Differential Expression Analysis for RNA-Seq: An Overview of Statistical Methods and Computational Software

Abstract

Deep sequencing has recently emerged as a powerful alternative to microarrays for the high-throughput profiling of gene expression. In order to account for the discrete nature of RNA sequencing data, new statistical methods and computational tools have been developed for the analysis of differential expression to identify genes that are relevant to a disease such as cancer. In this paper, it is thus timely to provide an overview of these analysis methods and tools. For readers with statistical background, we also review the parameter estimation algorithms and hypothesis testing strategies used in these methods.

Introduction

In the past decade, deep sequencing has emerged as a powerful alternative to microarrays for the high-throughput profiling of gene expression. Comparing with microarrays, RNA sequencing (RNA-seq) possesses a number of technological advantages such as a wider dynamic range and the freedom from predesigned probes.^1–3 It also comes with a unique data feature as discrete sequencing reads. In order to account for this unique data feature, statistical methodologies and computational algorithms have been developed based on various data distributional assumptions such as Poisson, negative binomial, beta binomial, (full or empirical) Bayesian, and nonparametric.^4–16,

For researchers who are new to the analysis of RNA-seq data, in this paper we provide an introductory overview of the methods and software available for the differential expression analysis (DEA) of RNA-seq data when the analysis goal is to identify genes that are relevant to a disease such as cancer.^1,17,18 In addition, for those who are interested in the statistical aspects of these methods, we also provide an overview of their parameter estimation algorithms and hypothesis testing strategies. The overview of these statistical aspects in our paper provides a unique contribution to the review literature on RNA-seq DEA methods.^3,18–23 For readers who are interested in a performance comparison of RNA-seq DEA methods, they can refer to a large body of such papers in the literature.^20–23

The rest of the paper is organized as follows. In the Notation and Normalization Methods section, we introduce the unified notations used for the methods reviewed in our paper and touch on the normalization methods typically used to preprocess RNA-seq data before DEA. In the Statistical Modeling of RNA-seq Data section, we review the statistical modeling RNA-seq DEA categorized by the distributional assumptions such as Poisson,^4–6 negative binomial,^7–10 beta binomial,^11,12 Bayesian,^13,14 and nonparametric.^15,16 All reviewed methods directly work with gene-level count data for DEA and have available R packages. For interested readers with advanced statistical knowledge, the parameter estimation algorithm for each method is presented separately in a text box, and the typical statistical testing frameworks that have been proposed for RNA-seq DEA are reviewed in the Statistical Testing section. Finally, computational tools implemented for the reviewed methods are summarized in Table 2. We note that the methods reviewed in this paper are not an exhaustive collection of available methods in the literature. Rather, we reviewed a list of most commonly used categories of modeling assumptions and included a few representative methods for each category, to help researchers who are new to the field orientated and started in the still evolving literature on this topic.

Table 1.

List of sequencing depth normalization methods and reference papers.

METHODS	RELEVANT REFERENCES
RPKM	Mortazavi et al. 26
Upper-quartile, Median	Bullard et al. 18
TMM	Robinson et al. 7
DESeq	Anders and Huber 8
Quantile	Bolstad et al. 32

Table 2.

RNA-seq count data DEA statistical methods and software.

MODEL	SOFTWARE	REFERENCES	DATA TYPE	TESTING STRATEGY	NOTES	LIMITATIONS	DATA USED
Poisson	DEGseq	Wang et al. 4	RNA-seq data	Fisher's exact test Likelihood ratio test	Support raw read counts or normalized gene expression values, identify DE of exons or transcripts	Ignore biological variation	Marioni RNA-seq data
	Myrna	Langmead et al. 5	RNA-seq data	Likelihood ratio test Parallelized permutation test	Handle dataset with over 1 billion rows, computationally efficient	Ignore biological variation, signal loss due to junction or repetitive reads, inconvenient cloud data transfer	HapMap expression data
	PoissonSeq	Li et al. 6	RNA-seq data Tag-seq data	Score test	Accommodate multiple covariate types, computationally efficient	Transformation power depends only on gene expression, libraries are totally exchangeable	Marioni RNA-seq data Tag-seq data
Negative binomial	edgeR	Robinson et al. 7	SAGE data RNA-seq data	Exact test Likelihood ratio test	Separate biological from technical variations	Limited to pairwise comparison	SAGE data Fly RNA-seq data
	DESeq	Anders and Huber 8	Tag-seq data RNA-seq data ChIP-seq data	Exact test Likelihood ratio test	Extend edgeR by allowing more general, data-driven relationship of mean and variance	Limited to pairwise comparison	Neural stem cell Tag-seq data Yeast RNA-seq data HapMap ChIP-seq data Fly RNA-seq data
	DESeq2	Love et al. 9	Tag-seq data RNA-seq data ChIP-seq data	Wald test Likelihood ratio test	Improve upon DESeq for better gene ranking, allow hypothesis tests above and below threshold	Limited to pairwise comparison	Fly RNA-seq data Mouse straiturn RNA-seq data
	NBPSeq	Di et al. 10	RNA-seq data	Adapted exact test	Introduce an additional parameter to allow the dispersion to depend on the mean	Assume all library sizes are equal	Arabidopsis RNA-seq data
Beta binomial	BBSeq	Zhou et al. 12	RNA-seq data	Wald test Likelihood ratio test	Handle outlier detection automatically	Sensitive to outliers of shrinkage or penalization methods	HapMap RNA-seq data
Bayesian and
Empirical Bayesian	ShrinkSeq	Van de Wiel et al. 13	RNA-seq data CAGE data	Evaluating posterior probability for inference	Provide joint shrink multiple parameters, allow for random effects, address multiplicity problems	Computationally intensive but allow parallelization	HapMap RNA-seq data CAGE data
	baySeq	Hardcastle and Kelly 14	Small RNAs data	Evaluating posterior probability for inference	Involve multiple comparison, accommodate different sample size	Computationally intensive but allow parallelization	Trans-acting small RNAs
