Hello,
I have 5 treatments with 3 biological replicates and 2 technical replicates for each one. My design looks like:
- DMEM (Basal condition)
- Vehicle (Basal condition + vehicles used for treatments 3-5)
- Drug
- Insult
drug_insult
After collapsing the technical reps my
sampleTablelooks like:
condition Replicates sex
Drug 1 male
DMEM 1 male
drug_insult 1 male
Insult 1 male
vehiculo 1 male
Drug 2 female
DMEM 2 female
drug_insult 2 female
Insult 2 female
vehiculo 2 female
Drug 3 female
DMEM 3 female
drug_insult 3 female
Insult 3 female
vehiculo 3 female
DMEM is the control only for Drug and Vehiculo is the control for treatments 3-5. My question is: how to account for these two controls in the formula for DESeq?
I normally would run
x <- DESeqDataSetFromTximport(tximport_object, sampleTable, ~ sex + condition + sex:condition)
x$condition <- relevel(x$condition, "vehiculo")
dds <- DESeq(x)
And then extract the contrast of interest using results(). For instance, if I want to compare Drug vs DMEM I would run:
res <- results(dds, contrast = c("condition", "drug", "DMEM"))
However, the control for that contrast, as specified with the relevel() call is actually Vehiculo. So, I presume this might be giving me different genes that it would otherwise give me if I had specified DMEM as the control for that contrast with:
x <- DESeqDataSetFromTximport(tximport_object, sampleTable, ~ sex + condition + sex:condition)
x$condition <- relevel(x$condition, "DMEM")
dds <- DESeq(x)
I have also checked other threads (1, 2), from which the closest one is 2, however, it is still not clear for me how could I encode the two controls in my formula. So any help will be appreciated,
Best,
Tain.
> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
>packageVersion("DESeq2")
[1] ‘1.24.0’
Hello @swbarnes2,
I am already using the
results()method to extract the contrast of interest, however, my concern is that I might be comparing, for instance, Drug vs DMEM using as control Vehicle (as this was the first level I specified in the formula) when in reality it must be Tibolone vs DMEM (using DMEM as the first level in the formula). The concern traces back to the vignette:Though I have 2 controls.
Best,
Tain.
Sorry, I misspoke in my answer. But your question is incomplete without saying exactly what comparison you are asking for in results, because the default is going to be to give you the genes which respond differently to treatment between sex, and it's not at all clear to me that that's what you really want.
Thank you for your suggestion. I have updated the question to state my problem more clearly.
Okay, so you observed this:
So you know that your setup is only looking at the difference between Drug and DMEM in the reference sex?
Yes, I noted that. I have no problem extracting the sex-specific differences using the
~groupdesign but instead trying to extract the general (i.e. no sex-specific) results for Drug vs DMEM.I think that what I am going to do is: 1. change my design to
~ sex + conditionset Vehiculo as control, 2. extract all comparisons except the Drug vs DMEM, 3. set DMEM as control and extract the Drug vs DMEM, 4. finally, extract the sex-specific differences by using the design~group.Would you advise such approach?
Thank you for your answers.