Hi everyone,
I am analyzing differential gene expression for a plant based on 3'cDNA tag libraries generated from total RNA. We are looking at a fungal treatment (n=6) compared to a sterile control as well as condition (Condition W vs Condition D).
I am using DESeq2 to determine DEGs of the plant when under a fungal treatment compared to the sterile control within the condition. Then I want to use RRPP to compare DEGs and determine which DEGs are unique for each treatment and condition.
I have never used RRPP and reading up on it. Currently, when I import my count data from tximport, I set the design as follows and create a DESeq object:
ddsTxi <- DESeqDataSetFromTximport(txi, sampleTable, design = ~condition + fungal_tr + condition:fungal_tr )
keep <- rowSums(counts(ddsTxi)) > 1
ddsTxi <- ddsTxi[keep,]
ddsTxi$fungal_tr <- relevel(ddsTxi$fungal_tr, ref = "CON")
ddsTxi$condition <- relevel(ddsTxi$condition, ref = "W")
dds_DESeq <- DESeq(ddsTxi)
So my starting questions are:
- Can I use the DESeq object and format for RRPP or do I need to transform via VST?
- Is my design incorrect for what I want to do with RRPP downstream?
Thank you and let me know what I need to clarify if need be!
Hi Michael, thanks for responding!
Randomization of residuals permutation procedure (RRPP) is an R package developed by Collyer and Adams but is not a Bioconductor package to my knowledge/checking online.
Our plan was to obtain a matrix of DEGs of the plant for each fungal treatment compared to the sterile control within each condition and then use RRPP to compare that data to determine global transcriptome changes between treatment and conditions. I'm worried I shouldn't be setting up the count matrix or design the way I have it for the DESeq2. Instead I should have one count matrix being for condition W and another being for condition D, and then each of these matrices have the design = ~fungal_tr.
The reason behind this concern may be due to mis-understanding the DESeq2 output. With the design I posted in my question my interpretation of the results was the following:
If this is correct, then I think that when I run RRPP I wouldn't be able to look at DEGs of FungalTreatment1W v.s. DEGs FungalTreatment1D, because I don't have the former.
Let me know if I need to clarify anything. Also side note, your blog and Bioconductor RNAseq workflow docs have been super helpful!
Take a look at the interaction section in the vignette. Your interpretation of 1-3 isn’t exactly correct but we have a diagram of the interactions in the section that may help. Otherwise I’d recommend speaking to a statistician to help with interpretation of the model.