Module 10: DE Testing

Objectives

• Generate tables of DE results
• Understand what a p-value represents.
• Understand multiple hypothesis correction application and importance

Differential Expression Workflow

Now we test for differential expression between our groups of interest and return a table of results.

---

Testing for DE

Before showing the code for generating differential expression results, let’s walk through some toy data and discuss our intuition, its limitations, and what we can do to formalize our thought around determining if a gene is differentially expressed.

Let’s get a sense for our intuition by looking at expression boxplots for a few “genes” where we’ve made up the data. In each plot, we are comparing the expression levels (on the y-axis) for samples (each point) representing a WT (red) and KO (teal) condition. The boxplot shows the 25% - 75% distribution, along with outliers, with a bar representing the median value, and a black-crossed-point representing the mean.

For the first example, we ask: Does the KO affect the expression of the gene?

It seems pretty clear that there is a large difference between the means of the two groups, and within-group spread is quite low. So the answer to the question is likely “Yes”.

Consider a second example:

Here there isn’t a large difference between the means of the two groups, and there is quite a bit of spread within-group. The answer here is likely “No”.

Finally, consider a third example:

The means are not so close to each other in this example, though there is still quite a bit of spread. This is example is perhaps less clear.

Looking at all three of the hypothetical genes together, we see them in relation to one another and how they span the range from No, to Maybe?, to Definitely.

Consider needing to make this decision for 20,000 genes. Even if all the genes behaved in a clear-cut manner, that would take a lot of time (or a lot of graduate students). But there are probably lots of Gene 2s out there.

We need a formal, reproducible, way to make this decision!

There are many tools to model expression between the groups, we have chosen to focus on DESeq2, but know that a tool like edgeR is out there. Both tools are statistically sound approachees to modeling and testing for differential expression.

With any statistical test, we need to clearly state what we are testing. For the toy data above, we assume that the KO doesn’t affect the expression of any particular gene. Statisticians would call this the “null hypothesis”. For the mouse data we’re working with, the null hypothesis is that “there is no difference between the expression of the iron deficient (minus) and normal (plus) mice”. The alternative is that there is a difference.

For each gene, DESeq2 computes a “Wald statistic” which is a single number encapsulating the difference in the means and the spread of the groups. However, this number alone doesn’t tell us how to decide if a gene is differentially expressed. We need a second number to give us an idea of how extreme that Wald statistic is among the distribution of Wald statistics, this is the “p-value”.

Imagine shuffling the KO and WT labels and recomputing the statistic over and over again. You’d get a distribution of statistics that would look similar to a normal curve. The p-value essentially tells you how likely you are to have seen the statistic you see by chance. So when we set a p-value = 0.05 as a threshold, we’re saying, “there is a 5% chance I’d see something this extreme when there was actually no effect”. So the evidence is strong, but not ironclad.

Generating DE Results

We can check what comparisons were automatically generated during fitting using the resultsNames() function.

resultsNames(dds)

[1] "Intercept"               "condition_minus_vs_plus"

There is only the one comparison in the results, so we will refer to it in the name parameter of the results() function, and assign the result as an object.

results_minus_vs_plus = results(dds, name = 'condition_minus_vs_plus')
head(results_minus_vs_plus)

log2 fold change (MLE): condition minus vs plus 
Wald test p-value: condition minus vs plus 
DataFrame with 6 rows and 6 columns
                      baseMean log2FoldChange     lfcSE      stat    pvalue
                     <numeric>      <numeric> <numeric> <numeric> <numeric>
ENSMUSG00000000001  1489.83039       0.297760  0.210310  1.415815  0.156830
ENSMUSG00000000028  1748.93544       0.226421  0.176795  1.280695  0.200301
ENSMUSG00000000031  2151.87715       0.457635  0.764579  0.598545  0.549476
ENSMUSG00000000037    24.91672       0.579130  0.561259  1.031840  0.302147
ENSMUSG00000000049     7.78377      -0.899483  1.553063 -0.579167  0.562476
ENSMUSG00000000056 19653.54030      -0.174048  0.203529 -0.855151  0.392468
                        padj
                   <numeric>
ENSMUSG00000000001  0.868573
ENSMUSG00000000028  0.902900
ENSMUSG00000000031  0.995391
ENSMUSG00000000037  0.950613
ENSMUSG00000000049  0.998043
ENSMUSG00000000056  0.982479

In the results table, the row names are gene identifiers (in theis case ENSEMBL IDs because that’s what the GTF we used in the call to RSEM+STAR used), and the columns have the following definitions:

1. baseMean is the average of the normalized count values, divided by size factors and taken over all samples, and can be interpreted as the relative expression level of that gene across all samples.
2. log2FoldChange is the log2 transformed ratio of the expression of the numerator group (first group) over the denominator group (second group after “vs”). Note that in our comparison, the log2FoldChange column compares the expression of the numerator group (minus) over the denominator group (plus). If the value is positive, that means the expression of that gene is greater across the minus samples than across the plus samples. If the value is negative, that means the expression of that gene is greater across the minus samples.
3. lfcSE is the standard error for the log2 fold change estimate.
4. stat is the calculated Wald statistic for that gene.
5. pvalue is the nominal significance for that gene.
6. padj is the adjusted p-value and is what we use for determining significantly differently expressed genes.

