Update single-cell simulation

Author: jaydu1

paper/simu_nb/Plot.ipynb

@@ -3,9 +3,11 @@ do
     Rscript simu_nb_data.R ${isimu}
     cp simu_nb_data.R data/simu_100${isimu}/
     cp simu_nb_fit.R data/simu_100${isimu}/
-    Rscript --no-save --no-restore --verbose simu_nb_fit.R ${isimu} &> out${isimu}_fit.txt
-    mv out${isimu}_fit.txt results/
-    
+    for seed in $(seq 0 49);
+    do
+        Rscript --no-save --no-restore --verbose simu_nb_fit.R ${isimu} ${seed} &> out${isimu}_fit_${seed}.txt
+        mv out${isimu}_fit_${seed}.txt results/
+    done
     python simu_nb_plot.py ${isimu} &> out${isimu}_plot.txt
     mv out${isimu}_plot.txt results/
-done
+done

paper/simu_nb/simu_nb.sh

@@ -1,4 +1,6 @@
 library(rhdf5)
+
+
 # Load Python and R functions for basic Wilcoxon Tests, DESeq2, and CocoA-Diff
 # and causarray, cinemaot
 require(reticulate)
@@ -15,8 +17,8 @@ n <- 200
 c = 0.1
 alpha = 0.1
 
