Question: model matrix DESeq2 1.10
gravatar for juliah
17 months ago by
juliah0 wrote:


I have a question concerning the use of interaction designs with DESeq2.

My experimental design looks like this:

3 treatments: control, low, high

4 sampling time points throughout a years season: 1,2,3,4

3 plots per treatment: control (a,b,c), nf (d,e,f), fert (g,h,i)

Each plot was sampled at each sampling time point, resulting in 36 samples in total.

It is 16S amplicon sequencing data, so instead of a RNA-Seq count-table i have a otu-table.

PCoA plots indicate that there is a separation of the samples by treatment and possibly a weak time effect overall samples (not specificly for one treatment). Additionally i have a strong plot effect: the 4 samples from each plot are closer to each other than to another sample of the same time from the same treatment. The plots are geographically close but they are random spread out over an experimental site, so not every time a control plot lies next to a low and a high plot-

One question is:can i implement this confounding plot factor in my design and how would that work?

i want to address these questions and find Differentially abundant otus for:

1. which are the significant otus causing a treatment effect?

2. is there a season effect and what are the otus driving the seasonal variation?

3. is there one treatment behaving different in its seasonal response than the others?

4. which otus are causing this seasonal behaviour in this treatment?

For the first question i would simply make a design ~treatment and test the three comparisons (control vs low), (control vs high) and (low vs high)- ignoring seasonal influences. or including seasonal influences with ~date+treatment but still test the same contrasts =(treatment, control, low)...?

For the second part i can again make a simple design ~date and make the three comparisons (1 vs 2), (1vs3),(1vs4) or do this with the design ~ treatment+date?

For question 3 i need a design with interactions ~treatment+date+treatment:date, with DESeq1.6 or 1.8 i would have received this results(Names)

> dds <- DESeqDataSetFromMatrix(
+ countData = data,
+ colData = data.frame(condition=conditions,
+ date=date,
+ treatment=treatment,
+ soilplot=SoilPlot),
+ design = ~ treatment+date+treatment:date)
> dds <- estimateSizeFactors(dds)
> dds <- estimateDispersions(dds)

> dds <- DESeq(dds)

> resultsNames(dds)
 [1] "Intercept"              "treatmentcontrol"         
 [3] "treatmenthigh"          "treatmentlow"           
 [5] "date1"                "date2"               
 [7] "date3"                "date4"              
 [9] "treatmentcontrol.date1"  "treatmenthigh.date1"
[11] "treatmentlow.date1"    "treatmentcontrol.date2"
[13] "treatmenthigh.date2"  "treatmentlow.date2"   
[15] "treatmentcontrol.date3"  "treatmenthigh.date3"
[17] "treatmentlow.date3"    "treatmentcontrol.date4"
[19] "treatmenthigh.date4" "treatmentlow.date4" 

and would have tested the interaction contrast=list("treatmenthigh.date1","treatmenthigh.date2")... to test if treatment high is behaving differently from control and low in its seasonal response. but with version 1.10 i dont know how to achieve the same as the same design results in this resultsNames(dds)

[1] "Intercept"              "treatment_high_vs_control"
 [3] "treatment_low_vs_control"   "date_2_vs_1"       
 [5] "date_3_vs_1"        "date_4_vs_1"      
 [7] "treatmenthigh.date2"  "treatmentlow.date2"   
 [9] "treatmenthigh.date3"  "treatmentlow.date3"   
[11] "treatmenthigh.date4" "treatmentlow.date4" 

is there a way that i can still receive the full model matrix also with betaprior FALSE settings for interaction designs?
or which contrast would i now set to achieve the previous comparison?

For the last question testing the seasonal effect only within f.e. treatment high, with the previous version 1.8 i would have formed this contrast=list(c(date1,treatmenthigh.date1),c(date2,treatmenthigh.date2)) and so on but again with 1.10 i dont know how to achieve this? otherwise i would separate the treatment high samples and only check ~date contrasts or form groups<-factor(dds$treatment,dds$date) and use the design ~group to make this contrast=c(group,high1,high2), this is still working the same.

thanks a lot

please tell me if i am thinking the wrong way and my testing doesnt make sense!






ADD COMMENTlink modified 17 months ago • written 17 months ago by juliah0
gravatar for Michael Love
17 months ago by
Michael Love12k
United States
Michael Love12k wrote:

