Hello,

I have a question to the DESeq2 contrast parameter of the results function. I have single end reads from 16 samples with 4 treatments (group) and 4 biological replicates (indi). However, the RNA of the replicates was isolated on different days. Can I still correct for the RNA-isolation day bias (iso)?

My meta table looks like this:

group indi iso
  T1   I1   A     
  T1   I2   A     
  T1   I3   B     
  T1   I4   B     
  T2   I1   A     
  T2   I2   A    
  T2   I3   B     
  T2   I4   B     
  T3   I1   A     
  T3   I2   A     
  T3   I3   B     
  T3   I4   B     
  T4   I1   A     
  T4   I2   A     
  T4   I3   B     
  T4   I4   B

I followed the instructions for such case from the DESeq2 manual:

ds_txi <- DESeqDataSetFromTximport(txi = txi_salmon,
                                   colData = meta,
                                   design = ~ indi+group)

ds_txi$indi_n <- c("I1","I2","I1","I2","I1","I2","I1","I2","I1","I2","I1","I2","I1","I2","I1","I2")

meta$indi_n <- c("I1","I2","I1","I2","I1","I2","I1","I2","I1","I2","I1","I2","I1","I2","I1","I2")

meta$indi_n <- as.factor(meta$indi_n)
ds_txi$indi_n <- as.factor(ds_txi$indi_n)

ds_txi <- DESeqDataSetFromTximport(txi = txi_salmon,
                                   colData = meta,
                                   design = ~ iso+ iso:indi_n + iso:group)

Resulting in following meta table:

group indi iso indi_n
  T1   I1    A    I1
  T1   I2    A    I2
  T1   I3    B    I1
  T1   I4    B    I2
  T2   I1    A    I1
  T2   I2    A    I2
  T2   I3    B    I1
  T2   I4    B    I2
  T3   I1    A    I1
  T3   I2    A    I2
  T3   I3    B    I1
  T3   I4    B    I2
  T4   I1    A    I1
  T4   I2    A    I2
  T4   I3    B    I1
  T4   I4    B    I2

With following contrast, I get the difference between treatment T1 and T2 within Batch of isolation date A:

dds<- DESeq(ds_txi)
res<- results(dds,contrast=list("isoA.groupT1","isoA.groupT2"), alpha= p_adjust_treshold,  lfcThreshold = L2FC_treshold)

But how can I get the general differences between treatment (group) T1 and T2 with elimination of the RNA-isolation date batch effect, if thats possible?

Could I maybe just do something like this:

res<- results(dds,contrast=list(c("isoA.groupT1","isoB.groupT1"),c("isoA.groupT2","isoB.groupT2")), alpha= p_adjust_treshold,  lfcThreshold = L2FC_treshold)

This is may session info as requested:

> sessionInfo()
R version 3.6.2 (2019-12-12)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    

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

other attached packages:
 [1] RColorBrewer_1.1-2          pheatmap_1.0.12             scatterplot3d_0.3-41       
 [4] edgeR_3.28.1                limma_3.42.2                ModCon_0.2.0               
 [7] data.table_1.12.8           readr_1.3.1                 vsn_3.54.0                 
[10] hexbin_1.28.1               DESeq2_1.26.0               SummarizedExperiment_1.16.1
[13] DelayedArray_0.12.3         BiocParallel_1.20.1         matrixStats_0.56.0         
[16] Biobase_2.46.0              GenomicRanges_1.38.0        GenomeInfoDb_1.22.1        
[19] IRanges_2.20.2              S4Vectors_0.24.4            BiocGenerics_0.32.0        
[22] rjson_0.2.20                tximport_1.14.2            

loaded via a namespace (and not attached):
 [1] bitops_1.0-6           bit64_0.9-7            tools_3.6.2            backports_1.1.7       
 [5] R6_2.4.1               affyio_1.56.0          rpart_4.1-15           Hmisc_4.4-0           
 [9] DBI_1.1.0              colorspace_1.4-1       nnet_7.3-14            tidyselect_1.1.0      
[13] gridExtra_2.3          bit_1.1-15.2           compiler_3.6.2         preprocessCore_1.48.0 
[17] htmlTable_1.13.3       scales_1.1.1           checkmate_2.0.0        genefilter_1.68.0     
[21] affy_1.64.0            stringr_1.4.0          digest_0.6.25          foreign_0.8-76        
[25] XVector_0.26.0         base64enc_0.1-3        jpeg_0.1-8.1           pkgconfig_2.0.3       
[29] htmltools_0.4.0        htmlwidgets_1.5.1      rlang_0.4.6            rstudioapi_0.11       
[33] RSQLite_2.2.0          farver_2.0.3           jsonlite_1.6.1         acepack_1.4.1         
[37] dplyr_0.8.5            RCurl_1.98-1.2         magrittr_1.5           GenomeInfoDbData_1.2.2
[41] Formula_1.2-3          Matrix_1.2-18          Rcpp_1.0.4.6           munsell_0.5.0         
[45] lifecycle_0.2.0        stringi_1.4.6          zlibbioc_1.32.0        grid_3.6.2            
[49] blob_1.2.1             crayon_1.3.4           lattice_0.20-41        splines_3.6.2         
[53] annotate_1.64.0        hms_0.5.3              locfit_1.5-9.4         knitr_1.28            
[57] pillar_1.4.4           geneplotter_1.64.0     XML_3.99-0.3           glue_1.4.1            
[61] latticeExtra_0.6-29    BiocManager_1.30.10    png_0.1-7              vctrs_0.3.0           
[65] gtable_0.3.0           purrr_0.3.4            assertthat_0.2.1       ggplot2_3.3.0         
[69] xfun_0.13              xtable_1.8-4           survival_3.1-12        tibble_3.0.1          
[73] AnnotationDbi_1.48.0   memoise_1.1.0          cluster_2.1.0          ellipsis_0.3.1

If I1 in group1 is the same as I1 in group 2 (you expect these to share a baseline effect), then you should just use a design of ~indiv + group.

However, if I1 in each group are not related to each other in any way, just arbitrary assignment of a number within a group, then you should use ~iso + group.

Nothing you said yet indicated the need for an interaction term with group.