bvieth / powsimR

Power analysis is essential to optimize the design of RNA-seq experiments and to assess and compare the power to detect differentially expressed genes. PowsimR is a flexible tool to simulate and evaluate differential expression from bulk and especially single-cell RNA-seq data making it suitable for a priori and posterior power analyses.
https://bvieth.github.io/powsimR/
Artistic License 2.0
103 stars 23 forks source link

estimateParam() fails with some but not all combinations of input dataframe columns #68

Open BenjaminATaylor opened 1 week ago

BenjaminATaylor commented 1 week ago

My input counts frame, counts.true.tmp, has dimensions (12319,30). Running estimateParam() using this frame as the input countData fails with "Error in ParamData$means : $ operator is invalid for atomic vectors". However, if I take a subset of the columns (anything from 3 to 19 columns), the function runs as expected. The problem in this case appears not to be tied to a specific column, i.e. it occurs when I use counts.true.tmp[,c(1:20)] but not if I use counts.true.tmp[,c(1:19)] or counts.true.tmp[,c(2:20)]. Example code:

estimateParam(countData = counts.true.tmp[c(1:20)],
                               readData = NULL,
                               batchData = NULL,
                               spikeData = NULL,
                               spikeInfo = NULL,
                               Lengths = NULL, 
                               MeanFragLengths = NULL,
                               RNAseq = 'bulk', Protocol = 'Read',
                               Distribution = 'NB', Normalisation = "MR",
                               GeneFilter = 0.1, SampleFilter = 3, 
                               sigma = 1.96, NCores = NULL, verbose = TRUE)

Because of the above, I initially thought this problem was tied to column number. However, when I used a subset of columns to produce a reproducible example for this issue, I found a different situation: counts.true.tmp[c(1:50),c(1:15,17:30)] works, but counts.true.tmp[c(1:50),c(1:30)] does not, indicating that the problem in this case is specifically with column 16. However, I cannot see anything obviously different about this column, and actually this column gives no problem when we use a smaller subset still, e.g. counts.true.tmp[c(1:50),c(10:16)]. So I'm very confused!

Here's the full counts frame, and session info below.

R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.4.1

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

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

time zone: Europe/London
tzcode source: internal

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

other attached packages:
 [1] powsimR_1.2.5               gamlss.dist_6.1-1           devtools_2.4.5              usethis_2.2.3              
 [5] DescTools_0.99.54           compcodeR_1.38.0            sm_2.2-6.0                  DESeq2_1.42.1              
 [9] SummarizedExperiment_1.32.0 Biobase_2.62.0              MatrixGenerics_1.14.0       matrixStats_1.3.0          
[13] GenomicRanges_1.54.1        GenomeInfoDb_1.38.8         IRanges_2.36.0              S4Vectors_0.40.2           
[17] BiocGenerics_0.48.1         lubridate_1.9.3             forcats_1.0.0               stringr_1.5.1              
[21] dplyr_1.1.4                 purrr_1.0.2                 readr_2.1.5                 tidyr_1.3.1                
[25] tibble_3.2.1                ggplot2_3.5.1               tidyverse_2.0.0            

