Three factor model matrix, contrasts with interactions - DESeq2 / edgeR
1
1
Entering edit mode
Alan Smith ▴ 150
@alan-smith-5987
Last seen 7.5 years ago
United States

Hello,
I'm working on differential expression analysis with the given sample data using DESeq2.

Strain=rep(rep(c("ST1","BT6","MT3","XT7"),each=3),4)
Treatment=rep(c("Low","High"),each=24,1)
Tissue=rep(c("Lvr","Hrt"),each=12,2)
targets<-data.frame(Tissue,Treatment,Strain)
rownames(targets)<-paste0("samp",c(1:nrow(targets)))
strain <- factor(targets$Strain, levels=c("ST1","BT6","MT3","XT7"))
treatment <- factor(targets$Treatment, levels=c("Low","High"))
tissue<-factor(targets$Tissue,levels=c("Lvr","Hrt"))

I'm using the following design:

designFormula= "~tissue+treatment+strain+tissue:treatment+tissue:strain+treatment:strain"
tmp<-model.matrix(as.formula(designFormula),targets)
colnames(tmp)

 [1] "(Intercept)"             "tissueHrt"
 [3] "treatmentHigh"           "strainBT6"
 [5] "strainMT3"               "strainXT7"
 [7] "tissueHrt:treatmentHigh" "tissueHrt:strainBT6"
 [9] "tissueHrt:strainMT3"     "tissueHrt:strainXT7"
[11] "treatmentHigh:strainBT6" "treatmentHigh:strainMT3"
[13] "treatmentHigh:strainXT7"

Questions

I) Is my design correct for addressing my questions below?

II) Could someone please confirm if I'm interpreting these coefficients correctly:
        2) tissueHrt = finds genes that are differentially expressed in Hrt vs Lvr comparison in strain ST1 and Low heat treatment?
                        OR does that give the genes that are DE in Hrt vs Lvr across all strains and treatments
        4) strainBT6 = genes DE between BT6 and ST1 across in only reference tissue and treatment or all
        7) tissueHrt:treatmentHigh = Interaction effect of High vs Low treatment in Htr tissue (OR Genes that respond to High heat treatment in Hrt tissue)
        8) tissueHrt:strainBT6 = Interaction effect of that strain BT6 has on Hrt tissue rather than on Lvr
        9) treatmentHigh:strainBT6 = Interaction effect that genotype BT6 has on treatment High heat rather than Low heat

III) Does this results(dds,contrast=list(c("tissueHrt","tissueHrt:strainBT6")),parallel=TRUE)
provide Htr tissue Vs. Lvr effect in Strain BT6?

IV) How to find genes that respond to the High heat treatment in any strain in DESeq2 (overall high heat treatment effect on all strains)
i.e., the equivalent of glmLRT(fit, coef=11:13) in edgeR?

VI) What will be the contrast that finds genes that respond differently to the heat treatment in BT6 vs MT3 strains?

V) What are the contrasts that find the overall interaction effect of the tissue and strain?

Thanks for your time.

Greatly appreciate your help.


Alan

 

sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X El Capitan 10.11.6

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] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] edgeR_3.16.5               limma_3.30.13              DESeq2_1.14.1              SummarizedExperiment_1.4.0 Biobase_2.34.0             GenomicRanges_1.26.4       GenomeInfoDb_1.10.3       
 [8] IRanges_2.8.2              S4Vectors_0.12.2           BiocGenerics_0.20.0       

loaded via a namespace (and not attached):
 [1] genefilter_1.56.0    locfit_1.5-9.1       splines_3.3.2        lattice_0.20-35      colorspace_1.3-2     htmltools_0.3.5      yaml_2.1.14          base64enc_0.1-3      survival_2.41-3     
[10] XML_3.98-1.6         foreign_0.8-67       DBI_0.6-1            BiocParallel_1.8.1   RColorBrewer_1.1-2   plyr_1.8.4           stringr_1.2.0        zlibbioc_1.20.0      munsell_0.4.3       
[19] gtable_0.2.0         htmlwidgets_0.8      evaluate_0.10        memoise_1.0.0        latticeExtra_0.6-28  knitr_1.15.1         geneplotter_1.52.0   AnnotationDbi_1.36.2 htmlTable_1.9       
[28] Rcpp_0.12.10         acepack_1.4.1        xtable_1.8-2         scales_0.4.1         backports_1.0.5      checkmate_1.8.2      Hmisc_4.0-2          annotate_1.52.1      XVector_0.14.1      
[37] gridExtra_2.2.1      BiocStyle_2.2.1      ggplot2_2.2.1        digest_0.6.12        stringi_1.1.3        grid_3.3.2           rprojroot_1.2        tools_3.3.2          bitops_1.0-6        
[46] magrittr_1.5         lazyeval_0.2.0       RCurl_1.95-4.8       tibble_1.3.0         RSQLite_1.1-2        Formula_1.2-1        cluster_2.0.6        Matrix_1.2-8         data.table_1.10.0   
[55] rmarkdown_1.4        rpart_4.1-10         nnet_7.3-12   
deseq2 edgeR multifactorial design limma voom • 1.8k views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 1 day ago
United States

hi Alan,

If removing a term from the design formula would remove all the columns for which you wish to perform an LRT, then you can just use the 'reduced' argument with that reduced design formula. Otherwise, the equivalent way to perform a likelihood ratio test of certain coefficients in the design matrix is to provide a 'full' and 'reduced' design matrix to the DESeq() function while specifying test="LRT". If you have your full model matrix:

full <- model.matrix(design, colData(dds))

Then you remove certain columns:

reduced <- full[, -idx]

Then:

dds <- DESeq(dds, full=full, reduced=reduced, test="LRT")
res <- results(dds)

You have a number of contrasts where you want to compare two groups, e.g. heart vs liver in ST1 and low treatment. You can follow the example in the DESeq2 vignette and just use a design of ~group which combines all the variables into one, and then use the 'contrast' argument to compare these two groups. This will make the comparisons much easier for you.

Otherwise, you should partner with a statistician who can help you to identify the interpretation of all the terms.

ADD COMMENT
0
Entering edit mode

Hi Mike,

I initially started with the "group" method to compare those individual contrasts. But in order to find the interaction effect of the strain and treatment, I attempted the 3 factor design. From the DESeq2 manual and workflows, I "understood" the main effects (which are specific to the reference level) and interaction effects in 2 factor design. As you can tell, I'm totally confused about the 3 factor design.

If I'm not asking too much, could you please answer two questions before I find and pair up with a statistician.

In this design

designFormula= "~tissue+treatment+strain+tissue:treatment+tissue:strain+treatment:strain"

1) Does coefficient "tissueHrt" finds genes that are differentially expressed in Hrt vs Lvr comparison in strain ST1 and Low heat treatment?

2) Does "treatmentHigh:strainBT6" finds the additional effect of high vs low heat treatment in genotype BT6 rather than in ST1 (the reference)

Thanks a lot for your help.

 

ADD REPLY
0
Entering edit mode

Got it. 

1) Yes

2) Yes

ADD REPLY

Login before adding your answer.

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