# Effect size estimation with apeglm

#### Anqi Zhu, Joseph G. Ibrahim, and Michael I. Love

#### 10/30/2019

Abstract

apeglm provides empirical Bayes shrinkage estimators for effect sizes for a variety of GLM models; apeglm stands for “Approximate Posterior Estimation for GLM”. apeglm package version: 1.8.0

# 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")

# Acknowledgments

Joshua Zitovsky contributed fast C++ code for the beta-binomial likelihood, demonstrated in a later section of the vignette.

We have benefited in the development of `apeglm` from feedback or contributions from the following individuals:

Wolfgang Huber, Cecile Le Sueur, Charlotte Soneson

# 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
    ##   28.53    0.00   28.54

    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
    ##    6.17    0.00    6.17

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
    ##    3.42    0.00    3.42

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.26    0.00    0.27

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.64565685883547  0.0607881068454746   0.229672330995979
    ## ENSG00000000419   6.23452224381112 -0.0849990629106255 -0.0151414731626152
    ## ENSG00000000457   5.45295046335281  0.0601670302479465 -0.0556310344838028
    ## ENSG00000000460   3.77047988584832   0.545203896896376   0.261356100282312
    ## ENSG00000000938 -0.314689549960007   -4.17748981935105  -0.539678670200385
    ## ...                            ...                 ...                 ...
    ## ENSG00000122756  -3.62985171762422    4.03910910257071   -2.00591343292522
    ## ENSG00000122778   5.23905962019254  0.0800614479132839   0.785718241416212
    ## ENSG00000122779   5.66860926529729 -0.0569747440555311  0.0501292974552218
    ## ENSG00000122783       6.0153472139  0.0368918236872043    -0.1578282988443
    ## ENSG00000122786   9.78251386174181  -0.182349284400559  0.0942307382677229
    ##                          cellN61311                dextrt
    ##                           <numeric>             <numeric>
    ## ENSG00000000003  -0.152039407728909    -0.264701484877937
    ## ENSG00000000419 -0.0476622243538528     0.115132140593628
    ## ENSG00000000457  0.0484622685322016   0.00907000392885538
    ## ENSG00000000460   0.349429848473312   -0.0425238102971824
    ## ENSG00000000938   -4.20084641409538  -0.00836574242418845
    ## ...                             ...                   ...
    ## ENSG00000122756    2.95799390215457 -0.000701743016021352
    ## ENSG00000122778   0.240098604814277    -0.285457153643093
    ## ENSG00000122779  -0.044131517147821    -0.109945276882692
    ## ENSG00000122783  -0.297399861535176    0.0690120376646303
    ## ENSG00000122786  -0.120280028801119     0.523514860344871

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
    ##    2.49    0.00    2.49

    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)

The s-value was proposed by Stephens (2016), as a statistic giving the aggregate false sign rate for tests with equal or lower s-value than the one considered. We recommend using a lower threshold on s-values than typically used for adjusted p-values, for example one might be interested in sets with 0.01 or 0.005 aggregate FSR.

    plot(res$padj, fit$svalue, col="blue",
         xlab="DESeq2 padj", ylab="apeglm svalue",
         xlim=c(0,.2), ylim=c(0,.02))

More scrutiny can be applied by using an LFC threshold greater than zero, and asking for the probability of a “false-sign-or-small” (FSOS) event: that the effect size is not further from zero in distance than the threshold amount. We can run the `svalue` function on these per-gene probabilities to produce s-values that bound the FSOS rate for sets of genes. By specifying `threshold=log(2) * 1` above, `apeglm` will then output a vector `thresh` in the results list that gives the per-gene probabilities of false-sign-or-small events.

    s.val <- svalue(fit$thresh)
    cols <- ifelse(s.val < .01, "red", "black")
    plot(res$baseMean, log2(exp(1)) * fit$map[,5], log="x", col=cols,
         xlab=xlab, ylab="LFC")
    abline(h=c(-1,0,1), col=rgb(1,0,0,.5), lwd=2)

# Modeling zero-inflated counts

We have created a separate GitHub repository giving an example of how the apeglm estimator can be used for Zero-Inflated Negative Binomial data. This approach uses the zinbwave method and Bioconductor package to estimate the probability of each zero belonging to the excess zero component. We compare using a Negative Binomial likelihood with the excess zeros down-weighted and using a Zero-Inflated Negative Binomial likelihood. These two approaches with apeglm perform similarly but we note that the first approach involves less additional code and is faster to compute.

