Interaction Factors Removing All Variability
Entering edit mode
2.2 years ago
tyak9569 • 0

Hi all!

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)

plotDispEsts(dds )

normed = counts(dds, normalized=TRUE)
write.csv(normed, file="all_samples_normCounts.csv")

res <- results(dds)
write.csv(, file="all_samples_results.csv")
DESeq2::plotMA(res, alpha = 0.1, ylim=c(-3,3), cex=.4)
abline(h=c(-1,1), col="dodgerblue", lwd=2)

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.


DESeq DESeq2 RNA-Seq R rna-seq • 650 views
Entering edit mode
2.0 years ago
antonioggsousa ★ 2.5k


I believe that when you're using a model formula like yours, with multiple factors and/or interactions, DESeq2 will provide you multiple pairwise comparisons, but the last contrast will be used to plot dispersion, retrieve results, etc, unless that you specify the contrast or name of the specific pairwise comparison that you want to retrieve

So, you can check which contrast is being used by ploMA() by doing resultsNames(dds). I would say that in both cases you're plotting different pairwise comparisons, that's why you're getting such a different result.

To call a specific pairwise comparison to plot the MA, you can do after running the resultsNames(dds):

plotMA(results(dds, name=<name-specific-pairwise-comparison-given-by-resultsNames(dds)> ))

To assess variability the best way is doing a PCA or other multidimensional reduction technique and do not rely on this MA plot, because this plot is a consequence of what you can see on PCA before, that is, if you don't have variability among the groups that you want to explore, you'll not find genes diferentially expressed.

I hope this helps to answer your question!



Login before adding your answer.

Traffic: 1880 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6