Hej

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!

Julia

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).