Deseq2 resultsNames() giving header instead of values
1
0
Entering edit mode
pturpin0 • 0
@pturpin0-23836
Last seen 5.4 years ago

I recently updated some packages, which I think is interfering with the resultsNames() function. I've tried to test it on premade datasets to make sure that it's not a personal data issue, and am running into the same problem.

The below function uses a phyloseq dataset, which is similar to my dataset and showing the same issue.

> library(phyloseq)
> data(soilrep)
> prevalence_subset = soilrep

> map = data.frame(sample_data(prevalence_subset))

> OTU_data = data.frame(otu_table(prevalence_subset))
> OTU_data = OTU_data + 1
> OTU_data = otu_table(OTU_data, taxa_are_rows = TRUE)
> otu_table(prevalence_subset) = OTU_data

> diagdds_1 = phyloseq_to_deseq2(prevalence_subset, ~Treatment)
converting counts to integer mode
> diagdds_1 = DESeq(diagdds_1, test="Wald", fitType="local", quiet = TRUE)
> current_RN = resultsNames(diagdds_1)

> current_RN 
[1] "Intercept"  "Treatment1" "Treatment2" "Treatment3"
> #I'd expect it to be call on the variable names; "UC" "UU" "WU" "WC"

I've also tried running the script available on the Bioconductor subpage, to make sure it wasn't a phyloseq issue. Everything ran without error until I reached line 151, where I ran into the same problem.

> ## ----lfcShrink----------------------------------------------------------------
> resultsNames(dds)
[1] "Intercept"  "condition1"
> resLFC <- lfcShrink(dds, coef="condition_treated_vs_untreated", type="apeglm")
Error in lfcShrink(dds, coef = "condition_treated_vs_untreated", type = "apeglm") : 
  coef %in% resultsNamesDDS is not TRUE

Here are the current package versions I'm using.

> sessionInfo()
R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/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] datasets  parallel  stats4    stats     graphics  grDevices utils     methods   base     

other attached packages:
 [1] phyloseq_1.30.0             airway_1.6.0                tximeta_1.4.5               tximportData_1.14.0        
 [5] readr_1.3.1                 tximport_1.14.2             DESeq2_1.26.0               vegan_2.5-6                
 [9] lattice_0.20-41             permute_0.9-5               pime_0.1.0                  stringr_1.4.0              
[13] ggplot2_3.3.2               SummarizedExperiment_1.16.1 DelayedArray_0.12.3         BiocParallel_1.20.1        
[17] matrixStats_0.56.0          Biobase_2.46.0              GenomicRanges_1.38.0        GenomeInfoDb_1.22.1        
[21] IRanges_2.20.2              S4Vectors_0.24.4            BiocGenerics_0.32.0        

loaded via a namespace (and not attached):
  [1] backports_1.1.8          Hmisc_4.4-0              BiocFileCache_1.10.2     plyr_1.8.6               igraph_1.2.5            
  [6] lazyeval_0.2.2           splines_3.6.2            digest_0.6.25            foreach_1.5.0            ensembldb_2.10.2        
 [11] htmltools_0.5.0          magrittr_1.5             checkmate_2.0.0          memoise_1.1.0            cluster_2.1.0           
 [16] doParallel_1.0.15        Biostrings_2.54.0        annotate_1.64.0          askpass_1.1              prettyunits_1.1.1       
 [21] jpeg_0.1-8.1             colorspace_1.4-1         blob_1.2.1               rappdirs_0.3.1           xfun_0.15               
 [26] dplyr_1.0.0              crayon_1.3.4             RCurl_1.98-1.2           jsonlite_1.7.0           genefilter_1.68.0       
 [31] survival_3.2-3           iterators_1.0.12         ape_5.4                  glue_1.4.1               gtable_0.3.0            
 [36] zlibbioc_1.32.0          XVector_0.26.0           Rhdf5lib_1.8.0           scales_1.1.1             DBI_1.1.0               
 [41] Rcpp_1.0.5               xtable_1.8-4             progress_1.2.2           htmlTable_2.0.1          foreign_0.8-72          
 [46] bit_1.1-15.2             Formula_1.2-3            htmlwidgets_1.5.1        httr_1.4.1               RColorBrewer_1.1-2      
 [51] acepack_1.4.1            ellipsis_0.3.1           pkgconfig_2.0.3          XML_3.99-0.3             farver_2.0.3            
 [56] nnet_7.3-14              dbplyr_1.4.4             locfit_1.5-9.4           tidyselect_1.1.0         labeling_0.3            
 [61] rlang_0.4.6              reshape2_1.4.4           AnnotationDbi_1.48.0     munsell_0.5.0            tools_3.6.2             
 [66] generics_0.0.2           RSQLite_2.2.0            ade4_1.7-15              ranger_0.12.1            biomformat_1.14.0       
 [71] knitr_1.29               bit64_0.9-7              purrr_0.3.4              randomForest_4.6-14      AnnotationFilter_1.10.0 
 [76] nlme_3.1-148             biomaRt_2.42.1           compiler_3.6.2           rstudioapi_0.11          curl_4.3                
 [81] png_0.1-7                tibble_3.0.2             geneplotter_1.64.0       stringi_1.4.6            GenomicFeatures_1.38.2  
 [86] forcats_0.5.0            ProtGenerics_1.18.0      Matrix_1.2-18            microbiome_1.8.0         multtest_2.42.0         
 [91] vctrs_0.3.1              pillar_1.4.4             lifecycle_0.2.0          BiocManager_1.30.10      data.table_1.12.8       
 [96] bitops_1.0-6             rtracklayer_1.46.0       R6_2.4.1                 latticeExtra_0.6-29      gridExtra_2.3           
