Fast effect size shrinkage software for beta-binomial models of allelic imbalance

Estimation methods

We evaluated three estimation methods on their ability to estimate allelic expression proportions (or equivalently, the effects of covariates on allelic expression proportions): maximum likelihood estimation (ML estimation or MLE) with the likelihood described below, approximate posterior estimation of GLM coefficients (apeglm) and adaptive shrinkage (ash). All analyses was done using R version 3.5.1¹². The first two methods mentioned are implemented in the apeglm v.1.7.5 package, while the last is implemented in the ashr v.2.2.32 package. When using the ash function in the latter package, we set the method parameter equal to "shrink". While there are many Bayesian estimation methods that can be used to quantify allelic imbalance, these allow for arbitrary design matrices. For instance, these methods can estimate differences in AI between groups while controlling for, or allowing interactions with, multiple additional variables, and can estimate the effects of continuous variables on AI.

For the g-th gene (1 ≤ g ≤ G), a beta-binomial GLM was fit to model allele-specific counts as follows. Let Y_ig be the read counts of the first of the two alleles (which allele is designated as the first allele is arbitrary) for the i-th subject, 1 ≤ i ≤ I. Investigators may designate the first and second alleles of a gene as the paternal and maternal alleles or as the alternate and reference alleles. It is assumed that Y_ig ∼ BetaBin(n_ig , p_ig ,ϕ_g), where n_ig is the equal to the total counts of both alleles for the i-th subject, p_ig is the probability of counts belonging to the first allele of the i-th subject, and ϕ_g is the overdispersion parameter. For the remainder of this paper, we will refer to the total allele-specific counts for both alleles of a particular gene and for a particular sample as the ‘total counts’ for that gene and sample. Furthermore, we will refer to the probability that counts for a particular gene belong to a particular allele for a particular sample as the ‘allelic proportion’ for that particular allele and sample. In this case, ϕ → ∞ implies no overdispersion beyond what would be seen in a binomial distribution and ϕ → 0 implies increasing variance. n_1g , ..., n_Ig are assumed to be fixed and known. As the beta-binomial probability density function has multiple forms and parameterizations, we specify our parametrization as:

where B specifies the beta function. Furthermore, let x_i be the i-th row of the design matrix X (matrix where columns are vectors of covariates of interest). Potential predictors include disease status for association studies, parent of origin for imprinting studies, and the presence of a SNP for eQTL linkage studies. We also assume that $p_{ig}={[1+\exp (-{\mathbf{\text{x}}}_i^T{\mathit{\beta }}_g)]}^{-1}$, or equivalently $\text{logit}(p_{ig})={\mathbf{\text{x}}}_i^T{\mathit{\beta }}_g,$ where β_g = (β_1g , ...,β_Kg)^T is a vector of coefficients representing the effect sizes for the predictors in the design matrix. For ML estimation, β_g is estimated via maximum likelihood. Constrained optimization is used for the nuisance parameter ϕ_g with a maximum of 500, so that genes with no overdispersion have finite estimated values of ϕ.

Apeglm additionally assumes a zero-centered Cauchy prior distribution for the effects of one of the predictors¹¹. For estimating the effect of the j-th predictor in our model, where 1 ≤ j ≤ K is chosen by the user, and for the g-th gene, we have:

Apeglm shrinks the effect of one chosen predictor at a time, across all genes. The scale parameter of the Cauchy prior, γ_j, is estimated by pooling information across genes. The posterior distribution of β_g is the product of the above Cauchy prior and beta-binomial likelihood, and apeglm provides Bayesian shrinkage estimates based on the mode of the posterior as well as standard errors. Genes with lower expression, smaller numbers of heterozygous subjects and higher dispersion in allelic proportions will have flatter likelihoods, which will lead to the prior having more influence and shrinkage being greater. Furthermore, if the ML estimates are tightly clustered about zero, the estimated scale parameter of the Cauchy prior will be smaller. This will lead to more peakedness in the prior and also cause shrinkage to be greater.

