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

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 catego-rized 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 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 testLikelihood 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 dataTag-seq data	Score test	Accommodate multiple covariate types, computationally efficient	Transformation power depends only on gene expression, libraries are totally exchangeable	Marioni RNA-seq dataTag-seq data
Negative binomial	edgeR	Robinson et al.7	SAGE dataRNA-seq data	Exact testLikelihood ratio test	Separate biological from technical variations	Limited to pairwise comparison	SAGE dataFly RNA-seq data
	DESeq	Anders and Huber8	Tag-seq dataRNA-seq dataChIP-seq data	Exact testLikelihood ratio test	Extend edgeR by allowing more general, data-driven relationship of mean and variance	Limited to pairwise comparison	Neural stem cell Tag-seq dataYeast RNA-seq dataHapMap ChIP-seq dataFly RNA-seq data
	DESeq2	Love et al.9	Tag-seq dataRNA-seq dataChIP-seq data	Wald testLikelihood ratio test	Improve upon DESeq for better gene ranking, allow hypothesis tests above and below threshold	Limited to pairwise comparison	Fly RNA-seq dataMouse 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 testLike lihood 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 dataCAGE 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 dataCAGE data
	baySeq	Hardcastle and Kelly14	Small RNAs data	Evaluating posterior probability for inference	Involve multiple comparison, accommodate different sample size	Computationally intensive but allow parallelization	Trans-acting small RNAs
