DESeq2: Model fitting and Hypothesis testing¶

Last updated: November 28, 2023

The final step in the DESeq2 workflow is taking the counts for each gene and fitting it to the model and testing for differential expression.

Generalized Linear Model¶

As described earlier, the count data generated by RNA-seq exhibits overdispersion (variance > mean) and the statistical distribution used to model the counts needs to account for this. As such, DESeq2 uses a negative binomial distribution to model the RNA-seq counts using the equation below:

The two parameters required are the size factor, and the dispersion estimate. Next, a generalized linear model (GLM) of the NB family is used to fit the data. Modeling is a mathematically formalized way to approximate how the data behaves given a set of parameters.

After the model is fit, coefficients are estimated for each sample group along with their standard error. The coefficients are the estimates for the log2 fold changes, and will be used as input for hypothesis testing.

Hypothesis testing¶

The first step in hypothesis testing is to set up a null hypothesis for each gene. In our case, the null hypothesis is that there is no differential expression across the two sample groups (LFC == 0). Notice that we can do this without observing any data, because it is based on a thought experiment. Second, we use a statistical test to determine if based on the observed data, the null hypothesis can be rejected.

Wald test¶

In DESeq2, the Wald test is the default used for hypothesis testing when comparing two groups. The Wald test is a test usually performed on parameters that have been estimated by maximum likelihood. In our case we are testing each gene model coefficient (LFC) which was derived using parameters like dispersion, which were estimated using maximum likelihood.

DESeq2 implements the Wald test by:

1. Taking the LFC and dividing it by its standard error, resulting in a z-statistic.
2. The z-statistic is compared to a standard normal distribution, and a p-value is computed reporting the probability that a z-statistic at least as extreme as the observed value would be selected at random.
3. If the p-value is small we reject the null hypothesis and state that there is evidence against the null (i.e. the gene is differentially expressed).

The model fit and Wald test were already run previously as part of the DESeq() function:

## DO NOT RUN THIS CODE

## Create DESeq2Dataset object
dds <- DESeqDataSetFromTximport(txi,
                                   colData = meta %>% column_to_rownames("sample"), 
                              design = ~ condition)
## Run analysis
dds <- DESeq(dds)

Likelihood ratio test (LRT)¶

An alternative to pair-wise comparisons is to analyze all levels of a factor at once. By default, the Wald test is used to generate the results table, but DESeq2 also offers the LRT which is used to identify any genes that show change in expression across the different levels. This type of test can be especially useful in analyzing time course experiments.

To use the LRT, we use the DESeq() function but this time adding two arguments:

1. specifying that we want to use the LRT test
2. the "reduced" model

# The full model was specified previously with the `design = ~ condition`:
# dds <- DESeqDataSetFromTximport(txi,
                                   # colData = meta %>% column_to_rownames("sample"), 
                              # design = ~ condition)

# Likelihood ratio test
dds_lrt <- DESeq(dds, test="LRT", reduced = ~ 1)

Since our "full" model only has one factor (sampletype), the "reduced" model (removing that factor) is just the intercept (~ 1). You will find that similar columns are reported for the LRT test. One thing to note is, even though there are fold changes present they are not directly associated with the actual hypothesis test.

full_model <- ~ genotype + treatment + time + treatment:time

reduced_model <- ~ genotype + treatment + time

dds_lrt <- DESeqDataSetFromMatrix(countData = data, colData = metadata, design = ~ genotype + treatment + time + treatment:time)

dds_lrt_time <- DESeq(dds_lrt, test="LRT", reduced = ~ genotype + treatment + time)

Multiple test correction¶

Regardless of whether we use the Wald test or the LRT, each gene that has been tested will be associated with a p-value. It is this result which we use to determine which genes are considered significantly differentially expressed. However, we cannot use the p-value directly.

What does the p-value mean?¶

A gene with a significance cut-off of p \< 0.05, means there is a 5% chance it is a false positive. For example, if we test 20,000 genes for differential expression, at p \< 0.05 we would expect to find 1,000 DE genes by chance. If we found 3000 genes to be differentially expressed total, roughly one third of our genes are false positives! We would not want to sift through our "significant" genes to identify which ones are true positives.

Since each p-value is the result of a single test (single gene). The more genes we test, the more we inflate the false positive rate. This is the multiple testing problem.

Correcting the p-value for multiple testing¶

There are a few common approaches for multiple test correction:

• Bonferroni: The adjusted p-value is calculated by: \(p-value*m\) (\(m\) = total number of tests). This is a very conservative approach with a high probability of false negatives, so is generally not recommended.
• FDR/Benjamini-Hochberg: Benjamini and Hochberg (1995) defined the concept of False Discovery Rate (FDR) and created an algorithm to control the expected FDR below a specified level given a list of independent p-values. More info about BH.
• Q-value / Storey method: The minimum FDR that can be attained when calling that feature significant. For example, if gene X has a q-value of 0.013 it means that 1.3% of genes that show p-values at least as small as gene X are false positives.