hi Julia,

Can you take a look at this answer I gave to a someone else about how to reformulate interactions designs in version 1.10 (without the expanded model matrices)?

A: Contrast All Interaction Effects in DESeq2 in current version (1.10.1)

In your case this would mean to remove the 'treatment' term and use '~date + date:treatment' to produce date-specific treatment effects for all dates (note that this gives the same model, just reformulated terms). You can then easily contrast these using the list style of 'contrast'. 

In general, to answer all of your questions, I think you should get help from a local statistician. There are a lot of issues here, and resolving them all is not so much a software issue, but key analysis decisions. With the interaction model, the terms which are appear will be identical to the ones you would get with standard linear regression in R.

The confounding plot issue is tricky, because there is no way to differentiate the individual plot effects and the treatment (even though you have them over time). Using limma, I believe it is possible to inform the model that samples within a single plot are correlated over time using duplicateCorrelation(), but with DESeq2 you can't use fixed effects to account for effects which are confounded with a comparison of interest (treatment).

ADD COMMENTlink written 17 months ago by Michael Love12k
gravatar for juliah
17 months ago by
juliah0 wrote:

Hej Michael

Thanks for your immediate reply!!!

I had seen the related post but had set up the design wrong ~treatment+treatment:date. So now i made it the way you suggested ~date+date:treatment and do get these resultsNames():

[1] "Intercept"              "date_2_vs_1"       
 [3] "date_3_vs_1"        "date_4_vs_1"      
 [5] "date1.treatmenthigh"  "date2.treatmenthigh" 
 [7] "date3.treatmenthigh"  "date4.treatmenthigh"
 [9] "date1.treatmentlow"    "date2.treatmentlow"   
[11] "date3.treatmentlow"    "date4.treatmentlow"  

Do i understand right now with list style to answer question three(is treatment high different in seasonal behaviour from control and treatment low)  it would be contrast=list(c("date1.treatmenthigh","date1.treatmentlow"),c("date2.treatmenthigh","date2.treatmentlow"))?with control as reference level. in the other post you mention sth with an average effect and 1/3- how would that look like?

is it possible in version 1.10 to use custom model matrices?

Also for correct understanding: If you have a design ~treatment+date the effect of treatment will be blocked/removed and only date will be considered?

Thanks for the advice with limma, i will try to use this for the plot effect if possible.

all the best


ADD COMMENTlink written 17 months ago by juliah0

That contrast doesn't make sense to me.

I can tell you what the terms in resultsNames(dds) mean with some example results tables, but I think you need to partner with someone if you want to explore testing of more hypotheses.

the treatment high vs control effect in date1:

results(dds, name="date1.treatmenthigh")

the treatment high vs low effect in date1:

results(dds, contrast=list("date1.treatmenthigh","date1.treatmentlow"))

is the high vs control effect different in date 2 vs date 1:

results(dds, contrast=list("date2.treatmenthigh","date1.treatmenthigh")) 
ADD REPLYlink written 17 months ago by Michael Love12k
gravatar for juliah
17 months ago by
juliah0 wrote:

Thanks again.

I was trying to formulate this contrast:


as it was possible using the interaction design in v 1.8, giving me the interaction term for date2 vs date1 in treatment high.

if this term is non-zero, then treatment high has a different date effect than the main effect for date

but couldnt see how to get it with the same interaction design in 1.10 as in my case now always treatment control is used as reference level also for the main effect of date. Using the last contrast you suggested will give the interaction term for date2 vs date1 in treatment high vs control. this tests if the date effect is different in treatment high compared to treatment control.

This is different right?




ADD COMMENTlink written 17 months ago by juliah0

It's true that this contrast is the interaction term, not the total effect of date 2 vs date 1 in treatment high. If you look over the current DESeq2 vignette, I've added a longer section on how interaction terms work. The important thing I diagram here is that, interaction terms for a particular group (here, treatment high) plus the main effect give the total effect in that group. So to compare date 2 vs date 1 in treatment high (not whether high vs control is different between these dates) you add the main effect and the interaction contrast. This looks like:

results(dds, contrast=list(c("date_2_vs_1","treatmenthigh.date2"), "treatmenthigh.date1"))
ADD REPLYlink written 17 months ago by Michael Love12k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 122 users visited in the last hour