Gamma-Poisson GLM - Seurat Blog

Introduction¶

The Gamma-Poisson (or Negative Binomial) distribution is a suitable model to describe variation accurately in count matrices of single-cell RNA-sequencing (scRNA-seq) data in both theory and practice Ahlmann-Eltze & Huber, 2021. It plays an essential role in normalizing count data in the SCTransform process before performing downstream analyses, such as principal component analysis (PCA) or clustering.

Although several packages implement this framework, glmGamPoi offers an outstanding and robust performance Ahlmann-Eltze & Huber, 2021.

Seurat command

This command is captured from the process of SCTransform.

fit <- glmGamPoi::glm_gp(data = umi,
                           design = as.formula(new_formula),
                           col_data = data,
                           offset = log_umi,
                           size_factors = FALSE)

Command explanation:

data: The count matrix, genes $\times$ samples.
design: The statistical model that fitted data. The model ~1, which treats all samples are in a same group, is used in this example.
col_data: The metadata of samples.
offset: Additional constants of each genes to adjust size factors. In this example, it is natural log of the total counts of each sample.
size_factors: Scaling factors to normalise counts across samples. In the commands, size factors are not applied.

Workflow¶

Overall, this workflow fits a raw count matrix to a Negative Binomial (NB) distribution. The primary target is to determine the model coefficients: $\beta$ , which represented as the log-scaled expected mean counts; and the dispersion parameter ( $\theta$ ) for each gene.

The workflow begins with an initial rough dispersion estimation, followed by estimating the model coefficients $\beta$ using the Maximum Likelihood method. After fixing $\beta$ , the maximum likelihood estimation of $\theta$ is conducted. These gene-specific dispersion estimates are then stabilized via a shrinkage process.

Rough dispersion estimation¶

\begin{cases} \begin{aligned} \xi &= \left( \frac{1}{N}\sum_{j=1}^{N}{\exp(o_j)} \right)^{-1} \\ \widehat{\theta}_{i} &= \max \left(\frac{var(C_i) - \xi mean(C_i)}{mean(C_i)^2}, 0 \right) \end{aligned} \end{cases}

(1)

Initially, a rough estimate of the dispersion for gene $i$ ( $\widehat{\theta}_{i}$ ) is computed using Equation (1). Here, $\xi$ is an offset correction factor, defined as the reciprocal of the average exponentiated offset ( $o_j$ ) across all $N$ samples. Since the offset in this context is the natural logarithm of the total counts for each sample, $\xi = (\overline{C})^{-1}$ , where $\overline{C}$ is the mean total count across samples.

The rough dispersion estimator is derived from the mean–variance relationship of the Gamma-Poisson distribution ( $\sigma^2 = \mu + \theta \mu^2$ ) Ahlmann-Eltze & Huber, 2021, adjusted by the factor $\xi$ and constrained to be non-negative.

Here, $C_i$ denotes the vector of counts for gene $i$ , while $var()$ and $mean()$ represent the variance and mean functions, respectively.

Beta estimation¶

\begin{cases} \begin{aligned} NC_{ij} &= \frac{C_{ij}}{\sum_{k=1}^{M}{C_{kj}}} \\ \widehat{\beta_i} &= \ln(mean(NC_i)) \end{aligned} \end{cases}

(2)

At the beginning, rough beta of each gene ( $\beta_i$ ) is estimated by Equation (2), defined as the logarithm of the gene-oriented mean of normalised counts ( $NC_{ij}$ ), which are proportion of count on total counts of each sample ( $M$ is the number of genes).

