DESeq2 workflow problem

[Sentinel156]

Hi everyone,

I'm investigating differentially expressed genes in my wild-type yeast vs a mutant strain. Each strain has been cultured in three different conditions (glucose, sucrose, glycerol). I've got two replicates of Illumina 100bp single end reads per condition and I'm analysing these with DESeq2 (v1.6.3).

Taking the glucose replicates for example. I've performed two analyses. The first using a counts table containing all replicates from all conditions. For the glucose comparison I specified a pairwise contrast like follows:

myresults <- results(ddsanalysis, contrast=c("condition","mutant_glucose","WT_glucose"))

Comparing this output to the second analysis in which the counts table only contained my glucose replicates, I see ~100 fewer DE genes returned from the first analysis containing all conditions and replicates. The subsets of genes returned are also different. For example, the first analysis suggests a gene (which I know to be upregulated in the mutant on glucose via qRT-PCR) is not among the significant DE genes; yet this gene is the most significant DE gene in my second analysis.

While I don't fully understand the modelling that DESeq2 conducts during the analysis, It seems to me that information is drawn from all replicates supplied via the counts table regardless of whether they are specified as part of a pairwise contrast or not and that this is influencing the determination of DE genes in the first analysis.

What would the community recommend as the best course of action for this analysis? To create separate counts tables for the pairwise comparisons I wish to perform? Or to continue using the contrasts method as I have described?

Thanks in advance,

[Michael Love]

"information is drawn from all replicates supplied via the counts table regardless of whether they are specified as part of a pairwise contrast"

You are correct, all samples in the DESeqDataSet are used to estimate the dispersion. Adding or removing samples will change the dispersion estimates, which will change p-values. So it is not surprising to see more or less genes passing an FDR threshold after adding or removing samples. p-values are tail probabilities which are highly sensitive to model parameters. Usually it's best to include all the samples for dispersion estimation. An exception to this is if you have plenty of biological replicates, and the within-group variance of the groups are very different. Look at a PCA plot (see vignette) and see if the different groups have similar or very different within-group spread. If they are very different and you have plenty of biological replicates, then you might want to run DESeq() on each group separately.

[Sentinel156]

Thanks for the reply Michael,

After looking over a PCA plot containing all of my replicates the glycerol condition replicates do not clearly cluster independently of the other conditions. It is possible that including these replicates in my initial analysis inflated the dispersion estimates and in turn, affected the DE gene predictions. Removing the glycerol condition results in a set of predicted DE genes more similar to that when only analysing the glucose condition. I will likely remove the glycerol replicates from my analysis and analyse them independently however there may be too much variability to resolve with only two biological replicates.