evalsim_sceps.py

import argparse, glob, sys
import pandas as pd 
import numpy as np 
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import re
import scipy as sp
import scipy.stats

from tqdm import tqdm


def robust_mean(vec):
    """
    calculate the mean with extreme ends removed
    """

    if vec.shape[0] == 0:
        return np.nan

    q1 = np.percentile(vec, 5)
    q3 = np.percentile(vec, 95)

    iqr = q3 - q1
    lower_bound = q1 - 1.5 * iqr
    upper_bound = q3 + 1.5 * iqr

    vec_ = vec[(vec >= lower_bound) & (vec <= upper_bound)]

    return np.mean(vec_)


def robust_se(vec):
    """
    calculate the mean with extreme ends removed
    """

    if vec.shape[0] == 0:
        return np.nan

    q1 = np.percentile(vec, 25)
    q3 = np.percentile(vec, 75)

    iqr = q3 - q1
    lower_bound = q1 - 1.5 * iqr
    upper_bound = q3 + 1.5 * iqr

    vec_ = vec[(vec >= lower_bound) & (vec <= upper_bound)]

    return np.std(vec_) / np.sqrt(vec_.shape[0])

def main():

    args = get_command_line()

    # load in all simulation resultsj
    df_all = []
    for j, prefix in enumerate(args.prefix):
        df_sim = []
        all_out = glob.glob('{}*.sceps.weighted.txt.gz'.format(prefix))
        for fnm in tqdm(all_out):
            sim_idx = fnm.split('/')[-1].replace('sim_out.', '')
            sim_idx = sim_idx.replace('.sceps.weighted.txt.gz', '')
            sim_idx = int(sim_idx)

            shrinkage = 0.0

            try:
                df = pd.read_table(fnm, sep='\t')
            except Exception as inst:
                continue
            df['IDX'] = sim_idx
            df['SHRINKAGE'] = shrinkage
            df_sim.append(df)
        if len(df_sim) > 0:
            df_sim = pd.concat(df_sim)
            df_sim['SIMULATED_DIFF'] = args.simulated_diff[j]
            df_all.append(df_sim)
    df_all = pd.concat(df_all)

    # get all shrinkage
    all_shrinkage = np.sort(pd.unique(df_all['SHRINKAGE']))

    # load scdrs result
    df_scdrs = pd.read_table('/gstore/project/statgen_methods/sceps/simulation/004_simulations_core/eval/scdrs.out.txt')

    # summarize results
    sim_sum = []
    avg_ratio = []
    for sim_diff in args.simulated_diff:

        for shrinkage in all_shrinkage:


            df_tmp = df_all[(df_all['SIMULATED_DIFF']==sim_diff) & (df_all['SHRINKAGE']==shrinkage)]
            prop_reject = np.sum(df_tmp['P_Z_WEIGHTED_DIFF']<0.05) / df_tmp.shape[0]

            mean_diff = np.mean(df_tmp['WEIGHTED_DIFF'].values)
            se_diff = np.std(df_tmp['WEIGHTED_DIFF'].values) / np.sqrt(df_tmp['WEIGHTED_DIFF'].shape[0])

            df_scdrs_tmp = df_scdrs[df_scdrs['barcode'].isin(df_tmp['CELL'])]
            #prop_reject_scdrs = np.sum(df_scdrs_tmp['pval']<0.05) / df_scdrs_tmp.shape[0]
            mean_scdrs = np.mean(df_scdrs_tmp['raw_score'])
            se_scdrs = np.std(df_scdrs_tmp['raw_score']) / np.sqrt(df_scdrs_tmp.shape[0])

            sim_sum.append((sim_diff, shrinkage, mean_diff, se_diff, prop_reject,
                            mean_scdrs, se_scdrs))
    sim_sum = pd.DataFrame(sim_sum)
    sim_sum.columns = ['SIM_DIFF', 'SHRINKAGE', 'MEAN', 'SE', 'PROP_REJECT',
                       'MEAN_scDRS', 'SE_scDRS']

    # sceps
    _, ax = plt.subplots(figsize=(6,4))
    num_shrinkage = all_shrinkage.shape[0] - 1
    width = np.median(args.simulated_diff)/15
    for i, shrinkage in enumerate(all_shrinkage):
        sim_sum_tmp = sim_sum[sim_sum['SHRINKAGE']==shrinkage]
        ax.errorbar(x=sim_sum_tmp['SIM_DIFF']+(i-num_shrinkage/2.0)*(2*width)/(num_shrinkage + 1e-16), \
                    y=sim_sum_tmp['MEAN'],  yerr=2.0*sim_sum_tmp['SE'],
                    ls='none', fmt='o', label='shrinkage={}'.format(shrinkage))
    for i in range(len(args.simulated_diff)):
        ax.plot([(args.simulated_diff[i]-width), (args.simulated_diff[i]+width)],
                [args.simulated_diff[i], args.simulated_diff[i]], 'k--', alpha=0.75, zorder=0)
    
    ax.set_xlabel('simulated scEPS statistics', fontsize=12)
    ax.set_ylabel('estimated scEPS statistics', fontsize=12)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    plt.savefig('./eval/scEPS_sim_eval_sceps_score.png', bbox_inches='tight', dpi=300)

    # scdrs
    """
    fig, ax = plt.subplots(figsize=(6,4))
    sim_sum_tmp = sim_sum[sim_sum['SHRINKAGE']==0.0]
    ax.errorbar(x=sim_sum['SIM_DIFF'], y=sim_sum['MEAN_scDRS'],  yerr=1.96*sim_sum['SE_scDRS'],
                ls='none', fmt='o', label='simulation results')
    ax.legend(loc='best')
    ax.set_xlabel('scDRS score', fontsize=14)
    ax.set_ylabel('estimated scEPS score', fontsize=14)
    plt.savefig('./eval/scEPS_sim_eval_scdrs_score.png', bbox_inches='tight', dpi=300)
    """

    df_pval = sim_sum[['SIM_DIFF', 'PROP_REJECT', 'SHRINKAGE']].melt(id_vars=['SIM_DIFF', 'SHRINKAGE'])

    fig, ax = plt.subplots(figsize=(6,4))
    sns.barplot(data=df_pval, x='SIM_DIFF', y='value', hue='SHRINKAGE')
    ax.set_xlabel('simulated scEPS score', fontsize=14)
    ax.set_ylabel('prop. of rejection (p < 0.05)', fontsize=14)
    plt.savefig('./eval/scEPS_sim_eval_pval.png', bbox_inches='tight', dpi=300)


def get_command_line():
 
    parser = argparse.ArgumentParser(description="Generate simulated data")

    # Create the parser
    parser = argparse.ArgumentParser(description="Generate simulated data")

    # Add the arguments
    parser.add_argument('--prefix', type=str, required=False, nargs='+',
                        help='Output directory')

    parser.add_argument('--simulated-diff', type=float, required=False, nargs='+',
                        help='Output directory')

    parser.add_argument('--out', type=str, required=False,
                        help='Output directory')

    # Execute the parse_args() method
    args = parser.parse_args()
    
    return args


if __name__ == '__main__':
    main()