The original apeglm package estimated regression coefficients using C++ for negative binomial GLMs, while GLMs with other likelihoods, such as the beta-binomial, were fit completely in R. To improve scalability for large data sets with beta-binomial GLMs, we wrote fast C++ code for calculating maximum likelihood and apeglm shrinkage estimates of beta-binomial regression coefficients. We also changed the source code to speed up computation of the standard errors (though such computations were still done in R) and prevent convergence issues. Details can be found in the Supplementary Methods section¹³.

Ash is a general Empirical Bayes shrinkage estimator for hypothesis testing and measuring uncertainty in a vector of effects of interest, such as a set of log fold changes in gene expression between biological conditions¹⁰. Suppose again that one is interested in the effect sizes of the j-th predictor, β_.j = (β_j1, ..., β_jG), where 1 ≤ j ≤ K. Ash takes as input a vector of ML estimated effects ${\widehat{\mathit{\beta }}}_{\mathit{.j}}=({\widehat{\beta }}_{j1},\mathrm{...},{\widehat{\beta }}_{jG})$ and corresponding estimated standard errors σ_β.j = (σ_{β j1} , ..., σ_{β jG}). Here we take the estimated standard errors to be the true standard errors as suggested in the original methodology for ash, though the developers of ash have recently proposed an extension to their method that allows for random errors¹⁴. For all 1 ≤ g ≤ G, it is assumed that ${\widehat{\beta }}_{jg}|{\beta }_{jg}\sim N({\beta }_{jg},{\sigma }_{\beta jg})$ and that β_jg ~ h_j , where h_j is some unimodal, zero-mode prior distribution. h_j is estimated from the ML estimates using mixtures of uniforms and a point-mass at zero, a choice guided by the author’s claim that any unimodal distribution can be approximated as a mixture of uniforms with arbitrary accuracy. The posterior is ${\beta }_{jg}|{\widehat{\beta }}_{jg}\sim N({\beta }_{jg},{\sigma }_{\beta jg})\times h_j,$ and ash provides Empirical Bayes shrinkage estimates using the mean of the posterior as well as standard errors. Genes with larger standard errors for their ML estimates will have a flatter likelihood that will be less impactful on the estimation. Thus, estimates for these genes will be shrunk more. Like apeglm, ash can only shrink estimates for one covariate at a time.

Datasets and simulations

We compared the three estimation methods using the data set from the allelic expression study by Crowley et al.^15,16 The study took mice from three divergent inbred strains (CAST/EiJ, PWK/PhJ and WSB/EiJ) and performed a diallel cross. The data set contains ASE counts for 72 mice and 23,297 genes in the resulting cross, with 12 mice of each possible parent combination (e.g. CAST/EiJ as mother and PWK/PhJ as father is one parent combination, and PWK/PhJ as mother and CAST/EiJ as father is another), and an equal number of males and females within each parent combination. Sequencing was performed with the Illumina HiSeq 2000 platform to generate 100-bp paired-end reads and following the TruSeq RNA Sample Preparation v2 protocol. To assure that the mice all had the same alleles, we chose one genotype to focus on, namely the genotype resulting from the cross with CAST/EiJ and PWK/PhJ. The resulting data set, which we will refer to for the remainder of this paper as the ‘mouse data set’, had 24 mice, 12 of each sex and 12 of each parent of origin, and each mouse had nearly the same nuclear genetic composition as a result of the cross.

To evaluate the estimators on estimating effect sizes of predictors when the truth is known, we first fit an intercept-only beta-binomial model on each gene for the mouse data set. ϕ = [ϕ_g] is the vector of ML estimates of the overdispersion parameter from each model, and µ = [µ_g] is the vector of ML estimates of allelic proportions (which ranges between 0 and 1). 8 mice were then selected from the data set. Denote N_I×G = [n_ig] as the matrix of total ASE counts for the 8 mice. Finally, a matrix of counts from one of the alleles Y_I×G = [y_ig] was simulated for a sample size of 4 vs. 4, where y_ig was simulated from BetaBin(n_ig, p_ig, ϕ_g), logit(p_ig) = µ_g + β_gx, β_g was simulated from a standard normal distribution, and x splits the mice into two groups of size four (x = 1 if a mouse is in the first group and 0 otherwise). Samples were drawn from the beta-binomial distribution using the emdbook v1.3.11 package¹⁷. We refer to this simulation throughout the paper as the ‘standard normal simulation’, reflecting the distribution of the true effect sizes.

