- Software
- Open Access

# ShrinkBayes: a versatile R-package for analysis of count-based sequencing data in complex study designs

- Mark A van de Wiel
^{1, 2}Email author, - Maarten Neerincx
^{3}, - Tineke E Buffart
^{3}, - Daoud Sie
^{4}and - Henk MW Verheul
^{3}

**15**:116

https://doi.org/10.1186/1471-2105-15-116

© van de Wiel et al.; licensee BioMed Central Ltd. 2014

**Received:**24 September 2013**Accepted:**11 April 2014**Published:**26 April 2014

## 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 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.

### 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 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.

### Conclusion

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.

## Keywords

- Differential expression
- Shrinkage
- Sequencing
- Bayesian analysis

## Background

Following the surge of count-based 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 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 [2]. 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.

## 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, 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 [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 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 [2].

### Setting

*j*= 1,…,

*n*denote independent samples,

*Y*

_{ ij }be the data for feature

*i*and sample

*j*,

*F*be the likelihood model (e.g. (zero-inflated) negative binomial) with mean

*μ*

_{ ij }and hyper-parameters γ

_{ i }and

*g*() a link-function. Here, γ

_{ i }contains distribution parameters that are not linked to covariates, e.g. zero-inflation and over-dispersion. 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, ${X}_{j}^{\alpha}$ (${X}_{j}^{\beta}$) denotes the *j*th 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 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 [2].

_{ i }=

*β*

_{ i }),

*H*

_{0 i}:

*β*

_{ i }= 0 against

*H*

_{1 i}:

*β*

_{ i }≠ 0 [1]. Here, we resolve this issue. In [1], a smooth non-parametric 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 multi-grid parameter estimation of such priors, and increased their flexibility by allowing non-equal mixture proportions for negative and positive effects. Moreover, we implemented a mixture of a spike and a smooth non-parametric component (SpNP prior). For the Results, we focus on the Spike-Gauss-Gauss (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 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 *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 half-plane. 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 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 *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 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.

### Multi-parameter inference

*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 ${\mathcal{\mathcal{M}}}_{1}=M\left({\mathit{\beta}}_{i}\right)$ when β

_{ i }is unconstrained and the null model ${\mathcal{\mathcal{M}}}_{0}=M\left(\mathbf{0}\right)$. 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 $\text{ML}({\mathbf{Y}}_{i};{\mathcal{\mathcal{M}}}_{0})$ and $\text{ML}({\mathbf{Y}}_{i};{\mathcal{\mathcal{M}}}_{1})$ are the marginal likelihoods under ${\mathcal{\mathcal{M}}}_{0}$ and ${\mathcal{\mathcal{M}}}_{1}$, 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 [6] we prefer to define threshold *t* using BFDR(*t*,**Y**_{
i
}) rather than lfdr. In any case, we need to compute $\text{ML}({\mathbf{Y}}_{i};{\mathcal{\mathcal{M}}}_{0}),\text{ML}({\mathbf{Y}}_{i};{\mathcal{\mathcal{M}}}_{1})$ and *p*_{0}.

*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 point-mass (*δ*(β= **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*${\widehat{p}}_{0}$ 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 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 [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 Savage-Dickey approximation [7]. 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.

## 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, area-under-the-curve (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 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 (${\text{B}}_{2000}^{2000}$, *p*_{0} = 0.64 and ${\text{B}}_{625}^{625}$, *p*_{0} = 0.9), because for the asymmetric cases the normalization procedure used in [1] 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 (${\text{S}}_{625}^{625}$). 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.

*n*= 5(10)] of the other methods that are based on count distributions.

**FDR results for target FDR=0.05**

Data set |
| SpGG | SpNP | NP |
---|---|---|---|---|

${\text{B}}_{2000}^{2000}$ | 5 | 0.085 | 0.078 | 0.29 |

${\text{B}}_{2000}^{2000}$ | 10 | 0.079 | 0.071 | 0.29 |

${\text{B}}_{625}^{625}$ | 5 | 0.115 | 0.115 | 0.37 |

${\text{B}}_{625}^{625}$ | 10 | 0.083 | 0.081 | 0.38 |

${\text{S}}_{625}^{625}$ | 5 | 0.111 | 0.108 | 0.38 |

${\text{S}}_{625}^{625}$ | 10 | 0.119 | 0.117 | 0.40 |

**Area-under-the-curves**

Data set |
| SpGG | SpNP | NP | Best |
---|---|---|---|---|---|

${\text{B}}_{2000}^{2000}$ | 5 | 0.897 | 0.898 | 0.85 | 0.87 |

${\text{B}}_{2000}^{2000}$ | 10 | 0.949 | 0.951 | 0.91 | 0.93 |

${\text{B}}_{625}^{625}$ | 5 | 0.874 | 0.879 | 0.82 | 0.87 |

${\text{B}}_{625}^{625}$ | 10 | 0.937 | 0.940 | 0.88 | 0.93 |

${\text{S}}_{625}^{625}$ | 5 | 0.866 | 0.871 | 0.81 | 0.85 |

${\text{S}}_{625}^{625}$ | 10 | 0.923 | 0.927 | 0.87 | 0.92 |

### 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[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 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 *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.

*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 $\hat{\text{FDR}}\le 0.1$ 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

*K*-group setting. Surprisingly, the Kruskal-Wallis test seems to be very competitive, although it also loses power for larger values of the FPR.

**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)**

Method | # Detections | True FDR |
---|---|---|

mean (sd) | Median (IQR) | |

ShrinkBayes | 37.4 (4.60) | 0.171 (0.072) |

ShrinkBayes, Old | 132.1 (15.3) | 0.509 (0.038) |

edgeR | 58.8 (12.9) | 0.258 (0.120) |

DESeq | 10.4 (3.75) | 0.191 (0.178) |

Kruskal-Wallis | 0 (0) | NA |

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

#### 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.

#### Analysis

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 ${\mathcal{\mathcal{M}}}_{0}$: 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 ${\mathcal{\mathcal{M}}}_{0}$, 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

*β*

_{ ik }-

*β*

_{ i ℓ },

*k*>

*ℓ*, which help to explain differential or non-differential 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, ${\widehat{p}}_{0}=0.92$, 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, multi-parameter 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

*K*-group comparisons only.

## 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 non-parametric component is important (see also [4]). When using an interval null-hypothesis, *H*_{0i}: |*β*_{
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[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 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 [2]. ShrinkBayes also applies to Gaussian data, like mRNA microarray data or high-throughput RNA interference screens [4]. Use is similar, as illustrated in the ShrinkBayes R-vignette, 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 **Operating system:** Platform independent (developed on Linux)**Programming language:** R**Other requirements:** R (>= 3.0.1); R-packages: INLA, from http://www.r-inla.org and snowfall, VGAM, mclust, logcondens, Iso from CRAN**Licence:** GNU GPL

## Declarations

### Acknowledgements

We thank Charlotte Soneson for providing simulated data, Andrea Riebler for discussions and Paul Eijk for his help with library preparation.

## Authors’ Affiliations

## References

- 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

## Copyright

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.