Nonparametric	SAMseq	Li and Tibshirani 15	RNA-seq data Tag-seq data miRNA-seq data	Wilcoxon test	Robust to outliers, remove Experimental effect, simplify test for feature effect, accommodate quantitative, survival and multiple group comparison	Overestimate FDR in some cases, relative low power for data with small sample size	Marioni RNA-seq data t'Hoen Tag-seq data Witten miRNA-seq data
	NOIseq	Tarazona et al. 16	RNA-seq data	Wilcoxon test	Robust and maintain a high true-positive rate	Not easy to identify true differential expression at a low count range, limited to pair-wise comparison	Marioni RNA-seq data

Notation and Normalization Methods

Notation

RNA-seq data for G genes and N samples can be described by a G x N matrix Y. Each entry y_gi (g = 1, …, G, i = 1, N) represents the count of sequencing reads for gene g in sample i. For a given g and i, y_gi is a nonnegative integer representing the number of reads mapped to gene g in sample i. For succinctness, we also use notations “·” for summations, eg, y_g = ${\sum }_{i=1}^N{y}_{gi}$ and ${\sum }_{g=1}^G{y}_{gi}.$ .

We use X to represent an N x P design matrix, where P is the number of covariates. For instance, x_ip can be an indicator variable of disease status, taking a value of 0 for a normal sample and a value of 1 for a tumor sample. When comparing K groups of samples, C_k represents the collection of indices of the samples in group k (k = 1, …, K), that is, C_k = {i: x_i = k}. Each sample can only belong to one group.

Normalization Methods

Similar to microarray data, RNA-seq data are also prone to nonbiological effects due to the experimental process. Consequently, these effects need to be adjusted before any further data analysis. ²⁴ One major source of nonbiological effects is sequencing depth, which can be adjusted by rescaling the sequencing counts with factors that mimic sequencing depth. ²⁵ Reads per kilobase per million reads (RPKM) is a simple adjustment that considers gene counts standardized by the gene length and the total number of reads in each library as expression values.^17,26 More sophisticated adjustment factors, including trimmed mean of M-values (TMM), ²⁷ DESeq size factor, ²⁸ and quantile-based normalizations such as upper quartile normalization, ¹⁸ are given in Table 1. Other sources of nonbiological effects for RNA-seq include gene length and GC-content,^21,29 whose effects are typically assumed to be consistent across samples for a given gene and hence cancel out in the analysis of differential expression. Interested readers can look up available normalization methods adjusting for gene length and GC-content in the publications such as Risso et al. ²⁹ , Benjamini and Speed, ³⁰ and Hansen et al. ³¹

Statistical Modeling of RNA-Seq Data

Poisson

Overview

Models for read counts originated from the idea that each read is sampled independently from a pool of reads and hence the number of reads for a given gene follows a binomial distribution, which can be approximated by a Poisson distribution. Based on the Poisson model assumption for repeated sequencings of a sample, Marioni et al. ¹⁷ proposed to use a log-linear model to model the mean difference between two samples and adopted the classical likelihood ratio test for calculating the P-values. Based on the same Poisson assumption, Bullard et al. ¹⁸ proposed to use two other test statistics, exact test statistics and score test statistics, in the generalized linear model (GLM) framework. Li et al. ⁶ proposed a method called PoissonSeq, which adapts a two-step procedure for fitting a Poisson model. The method first estimates sequencing depths using a Poisson goodness-of-fit statistic and then calculates a score statistic based on a log-linear model. In addition, Wang et al. ⁴ developed an R package, DEGseq, to identify differentially expressed (DE) genes with an MA-plot-based approach. Langmead et al. ⁵ incorporated cloud computing in their method called Myrna.

Modeling

In a Poisson model, one assumes that Y_gi,, the number of reads mapped to gene g in sample i, follows a Poisson distribution, y_gi ∼ Poisson(μ_gii). μ _gi is the rate parameter for gene g in sample i, which equals both the mean and the variance of the read counts. The probability mass function is: $$f(y_{gi}|{\mu }_{gi})=P(Y_{gi}|{\mu }_{gi})=\frac{{\mu }_{gi}^{y_{gi}}\exp (-{\mu }_{gi})}{y_{gi}!}$$ (3.1.1) and E(Y_gi) = μ _gi and Var(Y_gi) = μ _gi The association of μ _gi with the same sample group can be described by a log-linear model as follows: $$\log \left({\mu }_{gi}\right)=\log d_i+\log \beta +\sum _{k=1}^K{\gamma }_{gk}I\left(i\in C_k\right),$$ (3.1.2) where d_i represents the sequencing depth of sample i and ${\sum }_{i=1}^N{d}_i$ = 1 is assumed for generality. Let βg be the expression level of gene g and γ _g be the association of gene g with the covariate. For hypothesis testing, γ_g1 = … = γ _gK = 0 indicates that the expression of gene g is not associated with the sample group. In the case of two sample group comparison, if γ _g = 0, then gene g is not DE between the two sample groups.

Algorithm Overview 1: Li et al.'s ⁶ PoissonSeq

Li and others proposed PoissonSeq that assumes the hypotheses as follows. Under the null hypothesis where genes and covariates are not relevant, $$\log {\mu }_{gi}=\log d_i+\log {\beta }_g,$$ (3.1.a) where d_i is the sequencing depth in sample i and β _g is the expression of gene g. The model fit from Equation (3.1.a) is denoted as ${\text{N}}_{gi}^{\left(0\right)}$ in later equations: ${\text{N}}_{gi}^{\left(0\right)}=\exp \left(\log \left({\widehat{d}}_i\right)+\log \left({\widehat{\beta }}_g\right)\right)$ .

Under the alternative hypothesis where genes and covariates, ${\text{x}}_i^{*}$ are relevant, $$\log {\mu }_{gi}=\log d_i+\log {\beta }_g+{\gamma }_g{\text{x}}_i^{*}$$ (3.1.b) where ${\text{x}}_i^{*}$ would be $I_{\left(i\in c_k\right)}$ when comparing two or multiple sample groups. The authors suggested using the maximum likelihood to estimate ${\widehat{\beta }}_g,$ as a result ${\widehat{\beta }}_g=y_g$ . However, instead of using the maximum likelihood estimate of the sequencing depth in sample i, the authors sought for a set of genes, denoted by S, that are not DE to estimate sequencing depth in sample i: $${\widehat{d}}_i=\frac{{\sum }_{g\in S}{y}_{gi}}{{\sum }_{g\in S}{y}_g}.$$ (3.1.c)

