lfcshrink error DESeq2
1
0
Entering edit mode
Alina • 0
@36288a68
Last seen 17 months ago
Spain

Hello!

I'm having problems with lfcShrink in my DESeq2 workflow.

I'm trying to do a differential expression analysis (with only one comparison term: "MULTIseq_ID_call2") on my single-cell data. However when I do lfcShrink I get an error that I cannot interpret.

Can you help me?

dds <- DESeq(dds, test = "LRT", reduced=~1, useT=TRUE, minmu=1e-6,minReplicatesForReplace=Inf, fitType = "glmGamPoi")

res <- results(dds,
               alpha = 0.05)

res <- lfcShrink(dds, 
                 coef="MULTIseq_ID_call2iRBC_extract",
                  type="apeglm")

Error:

Error in if (objective(min.var) < 0) { :
missing value where TRUE/FALSE needed
3.
priorVar(mle)
2.
apeglm::apeglm(Y = Y, x = design, log.lik = log.lik, param = disps,
coef = coefNum, mle = mle, threshold = apeT, weights = weights,
offset = offset, method = apeMethod, ...)
1.
lfcShrink(dds, coef = "MULTIseq_ID_call2iRBC_extract", type = "apeglm",
res = res)

> sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Rocky Linux 8.7 (Green Obsidian)

Matrix products: default
BLAS/LAPACK: /g/easybuild/x86_64/Rocky/8/haswell/software/FlexiBLAS/3.2.0-GCC-11.3.0/lib64/libflexiblas.so.3.2

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] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ashr_2.2-54                 apeglm_1.18.0               glmGamPoi_1.8.0            
 [4] DESeq2_1.36.0               pheatmap_1.0.12             MAST_1.22.0                
 [7] ggrepel_0.9.1               NMF_0.25                    bigmemory_4.6.1            
[10] cluster_2.1.3               rngtools_1.5.2              registry_0.5-1             
[13] data.table_1.14.2           clusterProfiler_4.4.4       msigdbr_7.5.1              
[16] SingleR_1.10.0              celldex_1.6.0               writexl_1.4.0              
[19] bluster_1.6.0               scDblFinder_1.10.0          scran_1.24.1               
[22] scater_1.24.0               scuttle_1.6.3               SingleCellExperiment_1.18.1
[25] SummarizedExperiment_1.26.1 Biobase_2.56.0              GenomicRanges_1.48.0       
[28] GenomeInfoDb_1.32.2         IRanges_2.30.0              S4Vectors_0.34.0           
[31] BiocGenerics_0.42.0         MatrixGenerics_1.8.1        matrixStats_0.62.0         
[34] umap_0.2.9.0                reticulate_1.26             deMULTIplex_1.0.2          
[37] SeuratObject_4.1.3          Seurat_4.2.1                limma_3.52.2               
[40] forcats_0.5.1               stringr_1.4.0               dplyr_1.0.9                
[43] purrr_0.3.4                 readr_2.1.2                 tidyr_1.2.0                
[46] tibble_3.1.7                tidyverse_1.3.1             RColorBrewer_1.1-3         
[49] xlsx_0.6.5                  readxl_1.4.0                rstatix_0.7.0              
[52] ggpubr_0.4.0                ggplot2_3.4.0              

