Differential Gene Expression Analysis

We start with summarized measures of gene expression for the experiment, represented by a matrix of read or fragment counts. The rows of the matrix represents genes, $(g=1,\dots ,G)$ , and columns represent samples, $(i=1,\dots ,m)$ . Let Y_gi denote the count of RNA-seq fragments assigned to gene g in sample i. We assume that Y_gi follows a NB distribution with mean μ_gi and dispersion α_g, such that $\text{Var}(Y_{gi})={\mu }_{gi}+{\alpha }_g{\mu }_{gi}^2$ . The mean μ_gi is a product of a scaling factor s_gi and a quantity q_gi that is proportional to the expression level of the gene g. We follow the methods of Love et al. (2014) (link) to estimate α_g and s_gi sharing information across G genes, and consider estimates as fixed for the following. We fit a GLM to the count Y_gi for gene g and sample i,
$\begin{array}{l}{Y}_{gi}\sim \text{NB}({\mu }_{gi},{\alpha }_g)\\ {\mu }_{gi}=s_{gi}{q}_{gi}\\ \hspace{0.17em}\text{log}\hspace{0.17em}{q}_{gi}=X_{i,*}{\beta }_g\end{array}$
where X is the standard design matrix and β_g is the vector of regression coefficients specific to gene g. Usually X has one intercept column, and columns for covariates, e.g. indicators of the experimental conditions other than the reference condition, continuous covariates, or interaction terms. We consider design matrices where the first element of β_g is the intercept. For clarity, we partition the β_g into ${\beta }_g=({\beta }_{g0},{\beta }_{g1},\dots ,{\beta }_{gK})$ , where ${\beta }_{g0}$ is the intercept and β_gk, $k=1,\dots ,K$ is for kth covariate. The scaling factor s_gi accounts for the differences in library sizes, gene length (Soneson et al., 2015 (link)) or sample-specific experimental biases (Patro et al., 2017 (link)) between samples, and is used as an offset in our model.
In the GLM, we use the logarithmic link function. In the apeglm software, the estimated coefficients and corresponding SD estimates are reported on the same $\text{natural\hspace{0.17em}log}\hspace{0.17em}$ scale. The apeglm method can be easily called from DESeq2’s lfcShrink function, which provides LFC estimates on the ${\hspace{0.17em}\text{log}\hspace{0.17em}}_2$ scale. The apeglm method and software is generic for GLMs and can be used with other likelihoods. For example, it can be used for the Beta Binomial or zero-inflated NB model, as long as estimates for the additional parameters, e.g. dispersion or the zero component parameters, are provided. An example of apeglm applied to Beta Binomial counts, as could be used to detect differential allele-specific expression, is provided in the software package vignette.