They then estimated which genes belong to S by a Poisson goodness-of-fit statistic, ie, $$GO{F}_g=\sum _{i=1}^N\frac{{\left(y_{gi}-{\widehat{d}}_i{y}_{g.}\right)}^2}{{\widehat{d}}_i{y}_{g.}}$$ (3.1.d)

S is set to be the genes whose GOF _g values are in the(∊, 1 - ∊) quantile of all GOF _g values. Li and others used ∊ = 0.25 in their study. ⁶

The objective is to test H₀: $${\gamma }_{g1}=\mathrm{\dots }\cdot ={\gamma }_{gk}=0,$$ and score statistics were proposed to perform the testing. For a two-group or multiple-group covariate, the score statistic for gene g is $$\sum _{k=1}^K\frac{{\left[{\sum }_{i\in C_k}\left(y_{gi}-{\text{N}}_{gi}^{\left(0\right)}\right)\right]}^2}{{\sum }_{i\in C_k}{\text{N}}_{gi}^{\left(0\right)}}\sim X^2\left(K-1\right).$$ (3.1.e)

With accumulating empirical data (especially “with the data available for groups of multiple biological samples), researchers began to observe that in a group, the between-sample variation of sequencing reads for a gene often exceeds the mean.^17,23,33 This excessive variation that cannot be explained by the Poisson model is called overdispersion. Extensions of the classic Poisson model have been proposed in order to accommodate such overdispersion, including the two-stage Poisson models ³⁴ and the generalized Poisson model. ³⁵

Negative Binomial

Overview

A class of models based on the negative binomial distribution assumption has been developed in order to accommodate the overdispersion among biological replicate data.^8,9,33,36 Robinson and Smyth ³³ used the conditional maximum likelihood (CML) to estimate the dispersion parameter-a measure of the excessive variance that a Poisson model does not incorporate-when assuming a common dispersion parameter across genes. They compared the CML method with alternative estimation methods based on pseudolikelihood, quasi-likelihood, and conditional inference.^37–39 In a follow-up paper, ³⁶ they also extended the model to allow for gene-specific dispersion parameters and proposed to estimate the dispersion parameters by maximizing a weighted conditional likelihood with empirical Bayesian approximation. Details of their method, edgeR, can be found in Robinson and Smyth.^33,36 edgeRun is based on the same model as edgeR but it uses an unconditional exact test to achieve more power while paying the price of computational time. ⁴⁰ Anders and Huber ⁸ proposed a method called DESeq also under the negative binomial assumption. They advocated the use of a robust estimate of normalization factors for the estimation of dispersion parameter and a local regression to obtain smooth function for each group on the graphs of expected proportions vs sample variances. DESeq2 was developed in the study by Love et al. ⁹ as a successor of DESeq. It employs a number of new modeling features, such as the use of a shrunken fold change and a shrunken dispersion estimation method, to further improve the model performance. Di and others ¹⁰ proposed a method, NBPSeq, using a negative binomial power distribution instead of a regular negative binomial distribution. They hypothesized that $E\left(Y_{gi}\right)={\mu }_{gi},Var\left(Y_{gi}\right)={\mu }_{gi}\left(1+\varphi {\mu }_{gi}^{\alpha -1}\right),$ and ϕ is common across genes while a helps to accommodate the overdispersion. ϕ and α are estimated by maximizing conditional log-likelihood, ⁴¹ conditional on the total gene counts for each gene g. An exact test modified for negative binomial power distribution is used for hypothesis testing. More details can be found in the study by Di et al. ¹⁰

Modeling

The model setup for negative binomial is to assume y_gi ∼ negative binomial $({\mu }_{gi},{\varphi }_g)$ . The dispersion parameter, ϕ_g, accounts for the sample-to-sample variability, which is usually assumed to be common across samples. There are various estimation methods for this model assumption. More specifically, the negative binomial probability mass function is written as $$\begin{array}{c}f\left(y_{gi}|{\mu }_{gi},{\varphi }_g\right)=P\left(Y_{gi}|{\mu }_{gi},{\varphi }_g\right)\\ =\frac{\Gamma \left(y_{gi}+{\varphi }_h^{-1}\right)}{\Gamma \left({\varphi }_g^{-1}\right)\Gamma \left(y_{gi}+1\right)}{\left(\frac{1}{1+{\mu }_{gi}{\varphi }_g}\right)}^{{\varphi }_g^{-1}}{\left(\frac{{\mu }_{gi}}{{\varphi }_g^{-1}+{\mu }_{gi}}\right)}^{y_{gi}},\end{array}$$ (3.2.1) where $E\left(V_{gi}\right)={\mu }_{gi}$ and Var $\left(V_{gi}\right)={\mu }_{gi}+{\varphi }_g{\mu }_{gi}^2$ . Hypothesis testing is set up as H₀: no difference either between the expected normalized expression of gene g in groups or between the proportion of reads that are gene g in groups.

Algorithm Overview 2: Overdispersion

Negative binomial can be derived as a gamma-Poisson mixture model (subscripts g's and i's are omitted for brevity), under the assumption that technical replicates follow a Poisson distribution, and biological replicates follow a gamma distribution, with the latter accommodating the overdispersion observed in empirical data. $$\begin{array}{c}y\sim poisson\left(\mu \right),\mu \sim gamma\left(\alpha ,\beta \right)\\ P\left(y|\mu \right)=\frac{{\mu }^y\exp \left(-\mu \right)}{y!}\\ f\left(\mu \right)={\left(\Gamma \left(\alpha \right){\beta }^{\alpha }\right)}^{-1}\left({\mu }^{\alpha -1}\exp \left(-\mu /\beta \right)\right)\end{array}$$

Then, $$\begin{array}{c}P\left(y\right)={\int }_0^{\infty }P\left(y_{gi}|\mu \right)f\left(\mu \right)d\mu \\ ={\left(y!\Gamma \left(\alpha \right){\beta }^{\alpha }\right)}^{-1}{\int }_0^{\infty }{\mu }^{\left(y-a\right)-1}\exp \left(-\mu \left(1+1/\beta \right)\right)d\mu \\ =\frac{\Gamma \left(y+\alpha \right){\beta }^y}{y!\Gamma \left(\alpha \right){\left(1+\beta \right)}^{y+\alpha }}\end{array}$$

