Real Data with Simulated Signals: Part I

Introduction

Using David’s package seqgendiff, we are adding artefactual signals to the real GTEx Liver RNA-seq data.

mat = read.csv("../data/liver.csv")

The true signal comes from a mixture distribution

\[ g\left(\beta\right) = \pi_0\delta_0 + \left(1 - \pi_0\right)N\left(0, \sigma^2\right) \] The simulated data matrices are then fed into edgeR, limma pipeline. In the following simulations, we always use \(5\) vs \(5\).

Case 1: \(\pi_0 = 0.9\), \(\sigma^2 = 1\).

N = 100
nsamp = 10
pi0 = 0.9
sd = 1
system.time(ashvgdash <- N_simulations(N, mat, nsamp, pi0, sd))

     user    system   elapsed 
 6854.603   627.845 12855.081 

Case 2: \(\pi_0 = 0.9\), \(\sigma^2 = 4\)

N = 100
nsamp = 10
pi0 = 0.9
sd = 2
system.time(ashvgdash <- N_simulations(N, mat, nsamp, pi0, sd))

    user   system  elapsed 
5877.223  621.433 5256.693 

Case 3: \(\pi_0 = 0.5\), \(\sigma^2 = 4\)

N = 100
nsamp = 10
pi0 = 0.5
sd = 2
system.time(ashvgdash <- N_simulations(N, mat, nsamp, pi0, sd))

     user    system   elapsed 
 5269.488   574.588 15042.479 

Case 4: \(\pi_0 = 0.9\), \(\sigma^2 = 9\)

N = 100
nsamp = 10
pi0 = 0.9
sd = 3
system.time(ashvgdash <- N_simulations(N, mat, nsamp, pi0, sd))

     user    system   elapsed 
 5321.625   558.205 15356.345 

Session information

sessionInfo()

R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] Rmosek_7.1.3        PolynomF_0.94       cvxr_0.0.0.9400    
 [4] REBayes_0.85        Matrix_1.2-11       SQUAREM_2017.10-1  
 [7] EQL_1.0-0           ttutils_1.0-1       cate_1.0.4         
[10] sva_3.26.0          BiocParallel_1.12.0 genefilter_1.60.0  
[13] mgcv_1.8-22         nlme_3.1-131        seqgendiff_0.1.0   
[16] qvalue_2.10.0       edgeR_3.20.1        limma_3.34.0       
[19] ashr_2.1-27        

loaded via a namespace (and not attached):
 [1] Biobase_2.38.0       svd_0.4.1            bit64_0.9-7         
 [4] splines_3.4.2        foreach_1.4.3        stats4_3.4.2        
 [7] blob_1.1.0           yaml_2.1.14          RSQLite_2.0         
[10] backports_1.1.1      lattice_0.20-35      digest_0.6.12       
[13] colorspace_1.3-2     htmltools_0.3.6      plyr_1.8.4          
[16] XML_3.98-1.9         esaBcv_1.2.1         xtable_1.8-2        
[19] corpcor_1.6.9        scales_0.5.0         git2r_0.19.0        
[22] tibble_1.3.4         annotate_1.56.0      gmp_0.5-13.1        
[25] IRanges_2.12.0       ggplot2_2.2.1        BiocGenerics_0.24.0 
[28] lazyeval_0.2.1       survival_2.41-3      magrittr_1.5        
[31] memoise_1.1.0        evaluate_0.10.1      doParallel_1.0.11   
[34] MASS_7.3-47          truncnorm_1.0-7      tools_3.4.2         
[37] matrixStats_0.52.2   stringr_1.2.0        S4Vectors_0.16.0    
[40] munsell_0.4.3        locfit_1.5-9.1       AnnotationDbi_1.40.0
[43] compiler_3.4.2       rlang_0.1.4          grid_3.4.2          
[46] leapp_1.2            RCurl_1.95-4.8       iterators_1.0.8     
[49] bitops_1.0-6         rmarkdown_1.6        gtable_0.2.0        
[52] codetools_0.2-15     DBI_0.7              reshape2_1.4.2      
[55] ruv_0.9.6            knitr_1.17           bit_1.1-12          
[58] rprojroot_1.2        stringi_1.1.5        pscl_1.5.2          
[61] parallel_3.4.2       Rcpp_0.12.13