https://github.com/mikelove/zinbwave-apeglm

# Modeling ratios of counts

We also show an short example using an alternative likelihood to the negative binomial. Suppose we have allele-specific counts for n=20 vs 20 samples across 5000 genes. We can define a binomial model and test for the allelic balance across groups of samples.

Here we will simulate allele counts from our existing dataset for demonstration. We spike in 10 genes with strong allelic imbalance (instead of an allelic ratio close to 0.5, these will have a ratio of 0.75).

    library(emdbook)
    n <- 20
    f <- factor(rep(1:2,each=n))
    mu <- ifelse(res$baseMean > 50, res$baseMean, 50)
    set.seed(1)
    cts <- matrix(rnbinom(nrow(dds)*2*n,
                          mu=mu,
                          size=1/dispersions(dds)),
                  ncol=2*n)
    theta <- runif(nrow(cts),1,1000)
    prob <- rnorm(nrow(cts),.5,.05) # close to 0.5
    ase.cts <- matrix(rbetabinom(prod(dim(cts)), prob=prob,
                                 size=cts, theta=rep(theta,ncol(cts))),
                      nrow=nrow(cts))
    idx <- 1:10
    idx2 <- which(f == 2)
    theta[idx] <- 1000
    prob[idx] <- 0.75
    # the spiked in genes have an allelic ratio of 0.75
    ase.cts[idx,idx2] <- matrix(rbetabinom(length(idx)*length(idx2), prob=prob[idx],
                                           size=cts[idx,idx2], theta=theta[idx]),
                                nrow=length(idx))

We first need to estimate MLE coefficients and standard errors.

One option to run `apeglm` would be to define a beta-binomial likelihood function which uses the total counts as a parameter, and the logit function as a link. And then this function could be provided to the `log.lik` argument.

    betabinom.log.lik <- function(y, x, beta, param, offset) {
      xbeta <- x %*% beta
      p.hat <- (1+exp(-xbeta))^-1
      dbetabinom(y, prob=p.hat, size=param[-1], theta=param[1], log=TRUE)
    }

However, `apeglm` has faster C++ implementations for the beta-binomial, which were implemented and tested by Joshua Zitovsky. Here, using `method="betabinCR"` is 3 times faster than using the R-defined `log.lik`, and the `"betabinCR"` method also scales significantly better with more samples and more coefficients. As with the negative binomial, `method="betabinC"` can be used if the standard errors are not needed, and this will be 5 times faster than the R-defined `log.lik` approach on this dataset.

The following code performs two iterations of estimating the MLE coefficients, and then estimating the beta-binomial dispersion.

    theta.hat <- 100 # rough initial estimate of dispersion
    x <- model.matrix(~f)
    niter <- 2
    for (i in 1:niter) {
      param <- cbind(theta.hat, cts)
      system.time({
        fit.mle <- apeglm(Y=ase.cts, x=x, log.lik=NULL, param=param,
                          no.shrink=TRUE, log.link=FALSE, method="betabinCR")
      })
      theta.hat <- bbEstDisp(success=ase.cts, size=cts,
                             x=x, beta=fit.mle$map,
                             minDisp=1, maxDisp=500)
    }

We can then plot the MLE estimates over the mean:

    coef <- 2
    xlab <- "mean of normalized counts"
    plot(res$baseMean, fit.mle$map[,coef], log="x", xlab=xlab, ylab="log odds")
    points(res$baseMean[idx], fit.mle$map[idx,coef], col="dodgerblue", cex=3)

Now we run the posterior estimation, including a prior on the second coefficient:

    mle <- cbind(fit.mle$map[,coef], fit.mle$sd[,coef])
    param <- cbind(theta.hat, cts)
    system.time({
      fit2 <- apeglm(Y=ase.cts, x=x, log.lik=NULL, param=param,
                     coef=coef, mle=mle, threshold=0.5,
                     log.link=FALSE, method="betabinCR")
    })

    ##    user  system elapsed
    ##    5.20    0.00    5.22