A second simulation was also performed that was similar in setup to the first, but with modifications to the distribution of β_g and ϕ_g. In many studies, the effect sizes of a predictor will be zero for all but a handful of genes. Thus, β_g was simulated from t₃/10 (a Student’s t-distribution with 3 degrees of freedom scaled by 1/10), which gave effects mostly close to zero, but with moderate and large effects occasionally appearing (Supplementary Figure 1¹³). Furthermore, the distribution of ϕ_g from the mouse data appeared to be a mixture of two distributions: Genes without overdispersion had an obvious point mass at 500 with 70% proportion, and the remaining 30% genes had a distribution somewhat similar to an exponential with a mean of β = 179 (Supplementary Figure 2¹³). To get more over-dispersed allele-specific counts, ϕ_g was simulated from 0.5Exp(β=89) + 0.5(500), a mixture distribution where one component was exponential with a mean of 89 and had 50% proportion, and the other component was a point mass at 500 and had 50% proportion. We refer to this simulation throughout the paper as the ‘Student’s t simulation’, again reflecting the distribution of the true effect sizes. Note that these two simulations assume a data generating process, specifically the same data generating process as our assumed likelihood.

The estimators were then evaluated on real data with the focus on estimating mean, or gene-wide, allelic imbalance. From the mouse data set, random samples of size 6 were drawn, and this process was repeated 100 times. We will refer to these samples throughout the paper as the ‘random subsamples’. For each random subsample, the ML, apeglm and ash estimates of intercept-only models were calculated for the genes (where the intercept term was shrunk), and the MLE of the held-out 18 mice was taken to be the truth. Estimating the intercept in an intercept-only model for each gene is equivalent to estimating overall allelic imbalance for each gene.

Additional simulations were conducted for evaluating computational performance of our improvements to apeglm, to see how well they would scale to larger and more complicated data sets. Allele-specific counts were simulated in a similar manner as the apeglm vignette¹⁸. Briefly, we have Y_100×5000 = [y_ig] as our simulated count matrix for one allele with associated total count matrix N_100×5000 = [n_ig] where rows are samples and columns are genes, y_ig ~ BetaBin(n_ig, p_g, θ_g), θ_g ~ U (0, 1000), p_g ~ N (.5, 0.5²), n_ig ~ NB(µ_g, 1/ϕ_g), and µ_g, ϕ_g are based on the airway data set by Himes et al.¹⁹ To see how well our improvements scaled with increasing numbers of covariates, the data were split multiple times into differing numbers of groups of approximately equal size, where the number of groups ranged from 2 to 10. With K groups, the design matrix was X_100×K = [1 x₁ ... x_{K –1}], where x_j is an indicator variable for the (j + 1)-th group, or a row vector whose i-th element is 1 if the i-th sample is in the (j + 1)-th group and 0 otherwise. A simulation was also conducted to see how well apeglm would work with continuous predictors. This time, Y and N was kept the same, but with the design matrix X_100×4 = (1, x₁, x₂, x₃) = [x_ij], where x₁ = (1, 0, 1..., 1, 0)^T separates the samples into two equally sized groups and x_i2, x_i3 ~ N (0, 1). x₁ is the covariate whose effect size estimates are shrunk.

Data processing

Genes where at least three samples did not have at least 10 counts were removed, which we considered minimal filtering that shouldn’t decrease statistical power. Genes without at least one count for both alleles across all individuals were removed. Genes with a marginally significant sex or parent effect were removed, so that all samples could be assumed independent and identically distributed for all genes. Genes were removed from the mouse data set prior to conducting random sampling from the data set or simulations.

To determine whether sex or parent effects were significant, beta-binomial GLMs were estimated for each gene by maximum likelihood, with a design matrix that included a sex effect (an indicator that was 1 if male and 0 if female), a parent-of-origin effect (an indicator that was 1 if the mother was the CAST/EiJ strain and 0 if the father was the CAST/EiJ strain) and an interaction term. For each gene, if the p-value for the sex, parent-of-origin or interaction effect was less than 0.1, the effect was deemed marginally significant for that gene.