loaded via a namespace (and not attached):
  [1] igraph_2.0.3                plotly_4.10.4               Formula_1.2-5              
  [4] scater_1.30.1               maps_3.4.2                  zlibbioc_1.48.2            
  [7] fpc_2.2-12                  bit_4.0.5                   tidyselect_1.2.1           
 [10] doParallel_1.0.17           clue_0.3-65                 lattice_0.22-6             
 [13] rjson_0.2.21                blob_1.2.4                  urlchecker_1.0.1           
 [16] S4Arrays_1.2.1              parallel_4.3.2              png_0.1-8                  
 [19] plotrix_3.8-4               cli_3.6.3                   msir_1.3.3                 
 [22] outliers_0.15               DECENT_1.1.1                lavaan_0.6-18              
 [25] BiocIO_1.12.0               doSNOW_1.0.20               lpsymphony_1.30.0          
 [28] kernlab_0.9-32              bluster_1.12.0              BiocNeighbors_1.20.2       
 [31] curl_5.2.1                  zinbwave_1.24.0             mime_0.12                  
 [34] evaluate_0.24.0             stringi_1.8.4               nonnest2_0.5-7             
 [37] backports_1.5.0             fBasics_4032.96             desc_1.4.3                 
 [40] XML_3.99-0.16.1             Exact_3.2                   httpuv_1.6.15              
 [43] flexmix_2.3-19              AnnotationDbi_1.64.1        magrittr_2.0.3             
 [46] rappdirs_0.3.3              splines_4.3.2               mclust_6.1.1               
 [49] jpeg_0.1-10                 IHW_1.30.0                  DT_0.33                    
 [52] sctransform_0.4.1           rootSolve_1.8.2.4           ggbeeswarm_0.7.2           
 [55] sessioninfo_1.2.2           statip_0.2.3                DBI_1.2.3                  
 [58] genefilter_1.84.0           withr_3.0.0                 class_7.3-22               
 [61] rprojroot_2.0.4             brio_1.1.5                  rtracklayer_1.62.0         
 [64] compositions_2.0-8          htmlwidgets_1.6.4           fs_1.6.4                   
 [67] segmented_2.1-0             biomaRt_2.58.2              SingleCellExperiment_1.24.0
 [70] ggrepel_0.9.5               SCnorm_1.24.0               baySeq_2.36.2              
 [73] SparseArray_1.2.4           cellranger_1.1.0            DEoptimR_1.1-3             
 [76] mixtools_2.0.0              spatial_7.3-17              diptest_0.77-1             
 [79] annotate_1.80.0             truncnorm_1.0-9             lmom_3.0                   
 [82] reticulate_1.37.0           minpack.lm_1.2-4            hexbin_1.28.3              
 [85] zoo_1.8-12                  XVector_0.42.0              knitr_1.47                 
 [88] gmodels_2.19.1              modeest_2.4.0               timechange_0.3.0           
 [91] pbivnorm_0.6.0              foreach_1.5.2               fansi_1.0.6                
 [94] caTools_1.18.2              grid_4.3.2                  rhdf5_2.46.1               
 [97] data.table_1.15.4           timeDate_4032.109           vegan_2.6-6.1              