One substitutes back ${\mu }_{gi},y_{gi},\alpha ={\varphi }_g^{-1}$ , and $\beta ={\mu }_{gi}{\varphi }_g$ , a gamma–Poisson mixture can be viewed as a negative binomial, see Equation (3.2.1).

Algorithm Overview 3: Robinson and Smyth's^33,36 edgeR

In edgeR, ${\mu }_{gi}=m_i{\lambda }_{gk\left(i\right)}$ where m_i is the ith library size and ${\lambda }_{gk\left(i\right)}={\sum }_{i=1}^{C_k}{\lambda }_{gi}$ represents the proportion of the total reads that is gene g in group k and ${\lambda }_{gi}$ is the proportion of the total reads that is gene g in sample i.

Under the assumption of gene-wise (or tag-wise in the original paper) dispersion, ϕ _g is estimated by maximizing a weighted conditional log-likelihood, WL(ϕ _g ): $$WL\left({\varphi }_g\right)=l_g\left({\varphi }_g\right)+\alpha l_c\left({\varphi }_g\right)$$ (3.2.a) where α is the weight given to the common likelihood, l_C; the maximum estimator of WL(ϕ _g ) is denoted by ${\widehat{\varphi }}_g^{WL}$ . An α has to be chosen such that ${\widehat{\varphi }}_g^{WL}$ coincides with an empirical Bayesian solution, ${\widehat{\varphi }}_g^B$ , the Bayesian posterior mean estimator of ϕ _g where ${\widehat{\varphi }}_g|\varphi \sim N\left({\varphi }_g,{\tau }_g^2\right)$ and ${\widehat{\varphi }}_g\sim N\left({\varphi }_0,{\tau }_0^2\right)$ for g = 1, …, G. The approximation method is selected as a direct estimate of ϕ _g is difficult because of the lack of a conjugate prior for ϕ in negative binomial model. Details are given in the study by Robinson and Smyth. ³³

In the study by Robinson and Smyth, ³⁶ the overdispersion parameter is assumed to be common across all genes (ie, ϕ _g = ϕ). To estimate the shared dispersion parameter with and without equal library size, the authors proposed to use the CML and quantile-adjusted CML (qCML) as follows.

In a special case where m_i = m for i ∈ C_k where C_k = {i: k(i) = k}, ∼ negative binomial $\left({\mu }_{gi}=m{\lambda }_{gk},\varphi \right)$ in group k and Y_gi's evidently become identically distributed, and the maximum likelihood estimator (MLE) ${\widehat{\lambda }}_{gk\left(i\right)}$ becomes $\frac{{\sum }_{i\in C_k}{y}_{gi}}{{\sum }_{i\in C_k}{m}_i}$ in group k. CML function for dispersion ϕ given $z_k={\sum }_{i=1}^{n_k}{y}_{ki}$ was proposed. The function is as follows: $$\begin{array}{c}{l}_C\left(\varphi \right)=\sum _{g=1}^G{1}_g\left(\varphi \right)=\sum _{g=1}^g\sum _{k=1}^k[\sum _{j=1}^{n_k}\log \Gamma \left(y_{ki}+{\varphi }^{-1}\right)\\ +\log \Gamma \left(n_k{\varphi }^{-1}\right)\\ -\log \Gamma \left(Z_k+n_k{\varphi }^{-1}\right)-n_k\log \Gamma \left({\varphi }^{-1}\right)\end{array}$$ (3.2.b)

In the case of different m_i in group k, the MLE of ${\lambda }_{gk\left(i\right)}$ depends on ϕ (ie, maximum likelihood estimation of the two parameters proceeds jointly). As a result, an approximate approach called qCML was proposed to equate the library sizes. The quantile-adjusted pseudodata supposedly allows one to use a common likelihood $l_c\left(\varphi \right)$ to estimate an accurate estimate of ϕ. Specifically, let m^* = ${\left({\prod }_{i=1}^N{m}_i\right)}^{\frac{1}{N}}$ , where m^* is the geometric mean of the library sizes. Then, the observed data could be adjusted as if they were all sampled as identically distributed negative binomial $\left(m*\lambda ,\varphi \right)$ .

Hypothesis testing is set up as $H_0:{\lambda }_{g1}={\lambda }_{g2}$ in other words, no difference in proportion of gene g in samples between group 1 and group 2.

Algorithm Overview 4: Anders and Huber's ⁸ DESeq

The read count y_gi is modeled by a GLM of negative binomial distribution with a log link: $$\log \left({\lambda }_{gi}\right)=\sum _{p=1}^p{x}_{ip}{\beta }_{gp}$$ (3.2.c)

The mean μ _gi is the proportion of reads for gene g in sample i, ${\lambda }_{gk\left(i\right)}$ , scaled by a normalization factor, m_i. The variance ${\sigma }_{gi}^2$ is ${\mu }_{gi}+{\text{m}}_i^2{V}_{gk\left(i\right)},$ where $V_{gk(i)}$ is assumed to be a per gene raw variance, a smoothing function of λ _g and k. The use of the smoothing function can help stabilize the variance estimates especially when the number of samples is small. For the estimation of the normalization factor (which is referred to as the size factor by Anders and Huber), m, for each sample, the authors noted that highly DE genes are more likely to be influential on total count and so the median of the ratios of counts should be used for more robustness: $${\widehat{m}}_i=median\frac{y_{gi}}{{\left(\prod _{v-1}^N{y}_{gv}\right)}^{1/N}}$$ (3.2.d)

Since ${\lambda }_{gk\left(i\right)}$ is proportional to the expected value of the unknown proportion from gene g in group k, it is estimated by the average of counts from all samples in group k with a common scale.

$${\widehat{\lambda }}_{gk\left(i\right)}=\frac{1}{M_k}\sum _{i:k\left(i\right)=k}\frac{y_{gi}}{{\widehat{m}}_i},$$ (3.2.e) where M_k is the total number of replicates for group k. The sample variances with the common scale are calculated as: $$w_{gk}=\frac{1}{M_k-1}{\sum _{i:k\left(i\right)=k}\left(\frac{y_{gi}}{{\widehat{m}}_i}-{\lambda }_{gk\left(i\right)}\right)}^2$$ (3.2.f) $$z_{gk}=\frac{{\widehat{\lambda }}_{gk\left(i\right)}}{M_k}\sum _{i:k\left(i\right)}\frac{1}{{\widehat{m}}_i}$$ (3.2.g)

