DESeq2 test for synergistic vs additive effects
Entering edit mode
Jane • 0
Last seen 7 months ago
United States

I am trying to identify synergistic vs additive effects of various treatments on gene expression using DESeq2. I have data from 3 animals that were each treated with hormone A, hormone B, or their combination. Each animal received each of the treatments, so my original model design includes blocking by individual:

    dds <- DESeqDataSetFromTximport(txi, 
                                colData = samples, 
                                design = ~ Animal + Treatment)

    dds$Treatment <- relevel(dds$Treatment, ref = "control")
    dds <- DESeq(dds)
     [1] "Intercept"                 "Subject_2_vs_1"          "Subject_3_vs_1"          "Treatment_A_vs_control"     
     [5] "Treatment_B_vs_control"      "Treatment_AB_vs_control"  

resA <- results(dds, contrast = c("Treatment", "A", "control"))
resB <- results(dds, contrast = c("Treatment", "B", "control"))
resAB <- results(dds, contrast = c("Treatment", "AB", "control"))


sample            animal        treatment
con_rep1          1             control
con_rep2          2             control
cont_rep3         3             control
hormA_rep1        1             hormone A
hormA_rep2        2             hormone A
hormA_rep3        3             hormone A
hormB_rep1        1             hormone B
hormB_rep2        2             hormone B
hormB_rep3        3             hormone B
hormAB_rep1       1            hormone A + hormone B
hormAB_rep2       2            hormone A + hormone B
hormAB_rep3       3            hormone A + hormone B

I would like to test for synergy between treatments A and B (i.e., whether AB - A - B = 0 or not). Can someone please suggest a model design for this, as well as how to extract genes showing synergistic rather than additive effects using the results function?

I was trying to use the results function to do this as suggested in other posts:

res_syn <- results(dds, contrast = c("Treatment_AB_vs_control"  - "Treatment_A_vs_control" - "Treatment_B_vs_control"))

but it returns the error

Error in "Treatment_AB_vs_control" - "Treatment_A_vs_control" : 
  non-numeric argument to binary operator
synergistic additive DESeq2 • 912 views
Entering edit mode

Have you already checked the section on Interactions in the DESeq2 vignette, ?

Entering edit mode

Thanks, Wolfgang.

I saw that but thought that it referred to interactions between factors (e.g., in my case Animal*Treatment), but is there a way to test for interactions between specific levels of a factor?

Entering edit mode

Jane, you can (perhaps should) think of treatment with hormone A and hormone B as two different (binary: yes or no) factors, and then the concept of interactions between factors applies - just follow the steps I suggested previously.

I am not aware that "interactions between specific levels of a factor" is a thing.

Finally, note that the concepts of interactions is tricky in the context of drug or hormone treatments, and the biomedical intention or intuition of an interaction does not always directly map 1:1 to that of a statistical interaction. Namely, drug or hormone response phenotypes are often non-linear (e.g. sigmoidal) and thus, if you define a (statistical) interaction as "effect of A+B" is different from "effect of A" + "effect of B", then a drug can easily interact with itself (A=B): "effect of 2 x A" is different from 2 x "effect of A".

Entering edit mode
Last seen 11 hours ago
United States

No idea if this is possible using DESeq2, but it would be trivial to do so using edgeR. That said, algebraically what you are doing doesn't make sense. Note that each of your coefficients are already differences between a given treatment and control, so you are doing (here I substitute TrtA, etc for your long names)

(TrtAB - Control) - (TrtA - Control) - (TrtB - Control)
##cancel two of the controls and now
TrtAB - TrtA - TrtB + Control

## which isn't really a thing

But you have already blocked on animal, so the controls are extraneous anyway. So you could do

dglst <- DGEList(txi, genes = <add annotation data here>)
dglst <- calcNormFactors(dglst)
samples$trt <- factor(rep(c("Control","TrtA","TrtB","TrtAB"), each = 3), c("Control","TrtA","TrtB","TrtAB")
design <- model.matrix(~0 + Animal + trt, samples)
colnames(design) <- gsub("trt", "", colnames(design))
contrast <- makeContrasts(TrtAB - (TrtA + TrtB)/2, levels = design)

I leave the remainder of the analysis pipeline up to you to do.

Note that TrtAB - TrtA - TtrB isn't quite correct either, which is why I compared the double treatment to the mean of the two single treatments. As an numerical example, let's say there actually is a synergistic effect, and TrtAB increases a gene's exposure by 100% more than either TrtA or TrtB alone. For simplicity, let's assume no difference between TrtA and TrtB (and the logCPM for each is 10). That means TrtAB will have a logCPM of 11. So TrtAB - TrtA - TrtB = -9, which is obviously not what you want. Whereas TrtAB - (TrtA + TrtB)/2 = 1, which is what we would expect if it's synergistic.

Entering edit mode

In terms of design matrices and contrasts, I think the capabilities of DESeq2 and edgeR are pretty much equivalent. Differences exist on "details" of how things are fitted, and of course the interface.

To me, the real bottleneck here seems to be the formulation of the question / what Jane actually wants to achieve (also see comment above).

Entering edit mode

Why is TrtAB - TrtA - TrtB + Control "not really a thing"? i came to this forum answer when attempting to use EdgeR to answer a similar question, and this chapter seems to suggest TrtAB - TrtA - TrtB + Control (amongst other, non-averaging contrasts) is appropriate for this analysis. I'm just having trouble reconciling the two. I follow your numeric counterexample to TrtAB - TrtA - TtrB, which is also suggested in the reference in section 1.9.

Edit: Here is a reference that is germane to differential expression. It also raises the TrtAB - TrtA - TrtB + Control as the proper contrast.


Login before adding your answer.

Traffic: 380 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6