I am trying to perform power calculation using the powsimR. I ran the analysis as specified in the vignette but I keep encountering an error after running the simulateDE function:
#Estimation
estparam_gene <- estimateParam(countData = as.data.frame(GetAssayData (astrocyte_SCT_counts, assay = "SCT", slot = "count")),
readData = NULL,
batchData = NULL,
spikeData = NULL,
spikeInfo = NULL,
Lengths = NULL, MeanFragLengths = NULL,
RNAseq = 'singlecell', Protocol = 'UMI',
Distribution = 'ZINB', Normalisation = "sctransform",
GeneFilter = 0.05, SampleFilter = 5,
sigma = 1.96, NCores = NULL, verbose = TRUE)
# define log fold change
p.lfc <- function(x) sample(c(-1,1), size=x,replace=T)*rgamma(x, shape = 1, rate = 2)
# set up simulations
setupres <- Setup(ngenes = 10000, nsims = 25,
p.DE = 0.05, pLFC = p.lfc,
n1 = c(50, 100, 200, 400, 800), n2 = c(50, 100, 200, 400, 800),
Thinning = NULL, LibSize = 'equal',
estParamRes = estparam_gene,
estSpikeRes = NULL,
DropGenes = TRUE,
setup.seed = 5299, verbose = TRUE)
#Simulation
simres <- simulateDE(SetupRes = setupres,
Prefilter = "FreqFilter", Imputation = NULL,
Normalisation = 'sctransform', Label = 'none',
DEmethod = "DESeq2-ZINB-WaVE", DEFilter = FALSE,
NCores = NULL, verbose = TRUE)
Error in DESeq2::DESeqDataSetFromMatrix(countData, coldat, ~design, tidy = FALSE, :
ncol(countData) == nrow(colData) is not TRUE
(1) Although I have input the countData when I ran the estimateParam function at the start, I am unsure where to input the colData so that ncol(countData) == nrow(colData) would be TRUE. Or is there something I have not done correctly in the workflow?
(2) On a related note, is there way to include a model design since the error mentions "~ design"?
Thanks very much for your help.
R version 4.1.1 (2021-08-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.6 LTS
Matrix products: default
BLAS/LAPACK: /opt/OpenBLAS/lib/libopenblasp-r0.3.3.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
[4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=en_US.UTF-8 LC_ADDRESS=en_US.UTF-8
[10] LC_TELEPHONE=en_US.UTF-8 LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=en_US.UTF-8
attached base packages:
[1] grid stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] powsimR_1.2.5 gamlss.dist_6.0-3 MASS_7.3-54 POWSC_1.10.0
[5] MAST_1.20.0 RnaSeqSampleSize_2.12.0 RnaSeqSampleSizeData_1.34.0 edgeR_3.36.0
[9] Hmisc_4.7-1 Formula_1.2-4 survival_3.2-11 lattice_0.20-44
[13] sva_3.42.0 genefilter_1.76.0 mgcv_1.8-36 nlme_3.1-152
[17] purrr_1.0.1 muscat_1.8.2 reticulate_1.27 batchtools_0.9.15
[21] zinbwave_1.16.0 EnhancedVolcano_1.12.0 broom_1.0.0 tibble_3.2.1
[25] clusterProfiler_4.2.2 pathfindR_1.6.3 pathfindR.data_1.1.2 DEsingle_1.14.0
[29] vsn_3.62.0 variancePartition_1.24.1 BiocParallel_1.28.3 limma_3.50.3
[33] ggpubr_0.4.0 PCAtools_2.6.0 ggrepel_0.9.1 DESeq2_1.34.0
[37] Matrix.utils_0.9.8 scDblFinder_1.8.0 pheatmap_1.0.12 scRNAseq_2.8.0
[41] SingleR_1.8.1 biomaRt_2.50.3 org.Mm.eg.db_3.14.0 cowplot_1.1.1
[45] gridExtra_2.3 kableExtra_1.3.4 scran_1.22.1 mvoutlier_2.1.1
[49] sgeostat_1.0-27 RColorBrewer_1.1-3 ensembldb_2.18.4 AnnotationFilter_1.18.0
[53] GenomicFeatures_1.46.5 AnnotationDbi_1.56.2 AnnotationHub_3.2.2 BiocFileCache_2.2.1
[57] dbplyr_2.3.3 tidyr_1.3.0 stringr_1.5.0 lubridate_1.9.0
[61] timechange_0.1.1 Matrix_1.4-1 scater_1.22.0 ggplot2_3.4.4
[65] scuttle_1.4.0 SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0 Biobase_2.54.0
[69] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1 IRanges_2.28.0 S4Vectors_0.32.4
[73] BiocGenerics_0.40.0 MatrixGenerics_1.6.0 matrixStats_0.62.0 sp_1.5-0
[77] SeuratObject_4.1.0 Seurat_4.1.1 dendextend_1.16.0
loaded via a namespace (and not attached):
[1] abind_1.4-5 SAVER_1.1.3 ROTS_1.22.0
[4] fastICA_1.2-3 GO.db_3.14.0 dqrng_0.3.0
[7] affy_1.72.0 RANN_2.6.1 ellipsis_0.3.2
[10] memoise_2.0.1 emmeans_1.8.0 munsell_0.5.0
[13] scales_1.2.1 grr_0.9.5 locfit_1.5-9.6
[16] ggridges_0.5.3 segmented_1.6-0 xlsx_0.6.5
[19] listenv_0.9.0 blob_1.2.4 interactiveDisplayBase_1.32.0
[22] GOSemSim_2.20.0 fgsea_1.20.0 quantreg_5.93
[25] nloptr_2.0.3 truncnorm_1.0-9 DOSE_3.20.1
[28] DECENT_1.1.1 lme4_1.1-30 sparseMatrixStats_1.6.0
[31] mclust_5.4.10 jpeg_0.1-9 iCOBRA_1.30.0
[34] htmltools_0.5.5 vipor_0.4.5 doSNOW_1.0.20
[37] Rsamtools_2.10.0 progressr_0.13.0 ggplotify_0.1.0
[40] msir_1.3.3 fastmatch_1.1-3 leiden_0.4.2
[43] fastcluster_1.2.3 RUVSeq_1.28.0 XML_3.99-0.10
[46] foreach_1.5.2 aod_1.3.2 GEOquery_2.70.0
[49] enrichplot_1.14.2 NBPSeq_0.3.1 crayon_1.5.2
[52] maps_3.4.0 DelayedArray_0.20.0 tidyselect_1.2.0
[55] filelock_1.0.2 lpsymphony_1.30.0 withr_2.5.0
[58] amap_0.8-19 qlcMatrix_0.9.7 densityClust_0.3.2
[61] future.apply_1.11.0 DEoptimR_1.0-11 moments_0.14.1
[64] shiny_1.7.4.1 snow_0.4-4 ape_5.6-2
[67] BiocManager_1.30.20 tzdb_0.4.0 lmerTest_3.1-3
[70] RcppArmadillo_0.12.4.1.0 NOISeq_2.46.0 Rdpack_2.3.1
[73] SparseM_1.81 multcomp_1.4-20 graphlayouts_0.8.0
[76] SCnorm_1.24.0 pbkrtest_0.5.1 rappdirs_0.3.3
[79] circlize_0.4.15 doParallel_1.0.17 UpSetR_1.4.0
[82] nonnest2_0.5-6 diptest_0.76-0 DBI_1.1.3
[85] shape_1.4.6 VGAM_1.1-7 httr_1.4.6
[88] spatstat.data_3.0-1 parallelly_1.36.0 xml2_1.3.5
[91] DDRTree_0.1.5 mnormt_2.1.1 ZIM_1.1.0
[94] deldir_1.0-6 interp_1.1-3 fastmap_1.1.1
[97] ExperimentHub_2.2.1 spatstat.geom_3.2-1 glue_1.6.2
[100] flexmix_2.3-18 tidygraph_1.2.1 rlang_1.1.1
[103] viridis_0.6.2 rARPACK_0.11-0 bbmle_1.0.25
[106] maxLik_1.5-2 labeling_0.4.2 plotly_4.10.0
[109] gdata_2.19.0 spatstat.sparse_3.0-1 clue_0.3-61
[112] arm_1.13-1 mvtnorm_1.1-3 cli_3.6.1
[115] zingeR_0.1.0 minqa_1.2.4 hierarchicell_1.0.0
[118] fpc_2.2-9 ash_1.0-15 glmmTMB_1.1.4
[121] latticeExtra_0.6-30 RCurl_1.98-1.7 openxlsx_4.2.5.2
[124] polyclip_1.10-0 pillar_1.9.0 tensorA_0.36.2
[127] gplots_3.1.3 prettyunits_1.1.1 R6_2.5.1
[130] ROCR_1.0-11 bit64_4.0.5 plyr_1.8.8
[133] ggstance_0.3.6 gtable_0.3.1 estimability_1.4.1
[136] xfun_0.39 reshape2_1.4.4 restfulr_0.0.15
[139] tweenr_1.0.2 miniUI_0.1.1.1 gtools_3.9.3
[142] slam_0.1-50 base64url_1.4 scattermore_0.8
[145] lifecycle_1.0.3 gmodels_2.18.1.1 splines_4.1.1
[148] fdrtool_1.2.17 BSgenome_1.62.0 spatstat.core_2.4-4
[151] KernSmooth_2.23-20 testthat_3.1.10 scatterpie_0.1.7
[154] aroma.light_3.24.0 derfinderHelper_1.36.0 globals_0.16.2
[157] sctransform_0.3.3 spatstat.utils_3.0-3 backports_1.4.1
[160] viridisLite_0.4.1 utf8_1.2.3 jsonlite_1.8.4
[163] xlsxjars_0.6.1 blockmodeling_1.1.5 rpart_4.1-15
[166] GetoptLong_1.0.5 Rtsne_0.16 prabclus_2.3-2
[169] bayNorm_1.20.0 fitdistrplus_1.1-8 stringi_1.7.8
[172] generics_0.1.3 GenomicFiles_1.38.0 ggfun_0.0.6
[175] cobs_1.3-5 zlibbioc_1.40.0 htmlTable_2.4.1
[178] brew_1.0-8 EDASeq_2.28.0 shinydashboard_0.7.2
[181] hms_1.1.3 png_0.1-8 pbivnorm_0.6.0
[184] ellipse_0.4.3 DO.db_2.9 coda_0.19-4
[187] pbapply_1.5-0 rmarkdown_2.19 lazyeval_0.2.2
[190] fields_15.2 patchwork_1.1.2 numDeriv_2016.8-1.1
[193] sandwich_3.0-2 htmlwidgets_1.6.2 annotate_1.72.0
[196] scDD_1.26.0 R.oo_1.25.0 bluster_1.4.0
[199] beachmat_2.10.0 BiocNeighbors_1.12.0 bumphunter_1.36.0
[202] renv_0.15.5 rbibutils_2.2.8 caTools_1.18.2
[205] rsvd_1.0.5 beeswarm_0.4.0 BiocIO_1.4.0
[208] ggforce_0.3.3 zip_2.3.0 R.methodsS3_1.8.2
[211] spatstat.random_3.1-5 apcluster_1.4.11 fansi_1.0.4
[214] cluster_2.1.2 RhpcBLASctl_0.21-247.1 bit_4.0.5
[217] Linnorm_2.26.0 rjson_0.2.21 shadowtext_0.1.2
[220] systemfonts_1.0.4 robustbase_0.95-0 knitr_1.43
[223] lavaan_0.6-17 brio_1.1.3 bdsmatrix_1.3-6
[226] future_1.33.0 IHW_1.30.0 tools_4.1.1
[229] TMB_1.9.1 data.table_1.14.8 RcppAnnoy_0.0.19
[232] rhdf5filters_1.6.0 EBSeq_2.0.0 DrImpute_1.0
[235] dotCall64_1.1-1 statmod_1.4.36 plotrix_3.8-2
[238] GlobalOptions_0.1.2 qvalue_2.26.0 Rcpp_1.0.11
[241] penalized_0.9-52 ica_1.0-3 softImpute_1.4-1
[244] RcppEigen_0.3.3.9.4 combinat_0.0-8 bitops_1.0-7
[247] quadprog_1.5-8 digest_0.6.33 rngtools_1.5.2
[250] downloader_0.4 Rhdf5lib_1.16.0 TH.data_1.1-1
[253] yulab.utils_0.0.5 VariantAnnotation_1.40.0 ScaledMatrix_1.2.0
[256] rstudioapi_0.15.0 BPSC_0.99.2 parallel_4.1.1
[259] XVector_0.34.0 magrittr_2.0.3 sparsesvd_0.2-2
[262] svglite_2.1.0 doRNG_1.8.2 dplyr_1.1.2
[265] bayesm_3.1-6 rhdf5_2.38.1 ggbeeswarm_0.6.0
[268] affyio_1.64.0 geneplotter_1.72.0 gamlss_5.4-3
[271] CompQuadForm_1.4.3 ggraph_2.0.5 modeltools_0.2-23
[274] igraph_1.3.2 Rttf2pt1_1.3.10 scPower_1.0.4
[277] vctrs_0.6.3 uwot_0.1.11 zoo_1.8-12
[280] cachem_1.0.6 Biostrings_2.62.0 colorspace_2.0-3
[283] monocle_2.22.0 ShortRead_1.52.0 curl_5.0.1
[286] recount_1.28.0 httpuv_1.6.11 carData_3.0-5
[289] BiocVersion_3.14.0 GenomicAlignments_1.30.0 aplot_0.1.6
[292] gamlss.data_6.0-2 base64enc_0.1-3 class_7.3-19
[295] ProtGenerics_1.26.0 derfinder_1.36.0 DelayedMatrixStats_1.16.0
[298] permute_0.9-7 nnet_7.3-16 ComplexHeatmap_2.10.0
[301] pscl_1.5.5 later_1.3.1 extrafontdb_1.0
[304] MatrixModels_0.5-0 hwriter_1.3.2.1 tensor_1.5
[307] DT_0.23 spam_2.10-0 kernlab_0.9-31
[310] hexbin_1.28.3 rgeos_0.5-9 car_3.1-0
[313] xgboost_1.6.0.1 checkmate_2.1.0 rvest_1.0.3
[316] compiler_4.1.1 xtable_1.8-4 Rmagic_2.0.3
[319] rJava_1.0-6 farver_2.1.1 extrafont_0.18
[322] minpack.lm_1.2-4 codetools_0.2-18 rstatix_0.7.0
[325] gridGraphics_0.5-1 evaluate_0.19 yaml_2.3.6
[328] goftest_1.2-3 rentrez_1.2.3 preprocessCore_1.56.0
[331] ggtree_3.2.1 docopt_0.7.1 mixtools_1.2.0
[334] GenomeInfoDbData_1.2.7 RSpectra_0.16-1 promises_1.2.0.1
[337] tidytree_0.3.9 RSQLite_2.3.1 scone_1.26.0
[340] ggsignif_0.6.3 treeio_1.21.0 FNN_1.1.3.1
[343] matlab_1.0.4 lmtest_0.9-40 BiocSingular_1.10.0
[346] webshot_0.5.3 baySeq_2.36.0 heatmap3_1.1.9
[349] vegan_2.6-2 shinyBS_0.61.1 readr_2.1.4
[352] rtracklayer_1.54.0 boot_1.3-28 metapod_1.2.0
[355] R.utils_2.12.0 ggdendro_0.1.23 HSMMSingleCell_1.14.0
[358] ggrastr_1.0.1 KEGGREST_1.34.0 iterators_1.0.14
[361] proj4_1.0-11 blme_1.0-5 ggalt_0.4.0
[364] progress_1.2.2 compositions_2.0-7 irlba_2.3.5
[367] foreign_0.8-81 mime_0.12 miscTools_0.6-26
[370] outliers_0.15 pkgconfig_2.0.3
Hi,
I am trying to perform power calculation using the powsimR. I ran the analysis as specified in the vignette but I keep encountering an error after running the simulateDE function:
I encounter this error after running the simulateDE function:
Could I please get some help and advice about how to address this problem?
(1) Although I have input the countData when I ran the estimateParam function at the start, I am unsure where to input the colData so that ncol(countData) == nrow(colData) would be TRUE. Or is there something I have not done correctly in the workflow?
(2) On a related note, is there way to include a model design since the error mentions "~ design"?
Thanks very much for your help.
Session info