Non-parametri	SAMseq	Li and Tibshirani15	RNA-seq dataTag-seq datamiRNA-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 pairwise 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 × N matrix Y. Each entry y (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 is a nonnega-tive integer representing the number of reads mapped to gene g in sample i. For succinctness, we also use notations “·”for summations, eg, and .

We use X to represent an N × P design matrix, where P is the number of covariates. For instance, x 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 represents the collection of indices of the samples in group k (k = 1, …, K), that is, C = {i: x = 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.24 One major source of nonbiological effects is sequencing depth, which can be adjusted by rescaling the sequencing counts with factors that mimic sequencing depth.25 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),27 DESeq size factor,28 and quantile-based normalizations such as upper quartile normalization,18 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.29, Benjamini and Speed,30 and Hansen et al.31

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 Huber8
Quantile	Bolstad et al.32

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.17 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.18 proposed to use two other test statistics, exact test statistics and score test statistics, in the generalized linear model (GLM) framework. Li et al.6 proposed a method called Poisson-Seq, 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.4 developed an R package, DEGseq, to identify differentially expressed (DE) genes with an MA-plot-based approach. Langmead et al.5 incorporated cloud computing in their method called Myrna.

Modeling

In a Poisson model, one assumes that Y, the number of reads mapped to gene g in sample i, follows a Poisson distribution, y ∼ Poisson(µ). µ 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: and E(Y) = µ and Var(Y) = µ. The association of µ with the same sample group can be described by a log-linear model as follows: where d represents the sequencing depth of sample i and is assumed for generality. Let β be the expression level of gene g and γ be the association of gene g with the covariate. For hypothesis testing, γ1 = … = γ = 0 indicates that the expression of gene g is not associated with the sample group. In the case of two sample group comparison, if γ = 0, then gene g is not DE between the two sample groups.

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 models34 and the generalized Poisson model.35

Negative binomial

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 Smyth33 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,36 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.40 Anders and Huber8 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.9 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 others10 proposed a method, NBPSeq, using a negative binomial power distribution instead of a regular negative binomial distribution. They hypothesized that , and ϕ is common across genes while α helps to accommodate the overdispersion. ϕ and α are estimated by maximizing conditional log-likelihood,41 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.10

The model setup for negative binomial is to assume y ∼ negative binomial (µ, ϕ). The dispersion parameter, ϕ, 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 where E(Y) = µ and . Hypothesis testing is set up as H0: 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.

Beta binomial

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.11 to account for both between-library and within-library variations. The authors assumed that the true proportion of gene g within a library i, θ, is library-specific and follows a beta distribution: θ ∼ Beta(α, β), and that the count Y given θ follows binomial (m, θ). Zhou et al.12 proposed a method, BBSeq, which also assumes a beta-binomial distribution and models the proportions of gene g within sample i 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.

In a beta-binomial model, y is converted from the count of gene g in sample i, to proportion, θ where . The model is constructed as: where is a vector of the regression coefficients for sample cova-riates and is the parameter for hypothesis testing; 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. θ ∼ Beta with E(θ) = logit−1 (β) and var(θ) = ϕ(θ)(1 –E(θ)), where ϕ is the dispersion parameter. The hypothesis test is constructed as , where β denotes the estimated coefficient of the indicator variable with 1 for samples in group k and 0 otherwise.

Bayesian and empirical Bayesian

RNA-seq DEA can be modeled in Bayesian framework using various parametric and nonparametric priors. Van de Wiel et al.13 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 others13 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,14 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.

A Bayesian GLM for RNA-seq can be set as: where γ 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 ∼ Poisson(µ); µ follows a gamma: , where η and γ are hyperparameters and . x is the value of the pth covariate for sample i, such as β1 in a two-group comparison. With g(·) as a link function, µ = g−1(η). The conditional posterior distribution for β is proportional with its prior:

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 β ≤ prior under the null is used.

Initiate l = 0 and for b = 1, …, B.

Use INLA to estimate posteriors

Obtain for b = 1, …, B with ML′.

Nonparametric

In this section, we discuss two nonparametric methods for RNA-seq DEA by Li and Tibshirani15 and Tarazona et al.16 In SAMseq, Li and Tibshirani15 calculated a modified two-sample Wilcoxon statistic using the ranked counts for two-group comparison.44 The authors proposed two resampling strategies for producing equal sequencing depths of the samples: downsampling and Pois-son sampling, and also suggested that ties can be broken by inserting a small random number in resampling. NOISeq by Tarazona et al.16 first used pseudo-counts corrected by the library size mk 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.

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

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 λ() is the expression level of gene g in sample i belonging to sample group k. ϕ is the dispersion parameter. DE tests are proposed below for the null hypothesis (H0):

In parametric regime, one can employ classic log-likelihood ratio test.

In absence of overdispersion, where l denotes the log-likelihood function for the gth gene; and 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: where w is a known weight, is estimated by MLE under the null hypothesis, and is the variance function of µ.

Wilcoxon statistic is a rank-transformed version of t-statistics, used by the nonparametric method, SAMseq: where r is the rank of y across samples and (r0 is used to make E(W) = 0). W > 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: 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.36 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, for H0 vs H, each time when the outcome is permuted. For B permutations, the same procedure is applied to calculate null statistics for b = 1 ⋯ B. The permutation P-value is:

For Bayesian methods, since posterior probabilities are computed, Bayesian FDR or local FDR are conveniently used. Local false discovery rate (lFDR) is simply the posterior probability π0: where and , and Δ denotes prior. Bayesian false discovery rate (BFDR) is calculated as:

Note that I{π0 < t} = I{π1 ≥ 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,35 negative binomial power distribution,10 and beta-binomial distribution,11,12 as well as nonparametric models15,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.22, 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 Delorenzi21 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.23, 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.

Algorithm Overview 1: Li et al.’s6 PoissonSeq

Li and others proposed PoissonSeq that assumes the hypotheses as follows. Under the null hypothesis where genes and covariates are not relevant, logμ_gi=logd_i+logβ_g,(3.1.a)where di 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 Ngi(0) in later equations: Ngi(0)=exp(log(d^_i)+log(β^_g)).Under the alternative hypothesis where genes and covariates, xi*, are relevant, logμ_gi=logd_i+logβ_g+γ_gxi*(3.1.b)where xi* would be I_(i∈C_k) when comparing two or multiple sample groups. The authors suggested using the maximum like to estimate β^_g, as a result β^_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: d^_i=∑_g∈Sy_gi∑_g∈Sy_g⋅.(3.1.c)They then estimated which genes belong to S by a Poisson goodness-of-fit statistic, ie, GOF_g=∑i=1N(y_gi−d^_iy_g⋅)^2d^_iy_g⋅.(3.1.d)S is set to be the genes whose GOFg values are in the (ε, 1 − ε) quantile of all GOFg values. Li and others used ε = 0.25 in their study.6The objective is to test H0: γ_g1=…⋅=γ_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 ∑k=1K[∑_i∈C_k(y_gi−Ngi(0))]^2∑_i∈C_kNgi(0)~χ^2(K−1).(3.1.e)

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. y~Poisson(μ),μ~gamma(α,β)P(y|μ)=μ^yexp(−μ)y!f(μ)=(Γ(α)β^α)^−1(μ^α−1exp(−μ/β))Then, P(y)=∫0∞P(y_gi|μ)f(μ)dμ=(y!Γ(α)β^α)^−1∫0∞μ^(y−α)−1exp(−μ(1+1/β))dμ=Γ(y+α)β^yy!Γ(α)(1+β)^y+αOne substitutes back μgi, ygi, α=ϕg−1, and β = μgi ϕg, a gamma–Poisson mixture can be viewed as a negative binomial, see Equation (3.2.1).

Algorithm Overview 3: Robinson and Smyth’s33,36 edgeR

In edgeR, µgi = miλgk(i) where mi is the ith library size and λ_gk(i)=∑i=1C_kλ_gi represents the proportion of the total reads that is gene g in group k and λ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(ϕ_g)=l_g(ϕ_g)+αl_C(ϕ_g)(3.2.a)where α is the weight given to the common likelihood, lC; the maximum estimator of WL(ϕg) is denoted by ϕ^gWL. An α has to be chosen such that ϕ^gWL coincides with an empirical Bayesian solution, ϕ^gB, the Bayesian posterior mean estimator of ϕg where ϕ^_g|ϕ~N(ϕ_g,τg2) and ϕ_g~N(ϕ_0,τ02) 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.33In the study by Robinson and Smyth,36 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 mi = m for i ∈ Ck where Ck = {i: k(i) = k}, ygi ∼ negative binomial (µgi = mλgk, ϕ) in group k and Ygi’s evidently become identically distributed, and the maximum likelihood estimator (MLE) λ^_gk(i) becomes ∑_i∈C_ky_gi∑_i∈C_km_i in group k. CML function for dispersion ϕ given z_k=∑i=1n_ky_ki was proposed. The function is as follows: l_C(ϕ)=∑g=1G1_g(ϕ)=∑g=1G∑k=1K[∑j=1n_klogΓ(y_ki+ϕ^−1)+logΓ(n_kϕ^−1)−logΓ(z_k+n_kϕ^−1)−n_klogΓ(ϕ^−1)](3.2.b)In the case of different mi in group k, the MLE of λgk(i) 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 lC(ϕ) to estimate an accurate estimate of ϕ. Specifically, let m*=(∏i=1Nm_i)^1N, 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 (m*λ, ϕ).Hypothesis testing is set up as H0: λg1 = λ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’s8 DESeq

The read count ygi is modeled by a GLM of negative binomial distribution with a log link: log(λ_gi)=∑p=1Px_ipβ_gp(3.2.c)The mean µgi is the proportion of reads for gene g in sample i, λgk(i), scaled by a normalization factor, mi. The variance σgi2 is μ_gi+mi2ν_gk(i), where ν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), mi, 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: m^_i=medianiy_gi(∏v=1Ny_gv)^1/N(3.2.d)Since λgk(i) 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. λ^_gk(i)=1M_k∑i:k(i)=ky_gim^_i,(3.2.e)where Mk is the total number of replicates for group k. The sample variances with the common scale are calculated as: w_gk=1M_k−1∑i:k(i)=k(y_gim^_i−λ^_gk(i))^2(3.2.f) z_gk=λ^_gk(i)M_k∑i:k(i)=k1m^_i(3.2.g)In the case of a sufficiently large number of Mk, one can see wgk – zgk as the unbiased estimator of the raw variance νgk. In the case of a small number of Mk, local regression for a smooth function wk(λ) on the graph of (λ^_gk(i),w_gk) was suggested so that w_k(λ^_gk(i))−z_gk would be the estimate for the raw variance. More details are in the study by Anders and Huber.8

Algorithm Overview 5: Love et al.’s9 DESeq2

DESeq2 allows the normalization factors to be gene specific (mgi), rather than being fixed across genes (mi). The estimation of mgi is implemented in their new R packages.9When 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 userdefined false discovery rate. The authors claimed that this filtering step effectively controls the power of detecting DE genes. The null hypothesis becomes |βgp| ≤ η where β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, Cd. Samples are flagged with Cd ≥ 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.9 In conclusion, DESeq2 is recommended by its authors as an improved solution to perform differential analysis because it adopts many competitive features.

Algorithm Overview 6: Van de Wiel et al.’s13 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.43 The procedure of finding α is shown below and is said to be analogous to the EM algorithm: Initiate l = 0 and αb(0) for b = 1, …, B. Use INLA to estimate posteriors π_A^(l)(θ|Y_g) Obtain αb(l+1) for b = 1, …, B with ML′. Iterate from step 2 until convergence. Notes: let b be the number of informative priors and αb(l) be the bth element of A(l) at iteration l; let π_A^(l) be the posterior of θg condition on data Yg with A(l) as the current estimate of A. ML′ is α^ML^′,(l+1)=argmax_α∑s=1Slog(π_α(z_s,A^(l))) where this is the prior log-likelihood at z_A^(l) and s is a large independent sample set from π_A^(l)(θ)EmpBayes; ML′ has the same mechanism as the maximum likelihood.Dirac–Gaussian and Gaussian–Dirac–Gaussian mixture priors: π(β)=p_0δ_0+(1−p_0)N(β;0,τ^2),(3.4.a) π(β)=p_−1N(β;μ_−1,τ−12)+p_0δ_0+p_1N(β;μ_1,τ12),(3.4.b)The subscripts of p, ie, −1, 0, and 1, indicate the locations. For example, Dirac mass on 0 is denoted as δ0. Considering the p as probability where p−1, p0, and p1 sum up to 1, then p0 = 1 − p−1 – p1·µ−1, <0, µ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.13

Algorithm Overview 7: Hardcastle and Kelly’s14 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, {E1… El}. κ represents the set of parameters for each model, M, ie, {θ1 … θ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 (C1, C2, C3) together in a way that groups 1 and 2 are equivalently distributed and group 3 stands alone: M={A_i∈C_1,A_i∈C_2},{A_i∈C_3} where A is the sample. Dt is the data in tuple t:{{y_1t…y_it…y_n_tt},{m_1…m_i…m_n_t}}, which is the count in tuple t for sample i, mi is the library size. The posterior probability of model given data is: P(M|D_t)=P(D_t|M)P(M)P(D_t)(3.4.c) whereP(D_t|M)=∫P(D_t|κ,M)P(κ|M)dκ(3.4.d)Suppose that a sample Ai is in the set Eq where the count of this sample at a particular tuple t is yit, which follows a negative binomial(µit, ϕq) (θq = (λq, ϕq)). The mean count µit is a product of the library size scaling factor, mi, and the proportion of reads in set Eq, λq. We have: P(D_t|κ,M)=P(y_it|m_i,θ_q)=Γ(y_it+ϕq−1)Γ(ϕq−1)Γ(y_it+1)(11+μ_itϕ_q)^ϕq1(μ_itϕq−1+μ_it)^y_it(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.38 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.14 Hypothesis testing can be easily denoted with the tuple system, for instance a two-group case, H_0(non-DE):{A_i∈C_1,A_i∈C_2}H_a(DE):{A_i∈C_1}and{A_i∈C_2}

Algorithm Overview 8: Li and Tibshirani’s15 SAMseq

To use SAMseq, one ranks the counts of gene g across samples and denotes the ordered counts as yg1′⋯ygN′. 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 mmin, where mmin = min(m1, …, mN) and keeps this list of counts while resampling lists of counts for all other samples with the sequencing depth, mmin. Every count is randomly sampled with a success probability of mmin/mi and failure probability of its complement, ie, the resampled count is yi′_j~binomial(y_ij,m_minm_i).(3.5.a)In the case of an insufficient minimal sequencing depth, Li and Tibshirani15 introduced Poisson sampling strategy, wherein they employed the geometric mean of the sequencing depths for all samples: yij′~Poisson(m¯m_iN_ij),(3.5.b)where yij′ is resampled data and m¯=(Πi=1Nm_i)^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.15 In cases where mi 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.’s16 NOISeq

In NOISeq, for each ygk(i), the count of gene g in sample i from group k, the correction method for library size, mk(i), is the sum of counts over all genes for the ith sample replicate in condition k. Let mk(i) be simplified as mi. One would work with pseudocounts (after normalization) formulated as: y˜_gk(i)=y_gk(i)×10^6/m_i.With the pseudocounts, the log ratio (L) and the absolute value of difference (D) are calculated. y˜_gk is summarized over ith samples, a.k.a. y˜_gk=∑_i∈C_ky˜_gi. L_g=log_2(y˜_gC_1y˜_gC_2) and D_g=|y˜_gC_1−y˜_gC_2|, where C1 and C2 denote group 1 and 2, respectively. Zero counts are replaced by 0.5 or by mid(0, normalized minimum expression) when calculating Lg. 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: P(DE_g=1|y˜_gC_1−y˜_gC_2)=P(DE_g=1|L_g=l_g,D_g=d_g)=P(|L*|<|l_g|,D*<d_g)(3.5.c)DEg 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.16, for more details.