I'm quite confused about using DESeq2 to find the differential abundant taxa in microbiome studies, especially when there are more than two groups of the factor. I know DESeq2 was initially used for RNA-seq to detect the regulation of gene expressions. It's easy to understand when there are only two groups, e.g. treated vs. untreated. We can easily say which taxa was up-regulated by looking at the log2fold change (positive or negative).
BUT what about when there are three groups, control vs. treat1 vs. treat2? The DESeq2 can still handle this situation, but then I have no idea how to interpret the log2fold change. If we detected some taxa that were significantly different, how can we know in which group these taxa were up-regulated?
Although some people suggest to do the pairwise comparison, it's still unclear to me how to do it and how to interpret it. Does anyone have very good recommendation or idea about this?
Thanks in advance.
Hi Haitao Wang,
This is definitely confusing, and the meaning of the fold-change returned automatically is a bit convoluted. Ultimately, the reported fold-change depends on the models (full vs. reduced) you provide DESeq2 when performing the likelyhood ratio test (LRT), since the LRT goodness-of-fit tests the hypothesis that the data are _-times more likely to fit the full model (~treat1 + treat2) than the reduced (~1) model. This result has little biological meaning when examining multiple variables (though it is essentially the same as fold-change with Wald when comparing between two).
It is therefore more useful, once you have computed the p-value statistic, to report the fold-changes of treat1-vs-control and treat2-vs-control using the following to examine the possible iterations: resultsNames(dds). It will give you something like: "Intercept" "Treatment_treat1_vs_control" etc... You can then select the appropriate fold-change when you generate the data table:
resdf
Hi Haitao Wang,
This is definitely confusing, and the meaning of the fold-change returned automatically is a bit convoluted. Ultimately, the reported fold-change depends on the models (full vs. reduced) you provide DESeq2 when performing the likelyhood ratio test (LRT), since the LRT goodness-of-fit tests the hypothesis that the data are _-times more likely to fit the full model (~treat1 + treat2) than the reduced (~1) model. This result has little biological meaning when examining multiple variables (though it is essentially the same as fold-change with Wald when comparing between two).
It is therefore more useful, once you have computed the p-value statistic, to report the fold-changes of treat1-vs-control and treat2-vs-control using the following to examine the possible iterations: resultsNames(dds). It will give you something like: "Intercept" "Treatment_treat1_vs_control" etc... You can then select the appropriate fold-change when you generate the data table:
resdf
Dear Mark Pepin , thanks very much for your answer. Now I have a better understanding about the results. resultsNames(dds) really helps! The default result from results(dds) only grabs the last comparison of my pairwise comparisons. Let's say I have group A, B and C. But the resultsNames(dds) only returns "Intercept", "B_vs_A", and "C_vs_A". What does "Intercept" mean here? And how about "B_vs_C"?
Hi Haitao Wang,
The explanation for the "intercept" requires some knowledge of the Likelihood ratio test, which is again a goodness-of-fit with two models, the "full" and a "reduced" (AKA "intercept-only") model. So, the intercept takes the variable(s) used in your reduced model and creates a baseline expression differences (according to that "intercept" variable/model) that is hypothetically improved by the full model.
If you have a specific comparison that you want to report, you can use the "Contrast" function to get "B vs. C", as follows:
results(dds, contrast=c("treatment", "B", "C"))
- Badges
- Science topic
- Similar topics
- Mathematics
- Statistics