\beta_i^{(n+1)} = \beta_i^{(n)} + \frac{dl_{i}(\beta_i)}{dl'_{i}(\beta_i)}

(3)

Beta is estimated using the Newton-Raphson method Akram & Ann, 2015, which iterates Equation (3) until the value of beta converges. $dl_{i}$ is the derivative of likelihood function that count values of gene $i$ is of Negative Binomial (NB) distribution.

\begin{cases} \begin{aligned} \mu_{ij} &= \exp(\beta_i + o_j) \\ dl_{i}(\beta_i) &= \sum_{k=1}^{N}{\frac{C_{ik}-\mu_{ik}}{1 + \mu_{ik} \theta_i}} \\ dl'_{i}(\beta_i) &= - \sum_{k=1}^{N}{\frac{\mu_{ik}(1+C_{ik}\theta_i)}{(1 + \mu_{ik} \theta_i)^2}} \end{aligned} \end{cases}

(4)

Starting with an initial rough estimate of $\beta_i$ for gene $i$ , the sample-specific expected count $\mu_{ij}$ is updated in each iteration $n$ by combining the biological parameter ( $\beta_i$ ) and technical sampling offset ( $o_j$ ). Even though $\mu$ represents the theoretical mean of the NB distribution, incorporating these sample-specific offsets allows the model to account for variation in sequencing depth across the samples. This helps strip away library-size biases while computing the log-likelihood function and its derivative according to Equation (4).

Maximum likelihood of NB distribution

By assuming counts of each gene follows NB distribution, where variance quadratically increases by mean, data is utillized to estimated coefficents in the model by the maximum likelihood method.

\begin{cases} \begin{aligned} P(Y = y) &= \binom{y+r-1}{r-1} p^r (1-p)^y \\ p &= \frac{1}{1+\mu/r} \end{aligned} \end{cases}

(5)

Theoretically, the probability mass function (PMF) of NB distribution (5) estimates the probability of obtaining $y$ failures before archiving $r$ successes in a sequence of independent Bernoullian trials (two possible results, success or failure), with $p$ is the probability of succeeding any trial, and $\mu$ is the mean of number of failures Sinharay, 2010. In practical, since $y$ is a discrete value starts from 0, it is chosen to assign for $C_k$ without any biological meaning.

l(p,r) = \ln \left( \prod_{k=1}^{N}{ P(Y = C_k) } \right)

(6)

The log-likelihood function (6) converts the product of probabilities into a sum, making it analytically straightforward to differentiate. As mentioned earlier, the relationship between the dispersion parameter $\theta$ and the number of successes $r$ is $r = 1/\theta$ .

l(\mu,\theta) = \sum_{k=1}^N{\left( \ln \left[\binom{C_{k}+1/\theta-1}{1/\theta-1} \right] - \frac{1}{\theta}\ln(\mu\theta + 1) + C_k\ln(\frac{\mu\theta}{1+\mu\theta}) \right) }

(7)

Substituting $p$ according to (5), Equation (7) computes likelihood depending on the mean count of a gene. Maximum likelihood is found by rooting derivative of the likelihood function.

\frac{\partial l}{\partial \mu} = \sum_{k=1}^N{\frac{C_k - \mu}{\mu(1+\mu\theta)}}

(8)

Equation (8) provides the partial derivative of the log-likelihood with respect to $\mu$ .

\frac{\partial l}{\partial \beta} = \frac{\partial l}{\partial \mu} \cdot \frac{\partial \mu}{\partial \beta} = \sum_{k=1}^N{\frac{C_k - \mu}{1+\mu\theta}}

(9)

Hence the package glmGamPoi uses logarithm-scaled mean $\beta$ mentioned in (4) instead of count mean $\mu$ , it is transformed to paritial derivative by $\beta$ following Equation (9).

Overdispersion (theta) estimation¶

Hence using the likelihood function $l(\mu,\theta)$ to model NB distribution for each gene, the overdispersion parameter $\theta$ is estimated after determination of $\beta$ derived from $\mu$ . The vector of mean count $\mu_i$ of gene $i$ is deduced via the equation in (4), which stands for constant count mean across $N$ samples during optimizing $\theta$ .

l(\theta) = \sum_{k=1}^{N}{ \left[ \ln\Gamma(C_k + \theta^{-1}) - \ln\Gamma(\theta^{-1}) - \left( C_k +\theta^{-1} \right)\ln(\mu_k+\theta^{-1}) - \frac{1}{\theta} \ln(\theta) \right] }

(10)

The likelihood function $l(\theta)$ presented in Equation (10) is derived from general likelihood function of NB distribution (7). Here the factorial components in the combination term is substituted with Gamma functions to generalize the domain to continuous real numbers ( $\mathbb{R}$ ), while any additive terms independent of $\theta$ are omitted.

\frac{d l}{d \theta} = \frac{1}{\theta} \left[ \sum_{k=1}^{N}{ \left( -\frac{1}{\theta} \left( \psi(C_k + \theta^{-1}) - \psi(\theta^{-1}) \right) + \ln(1+\mu_k\theta) + \frac{C_k - \mu_k}{\mu_k + \theta^{-1}} \right) } \right]

(11)

To maximize the log-likelihood function, we first compute its derivative with respect to $\theta$ , as shown in Equation (11), where $\psi(x)$ denotes the digamma function (the first derivative of $\ln\Gamma(x)$ ).

- \underbrace{\frac{1}{\theta}\sum_{k=1}^{N}{ \left( \psi(C_k + \theta^{-1}) - \psi(\theta^{-1}) \right) } }_{D(\theta)} + \underbrace{\sum_{k=1}^{N}{ \left( \ln(1+\mu_k\theta) + \frac{C_k - \mu_k}{\mu_k + \theta^{-1}} \right) } }_{L(\theta)} = 0

(12)

Setting $dl/d\theta = 0$ yields the root-finding problem in Equation (12), which can be separated into two distinct components. The first part, $D(\theta)$ , accounts for the digamma terms and involves intensive computation. To eliminate redundant calculations for identical count values, glmGamPoi optimizes this step by utilizing a count frequency table Ahlmann-Eltze & Huber, 2021. Meanwhile, the remaining log-mean component, $L(\theta)$ , is more straightforward to compute.

G(\theta) = L(\theta) - D(\theta)

(13)

Structurally, the gradient $G(\theta)$ defined in Equation (13) governs the direction of the variation of $\theta$ . When $G(\theta) > 0$ , the log-likelihood slope is positive, indicating that $\theta$ must be increased. Conversely, when $G(\theta) < 0$ , the slope is negative, meaning $\theta$ must be decreased. A state of $G(\theta) = 0$ indicates that the log-likelihood has reached a local extremum.

To verify whether this stationary point corresponds to a local maximum or minimum. The second-derivative of the gradient is tested, if $d^{2}G/d\theta^2 < 0$ , the curvature is concave down, confirming that the extremum is a local maximum representing the Maximum Likelihood Estimate (MLE) of $\theta$ .

The Cox-Reid (CR) adjustment

The Cox-Reid (CR) adjustment Cox & Reid, 1987 is widely used in bulk and single-cell RNA-seq data analysis to correct for the negative bias that occurs when performing dispersion estimation via Maximum Likelihood Estimation (MLE) Love et al., 2014.

\begin{cases} \begin{aligned} w_{jj} &= \frac{1}{\theta + 1/\mu_j } \\ B &= M^TWM \\ b(\theta) &= \ln(\det(B)) \end{aligned} \end{cases}

(14)

The log-determinant of the Fisher information matrix for the linear model coefficients is denoted as $b(\theta)$ . Equation (14) outlines its calculation.

For a simple null model ~1, the design matrix ( $M$ ) contains only an intercept term (a column vector of ones). $W$ is a diagonal weight matrix whose elements $w_{jj}$ are determined by the mean $\mu_j$ and the dispersion $\theta$ as specified by Love et al. (2014).

l(\theta)_{\text{cr}} = l(\theta) - \frac{1}{2} b(\theta) F_{\text{cr}}

(15)

The Cox-Reid adjusted log-likelihood ( $l(\theta)_{\text{cr}}$ ) is defined in Equation (15), where the second term is the adjustment from the bias. This is the main likelihood function for further calculating and estimation. The scaling factor $F_{\text{cr}}$ (is set to 0.99) prevents the mathematical cancellation leading to exceptionally large values of $\theta$ in edge cases.

Shrinkage¶

Even though NB distribution is able to obtain flexible variance-mean relationship, it is proved virtually being uncertainty in parameter estimation. Quasi-likelihood (QL) comes as a solution to adjust dispersion Lund et al., 2012. The variance-mean equation of quasi-distribution by genes is shown as (16).

\sigma^2 = \theta_{QL}(\mu + \mu^2\theta_{\text{trend}})

(16)

Initially, the trend dispersion ( $\theta_{\text{trend}}$ ) is derived from the vector of likelihood estimated dispersion ( $\exp(1/\theta_{\text{ML}})$ ), which is performed local median regression to stablize dispersion. Specifically, the dispersion according to each mean count value ( $\mu_i$ ) would be performed weighted median Cormen et al., 2022 with neighbors (sorted by mean, analyzing 100 neighbors by default), where the weight range originates from the range of probability in the domain $[-3,3]$ of standard normal distribution.

\theta_{QL} = \frac{1+\mu\theta_{\text{ML}}}{1+\mu\theta_{\text{trend}}}

(17)

$\theta_{QL}$ is defined by (17), which derived from (16) and mean-variance equation of NB distribution, where likelihood dispersion $\theta_{\text{ML}}$ was smoothed into $\theta_{\text{trend}}$ .

\theta_{SQL} = \frac{\text{df}_0\tau^2_0 + \text{df}\theta_{QL}}{\text{df}_0 + \text{df}}

(18)

The shrunken quasi-likelihood overdispersion ( $\theta_{SQL}$ ) for each gene is estimated by Equation (18). Particularly, $\text{df}$ is the residual degree of freedom - DOF ( $\text{\#samples} - 1$ ), while $\text{df}_0$ and $\tau^2_0$ are prior parameters representing DOF and dispersion scaler respectively. These prior parameters are estimated by modeling the $\theta_{QL}$ values using a scaled Chi-squared distribution via a natural cubic spline framework. The process successfully yields a theoretical pseudo-sample size ( $\text{df}_0$ ) and a denoised central dispersion trend ( $\tau_0^2$ ). Ultimately, the resulting $\theta_{SQL}$ reflects the true amount of informative dispersion captured from the data.

Under the intercept-only model (~1), all samples are treated as a single group. The coefficient vector $\beta$ is re-estimated using the trend dispersion ( $\theta_\text{trend}$ ), which is then used to recalculate the expected mean counts.

Summary¶

The Gamma-Poisson model is a good fit for overdispersed scRNA-seq data. Ultimately, the function glm_gp() of the glmGamPoi package, which models NB distribution, targets essential regression coefficients (relative log-scaled count mean $\beta$ and overdispersion $\theta$ ), by maximizing their likelihood estimations.

References¶

Ahlmann-Eltze, C., & Huber, W. (2021). glmGamPoi: Fitting Gamma-Poisson Generalized Linear Models on Single Cell Count Data. Bioinformatics, 36(24), 5701–5702. 10.1093/bioinformatics/btaa1009
Akram, S., & Ann, Q. U. (2015). Newton raphson method. International Journal of Scientific & Engineering Research, 6(7), 1748–1752.
Sinharay, S. (2010). Discrete Probability Distributions. In International Encyclopedia of Education (pp. 132–134). Elsevier. 10.1016/B978-0-08-044894-7.01721-8
Cox, D. R., & Reid, N. (1987). Parameter Orthogonality and Approximate Conditional Inference. Journal of the Royal Statistical Society Series B: Statistical Methodology, 49(1), 1–18. 10.1111/j.2517-6161.1987.tb01422.x
Love, M. I., Huber, W., & Anders, S. (2014). Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2. Genome Biology, 15(12), 550. 10.1186/s13059-014-0550-8
Lund, S. P., Nettleton, D., McCarthy, D. J., & Smyth, G. K. (2012). Detecting Differential Expression in RNA-sequence Data Using Quasi-likelihood with Shrunken Dispersion Estimates. Statistical Applications in Genetics and Molecular Biology, 11(5). 10.1515/1544-6115.1826
Cormen, T. H., Leiserson, C. E., Rivest, R. L., & Stein, C. (2022). Introduction to algorithms. MIT press.