Correct way to compare multiple treaments between RNA-Seq samples using edgeR?

Question

[Guille]

Hello,

I'm doing my masters on Bioinformatics and, as expected, I've hit a sort of roadblock of which I'm unsure how to proceed (my PI isn't very well versed on bioinformatics).

We are researching different treatments for a DSS-induced colitis model in mice. We usually go about this by quantifying gene expression by means of RT-qPCR in which logFC for each gene is calculated by averaging normalized expression vs the control group and comparing those values by running ANOVA or t-test on them.

With that in mind, the groups for this experiment (n=3) were as follows:

1. Control group.
2. DSS-induced colitis.
3. DSS + Treatment 1.
4. DSS + Treatment 2.
5. DSS + Treatment 1 + Treatment 2.

I'm using glmQLFTest to make the pertinent comparisons (and I've had no issue so far). My IP wants to find out whether adding Treatment 2 to Treatment 1 has a different enough response in colitis treatment.

Normally, for RT-qPCR data I'd just run an ANOVA test and compare every group to each other and find out which groups are different enough. However, I'm second-guessing myself on whether the comparisons I'm using are adequate.

Most of these comparisons are very straightforward when comparing to the Control or Disease groups, but when comparing group 5 to group 3 (as we want to know how treatment 2 can influence treatment 1), I'm getting fairly low p-values (that fade into nothingness after adjusting with topTags.

My code is as follows

contrast = makeContrasts(
  "DSS-CONTROL", #1
  "DSS.1-DSS", #2
  "DSS.2-DSS", #3
  "DSS.1.2-DSS", #4
  "DSS.1.2-DSS.1", #5
  levels=design)

res = glmQLFTest(fit, contrast = contrast[,5]) # I change this around to select a comparison
res_corrected = topTags(res, n = Inf)

Would this be the correct way to go about it? In this snippet, I only want to get the p-value for DEGs between the aforementioned groups 3 and 5, but the FDR is fairly high (it only goes as low as 0.10 for this single comparison).

Thank you for your time. I just might be an idiot and the differences caused by the treatment aren't biologically meaningful (which might spell trouble for me).

Answer 1

[Gordon Smyth]

The comparison you've done seems exactly right to test whether treatment 2 modifies or enhances the effect of treatment 1.

Have you examined an MDS plot to see whether the different groups appear separate?

[Guille]

Thank you for your answer!

The MDS plot seems to group most treatment groups close-by, so I'm starting to suspect the effect of treatment 2 may very well be masked by treatment 1 (we are not really testing for synergy, more like alleviating side effects of treatment 1). Biological constants like body weight and disease activity index do differ noticeably between groups, so it may be that changes are more subtle and conducting a bulk RNA-Seq experiment might be casting a far too wide web. I'll just consult my PI on whether using the unadjusted p-values may be valid (they really get squashed to smithereens but some of the top DEGs are ones we have confirmed to be diferentially expressed by RT-qPCR).