Typical RNA-seq call from DESeq2

Note: the typical RNA-seq workflow for users would be to call apeglm estimation from within the lfcShrink function from the DESeq2 package. The unevaluated code chunk shows how to obtain apeglm shrinkage estimates after running DESeq. See the DESeq2 vignette for more details. The lfcShrink wrapper function takes care of many details below, and unifies the interface for multiple shrinkage estimators. The coefficient to shrink can be specified either by name or by number (following the order in resultsNames(dds)). Be aware that DESeq2’s lfcShrink interface provides LFCs on the log2 scale, while apeglm provides coefficients on the natural log scale.

res <- lfcShrink(dds, coef=2, type="apeglm")

Example RNA-seq analysis

Here we show example code which mimics what will happen inside the lfcShrink function when using the apeglm method (Zhu, Ibrahim, and Love 2018).

Load a prepared SummarizedExperiment:

library(airway)
data(airway)
head(assay(airway))

##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003        679        448        873        408       1138
## ENSG00000000005          0          0          0          0          0
## ENSG00000000419        467        515        621        365        587
## ENSG00000000457        260        211        263        164        245
## ENSG00000000460         60         55         40         35         78
## ENSG00000000938          0          0          2          0          1
##                 SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003       1047        770        572
## ENSG00000000005          0          0          0
## ENSG00000000419        799        417        508
## ENSG00000000457        331        233        229
## ENSG00000000460         63         76         60
## ENSG00000000938          0          0          0

For demonstration, we will use 5000 genes of the airway dataset, the first from those genes with at least 10 counts across all samples.

keep <- head(which(rowSums(assay(airway)) > 0), 5000)
airway <- airway[keep,]

First run a DESeq2 differential expression analysis (Love, Huber, and Anders 2014) (size factors, and dispersion estimates could similarly be estimated using edgeR):

library(DESeq2)
dds <- DESeqDataSet(airway, ~cell + dex)
dds$dex <- relevel(dds$dex, "untrt")
dds <- DESeq(dds)
res <- results(dds)

Defining data and parameter objects necessary for apeglm. We must multiply the coefficients from DESeq2 by a factor, because apeglm provides natural log coefficients. Again, this would be handled inside of lfcShrink in DESeq2 for a typical RNA-seq analysis.

x <- model.matrix(design(dds), colData(dds))
param <- dispersions(dds)
mle <- log(2) * cbind(res$log2FoldChange, res$lfcSE)
offset <- matrix(log(sizeFactors(dds)),
                 ncol=ncol(dds),
                 nrow=nrow(dds),byrow=TRUE)

Running apeglm

Here apeglm on 5000 genes takes less than a minute on a laptop. It scales with number of genes, the number of samples and the number of variables in the design formula, where here we have 5 coefficients (one for the four cell cultures and one for the difference due to dexamethasone treatment).

We provide apeglm with the SummarizedExperiment although the function can also run on a matrix of counts or other observed data. We specify a coef as well as a threshold which we discuss below. Note that we multiple the threshold by log(2) to convert from log2 scale to natural log scale.

library(apeglm)
system.time({
  fit <- apeglm(Y=airway, x=x, log.lik=logLikNB, param=param, coef=ncol(x),
                threshold=log(2) * 1, mle=mle, offset=offset)
})

##    user  system elapsed 
##  24.816   0.000  24.838

str(fit$prior.control)

## List of 7
##  $ no.shrink            : int [1:4] 1 2 3 4
##  $ prior.mean           : num 0
##  $ prior.scale          : num 0.299
##  $ prior.df             : num 1
##  $ prior.no.shrink.mean : num 0
##  $ prior.no.shrink.scale: num 15
##  $ prior.var            : num 0.0894

There are faster implementations of apeglm specifically for negative binomial likelihoods. The version nbinomR is ~5 times faster than the default general.

system.time({
  fitR <- apeglm(Y=assay(airway), x=x, log.lik=NULL, param=param, coef=ncol(x),
                 threshold=log(2) * 1, mle=mle, offset=offset, method="nbinomR")
})

##    user  system elapsed 
##   5.356   0.000   5.362

The version nbinomCR is ~10 times faster than the default general.

system.time({
  fitCR <- apeglm(Y=assay(airway), x=x, log.lik=NULL, param=param, coef=ncol(x),
                 threshold=log(2) * 1, mle=mle, offset=offset, method="nbinomCR")
})

##    user  system elapsed 
##   2.788   0.004   2.793