[101] codetools_0.2-16         MASS_7.3-51.6            assertthat_0.2.1         rhdf5_2.30.1             openssl_1.4.2           
[106] withr_2.2.0              GenomicAlignments_1.22.1 Rsamtools_2.2.3          GenomeInfoDbData_1.2.2   mgcv_1.8-31             
[111] hms_0.5.3                grid_3.6.2               rpart_4.1-15             tidyr_1.1.0              Rtsne_0.15              
[116] base64enc_0.1-3         

Thank you so much for your time!

deseq2 • 928 views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 2 days ago
United States

Here is my take on the above code snippets:

First snippet: resultsNames will be based on the colData of dds and on the design. Here the creation of the DESeqDataSet is actually being done within another function from phyloseq, so you'd have to see how they are converting your sample information (we call colData). But you can check what is being provided to DESeq2 by looking at colData(dds). This is the information DESeq() is working with and that is what determines the resultsNames.

Second snippet: when you specify coef it has to be the name of a character string in resultsNames(dds). The printed error gives you a hint to this. Since "condition_treated_vs_untreated" is not present in the vector: c("Intercept", "condition1"), it gives an error.

One issue is that I can't see any lines where DESeqDataSet is created (because they are within other functions, or not shown above), so I can't help resolve the DESeq2 side of things.

ADD COMMENT
0
Entering edit mode

In regards to the second snippet, I used the code from the Deseq2 Bioconductor site: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.R.

I'm sorry I wasn't clearer, I just didn't want to have to paste the whole thing in since it's so long, but I'm happy to. I've successfully run the entirety of the provided script, until line 151, as shown above.

ADD REPLY
0
Entering edit mode

I'm not sure why you're not seeing the same resultsNames(dds) as in the vignette for Bioconductor 3.10.

https://bioconductor.org/packages/3.10/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#log-fold-change-shrinkage-for-visualization-and-ranking

What do you get for:

BiocManager::valid()

ADD REPLY
0
Entering edit mode
> BiocManager::valid()

* sessionInfo()

R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/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] datasets  parallel  stats4    stats     graphics  grDevices utils     methods   base     

other attached packages:
 [1] pasilla_1.14.0              phyloseq_1.30.0             airway_1.6.0                tximeta_1.4.5              
 [5] tximportData_1.14.0         readr_1.3.1                 tximport_1.14.2             DESeq2_1.26.0              
 [9] vegan_2.5-6                 lattice_0.20-41             permute_0.9-5               pime_0.1.0                 
[13] stringr_1.4.0               ggplot2_3.3.2               SummarizedExperiment_1.16.1 DelayedArray_0.12.3        
[17] BiocParallel_1.20.1         matrixStats_0.56.0          Biobase_2.46.0              GenomicRanges_1.38.0       
[21] GenomeInfoDb_1.22.1         IRanges_2.20.2              S4Vectors_0.24.4            BiocGenerics_0.32.0        