loaded via a namespace (and not attached):
  [1] rsvd_1.0.5                    ica_1.0-2                     Rsamtools_2.12.0             
  [4] foreach_1.5.2                 lmtest_0.9-40                 crayon_1.5.1                 
  [7] MASS_7.3-57                   nlme_3.1-158                  backports_1.4.1              
 [10] reprex_2.0.1                  GOSemSim_2.22.0               rlang_1.0.6                  
 [13] XVector_0.36.0                ROCR_1.0-11                   irlba_2.3.5                  
 [16] filelock_1.0.2                xgboost_1.6.0.1               BiocParallel_1.30.4          
 [19] rjson_0.2.21                  bit64_4.0.5                   glue_1.6.2                   
 [22] mixsqp_0.3-43                 sctransform_0.3.5             parallel_4.2.1               
 [25] vipor_0.4.5                   spatstat.sparse_3.0-0         AnnotationDbi_1.58.0         
 [28] DOSE_3.22.0                   spatstat.geom_3.0-3           haven_2.5.0                  
 [31] tidyselect_1.1.2              fitdistrplus_1.1-8            XML_3.99-0.10                
 [34] zoo_1.8-10                    GenomicAlignments_1.32.0      xtable_1.8-4                 
 [37] magrittr_2.0.3                evaluate_0.15                 cli_3.4.1                    
 [40] zlibbioc_1.42.0               rstudioapi_0.13               miniUI_0.1.1.1               
 [43] sp_1.5-1                      fastmatch_1.1-3               treeio_1.20.1                
 [46] shiny_1.7.1                   BiocSingular_1.12.0           xfun_0.31                    
 [49] askpass_1.1                   tidygraph_1.2.1               KEGGREST_1.36.2              
 [52] interactiveDisplayBase_1.34.0 ape_5.6-2                     listenv_0.8.0                
 [55] xlsxjars_0.6.1                Biostrings_2.64.0             png_0.1-7                    
 [58] future_1.26.1                 withr_2.5.0                   bitops_1.0-7                 
 [61] ggforce_0.3.3                 plyr_1.8.7                    cellranger_1.1.0             
 [64] coda_0.19-4                   dqrng_0.3.0                   pillar_1.7.0                 
 [67] cachem_1.0.6                  fs_1.5.2                      DelayedMatrixStats_1.18.0    
 [70] vctrs_0.5.1                   ellipsis_0.3.2                generics_0.1.2               
 [73] tools_4.2.1                   beeswarm_0.4.0                munsell_0.5.0                
 [76] tweenr_1.0.2                  fgsea_1.22.0                  DelayedArray_0.22.0          
 [79] fastmap_1.1.0                 compiler_4.2.1                abind_1.4-5                  
 [82] httpuv_1.6.5                  rtracklayer_1.56.1            ExperimentHub_2.4.0          
 [85] plotly_4.10.1                 rJava_1.0-6                   GenomeInfoDbData_1.2.8       
 [88] gridExtra_2.3                 edgeR_3.38.1                  lattice_0.20-45              
 [91] deldir_1.0-6                  utf8_1.2.2                    later_1.3.0                  
 [94] BiocFileCache_2.4.0           jsonlite_1.8.0                scales_1.2.0                 
 [97] ScaledMatrix_1.4.0            tidytree_0.3.9                pbapply_1.5-0                