Technical details of evaluations

For each gene, we define the shrinkage score as movement from MLE to zero. We define a gene as (noticeably) shrunk if shrinkage exceeds 0.1, and substantially or most shrunk if shrinkage is greater than max $(1,|{\widehat{\beta }}_{\text{MLE}}|/4).$ For instance, if an apeglm estimate for a gene is 0.15 closer to zero than the MLE, then the shrinkage score is 0.15 and the gene is noticeably shrunk but not substantially shrunk by apeglm.

Concordance at the top (CAT) plots²⁰ were used to determine which estimation method could best find the most important genes (the genes with the greatest allelic imbalance or largest effect size). For an estimation method, concordance at the top takes the top genes according to the true ranking and compares it to the top genes according to the estimates, where the top genes are the genes with the largest true or estimated effect sizes in absolute value. For instance, a concordance at the top 10 of 90% means that the top 10 genes according to the estimation method and the top 10 genes according to the truth agree for 9 out of 10 genes.

For evaluating the performance of the three methods in estimating intervals, we calculated normality-based 95% confidence and credible intervals (both of which we will abbreviate as CIs) of the ML and apeglm estimators using their standard errors, or intervals based on the Laplace approximation of the likelihood and posterior. Such normality-based intervals are the default and suggested method for the apeglm package. Credible intervals in the ashr package were calculated from directly estimating tail probabilities of the posterior.

For each of the design matrices posited in our computation simulation, computational performance of apeglm estimation was compared between the old and new apeglm code. From apeglm v1.7.5, we set the method parameter equal to “betabinCR" to run the new C++ code, and set the log.lik parameter equal to a beta-binomial log-likelihood function to run the old code from before our improvements were introduced (version 1.6.0 of the package). Details can be found in the vignette¹⁸. Computational performance of ML estimation was also compared between our improved apeglm package and the following packages: aod v1.3.1²¹, VGAM v1.1²², aods3 v0.4²³, gamlss v5.1²⁴ and HRQoL v1.0²⁵. Computational performance was evaluated using the microbenchmark v1.4.6 package²⁶ for estimation of a single gene and elapsed time for estimation of all 5000 genes, on a 2012 15-inch MacBook Pro with an Intel Core i7-3720QM processor.

Determining the optimal filtering rule

In addition to comparing the three estimation methods described above, maximum likelihood estimation paired with optimal filtering criteria was also assessed via concordance at the top. CAT was chosen over other benchmark metrics, such as mean absolute error, as the different number of genes after filtering would make comparisons between filtered MLE and the three unfiltered methods biased. Furthermore, as we were primarily interested in whether a good filtering rule even existed, the true ranking of genes was used to determine the filtering rule. We looked at three rules: 1) removing genes where less than half the samples had a minimum total count threshold, 2) removing genes where less than all the samples had a minimum total count threshold, and 3) removing genes where the sum of total counts across samples was less than a certain threshold. For the remainder of the paper, we will refer to the sum of total counts across samples as the ‘summed counts’ of a gene. For each rule, various different thresholds were looked at: {0, 10, ..., 200} were potential thresholds for rule 1, {0, 10, ..., 100} were potential thresholds for rule 2, and {0, 50, ..., 1000} were potential thresholds for rule 3. For each rule and threshold, the MLE was calculated and concordance among the top 50, 100, 200, 300, 400 and 500 genes were averaged. We will refer to the rule and threshold that had the best concordance as the ‘optimal filtering rule’.

Standard normal simulation

We began by looking at a simulation where allelic counts came from known beta-binomial distributions and effect sizes came from a standard normal distribution. In this simulation, apeglm and ash successfully shrunk erroneously large estimates and reduced estimation error, particularly for genes that were noticeably shrunk (see Table 1 and Figure 1).

Figure 1. Truth vs. estimate and CAT Plots for normal simulation.

