Hi:
I am new to RNA-seq analysis, R and DESeq. I'm trying to analyze 12 samples divided into four groups, as follows:
| Samples | Group | Treatment |
|---|---|---|
| S1 | Control | Untreated |
| S2 | Control | Untreated |
| S3 | Control | Untreated |
| S4 | PKD | Untreated |
| S5 | PKD | Untreated |
| S6 | PKD | Untreated |
| S7 | PKD | Treated |
| S8 | PKD | Treated |
| S9 | PKD | Treated |
| S10 | Control | Treated |
| S11 | Control | Treated |
| S12 | Control | Treated |
Q: I'm interesting in genes which are unique DEGs for PKD. To achieve this, I followed the following code:
library(DESeq2)
counts <- read.table("count.txt", header=TRUE, row.names=1)
samples <- read.table("above edf file.txt", header=TRUE)
ds <- DESeqDataSetFromMatrix(countData=counts, colData=samples, design=~group+treatment+group:treatment)
ds$group <- relevel(ds$group, "PKD")
ds$treatment <- relevel(ds$treatment, "untreated")
ds=DESeq(ds)
resultsNames(ds)
> resultsNames(ds)
[1] "Intercept" "group_PKD_vs_Control" "treatment_Treated_vs_Untreated" "groupPKD.treatmentTreated"
Which of the resultNames contain my genes of interest (unique to PKDs)? How can I achieve that? please help.
Thank you,
Sofia
Michael:
Thank you for your reply. I'm interesting in those genes which are differential expressed in PKD i.e.
1. PKD-untreated vs Control-untreated
2. PKD-treated vs Control-treated
But I like to combine these aforementioned comparisons (providing entire model but not separated) to get a single result file, containing only those genes that are differentially expressed during PKD (untreated + Treated). So my end point is to find out the genes which are significant in PKD (treated+untreated) vs Control (treated+Untreated).
Thank you for your time and help.
Sofia
If you need them to be DE in both groups, you can't get a single p-value for that. You can make the two comparisons, and then look at the set of genes with adjusted p-value < x in both groups. To make the two comparisons, I'd recommend following the advice in the vignette (the very beginning of the section I linked above) about combining the two variables into one, called "group", and then using the 'contrast' argument to make comparisons. It will look something like
res.trt <- results(dds, contrast=c("group","PKDTreated","ControlTreated"))Michael:
Thank you for the suggestion of comparing with a cut-off of pvalue. Am I missing something (other comparison) from the above experiment which I should also compare? Any suggestion on this will be greatly appreciated.
Thank you very much!
Sofia
You mentioned two comparisons, and I showed example code for how to produce one of those. You can follow the pattern for how to get the other.
Thank you Michael.
Sorry to trouble you again but after reading interactions, I still have some doubts. I would appreciate if you can provide a brief explanation on the following results:
What do "group_PKD_vs_Control", "treatment_Treated_vs_Untreated" and "groupPKD.treatmentTreated" stand for in terms of comparisons?
Thanks
Sofia
If the vignette doesn't help clarify, I'd recommend speaking with a statistician who can better explain interactions in linear models in person. The vignette is my best attempt at explain the concept without being able to speak to a person face-to-face. My rephrasing the words from the vignette here probably won't help much, but speaking to someone in person probably will help a lot.
My suggestion for your code though is to not use a design with interactions for this model, but to combine the two variables into one. It accomplishes the same thing, while avoiding the confusion users have with the meaning of terms produced by the design with interactions.