In the case of a sufficiently large number of M_k, one can see $w_{gk}-z_{gk}$ as the unbiased estimator of the raw variance v_gk. In the case of a small number of M_k, local regression for a smooth function w_k(λ) on the graph of $\left({\widehat{\lambda }}_{gk\left(i\right)},w_{gk}\right)$ was suggested so that $w_k\left({\widehat{\lambda }}_{gk\left(i\right)}\right)-z_{gk}$ would be the estimate for the raw variance. More details are in the study by Anders and Huber. ⁸

Algorithm Overview 5: Love et al.'s ⁹ DESeq2

DESeq2 allows the normalization factors to be gene specific (m_gi), rather than being fixed across genes (m_i). The estimation of m_gi is implemented in their new R packages. ⁹

When modeling dispersion parameters, a large variation in estimates usually arises because of small sample sizes. DESeq2 proposed to pool genes with similar average expression together for the estimation of dispersions. To do this, one first separately estimates dispersion with maximum likelihood. Then, one identifies a location parameter for the distribution of the estimates by fitting a smooth curve dependent on average normalized expressions, before finally shrinking gene-specific dispersions to the fitted curve using an empirical Bayesian approach. The authors stated that this procedure is more superior than DESeq.

In order to avoid identifying differential expressions in genes of small average expression, fold change estimation is shrunken toward 0 for genes with insufficient information by employing an empirical Bayesian shrinkage. The procedure is as follows: (1) obtain the maximum likelihood estimates for the log fold changes from the GLM fit, then (2) fit a normal distribution with mean 0 to the estimates, and (3) use that as the prior for a second GLM fit. The maximum a posterior and the standard error for each estimate are the products of this procedure and will be used for the calculation of Wald statistics for DEA.

DESeq2 computes a threshold, η, to filter genes based on their average normalized expressions. The threshold is calculated for maximizing the number of genes with a user-defined false discovery rate. The authors claimed that this filtering step effectively controls the power of detecting DE genes. The null hypothesis becomes $\left|{\beta }_{gp}\right|\le \eta$ where ${\beta }_{gp}$ is the shrunken log fold change.

Finally, the method provides a way to diagnose outliers using the Cook's distance from the GLM within each gene,C_d Samples are flagged with C_d  99% quantile of an F distribution with degrees of freedom as the number of parameter, P, and the difference in the number of samples and the number of parameter, N - P. When there is a large number of replicates available, influential data can be removed without removing the whole gene; however, when there is a small number of replicates, the entire gene with influential points should be removed from the analysis to preclude bias. More details on DESeq2's features can be found in the study by Love et al. ⁹ In conclusion, DESeq2 is recommended by its authors as an improved solution to perform differential analysis because it adopts many competitive features.

Beta Binomial

Overview

A beta-binomial model is another alternative distribution to accommodate overdispersion.^11,12,42 The beta-binomial distribution has been used in the study by Baggerly et al. ¹¹ to account for both between-library and within-library variations. The authors assumed that the true proportion of gene g within a library i, ${\theta }_{gi},$ is library-specific and follows a beta distribution: ${\theta }_{gi}$ ∼ Beta(α, β), and that the count Y_gi given θ _gi follows binomial (m_i, θ _i ). Zhou et al. ¹² proposed a method, BBSeq, which also assumes a beta-binomial distribution and models the proportions of gene g within sample with a logistic regression. To estimate overdispersion parameters, BBSeq either treats the parameter as free and maximizing likelihood directly, or estimates the parameter through modeling the mean-overdispersion relationship.

Modeling

In a beta-binomial model, y_g is converted from the count of gene g in sample i, to proportion, θ _gi where ${\theta }_{gi}=\frac{y_{gi}}{\sum y_{gi}}$ . The model is constructed as: $$\log it\left(E{\theta }_{g.}\right))=\log \left(\frac{E\left({\theta }_{g.}\right)}{1-E\left({\theta }_{g.}\right)}\right)=X{\beta }_{g.}$$ (3.3.1) where β _g is a vector of the regression coefficients for sample covariates and is the parameter for hypothesis testing; θ _g. is a vector consisting of the proportion of gene g for sample i through N. With the beta-binomial distribution, we are no longer working with a log link but a logit link. θ _gi ∼ Beta with E(θ _gi ) = logit⁻¹(Xβ _g ) and var(θ _gi ) = ${\varphi }_{g}E\left({\varphi }_{gi}\right)\left(1-E\left({\varphi }_{gi}\right)\right)$ , where ϕ _g is the dispersion parameter. The hypothesis test is constructed as $H_0:{\beta }_{g{C}_1}=\dots ={\beta }_{g{C}_K}$ , where ${\beta }_{g{C}_K}$ denotes the estimated coefficient of the indicator variable with 1 for samples in group k and 0 otherwise.

Bayesian and Empirical Bayesian

Overview

RNA-seq DEA can be modeled in Bayesian framework using various parametric and nonparametric priors. Van de Wiel et al. ¹³ proposed a Bayesian method, ShrinkSeq, which either assumes an informative prior for the overdispersion such as the Dirac–Gaussian prior or estimates one with the empirical Bayesian approach. An empirical Bayesian approach differs from a fully Bayesian approach in that it borrows information from data to elicit priors for overdispersion parameters. For estimating posteriors, Van de Wiel and others ¹³ adapted the use of integrated nested Laplace approximations, a method that only considers marginal posteriors, but adds a direct maximization of marginal likelihood to allow information sharing from joint posteriors. They further suggested that the use of informative priors for shrinkage, as in ShrinkSeq, can ensure stability and accommodate multiplicity correction. They also suggested that shrinkage should be applied not only to overdispersion parameters but also to the regression coefficient parameters. baySeq, proposed by Hardcastle and Kelly, ¹⁴ constructs the data with tuples grouping genes together based on the study of interest. The distribution of a tuple shares the parameters of some prior distribution so that one can consider many hypotheses for testing beyond two group comparison. The method assumes a negative binomial distribution from the data. baySeq first estimates the empirical distribution on the set of parameters for null and alternative models with the quasi-likelihood approach. Then, it estimates the prior probabilities starting from a prior followed by an iterative process updating the priors until convergence. The authors suggested using a log posterior probability ratio of DE for DEA and noted that the posterior probability of DE for each individual model can be conveniently summed up for hypothesis testing.

