Interpretation of DESeq2 results with 3 factorial design with interaction terms.
Entering edit mode
s-rhussien • 0
Last seen 2.2 years ago

Assalam Alykom,

I'm doing a 3 factorial DESeq2 analysis with interaction terms, I want to extract results to answer specific biological questions. Here is the code for the design:

dds<-DESeqDataSetFromMatrix(countData = CountData, colData = metaData, design = ~ age + gender + viralLoad + age:viralLoad + gender:viralLoad)
dds <- DESeq(dds)

# [1] "Intercept"               "age_old_vs_young"        "gender_F_vs_M"           "viralLoad_high_vs_N.A"  
# [5] "viralLoad_low_vs_N.A"    "viralLoad_medium_vs_N.A" "ageold.viralLoadhigh"    "ageold.viralLoadlow"    
# [9] "ageold.viralLoadmedium"  "genderF.viralLoadhigh"   "genderF.viralLoadlow"    "genderF.viralLoadmedium"

If I want to answer the following question: 1- For females, what are the DEGs in high viral load vs. control? Reading the example here, I drew the following design graph (reference levels are: age: young, viral load: N.A, gender: M): design graph I use the following coefficients: viralLoad_high_vs_N.A + genderF.viralLoadhigh. I use the following code to extract this data:

res <- results(dds, contrast=c(0,0,0,1,0,0,0,0,0,1,0,0), alpha=0.05)

But what if I want to answer: 2- For males, what are the DEGs in high viral load vs. control? How should I do this? as there is no interaction term such as: genderM.viralLoadhigh in the list of resultsNames(dds). Could I use something like: viralLoad_high_vs_N.A - gender_F_vs_M - genderF.viralLoadhigh ?

sessionInfo( )

R version 4.0.4 (2021-02-15)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

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

other attached packages:
 [1] pheatmap_1.0.12             RColorBrewer_1.1-2          EnhancedVolcano_1.8.0      
 [4] ggrepel_0.9.1               forcats_0.5.1               stringr_1.4.0              
 [7] purrr_0.3.4                 readr_2.0.1                 tidyr_1.1.3                
[10] tibble_3.1.3                ggplot2_3.3.5               tidyverse_1.3.1            
[13] DESeq2_1.30.1               SummarizedExperiment_1.20.0 MatrixGenerics_1.2.1       
[16] matrixStats_0.60.0          dplyr_1.0.7                 GEOquery_2.58.0            
[19] Biobase_2.50.0              rtracklayer_1.49.5          GenomicRanges_1.42.0       
[22] GenomeInfoDb_1.26.7         Biostrings_2.58.0           XVector_0.30.0             
[25] IRanges_2.24.1              S4Vectors_0.28.1            AnnotationHub_2.22.1       
[28] BiocFileCache_1.14.0        dbplyr_2.1.1                BiocGenerics_0.36.1        

loaded via a namespace (and not attached):
  [1] ggbeeswarm_0.6.0              colorspace_2.0-2              ellipsis_0.3.2               
  [4] fs_1.5.0                      rstudioapi_0.13               bit64_4.0.5                  
  [7] interactiveDisplayBase_1.28.0 AnnotationDbi_1.52.0          fansi_0.5.0                  
 [10] lubridate_1.7.10              xml2_1.3.2                    splines_4.0.4                
 [13] extrafont_0.17                cachem_1.0.5                  geneplotter_1.68.0           
 [16] knitr_1.33                    jsonlite_1.7.2                Rsamtools_2.6.0              
 [19] Rttf2pt1_1.3.9                broom_0.7.9                   annotate_1.68.0              
 [22] shiny_1.6.0                   BiocManager_1.30.16           compiler_4.0.4               
 [25] httr_1.4.2                    backports_1.2.1               assertthat_0.2.1             
 [28] Matrix_1.3-4                  fastmap_1.1.0                 limma_3.46.0                 
 [31] cli_3.0.1                     later_1.3.0                   htmltools_0.5.1.1            
 [34] tools_4.0.4                   gtable_0.3.0                  glue_1.4.2                   
 [37] GenomeInfoDbData_1.2.4        maps_3.3.0                    rappdirs_0.3.3               
 [40] tinytex_0.33                  Rcpp_1.0.7                    cellranger_1.1.0             
 [43] vctrs_0.3.8                   ggalt_0.4.0                   extrafontdb_1.0              
 [46] xfun_0.23                     rvest_1.0.1                   mime_0.11                    
 [49] lifecycle_1.0.0               XML_3.99-0.6                  MASS_7.3-54                  
 [52] zlibbioc_1.36.0               scales_1.1.1                  vroom_1.5.4                  
 [55] hms_1.1.0                     promises_1.2.0.1              proj4_1.0-10.1               
 [58] yaml_2.2.1                    curl_4.3.2                    memoise_2.0.0                
 [61] ggrastr_0.2.3                 stringi_1.7.4                 RSQLite_2.2.7                
 [64] BiocVersion_3.12.0            genefilter_1.72.1             BiocParallel_1.24.1          
 [67] rlang_0.4.11                  pkgconfig_2.0.3               bitops_1.0-7                 
 [70] evaluate_0.14                 lattice_0.20-44               GenomicAlignments_1.26.0     
 [73] bit_4.0.4                     tidyselect_1.1.1              magrittr_2.0.1               
 [76] R6_2.5.1                      generics_0.1.0                DelayedArray_0.16.3          
 [79] DBI_1.1.1                     pillar_1.6.2                  haven_2.4.3                  
 [82] withr_2.4.2                   ash_1.0-15                    survival_3.2-11              
 [85] RCurl_1.98-1.3                modelr_0.1.8                  crayon_1.4.1                 
 [88] KernSmooth_2.23-20            utf8_1.2.2                    tzdb_0.1.2                   
 [91] rmarkdown_2.10                locfit_1.5-9.4                grid_4.0.4                   
 [94] readxl_1.3.1                  blob_1.2.2                    reprex_2.0.1                 
 [97] digest_0.6.27                 xtable_1.8-4                  httpuv_1.6.2                 
[100] munsell_0.5.0                 beeswarm_0.4.0                vipor_0.4.5
r DESeq DESeq2 • 1.4k views
Entering edit mode
Last seen 12 hours ago
United States

Unfortunately, I only have time to answer software related questions on the support site. For statistical consultation about design and contrasts, I recommend meeting with a local statistician at your institute.

Entering edit mode

Thanks very much

Entering edit mode
swbarnes2 ★ 1.3k
Last seen 4 days ago
San Diego

I think the simplest way to do your first question is to make a new column that joins sex and load, make your design sex_load + age, and then use contrasts to specify that you want to compare F_high to F_control.

See here. I think your questions fall in the category of "Can be done simpler without interactions"

Entering edit mode

You mean that if I have a factor with more than 2 levels, I can compare any 2 levels of that factor using contrast=c(factor_name, 1st_desired_level, 2nd_desired_level)?

Entering edit mode

Yes, that's just how the vignette says to use contrasts.

Entering edit mode

Thanks very much


Login before adding your answer.

Traffic: 885 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6