Differential expression analysis with DESeq2

A basic task in the analysis of RNA-seq count data is the detection of differentially expressed genes. The count data are presented as a table which reports, for each sample, the number of sequence fragments that have been assigned to each gene. An important analysis question is the quantification and statistical inference of systematic changes between conditions, as compared to within-condition variability.

We start by loading the DESeq2 package, a very popular method for analysing differential expression of bulk RNA-seq data.

library("DESeq2")
Warning: package 'DESeq2' was built under R version 4.3.3
Warning: package 'GenomeInfoDb' was built under R version 4.3.3
library(tidyverse)

DESeq2 requires count data like that in the SummarizedExperiment we have been working with.

The airway experimental data package contains an example dataset from an RNA-Seq experiment of read counts per gene for airway smooth muscles. These data are stored in a RangedSummarizedExperiment object which contains 8 different experimental samples and assays 64,102 gene transcripts.

library(airway)
data(airway)
se <- airway
se
class: RangedSummarizedExperiment 
dim: 63677 8 
metadata(1): ''
assays(1): counts
rownames(63677): ENSG00000000003 ENSG00000000005 ... ENSG00000273492
  ENSG00000273493
rowData names(10): gene_id gene_name ... seq_coord_system symbol
colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
colData names(9): SampleName cell ... Sample BioSample
rowRanges(se)
GRangesList object of length 63677:
$ENSG00000000003
GRanges object with 17 ranges and 2 metadata columns:
       seqnames            ranges strand |   exon_id       exon_name
          <Rle>         <IRanges>  <Rle> | <integer>     <character>
   [1]        X 99883667-99884983      - |    667145 ENSE00001459322
   [2]        X 99885756-99885863      - |    667146 ENSE00000868868
   [3]        X 99887482-99887565      - |    667147 ENSE00000401072
   [4]        X 99887538-99887565      - |    667148 ENSE00001849132
   [5]        X 99888402-99888536      - |    667149 ENSE00003554016
   ...      ...               ...    ... .       ...             ...
  [13]        X 99890555-99890743      - |    667156 ENSE00003512331
  [14]        X 99891188-99891686      - |    667158 ENSE00001886883
  [15]        X 99891605-99891803      - |    667159 ENSE00001855382
  [16]        X 99891790-99892101      - |    667160 ENSE00001863395
  [17]        X 99894942-99894988      - |    667161 ENSE00001828996
  -------
  seqinfo: 722 sequences (1 circular) from an unspecified genome

...
<63676 more elements>
colData(se)
DataFrame with 8 rows and 9 columns
           SampleName     cell      dex    albut        Run avgLength
             <factor> <factor> <factor> <factor>   <factor> <integer>
SRR1039508 GSM1275862  N61311     untrt    untrt SRR1039508       126
SRR1039509 GSM1275863  N61311     trt      untrt SRR1039509       126
SRR1039512 GSM1275866  N052611    untrt    untrt SRR1039512       126
SRR1039513 GSM1275867  N052611    trt      untrt SRR1039513        87
SRR1039516 GSM1275870  N080611    untrt    untrt SRR1039516       120
SRR1039517 GSM1275871  N080611    trt      untrt SRR1039517       126
SRR1039520 GSM1275874  N061011    untrt    untrt SRR1039520       101
SRR1039521 GSM1275875  N061011    trt      untrt SRR1039521        98
           Experiment    Sample    BioSample
             <factor>  <factor>     <factor>
SRR1039508  SRX384345 SRS508568 SAMN02422669
SRR1039509  SRX384346 SRS508567 SAMN02422675
SRR1039512  SRX384349 SRS508571 SAMN02422678
SRR1039513  SRX384350 SRS508572 SAMN02422670
SRR1039516  SRX384353 SRS508575 SAMN02422682
SRR1039517  SRX384354 SRS508576 SAMN02422673
SRR1039520  SRX384357 SRS508579 SAMN02422683
SRR1039521  SRX384358 SRS508580 SAMN02422677

The package requires count data like that in the SummarizedExperiment we have been working with, in addition to a formula describing the experimental design. We use the cell line as a covariate, and dexamethazone treatment as the main factor that we are interested in.

dds <- DESeqDataSet(se, design = ~ cell + dex)
dds
class: DESeqDataSet 
dim: 63677 8 
metadata(2): '' version
assays(1): counts
rownames(63677): ENSG00000000003 ENSG00000000005 ... ENSG00000273492
  ENSG00000273493
rowData names(10): gene_id gene_name ... seq_coord_system symbol
colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
colData names(9): SampleName cell ... Sample BioSample

The dds object can be manipulated very much like a SummarizedExperiment (in fact: it is a SummarizedExperiment).

There are two reasons which make pre-filtering useful: by removing genes with only few reads across samples, we reduce the size of the dds data object, and thus increase the speed of the transformation and testing functions within DESeq2.

Here we perform a minimal pre-filtering to keep only rows that have at least 10 reads total.

keep <- rowSums(counts(dds)) >= 10
table(keep)
keep
FALSE  TRUE 
41308 22369 
dds <- dds[keep,]

The DESeq workflow is summarized by a single function call, which performs statistical analysis on the data in the dds object.

dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

A table summarizing measures of differential expression can be extracted from the object, and visualized or manipulated using commands we learned earlier.

res <- results(dds)
res
log2 fold change (MLE): dex untrt vs trt 
Wald test p-value: dex untrt vs trt 
DataFrame with 22369 rows and 6 columns
                 baseMean log2FoldChange     lfcSE       stat      pvalue
                <numeric>      <numeric> <numeric>  <numeric>   <numeric>