Modeling

A Bayesian GLM for RNA-seq can be set as: $$Y_{gi}{}^{\underset{\_}{\underset{\_}{d}}}{F}_{{\mu }_{gi},{\gamma }_g},$$ (3.4.1) where γ _g is a vector of parameters not in the regression. The model is in fact flexible in that F can be negative binomial or other distributions. Suppose F follows a negative binomial distribution, then y_gi ∼ Poisson $\left({\mu }_{gi}\right);{\mu }_{gi}$ follows a gamma: ${\mu }_{gi}$ ∼ Gamma( $e^{{\eta }_{gi}},{\gamma }_g$ ), where ${\eta }_{gi}$ and γ _g are hyperparameters and ${\eta }_{gi}=X{\beta }_g={\beta }_{g0}+\sum _{p=1}^P{\beta }_{gp}{x}_{ip}.x_{ip}$ is the value of the pth covariate for sample i, such as ${\beta }_{g1}$ in a two-group comparison. With g(·) as a link function, ${\mu }_{gi}=g^{-1}\left({\eta }_{gi}\right)$ . The conditional posterior distribution for β is proportional with its prior: $$P\left(\beta |{\gamma }_g,y_{gi}\right)\propto P\left(\beta \right)\prod \frac{\exp {\left(X\beta \right)}^{y_{gi}}}{1+\exp {\left(X\beta \right)}^{y_{gi}+{\gamma }_g}}$$ (3.4.2)

Each parameter has its respective informative prior and one has to specify priors conditional on the model of interest as well as the prior itself to reach the posterior probability. For testing, a null hypothesis of β _g ≤ prior under the null is used.

Algorithm Overview 6: Van de Wiel et al.'s ¹³ ShrinkSeq

ShrinkSeq assumes that α is the unknown hyperparameter from a collection of all unknown hyperparameter vectors A. It uses a direct maximization of the marginal likelihood method for the estimation of A; this method is a modified version of INLA. ⁴³ The procedure of finding a is shown below and is said to be analogous to the EM algorithm:

1. 1.

   Initiate l = 0 and α⁽⁰⁾ _b for b = 1, …, B.

2. 2.

   Use INLA to estimate posteriors ${\pi }_{A\left(1\right)}\left(\varphi |Y_g\right).$

3. 3.

   Obtain ${\alpha }_b^{\left(l+1\right)}$ for b = 1, …, B with ML'.

4. 4.

   Iterate from step 2 until convergence.

Notes: let b be the number of informative priors and α^(l) _b be the bth element of A^(l) at iteration l; let ${\pi }_{A\left(1\right)}$ be the posterior of θg condition on data Y _g with A^(l) as the current estimate of A. ML' is ${\alpha }^{ML:,\left(l+1\right)}$ = argmax_α $\sum _{S=1}^S\log \left({\pi }_{\alpha }\left(z_{s,A^{\left(1\right)}}\right)\right)$ , where this is the prior log-likelihood at $z_{A^{\left(1\right)}}$ and s is a large independent sample set from ${\pi }_{A^{\left(1\right)}}^{EmpBayes}(\theta )$ ; ML' has the same mechanism as the maximum likelihood.

Dirac–Gaussian and Gaussian–Dirac–Gaussian mixture priors: $$\pi \left(\beta \right)=p_0{\delta }_0+\left(1-p_0\right)N\left(\beta ;0,{\tau }^2\right),$$ (3.4.a) $$\pi \left(\beta \right)=p_{-1}N\left(\beta ;{\mu }_{-1},{\tau }_{-1}^2\right)+p_0{\delta }_0+p_{1}N\left(\beta ;{\mu }_1,{\tau }_1^2\right),$$ (3.4.b)

The subscripts of p, ie, ™1, 0, and 1, indicate the locations. For example, Dirac mass on 0 is denoted as δ₀. Considering the p as probability where p_-1,p₀, and p₁ sum up to 1, then $p_0=1-p_{-1}-p_1.{\mu }_{-1}<0,{\mu }_1>0.$ Priors with positive mass on zero were intentionally selected because it reflects the non-DE condition. For more details on priors, please refer to the study by Van de Wiel et al. ¹³

Algorithm Overview 7: Hardcastle and Kelly's ¹⁴ baySeq

The tuple system in baySeq is as follows. Let a model be denoted as M. E refers to a set of models described by the data, {E₁ … E_l}. κ represents the set of parameters for each model, M, ie, {θ₁ … θ _l }. Let q be the index of each underlying distribution for model 1, …, l. An example would be that samples in groups 1, 2, and 3 (C₁, C₂, C₃) are grouped together in a way that groups 1 and 2 are equivalently distributed and group 3 stands alone: M = $\{A_{i\in C_1},A_{i\in C_2}\},\left\{A_{i\in C_3}\right\}$ where A is the sample. D_t is the data in tuple $t:\left\{\right\{y_{1t}\mathrm{\dots }{y}_{it}\mathrm{\dots }{y}_{n_{t}t}\},\{m_1\mathrm{\dots }{m}_i\mathrm{\dots }{m}_{n_t}\left\}\right\},$ which is the count in tuple t for sample i, m_i is the library size. The posterior probability of model given data is: $$P\left(M|D_t\right)=\frac{P\left(D_t|M\right)P\left(M\right)}{P\left(D_t\right)}$$ (3.4.c) $$P\left(\left|D_t\right|M\right)=\int P\left(D_t\left|\kappa \right|M\right)d\kappa$$ (3.4.d)

