Hi,

I have an RNA-seq experiment where there are 2 conditions and 2 genotypes. I am trying to figure out how to output the 2 conditions with 2 genotypes from the dds object. I have read online resources, however, it is still not clear what is extracted. I followed and compared the results from the example within the R Deseq2 documentation, however, depending on how the results are saved, the outcome is wildly different.

In the below code, within the dds design option, the 'variant' class contains "wild type" and "mutant" genotypes, and the "Type" class contains "genomic" and "Total".

I would like to compare 'Total' versus 'genomic' variants for each "wild type" and "mutant" genotypes, and then compare total/genomic logFC of 'mutant' versus 'wildtype' genotypes.

How should I write out the output for this comparison? I have 2 output examples in the code below, which are based on the examples, however it is not clear what exact comparison is being done. Any feedback will be appreciated.

I am using the code below:

#Import data
counts <-  read.table(file= "counts_for_DEseq.csv", sep=",", header=T, row.names=1)
metadata <-  read.table(file="metadata_for_DEseq.csv", sep=",", header=T)

#Note the design option 
dds <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~ variant+Type+variant:Type ) 

#Pre-filtering : min 10 reads in at least 3 samples
smallestGroupSize <- 3
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep,]

#specifying the level that represents the control group
dds$condition <- relevel(dds$Type, ref = "genomic")

# estimate size factors
dds <- estimateSizeFactors(dds)

# calculate normalized counts
norm.counts <- counts(dds, normalized=TRUE)

#Run DEseq with an interaction term
dds <- DESeq(dds)

#Output #1
MutantEffect <- results(dds, list( c("Type_total_vs_genomic","Mutant.Typetotal") ))
write.csv(MutantEffect, file= "MutantversusWildtype_total_genomic_DESeqRes.csv")

#Output #2
Cond_Eff_Across_Gen <- results(dds, name="variantmutant.Typetotal")
write.csv(Cond_Eff_Across_Gen, file= "Cond_Eff_Across_Gen_interactionTerm_DEseqRes.csv")



# include your problematic code here with any corresponding output 
# please also include the results of running the following in an R session 

sessionInfo( )
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

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.0/Resources/lib/libRlapack.dylib

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

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

other attached packages:
 [1] DESeq2_1.28.1               SummarizedExperiment_1.18.2
 [3] DelayedArray_0.14.1         matrixStats_0.61.0         
 [5] Biobase_2.48.0              GenomicRanges_1.40.0       
 [7] GenomeInfoDb_1.24.2         IRanges_2.22.2             
 [9] S4Vectors_0.26.1            BiocGenerics_0.34.0        
[11] shinyBS_0.61                shinydashboard_0.7.2       
[13] shinyjs_2.0.0               jsonlite_1.7.2             
[15] shiny_1.7.1                 plotly_4.10.0              
[17] debrowser_1.16.3            pheatmap_1.0.12            
[19] RColorBrewer_1.1-2          ggplot2_3.3.5              
[21] gplots_3.1.1               

