ShrinkBayes: a versatile R-package for analysis of count-based sequencing data in complex study designs
© van de Wiel et al.; licensee BioMed Central Ltd. 2014
Received: 24 September 2013
Accepted: 11 April 2014
Published: 26 April 2014
Complex designs are common in (observational) clinical studies. Sequencing data for such studies are produced more and more often, implying challenges for the analysis, such as excess of zeros, presence of random effects and multi-parameter inference. Moreover, when sample sizes are small, inference is likely to be too liberal when, in a Bayesian setting, applying a non-appropriate prior or to lack power when not carefully borrowing information across features.
We show on microRNA sequencing data from a clinical cancer study how our software ShrinkBayes tackles the aforementioned challenges. In addition, we illustrate its comparatively good performance on multi-parameter inference for groups using a data-based simulation. Finally, in the small sample size setting, we demonstrate its high power and improved FDR estimation by use of Gaussian mixture priors that include a point mass.
ShrinkBayes is a versatile software package for the analysis of count-based sequencing data, which is particularly useful for studies with small sample sizes or complex designs.
KeywordsDifferential expression Shrinkage Sequencing Bayesian analysis
Following the surge of count-based sequencing data, a plethora of software packages for differential expression analysis of such data has emerged . Many of these methods are limited in use due to restrictions on the study design, the model and inference like a) 2- or K - group comparisons only; b) no random effects; c) no explicit solution for excess of zeros and d) no multi-parameter inference. We introduced ShrinkBayes as a versatile analysis method which allows generalized linear mixed models and zero-inflation and with, due to its multi-parameter shrinkage options, good reproducibility and power characteristics . This paper illustrates the R-package ShrinkBayes on a challenging microRNA sequencing (miRseq) colon tumor-plus-metastasis study. In addition, we automated the use of mixture priors containing a spike, leading to improved FDR-based inference. Finally, we extend the class of admitted priors with mixtures of a multivariate point mass and a Gaussian product density to allow for powerful multi-parameter inference.
ShrinkBayes applies Integrated Nested Laplace Approximation, INLA, in combination with Empirical Bayes principles to provide shrunken parameter estimates and inference. In a Bayesian setting, multi-parameter shrinkage is effectuated by estimating hyper-parameters of priors. The core of ShrinkBayes is iterative estimation of priors: each prior is fit to the point-wise empirical mean of the marginal posteriors of those parameters θ i ,i = 1,…,p = # features, that correspond to the prior . Shrinkage is known to be potentially beneficial for dispersion parameters, but may be as important for parameters of interest to accomplish better inference  and for nuisance parameters to reduce their impact when unimportant .
A typical ShrinkBayes analysis consists of the following modules: a) Iterative Empirical Bayes estimation of multiple priors which need to obey the parametric forms included in INLA; b) Fitting of the full model and the null model; c) Updating one prior resulting from a) to a non-parametric or mixture prior to allow for for more flexibility and/or better inferential properties; d) Updating the posteriors of the corresponding parameters; e) Computing summary statistics including estimates of lfdr and (B)FDR. The steps are detailed in the Example section. Below we discuss novel implementations and methods with respect to .
where β i = (βi1,…,β iK ) denotes the parameter(s) for which (joint) inference is desired, while α i contains all the other regression parameters, including the intercept. In addition, () denotes the jth row of the design matrix restricted to those columns of this matrix that are relevant for α i (β i ).
ShrinkBayes inherits much of its flexibility from the INLA R-package, including its ability to deal with arbitrary designs and random effects. INLA, however, requires use of specific parametric priors. Since the prior may be crucial for inference in a multiple testing setting, we extended the class of admissible priors to non-parametric and parametric mixture priors .
where δ () is the dirac delta function, i.e. a spike. The spike is essential, because it allows the posteriors to have non-zero mass on the null-hypothesis, β i = 0, hence accommodating selection. The smooth parts of both these priors allow asymmetry between under and overexpression. All parameters are determined by maximizing the total (log-) marginal likelihood (i.e. the sum of marginal likelihoods over all features). This maximization is explicit for the parametric SpGG prior, whereas FNP is obtained by the iterative marginal procedure  with the restriction that it contains maximally one mode on both the negative and positive half-plane. The restriction helps to identify FNP together with p0. In words, given a current proposal for p0 and FNP the iterative procedure proposes a new estimate of p0 and FNP by fitting the SpNP prior to the point-wise empirical mean (over features i=1,…,p) of the current posteriors π(β i |Y i ), where the fit needs to respect the aforementioned restriction. Any reasonable starting value of p0 (we use 0.8) and FNP (we use a sufficiently vague central Gaussian, e.g. N(0,5)) can be used and convergence is checked by assessing the total (log-)marginal likelihood.
ShrinkBayes allows for other parametric priors, such as the Spike-Gauss’ (SpG) and the Spike-and-Slab’ (SpSlab). Both are mixtures of a point mass and a central Gaussian distribution, but the first has a data-adaptive variance fitted with the same direct maximization procedure as for the SpGG prior, whereas the latter has a prescribed large variance. Both alternatives are discussed in more detail in the Additional file 1.
where and are the marginal likelihoods under and , respectively. On its turn, lfdr determines BFDR(t,Y i ) = E[lfdr|lfdr < t] : the mean of all local fdrs smaller than t. Given its analogous interpretation to ordinary FDR  we prefer to define threshold t using BFDR(t,Y i ) rather than lfdr. In any case, we need to compute and p0.
hence a mixture of a multivariate point-mass (δ(β= 0)) and a Gaussian product density for the regression parameters β i = (βi 1,…,β iK ). In particular when the true p0 is large, the total (log-)marginal likelihood may contain ridges and/or multiple modalities with respect to the parameters of (5). For example, when the true p0 is large a prior (5) with small and small values of σ k may also fit rather well. To counter this, we use the constraint p0 ≥ 0.5 (which is realistic in most cases) and use a large default starting value of p0 (0.8). Moreover, iteration is stopped when the total (log-)marginal likelihood decreases by less than 0.1% to avoid walking on a ridge’.
In addition to the improved implementation of spike-priors and the multi-parameter inference, ShrinkBayes versions 2.3 and higher contain a number of novelties and changes compared to version 1.6, which corresponds to . In particular, it is faster, because convergence of the parameters of the prior(s) is assessed in terms of total marginal likelihood instead of on the separate parameters. The new version also allows to approximate marginal likelihood for a null model from the results of the full model using the Savage-Dickey approximation . This is particulary convenient for contrasts for which a null-model can not be defined without the use of constraints. Additional file 1, Section 2, contains more details and a full list of changes.
To study which of the priors performs best in terms of FDR estimation and power, we compared them on simulated data sets, including those in .
Results on simulations for various effect size distributions
The true effect size distribution, i.e. the true generating distribution of the parameter of interest, may have impact on what prior performs best. Hence, we study several effect size distributions, including a Gamma, t, Uniform and Gaussian mixture (see Additional file 1, Section 1). We compared performance of the SpGG, SpNP, SpG and SpSlab priors in terms of accuracy of FDR estimation, area-under-the-curve (AUC), number of detections and absence of detections when H0iis true for all features (p0 = 1). From the results (Additional file 1, Section 1) we conclude that SpGG and SpNP lead to accurate estimates of FDR and are very competitive in terms of power, whereas SpSlab is often too conservative; SpG generally performs well except for the (asymmetric) Gamma distribution for which it is less powerful than SpGG and SpNP. In the case p0 = 1, none of the prior returns a significant result at BFDR≤0.1, but the SpGG prior performs best in the sense that it produces the highest BFDRs.
Results on simulations in
Next, we report results of ShrinkBayes with the SpGG and SpNP priors on simulations in , which compared several methods, including ShrinkBayes (referred to as ShrinkSeq), on a variety of data sets. ShrinkBayes was used with a smooth non-parametric prior (NP), so not containing a spike. The number of features equals 12500. We focus on data sets where counts are exclusively generated from the negative binomial. Moreover, we report results on the symmetric cases (in terms of up- and down-regulation) only (, p0 = 0.64 and , p0 = 0.9), because for the asymmetric cases the normalization procedure used in  introduces artificial differential signal for the non-differential features. We do include a case with outliers which contains, on average, 5% outliers for 10% of the features (). For sample sizes we focus on n = N/2 = 5,10, because the ShrinkBayes results reported in  were relatively worse for those sample sizes.
FDR results for target FDR=0.05
Multi-parameter inference: data-based simulation
We compare our solution for multi-parameter inference to the likelihood-ratio tests that are implemented in the popular RNAseq data analysis programs edgeR and DESeq. We believe such a comparison is most meaningful and fair when the data is simulated in a relevant and realistic way, preferably avoiding distributional assumptions as much as possible. Therefore, we generated the data in three steps. First, we create a realistic null data set: we simply re-sample 3*5 observations’ from our miRseq data set, independently for each of the 2060 features. Hence, per feature 5 observations are generated from the same empirical distribution for each of the 3 groups. Next, modest filtering on the number of non-zeros is applied, because this is recommended for the use of edgeR and DESeq: at least 3 non-zeros should be present. Finally, we need a realistic effect size distribution for the features. To avoid parametric assumptions this is estimated by FNP, the smooth component of the SpNP prior (3), for the groups in the miRseq study (organs). We create 20% differential features by sampling independently from FNP for groups 2 and 3 and multiplying the respective counts by the exponentiated sampled effect sizes. This entire simulation was repeated 10 times.
Number of detections (mean and standard deviation) at target FDR = 0.1 and true FDR for the set of detections (median and IQR: interquartile range)
Note that ShrinkBayes is still liberal in the sense that it underestimates true FDR. This is probably due to the data not being generated from a specific parametric distribution. In particular, we observed that the data contains outliers for some features. Dedicated detection of such outliers can certainly reduce the number of false positives. A simple, heuristic, practical alternative is to additionally require for selection the corresponding uncorrected Kruskal-Wallis p-value to be smaller than 0.05. Then, power of a parametric approach like ShrinkBayes, which is essential in a multiple testing setting, is combined with the robustness of a nonparametric test. In this case, the median true FDR drops from 0.171 to 0.134 (target equals 0.1), while detecting 32 features on average instead of 37.4.
Example: analysis of miRseq count data
We applied ShrinkBayes to a challenging data set. The data set contains miRseq counts of 2060 miRNAs (3p- and 5p-variants) for 55 resections from primary colon tumors (P) and corresponding metastases (M) coded by the covariate PM. In addition, several other covariates are available: indiv: most individuals correspond to 2 samples (one for P, M), but some have multiple measurements for M, because the metastasis occurred at multiple locations; organ: organ where the metastasis occurred; time: binary, indicating whether resections of the primary tumor and the metastasis were at different dates; chemo: binary, indicating whether chemotherapy was applied in between the resections. In addition to other software, ShrinkBayes provides two important extra features to correctly analyze these data: it explicitly accounts for excess of zeros and allows for random effects (here indiv). Both are important for appropriate inference. In addition, we demonstrate here that joint inference for related parameters like those corresponding to organ is feasible. Note that separate inference for each organ has limited power due to the small number of samples per organ. We focus on the statistical analysis. Preprocessing is described in the Additional file 1, Section 3, which also contains annotated R-code for the entire analysis, including inferences for organ and the P-M contrast.
The analysis consists of the following steps: 1) Likelihood specification for the counts, here the zero-inflated negative binomial one; 2a) Specification of the regression model. Here, the model is the linear model with fixed effects PM, time, chemo and organ plus random effect indiv; 2b) Specification of the null-model : as , without organ; 3) Choice of parameters to shrink. Here, all fixed parameters plus the over-dispersion parameter of the negative binomial.
4) Estimation of priors for the purpose of shrinkage. Standard priors (Gaussian and inverse-Gamma) are used for all parameters, except for the inferential variable, organ, for which the multivariate mixture prior (5) is used; 5) Computation of posteriors under models and , given the prior parameters; 6) Combination of the two posteriors to one given the parameters of the mixture prior; 7) Compute local and Bayesian false discovery rates (lfdr; BFDR). The most complex steps, 4) to 7), are completely automated including setting of tested defaults, which allows users with little experience in Bayesian computing to apply ShrinkBayes. The joint mixture prior is discussed above; other technical details are given in .
For the choice of prior, we recommend to use the SpGG prior when inference on a parameter equalling zero is desired, because of its uniformly good performance in terms of FDR estimation and power. The SpNP prior is a good alternative which may be attractive in extremely small sample size settings for which the flexible shape of the non-parametric component is important (see also ). When using an interval null-hypothesis, H0i: |β i | < δ, inclusion of a spike is less relevant, so smooth (non-parametric) priors generally suffice.
Given the good performance of the SpGG prior in a univariate setting, it may be good to extend (5) to the multivariate analogue of the SpGG prior: a mixture of a multivariate point mass and a two-component Gaussian mixture product density. However, while this is conceptually feasible, it may be computationally cumbersome, because it would require combining several different fits from INLA under combinations of the components of the mixture.
Although ShrinkBayes is much more efficient that MCMC-based methods, it is computationally more demanding than frequentist counterparts like edgeR and DESeq. As an indication: the data example above (on approx 2,000 features) runs in approximately 30 minutes on 6 cpus of a Linux-cluster, whereas approximately 6 hours would be required for 100,000 features. For extremely large data sets, ShrinkBayes provides quick pre-screen functions, application of which potentially reduces computing time by a large factor.
We focused on sequencing count data for fairly complex designs. To our knowledge, extensively validated data are still not available for such studies, which hampers a thorough comparison between methods. Even when such a data set would be available, it is uncertain to what extent conclusions from one data set could be extrapolated to others, because the relative performance of a method may depend on many aspects such as the proportion of outliers and zero counts and/or the presence of multiple noise levels (e.g. within and between individuals). We emphasize that ShrinkBayes is currently the only RNAseq analysis method that can deal with the latter, by allowing random and mixed effects models, concepts that are widely accepted and used in other fields of statistical data analysis.
For simple designs, ShrinkBayes can be useful as well, in particular due to its good reproducibility, as shown for publicly available RNAseq data in . ShrinkBayes also applies to Gaussian data, like mRNA microarray data or high-throughput RNA interference screens . Use is similar, as illustrated in the ShrinkBayes R-vignette, which also contains additional examples on count data.
We illustrated the versatility of ShrinkBayes on a data set which reflects a level of complexity that is common in clinical practice. With the decrease of costs for sequencing, we are likely to encounter such complex data sets frequently in the near future and ShrinkBayes provides the means and power to analyze these.
Availability and requirements
Project name: ShrinkBayes Project home page: http://www.few.vu.nl/~mavdwiel/ShrinkBayes.html Operating system: Platform independent (developed on Linux)Programming language: ROther requirements: R (>= 3.0.1); R-packages: INLA, from http://www.r-inla.org and snowfall, VGAM, mclust, logcondens, Iso from CRANLicence: GNU GPL
We thank Charlotte Soneson for providing simulated data, Andrea Riebler for discussions and Paul Eijk for his help with library preparation.
- Soneson C, Delorenzi M: A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics. 2013, 14: 91-10.1186/1471-2105-14-91.View ArticlePubMed CentralPubMedGoogle Scholar
- Van de Wiel MA, Leday GGR, Pardo L, Rue H, van der Vaart AW, van Wieringen WN: Bayesian analysis of RNA sequencing data by estimating multiple shrinkage priors. Biostatistics. 2012, 14: 113-128.View ArticlePubMedGoogle Scholar
- Rue H, Martino S, Chopin N: Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations (with discussion). J R Stat Soc B. 2009, 71: 319-392. 10.1111/j.1467-9868.2008.00700.x.View ArticleGoogle Scholar
- Van de Wiel MA, de Menezes RX, Siebring-van Olst E, van Beusechem VW: Analysis of small-sample clinical genomics studies using multi-parameter shrinkage: application to high-throughput RNA interference screening. BMC Med Genom. 2013, 6: 1-10.1186/1755-8794-6-1.View ArticleGoogle Scholar
- Efron B: Large-scale Inference. Institute of Mathematical Statistics Monographs. 2010, Cambridge: Cambridge University PressView ArticleGoogle Scholar
- Ventrucci M, Scott EM, Cocchi D: Multiple testing on standardized mortality ratios: a Bayesian hierarchical model for FDR estimation. Biostatistics. 2011, 12: 51-67. 10.1093/biostatistics/kxq040.View ArticlePubMedGoogle Scholar
- Wetzels R, Grasman RPPP, Wagenmakers E-J: An encompassing prior generalization of the Savage-Dickey density ratio. Comp Stat Data Anal. 2010, 54: 2094-2102. 10.1016/j.csda.2010.03.016.View ArticleGoogle Scholar
- Robinson MD, McCarthy DJ, Smyth GK: edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010, 26: 139-140. 10.1093/bioinformatics/btp616.View ArticlePubMed CentralPubMedGoogle Scholar
- Anders S, Huber W: Differential expression analysis for sequence count data. Genome Biol. 2010, 11: 106-10.1186/gb-2010-11-2-106.View ArticleGoogle Scholar
- Si Y, Liu P: An optimal test with maximum average power while controlling FDR with application to RNA-seq data. Biometrics. 2013, 69: 594-605. 10.1111/biom.12036.View ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.