Differential Expression Analysis for RNA-Seq: An Overview of Statistical Methods and Computational Software <s> Algorithm overview 5: Love et al.'s 9 DESeq2 <s> Background ::: High throughput sequencing has become an important technology for studying expression levels in many types of genomic, and particularly transcriptomic, data. One key way of analysing such data is to look for elements of the data which display particular patterns of differential expression in order to take these forward for further analysis and validation. <s> BIB001 </s> Differential Expression Analysis for RNA-Seq: An Overview of Statistical Methods and Computational Software <s> Algorithm overview 5: Love et al.'s 9 DESeq2 <s> Motivation: A number of penalization and shrinkage approaches have been proposed for the analysis of microarray gene expression data. Similar techniques are now routinely applied to RNA sequence transcriptional count data, although the value of such shrinkage has not been conclusively established. If penalization is desired, the explicit modeling of mean–variance relationships provides a flexible testing regimen that ‘borrows’ information across genes, while easily incorporating design effects and additional covariates. ::: ::: Results: We describe BBSeq, which incorporates two approaches: (i) a simple beta-binomial generalized linear model, which has not been extensively tested for RNA-Seq data and (ii) an extension of an expression mean–variance modeling approach to RNA-Seq data, involving modeling of the overdispersion as a function of the mean. Our approaches are flexible, allowing for general handling of discrete experimental factors and continuous covariates. We report comparisons with other alternate methods to handle RNA-Seq data. Although penalized methods have advantages for very small sample sizes, the beta-binomial generalized linear model, combined with simple outlier detection and testing approaches, appears to have favorable characteristics in power and flexibility. ::: ::: Availability: An R package containing examples and sample datasets is available at http://www.bios.unc.edu/research/genomic_software/BBSeq ::: ::: Contact:yzhou@bios.unc.edu; fwright@bios.unc.edu ::: ::: Supplementary information:Supplementary data are available at Bioinformatics online. <s> BIB002 </s> Differential Expression Analysis for RNA-Seq: An Overview of Statistical Methods and Computational Software <s> Algorithm overview 5: Love et al.'s 9 DESeq2 <s> Next generation sequencing is quickly replacing microarrays as a technique to probe different molecular levels of the cell, such as DNA or RNA. The technology provides higher resolution, while reducing bias. RNA sequencing results in counts of RNA strands. This type of data imposes new statistical challenges. We present a novel, generic approach to model and analyze such data. Our approach aims at large flexibility of the likelihood (count) model and the regression model alike. Hence, a variety of count models is supported, such as the popular NB model, which accounts for overdispersion. In addition, complex, non-balanced designs and random effects are accommodated. Like some other methods, our method provides shrinkage of dispersion-related parameters. However, we extend it by enabling joint shrinkage of parameters, including those for which inference is desired. We argue that this is essential for Bayesian multiplicity correction. Shrinkage is effectuated by empirically estimating priors. We discuss several parametric (mixture) and non-parametric priors and develop procedures to estimate (parameters of) those. Inference is provided by means of local and Bayesian false discovery rates. We illustrate our method on several simulations and two data sets, also to compare it with other methods. Model- and data-based simulations show substantial improvements in the sensitivity at the given specificity. The data motivate the use of the ZI-NB as a powerful alternative to the NB, which results in higher detection rates for low-count data. Finally, compared with other methods, the results on small sample subsets are more reproducible when validated on their large sample complements, illustrating the importance of the type of shrinkage. <s> BIB003
The read count y gi is modeled by a GLM of negative bino mial distribution with a log link: The mean µ gi is the proportion of reads for gene g in sample i, λ gk(i) , scaled by a normalization factor, m i . The variance σ gi 2 is µ ν , 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 stabi lize 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 i , 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: 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. where M k is the total number of replicates for group k. The sample variances with the common scale are calculated as: 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 ν gk . In the case of a small number of M k , local regression for a smooth function w k (λ) on the graph of ( , 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 varia tion 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 genespecific 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 proce dure is as follows: (1) obtain the maximum likelihood esti mates 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 calcula tion 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 |β gp | # η where β gp is the shrunken log fold change. Finally, the method provides a way to diagnose outli ers 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 num ber 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 betabinomial model is another alternative distribution to accommodate overdisper sion.  BIB002  The betabinomial distribution has been used in the study by Baggerly et al.  to account for both between library and withinlibrary variations. The authors assumed that the true proportion of gene g within a library i, θ gi , is libraryspecific and follows a beta distribution: θ gi ∼ Beta(α, β), and that the count Y gi given θ gi follows binomial (m i , θ i ). Zhou et al. BIB002 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 overdis persion parameters, BBSeq either treats the parameter as free and maximizing likelihood directly, or estimates the para meter through modeling the meanoverdispersion relationship. Modeling. In a betabinomial model, y gi is converted from the count of gene g in sample i, to proportion, θ gi where The model is constructed as: where β g is a vector of the regression coefficients for sample cova riates 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 betabinomial distribution, we are no longer working with a log link but a logit link. θ gi ∼ Beta with E(θ gi ) = logit -1 (X β g ) and var(θ gi ) = φ g E(θ gi )(1 -E(θ gi )), where φ g is the dispersion parameter. The hypothesis test is constructed as H 0 : , where β gC K denotes the estimated coefficient of the indica tor 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. BIB003 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 informa tion from data to elicit priors for overdispersion param eters. For estimating posteriors, Van de Wiel and others BIB003 adapted the use of integrated nested Laplace approxima tions, a method that only considers marginal posteriors, but adds a direct maximization of marginal likelihood to allow information sharing from joint posteriors. They further sug gested that the use of informative priors for shrinkage, as in ShrinkSeq, can ensure stability and accommodate multi plicity 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, BIB001 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 quasilikelihood 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 RNAseq can be set as: 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.