The version nbinomC returns only the MAP coefficients and can be ~50-100 times faster than the default general. The MAP coefficients are the same as returned by nbinomCR above, we just skip the calculation of posterior SD. A variant of nbinomC is nbinomC* which includes random starts.

system.time({
  fitC <- apeglm(Y=assay(airway), x=x, log.lik=NULL, param=param, coef=ncol(x),
                 threshold=log(2) * 1, mle=mle, offset=offset, method="nbinomC")
})

##    user  system elapsed 
##   0.188   0.000   0.189

Among other output, we have the estimated coefficients attached to the ranges of the SummarizedExperiment used as input:

class(fit$ranges)

## [1] "CompressedGRangesList"
## attr(,"package")
## [1] "GenomicRanges"

mcols(fit$ranges, use.names=TRUE)

## DataFrame with 5000 rows and 5 columns
##                       X.Intercept.         cellN061011         cellN080611
##                          <numeric>           <numeric>           <numeric>
## ENSG00000000003     6.645656858832  0.0607881068431617   0.229672330993859
## ENSG00000000419   6.23452224381182  -0.084999062909147 -0.0151414731605632
## ENSG00000000457   5.45295046291592   0.060167029688361 -0.0556310344682003
## ENSG00000000460   3.77047988584934   0.545203896906746   0.261356100290922
## ENSG00000000938 -0.314689549938794   -4.17748981892207  -0.539678670081039
## ...                            ...                 ...                 ...
## ENSG00000122756  -3.62985171699038    4.03910910191739   -2.00591343576517
## ENSG00000122778   5.23905962019645  0.0800614479140051    0.78571824141725
## ENSG00000122779   5.66860926529946 -0.0569747440556599  0.0501292974549099
## ENSG00000122783   6.01534721390051  0.0368918236863277  -0.157828298844423
## ENSG00000122786   9.78251386175167   -0.18234928440018  0.0942307382694629
##                          cellN61311                dextrt
##                           <numeric>             <numeric>
## ENSG00000000003  -0.152039407731068    -0.264701484867069
## ENSG00000000419 -0.0476622243521883     0.115132140591595
## ENSG00000000457  0.0484622689335007   0.00907000351291652
## ENSG00000000460   0.349429848485217   -0.0425238103250921
## ENSG00000000938   -4.20084641372255  -0.00836574242756413
## ...                             ...                   ...
## ENSG00000122756     2.9579939014709 -0.000701743005250805
## ENSG00000122778   0.240098604814436    -0.285457153652284
## ENSG00000122779 -0.0441315171481841    -0.109945276886632
## ENSG00000122783   -0.29739986153585    0.0690120376649482
## ENSG00000122786   -0.12028002880088     0.523514860323912

We can compare the coefficients from apeglm with those from DESeq2. apeglm provides coefficients on the natural log scale, so we must convert to log2 scale by multiplying by log2(exp(1)). Note that DESeq2’s lfcShrink function converts apeglm coefficients to the log2 scale internally.

system.time({
  res.shr <- lfcShrink(dds, coef=5)
})

##    user  system elapsed 
##   1.920   0.004   1.930

DESeq2.lfc <- res.shr$log2FoldChange
apeglm.lfc <- log2(exp(1)) * fit$map[,5]

Here we plot apeglm estimators against DESeq2:

plot(DESeq2.lfc, apeglm.lfc)
abline(0,1)

Here we plot MLE, DESeq2 and apeglm estimators against the mean of normalized counts:

par(mfrow=c(1,3))
lims <- c(-8,8)
hline <- function() abline(h=c(-4:4 * 2),col=rgb(0,0,0,.2))
xlab <- "mean of normalized counts"
plot(res$baseMean, res$log2FoldChange, log="x",
     ylim=lims, main="MLE", xlab=xlab)
hline()
plot(res$baseMean, DESeq2.lfc, log="x",
     ylim=lims, main="DESeq2", xlab=xlab)
hline()
plot(res$baseMean, apeglm.lfc, log="x",
     ylim=lims, main="apeglm", xlab=xlab)
hline()

Specific coefficients

Note that p-values and FSR define different events, and are not on the same scale. An FSR of 0.5 means that the estimated sign is as bad as random guess.

par(mfrow=c(1,2),mar=c(5,5,1,1))
plot(res$pvalue, fit$fsr, col="blue",
     xlab="DESeq2 pvalue", ylab="apeglm local FSR",
     xlim=c(0,1), ylim=c(0,.5))
abline(0,1)
plot(-log10(res$pvalue), -log10(fit$fsr),
     xlab="-log10 DESeq2 pvalue", ylab="-log10 apeglm local FSR",
     col="blue")
abline(0,1)