Suppose that a sample A_i is in the set E_q where the count of this sample at a particular tuple t is y_it, which follows a negative binomial $\left({\mu }_{it},{\phi }_q\right)\left({\theta }_q=\left({\lambda }_q,{\phi }_q\right)\right)$ . The mean count μ _it is a product of the library size scaling factor, m_i, and the proportion of reads in set E _q , λ _q . We have: $$\begin{array}{l}P\left(D_t|\kappa ,M\right)=P\left(y_{it}|m_i,{\theta }_q\right)\\ =\frac{\Gamma \left(y_{it}+{\varphi }_q^{-1}\right)}{\Gamma \left({\varphi }_q^{-1}\right)\Gamma \left(y_{it}+1\right)}{\left(\frac{1}{1+{\mu }_{it}{\varphi }_q}\right)}^{{\varphi }_q^{-1}}{\left(\frac{{\mu }_{it}}{{\varphi }_q^{-1}+{\mu }_{it}}\right)}^{y_{it}}\end{array}$$ (3.4.e) baySeq first estimates the empirical distribution on the set of parameters for null and alternative models through sampling from a negative binomial distribution and a quasi-likelihood approach. ³⁸ Then, it estimates the prior probabilities starting from a prior followed by an iterative process updating the priors until convergence. For detailed steps, please refer Hardcastle and Kelly. ¹⁴ Hypothesis testing can be easily denoted with the tuple system, for instance a two-group case, $$\begin{array}{l}{H}_0\left(non-DE\right):\left\{A_{i\in C_1,}{A}_{i\in C_2}\right\}\\ H_0\left(DE\right):\left\{A_{i\in C_1}\right\}and\left\{A_{i\in C_2}\right\}\end{array}$$

Nonparametric

Overview

In this section, we discuss two nonparametric methods for RNA-seq DEA by Li and Tibshirani ¹⁵ and Tarazona et al. ¹⁶ In SAMseq, Li and Tibshirani ¹⁵ calculated a modified two-sample Wilcoxon statistic using the ranked counts for two-group comparison. ⁴⁴ The authors proposed two resampling strategies for producing equal sequencing depths of the samples: downsampling and Poisson sampling, and also suggested that ties can be broken by inserting a small random number in resampling. NOISeq by Tarazona et al. ¹⁶ first used pseudo-counts corrected by the library size under two conditions (K=2) to calculate log-ratio (M) and absolute value of difference (D). Then, a test statistic is derived from M and D with a null hypothesis of no differential expression; in other words, M and D are no different than random variables either estimated from the real or simulated data.

Modeling

The two nonparametric methods discussed here are explained separately in the test boxes, as they each has a unique model setup.

Algorithm Overview 8: Li and Tibshirani's ¹⁵ SAMseq

To use SAMseq, one ranks the counts of gene g across samples and denotes the ordered counts as y'_g1 … y' _gN . If needed, resampling strategy may be used to fulfill the requirement of equal sequencing depths of samples in Wilcoxon test.

In the case of a sufficient minimal sequencing depth, the authors proposed a downsampling strategy where one first identifies the smallest sequencing depth, denoted as m_min, where m_min = min(m₁, …, m_N) and keeps this list f counts while resampling lists of counts for all other samples with the sequencing depth, m_min. Every count is randomly sampled with a success probability of m_min/m_i and failure probability of its complement, ie, the resampled count is $${y^{\prime }}_{ij}\sim binomial\left(y_{it},\frac{m_{\min }}{m_i}\right).$$ (3.5.a)

In the case of an insufficient minimal sequencing depth, Li and Tibshirani ¹⁵ introduced Poisson sampling strategy, wherein they employed the geometric mean of the sequencing depths for all samples: $${y^{\prime }}_{ij}\sim Poisson\left(\frac{\overline{m}}{m_i}{N}_{ij}\right),$$ (3.5.b) where y'_iJ is resampled data and $\overline{m}={\left(\prod _{i=1}^N{m}_i\right)}^{1/N}$ . Small random numbers are introduced into the resampling process to break ties, as well as multiple resampling to ensure stability. Poisson sampling is generally preferred based on the simulation. ¹⁵ In cases where m_i is unknown, one could use normalization methods to estimate. Differential expression of gene g is identified based on a comparison of the ranks of gene g between the two sample groups.

Algorithm Overview 9: Tarazona et al.'s ¹⁶ NOISeq

In NOISeq, for each $y_{gk\left(1\right)},$ the count of gene g in sample i from group k, the correction method for library size, m_k(i), is the sum of counts over all genes for the ith sample replicate in condition k. Let m_k(i) be simplified as m _i . One would work with pseudocounts (after normalization) formulated $as:{\tilde{y}}_{gk\left(i\right)}\times {10}^6/m_i.$

With the pseudocounts, the log ratio (L) and the absolute value of difference (D) are calculated. ${\tilde{y}}_{gk}$ is summarized over ith samples, a.k.a. ${\tilde{y}}_{gk}={{\sum }^{\text{}}}_{i\in C_k}{\tilde{y}}_{gi}.L_g={\log }_2(\frac{{\tilde{y}}_{g{C}_1}}{{\tilde{y}}_{g{C}_2}})$ and $D_g=|{\tilde{y}}_{g{C}_1}-{\tilde{y}}_{g{C}_2}|$ , where C₁ and C₂ denote group 1 and 2, respectively. Zero counts are replaced by 0.5 or by mid(0, normalized minimum expression) when calculating L_g. Samples with only zeros are dropped.

Null hypothesis: L and D values are no different than noise if no DE. Probability distribution for random variables L^* and D^* are either estimated from real data or simulated data and are used for the noises. One then obtains the probability of DE as: $$\begin{array}{c}P\left(D{E}_g=1|{\tilde{y}}_{g{C}_1}-{\tilde{y}}_{g{C}_2}\right)\\ =P\left(D{E}_g=1|L_g=l_g,D_g=d_g\right)\\ =P\left(|L*|<\left|l_g\right|,D*<d_g\right)\end{array}$$ (3.5.c)

DE _g equals 1 when gene g is DE. Note that log ratio is in absolute term because either direction indicates DE. See the study by Tarazona et al. ¹⁶ , for more details.

Statistical Testing

After performing parameter estimation for a statistical model, significance of differential expression can be assessed comparing the expression of gene g among K groups. Assume that ${\lambda }_{gk\left(i\right)}$ is the expression level of gene g in sample i belonging to sample group k. ϕ _g is the dispersion parameter. DE tests are proposed below for the null hypothesis (H₀): $${\lambda }_{g1}=\mathrm{\dots }={\lambda }_{gK}.$$