a) truth vs. estimate plot for MLE. Blue points represent genes substantially shrunk by apeglm only, orange points represent genes substantially shrunk by ash only and green points represent genes substantially shrunk by both ash and apeglm. b) truth vs. estimate plots for apeglm. c) truth vs. estimate plots for ash. d) CAT plot for the three methods as well as for MLE after filtering. CAT: Concordance at the Top, MLE: Maximum Likelihood Estimation, apeglm: Approximate Posterior Estimation of Generalized Linear Model Coefficients, ash: Adaptive Shrinkage.

All three estimation methods gave similar mean absolute error (MAE), as many genes did not differ much between the methods (Table 1). In exploring the behavior of shrinkage estimators, we were most interested in genes where shrinkage was high, and thus where estimates would be much closer to or much farther from the truth for one estimation method than for another. Thus, in addition to overall MAE, we also calculated MAE among genes that were noticeably shrunk by apeglm and genes that were noticeably shrunk by ash, to determine whether there was substantial improvement on average when apeglm or ash did noticeably shrink a gene. Among genes that were shrunk by apeglm, apeglm decreased the mean absolute error by 18.1%, and among genes that were shrunk by ash, ash decreased the mean absolute error by 21.1%. Moreover, from Figure 1a–c, it can be seen that apeglm shrunk most ML estimates that were inflated, bringing them closer to the truth, and mostly left truly large effects alone. Ash also shrunk ML estimates that were inflated, including some inflated estimates missed by apeglm. However, ash also had a tendency to incorrectly and excessively shrink: some genes with estimates close to the truth were severely shrunk, and several genes with truly large effects were shrunk to zero. Because of this tendency to over-shrink, ash performed worse among genes with large effects than among genes with small effects. For instance, among genes with effect sizes greater than two in absolute value, ash estimates had a higher mean absolute error than the MLE.

Ash and apeglm also performed better than the MLE in determining the most important genes, where concordance at the top was higher regardless of the number of genes being considered (Figure 1d). Apeglm performed slightly better than ash in concordance at the top 100 genes, but otherwise they performed about the same. Concordance at the top for the MLE was optimized when filtering out genes with summed counts less than 350. Using this filtering, we were able to get CAT results better than that of the shrinkage estimates, even if only by a very small amount. Thus, for this simulation, it was possible to outperform both apeglm and ash with filtering alone (provided that the true ranking of genes was known, and used to determine the optimal filtering rule).

With regard to the extent of shrinkage, both apeglm and ash mainly exhibited shrinkage for genes that had very low counts (Supplementary Figure 3¹³). This is not too surprising for this particular simulation, as after filtering out lowly-expressed genes, the remaining ML estimates were much closer to the truth (Supplementary Figure 4¹³). When comparing shrinkage scores between apeglm and ash, we found that there was a clear upward shift of shrinkage scores for ash (Supplementary Table 1¹³), further showing that ash had more extreme shrinkage than apeglm for this dataset. Though all three methods gave intervals that were similar in coverage probability, average interval width was smaller for apeglm and ash compared to the MLE (Table 1).

Student’s t Simulation

We also investigated the performance of the estimators when most of the effect sizes were close to zero and overdispersion was large. Here the shrinkage estimates had even more marked improvement over the ML estimates (see Table 2 and Figure 2).

Figure 2. Truth vs. estimate and CAT Plots for Student’s t Simulation.

a) truth vs. estimate plot for MLE. Orange points represent genes substantially shrunk by ash only and green points represent genes substantially shrunk by both ash and apeglm. All genes substantially shrunk by apeglm were shrunk by practically the same amount or more by ash. b) truth vs. estimate plots for apeglm. c) truth vs. estimate plots for ash. d) CAT plot for the three methods as well as for ML after filtering. MLE: Maximum Likelihood Estimation, apeglm: Approximate Posterior Estimation of Generalized Linear Model Coefficients, ash: Adaptive Shrinkage.

Apeglm improved mean absolute error by 49.5% among all genes, and by 65.4% among noticeably shrunk genes specifically (Table 2). Ash improved mean absolute error by 52.2% among all genes and by 65.8% among noticeably shrunk genes specifically. These improvements were greater than that seen from the standard normal simulation. Figure 2a–c show that ash successfully shrunk inflated ML estimates closer to the truth while leaving truly large effects mostly unchanged. Apeglm brought many inflated ML estimates closer to the truth as well, but not as many as ash.

