reczko@estia:/data/images/proton/run164$ /home/moulos/Rbase/bin/R library(metaseqR) the.path="/data/images/proton/Cereghini/run142" #26042015 metaseqr( sample.list=file.path(the.path,"targets_run142_woSCR22.txt.txt"), contrast=c("WT_vs_HET"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_wt_vs_het_run142_woSCR22"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.3, fig.format=c("png","pdf") ) #27022015 metaseqr( sample.list=file.path(the.path,"targets_run142.txt"), contrast=c("WT_vs_HET"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_wt_vs_het_run142"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.3, fig.format=c("png","pdf") ) metaseqr( sample.list=file.path(the.path,"targets_run128.txt"), contrast=c("WT_vs_MUT"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_wt_vs_mut_run128"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, fig.format=c("png","pdf") ) metaseqr( sample.list=file.path(the.path,"targets_run107.txt"), contrast=c("WT_vs_HET"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_wt_vs_het_run107"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.3, fig.format=c("png","pdf") ) metaseqr( sample.list=file.path(the.path,"targets_run20.txt"), contrast=c("WT_vs_TR"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_wt_vs_tr_run20"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.3, fig.format=c("png","pdf") ) reczko@estia:/data/images/proton/Cereghini$ ls /data/images/proton/Cereghini/metaseq*/lists/metaseqr_sig_out*.txt.gz r107=read.table("/data/images/proton/Cereghini/metaseqr_wt_vs_het_run107/lists/metaseqr_sig_out_WT_vs_HET.txt.gz",header=TRUE) r128=read.table("/data/images/proton/Cereghini/metaseqr_wt_vs_mut_run128/lists/metaseqr_sig_out_WT_vs_MUT.txt.gz",header=TRUE) r142=read.table("/data/images/proton/Cereghini/metaseqr_wt_vs_mut_run142/lists/metaseqr_sig_out_WT_vs_HET.txt.gz",header=TRUE) r20=read.table("/data/images/proton/Cereghini/metaseqr_wt_vs_tr_run20/lists/metaseqr_sig_out_WT_vs_TR.txt.gz",header=TRUE) ar107=read.table("/data/images/proton/Cereghini/run142/metaseqr_wt_vs_het_run107/lists/metaseqr_all_out_WT_vs_HET.txt.gz",header=TRUE) ar128=read.table("/data/images/proton/Cereghini/run142/metaseqr_wt_vs_mut_run128/lists/metaseqr_all_out_WT_vs_MUT.txt.gz",header=TRUE) ar142=read.table("/data/images/proton/Cereghini/run142/metaseqr_wt_vs_mut_run142/lists/metaseqr_all_out_WT_vs_HET.txt.gz",header=TRUE) ar20=read.table("/data/images/proton/Cereghini/run142/metaseqr_wt_vs_tr_run20/lists/metaseqr_all_out_WT_vs_TR.txt.gz",header=TRUE) genes=levels(rbind(r20[,1:10],r107[,1:10],r128[,1:10],r142[,1:10])) x=merge(r20[,1:23], r107[,1:23], by='gene_name', all=T,suffixes = c(".r20",".r107")) y=merge(x, r142[,1:23], by='gene_name', all=T,suffixes = c("",".r142")) x=merge(y, r128[,1:23], by='gene_name', all=T,suffixes = c("",".r128")) write.table(x,"foo") str(x) 'data.frame': 2320 obs. of 89 variables: $ gene_name : Factor w/ 2320 levels "Aadat","Abhd2",..: 1 2 3 4 5 6 7 8 9 10 ... $ chromosome.r20 : Factor w/ 20 levels "chr1","chr10",..: 18 17 19 3 17 11 15 15 7 11 ... $ start.r20 : int 62984920 86418151 119057160 105829258 126705221 34314825 90947974 90919739 72932054 9063773 ... $ end.r20 : int 63024474 86506487 119066211 105851278 126744208 34329863 90982570 90937933 73015377 9093686 ... $ gene_id.r20 : Factor w/ 212 levels "NM_001004366",..: 102 111 177 209 175 36 168 37 180 75 ... $ gc_content.r20 : num 0.399 0.474 0.558 0.526 0.425 ... $ strand.r20 : Factor w/ 2 levels "-","+": 2 2 1 2 2 1 2 2 1 2 ... $ biotype.r20 : Factor w/ 3 levels "nonsense_mediated_decay",..: 2 2 2 2 2 2 2 2 2 2 ... $ p.value_edger.r20 : num 7.35e-12 4.96e-02 1.82e-02 3.47e-05 2.00e-09 ... $ FDR_edger.r20 : num 4.38e-09 1.00 1.00 6.27e-03 8.45e-07 ... $ log2_normalized_fold_change_WT_vs_TR : num -3.185 0.836 -0.997 -1.783 -2.724 ... $ log2_normalized_mean_counts_WT.r20 : num 9.8 8.29 9.91 10.15 9.69 ... $ log2_normalized_median_counts_WT.r20 : num 9.8 8.29 9.91 10.15 9.69 ... $ log2_normalized_sd_counts_WT.r20 : logi NA NA NA NA NA NA ... $ log2_normalized_mad_counts_WT.r20 : int 0 0 0 0 0 0 0 0 0 0 ... $ log2_normalized_cv_counts_WT.r20 : logi NA NA NA NA NA NA ... $ log2_normalized_rcv_counts_WT.r20 : int 0 0 0 0 0 0 0 0 0 0 ... $ log2_normalized_mean_counts_TR : num 6.63 9.12 8.91 8.37 6.98 ... $ log2_normalized_median_counts_TR : num 6.63 9.12 8.91 8.37 6.98 ... $ log2_normalized_sd_counts_TR : logi NA NA NA NA NA NA ... $ log2_normalized_mad_counts_TR : int 0 0 0 0 0 0 0 0 0 0 ... $ log2_normalized_cv_counts_TR : logi NA NA NA NA NA NA ... $ log2_normalized_rcv_counts_TR : int 0 0 0 0 0 0 0 0 0 0 ... $ chromosome.r107 : Factor w/ 20 levels "chr1","chr10",..: NA NA NA 3 NA NA NA NA NA NA ... $ start.r107 : int NA NA NA 105829258 NA NA NA NA NA NA ... $ end.r107 : int NA NA NA 105851278 NA NA NA NA NA NA ... $ gene_id.r107 : Factor w/ 440 levels "NM_001001130",..: NA NA NA 423 NA NA NA NA NA NA ... $ gc_content.r107 : num NA NA NA 0.526 NA ... $ strand.r107 : Factor w/ 2 levels "-","+": NA NA NA 2 NA NA NA NA NA NA ... $ biotype.r107 : Factor w/ 10 levels "antisense","lincRNA",..: NA NA NA 6 NA NA NA NA NA NA ... $ p.value_edger.r107 : num NA NA NA 0.0386 NA ... $ FDR_edger.r107 : num NA NA NA 1 NA ... $ log2_normalized_fold_change_WT_vs_HET : num NA NA NA -1.48 NA ... $ log2_normalized_mean_counts_WT.r107 : num NA NA NA 8.01 NA ... $ log2_normalized_median_counts_WT.r107 : num NA NA NA 8.01 NA ... $ log2_normalized_sd_counts_WT.r107 : num NA NA NA 1.77 NA ... $ log2_normalized_mad_counts_WT.r107 : num NA NA NA 1.86 NA ... $ log2_normalized_cv_counts_WT.r107 : num NA NA NA 0.221 NA ... $ log2_normalized_rcv_counts_WT.r107 : num NA NA NA 0.232 NA ... $ log2_normalized_mean_counts_HET : num NA NA NA 7 NA ... $ log2_normalized_median_counts_HET : num NA NA NA 7 NA ... $ log2_normalized_sd_counts_HET : num NA NA NA 0.436 NA ... $ log2_normalized_mad_counts_HET : num NA NA NA 0.458 NA ... $ log2_normalized_cv_counts_HET : num NA NA NA 0.0624 NA ... $ log2_normalized_rcv_counts_HET : num NA NA NA 0.0654 NA ... $ chromosome : Factor w/ 21 levels "chr1","chr10",..: NA NA 19 NA 17 NA 15 15 NA NA ... $ start : int NA NA 119057160 NA 126705221 NA 90947974 90919739 NA NA ... $ end : int NA NA 119066211 NA 126744208 NA 90982570 90937933 NA NA ... $ gene_id : Factor w/ 414 levels "NM_001001806",..: NA NA 364 NA 362 NA 348 64 NA NA ... $ gc_content : num NA NA 0.558 NA 0.425 ... $ strand : Factor w/ 2 levels "-","+": NA NA 1 NA 2 NA 2 2 NA NA ... $ biotype : Factor w/ 7 levels "antisense","lincRNA",..: NA NA 4 NA 4 NA 4 4 NA NA ... $ p.value_edger : num NA NA 6.56e-07 NA 5.77e-08 ... $ FDR_edger : num NA NA 3.20e-04 NA 4.12e-05 ... $ log2_normalized_fold_change_WT_vs_HET.r142: num NA NA -1.79 NA -2.12 ... $ log2_normalized_mean_counts_WT : num NA NA 11.6 NA 13.7 ... $ log2_normalized_median_counts_WT : num NA NA 11.6 NA 13.7 ... $ log2_normalized_sd_counts_WT : num NA NA 0.00261 NA 0.28602 ... $ log2_normalized_mad_counts_WT : num NA NA 0.00273 NA 0.29985 ... $ log2_normalized_cv_counts_WT : num NA NA 0.000224 NA 0.02089 ... $ log2_normalized_rcv_counts_WT : num NA NA 0.000235 NA 0.021901 ... $ log2_normalized_mean_counts_HET.r142 : num NA NA 9.8 NA 11.6 ... $ log2_normalized_median_counts_HET.r142 : num NA NA 9.8 NA 11.6 ... $ log2_normalized_sd_counts_HET.r142 : num NA NA 0.358 NA 0.191 ... $ log2_normalized_mad_counts_HET.r142 : num NA NA 0.375 NA 0.2 ... $ log2_normalized_cv_counts_HET.r142 : num NA NA 0.0365 NA 0.0165 ... $ log2_normalized_rcv_counts_HET.r142 : num NA NA 0.0383 NA 0.0173 ... $ chromosome.r128 : Factor w/ 20 levels "chr1","chr10",..: 18 NA 19 3 17 NA 15 15 NA NA ... $ start.r128 : int 62984920 NA 119057160 105829258 126705221 NA 90947974 90919739 NA NA ... $ end.r128 : int 63024474 NA 119066211 105851278 126744208 NA 90982570 90937933 NA NA ... $ gene_id.r128 : Factor w/ 1703 levels "NM_001001144",..: 677 NA 1424 1656 1418 NA 1360 281 NA NA ... $ gc_content.r128 : num 0.399 NA 0.558 0.526 0.425 ... $ strand.r128 : Factor w/ 2 levels "-","+": 2 NA 1 2 2 NA 2 2 NA NA ... $ biotype.r128 : Factor w/ 14 levels "ambiguous_orf",..: 8 NA 8 8 8 NA 8 8 NA NA ... $ p.value_edger.r128 : num 4.21e-30 NA 7.50e-09 1.73e-25 1.09e-33 ... $ FDR_edger.r128 : num 2.20e-27 NA 6.54e-07 5.66e-23 7.14e-31 ... $ log2_normalized_fold_change_WT_vs_MUT : num -4.12 NA -1.61 -2.77 -3.13 ... $ log2_normalized_mean_counts_WT.r128 : num 10.67 NA 9.63 11.09 10.52 ... $ log2_normalized_median_counts_WT.r128 : num 10.67 NA 9.63 11.09 10.52 ... $ log2_normalized_sd_counts_WT.r128 : num 0.351 NA 0.135 0.174 0.327 ... $ log2_normalized_mad_counts_WT.r128 : num 0.368 NA 0.142 0.182 0.342 ... $ log2_normalized_cv_counts_WT.r128 : num 0.0329 NA 0.0141 0.0156 0.031 ... $ log2_normalized_rcv_counts_WT.r128 : num 0.0344 NA 0.0147 0.0164 0.0325 ... $ log2_normalized_mean_counts_MUT : num 6.44 NA 8.02 8.26 7.39 ... $ log2_normalized_median_counts_MUT : num 6.44 NA 8.02 8.26 7.39 ... $ log2_normalized_sd_counts_MUT : num 0.931 NA 0.0471 0.6224 0.374 ... $ log2_normalized_mad_counts_MUT : num 0.976 NA 0.0494 0.6525 0.3921 ... $ log2_normalized_cv_counts_MUT : num 0.14458 NA 0.00587 0.07535 0.05058 ... $ log2_normalized_rcv_counts_MUT : num 0.15157 NA 0.00615 0.079 0.05303 ... y=cbind(x$FDR_edger.r20,x$log2_normalized_fold_change_WT_vs_TR,x$log2_normalized_mean_counts_WT.r20,x$log2_normalized_mean_counts_TR,x$log2_normalized_fold_change_WT_vs_HET,x$log2_normalized_mean_counts_WT.r107,x$log2_normalized_mean_counts_HET,x$log2_normalized_fold_change_WT_vs_HET.r142,x$log2_normalized_mean_counts_WT,x$log2_normalized_mean_counts_HET.r142,x$log2_normalized_fold_change_WT_vs_MUT,x$log2_normalized_mean_counts_WT.r128,x$log2_normalized_mean_counts_MUT) # diff.exp correlation png("diff_exp_e15.5_VS_diff_exp_e14.5.png") plot(x$log2_normalized_fold_change_WT_vs_TR,x$log2_normalized_fold_change_WT_vs_HET,xlab="WT_vs_HET_e15.5",ylab="WT_vs_HET_e14.5") dev.off() cor.test(x$log2_normalized_fold_change_WT_vs_TR,x$log2_normalized_fold_change_WT_vs_HET) Pearson's product-moment correlation data: x$log2_normalized_fold_change_WT_vs_TR and x$log2_normalized_fold_change_WT_vs_HET t = 3.5537, df = 31, p-value = 0.00124 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.2388222 0.7439209 sample estimates: cor 0.5380191 png("diff_exp_e15.5_VS_diff_exp_RCAD_P1.png") plot(x$log2_normalized_fold_change_WT_vs_TR,x$log2_normalized_fold_change_WT_vs_HET.r142,xlab="WT_vs_HET_e15.5",ylab="WT_vs_HET_RCAD_P1") dev.off() cor.test(x$log2_normalized_fold_change_WT_vs_TR,x$log2_normalized_fold_change_WT_vs_HET.r142) Pearson's product-moment correlation data: x$log2_normalized_fold_change_WT_vs_TR and x$log2_normalized_fold_change_WT_vs_HET.r142 t = 3.6153, df = 76, p-value = 0.000537 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.1754955 0.5580284 sample estimates: cor 0.3830663 png("diff_exp_e15.5_VS_diff_exp_HoxB7cre.png") plot(x$log2_normalized_fold_change_WT_vs_TR,x$log2_normalized_fold_change_WT_vs_MUT,xlab="WT_vs_HET_e15.5",ylab="WT_vs_HET_HoxB7cre") dev.off() cor.test(x$log2_normalized_fold_change_WT_vs_TR,x$log2_normalized_fold_change_WT_vs_MUT) Pearson's product-moment correlation data: x$log2_normalized_fold_change_WT_vs_TR and x$log2_normalized_fold_change_WT_vs_MUT t = 28.558, df = 145, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.8927139 0.9427031 sample estimates: cor 0.9214371 png("diff_exp_e14.5_VS_diff_exp_RCAD_P1.png") plot(x$log2_normalized_fold_change_WT_vs_HET,x$log2_normalized_fold_change_WT_vs_HET.r142,xlab="WT_vs_HET_e14.5",ylab="WT_vs_HET_RCAD_P1") dev.off() cor.test(x$log2_normalized_fold_change_WT_vs_HET,x$log2_normalized_fold_change_WT_vs_HET.r142) Pearson's product-moment correlation data: x$log2_normalized_fold_change_WT_vs_HET and x$log2_normalized_fold_change_WT_vs_HET.r142 t = -0.0977, df = 26, p-value = 0.9229 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.3894473 0.3564709 sample estimates: cor -0.0191532 png("diff_exp_e14.5_VS_diff_exp_HoxB7cre.png") plot(x$log2_normalized_fold_change_WT_vs_HET,x$log2_normalized_fold_change_WT_vs_MUT,xlab="WT_vs_HET_e14.5",ylab="WT_vs_HET_HoxB7cre") dev.off() cor.test(x$log2_normalized_fold_change_WT_vs_HET,x$log2_normalized_fold_change_WT_vs_MUT) Pearson's product-moment correlation data: x$log2_normalized_fold_change_WT_vs_HET and x$log2_normalized_fold_change_WT_vs_MUT t = 4.677, df = 118, p-value = 7.812e-06 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.2327083 0.5366569 sample estimates: cor 0.3954548 png("diff_exp_HoxB7cre_VS_diff_exp_RCAD_P1.png") plot(x$log2_normalized_fold_change_WT_vs_MUT,x$log2_normalized_fold_change_WT_vs_HET.r142,xlab="WT_vs_HET_HoxB7cre",ylab="WT_vs_HET_RCAD_P1") dev.off() cor.test(x$log2_normalized_fold_change_WT_vs_MUT,x$log2_normalized_fold_change_WT_vs_HET.r142) Pearson's product-moment correlation data: x$log2_normalized_fold_change_WT_vs_MUT and x$log2_normalized_fold_change_WT_vs_HET.r142 t = 9.6145, df = 166, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.4911295 0.6872438 sample estimates: cor 0.5980661 # WT.exp correlation png("WT_e15.5_VS_WT_e14.5.png") plot(x$log2_normalized_mean_counts_WT.r20,x$log2_normalized_mean_counts_WT.r107,xlab="WT_e15.5",ylab="WT_e14.5") dev.off() cor.test(x$log2_normalized_mean_counts_WT.r20,x$log2_normalized_mean_counts_WT.r107) Pearson's product-moment correlation data: x$log2_normalized_mean_counts_WT.r20 and x$log2_normalized_mean_counts_WT.r107 t = 3.8234, df = 31, p-value = 0.0005952 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.2765152 0.7614188 sample estimates: cor 0.5660856 png("WT_e15.5_VS_WT_r142.png") plot(x$log2_normalized_mean_counts_WT.r20,x$log2_normalized_mean_counts_WT,xlab="WT_e15.5",ylab="WT_r142") dev.off() cor.test(x$log2_normalized_mean_counts_WT.r20,x$log2_normalized_mean_counts_WT) Pearson's product-moment correlation data: x$log2_normalized_mean_counts_WT.r20 and x$log2_normalized_mean_counts_WT t = 6.5386, df = 76, p-value = 6.41e-09 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.4356562 0.7256566 sample estimates: cor 0.6000171 png("WT_e15.5_VS_WT_r128.png") plot(x$log2_normalized_mean_counts_WT.r20,x$log2_normalized_mean_counts_WT.r128,xlab="WT_e15.5",ylab="WT_r128") dev.off() cor.test(x$log2_normalized_mean_counts_WT.r20,x$log2_normalized_mean_counts_WT.r128) Pearson's product-moment correlation data: x$log2_normalized_mean_counts_WT.r20 and x$log2_normalized_mean_counts_WT.r128 t = 21.5004, df = 145, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.8274677 0.9063538 sample estimates: cor 0.872482 png("WT_e14.5_VS_WT_r142.png") plot(x$log2_normalized_mean_counts_WT.r107,x$log2_normalized_mean_counts_WT,xlab="WT_e14.5",ylab="WT_r142") dev.off() cor.test(x$log2_normalized_mean_counts_WT.r107,x$log2_normalized_mean_counts_WT) Pearson's product-moment correlation data: x$log2_normalized_mean_counts_WT.r107 and x$log2_normalized_mean_counts_WT t = 1.0602, df = 26, p-value = 0.2988 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.1834392 0.5359430 sample estimates: cor 0.2035694 png("WT_e14.5_VS_WT_r128.png") plot(x$log2_normalized_mean_counts_WT.r107,x$log2_normalized_mean_counts_WT.r128,xlab="WT_e14.5",ylab="WT_r128") dev.off() cor.test(x$log2_normalized_mean_counts_WT.r107,x$log2_normalized_mean_counts_WT.r128) Pearson's product-moment correlation data: x$log2_normalized_mean_counts_WT.r107 and x$log2_normalized_mean_counts_WT.r128 t = 9.6881, df = 118, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.5522460 0.7547943 sample estimates: cor 0.6656025 png("WT_r128_VS_WT_r142.png") plot(x$log2_normalized_mean_counts_WT.r128,x$log2_normalized_mean_counts_WT,xlab="WT_e14.5",ylab="WT_r142") dev.off() cor.test(x$log2_normalized_mean_counts_WT.r128,x$log2_normalized_mean_counts_WT) Pearson's product-moment correlation data: x$log2_normalized_mean_counts_WT.r128 and x$log2_normalized_mean_counts_WT t = 11.1059, df = 166, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.5565055 0.7319536 sample estimates: cor 0.6529015 # all genes: ar107=read.table("/data/images/proton/Cereghini/run142/metaseqr_wt_vs_het_run107/lists/metaseqr_all_out_WT_vs_HET.txt.gz",header=TRUE) ar128=read.table("/data/images/proton/Cereghini/run142/metaseqr_wt_vs_mut_run128/lists/metaseqr_all_out_WT_vs_MUT.txt.gz",header=TRUE) ar142=read.table("/data/images/proton/Cereghini/run142/metaseqr_wt_vs_mut_run142/lists/metaseqr_all_out_WT_vs_HET.txt.gz",header=TRUE) ar20=read.table("/data/images/proton/Cereghini/run142/metaseqr_wt_vs_tr_run20/lists/metaseqr_all_out_WT_vs_TR.txt.gz",header=TRUE) genes=levels(rbind(ar20[,1:10],ar107[,1:10],ar128[,1:10],ar142[,1:10])) x=merge(ar20[,1:23], ar107[,1:23], by='gene_name', all=T,suffixes = c(".ar20",".ar107")) y=merge(x, ar142[,1:23], by='gene_name', all=T,suffixes = c("",".ar142")) x=merge(y, ar128[,1:23], by='gene_name', all=T,suffixes = c("",".ar128")) write.table(x,"foo") png("WT_e15.5_VS_WT_e14.5.png") plot(x$log2_normalized_mean_counts_WT.ar20,x$log2_normalized_mean_counts_WT.ar107,xlab="WT_e15.5",ylab="WT_e14.5") dev.off() cor.test(x$log2_normalized_mean_counts_WT.ar20,x$log2_normalized_mean_counts_WT.ar107) png("WT_e15.5_VS_WT_RCAD_P1.png") plot(x$log2_normalized_mean_counts_WT.ar20,x$log2_normalized_mean_counts_WT,xlab="WT_e15.5",ylab="WT_RCAD_P1") dev.off() cor.test(x$log2_normalized_mean_counts_WT.ar20,x$log2_normalized_mean_counts_WT) png("WT_e15.5_VS_WT_HoxB7cre.png") plot(x$log2_normalized_mean_counts_WT.ar20,x$log2_normalized_mean_counts_WT.ar128,xlab="WT_e15.5",ylab="WT_HoxB7cre") dev.off() cor.test(x$log2_normalized_mean_counts_WT.ar20,x$log2_normalized_mean_counts_WT.ar128) png("WT_e14.5_VS_WT_RCAD_P1.png") plot(x$log2_normalized_mean_counts_WT.ar107,x$log2_normalized_mean_counts_WT,xlab="WT_e14.5",ylab="WT_RCAD_P1") dev.off() cor.test(x$log2_normalized_mean_counts_WT.ar107,x$log2_normalized_mean_counts_WT) png("WT_e14.5_VS_WT_HoxB7cre.png") plot(x$log2_normalized_mean_counts_WT.ar107,x$log2_normalized_mean_counts_WT.ar128,xlab="WT_e14.5",ylab="WT_HoxB7cre") dev.off() cor.test(x$log2_normalized_mean_counts_WT.ar107,x$log2_normalized_mean_counts_WT.ar128) png("WT_HoxB7cre_VS_WT_RCAD_P1.png") plot(x$log2_normalized_mean_counts_WT.ar128,x$log2_normalized_mean_counts_WT,xlab="WT_HoxB7cre",ylab="WT_RCAD_P1") dev.off() cor.test(x$log2_normalized_mean_counts_WT.ar128,x$log2_normalized_mean_counts_WT) #@ rerun 09062017 library(metaseqR) the.path="/data/images/proton/Cereghini/run142" metaseqr( sample.list=file.path(the.path,"targets_run142_rerun09062017.txt"), contrast=c("WT_vs_HET"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_wt_vs_het_run142_rerun09062017"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.3, fig.format=c("png","pdf") )