ENSG00000000003  708.5979      0.3812272 0.1007023   3.785685 1.53286e-04
ENSG00000000419  520.2963     -0.2068403 0.1121077  -1.845013 6.50356e-02
ENSG00000000457  237.1621     -0.0379542 0.1428231  -0.265743 7.90437e-01
ENSG00000000460   57.9324      0.0885314 0.2849344   0.310708 7.56023e-01
ENSG00000000971 5817.3108     -0.4264245 0.0888056  -4.801774 1.57267e-06
...                   ...            ...       ...        ...         ...
ENSG00000273483   2.68955      -0.849208  1.253365 -0.6775424    0.498062
ENSG00000273485   1.28646       0.123614  1.588250  0.0778301    0.937963
ENSG00000273486  15.45244       0.150430  0.482098  0.3120312    0.755017
ENSG00000273487   8.16327      -1.045638  0.693057 -1.5087328    0.131367
ENSG00000273488   8.58437      -0.108945  0.632300 -0.1722990    0.863203
                       padj
                  <numeric>
ENSG00000000003 1.28920e-03
ENSG00000000419 1.94929e-01
ENSG00000000457 9.09901e-01
ENSG00000000460 8.92994e-01
ENSG00000000971 2.06392e-05
...                     ...
ENSG00000273483          NA
ENSG00000273485          NA
ENSG00000273486    0.892518
ENSG00000273487    0.323299
ENSG00000273488    0.943414

Task:

Use the contrast argument of the results function to compare trt vs. untrt groups instead of untrt vs. trt (changes the direction of the fold change).

res_trt_untrt <- results(dds, contrast = c("dex","trt","untrt"))

Volcano plot

A useful illustration of differential expression results is to plot the fold change against the p-value in a volcano plot. This allows to inspect direction and magnitude (fold change) as well as the statistical significance (p-value) of the expression change.

library(ggplot2)
ggplot(as.data.frame(res), 
       aes(x = log2FoldChange, y = -log10(padj))) + 
  geom_point() +
  geom_hline(yintercept = -log10(0.05), col = "red") +
  geom_vline(xintercept = -1, col = "red") +
  geom_vline(xintercept = 1, col = "red")
Warning: Removed 4337 rows containing missing values or values outside the scale range
(`geom_point()`).

We can get more advanced with creating volcano plots, labeling cells with things like the

library(ggrepel) #This is a good library for displaying text without overlap
library(biomaRt, quietly = TRUE) #for ID mapping

res_df <- as.data.frame(res)
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
genes <- rownames(res_df)
gene_map <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id","hgnc_symbol"),values=genes,mart= mart)
ind <- match(rownames(res_df), gene_map$ensembl_gene_id)
res_df$gene <- gene_map$hgnc_symbol[ind]

res_df <- mutate(res_df, sig = ((padj < 0.05) & abs(log2FoldChange) > 2)) 
  ggplot(res_df, aes(x = log2FoldChange, y = -log10(padj), col=sig)) +
       geom_point() +
       geom_vline(xintercept=c(-2, 2), col="red") +
       geom_hline(yintercept=-log10(0.05), col="red") +
       geom_text_repel(data=filter(res_df, sig), aes(label=gene))
Warning: Removed 4337 rows containing missing values or values outside the scale range
(`geom_point()`).
Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text_repel()`).
Warning: ggrepel: 209 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

MA plot

Another useful illustration of differential expression results is to plot the fold changes as a function of the mean of the expression level (normalized counts) across samples in an MA plot.

Points will be colored if the adjusted p-value is less than a defined significance threshold (default: 0.1). Points which fall out of the window are plotted as open triangles pointing either up or down.

plotMA(res)

The DESeq2 vignette also describes several other useful result exploration and data quality assessment plots.

Exercises

Basic

  1. Use the coef function to examine the actual coefficients of the model.

  2. Get the number of genes considered significantly expressed at the alpha level of 0.1 (for the adjusted p value).

  3. Now see how many genes would be considered differentially expressed at an alpha level of 0.05 and a log2 fold change cutoff of at least 1. Note: Take a look again at the arguments of the results function. Are there any you should change?

#1.  Use the `coef` function to examine the actual coefficients of the model.
coef(dds)

#2.  
sum(res$padj<0.1, na.rm = TRUE)

#3. 
res_05_1 <- results(dds, alpha=0.05,lfcThreshold = 1)
sum((res_05_1$padj<0.1)&(res_05_1$log2FoldChange>=1), na.rm = TRUE)
Advanced

Let’s imagine that, instead of all being untreated, half of the samples had been treated with albuterol:

fake_se <- se
#Currently the factor only has the untrt level, so we need to add another
levels(colData(fake_se)$albut) <- c(levels(colData(fake_se)$albut), "trt")
colData(fake_se)$albut[c(1,3,4,8)] <- "trt"

Remake the dds object such that the albut column is an additional covartiate in the experimental design.

If there is time, compare the number of significant results (for comparing cell and dex) when albut is and is not accounted for. Does it make a difference?

dds_fake <- DESeqDataSet(fake_se, design = ~ albut + cell + dex)
keep <- rowSums(counts(dds_fake)) >= 4
dds_fake <- dds_fake[keep,]
dds_fake <- DESeq(dds_fake)
res_fake <- results(dds_fake)

sum(res$padj<0.1, na.rm = TRUE)
sum(res_fake$padj<0.1, na.rm = TRUE)
Bonus

Try installing the Glimma package from bioconductor. Use it to create an interactive multidimensional scaling (MDS) plot of your results.

library(Glimma)
glimmaMDS(dds)

This lesson was adapted from materials created by Ludwig Geistlinger