-ind <- '_d_1_r_4_noise_0.1'
-num_r <- 1
+ind <- '_d_1_r_4_noise_0.2'
+num_r <- 1; d <- 1
 
 args = commandArgs(trailingOnly=TRUE)
 if(length(args)==0){
@@ -38,17 +40,15 @@ if((ind != '') && grepl( 'r', ind, fixed = TRUE)){
 }
 
 
-
+seeds <- as.integer(args[2])
 for(n in c(100, 500, 1000, 5000)){
     path_result <- sprintf(paste0(path_base,'results/simu_%d%s/'), n, ind)
     dir.create(path_result, recursive=TRUE, showWarnings = FALSE)
 
-    for (seed in 0:49) {
+    for (seed in c(seeds)) {
         path_data <- sprintf(paste0(path_base,'data/simu_%d%s/simu_data_%d.h5'), n, ind, seed)
 
-        if(file.exists(sprintf('%scausarray_r_%d_%d.csv', path_result, 6, seed))){
-            next
-        }
+        cat(n, seed, '...\n')
 
         Y <- t(h5read(path_data, '/Y'))
         metadata <- t(h5read(path_data, '/metadata'))
@@ -72,12 +72,13 @@ for(n in c(100, 500, 1000, 5000)){
         possibleError <- tryCatch(
             {
         # Wilcoxon Tests ----
-        res.wilc <- run_wilcoxon(Y, metadata=scaleW, raw=T)
+        res.wilc <- run_wilcoxon(Y, metadata=scaleW)
 
         # DESeq ----
         res.DESeq <- run_DESeq(Y, metadata=scaleW, return_ds=TRUE)
         dds <- res.DESeq[[2]]
         res.DESeq <- res.DESeq[[1]]
+        # causarray$comp_stat(theta, res.DESeq$padj<0.1, 0.1)
 
         # CocoA-Diff ----
         cocoa <- run_cocoa(sc = Y, indvs=metadata[,2], metadata=scaleW, cocoAWriteName=sprintf('%stmp_cocoa/tmp_%d', path_result, seed))
@@ -97,20 +98,34 @@ for(n in c(100, 500, 1000, 5000)){
         W.cinemaotw <- res.cinemaotw[[2]]$W
         res.cinemaotw <- res.cinemaotw[[1]]
 
+
         # Save confounder estimation results
         write.csv(cf.cocoa, sprintf('%scocoa_cf_%d.csv', path_result, seed))        
         write.csv(cf.cinemaot, sprintf('%scinemaot_cf_%d.csv', path_result, seed))
         write.csv(W.cinemaot, sprintf('%scinemaot_W_%d.csv', path_result, seed))
         write.csv(cf.cinemaotw, sprintf('%scinemaotw_cf_%d.csv', path_result, seed))
         write.csv(W.cinemaotw, sprintf('%scinemaotw_W_%d.csv', path_result, seed))
 
+        
         # Save test results
         write.csv(res.wilc, sprintf('%swilc_%d.csv', path_result, seed))
         write.csv(res.DESeq, sprintf('%sDESeq_%d.csv', path_result, seed))
         write.csv(res.cocoa, sprintf('%scocoa_%d.csv', path_result, seed))
         write.csv(res.cinemaot, sprintf('%scinemaot_%d.csv', path_result, seed))
         write.csv(res.cinemaotw, sprintf('%scinemaotw_%d.csv', path_result, seed))
 
+
+        # mixscape ----
+        res.mixscape <- run_mixscape(Y, A, raw=TRUE)
+        cf.mixscape <- res.mixscape[[2]]$Y_hat_0
+        W.mixscape <- res.mixscape[[2]]$W
+        res.mixscape <- res.mixscape[[1]]
+
+        write.csv(cf.mixscape, sprintf('%smixscape_cf_%d.csv', path_result, seed))
+        write.csv(W.mixscape, sprintf('%smixscape_W_%d.csv', path_result, seed))
+        write.csv(res.mixscape, sprintf('%smixscape_%d.csv', path_result, seed))
+
+
         for(r_hat in c(2,4,6)){
             # RUV
             ruv <- run_ruv(Y, metadata=scaleW, r_hat)
@@ -140,7 +155,7 @@ for(n in c(100, 500, 1000, 5000)){
 
             res.causarray <- run_causarray(Y, scaleW[,-ncol(scaleW)], A, 
                 fdx=T, r=r_hat, glm_alpha=.5, shrinkage=T)
-            cf.causarray <- log1p(res.causarray[[2]]$Y_hat_0)
+            cf.causarray <- log1p(res.causarray[[2]]$Y_hat[,,1,1])
             W.causarray <- res.causarray[[2]]$W
             res.causarray <- res.causarray[[1]]
 

paper/simu_nb/simu_nb_fit.R

@@ -31,7 +31,7 @@ def legend_title_left(leg):
 else:
     ind = ''
 r_list = [2,4,6]
-method_list = ['wilc', 'DESeq', 'cocoa', 'cinemaot', 'cinemaotw'] \
+method_list = ['wilc', 'DESeq', 'cocoa', 'cinemaot', 'cinemaotw', 'mixscape'] \
     + ['ruv_r_{}'.format(r) for r in r_list] \
     + ['ruv3nb_r_{}'.format(r) for r in r_list] \
     + ['causarray_r_{}'.format(r) for r in r_list]
@@ -83,6 +83,27 @@ def legend_title_left(leg):
 df_res.reset_index(drop=True, inplace=True)
 df_res.to_csv(path_base+'results/result{}_test.csv'.format(ind))
 
+
+sns.set(font_scale=1.2)
+fig, axes = plt.subplots(1,2, figsize=(10,4), sharex=True, sharey=False)
+for j, metric in enumerate(['FDR', 'power']):
+    sns.boxplot(data=df_res, x='n', y=metric, hue='method', ax=axes[j])
+
+axes[0].axhline(alpha, color='r', linestyle='--')
+lines_labels = [ax.get_legend_handles_labels() for ax in [axes[1]]]
+handles, labels = [sum(lol, []) for lol in zip(*lines_labels)]
+
+axes[0].get_legend().remove()
+axes[1].get_legend().remove()
+legend = fig.legend(handles=handles, labels=labels,
+                    loc=9, ncol=5, title=None, frameon=False)
+legend_title_left(legend)
+
+fig.tight_layout()
+fig.subplots_adjust(top=0.85)
+
+plt.savefig(path_base+'results/simu{}_res.pdf'.format(ind), bbox_inches='tight', pad_inches=0, dpi=300)
+
 print(df_res.groupby(['n','method'])[['typeI_err', 'FDR', 'power', 'FDX', 'num_dis']].median())
 
 
@@ -93,7 +114,7 @@ def legend_title_left(leg):
 
 
 r_list = [2,4,6]
-method_list = ['cocoa', 'cinemaot', 'cinemaotw']  \
+method_list = ['cocoa', 'cinemaot', 'cinemaotw', 'mixscape']  \
     + ['ruv_r_{}'.format(r) for r in r_list] \
     + ['ruv3nb_r_{}'.format(r) for r in r_list] \
     + ['causarray_r_{}'.format(r) for r in r_list]
@@ -146,4 +167,67 @@ def legend_title_left(leg):
     df_res = pd.concat([df_res, df], axis=0)
 
 df_res.reset_index(drop=True, inplace=True)
-df_res.to_csv(path_base+'results/result{}_deconfound.csv'.format(ind))
+df_res.to_csv(path_base+'results/result{}_deconfound.csv'.format(ind))
+
+
+
+
+
+df_test = pd.read_csv(path_base+'results/result{}_test.csv'.format(ind)).rename({'FDR':'FPR', 'power':'TPR'}, axis=1)
+df_cf = pd.read_csv(path_base+'results/result{}_deconfound.csv'.format(ind))
+
+r_list = [2,4,6]
+method_name = {
+    'wilc':'Wilcoxon', 'DESeq':'DESeq2', 'cocoa':'CoCoA', 'cinemaot':'CINEMA-OT', 'cinemaotw':'CINEMA-OT-W', 'mixscape':'Mixscape',
+    }
+    
+method_name = reduce(lambda a, b: dict(a, **b), 
+    [{'ruv_r_{}'.format(r):'RUV $r={}$'.format(r) for r in r_list}, 
+      {'ruv3nb_r_{}'.format(r):'RUV-III-NB $r={}$'.format(r) for r in r_list}, 
+      {'causarray_r_{}'.format(r):'causarray $r={}$'.format(r) for r in r_list}
+    ])
+
+df_test = df_test[df_test['method'].isin(method_name.keys())]
+df_cf = df_cf[df_cf['method'].isin(method_name.keys())]
+df_test['method'] = df_test['method'].map(method_name)
+df_cf['method'] = df_cf['method'].map(method_name)
+
+df_test = df_test[df_test['n'].isin(n_list)]
+df_cf = df_cf[df_cf['n'].isin(n_list)]
+
+method_list = method_name.values()#df_test['method'].unique()
+# palette = sns.color_palette()[:len(method_list)]
+palette = reduce(lambda l1, l2: l1+l2, [sns.color_palette(name)[:len(r_list)*2:2] for name in ['Reds', 'Greens', 'Blues']])
+hue_order = {i:c for i,c in zip(method_list, palette) }
+
+sns.set(font_scale=1.3)
+fig, axes = plt.subplots(1,4, figsize=(16,5), sharex=False, sharey=False)
+for j, metric in enumerate(['FPR', 'TPR']):
+    sns.boxplot(data=df_test, x='n', y=metric, hue='method', hue_order=hue_order,
+        ax=axes[j+2], palette=palette, showfliers=False)
+
+for j, metric in enumerate(['ARI', 'ASW']):
+    sns.boxplot(data=df_cf, x='n', y=metric, hue='method', hue_order=hue_order,
+        ax=axes[j], palette=palette, showfliers=False)
+
+axes[2].axhline(0.1, color='r', linestyle='--')
+lines_labels = [ax.get_legend_handles_labels() for ax in [axes[1]]]
+handles, labels = [sum(lol, []) for lol in zip(*lines_labels)]
+handles = [mlines.Line2D([], [], linestyle='None')] * 3 + handles[::3] + handles[1::3] + handles[2::3]
+labels = ['RUV', 'RUV-III-NB', 'causarray',
+            '$r=2$', '$r=2$', '$r=2$',
+            '$r=4$', '$r=4$', '$r=4$',
+            '$r=6$', '$r=6$', '$r=6$']
+
+for j in range(4):
+    axes[j].get_legend().remove()
+    axes[j].tick_params(axis='both', which='major', labelsize=10)
+    axes[j].set_xlabel('Sample size $n$')
+legend = fig.legend(handles=handles, labels=labels,
+                    loc=9, ncol=4, title=None, frameon=False)           
+legend_title_left(legend)
+
+fig.tight_layout()
+fig.subplots_adjust(top=0.78)
+
+plt.savefig(path_base + 'results/simu_nb_r{}.pdf'.format(ind), bbox_inches='tight', pad_inches=0, dpi=300)