[100] carData_3.0-5                 sparseMatrixStats_1.8.0       genefilter_1.78.0            
[103] lazyeval_0.2.2                promises_1.2.0.1              doParallel_1.0.17            
[106] car_3.1-0                     goftest_1.2-3                 spatstat.utils_3.0-1         
[109] rmarkdown_2.14                cowplot_1.1.1                 statmod_1.4.36               
[112] Rtsne_0.16                    downloader_0.4                uwot_0.1.14                  
[115] igraph_1.3.2                  numDeriv_2016.8-1.1           survival_3.3-1               
[118] yaml_2.3.5                    SQUAREM_2021.1                htmltools_0.5.2              
[121] memoise_2.0.1                 BiocIO_1.6.0                  locfit_1.5-9.5               
[124] graphlayouts_0.8.0            viridisLite_0.4.0             digest_0.6.29                
[127] assertthat_0.2.1              mime_0.12                     rappdirs_0.3.3               
[130] emdbook_1.3.12                bigmemory.sri_0.1.3           RSQLite_2.2.14               
[133] yulab.utils_0.0.5             future.apply_1.9.0            blob_1.2.3                   
[136] labeling_0.4.2                splines_4.2.1                 AnnotationHub_3.4.0          
[139] RCurl_1.98-1.7                broom_0.8.0                   hms_1.1.1                    
[142] modelr_0.1.8                  colorspace_2.0-3              BiocManager_1.30.18          
[145] ggbeeswarm_0.6.0              aplot_0.1.6                   Rcpp_1.0.8.3                 
[148] RANN_2.6.1                    mvtnorm_1.1-3                 enrichplot_1.16.2            
[151] fansi_1.0.3                   tzdb_0.3.0                    truncnorm_1.0-8              
[154] parallelly_1.32.0             R6_2.5.1                      grid_4.2.1                   
[157] ggridges_0.5.3                lifecycle_1.0.3               curl_4.3.2                   
[160] ggsignif_0.6.3                leiden_0.4.2                  DO.db_2.9                    
[163] Matrix_1.5-1                  qvalue_2.28.0                 RcppAnnoy_0.0.19             
[166] iterators_1.0.14              spatstat.explore_3.0-3        htmlwidgets_1.5.4            
[169] beachmat_2.12.0               polyclip_1.10-0               shadowtext_0.1.2             
[172] gridGraphics_0.5-1            rvest_1.0.2                   globals_0.15.0               
[175] openssl_2.0.2                 patchwork_1.1.2               spatstat.random_3.0-1        
[178] bdsmatrix_1.3-6               progressr_0.10.1              invgamma_1.1                 
[181] codetools_0.2-18              lubridate_1.8.0               GO.db_3.15.0                 
[184] metapod_1.4.0                 dbplyr_2.2.0                  gridBase_0.4-7               
[187] RSpectra_0.16-1               gtable_0.3.0                  DBI_1.1.3                    
[190] ggfun_0.0.6                   tensor_1.5                    httr_1.4.3                   
[193] KernSmooth_2.23-20            stringi_1.7.6                 uuid_1.1-0                   
[196] reshape2_1.4.4                farver_2.1.0                  annotate_1.74.0              
[199] viridis_0.6.2                 ggtree_3.4.1                  xml2_1.3.3                   
[202] bbmle_1.0.25                  BiocNeighbors_1.14.0          restfulr_0.0.15              
[205] geneplotter_1.74.0            ggplotify_0.1.0               scattermore_0.8              
[208] BiocVersion_3.15.2            bit_4.0.4                     scatterpie_0.1.7             
[211] spatstat.data_3.0-0           ggraph_2.0.5                  pkgconfig_2.0.3              
[214] babelgene_22.3                knitr_1.39
DESeq2 lfcshrink DESe • 2.4k views
ADD COMMENT •
0
Entering edit mode

I bet this is because apeglm required a lfcSE which is NA though when using apeglm, right Michael Love ?

ADD REPLY •
0
Entering edit mode

yes, lfcSE is NA. Can I change that?

ADD REPLY •
0
Entering edit mode

This case is now addressed (gives an informative error) in devel:

https://github.com/mikelove/DESeq2/issues/73

You can test that the error works with devtools::install_github("mikelove/DESeq2")

ADD REPLY •
0
Entering edit mode
@mikelove
Last seen 7 hours ago
United States

What is resultsNames(dds)
ADD COMMENT •
0
Entering edit mode
> resultsNames(dds)
[1] "Intercept"                     "MULTIseq_ID_call2iRBC_extract"
ADD REPLY •
0
Entering edit mode

Is the general recommendation to always use "LRT" on single-cell data? Even if there is only one comparison factor?

ADD REPLY •
0
Entering edit mode

You can shrink coefficients listed in resultsNames, see ?lfcShrink

Yes, LRT recommended for scRNA-seq

ADD REPLY •
0
Entering edit mode

Thank you!

So it should be possible with coef="MULTIseq_ID_call2iRBC_extract"?

Is the problem that my lfcSE column is NA?

ADD REPLY •
0
Entering edit mode

Yeah, lfcSE needs to not be NA.

You can use a different fitType (the default for example).

ADD REPLY •
0
Entering edit mode

^ my advice is here.

ADD REPLY •
0
Entering edit mode

Thanks a lot!

ADD REPLY •
0
Entering edit mode

It's a reproducible error with glmGamPoi+apeglm causing it:

