--- a
+++ b/exseek/scripts/report.py
@@ -0,0 +1,640 @@
+#! /usr/bin/env python
+from __future__ import print_function
+import argparse, sys, os, errno
+import logging
+logging.basicConfig(level=logging.INFO, format='[%(asctime)s] [%(levelname)s] %(name)s: %(message)s')
+
+from scipy import interp
+from sklearn.metrics import roc_curve, roc_auc_score
+from sklearn.preprocessing import RobustScaler
+import numpy as np
+import pandas as pd
+
+command_handlers = {}
+def command_handler(f):
+    command_handlers[f.__name__] = f
+    return f
+
+def _compare_feature_selection_params(input_dir):
+    from tqdm import tqdm
+    import pandas as pd
+    import matplotlib.pyplot as plt
+
+    records = []
+    pbar = tqdm(unit='directory')
+    for compare_group in os.listdir(input_dir):
+        for path in os.listdir(os.path.join(input_dir, compare_group)):
+            classifier, n_features, selector, resample_method  = path.split('.')
+            record = {
+                'compare_group': compare_group,
+                'classifier': classifier,
+                'n_features': n_features,
+                'selector': selector,
+                'resample_method': resample_method
+            }
+            metrics = pd.read_table(os.path.join(input_dir, compare_group, path, 'metrics.{}.txt'.format(resample_method)))
+            record['test_roc_auc_mean'] = metrics['test_roc_auc'].mean()
+            if resample_method == 'leave_one_out':
+                record['test_roc_auc_std'] = 0
+            elif resample_method == 'stratified_shuffle_split':
+                record['test_roc_auc_std'] = metrics['test_roc_auc'].std()
+            pbar.update(1)
+            records.append(record)
+    pbar.close()
+    records = pd.DataFrame.from_records(records)
+    records.loc[:, 'n_features'] = records.loc[:, 'n_features'].astype('int')
+    compare_groups = records.loc[:, 'compare_group'].unique()
+
+    figsize = 3.5
+    # Compare resample methods
+    fig, axes = plt.subplots(1, len(compare_groups), 
+                             figsize=(figsize*len(compare_groups), figsize),
+                             sharey=True, sharex=False)
+    for i, compare_group in enumerate(compare_groups):
+        if len(compare_groups) > 1:
+            ax = axes[i]
+        else:
+            ax = axes
+        sub_df = records.query('compare_group == "{}"'.format(compare_group))
+        pivot = sub_df.pivot_table(index=['classifier', 'n_features', 'selector'], 
+                  columns=['resample_method'], 
+                  values='test_roc_auc_mean')
+        ax.scatter(pivot.loc[:, 'leave_one_out'], pivot.loc[:, 'stratified_shuffle_split'], s=12)
+        ax.set_xlabel('AUROC (leave_one_out)')
+        ax.set_ylabel('AUROC (stratified_shuffle_split)')
+        ax.set_xlim(0.5, 1)
+        ax.set_ylim(0.5, 1)
+        ax.plot([0.5, 1], [0.5, 1], linestyle='dashed', color='gray', linewidth=0.8)
+        ax.set_title(compare_group)
+
+    # Compare classifiers
+    fig, axes = plt.subplots(1, len(compare_groups), 
+                             figsize=(figsize*len(compare_groups), figsize),
+                             sharey=True, sharex=False)
+    for i, compare_group in enumerate(compare_groups):
+        if len(compare_groups) > 1:
+            ax = axes[i]
+        else:
+            ax = axes
+        sub_df = records.query('compare_group == "{}"'.format(compare_group))
+        pivot = sub_df.pivot_table(index=['resample_method', 'n_features', 'selector'], 
+                  columns=['classifier'], 
+                  values='test_roc_auc_mean')
+        ax.scatter(pivot.loc[:, 'logistic_regression'], pivot.loc[:, 'random_forest'], s=12)
+        ax.set_xlabel('AUROC (logistic_regression)')
+        ax.set_ylabel('AUROC (random_forest)')
+        ax.set_xlim(0.5, 1)
+        ax.set_ylim(0.5, 1)
+        ax.plot([0.5, 1], [0.5, 1], linestyle='dashed', color='gray', linewidth=0.8)
+        ax.set_title(compare_group)
+
+    # Compare number of features
+    fig, axes = plt.subplots(1, len(compare_groups), 
+                             figsize=(figsize*len(compare_groups), figsize),
+                             sharey=False, sharex=False)
+    for i, compare_group in enumerate(compare_groups):
+        if len(compare_groups) > 1:
+            ax = axes[i]
+        else:
+            ax = axes
+        sub_df = records.query('compare_group == "{}"'.format(compare_group))
+        pivot = sub_df.pivot_table(index=['classifier', 'selector', 'resample_method'], 
+                  columns=['n_features'], 
+                  values='test_roc_auc_mean')
+        ax.plot(np.repeat(pivot.columns.values.reshape((-1, 1)), pivot.shape[0], axis=1),
+                pivot.values.T)
+        ax.set_ylim(0.5, 1)
+        ax.set_xlabel('Number of features')
+        ax.set_ylabel('AUROC')
+        ax.set_title(compare_group)
+
+    # Compare feature selection methods
+    fig, axes = plt.subplots(1, len(compare_groups), 
+                             figsize=(figsize*len(compare_groups), figsize),
+                             sharey=True, sharex=False)
+    for i, compare_group in enumerate(compare_groups):
+        if len(compare_groups) > 1:
+            ax = axes[i]
+        else:
+            ax = axes
+        sub_df = records.query('compare_group == "{}"'.format(compare_group))
+        pivot = sub_df.pivot_table(index=['classifier', 'n_features', 'resample_method'], 
+                  columns=['selector'], 
+                  values='test_roc_auc_mean')
+        ax.plot(np.repeat(pivot.columns.values.reshape((-1, 1)), pivot.shape[0], axis=1),
+                pivot.values.T)
+        ax.set_ylim(0.5, 1)
+        ax.set_xlabel('Feature selection method')
+        ax.set_ylabel('AUROC')
+        ax.set_title(compare_group)
+    return records
+
+@command_handler
+def compare_feature_selection_params(args):
+    _compare_feature_selection_params(args.input_dir)
+    logger.info('save plot: ' + args.output_file)
+    plt.savefig(args.output_file)
+
+def _compare_features(input_dir, datasets):
+    import pandas as pd
+    from tqdm import tqdm
+    import seaborn as sns
+
+    pbar = tqdm(unit='directory')
+    records = []
+    feature_matrices = {}
+    #feature_support_matrices = {}
+    feature_indicator_matrices = {}
+    for dataset in datasets:
+        cpm = pd.read_table('output/cpm_matrix/{}.txt'.format(dataset), index_col=0)
+        for compare_group in os.listdir(os.path.join(input_dir, dataset)):
+            feature_lists = {}
+            #feature_supports = {}
+            for path in os.listdir(os.path.join(input_dir, dataset, compare_group)):
+                classifier, n_features, selector, resample_method  = path.split('.')
+                if int(n_features) > 10:
+                    continue
+                if (classifier != 'random_forest') or (selector != 'robust'):
+                    continue
+                if resample_method != 'stratified_shuffle_split':
+                    continue
+                record = {
+                    'compare_group': compare_group,
+                    'classifier': classifier,
+                    'n_features': n_features,
+                    'selector': selector,
+                    'resample_method': resample_method,
+                    'dataset': dataset
+                }
+                # feature importance
+                feature_lists[n_features] = pd.read_table(os.path.join(input_dir, dataset, compare_group,
+                    path, 'feature_importances.txt'), header=None, index_col=0).iloc[:, 0]
+                feature_lists[n_features].index = feature_lists[n_features].index.astype('str')
+                # feature support
+                #with h5py.File(os.path.join(input_dir, dataset, compare_group,
+                #    path, 'evaluation.{}.h5'.format(resample_method)), 'r') as f:
+                #    feature_support = np.mean(f['feature_selection'][:], axis=0)
+                #    feature_support = pd.Series(feature_support, index=cpm.index.values)
+                #    feature_support = feature_support[feature_lists[n_features].index.values]
+                #    feature_supports[n_features] = feature_support
+                # read metrics
+                metrics = pd.read_table(os.path.join(input_dir, dataset, compare_group, 
+                    path, 'metrics.{}.txt'.format(resample_method)))
+                record['test_roc_auc_mean'] = metrics['test_roc_auc'].mean()
+                if resample_method == 'leave_one_out':
+                    record['test_roc_auc_std'] = 0
+                elif resample_method == 'stratified_shuffle_split':
+                    record['test_roc_auc_std'] = metrics['test_roc_auc'].std()
+                pbar.update(1)
+                records.append(record)
+            # feature union set
+            feature_set = reduce(np.union1d, [a.index.values for a in feature_lists.values()])
+            # build feature importance matrix
+            feature_matrix = pd.DataFrame(np.zeros((len(feature_set), len(feature_lists))),
+                                          index=feature_set, columns=list(feature_lists.keys()))
+            for n_features, feature_importance in feature_lists.items():
+                feature_matrix.loc[feature_importance.index.values, n_features] = feature_importance.values
+            feature_matrix.columns = feature_matrix.columns.astype('int')
+            feature_matrix.index = feature_matrix.index.astype('str')
+            feature_matrix = feature_matrix.loc[:, feature_matrix.columns.sort_values().values]
+                    
+            feature_matrices[(dataset, compare_group)] = feature_matrix
+            # build feature indicator matrix
+            feature_indicator_matrix = pd.DataFrame(np.zeros((len(feature_set), len(feature_lists))),
+                                          index=feature_set, columns=list(feature_lists.keys()))
+            for n_features, feature_importance in feature_lists.items():
+                feature_indicator_matrix.loc[feature_importance.index.values, n_features] = 1
+            feature_indicator_matrix.columns = feature_indicator_matrix.columns.astype('int')
+            feature_indicator_matrix = feature_indicator_matrix.loc[:, feature_indicator_matrix.columns.sort_values().values]
+            feature_indicator_matrices[(dataset, compare_group)] = feature_indicator_matrix
+            
+            if dataset in feature_fields:
+                feature_meta = feature_matrix.index.to_series().str.split('|', expand=True)
+                feature_meta.columns = feature_fields[dataset]
+                if 'transcript_id' in feature_fields[dataset]:
+                    feature_matrix.insert(
+                        0, 'gene_type', 
+                        transcript_table_by_transcript_id.loc[feature_meta['transcript_id'].values, 'gene_type'].values)
+                    feature_matrix.insert(
+                        0, 'gene_name', 
+                        transcript_table_by_transcript_id.loc[feature_meta['transcript_id'].values, 'gene_name'].values)
+                    
+                elif 'gene_id' in feature_fields[dataset]:
+                    feature_matrix.insert(
+                        0, 'gene_type', 
+                        transcript_table_by_gene_id.loc[feature_meta['gene_id'].values, 'gene_type'].values)
+                    feature_matrix.insert(
+                        0, 'gene_name', 
+                        transcript_table_by_gene_id.loc[feature_meta['gene_id'].values, 'gene_name'].values)
+                elif 'transcript_name' in feature_fields[dataset]:
+                    feature_matrix.insert(
+                        0, 'gene_type', 
+                        transcript_table_by_transcript_name.loc[feature_meta['transcript_name'].values, 'gene_type'].values)
+                    feature_matrix.insert(
+                        0, 'gene_name', 
+                        transcript_table_by_transcript_name.loc[feature_meta['transcript_name'].values, 'gene_name'].values)
+                    
+                feature_indicator_matrix.index = feature_matrix.loc[:, 'gene_name'].values + '|' + feature_matrix.loc[:, 'gene_type'].values
+            # build feature support matrix
+            #feature_support_matrix = pd.DataFrame(np.zeros((len(feature_set), len(feature_lists))),
+            #                              index=feature_set, columns=list(feature_lists.keys()))
+            #for n_features, feature_support in feature_supports.items():
+            #    feature_support_matrix.loc[feature_support.index.values, n_features] = feature_support.values
+            #feature_support_matrix.columns = feature_support_matrix.columns.astype('int')
+            #feature_support_matrix = feature_matrix.loc[:, feature_support_matrix.columns.sort_values().values]
+            #feature_support_matrices[(dataset, compare_group)] = feature_support_matrix
+            fig, ax = plt.subplots(figsize=(6, 8))
+            sns.heatmap(feature_indicator_matrix,
+                        cmap=sns.light_palette('green', as_cmap=True), cbar=False, ax=ax, linewidth=1)
+            ax.set_xlabel('Number of features')
+            ax.set_ylabel('Fetures')
+            ax.set_title('{}, {}'.format(dataset, compare_group))
+
+            display(feature_matrix.style\
+                .background_gradient(cmap=sns.light_palette('green', as_cmap=True))\
+                .set_precision(2)\
+                .set_caption('{}, {}'.format(dataset, compare_group)))
+
+    pbar.close()
+    metrics = pd.DataFrame.from_records(records)
+    return metrics, feature_matrices, feature_indicator_matrices
+
+def plot_roc_curve_ci(y, is_train, predicted_scores, ax, title=None):
+    # ROC curve
+    n_splits = is_train.shape[0]
+    all_fprs = np.linspace(0, 1, 100)
+    roc_curves = np.zeros((n_splits, len(all_fprs), 2))
+    roc_aucs = np.zeros(n_splits)
+    for i in range(n_splits):
+        fpr, tpr, thresholds = roc_curve(y[~is_train[i]], predicted_scores[i, ~is_train[i]])
+        roc_aucs[i] = roc_auc_score(y[~is_train[i]], predicted_scores[i, ~is_train[i]])
+        roc_curves[i, :, 0] = all_fprs
+        roc_curves[i, :, 1] = interp(all_fprs, fpr, tpr)
+    roc_curves = pd.DataFrame(roc_curves.reshape((-1, 2)), columns=['fpr', 'tpr'])
+    sns.lineplot(x='fpr', y='tpr', data=roc_curves, ci='sd', ax=ax,
+                 label='Average ROAUC = {:.4f}'.format(roc_aucs.mean()))
+    #ax.plot(fpr, tpr, label='ROAUC = {:.4f}'.format(roc_auc_score(y_test, y_score[:, 1])))
+    #ax.plot([0, 1], [0, 1], linestyle='dashed')
+    ax.set_xlabel('False positive rate')
+    ax.set_ylabel('True positive rate')
+    ax.plot([0, 1], [0, 1], linestyle='dashed', color='gray')
+    if title:
+        ax.set_title(title)
+    ax.legend()
+
+def _plot_10_features(input_dir, datasets, use_log=False, scale=False, title=None):
+    pbar = tqdm_notebook(unit='directory')
+    for dataset in datasets:
+        sample_classes = pd.read_table('metadata/sample_classes.{}.txt'.format(groups[dataset]),
+                                       header=None, index_col=0).iloc[:, 0]
+        cpm = pd.read_table('output/cpm_matrix/{}.txt'.format(dataset), index_col=0)
+        if use_log:
+            cpm = np.log2(cpm + 0.001)
+        if scale:
+            X = RobustScaler().fit_transform(cpm.T.values).T
+        X = cpm.values
+        X = pd.DataFrame(X, index=cpm.index.values, columns=cpm.columns.values)
+        for compare_group in os.listdir(os.path.join(input_dir, dataset)):
+            for path in os.listdir(os.path.join(input_dir, dataset, compare_group)):
+                classifier, n_features, selector, resample_method  = path.split('.')
+                if int(n_features) != 10:
+                    continue
+                if (classifier != 'random_forest') or (selector != 'robust'):
+                    continue
+                if resample_method != 'stratified_shuffle_split':
+                    continue
+                record = {
+                    'compare_group': compare_group,
+                    'classifier': classifier,
+                    'n_features': n_features,
+                    'selector': selector,
+                    'resample_method': resample_method,
+                    'dataset': dataset
+                }
+                result_dir = os.path.join(input_dir, dataset, compare_group, path)
+                with h5py.File(os.path.join(result_dir, 'evaluation.{}.h5'.format(resample_method))) as f:
+                    train_index = f['train_index'][:]
+                    predicted_scores = f['predictions'][:]
+                    labels = f['labels'][:]
+                fig, ax = plt.subplots(figsize=(8, 8))
+                plot_roc_curve_ci(labels, train_index, predicted_scores, ax, 
+                                  title='{}, {}'.format(dataset, compare_group))
+                
+                features = pd.read_table(os.path.join(result_dir, 'features.txt'), header=None).iloc[:, 0].values
+                pbar.update(1)
+
+    pbar.close()
+
+
+def _evaluate_preprocess_methods(input_dirs, preprocess_methods, title=None):
+    records = []
+    pbar = tqdm_notebook(unit='directory')
+    for preprocess_method, input_dir in zip(preprocess_methods, input_dirs):
+        for compare_group in os.listdir(input_dir):
+            for path in os.listdir(os.path.join(input_dir, compare_group)):
+                classifier, n_features, selector, resample_method  = path.split('.')
+                if int(n_features) > 50:
+                    continue
+                if (classifier != 'random_forest') or (selector != 'robust'):
+                    continue
+                if resample_method != 'stratified_shuffle_split':
+                    continue
+                record = {
+                    'compare_group': compare_group,
+                    'classifier': classifier,
+                    'n_features': n_features,
+                    'selector': selector,
+                    'resample_method': resample_method,
+                    'preprocess_method': preprocess_method
+                }
+                metrics = pd.read_table(os.path.join(input_dir, compare_group, path, 'metrics.{}.txt'.format(resample_method)))
+                record['test_roc_auc_mean'] = metrics['test_roc_auc'].mean()
+                if resample_method == 'leave_one_out':
+                    record['test_roc_auc_std'] = 0
+                elif resample_method == 'stratified_shuffle_split':
+                    record['test_roc_auc_std'] = metrics['test_roc_auc'].std()
+                pbar.update(1)
+                records.append(record)
+    pbar.close()
+    records = pd.DataFrame.from_records(records)
+    records['n_features'] = records.loc[:, 'n_features'].astype(np.int32)
+    for compare_group, sub_df in records.groupby('compare_group'):
+        pivot = sub_df.pivot_table(
+            index='preprocess_method', columns='n_features', values='test_roc_auc_mean')
+        #print(pivot.iloc[:, 0])
+        #print(np.argsort(np.argsort(pivot.values, axis=0), axis=0)[:, 0])
+        mean_ranks = np.mean(pivot.shape[0] - np.argsort(np.argsort(pivot.values, axis=0), axis=0), axis=1)
+        mean_ranks = pd.Series(mean_ranks, index=pivot.index.values)
+        mean_ranks = mean_ranks.sort_values()
+        rename_index = ['{} (rank = {:.1f})'.format(name, value) for name, value in zip(mean_ranks.index, mean_ranks.values)]
+        rename_index = pd.Series(rename_index, index=mean_ranks.index.values)
+        sub_df = sub_df.copy()
+        sub_df['preprocess_method'] = rename_index[sub_df['preprocess_method'].values].values
+        sub_df['n_features'] = sub_df['n_features'].astype('int')
+        sub_df = sub_df.sort_values(['preprocess_method', 'n_features'], ascending=True)
+        sub_df['n_features'] = sub_df['n_features'].astype('str')
+        fig, ax = plt.subplots(figsize=(8, 8))                      
+        #sns.lineplot('n_features', 'test_roc_auc_mean', hue='preprocess_method', data=sub_df, 
+        #          ci=None, ax=ax, markers='o', hue_order=rename_index.values, sort=False)
+        for preprocess_method in rename_index.values:
+            tmp_df = sub_df[sub_df['preprocess_method'] == preprocess_method]
+            ax.plot(np.arange(tmp_df.shape[0]) + 1, tmp_df['test_roc_auc_mean'], label=preprocess_method)
+            ax.set_xticks(np.arange(tmp_df.shape[0]) + 1)
+            ax.set_xticklabels(tmp_df['n_features'])
+        ax.set_xlabel('Number of features')
+        ax.set_ylabel('Average AUROC')
+        if len(preprocess_methods) > 1:
+            ax.legend(title='Preprocess method', bbox_to_anchor=(1.04,0.5), 
+                      loc="center left", borderaxespad=0)
+        ax.set_ylim(0.5, 1)
+        if title:
+            ax.set_title(title + ', ' + compare_group)
+
+@command_handler
+def evaluate_preprocessing_methods(args):
+    _evaluate_preprocess_methods(args.input_dirs, args.precessing_methods)
+
+def bigwig_fetch(filename, chrom, start, end, dtype='float'):
+    import subprocess
+    p = subprocess.Popen(['bigWigToBedGraph', filename, 'stdout',
+                      '-chrom={}'.format(chrom), '-start={}'.format(start), '-end={}'.format(end)],
+                    stdout=subprocess.PIPE)
+    data = np.zeros(end - start, dtype=dtype)
+    for line in p.stdout:
+        line = str(line, encoding='ascii')
+        c = line.strip().split('\t')
+        data[(int(c[1]) - start):(int(c[2]) - start)] = float(c[3])
+    return data
+    
+
+def extract_feature_sequence(feature, genome_dir):
+    from pyfaidx import Fasta
+    from Bio.Seq import Seq
+
+    feature = line.split('\t')[0]
+    gene_id, gene_type, gene_name, domain_id, transcript_id, start, end = feature.split('|')
+    start = int(start)
+    end = int(end)
+    if gene_type == 'genomic':
+        gene_type = 'genome'
+    fasta = Fasta(os.path.join(args.genome_dir, 'fasta', gene_type + '.fa'))
+    if gene_type == 'genome':
+        chrom, gstart, gend, strand = gene_id.split('_')
+        gstart = int(gstart)
+        gend = int(gend)
+        seq = fasta[chrom][gstart:gend].seq
+        if strand == '-':
+            seq = str(Seq(seq).reverse_complement())
+    else:
+        seq = fasta[transcript_id][start:end].seq
+    seq = seq.upper()
+
+
+@command_handler
+def visualize_domains(args):
+    import numpy as np
+    import matplotlib
+    matplotlib.use('Agg')
+    import matplotlib.pyplot as plt
+    from matplotlib.backends.backend_pdf import PdfPages
+    from matplotlib.gridspec import GridSpec
+    import seaborn as sns
+    sns.set_style('white')
+    import pandas as pd
+    plt.rcParams['figure.dpi'] = 96
+    from tqdm import tqdm
+    from pykent import BigWigFile
+    from scipy.cluster.hierarchy import linkage, dendrogram
+    from pyfaidx import Fasta
+    from Bio.Seq import Seq
+    from call_peak import call_peaks
+
+    # read sample ids
+    #logger.info('read sample ids: ' + args.sample_ids)
+    #sample_ids = open(args.sample_ids_file, 'r').read().split()
+    logger.info('reads sample classes: ' + args.sample_classes)
+    sample_classes =  pd.read_table(args.sample_classes, sep='\t', index_col=0).iloc[:, 0]
+    sample_classes = sample_classes.sort_values()
+    sample_ids = sample_classes.index.values
+
+    # read features
+    features = pd.read_table(args.features, header=None).iloc[:, 0]
+    feature_info = features.str.split('|', expand=True)
+    feature_info.columns = ['gene_id', 'gene_type', 'gene_name', 'domain_id', 'transcript_id', 'start', 'end']
+    feature_info.index = features.values
+    # read count matrix to get read depth
+    #counts = pd.read_table(args.count_matrix, index_col=0)
+    #read_depth = counts.sum(axis=0)
+    #del counts
+    # read chrom sizes
+    #chrom_sizes = pd.read_table(args.chrom_sizes, sep='\t', index_col=0, header=None).iloc[:, 0]
+
+    with PdfPages(args.output_file) as pdf:
+        for feature_name, feature in tqdm(feature_info.iterrows(), unit='feature'):
+            #logger.info('plot feature: {}'.format(feature_name))
+            if feature['gene_type'] == 'genomic':
+                chrom, start, end, strand = feature['gene_id'].split('_')
+                start = int(start)
+                end = int(end)
+                bigwig_file = os.path.join(args.output_dir, 'bigwig', '{{0}}.genome.{0}.bigWig'.format(strand))
+            elif feature['gene_type'] in ('piRNA', 'miRNA'):
+                continue
+            else:
+                start = int(feature['start'])
+                end = int(feature['end'])
+                chrom = feature.transcript_id
+                bigwig_file = os.path.join(args.output_dir, 'tbigwig_normalized', '{0}.transcriptome.bigWig')
+            # read coverage from BigWig files
+            coverage = None
+            for i, sample_id in enumerate(sample_ids):
+                bwf = BigWigFile(bigwig_file.format(sample_id))
+                if coverage is None:
+                    # interval to display coverage
+                    chrom_size = bwf.get_chrom_size(feature['transcript_id'])
+                    if chrom_size == 0:
+                        raise ValueError('cannot find transcript id {} in bigwig'.format(feature['transcript_id']))
+                    view_start = max(start - args.flanking, 0)
+                    view_end = min(end + args.flanking, chrom_size)
+                    coverage = np.zeros((len(sample_ids), view_end - view_start), dtype='float')
+                    logger.info('create_coverage_matrix: ({}, {})'.format(*coverage.shape))
+                #logger.info('bigWigQuery: {}:{}-{}'.format(chrom, view_start, view_end))
+                values = bwf.query(chrom, view_start, view_end, fillna=0)
+                del bwf
+                if values is not None:
+                    coverage[i] = values
+                #coverage[i] = bigwig_fetch(bigwig_file.format(sample_id), chrom, view_start, view_end, dtype='int')
+                # normalize coverage by read depth
+                #coverage[i] *= 1e6/read_depth[sample_id]
+                # log2 transformation
+                coverage[i] = np.log2(coverage[i] + 1)
+            
+            # get sequence
+            gene_type = feature['gene_type']
+            if gene_type == 'genomic':
+                gene_type = 'genome'
+            fasta = Fasta(os.path.join(args.genome_dir, 'fasta', gene_type + '.fa'))
+            seq = fasta[feature['transcript_id']][view_start:view_end].seq
+            if (gene_type == 'genome') and (strand == '-'):
+                seq = str(Seq(seq).reverse_complement())
+            seq = seq.upper()
+
+            # draw heatmap
+            '''
+            plot_data = pd.DataFrame(coverage)
+            cmap = sns.light_palette('blue', as_cmap=True, n_colors=6)
+            g = sns.clustermap(plot_data, figsize=(20, 8), col_cluster=False, row_colors=None, cmap='Blues')
+            g.ax_heatmap.set_yticklabels([])
+            g.ax_heatmap.set_yticks([])
+            xticks = np.arange(0, coverage.shape[1], 10)
+            g.ax_heatmap.set_xticks(xticks)
+            g.ax_heatmap.set_xticklabels(xticks, rotation=0)
+            g.ax_heatmap.vlines(x=domain_start, ymin=0, ymax=g.ax_heatmap.get_ylim()[0], linestyle='dashed', linewidth=1.0)
+            g.ax_heatmap.vlines(x=domain_end, ymin=0, ymax=g.ax_heatmap.get_ylim()[0], linestyle='dashed', linewidth=1.0)
+            g.ax_heatmap.set_title(feature_name)
+            '''
+            # hierarchical clustering
+            order = np.arange(coverage.shape[0], dtype='int')
+            for label in np.unique(sample_classes):
+                mask = (sample_classes == label)
+                Z = linkage(coverage[mask], 'single')
+                R = dendrogram(Z, no_plot=True, labels=order[mask])
+                order[mask] = R['ivl']
+            sample_classes = sample_classes.iloc[order]
+            coverage = coverage[order]
+
+            plt.rcParams['xtick.minor.visible'] = True
+            plt.rcParams['xtick.minor.size'] = 4
+            plt.rcParams['xtick.bottom'] = True
+            plt.rcParams['xtick.labelsize'] = 8
+            fig = plt.figure(figsize=(20, 6))
+            gs = GridSpec(4, 3, figure=fig, width_ratios=[0.95, 0.03, 0.02], height_ratios=[0.6, 0.15, 0.1, 0.1], hspace=0.2, wspace=0.15)
+            #fig, axes = plt.subplots(1, 2, figsize=(20, 3), sharey=True, 
+            #    gridspec_kw={'width_ratios': [0.98, 0.02], 'hspace': 0})
+            ax_heatmap = plt.subplot(gs[0, 0])
+            ax_colorbar = plt.subplot(gs[0, 1])
+            ax_label = plt.subplot(gs[0, 2])
+            ax_line = plt.subplot(gs[1, 0])
+            ax_domain = plt.subplot(gs[2, 0])
+            ax_refined_domain = plt.subplot(gs[3, 0])
+
+            p = ax_heatmap.pcolormesh(coverage, cmap='Blues')
+            ax_heatmap.set_xticks(np.arange(coverage.shape[1]) + 0.5, minor=True)
+            ax_heatmap.set_xlim(0, coverage.shape[1])
+            xticks = ax_heatmap.get_xticks()
+            ax_heatmap.set_xticks(xticks + 0.5)
+            ax_heatmap.set_xticklabels(xticks.astype('int'))
+            ax_heatmap.set_title(feature_name)
+
+            fig.colorbar(p, cax=ax_colorbar, use_gridspec=False, orientation='vertical')
+
+            for label in sample_classes.unique():
+                ax_label.barh(y=np.arange(coverage.shape[0]), width=(sample_classes == label).astype('int'), height=1,
+                    edgecolor='none', label=label)
+            ax_label.set_xlim(0, 1)
+            ax_label.set_ylim(0, coverage.shape[0])
+            ax_label.tick_params(labelbottom=False, bottom=False)
+            ax_label.set_xticks([])
+            ax_label.set_yticks([])
+            ax_label.legend(title='Class', bbox_to_anchor=(1.1, 0.5), loc="center left", borderaxespad=0)
+
+            ax_line.fill_between(np.arange(coverage.shape[1]), coverage.mean(axis=0), step='pre', alpha=0.9)
+            ax_line.set_xlim(0, coverage.shape[1])
+            #ax_line.set_xticks(np.arange(coverage.shape[1]) + 0.5, minor=True)
+            #xticks = ax_line.get_xticks()
+            #ax_line.set_xticks(xticks + 0.5)
+            #ax_line.set_xticklabels(xticks.astype('int'))
+            ax_line.set_xticks(np.arange(coverage.shape[1]) + 0.5)
+            ax_line.set_xticks([], minor=True)
+            ax_line.set_xticklabels(list(seq))
+            ax_line.set_ylim(0, ax_line.get_ylim()[1])
+            #ax_line.vlines(x=start - view_start + 0.5, ymin=0, ymax=ax_line.get_ylim()[1], linestyle='dashed', linewidth=1.0)
+            #ax_line.vlines(x=end - view_start + 0.5, ymin=0, ymax=ax_line.get_ylim()[1], linestyle='dashed', linewidth=1.0)
+
+            ax_domain.hlines(y=0.5, xmin=start - view_start, xmax=end - view_start, linewidth=5, color='C0')
+            ax_domain.set_ylim(0, 1)
+            ax_domain.set_ylabel('Domain')
+            ax_domain.set_yticks([])
+            ax_domain.set_xticks([])
+            ax_domain.set_xticks(np.arange(coverage.shape[1]) + 0.5, minor=True)
+            ax_domain.set_xlim(0, coverage.shape[1])
+            ax_domain.spines['top'].set_visible(False)
+            ax_domain.spines['right'].set_visible(False)
+
+            coverage_mean = coverage.mean(axis=0)
+            for _, peak_start, peak_end in call_peaks([coverage_mean], min_length=10):
+                ax_refined_domain.hlines(y=0.5, xmin=peak_start, xmax=peak_end, linewidth=5, color='C0')
+            ax_refined_domain.set_ylim(0, 1)
+            ax_refined_domain.set_ylabel('Refined')
+            ax_refined_domain.set_yticks([])
+            ax_refined_domain.set_xticks([])
+            ax_refined_domain.set_xticks(np.arange(coverage.shape[1]) + 0.5, minor=True)
+            ax_refined_domain.set_xlim(0, coverage.shape[1])
+            ax_refined_domain.spines['top'].set_visible(False)
+            ax_refined_domain.spines['right'].set_visible(False)
+
+            fig.tight_layout()
+            # save plot
+            pdf.savefig(fig)
+            plt.close()
+
+
+if __name__ == '__main__':
+    main_parser = argparse.ArgumentParser(description='Preprocessing module')
+    subparsers = main_parser.add_subparsers(dest='command')
+
+    parser = subparsers.add_parser('visualize_domains', help='plot read coverage of domains as heatmaps')
+    parser.add_argument('--sample-classes', type=str, required=True, help='e.g. {data_dir}/sample_classes.txt')
+    parser.add_argument('--output-dir', type=str, required=True, help='e.g. output/scirep')
+    parser.add_argument('--features', type=str, required=True, help='list of selected features')
+    #parser.add_argument('--count-matrix', type=str, required=True, help='count matrix')
+    parser.add_argument('--output-file', '-o', type=str, required=True, help='output PDF file')
+    parser.add_argument('--flanking', type=int, default=20, help='flanking length for genomic domains')
+    parser.add_argument('--genome-dir', type=str, required=True, help='e.g. genome/hg38')
+    
+    args = main_parser.parse_args()
+    if not args.command:
+        main_parser.print_help()
+        sys.exit(1)
+    logger = logging.getLogger('report.' + args.command)
+
+    command_handlers.get(args.command)(args)
\ No newline at end of file