loaded via a namespace (and not attached):
  [1] backports_1.1.8          Hmisc_4.4-0              BiocFileCache_1.10.2     plyr_1.8.6               igraph_1.2.5            
  [6] lazyeval_0.2.2           splines_3.6.2            digest_0.6.25            foreach_1.5.0            ensembldb_2.10.2        
 [11] htmltools_0.5.0          magrittr_1.5             checkmate_2.0.0          memoise_1.1.0            cluster_2.1.0           
 [16] doParallel_1.0.15        Biostrings_2.54.0        annotate_1.64.0          askpass_1.1              prettyunits_1.1.1       
 [21] jpeg_0.1-8.1             colorspace_1.4-1         blob_1.2.1               rappdirs_0.3.1           xfun_0.15               
 [26] dplyr_1.0.0              crayon_1.3.4             RCurl_1.98-1.2           jsonlite_1.7.0           genefilter_1.68.0       
 [31] survival_3.2-3           iterators_1.0.12         ape_5.4                  glue_1.4.1               gtable_0.3.0            
 [36] zlibbioc_1.32.0          XVector_0.26.0           Rhdf5lib_1.8.0           scales_1.1.1             DBI_1.1.0               
 [41] Rcpp_1.0.5               xtable_1.8-4             progress_1.2.2           htmlTable_2.0.1          foreign_0.8-72          
 [46] bit_1.1-15.2             Formula_1.2-3            htmlwidgets_1.5.1        httr_1.4.1               RColorBrewer_1.1-2      
 [51] acepack_1.4.1            ellipsis_0.3.1           pkgconfig_2.0.3          XML_3.99-0.3             farver_2.0.3            
 [56] nnet_7.3-14              dbplyr_1.4.4             locfit_1.5-9.4           tidyselect_1.1.0         labeling_0.3            
 [61] rlang_0.4.6              reshape2_1.4.4           AnnotationDbi_1.48.0     munsell_0.5.0            tools_3.6.2             
 [66] generics_0.0.2           RSQLite_2.2.0            ade4_1.7-15              ranger_0.12.1            biomformat_1.14.0       
 [71] knitr_1.29               bit64_0.9-7              purrr_0.3.4              randomForest_4.6-14      AnnotationFilter_1.10.0 
 [76] nlme_3.1-148             biomaRt_2.42.1           compiler_3.6.2           rstudioapi_0.11          curl_4.3                
 [81] png_0.1-7                tibble_3.0.2             geneplotter_1.64.0       stringi_1.4.6            GenomicFeatures_1.38.2  
 [86] forcats_0.5.0            ProtGenerics_1.18.0      Matrix_1.2-18            microbiome_1.8.0         multtest_2.42.0         
 [91] vctrs_0.3.1              pillar_1.4.5             lifecycle_0.2.0          BiocManager_1.30.10      data.table_1.12.8       
 [96] bitops_1.0-6             rtracklayer_1.46.0       R6_2.4.1                 latticeExtra_0.6-29      gridExtra_2.3           
[101] codetools_0.2-16         MASS_7.3-51.6            assertthat_0.2.1         rhdf5_2.30.1             openssl_1.4.2           
[106] withr_2.2.0              GenomicAlignments_1.22.1 Rsamtools_2.2.3          GenomeInfoDbData_1.2.2   mgcv_1.8-31             
[111] hms_0.5.3                grid_3.6.2               rpart_4.1-15             tidyr_1.1.0              Rtsne_0.15              
[116] base64enc_0.1-3         

Bioconductor version '3.10'

  * 2 packages out-of-date
  * 0 packages too new

create a valid installation with

  BiocManager::install(c(
    "quantreg", "rlang"
  ), update = TRUE, ask = FALSE)

more details: BiocManager::valid()$too_new, BiocManager::valid()$out_of_date

Warning message:
2 packages out-of-date; 0 packages too new 

When I attempted to update the out-of-date packages ("quantreg", "rlang") I run into the errors, which I can try to iron out if you think that's the problem.

ADD REPLY
0
Entering edit mode

Let's do a smaller example:

> cts <- matrix(rpois(1000,lambda=10),ncol=10)
> cdata <- data.frame(condition=factor(rep(1:2,each=5)))
> dds <- DESeqDataSetFromMatrix(cts, cdata, ~condition)
> dds <- DESeq(dds, fitType="mean", quiet=TRUE)
> resultsNames(dds)
[1] "Intercept"        "condition_2_vs_1"
ADD REPLY
0
Entering edit mode

Here is my output for the above:

> cts <- matrix(rpois(1000,lambda=10),ncol=10)
> cdata <- data.frame(condition=factor(rep(1:2,each=5)))
> dds <- DESeqDataSetFromMatrix(cts, cdata, ~condition)
> dds <- DESeq(dds, fitType="mean", quiet=TRUE)

> resultsNames(dds)
[1] "Intercept"  "condition1"

> dds
class: DESeqDataSet 
dim: 100 10 
metadata(1): version
assays(4): counts mu H cooks
rownames: NULL
rowData names(22): baseMean baseVar ... deviance maxCooks
colnames: NULL
colData names(2): condition sizeFactor

> cdata
   condition
1          1
2          1
3          1
4          1
5          1
6          2
7          2
8          2
9          2
10         2
ADD REPLY
0
Entering edit mode

Do you have anything in ~/.Rprofile that would modify the default behavior of model.matrix?

I haven't seen this before...

ADD REPLY

Login before adding your answer.

Traffic: 943 users visited in the last hour
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