Skip to content
Bench to Insights
Blog

Differential Expression

Differential expression asks, for each gene, whether its count differs between conditions by more than the biological and technical variation would explain on its own. DESeq2 answers with a negative binomial model fitted to each gene. The steps it runs to get there explain the numbers in the results table.

What DESeq2 does

Section titled “What DESeq2 does”
  1. Size factor normalization accounts for sequencing depth between samples.
  2. Dispersion estimation models the gene-wise variance, shrinking each gene’s estimate toward a fit from all genes.
  3. A generalized linear model is fit for each gene.
  4. A Wald test gives a p-value for the coefficient of interest.
  5. Log fold change shrinkage pulls noisy fold changes from low-count genes toward zero.

Shrinkage is the step that changes the numbers you read. A gene with 5 reads in one condition and 50 in the other has a large raw fold change and a wide confidence interval, and shrinkage moves the estimate toward zero. The shrunken fold change is the one to plot and to rank on.

Run DESeq2 on the airway data

Section titled “Run DESeq2 on the airway data”

The model ~ cell + dex estimates the treatment effect after accounting for the donor cell line. With the untreated level as the reference, the coefficient dex_trt_vs_untrt is the effect of treatment.

library(DESeq2)
library(airway)


data(airway)
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds$dex <- relevel(dds$dex, ref = "untrt")


# Filter low-count genes.
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]


# Run the full pipeline: size factors, dispersion, model, test.
dds <- DESeq(dds)

Extract results

Section titled “Extract results”

results() returns the table. lfcShrink() takes the fitted object and returns the same table with the log fold change column replaced by shrunken estimates. Those are the ones to use for ranking and plotting.

# Shrunken results for the treatment coefficient.
res <- lfcShrink(dds, coef = "dex_trt_vs_untrt", type = "apeglm")


# The first rows.
head(res)


# Counts of up, down, and unchanged genes at padj < 0.1.
summary(res)
log2 fold change (MAP): dex trt vs untrt
Wald test p-value: dex trt vs untrt
DataFrame with 6 rows and 5 columns
                 baseMean log2FoldChange     lfcSE      pvalue        padj
                <numeric>      <numeric> <numeric>   <numeric>   <numeric>
ENSG00000000003  708.5979     -0.3640838 0.1000795 1.53286e-04 1.28920e-03
ENSG00000000419  520.2963      0.1864336 0.1081165 6.50354e-02 1.94930e-01
ENSG00000000457  237.1621      0.0310729 0.1300506 7.90437e-01 9.09900e-01
ENSG00000000460   57.9324     -0.0471496 0.2123688 7.56024e-01 8.92994e-01
ENSG00000000971 5817.3108      0.4025264 0.0886677 1.57266e-06 2.06391e-05
ENSG00000001036 1282.1007     -0.2279666 0.0873527 6.75809e-03 3.34601e-02


out of 22369 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 2610, 12%
LFC < 0 (down)     : 2224, 9.9%
outliers [1]       : 0, 0%
low counts [2]     : 4337, 19%
(mean count < 5)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

The summary counts genes with an adjusted p-value below 0.1, split by the direction of the fold change. The 0.1 is DESeq2’s default, a choice rather than a truth. The visualization page uses 0.05 for the figures. For your own study, choose the threshold from the cost of a false positive in the downstream experiment.

The results table

Section titled “The results table”

Each row is a gene. The columns are the base mean, the log fold change, its standard error, the Wald p-value, and the adjusted p-value.

# Convert to a data frame for inspection.
res_df <- as.data.frame(res)
head(res_df[order(res_df$padj), ], 5)
                  baseMean log2FoldChange     lfcSE        pvalue          padj
ENSG00000152583   997.4447       4.559848 0.1858830 4.110667e-136 7.412355e-132
ENSG00000165995   495.0957       3.280119 0.1332056 4.463384e-135 4.024187e-131
ENSG00000120129  3409.0384       2.936591 0.1225024 3.033840e-129 1.823540e-125
ENSG00000101347 12703.4128       3.754342 0.1576101 7.682657e-129 3.463342e-125
ENSG00000189221  2341.7807       3.341300 0.1432340 5.212706e-123 1.879910e-119

The adjusted p-value corrects for testing every gene at once. With 22369 genes under test, a raw p-value of 0.01 is expected by chance in hundreds of them. The adjusted value, Benjamini-Hochberg by default, is the one to filter on.

The MA plot

Section titled “The MA plot”

The MA plot puts the log fold change against the mean expression, and it is where the effect of shrinkage shows. The low-count genes on the left had extreme raw fold changes. After shrinkage they sit close to zero.

# MA plot data: mean expression against log fold change.
ma_data <- as.data.frame(res)
head(ma_data[, c("baseMean", "log2FoldChange", "padj")], 3)
                baseMean log2FoldChange        padj
ENSG00000000003 708.5979     -0.3640838 0.001289205
ENSG00000000419 520.2963      0.1864336 0.194929522
ENSG00000000457 237.1621      0.0310729 0.909899502

The visualization page draws the MA and volcano plots as figures. The edgeR and limma page repeats the comparison with two other packages, and the numbers agree.