R formula in RNA-Seq differential analysis with multiple parameters

Question

[Thibault D.]

Dear all,

I'm working on differential genes/transcripts analysis in which I am interested in the effect of a treatment between patients, and their reaction across two age groups. I have a large number of patients in that study.

Here, both age and treatment factors are interesting. However, I am experiencing doubts about the statistical formula to be written withing R (for instance in DESeq2 or in Sleuth).

For now, I have counted reads over the transcriptome and performed classic clustered heatmap and PCA. They are showing the effect of the age is the most powefull effect over the experimental design I am analyzing.

For the differential expression analysis, which of the the following formulas should I consider :

1. ~treatment+age (treatment is considered a batch effect, according to DESeq2 documentation. I don't feel like this is correct here.)
2. ~age (forget the treatment)
3. ~treatment:age (consider interactions between age and treatment)
4. ~treatment*age (consider both effect and interaction of age and treatment, it is similar to ~treatment+age+treatment:age)

I am a bit confused, and I have always been confronted to the classic "condition-and-batch-effect" type of experimental design.

I have tried all of those formulas, and I have some genes/transcripts in common, other are diverging, adjusted p-values and q-values are sometimes very diverging. Most of my fellow biologist coworkers' hypothesis are answered by those formulas, but I guess that only one is correct.

Thanks!

Answer 1

[Devon Ryan]

1. treatment is not so much a batch effect as a "main" effect in the normal parlance (the same as age). This is one of the two possible model choices for you.
2. Egad don't do this!
3. Egad don't do this!
4. This is the other option (it is identical to ~treatment + age + treatment:age)

I wouldn't be surprised if there's an interaction between age and the size of the effect, so my knee-jerk reaction is to favor model 4 over model 1. I presume that the biologists also presumed this could happen, since otherwise there wouldn't have been multiple ages.