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:
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:
and using the following design formulae to analyse with DESeq2:
Using this design matrix and formulae, I got the following output from the resultsNames function:
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:
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:
Thanks in advance for any help.
For information, I use R version 3.3.0 (2016-05-03) and DESeq2_1.12.4.
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:
Code:
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
Code:
> 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
Code:
> df ~ Disease + Disease:PN + Disease:Treatment
Code:
> 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"
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:
Code:
> 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"
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:
Code:
(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 ...
For information, I use R version 3.3.0 (2016-05-03) and DESeq2_1.12.4.
Comment