DESeq2 helps reduce the number of genes tested by removing those genes unlikely to be significantly DE prior to testing, such as those with low number of counts and outlier samples (see below). However, multiple test correction is also implemented to reduce the False Discovery Rate using an interpretation of the Benjamini-Hochberg procedure.

Exploring Results (Wald test)¶

By default DESeq2 uses the Wald test to identify genes that are differentially expressed between two sample groups. Given the factor(s) used in the design formula, and how many factor levels are present, we can extract results for a number of different comparisons. Here, we will walk you through how to obtain results from the dds object and provide some explanations on how to interpret them.

Specifying contrasts¶

In our dataset, we have three sample groups so we can make three possible pairwise comparisons:

1. Control vs. Vampirium
2. Garlicum vs. Vampirium
3. Garlicum vs. Control

We are really only interested in 1 and 2 from above. When we initially created our dds object we had provided ~ condition as our design formula, indicating that sampletype is our main factor of interest.

To indicate which two sample groups we are interested in comparing, we need to specify contrasts. The contrasts are used as input to the DESeq2 results() function to extract the desired results.

Contrasts can be specified in three different ways:

1.Contrasts can be supplied as a character vector with exactly three elements: the name of the factor (of interest) in the design formula, the name of the two factors levels to compare. The factor level given last is the base level for the comparison. The syntax is given below:

# DO NOT RUN!
contrast <- c("condition", "level_to_compare", "base_level")
results(dds, contrast = contrast)

2.Contrasts can be given as a list of 2 character vectors: the names of the fold changes for the level of interest, and the names of the fold changes for the base level. These names should match identically to the elements of resultsNames(object). This method can be useful for combining interaction terms and main effects.

# DO NOT RUN!
resultsNames(dds) # to see what names to use
contrast <- list(resultsNames(dds)[1], resultsNames(dds)[2])
results(dds, contrast = contrast)

3.One of the results from resultsNames(dds) and the name argument. This one is the simplest but it can also be very restrictive:

# DO NOT RUN!
resultsNames(dds) # to see what names to use
results(dds, name = resultsNames(dds)[2])

Alternatively, if you only had two factor levels you could do nothing and not worry about specifying contrasts (i.e. results(dds)). In this case, DESeq2 will choose what your base factor level based on alphabetical order of the levels.

To start, we want to evaluate expression changes between the control samples and the vampirium samples. As such we will use the first method for specifying contrasts and create a character vector:

## your code here
#contrast_cont <- 

meta

contrast_cont <- c("condition", "control", "vampirium")

The results table¶

Now that we have our contrast created, we can use it as input to the results() function. Let's take a quick look at the help manual for the function:

?results

You will see we have the option to provide a wide array of arguments and tweak things from the defaults as needed. As we go through the lesson we will keep coming back to the help documentation to discuss some arguments that are good to know about.

## Extract results for Control vs Vampirium samples
res_tableCont <- results(dds, contrast=contrast_cont, alpha = 0.05)

The results table that is returned to us is a DESeqResults object, which is a simple subclass of DataFrame. In many ways it can be treated like a dataframe (i.e when accessing/subsetting data), however it is important to recognize that there are differences for downstream steps like visualization.

# Check what type of object is returned
class(res_tableCont)

Now let's take a look at what information is stored in the results:

# What is stored in results?
res_tableCont %>% 
  data.frame() %>% 
  head()

We have six columns of information reported for each gene (row). We can use the mcols() function to extract information on what the values stored in each column represent:

# Get information on each column in results
mcols(res_tableCont, use.names=T)

• baseMean: mean of normalized counts for all samples
• log2FoldChange: log2 fold change
• lfcSE: standard error
• stat: Wald statistic
• pvalue: Wald test p-value
• padj: BH adjusted p-values

p-values¶

The p-value is a probability value used to determine whether there is evidence to reject the null hypothesis. A smaller p-value means that there is stronger evidence in favor of the alternative hypothesis.

The padj column in the results table represents the p-value adjusted for multiple testing, and is the most important column of the results. Typically, a threshold such as padj \< 0.05 is a good starting point for identifying significant genes. The default method for multiple test correction in DESeq2 is the FDR. Other methods can be used with the pAdjustMethod argument in the results() function.

Gene-level filtering¶

Let's take a closer look at our results table. As we scroll through it, you will notice that for selected genes there are NA values in the pvalue and padj columns. What does this mean?

The missing values represent genes that have undergone filtering as part of the DESeq() function. Prior to differential expression analysis it is beneficial to omit genes that have little or no chance of being detected as differentially expressed. This will increase the power to detect differentially expressed genes. DESeq2 does not physically remove any genes from the original counts matrix, and so all genes will be present in your results table. The genes omitted by DESeq2 meet one of the three filtering criteria outlined below:

1. Genes with zero counts in all samples

