Hi Mike and community,

I’m using theDESeq2 package for my RNA-seq analysis. I would like to get feedback on the choice of my models and applied contrasts. If there is a better model to answer our questions, I would appreciate any suggestions. 

I use a 3-level factor analysis:

Cultivar factor with 2 levels (cult1, cult2)

Time factor with 3 levels (T1, T2, T3)

Treatment factor with 2 levels (Infected, Control) with 3 biological replicates each.

Aim:

1) Effect of treatment on cultivar

2) Effect of time on cultivar 

Complication: In order to compare the effect of treatment within and across cultivars, we have to reduce infected samples with control samples before making any comparison for each time point and cultivar. Therefore, I have nested cultivar with time point for my first question, and run full model using Wald test where I have used only interaction terms to compare effect of treatment within and between time points for each cultivar.

I have been wondering if using main effect + interaction term could give me the same control over treatment within each cultivar and time point without nesting?

1) Effect of treatment on cultivar:

My reference levels are in alphabetical order. cult1, Control. 

full~Cultivar + Cultivar:Treatment

In order to compare across cultivars while controlling for specific time points, I applied this contrast:

contrast = list(c("Cultivarcult1_t1.TreatmentInf","Cultivarcult2_t1.TreatmentInf")), test = 'Wald', alpha = 0.01)

2) Effect of time on cultivar:

I want to detect genes that are differentially expressed due to time within cultivar and across cultivars, I didn’t nest any variables in the model, I have followed this post: https://support.bioconductor.org/p/101002/.  

full ~Cultivar + Treatment + Time + Cultivar:Time

reduced ~Cultivar + Treatment + Time, test="LRT"

Thank you!

Iryna

1. Changes were made, ignore this submission.

hi Iryna

I think you can model everything in one go, and extract contrasts afterward. I tried to mock up your experiment (minus the replication) with this code:

cond <- gl(2,3,12) # control vs treated
cltv <- gl(2,6) # cultivar 1 or 2
time <- gl(3,1,12) # three time points

Then the following design, ~cltv + cltv:time + cond:cltv:time, gives you a cultivar difference at time point, differences for each cultivar at each time point, and a potential effect of treatment for each cultivar and time point. This sounds like the correct design for your experiment. You can test individual effects in this design using a Wald test, or if you want to test multiple coefficients, you can remove the columns from the design matrix and use a LRT. Are you familiar with how to build the model matrix:

mm <- model.matrix(~cltv + cltv:time + cond:cltv:time, data=colData(dds))

If you wanted to do the LRT, you can provide mm to DESeq() as the full argument, and then remove specific columns from mm and provide this reduced matrix to the reduced argument.