Question

DESeq2 test for synergistic vs additive effects

1

Entering edit mode

Jane ▴ 10

@jkhudyakov-23010

Last seen 18 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)
    resultsNames(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"))

colData:

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 • 1.9k views

ADD COMMENT • link updated 8 months ago by Wolfgang Huber ★ 13k • written 18 months ago by Jane ▴ 10

0

Entering edit mode

Have you already checked the section on Interactions in the DESeq2 vignette, http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions ?

ADD REPLY • link 18 months ago Wolfgang Huber ★ 13k

0

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?

ADD REPLY • link 18 months ago Jane ▴ 10

0

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

ADD REPLY • link 17 months ago Wolfgang Huber ★ 13k

0

Entering edit mode

Hi, Wolfgang.

I have a similar experiment design, so when looking for help in the design matrix, I came across this post which I found really useful. In my case I have data from a cell line with two treatments, A and B, the combined one, AB, and a Control treatment. I also have a considerable batch effect between replicates (two replicates).

I have read the DESeq2 vignette as well as several articles explaining this type of designs. They helped me a lot but now I have some questions.

Originally I went for this design: ~ batch + treatment. In the treatment factor I separated the four treatments as CTRL, A, B, AB. This allowed me to test all the comparisons (A vs CTRL, B vs CTRL, AB vs CTRL, B vs A, AB vs A, AB vs B), while also accounting for the batch effect.

Now, I am rethinking this analysis including an interaction between the treatments. For that I created dummy variables for the two treatments A and B. So now, my design looks like this: ~ batch + A + B + A:B. This design allows me to test for the interaction and the main effects while accounting for the batch effect. Using the contrasts I can test:

A vs CTRL, using the contrast=c("A_YES_vs_NO") because in the main effect of A, B=NO
B vs CTRL, using the contrast=c("B_YES_vs_NO") because in the main effect of B, A=NO
AB vs A, using the contrast=list(c("B_YES_vs_NO", "BYES.AYES"))
AB vs B, using the contrast=list(c("A_YES_vs_NO", "BYES.AYES"))

But I can't figure out if the following contrasts can be done:

A=YES & B=YES vs A=NO and B=NO, which would be something like doing AB vs CTRL in my original design.
B=YES vs A=YES, which would be something like doing A vs B in my original design.

1) Can they be done? How?

2) Also, is the batch effect really corrected in the contrasts mentioned before? I mean, because I didn't include the interaction between the treatments and the batch, I assume the main effects of the tratments and the interaction have been corrected by the batch. Is this rigth? Or the main effects of the tratments and the interaction are relative to the reference level of the batch?

Last but not least, I found really interesting what you said in the previous comment about biomedical and statistical interaction. My model is a breast cancer human cell line with two different hormones. So, taking what you said into account...

3) Is this approach correct? Or is there a better way of analyzing RNA-seq data to test for additive/synergy/antagonism effects? What do you reccomend?

Thank you in advance!

ADD REPLY • link 8 months ago eve.olszanowski • 0

0

Entering edit mode

Dear Eve

Thank you for your detailed questions.

(1) Yes. Any contrast can be computed, these are just linear combinations of estimated effects. See https://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts and perhaps also https://github.com/tavareshugo/tutorial_DESeq2_contrasts

(2) The way you describe the data, I would, in addition to ~ batch + A + B + A:B (which is the same as ~ batch + A B) also do just ~ A B separately for each replicate (batch) and then compare the results manually, e.g. plot them against each other. This could help in particular if one of the two batches catastrophically failed. Vice versa, seeing good correlation between estimated effects in the two replicates is a useful quality assessment.

(3) I'm tempted to say "both". The approach is correct. But I also suspect that this is an exploratory analysis, i.e. one that is followed up by in-depth investigation of findings and follow-up experiments - rather than a definite, final result and the very last step in your project, I would put emphasis on discovery, priorisation and visualisation, rather than on inference ("p values").

In particular, since we are working essentially on a log scale, "no interaction" in the above model means that the effects of A and B on a gene's expression [as measured in something like mRNA abundance or concentration] are multiplicative. This is mathematically convenient, but perhaps there are also biological situations in which addition would be the more natural "no interaction" concept.

ADD REPLY • link 8 months ago Wolfgang Huber ★ 13k

score 0 · Answer 1 · 2023-07-03

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.