Correctly set contrast with DESeq2 for not full rank design matrix
1
0
Entering edit mode
llooll1 • 0
@llooll1-23715
Last seen 4.1 years ago

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
DESeq2 DGE contrast • 1.5k views
ADD COMMENT
0
Entering edit mode

We discourage cross-posting to biostars and BioC without explicitly linking the posts:

https://www.biostars.org/p/444660/

A lot of the advice you had already received on biostars is appropriate here.

ADD REPLY
0
Entering edit mode

Hi, unfortunately there was no advice for correcting the bias, thats why I tried again with a better phrasing of the question to avoid confusion.

ADD REPLY
1
Entering edit mode

Thanks for calling a one-day endeavour involving two experienced users who tried to provide you with help "no advise".

ADD REPLY
0
Entering edit mode

I appreciate the answers very much, but the question was closed on biostars by you for "not fitting the main topic of this site", before I knew how to proceed with my analysis. So I asked the question again slightly rephrased here. Now I know, that I should just ask here for futur similar topics. Sorry, I didn't want to offend you.

ADD REPLY
0
Entering edit mode

No, I closed the one on Boostars after this one was posted here to avoid that even more users invest double-effort. The ...does not fit the main topic... is a phrase that is automatically being added when a topic is closed. That all is spilled milk under the bridge now so I suggest to forget about it and proceed with our daily duties. You are always welcome to post at any community you like but it is good practice to avoid crossposting.

ADD REPLY
3
Entering edit mode
@mikelove
Last seen 19 hours ago
United States

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.

ADD COMMENT
0
Entering edit mode

Thank you for your answer. I1 are samples from the same individual in every group.

However, the RNA-isolation for all samples of 2 individuals was done seperately from the RNA-isolation for samples of the other 2 individuals and someone pointed out, that this is a bias which should better be corrected, if possible.

Can I still use indiv + group?

ADD REPLY
1
Entering edit mode

Yes still use indiv + group. It will correct for A/B differences as part of the indiv baseline. If you have further questions I’d recommend consulting with a statistician.

ADD REPLY
0
Entering edit mode

Ok, thank you very much for the help and for the tool!

ADD REPLY

Login before adding your answer.

Traffic: 698 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