Abstract
Background
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 multiparameter inference. Moreover, when sample sizes are small, inference is likely to be too liberal when, in a Bayesian setting, applying a nonappropriate prior or to lack power when not carefully borrowing information across features.
Results
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 multiparameter inference for groups using a databased 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.
Conclusion
ShrinkBayes is a versatile software package for the analysis of countbased sequencing data, which is particularly useful for studies with small sample sizes or complex designs.
Keywords:
Differential expression; Shrinkage; Sequencing; Bayesian analysisBackground
Following the surge of countbased sequencing data, a plethora of software packages for differential expression analysis of such data has emerged [1]. 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 multiparameter inference. We introduced ShrinkBayes as a versatile analysis method which allows generalized linear mixed models and zeroinflation and with, due to its multiparameter shrinkage options, good reproducibility and power characteristics [2]. This paper illustrates the Rpackage ShrinkBayes on a challenging microRNA sequencing (miRseq) colon tumorplusmetastasis study. In addition, we automated the use of mixture priors containing a spike, leading to improved FDRbased 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 multiparameter inference.
Implementation
Shrinkage
ShrinkBayes applies Integrated Nested Laplace Approximation, INLA[3], in combination with Empirical Bayes principles to provide shrunken parameter estimates and inference. In a Bayesian setting, multiparameter shrinkage is effectuated by estimating hyperparameters of priors. The core of ShrinkBayes is iterative estimation of priors: each prior is fit to the pointwise empirical mean of the marginal posteriors of those parameters θ_{i},i = 1,…,p = # features, that correspond to the prior [2]. Shrinkage is known to be potentially beneficial for dispersion parameters, but may be as important for parameters of interest to accomplish better inference [2] and for nuisance parameters to reduce their impact when unimportant [4].
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 nonparametric 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 [2].
Setting
The setting is a generalized linear model. Let j = 1,…,n denote independent samples, Y_{ij }be the data for feature i and sample j, F be the likelihood model (e.g. (zeroinflated) negative binomial) with mean μ_{ij }and hyperparameters γ_{i }and g () a linkfunction. Here, γ_{i }contains distribution parameters that are not linked to covariates, e.g. zeroinflation and overdispersion. Then,
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}).
Priors
ShrinkBayes inherits much of its flexibility from the INLA Rpackage, 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 nonparametric and parametric mixture priors [2].
ShrinkBayes was praised for its power and versatility, but also criticized for its poor FDR estimation in case of a point nullhypothesis for one parameter (so β_{i }= β_{i}), H_{0 i }: β_{i }= 0 against H_{1 i }: β_{i }≠ 0 [1]. Here, we resolve this issue. In [1], a smooth nonparametric prior was used for β_{i}, which does not suit H _{0 i}. To promote more suitable priors, we simplified application of parametric mixture priors with a spike on zero by automating multigrid parameter estimation of such priors, and increased their flexibility by allowing nonequal mixture proportions for negative and positive effects. Moreover, we implemented a mixture of a spike and a smooth nonparametric component (SpNP prior). For the Results, we focus on the SpikeGaussGauss (SPGG) and SpNP priors:
where δ () is the dirac delta function, i.e. a spike. The spike is essential, because it allows the posteriors to have nonzero mass on the nullhypothesis, β_{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 F_{NP} is obtained by the iterative marginal procedure [2] with the restriction that it contains maximally one mode on both the negative and positive halfplane. The restriction helps to identify F_{NP }together with p_{0}. In words, given a current proposal for p_{0 }and F_{NP }the iterative procedure proposes a new estimate of p_{0 }and F_{NP }by fitting the SpNP prior to the pointwise 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 p_{0 }(we use 0.8) and F_{NP }(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 SpikeGauss’ (SpG) and the SpikeandSlab’ (SpSlab). Both are mixtures of a point mass and a central Gaussian distribution, but the first has a dataadaptive 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.
Additional file 1. Supplementary Material. It contains: additional simulation results, a list of changes of the software with respect to previous versions, details on the SavageDickey approximation for marginal likelihood and extensive Rcode for the miRseq example.
Format: PDF Size: 273KB Download file
This file can be viewed with: Adobe Acrobat Reader
Multiparameter inference
Multiparameter inference is desirable when the parameters represent multiple groups or covariates with a similar interpretation. In a frequentist setting, this is often done by likelihoodratio tests. Below we discuss the Bayesian counterpart. Suppose one aims at testing H_{0i }: β_{i }= 0 against H_{1 i }: β_{i }≠ 0 in a linear model M (β_{i}), which also includes response Y_{i}, covariates X and, possibly, additional parameters λ_{i}. Refer to the full model when β_{i }is unconstrained and the null model . Traditionally, comparison of two models is done by computation of the Bayes Factor (BF). However, in a multiple testing setting a good threshold for BF requires knowing p_{0}, the proportion of true null models (see [5], Ch. 5). Then, thresholding for BF is directly linked to local fdr, which simply equals
where and are the marginal likelihoods under and , respectively. On its turn, lfdr determines BFDR(t,Y_{i}) = E[lfdrlfdr < t] : the mean of all local fdrs smaller than t. Given its analogous interpretation to ordinary FDR [6] we prefer to define threshold t using BFDR(t,Y_{i}) rather than lfdr. In any case, we need to compute and p_{0}.
The marginal likelihoods and are conveniently supplied by INLA from the two separate fits of the models and . Finally, p_{0 }is determined by our iterative joint procedure[2], which determines the value of p_{0 }(along with other parameters) that maximizes the total (log)marginal likelihood with respect to prior:
hence a mixture of a multivariate pointmass (δ(β = 0)) and a Gaussian product density for the regression parameters β_{i }= (β_{i 1},…,β_{iK }). In particular when the true p_{0} 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 p_{0 }is large a prior (5) with small and small values of σ_{k }may also fit rather well. To counter this, we use the constraint p_{0 }≥ 0.5 (which is realistic in most cases) and use a large default starting value of p_{0 }(0.8). Moreover, iteration is stopped when the total (log)marginal likelihood decreases by less than 0.1% to avoid walking on a ridge’.
Additional changes
In addition to the improved implementation of spikepriors and the multiparameter inference, ShrinkBayes versions 2.3 and higher contain a number of novelties and changes compared to version 1.6, which corresponds to [2]. 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 SavageDickey approximation [7]. This is particulary convenient for contrasts for which a nullmodel can not be defined without the use of constraints. Additional file 1, Section 2, contains more details and a full list of changes.
Results
Priors
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 [1].
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, areaunderthecurve (AUC), number of detections and absence of detections when H_{0i }is true for all features (p_{0 }= 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 p_{0 }= 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 [1]
Next, we report results of ShrinkBayes with the SpGG and SpNP priors on simulations in [1], which compared several methods, including ShrinkBayes (referred to as ShrinkSeq), on a variety of data sets. ShrinkBayes was used with a smooth nonparametric 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 downregulation) only (, p_{0 }= 0.64 and , p_{0 }= 0.9), because for the asymmetric cases the normalization procedure used in [1] introduces artificial differential signal for the nondifferential 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 [1] were relatively worse for those sample sizes.
Table 1 contains the results on FDR estimation. Note that the target FDR equals 0.05 here, in order to be consistent with [1]. We observe that ShrinkBayes with SpGG or SpNP is still liberal, but the results are much better than those for the NP prior. In fact, when comparing the results of Table 1 with those in Figure four of [1], we observe that ShrinkBayes has improved from the worse to at least average in terms of FDR estimation. In particular, for the data sets with outliers it outperforms 5/6 (4/6) [ n = 5(10)] of the other methods that are based on count distributions.
Table 1. FDR results for target FDR=0.05
Table 2 contains the results on AUC. Again, we observe a uniform improvement when using ShrinkBayes with SpGG or SpNP instead of NP. Strikingly, ShrinkBayes with both SpGG and SpNP generally outperforms all the other methods reported in [1] when it comes to AUC.
Table 2. Areaunderthecurves
Multiparameter inference: databased simulation
We compare our solution for multiparameter inference to the likelihoodratio tests that are implemented in the popular RNAseq data analysis programs edgeR[8] and DESeq[9]. 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 resample 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 nonzeros is applied, because this is recommended for the use of edgeR and DESeq: at least 3 nonzeros should be present. Finally, we need a realistic effect size distribution for the features. To avoid parametric assumptions this is estimated by F_{NP}, 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 F_{NP} for groups 2 and 3 and multiplying the respective counts by the exponentiated sampled effect sizes. This entire simulation was repeated 10 times.
We analyzed the simulated data sets using ShrinkBayes, edgeR, DESeq and a simple nonparametric KruskalWallis test. In addition, the old version of ShrinkBayes was applied with a smooth nonparametric prior and an a posteriori multiple comparison of the 3 groups, as suggested in [2]. Figure 1 shows the ROC curves, as averaged over the 10 repeats, for False Positive Rate (FPR) smaller than 0.05. We focus on this FPR range, because when using as a selection criterion, all 5 methods produce sets of significant features with FPR≤0.05. ShrinkBayes seems somewhat superior to edgeR across the entire range, while it is competitive with DESeq. Possibly due to the smoothness of the prior ShrinkBayes,Old performs a little bit better in terms of ranking than ShrinkBayes for very small FPR, but becomes inferior for larger values. The latter may be caused by loss of power when using a multiple comparison approach in a Kgroup setting. Surprisingly, the KruskalWallis test seems to be very competitive, although it also loses power for larger values of the FPR.
Figure 1. ROC curves for multiparameter inference: mean False Positive Rate (FPR; xaxis) vs mean True Positive Rate (TPR; yaxis), as averaged over 10 repeats of the databased simulation, which consists of 3 groups with 5 counts for ≈2000 features.
ROC curves, however, only allow comparison of the rankings. In practice, the actual selection is most important. Table 3 shows the results summarized over the 10 repeats when using as selection criterion. Note that for all pvaluebased methods we use the BenjaminiHochberg FDR correction, which is appropriate here given the independent sampling per feature in our simulated data set. BFDR is used as an estimate of FDR in the ShrinkBayes setting. True FDR is evaluated on the selected sets by simply dividing the number of false positives by the total number of positives. Here, the differences are much clearer: the KruskalWallis test is useless in this setting, because it does not select anything. ShrinkBayes,Old selects too much at a too high true FDR, probably due to the smooth prior, as discussed before. DESeq and ShrinkBayes produce better true FDRs (with the DESeq ones more variable), but, on average, ShrinkBayes detects almost four times as many features. edgeR selects more, but is both more liberal and more variable. In fact, as can be inferred from the ROC curves, ShrinkBayes would achieve a smaller true FDR with the same number of detections as edgeR.
Table 3. 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 KruskalWallis pvalue 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
Data
We applied ShrinkBayes to a challenging data set. The data set contains miRseq counts of 2060 miRNAs (3p and 5pvariants) 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 Rcode for the entire analysis, including inferences for organ and the PM contrast.
Analysis
The analysis consists of the following steps: 1) Likelihood specification for the counts, here the zeroinflated 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 nullmodel : as , without organ; 3) Choice of parameters to shrink. Here, all fixed parameters plus the overdispersion parameter of the negative binomial.
4) Estimation of priors for the purpose of shrinkage. Standard priors (Gaussian and inverseGamma) 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 [2].
Discoveries
At BFDR = 0.10, we discovered 43 miRs for which organ is associated to expression in the metastasis. Figure 2 shows two posteriors of contrasts β_{ik } β_{i ℓ},k > ℓ, which help to explain differential or nondifferential miR expression. For example, for the significant differential miR, which corresponds to the left display of Figure 2, differences are largest between organs 0 and 3 on one side and organs 1 and 2 on the other. To accommodate users, ShrinkBayes contains functions to easily produce such posterior plots and also summary tables. Importantly, the estimate of p_{0 }in (5) is large, , which implies strong shrinkage of organ effects towards zero, rendering more degrees of freedom’ and hence more power for other inferences. This is another strong aspect of ShrinkBayes: in studies with relatively few samples, multiparameter shrinkage helps to increase power for a particular parameter of interest [4]. The idea of jointly shrinking multiple parameters was recently also adopted in [10], although their approach currently applies to Kgroup comparisons only.
Figure 2. Posterior densities and joint nullprobability (pi0All) of 6 contrastsβ_{ik } β_{i ℓ}, k>ℓ, representing logefold expression differences (xaxis) between 4 organs, for a significant miR (left) and nonsignificant miR (right).
Discussion
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 nonparametric component is important (see also [4]). When using an interval nullhypothesis, H_{0i }: β_{i} < δ, inclusion of a spike is less relevant, so smooth (nonparametric) 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 twocomponent 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 MCMCbased methods, it is computationally more demanding than frequentist counterparts like edgeR[8] and DESeq[9]. As an indication: the data example above (on approx 2,000 features) runs in approximately 30 minutes on 6 cpus of a Linuxcluster, whereas approximately 6 hours would be required for 100,000 features. For extremely large data sets, ShrinkBayes provides quick prescreen 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 [2]. ShrinkBayes also applies to Gaussian data, like mRNA microarray data or highthroughput RNA interference screens [4]. Use is similar, as illustrated in the ShrinkBayes Rvignette, which also contains additional examples on count data.
Conclusion
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 webciteOperating system: Platform independent (developed on Linux)Programming language: ROther requirements: R (>= 3.0.1); Rpackages: INLA, from http://www.rinla.org webcite and snowfall, VGAM, mclust, logcondens, Iso from CRANLicence: GNU GPL
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
MAvdW developed the methodology and the software. MAvdW and MN wrote the manuscript. HMWV and TEB conceived of the miRseq study, which was carried out by MN. DS facilitated the sequencing and performed the initial analysis of the data. All authors read, discussed and approved the manuscript.
Acknowledgements
We thank Charlotte Soneson for providing simulated data, Andrea Riebler for discussions and Paul Eijk for his help with library preparation.
References

