So recently we have got sample back with the following design
| Condition | Treatment | Lane | Batch | RIN |
| Case | Treated | 1 | 1 | 7.7 |
| Case | Treated | 2 | 1 | 7.7 |
| Case | Treated | 1 | 4 | 7.6 |
| Control | Treated | 2 | 4 | 8.1 |
| Control | Treated | 1 | 1 | 7.8 |
| Control | Treated | 2 | 1 | 7.7 |
| Case | Untreated | 1 | 4 | 7.5 |
| Case | Untreated | 2 | 3 | 7.8 |
| Case | Untreated | 2 | 2 | 7.9 |
| Control | Untreated | 2 | 1 | 7.4 |
| Control | Untreated | 1 | 3 | 8 |
| Control | Untreated | 1 | 2 | 7.8 |
data.frame(Condition=rep(rep(c("Case","Control"), each=3),2), Treatment=rep(c("Treated", "Untreated"),each=6), Lane=c(1,2,1,2,1,2,1,2,2,2,1,1), Batch=c(1,1,4,4,1,1,4,3,2,1,3,2), RIN=c(7.7,7.7,7.6,8.1,7.8,7.7,7.5,7.8,7.9,7.4,8,7.8))
What we would like to do is to:
1. Compare effect of treatment in case
2. Compare effect of treatment in control
3. Compare and contrast the effect of treatment in case when compared to control.
From previous answers, it seems like there are two different model that we can use:
1. ~Batch+Lane+RIN+Condition+Condition:Treatment
2. ~Batch+Lane+RIN+Condition+Group
Where group is <Treatment><Condition>
Take for example, when we try to compare the difference of condition when there is no treatment
e.g. Untreated Case vs Untreated Control
For 1, we use
dds = DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design = ~Batch+RIN+Lane+Condition+Condition:Treatment)
dds$Condition = relevel(dds$Condition, ref="Control")
dds$Treatment = relevel(dds$Treatment, ref="Untreated")
dds = DESeq(dds)
results(dds, contrast=c("Condition", "Control", "Case"))
For 2 we use
coldata$Group = paste0(coldata$Treatment, "-", coldata$Condition)
dds = DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design = ~Batch+RIN+Lane+Group)
dds = DESeq(dds)
results(dds, contrast=c("Group", "Untreated-Control", "Untreated-Case"))
However, I found that the two design actually give very different results.
Using the interaction term model, we will get around 438 genes with padj < 0.05 whereas none of the genes can anywhere close to significance when we use the group model.
Is there something that I did wrongly? Which one is the more appropriate method? And most importantly, why was there such a large difference?
Thank you Michael,
I am currently using DESeq2 version 1.10.0. I was adding the the RIN just to see if the RNA degradation might affect the results. I have tried to get the correlation of the p-value (which might not be correct but was only try to get an idea) of all possible combinations and the results from model 1 and model 2 doesn't agrees (the correlation file is here). So I am uncertain whether if there is some problem with my scripts...
(Forgot to mention, the rows with Interact as prefix are those that uses the interaction model).
Below is my R sessionInfo():
Repeat the analysis but definitely do not include RIN in the design. Then compare 1 and 2.
It is slightly better in that both doesn't return any significance, however, the p-value still differs in that their correlation is only around 0.6827998
Really thank you for your help Michael
Correlation on p-values isn't a great comparison. Better is a cross tabulation:
table(res1$padj < .1, res2$padj < .1)
Remember though, that they are not identical methods so the results are not identical (2 uses LFC shrinkage and the 1 does not). My advice is if you are interested in testing the interaction term use 1, if you are interested in comparing groups use 2.
Thank you Michael. None of them return any significance so I tried to use pvalue < 0.05 instead
If that is normal, maybe I will stick with model 2 just because it is much simpler to understand. Thanks again for your timely feedback!