In the `apeglm` plot, we color in red the genes with a low aggregate probability of false-sign-or-small (FSOS) events (s-value < .01), where we’ve again defined “small” on the log odds scale using the `threshold` argument above.

    par(mfrow=c(1,2))
    ylim <- c(-1,1.5)
    s.val <- svalue(fit2$thresh) # small-or-false-sign value
    plot(res$baseMean, fit.mle$map[,coef], main="MLE",
         log="x", xlab=xlab, ylab="log odds", ylim=ylim)
    points(res$baseMean[idx], fit.mle$map[idx,coef], col="dodgerblue", cex=3)
    abline(h=0,col=rgb(1,0,0,.5))
    cols <- ifelse(s.val < .01, "red", "black")
    plot(res$baseMean, fit2$map[,coef], main="apeglm",
         log="x", xlab=xlab, ylab="log odds", col=cols, ylim=ylim)
    points(res$baseMean[idx], fit2$map[idx,coef], col="dodgerblue", cex=3)
    abline(h=0,col=rgb(1,0,0,.5))

    logit <- function(x) log(x/(1-x))
    logit(.75)

    ## [1] 1.098612

    table(abs(logit(prob[s.val < .01])) > .5)

    ##
    ## TRUE
    ##   10

# Session Info

    sessionInfo()

    ## R version 3.6.1 (2019-07-05)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] emdbook_1.3.11              apeglm_1.8.0
    ##  [3] DESeq2_1.26.0               airway_1.5.3
    ##  [5] SummarizedExperiment_1.16.0 DelayedArray_0.12.0
    ##  [7] BiocParallel_1.20.0         matrixStats_0.55.0
    ##  [9] Biobase_2.46.0              GenomicRanges_1.38.0
    ## [11] GenomeInfoDb_1.22.0         IRanges_2.20.0
    ## [13] S4Vectors_0.24.0            BiocGenerics_0.32.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] bit64_0.9-7            splines_3.6.1          Formula_1.2-3
    ##  [4] assertthat_0.2.1       latticeExtra_0.6-28    blob_1.2.0
    ##  [7] GenomeInfoDbData_1.2.2 yaml_2.2.0             numDeriv_2016.8-1.1
    ## [10] pillar_1.4.2           RSQLite_2.1.2          backports_1.1.5
    ## [13] lattice_0.20-38        glue_1.3.1             bbmle_1.0.20
    ## [16] digest_0.6.22          RColorBrewer_1.1-2     XVector_0.26.0
    ## [19] checkmate_1.9.4        colorspace_1.4-1       plyr_1.8.4
    ## [22] htmltools_0.4.0        Matrix_1.2-17          XML_3.98-1.20
    ## [25] pkgconfig_2.0.3        genefilter_1.68.0      zlibbioc_1.32.0
    ## [28] purrr_0.3.3            xtable_1.8-4           scales_1.0.0
    ## [31] htmlTable_1.13.2       tibble_2.1.3           annotate_1.64.0
    ## [34] ggplot2_3.2.1          nnet_7.3-12            lazyeval_0.2.2
    ## [37] survival_2.44-1.1      magrittr_1.5           crayon_1.3.4
    ## [40] memoise_1.1.0          evaluate_0.14          MASS_7.3-51.4
    ## [43] foreign_0.8-72         tools_3.6.1            data.table_1.12.6
    ## [46] stringr_1.4.0          locfit_1.5-9.1         munsell_0.5.0
    ## [49] cluster_2.1.0          AnnotationDbi_1.48.0   compiler_3.6.1
    ## [52] rlang_0.4.1            grid_3.6.1             RCurl_1.95-4.12
    ## [55] rstudioapi_0.10        htmlwidgets_1.5.1      bitops_1.0-6
    ## [58] base64enc_0.1-3        rmarkdown_1.16         gtable_0.3.0
    ## [61] DBI_1.0.0              R6_2.4.0               gridExtra_2.3
    ## [64] knitr_1.25             dplyr_0.8.3            zeallot_0.1.0
    ## [67] bit_1.1-14             Hmisc_4.2-0            stringi_1.4.3
    ## [70] Rcpp_1.0.2             geneplotter_1.64.0     vctrs_0.2.0
    ## [73] rpart_4.1-15           acepack_1.4.1          coda_0.19-3
    ## [76] tidyselect_0.2.5       xfun_0.10

# References

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology 15 (12):550. https://doi.org/10.1186/s13059-014-0550-8.

Stephens, Matthew. 2016. “False Discovery Rates: A New Deal.” Biostatistics 18 (2). https://doi.org/10.1093/biostatistics/kxw041.

Zhu, Anqi, Joseph G. Ibrahim, and Michael I. Love. 2018. “Heavy-Tailed Prior Distributions for Sequence Count Data: Removing the Noise and Preserving Large Differences.” Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895.