In parametric regime, one can employ classic log-likelihood ratio test. $$L{R}_g=\frac{2\left(l_g\left(\widehat{\lambda },Y_g\right)-l_g\left({\widehat{\lambda }}^0,Y_g\right)\right)}{{\widehat{\varphi }}_g}~F_{K-1,N-K}$$ (4.1)

In absence of overdispersion, $$L{R}_g=2\left(l_g\left(\widehat{\lambda },Y_g\right)-l_g\left({\widehat{\lambda }}^0,Y_g\right)\right)~X^2\left(K-1\right)$$ (4.2) where l_g denotes the log-likelihood function for the gth gene; $l_g\left(\widehat{\lambda },Y_g\right)$ and $l_g\left({\widehat{\lambda }}^0,Y_g\right)$ denote the MLE of biological and experimental effects under the full model and null model, respectively.

An exact test for negative binomial, analogous to the Fisher's exact test, is used by methods, such as edgeR and DESeq. By conditioning on the total sum, one can calculate the probability of observing counts as extreme or more extreme than what is really obtained, resulting in an exact P-value. Note that a sum of gene counts from all replicates in each group that is either too large or too small indicates a differential expression, so a two-sided test is used.

A score statistic is used by PoissonSeq, which tests for the significance of the association of gene g with expression of groups. In the context of gene count with unknown dispersion parameters, a score test is as follows: $$S_g=\sum _{k=1}^K\sum _{i\in C_k}\frac{w_g{\left(y_{gi}-{\widehat{\mu }}_{gi}\right)}^2}{{\varphi }_{g}v\left({\widehat{\mu }}_{gi}\right)}~F_{K-1,N-K}$$ (4.3)

where w_g is a known weight, ${\widehat{\mu }}_{gi}$ is estimated by MLE under the null hypothesis, and $v\left({\widehat{\mu }}_{gi}\right)$ is the variance function of ${\mu }_{gi}.$ .

Wilcoxon statistic is a rank-transformed version of t-statistics, used by the nonparametric method, SAMseq: $$W_g=\sum _{i\in C_k}{r}_{gi}-r_0,$$ (4.4) where r_gi is the rank of y_gi across samples and r₀ = $\left(\sum I_{\left(i\in C_k\right)}\right)\left(N+1\right)/2$ (r₀ is used to make E(W_g) = 0). W_g > 0 identifies that gene g is overly expressed in group k.

Under a Bayesian or empirical Bayesian framework, methods like baySeq use posterior likelihood of the DE model per gene to identify differential expression: $$P\left(M_{H_0}|Y_g\right)=\frac{P\left(Y_g|M_{H_0}\right)P\left(M_{H_0}\right)}{P\left(Y_g\right)},$$ (4.5) where M denotes a model. Posterior probability of DE to non-DE ratio is often used.

The choice of a testing strategy is a decision that often depends on the chosen method and other factors such as sample size. With a small sample size, the large-sample approximations based on the Wald test, score test, and likelihood ratio test are questionable and an exact test is usually preferred. ³⁶ We summarize testing strategies that are plausible for each method in Table 2.

Finally, almost all the methods we mentioned in this paper use standard approaches for multiple hypothesis correction to control false discoveries.^45,46 PoissonSeq is an exception that builds its own estimation of false discovery rate (FDR) from a permutation test. Permutation test calculates a score test per gene, S_g, for H_0g vs H_ag, each time when the outcome is permuted. For B permutations, the same procedure is applied to calculate null statistics ${\text{S}}_g^{0b}$ for b = 1 ··· B. The permutation P-value is: $$p_g=\sum _{b=1}^B\sum _{i=1}^N\frac{I\{{\text{S}}_g^{0b}>S_g\}+1}{N\times B+1}$$ (4.6)

For Bayesian methods, since posterior probabilities are computed, Bayesian FDR or local FDR are conveniently used. Local false discovery rate (1FDR_g) is simply the posterior probability ${\pi }_{0g}:$ : $$lFD{R}_g=P\left(M_{H_0}|Y_g\right)=P_0/\left(P_0+P_{\alpha }\right),$$ (4.7) where $P_0={\int }_{-\infty }^{△}P\left(Y_g|{\lambda }_g=\lambda \right)\pi \left(\lambda \right)d\lambda ,$ and $P_1={\int }_{△}^{\infty }P\left(Y_g|{\lambda }_g=\lambda \right)\pi \left(\lambda \right)d\lambda ,$ δ denotes prior. Bayesian false discovery rate (BFDR) is calculated as: $$BFDR\left(t\right)=\frac{\sum _{g=1}^{G}lFD{R}_g\times I\{{\pi }_{0_g}<t\}}{\sum _{g=1}^{G}I\{{\pi }_{0_g}<t\}}.$$ (4.8)

Note that $I\{{\pi }_{0_g}<t\}=I\{{\pi }_{1_g}\ge t\}$ for small t of interest.

Conclusion

RNA-seq data analysis is a relatively new and rapidly growing research area. The statistical model used for sequencing data has been evolving. The first proposed Poisson distribution has become obsolete because it fails to accommodate commonly-observed overdispersion in RNA-seq data. In a parametric framework, the negative binomial distribution is the most common assumption for modeling the marginal distribution due to the technical and biological variations.^8,9,33,36 Other available methods that account for overdispersions include the generalized Poisson distribution, ³⁵ negative binomial power distribution, ¹⁰ and beta-binomial distribution,^11,12 as well as nonparametric models^15,16 and Bayesian methods.^13,14 Table 2 summarizes all the reviewed methods in this paper.

For readers who are interested in the performance evaluation and method comparison of the available methods, they can refer to the original paper as well as the body of literature on this issue. For instance, in the study by Seyednasrollah et al. ²² , DESeq has been recommended as one of the most robust methods and caution is advised when dealing with a small number of replicates regardless of which method is being used. Similarly, Soneson and Delorenzi ²¹ also advise caution when interpreting results drawn from a small number of replicates and show that SAMseq surpasses many other reviewed methods. In the study by Rapaport et al. ²³ , DESeq, edgeR, and baySeq, which all assume a negative binomial model, have better specificity, sensitivity, and control of false positive errors than other nonnegative binomial models. As the technology continues to improve and the empirical data accumulate, more compelling statistical modeling for RNA-seq data can be expected.