If within a row, all samples have zero counts there is no expression information and therefore these genes are not tested. In our case, there are no genes that fulfill this criteria, since we have already filtered out these genes ourselves when we created our dds object.

# Show genes with zero expression
res_tableCont %>%
  as_tibble(rownames = "gene") %>% 
  dplyr::filter(baseMean==0) %>%
  head()

2. Genes with an extreme count outlier

The DESeq() function calculates, for every gene and for every sample, a diagnostic test for outliers called Cook's distance. Cook's distance is a measure of how much a single sample is influencing the fitted coefficients for a gene, and a large value of Cook's distance is intended to indicate an outlier count. Genes which contain a Cook's distance above a threshold are flagged, however at least 3 replicates are required for flagging, as it is difficult to judge which sample might be an outlier with only 2 replicates. We can turn off this filtering by using the cooksCutoff argument in the results() function.

# Show genes that have an extreme outlier
res_tableCont %>% 
  as_tibble(rownames = "gene") %>% 
  dplyr::filter(is.na(pvalue) & is.na(padj) & baseMean > 0) %>%
  head()

It seems that we have some genes with outliers!

3. Genes with a low mean normalized counts

DESeq2 defines a low mean threshold, that is empirically determined from your data, in which the fraction of significant genes can be increased by reducing the number of genes that are considered for multiple testing. This is based on the notion that genes with very low counts are not likely to see significant differences typically due to high dispersion.

Image courtesy of slideshare presentation from Joachim Jacob, 2014.

At a user-specified value (alpha = 0.05), DESeq2 evaluates the change in the number of significant genes as it filters out increasingly bigger portions of genes based on their mean counts, as shown in the figure above. The point at which the number of significant genes reaches a peak is the low mean threshold that is used to filter genes that undergo multiple testing. There is also an argument to turn off the filtering off by setting independentFiltering = F.

# Filter genes below the low mean threshold
res_tableCont %>% 
  as_tibble(rownames = "gene") %>% 
  dplyr::filter(!is.na(pvalue) & is.na(padj) & baseMean > 0) %>%
  head()

Fold changes¶

Another important column in the results table, is the log2FoldChange. With large significant gene lists it can be hard to extract meaningful biological relevance. To help increase stringency, one can also add a fold change threshold. Keep in mind when setting that value that we are working with log2 fold changes, so a cutoff of log2FoldChange \< 1 would translate to an actual fold change of 2.

res_tableCont_LFC1 <- results(dds, contrast=contrast_cont, alpha = 0.05, lfcThreshold = 1)

The fold changes reported in the results table are the coefficients calculated in the GLM mentioned in the previous section.

Summarizing results¶

To summarize the results table, a handy function in DESeq2 is summary(). Confusingly it has the same name as the function used to inspect data frames. This function, when called with a DESeq results table as input, will summarize the results using a default threshold of padj \< 0.1. However, since we had set the alpha argument to 0.05 when creating our results table threshold: FDR \< 0.05 (padj/FDR is used even though the output says p-value < 0.05). Let's start with the OE vs control results:

## Summarize results
summary(res_tableCont, alpha = 0.05)

In addition to the number of genes up- and down-regulated at the default threshold, the function also reports the number of genes that were tested (genes with non-zero total read count), and the number of genes not included in multiple test correction due to a low mean count.

Extracting significant differentially expressed genes¶

Let's first create variables that contain our threshold criteria. We will only be using the adjusted p-values in our criteria:

### Set thresholds
padj.cutoff <- 0.05

We can easily subset the results table to only include those that are significant using the dplyr::filter() function, but first we will convert the results table into a tibble:

# Create a tibble of results and add gene IDs to new object
res_tableCont_tb <- res_tableCont %>%
  as_tibble(rownames = "gene") %>%
  relocate(gene, .before = baseMean)

head(res_tableCont_tb)

Now we can subset that table to only keep the significant genes using our pre-defined thresholds:

# Subset the tibble to keep only significant genes
sigCont <- res_tableCont_tb %>%
  dplyr::filter(padj < padj.cutoff)

# Take a quick look at this tibble
sigCont

Now that we have extracted the significant results, we are ready for visualization!

contrast_gar <- c("condition", "garlicum", "vampirium")

res_tableGar <- results(dds, contrast = contrast_gar)

padj.cutoff <- 0.05
sigGar <- res_tableGar %>% as_tibble(rownames = "gene") %>%
  relocate(gene, .before = baseMean) %>% dplyr::filter(padj < padj.cutoff)

sigGar

nrow(sigGar) #number of significant genes in Garlicum vs Vampirium
nrow(sigCont) #number of significant genes in Control vs Vampirium

nrow(sigCont) - nrow(sigGar)

This lesson was originally developed by members of the teaching team (Mary Piper, Meeta Mistry, Radhika Khetani) at the Harvard Chan Bioinformatics Core (HBC).

Some materials and hands-on activities were adapted from RNA-seq workflow on the Bioconductor website