[100] R.oo_1.26.0                 quantreg_5.98               RSpectra_0.16-1            
[103] irlba_2.3.5.1               ellipsis_0.3.2              lazyeval_0.2.2             
[106] yaml_2.3.8                  aroma.light_3.32.0          survival_3.7-0             
[109] crayon_1.5.3                tensorA_0.36.2.1            RColorBrewer_1.1-3         
[112] later_1.3.2                 codetools_0.2-20            base64enc_0.1-3            
[115] profvis_0.3.8               KEGGREST_1.42.0             HSMMSingleCell_1.22.0      
[118] Rtsne_0.17                  fastICA_1.2-4               gdata_3.0.0                
[121] limma_3.58.1                filelock_1.0.3              Rsamtools_2.18.0           
[124] shinyBS_0.61.1              DrImpute_1.0                foreign_0.8-86             
[127] pkgconfig_2.0.3             NBPSeq_0.3.1                xml2_1.3.6                 
[130] ggpubr_0.6.0                GenomicAlignments_1.38.2    rmutil_1.1.10              
[133] viridisLite_0.4.2           xtable_1.8-4                interp_1.1-6               
[136] hwriter_1.3.2.1             car_3.1-2                   SAVER_1.1.3                
[139] fastcluster_1.2.6           plyr_1.8.9                  httr_1.4.7                 
[142] tools_4.3.2                 globals_0.16.3              prabclus_2.3-3             
[145] pkgbuild_1.4.4              beeswarm_0.4.0              htmlTable_2.4.2            
[148] broom_1.0.6                 checkmate_2.3.1             nlme_3.1-165               
[151] EBSeq_2.0.0                 dbplyr_2.5.0                MatrixModels_0.5-3         
[154] lme4_1.1-35.4               digest_0.6.36               bayesm_3.1-6               
[157] permute_0.9-7               Matrix_1.6-5                tzdb_0.4.0                 
[160] reshape2_1.4.4              viridis_0.6.5               rpart_4.1.23               
[163] BiocFileCache_2.10.2        glue_1.7.0                  cachem_1.1.0               
[166] UpSetR_1.4.0                Biostrings_2.70.3           Hmisc_5.1-3                
[169] generics_0.1.3              CompQuadForm_1.4.3          mvtnorm_1.2-5              
[172] amap_0.8-19                 parallelly_1.37.1           pkgload_1.3.4              
[175] mnormt_2.1.1                statmod_1.5.0               pscl_1.5.9                 
[178] arm_1.14-4                  here_1.0.1                  Linnorm_2.26.0             
[181] blockmodeling_1.1.5         ScaledMatrix_1.10.0         carData_3.0-5              
[184] minqa_1.2.7                 timeSeries_4032.109         fields_15.2                
[187] BPSC_0.99.2                 spam_2.10-0                 bayNorm_1.20.0             
[190] ggstance_0.3.7              dqrng_0.4.1                 snow_0.4-4                 
[193] utf8_1.2.4                  gtools_3.9.5                readxl_1.4.3               
[196] softImpute_1.4-1            ggsignif_0.6.4              RcppEigen_0.3.4.0.0        
[199] gridExtra_2.3               shiny_1.8.1.1               GenomeInfoDbData_1.2.11    
[202] apcluster_1.4.13            rhdf5filters_1.14.1         R.utils_2.12.3             
[205] RCurl_1.98-1.14             memoise_2.0.1               rmarkdown_2.27             
[208] pheatmap_1.0.12             R.methodsS3_1.8.2           scales_1.3.0               
[211] Rmagic_2.0.3                stabledist_0.7-1            gld_2.6.6                  
[214] future_1.33.2               RANN_2.6.1                  rstudioapi_0.16.0          
[217] cluster_2.1.6               cobs_1.3-8                  hms_1.1.3                  
[220] fitdistrplus_1.1-11         munsell_0.5.1               vioplot_0.4.0              
[223] fdrtool_1.2.17              cowplot_1.1.3               colorspace_2.1-0           
[226] ZIM_1.1.0                   ellipse_0.5.0               rlang_1.1.4                
[229] ggdendro_0.2.0              quadprog_1.5-8              iCOBRA_1.30.0              
[232] DelayedMatrixStats_1.24.0   sparseMatrixStats_1.14.0    dotCall64_1.1-1            
[235] shinydashboard_0.7.2        DDRTree_0.1.5               scuttle_1.12.0             
[238] mgcv_1.9-1                  xfun_0.45                   NOISeq_2.46.0              
[241] coda_0.19-4.1               e1071_1.7-14                rARPACK_0.11-0             
[244] remotes_2.5.0               iterators_1.0.14            modeltools_0.2-23          
[247] abind_1.4-5                 Rhdf5lib_1.24.2             bitops_1.0-7               
[250] ps_1.7.6                    promises_1.3.0              RSQLite_2.3.7              
[253] leidenbase_0.1.27           qvalue_2.34.0               sandwich_3.1-0             
[256] DelayedArray_0.28.0         proxy_0.4-27                compiler_4.3.2             
[259] prettyunits_1.2.0           boot_1.3-30                 beachmat_2.18.1            
[262] SparseM_1.83                listenv_0.9.1               Rcpp_1.0.12                
[265] edgeR_4.0.16                BiocSingular_1.18.0         progress_1.2.3             
[268] ROTS_1.30.0                 MASS_7.3-60.0.1             BiocParallel_1.36.0        
[271] MAST_1.28.0                 stable_1.1.6                monocle_2.30.1             
[274] R6_2.5.1                    RcppArmadillo_0.12.8.4.0    fastmap_1.2.0              
[277] rstatix_0.7.2               EDASeq_2.36.0               vipor_0.4.7                
[280] ROCR_1.0-11                 rsvd_1.0.5                  nnet_7.3-19                
[283] gtable_0.3.5                KernSmooth_2.23-24          latticeExtra_0.6-30        
[286] miniUI_0.1.1.1              zingeR_0.1.0                deldir_2.0-4               
[289] htmltools_0.5.8.1           bit64_4.0.5                 scDD_1.26.0                
[292] penalized_0.9-52            lifecycle_1.0.4             scone_1.26.0               
[295] processx_3.8.4              RUVSeq_1.36.0               restfulr_0.0.15            
[298] nloptr_2.1.0                callr_3.7.6                 vctrs_0.6.5                
[301] testthat_3.2.1.1            VGAM_1.1-11                 slam_0.1-50                
[304] robustbase_0.99-2           scran_1.30.2                future.apply_1.11.2        
[307] pillar_1.9.0                GenomicFeatures_1.54.4      moments_0.14.1             
[310] gplots_3.1.3.1              ShortRead_1.60.0            metapod_1.10.1             
[313] locfit_1.5-9.10             combinat_0.0-8              jsonlite_1.8.8             
[316] expm_0.999-9                markdown_1.13