Concordance at the top was better for the shrinkage estimates than for the ML estimates, regardless of the number of top genes in question. Furthermore, similar to mean absolute error, the improvements seen from the shrinkage estimates over the MLE was larger than those seen from the standard normal simulation. Ash performed better than apeglm in concordance at the top 50 and 100 genes, though performance was similar when looking at larger number of genes. Concordance at the top for the MLE was optimized when filtering out genes where less than half the samples had at least 110 counts. Though this improved CAT by quite a lot, performance was still much lower than apeglm and ash Thus, unlike in the standard normal simulation, the performance in CAT obtained by shrinkage could not be matched with filtering, even when using the true gene ranking to determine the optimal filtering rule.

In this simulation, due to the increased overdispersion, there were many effects that were overestimated or underestimated by ML, even among genes with large counts. Because of this, both ash and apeglm exhibited shrinkage for effects across the dynamic range of summed counts, as opposed to only shrinking effects with small counts (Supplementary Figure 5¹³). Together with the true vs. estimate plots, this shows that both apeglm and ash can correctly shrink falsely large effects even when the summed counts are large. Ash had larger shrinkage scores than apeglm on average, indicating that ash tended to shrink estimates more than apeglm (Supplementary Table 2¹³). All methods had coverage slightly less than nominal (95%), ranging from 92 to 94%. However, both apeglm and ash had half the average interval width compared to maximum likelihood, despite both having slightly higher coverage rates.

We also conducted simulations similar to that of the standard normal and Student’s t, but with 5 vs. 5 samples. Like the 4 vs. 4 case, both apeglm and ash had lower average estimation error and higher concordance at the top than the MLE (results not shown).

Sampling from the mouse dataset

To evaluate performance on real data, we took 100 random subsamples of 6 mice from the mouse data set and averaged various performance metrics across the random subsamples. Similar to the simulations, apeglm appeared to improve estimation accuracy and shrink erroneously large genes. Ash, on the other hand, appeared to perform worse than the MLE according to mean absolute error, concordance at the top and interval coverage (see Table 3 and Figure 3).

Figure 3. Truth vs. estimate and CAT plots for random subsamples.

a) through c) is based on one of the 100 random subsamples used in the mouse data benchmarking, and plots the ML estimates of the associated held-out set against the MLE, apeglm and ash estimates from the random subsample, respectively. Orange points represent genes substantially shrunk by ash only and green points represent genes substantially shrunk by both ash and apeglm. All genes substantially shrunk by apeglm were shrunk by practically the same amount or more by ash. d) plots concordance at the top averaged across the 100 random subsamples for each method. MLE: Maximum Likelihood Estimation, apeglm: Approximate Posterior Estimation of Generalized Linear Model Coefficients, ash: Adaptive Shrinkage.

Among genes that were shrunk by apeglm, mean absolute error was 14.2% lower on average for apeglm (Table 1). On the other hand, average MAE rose by 34.1% for ash among genes that were shrunk. Moreover, the minimum MAE obtained by ash across all 100 random subsamples was larger than the maximum MAE obtained by apeglm (results not shown). From Figure 3a–c, ash appears to be over-shrinking, and some of the genes with the largest held-out effect estimates were shrunk to zero. Though some genes also appeared to have been incorrectly or overly shrunk by apeglm, apeglm mainly was observed to shrink genes with inflated estimates and over-shrinkage was normally less severe when it occurred.

Both apeglm and MLE had universally higher concordance at the top than ash (Figure 1d). While apeglm performed slightly better than the MLE in concordance at the top 50 genes, performance was identical when looking at larger numbers of genes. We found that any filtering only decreased concordance at the top, as many top genes had low counts (i.e. the optimal filtering rule was no filtering). The most likely reason for this is that for each random subsample, we are treating the MLE of the held-out set as the truth. Thus, estimation error in the face of low-count genes would affect the held-out effect estimates and bias CAT results to some degree, even though the held-out sets have larger numbers of samples and performance metrics are averaged over many random subsamples. However, because genes with very large held-out effect estimates are more likely to have low counts, metrics that average across all genes, such as mean absolute error, would not be biased as much by estimation error.