suppressMessages(library(DESeq2))
dds <- makeExampleDESeqDataSet()
dds <- DESeq(dds, test="LRT", reduced=~1, fitType="glmGamPoi", quiet=TRUE)
res <- lfcShrink(dds, coef=resultsNames(dds)[2], type="apeglm", quiet=TRUE)
#> Error in if (objective(min.var) < 0) {: Fehlender Wert, wo TRUE/FALSE nötig ist
sessionInfo()
#> R version 4.2.2 (2022-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22000)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8   
#> [3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C                   
#> [5] LC_TIME=German_Germany.utf8    
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] DESeq2_1.38.3               SummarizedExperiment_1.28.0
#>  [3] Biobase_2.58.0              MatrixGenerics_1.10.0      
#>  [5] matrixStats_0.63.0          GenomicRanges_1.50.1       
#>  [7] GenomeInfoDb_1.34.9         IRanges_2.32.0             
#>  [9] S4Vectors_0.36.1            BiocGenerics_0.44.0        
#> 
#> loaded via a namespace (and not attached):
#>  [1] httr_1.4.4                bit64_4.0.5              
#>  [3] DelayedMatrixStats_1.20.0 highr_0.9                
#>  [5] blob_1.2.3                GenomeInfoDbData_1.2.9   
#>  [7] yaml_2.3.6                numDeriv_2016.8-1.1      
#>  [9] pillar_1.9.0              RSQLite_2.2.19           
#> [11] lattice_0.20-45           glue_1.6.2               
#> [13] bbmle_1.0.25              digest_0.6.30            
#> [15] RColorBrewer_1.1-3        XVector_0.38.0           
#> [17] colorspace_2.0-3          plyr_1.8.8               
#> [19] htmltools_0.5.4           Matrix_1.5-1             
#> [21] XML_3.99-0.13             pkgconfig_2.0.3          
#> [23] emdbook_1.3.12            zlibbioc_1.44.0          
#> [25] mvtnorm_1.1-3             xtable_1.8-4             
#> [27] scales_1.2.1              apeglm_1.20.0            
#> [29] BiocParallel_1.32.4       tibble_3.2.1             
#> [31] annotate_1.76.0           KEGGREST_1.38.0          
#> [33] generics_0.1.3            ggplot2_3.4.2            
#> [35] cachem_1.0.6              withr_2.5.0              
#> [37] cli_3.6.1                 magrittr_2.0.3           
#> [39] crayon_1.5.2              memoise_2.0.1            
#> [41] evaluate_0.18             fs_1.5.2                 
#> [43] fansi_1.0.3               MASS_7.3-58.1            
#> [45] tools_4.2.2               lifecycle_1.0.3          
#> [47] stringr_1.5.0             munsell_0.5.0            
#> [49] reprex_2.0.2              locfit_1.5-9.6           
#> [51] DelayedArray_0.24.0       AnnotationDbi_1.60.0     
#> [53] Biostrings_2.66.0         compiler_4.2.2           
#> [55] rlang_1.1.0               grid_4.2.2               
#> [57] RCurl_1.98-1.9            rstudioapi_0.14          
#> [59] glmGamPoi_1.10.2          bitops_1.0-7             
#> [61] rmarkdown_2.18            gtable_0.3.1             
#> [63] codetools_0.2-18          DBI_1.1.3                
#> [65] R6_2.5.1                  bdsmatrix_1.3-6          
#> [67] knitr_1.41                dplyr_1.1.2              
#> [69] fastmap_1.1.0             bit_4.0.5                
#> [71] utf8_1.2.2                stringi_1.7.8            
#> [73] parallel_4.2.2            Rcpp_1.0.9               
#> [75] vctrs_0.6.2               geneplotter_1.76.0       
#> [77] png_0.1-8                 coda_0.19-4              
#> [79] tidyselect_1.2.0          xfun_0.35                
#> [81] sparseMatrixStats_1.10.0
Created on 2023-05-09 with reprex v2.0.2
ADD REPLY •
1
Entering edit mode

Ah, I can just add a check for NA standard errors here. Will put it on the list.

ADD REPLY •
1
Entering edit mode

This case is now addressed (gives an informative error) in devel:

https://github.com/mikelove/DESeq2/issues/73

ADD REPLY •
0
Entering edit mode

Thank you!!

So should I for now just not use glmGamPoi? Or is there anything else I could do? Which fit type would you recommend for my data instead?

ADD REPLY •
0
Entering edit mode

I have exactly the same issue as Alina. Which fit type should we use?

ADD REPLY •
0
Entering edit mode

Answer is above in another thread.

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 637 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6