loaded via a namespace (and not attached):
  [1] fastmatch_1.1-3        plyr_1.8.6            
  [3] igraph_1.2.8           lazyeval_0.2.2        
  [5] splines_4.0.2          Harman_1.16.0         
  [7] BiocParallel_1.22.0    pathview_1.28.1       
  [9] sva_3.36.0             urltools_1.7.3        
 [11] digest_0.6.28          foreach_1.5.1         
 [13] yulab.utils_0.0.4      htmltools_0.5.2       
 [15] GOSemSim_2.14.2        viridis_0.6.2         
 [17] GO.db_3.11.4           fansi_0.5.0           
 [19] magrittr_2.0.1         memoise_2.0.0         
 [21] limma_3.44.3           Biostrings_2.56.0     
 [23] annotate_1.66.0        graphlayouts_0.7.2    
 [25] enrichplot_1.8.1       prettyunits_1.1.1     
 [27] colorspace_2.0-2       blob_1.2.2            
 [29] ggrepel_0.9.1          dplyr_1.0.7           
 [31] crayon_1.4.2           RCurl_1.98-1.5        
 [33] graph_1.66.0           scatterpie_0.1.7      
 [35] org.Mm.eg.db_3.11.4    genefilter_1.70.0     
 [37] survival_3.2-13        iterators_1.0.13      
 [39] glue_1.5.0             polyclip_1.10-0       
 [41] registry_0.5-1         gtable_0.3.0          
 [43] zlibbioc_1.34.0        XVector_0.28.0        
 [45] webshot_0.5.2          Rgraphviz_2.32.0      
 [47] scales_1.1.1           DOSE_3.14.0           
 [49] edgeR_3.30.3           DBI_1.1.1             
 [51] miniUI_0.1.1.1         Rcpp_1.0.7            
 [53] viridisLite_0.4.0      xtable_1.8-4          
 [55] progress_1.2.2         gridGraphics_0.5-1    
 [57] bit_4.0.4              europepmc_0.4.1       
 [59] DT_0.20                htmlwidgets_1.5.4     
 [61] httr_1.4.2             fgsea_1.14.0          
 [63] ellipsis_0.3.2         pkgconfig_2.0.3       
 [65] XML_3.99-0.8           farver_2.1.0          
 [67] locfit_1.5-9.4         utf8_1.2.2            
 [69] ggplotify_0.1.0        tidyselect_1.1.1      
 [71] rlang_0.4.12           reshape2_1.4.4        
 [73] later_1.3.0            AnnotationDbi_1.50.3  
 [75] munsell_0.5.0          tools_4.0.2           
 [77] cachem_1.0.6           downloader_0.4        
 [79] generics_0.1.1         RSQLite_2.2.8         
 [81] ggridges_0.5.3         stringr_1.4.0         
 [83] fastmap_1.1.0          heatmaply_1.3.0       
 [85] org.Hs.eg.db_3.11.4    bit64_4.0.5           
 [87] tidygraph_1.2.0        caTools_1.18.2        
 [89] purrr_0.3.4            KEGGREST_1.28.0       
 [91] dendextend_1.15.2      ggraph_2.0.5          
 [93] nlme_3.1-153           mime_0.12             
 [95] KEGGgraph_1.48.0       DO.db_2.9             
 [97] xml2_1.3.2             rstudioapi_0.13       
 [99] compiler_4.0.2         png_0.1-7             
[101] tibble_3.1.6           tweenr_1.0.2          
[103] geneplotter_1.66.0     stringi_1.7.5         
[105] lattice_0.20-45        Matrix_1.3-4          
[107] vctrs_0.3.8            pillar_1.6.4          
[109] lifecycle_1.0.1        BiocManager_1.30.16   
[111] triebeard_0.3.0        data.table_1.14.2     
[113] cowplot_1.1.1          bitops_1.0-7          
[115] seriation_1.3.1        httpuv_1.6.3          
[117] qvalue_2.20.0          R6_2.5.1              
[119] promises_1.2.0.1       TSP_1.1-11            
[121] KernSmooth_2.23-20     gridExtra_2.3         
[123] codetools_0.2-18       colourpicker_1.1.1    
[125] MASS_7.3-54            gtools_3.9.2          
[127] assertthat_0.2.1       withr_2.4.2           
[129] GenomeInfoDbData_1.2.3 mgcv_1.8-38           
[131] hms_1.1.1              clusterProfiler_3.16.1
[133] ggfun_0.0.4            grid_4.0.2            
[135] tidyr_1.1.4            rvcheck_0.2.1         
[137] ggforce_0.3.3

I would like to compare 'Total' versus 'genomic' variants for each "wild type" and "mutant" genotypes

The usual and far easier to read way to do that is to make a new column where you concatenate genotype and treatment, and make that column your design, and use contrasts to pick what subgroup to compare to what.

and then compare total/genomic logFC of 'mutant' versus 'wildtype' genotypes.

I believe that's what your second one will do. But you can confirm this by doing the subgroup comparisons first, and see if the LFC you get by manually subtracting the two LFC's matches what you get from the design with interactions.