A large amount of variability in the ML estimates was discernible for genes with low counts (Supplementary Figure 6¹³). Like in the standard normal simulation, the low-count genes were mainly the ones shrunk by apeglm and ash. As the truth vs. estimate plots suggest, ash had larger shrinkage scores than apeglm (indicating more extreme shrinkage), and with the difference in shrinkage between the two methods being larger than in the simulations (Supplementary Table 4¹³). Though apeglm intervals appeared to have smaller coverage than the ML intervals, the difference in coverage was very small, and average interval width was also 26.8% smaller for apeglm than that of maximum likelihood. Ash intervals were slightly more narrow than apeglm, with average interval width 32.5% smaller than that of maximum likelihood, but coverage was also lower.

As the mouse data set only had 24 samples, we determined that we didn’t have the sufficient sample size to evaluate our methods on estimating effect sizes of predictors, or models with more than just an intercept term. For instance, even if we wanted to look at the performance of our methods on estimating a group effect with only four samples in each group, each held-out set would only have eight samples in each group. Thus, the ML estimates of the held-out sets would have a lot of variance and could be far from the truth.

Computational performance of Apeglm

To evaluate the computational performance of our package on larger datasets, we simulated allelic counts for 5000 genes and 100 samples, and randomly divided the samples into differing numbers of groups. apeglm with our improvements had very fast running times for both ML and apeglm estimation and scaled well with the number of covariates (see Figure 4 and Figure 5).

Figure 4. Comparisons in estimation time for one gene (in milliseconds).

Figure 5. Comparisons in estimation time for all genes.

a) computational time of ML estimation (in seconds) for the apeglm and aod packages by the number of groups (covariates). b) computational time of apeglm estimation for the new and old apeglm packages by number of groups (covariates).

Estimation times per gene for ML estimation was substantially faster for apeglm than all other packages (Figure 4). The next best package, aods3, took 5 to 11 times longer than apeglm and did not scale as well with the number of groups. Furthermore, the aods3. gamlss and HRQoL packages occasionally produced errors and could not fit beta-binomial models for all the simulated genes.

For estimating all genes in the simulation via maximum likelihood, apeglm took 24 seconds for two groups and added only 1–2 seconds of computational time for every group added (Figure 5a). The next fastest package that could fit beta-binomial models for all the genes, aod, took seven times longer for two groups and grew 80 times as much for every group added. Comparisons in apeglm estimation between our improved apeglm package and the original package gave similar conclusions. Furthermore, unlike the new apeglm package, which grew roughly linearly with the number of groups in the range we assessed, the order of growth from the original package was not linear: the greater the number of groups already in the model, the greater the computational time increased for adding additional groups. At 10 groups, our improvements made apeglm 27 times faster than aod for ML estimation and 33 times faster than the old package for apeglm estimation. Our improvements also performed quite favorably when fitting beta-binomial models with two groups and two numerical controls. Elapsed time was 31 seconds for ML estimation and 43 seconds for apeglm estimation with the new apeglm package. In contrast, ML estimation took over nine minutes for aod and apeglm estimation took over seven minutes for the old apeglm package. Introducing multicollinearity into the design matrix did not substantially change computational performance for any package (results not shown).

Underlying data

Zenodo: RNA-seq Dataset from Crowley et al. 2015. http://doi.org/10.5281/zenodo.3404689¹⁶.

This project contains the following underlying data:

This file contains the Crowley et al. mouse dataset which was was obtained from http://csbio.unc.edu/gecco/data/fullGeccoRnaDump.csv.gz^15,30. We uploaded the dataset to Zenodo on the authors’ behalf with their permission, due to the fact that the original dataset is not currently hosted in a stable repository.

The dataset from this repository is available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).

Extended data

Zenodo: Supplementary Material for Zitovsky and Love 2019. http://doi.org/10.5281/zenodo.3404697¹³.

This project contains the following extended data:

Data are available under the terms of the CC-BY 4.0 license.