Note: results() defaults

If no arguments are passed to results(), then the log2 fold changes and Wald test p-value will be for the last variable in the design formula, and if this is a factor, the comparison will be the last level over the reference level. If you specify name, as we did above, then the behavior is given by the name used from resultsNames().

There are multiple ways to specify the test to be done using the results() function. It is especially helpful to know this when fitting more complex models and testing more complex contrasts. To demonstrate this, consider this description from the help for results():

contrast: a character vector with exactly three elements: the name of a factor in the design formula, the name of the numerator level for the fold change, and the name of the denominator level for the fold change

So an alternative way to test the same contrast as above (i.e. plus / minus) is:

alt_results_minus_vs_plus = results(dds, contrast = c('condition', 'minus', 'plus'))
head(alt_results_minus_vs_plus)

log2 fold change (MLE): condition minus vs plus 
Wald test p-value: condition minus vs plus 
DataFrame with 6 rows and 6 columns
                      baseMean log2FoldChange     lfcSE      stat    pvalue
                     <numeric>      <numeric> <numeric> <numeric> <numeric>
ENSMUSG00000000001  1489.83039       0.297760  0.210310  1.415815  0.156830
ENSMUSG00000000028  1748.93544       0.226421  0.176795  1.280695  0.200301
ENSMUSG00000000031  2151.87715       0.457635  0.764579  0.598545  0.549476
ENSMUSG00000000037    24.91672       0.579130  0.561259  1.031840  0.302147
ENSMUSG00000000049     7.78377      -0.899483  1.553063 -0.579167  0.562476
ENSMUSG00000000056 19653.54030      -0.174048  0.203529 -0.855151  0.392468
                        padj
                   <numeric>
ENSMUSG00000000001  0.868573
ENSMUSG00000000028  0.902900
ENSMUSG00000000031  0.995391
ENSMUSG00000000037  0.950613
ENSMUSG00000000049  0.998043
ENSMUSG00000000056  0.982479

This way of calling results() is especially helpful when the levels of the column of interest contain more than two levels because you can specify exactly which levels to test with little confusion.

Question

Why should we use values from padj instead of the pvalue? Post in the Slack thread.

Multiple hypothesis testing and FDR correction

Each p-value is the result of a single test for a single gene. With a p-value < 0.05 significance cut-off, there is a 5% chance it is a false positive. The more genes we test, the greater chance we have of seeing a significant results by chance. So if we are testing 20,000 genes for differential expression, we would expect to see ~1,000 significant genes just by chance.

To address this we will correct for multiple hypothesis testing to reduce the number of false positives. While there are a few approaches, the default method is the False Discovery Rate (FDR) (Benjamini and Hochberg (1995)).

The default FDR rate cutoff for our analyses is 0.05, meaning the proportion of false positives amongst our differentially expressed genes is controlled to 5%. So if we call 500 genes as differentially expressed with this FDR cutoff, we expect only 25 of them to be false positives. DESeq2 vignette’s includes a further discussion of filtering and multiple testing.

Note on padj values set to NA

As discussed in the HBC tutorial as well as the DESeq2 vignette, DESeq2 reduces the number of genes that will be tested by removing genes with low number of counts and outlier samples.

• If within a row, all samples have zero counts, the baseMean column will be zero, and the log2 fold change estimates, p-value and adjusted p-value will all be set to NA.
• If a row contains a sample with an extreme count outlier then the p-value and adjusted p-value will be set to NA. These outlier counts are detected by Cook’s distance.
• If a row is filtered by automatic independent filtering, e.g. for having a low mean normalized count, then only the adjusted p-value will be set to NA.

Now that we’ve generated our differential comparisons and have an understanding of our results, including multiple hypothesis correction, we can proceed with generating summary figures and tables for our differential expression analysis.

---

Summary

In this section, we:

• Performed statistical tests for comparisons of interest
• Generated tables of differential expression results - i.e. fold changes and adjusted pvalues for each gene in dataset
• Discussed importance and application of multiple hypothesis correction

Now that we’ve generated our differential comparisons and have an understanding of our results, including multiple hypothesis correction, we can proceed with generating summary figures and tables for our differential expression analysis.

---

Sources

• HBC DGE training module, part 1: https://hbctraining.github.io/DGE_workshop/lessons/04_DGE_DESeq2_analysis.html
• HBC DGE training module, part 2: https://hbctraining.github.io/DGE_workshop/lessons/05_DGE_DESeq2_analysis2.html
• DESeq2 vignette: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#differential-expression-analysis

---

---

These materials have been adapted and extended from materials listed above. These are open access materials distributed under the terms of the Creative Commons Attribution license (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.