Extract contrast from DESeq2 advanced design

[Michel Meyer]

Dear all,

I wish to analyze RNAseq experiment data with complex matrix design. The matrix I wish to analyze is the same as displayed in the edgeR documentation section 3.5:

   Disease  Patient Treatment
1  Healthy  1       None
2  Healthy  1       Hormone
3  Healthy  2       None
4  Healthy  2       Hormone
5  Healthy  3       None
6  Healthy  3       Hormone
7  Disease1 4       None
8  Disease1 4       Hormone
9  Disease1 5       None
10 Disease1 5       Hormone
11 Disease1 6       None
12 Disease1 6       Hormone
13 Disease2 7       None
14 Disease2 7       Hormone
15 Disease2 8       None
16 Disease2 8       Hormone
17 Disease2 9       None
18 Disease2 9       Hormone

So, to analyse this kind of matrix design with DESeq2, I have followed the instructions in the DESeq2 vignette section 3.12.1. Creating a "nested" Patient variable (called PN) like that:

> dm
   Disease  PN Treatment
1  Healthy  P1 None
2  Healthy  P1 Hormone
3  Healthy  P2 None
4  Healthy  P2 Hormone
5  Healthy  P3 None
6  Healthy  P3 Hormone
7  Disease1 P1 None
8  Disease1 P1 Hormone
9  Disease1 P2 None
10 Disease1 P2 Hormone
11 Disease1 P3 None
12 Disease1 P3 Hormone
13 Disease2 P1 None
14 Disease2 P1 Hormone
15 Disease2 P2 None
16 Disease2 P2 Hormone
17 Disease2 P3 None
18 Disease2 P3 Hormone

and using the following design formulae to analyse with DESeq2:

> df
~ Disease + Disease:PN + Disease:Treatment

Using this design matrix and formulae, I got the following output from the resultsNames function:

> dds <- DESeqDataSetFromMatrix(countData=rawData, colData=dm, design=df) 
> dds <- DESeq(dds, quiet=TRUE)
> resultsNames(dds)
[1] "Intercept"             "Disease_Disease2_vs_Disease1"  "Disease_Healthy_vs_Disease1"   "DiseaseDisease1.PNP2"         
[5] "DiseaseDisease2.PNP2"  "DiseaseHealthy.PNP2"           "DiseaseDisease1.PNP3"          "DiseaseDisease2.PNP3"         
[9] "DiseaseHealthy.PNP3"   "DiseaseDisease1.TraitmentNone" "DiseaseDisease2.TraitmentNone" "DiseaseHealthy.TraitmentNone"

So, now I can extract comparisons DiseaseDisease1.TraitmentNone vs DiseaseHealthy.TraitmentNone and DiseaseDisease2.TraitmentNone vs DiseaseHealthy.TraitmentNone, using contrast option of results function.

However, I'm also interested in comparisons DiseaseDisease1.TraitmentHormone vs DiseaseHealthy.TraitmentHormone and DiseaseDisease2.TraitmentHormone vs DiseaseHealthy.TraitmentHormone. To get these comparisons, I have rerun a DESeq2 analyses after inverting the factor levels from the Treatment condition. So, now I can access to the other ones:

> dm$Traitment
[1] None Hormone None Hormone None Hormone None Hormone None Hormone None Hormone None Hormone None Hormone None Hormone
Levels: Hormone None
> dm$Traitment <- relevel(dm$Traitment,'None')
> dm$Traitment
[1] None Hormone None Hormone None Hormone None Hormone None Hormone None Hormone None Hormone None Hormone None Hormone
Levels: None Hormone
> dds <- DESeqDataSetFromMatrix(countData=rawdata, colData=dm, design=df) 
> dds <- DESeq(dds, quiet=TRUE)
> resultsNames(dds)
[1] "Intercept"             "Disease_Disease2_vs_Disease1"     "Disease_Healthy_vs_Disease1"      "DiseaseDisease1.PNP2"            
[5] "DiseaseDisease2.PNP2"  "DiseaseHealthy.PNP2"              "DiseaseDisease1.PNP3"             "DiseaseDisease2.PNP3"            
[9] "DiseaseHealthy.PNP3"   "DiseaseDisease1.TraitmentHormone" "DiseaseDisease2.TraitmentHormone" "DiseaseHealthy.TraitmentHormone"

I wish to know if it is the right way to get the requested comparisons?
Because I have seen that the results between (1) DiseaseDisease1.TraitmentNone vs DiseaseHealthy.TraitmentNone and (2) DiseaseDisease2.TraitmentHormone vs DiseaseHealthy.TraitmentHormone are nearly the same:

 (1)        (2)
 0,268438  -0,268438
 0,141639  -0,141639
 0,000000   0,000000
-0,696627   0,696634
 0,000000  -0,973195
 0,034611  -0,034611
 0,530830  -0,530829
-0,118112   0,118113
-0,016033   0,016033
 0,270009  -0,270009
 3,021271  -3,021392
 0,000000   0,000000
 0,000000   0,000000
 1,617347  -1,617340
-0,491919   0,491938
-0,380177   0,380001
 0,528211  -0,528204
-0,052266   0,052266
 0,204210  -0,204210
 0,267654  -0,267591
 0,012433  -0,012432
-0,046289   0,046290
 ...

Thanks in advance for any help.
For information, I use R version 3.3.0 (2016-05-03) and DESeq2_1.12.4.

[shunyip]

Hello Michel,

While your way of doing it can be correct, maybe you can use a much simpler way to perform this analysis.

For example, you can extract DiseaseDisease1.TreatmentNone and DiseaseHealthy.TreatmentNone's samples into one matrix. Then, create a much less complicated design matrix to compare them. (1 condition vs 1 condition) This way, you will definitely know that you are doing the right thing.

Best,

[Michel Meyer]

Thank you for your answer.