Soneson C, Delorenzi M: A comparison of methods for differential expression analysis of RNAseq data.
BMC Bioinformatics 2013, 14:91. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

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:113128. PubMed Abstract  Publisher Full Text

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:319392. Publisher Full Text

Van de Wiel MA, de Menezes RX, Siebringvan Olst E, van Beusechem VW: Analysis of smallsample clinical genomics studies using multiparameter shrinkage: application to highthroughput RNA interference screening.
BMC Med Genom 2013, 6:1. BioMed Central Full Text

Efron B: Largescale Inference. Institute of Mathematical Statistics Monographs. Cambridge: Cambridge University Press; 2010.

Ventrucci M, Scott EM, Cocchi D: Multiple testing on standardized mortality ratios: a Bayesian hierarchical model for FDR estimation.
Biostatistics 2011, 12:5167. PubMed Abstract  Publisher Full Text

Wetzels R, Grasman RPPP, Wagenmakers EJ: An encompassing prior generalization of the SavageDickey density ratio.
Comp Stat Data Anal 2010, 54:20942102. Publisher Full Text

Robinson MD, McCarthy DJ, Smyth GK: edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.
Bioinformatics 2010, 26:139140. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Anders S, Huber W: Differential expression analysis for sequence count data.
Genome Biol 2010, 11:106. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Si Y, Liu P: An optimal test with maximum average power while controlling FDR with application to RNAseq data.
Biometrics 2013, 69:594605. PubMed Abstract  Publisher Full Text