I've been running into some wackiness when it comes to the creation of my DESeq object during the analysis of my RNA-Seq data. To give you some background, we have 8 total grouping, with two genotypes (transgenic and control), two treatment conditions, and two time points. For my masters' project, I put in our three independent variables into my DESeq object like so:
dds = DESeqDataSetFromMatrix(countData = counts, colData = coldata, design= ~treatment + time + genotype)
I ended up with differentially expressed genes with an effect of each variable besides treatment. In a previous experiment with a similar paradigm, we noted a significant genotype:time interaction, and I wanted to include that in the model just to see if anything new popped up. What I found was that including the interaction into my DESeq object completely removed any variability in my RNASeq data set, completely flattening my MA plots and not identifying any genes of interest. My code is as follows:
dds = DESeqDataSetFromMatrix(countData = counts1, colData = coldata, design = ~ time:genotype + treatment) keep <- rowSums(counts(dds)) >= 20 dds <- dds[keep,] dds$treatment <- relevel(dds$treatment, ref = "same") dds$genotype <- relevel(dds$genotype, ref = "ctrl") dds = DESeq(dds) pdf("all_samples_dispersionPlot.pdf") plotDispEsts(dds ) dev.off() normed = counts(dds, normalized=TRUE) write.csv(normed, file="all_samples_normCounts.csv") res <- results(dds) write.csv(as.data.frame(res), file="all_samples_results.csv") pdf("all_samples_MAPlot.pdf") DESeq2::plotMA(res, alpha = 0.1, ylim=c(-3,3), cex=.4) abline(h=c(-1,1), col="dodgerblue", lwd=2) dev.off()
Here's what my MA plot looked like previously:
Here's what it looks like with the interaction thrown in:
If you guys have any suggestions or explanation that I'm not aware of I'd love to hear it. My initial conclusion is that all of the variability in gene expression is due to the interaction, and that my data is basically useless, and I'm just hoping I'm wrong.