SCTransform workflow

Introduction¶

Normalization and preprocessing are the key challenges affecting directly to downstream scRNA-seq results, in which technically sequencing bias is cut down while biological variation is reserved. An effective workflow would eliminate technical bias among cells/samples and keep genewise heterogeneity without being overwhelmed by dominant genes Hafemeister & Satija, 2019Choudhary & Satija, 2022.

Based on that, SCTransform is introduced as a probabilistic approach stabilizing variation in the count matrix Hafemeister & Satija, 2019. When there is a significant difference in cell expression originally, the SCTransform method is able to replace Seurat standard preprocessing workflow (normalization, feature selection, and scaling), mentioned in pbmc3k_tutorial.

Workflow¶

Overall, the raw count matrix is utilized for Negative Binomial (NB) model estimate to determine the expected mean and overdispersion parameters for each gene. After the estimating, the raw likelihood parameters are highly variable across genes, a regulation step is applied to stabilize variance and reduce the sampling noise.

Next, the contribution of each element in the count matrix to the Chi-square dependence between genes and samples is quantified via Pearson residuals. Subsequently, the count matrix is corrected by eliminating technical noise across samples, and reconstructing normalized expression counts from the stabilized residuals.

Finally, the Pearson residuals undergo post-processing, where highly variable features are identified, and confounding biological or technical covariates (such as cell viability indicated by the percentage of mitochondrial reads) can be regressed out, yielding a polished residual matrix ready for downstream analysis (e.g., PCA and clustering).

Fitting model¶

At the beginning, only genes obtaining an overdispersion factor ($\sigma^2 > \mu$) are used to train. A number of genes is randomly selected (default is 2,000) across their expression levels, and the total count of selected genes must be greater than 5 by default. Subsequently, their raw count matrix are modeled by Gamma-Poisson general linear model (Gamma-Poisson GLM), i.e Negative Binomial (NB) model, to properly obtain overdispersion. See Gamma-Poisson GLM for more detail about the model-fitting process.

Briefly, it uses the Maximum Likelihood Estimates (MLE) method to estimate coefficients of the model according to the count values of each gene. The coefficients include relative log-scaled count mean $\beta$, and overdispersion $\theta$. Note that the definition of overdispersion between sctransform and Gamma-Poisson GLM are reciprocal.

Next, theoretical overdispersion $\hat{\theta}$, which is derived from mean-variance relationship of NB model, shown as Equation (1). $\alpha$ is the ratio of expected overdispersion ($\hat{\theta}$) over observed overdispersion ($\theta$), which originates from estimation above. If $\alpha < 0.001$, the model for that gene is assumed to follow the Poisson distribution ($\sigma^2 = \mu$, and $\theta \rightarrow -\infty$).

Regularizing model¶

Geometric mean¶

For each gene, logarithmic geometric mean ($\mu_{\text{g}}$) computed by Equation (2), represents as a centric count without affected by outlier samples. With $C_j$ is the count value of sample $j$ in total $N$ samples, and $\epsilon$ (default is 1) is a small fixed number to avoid $\ln(0)$ Hafemeister & Satija, 2019.

Local outlier detection¶

Local outlier genes along the $\mu_{g}$-axis are detected based on the matrices of the log-scaled overdispersion factor ($F = 1+\mu_g/\theta$) and $\beta$ estimated from Gamma-Poisson GLM. These parameters serve for variance and mean outlier indicators, respectively.

Initially, the genes are assigned to bins according to their $\mu_{g}$ value. The binwidth is determined using a heuristic rescale from the optimal bandwidth featured by Sheather & Jones (1991) to ensure data is fairly distributed across bins (satijalab/sctransform#214).

For each evaluated matrix ($Y$ representing either $F$ or $\beta$), the outlier scores ($S$) are computed for the points within each bin following Equation (3), which measures the distance to the median and eliminates deviation by Median Absolute Deviation (MAD) to enable comparability. For enhancing reliability, the measurement is processed on two overlapping grids, offset by half a binwidth. A gene is flagged as an outlier within a specific matrix when its absolute score on both grids exceeds a threshold, which defaults to 10.

Ultimately, a gene is marked as an outlier if it is flagged by either of the evaluated matrices (variance or mean).

Kernel smoothing¶

The estimated parameters ($F$ and $\beta$) are sequentially smoothed across gene expression levels (log geometric mean) using the Nadaraya-Watson kernel regression estimator Nadaraya, 1964. To avoid overfitting, this smoothing is anchored by the subset comprising overdispersion genes which exclude the previously flagged outliers and Poisson-like genes.

The smoothed value $\overline{y}_j$ (in $F$ or $\beta$ matrix) of any gene $j$ is computed by weighted average of anchored genes, as defined in Equation (4). Particularly, the weight ($w_{ij}$) is the density determined by the Gaussian kernel function $K$ at the distance between the current gene $j$ and anchored gene $i$ in the $\mu_g$ level. In $m$ overdispersion genes, the employed subset starts from anchored gene $i_{\text{min}}$, which is nearest to the boundary, four times the bandwidth h, from the targeted $\mu_g$ range.

The interquartile of the distance distribution is assumed to be the range $[-0.25h, +0.25h]$, where $h$ is the optimal bandwidth of overdispersion genes by the method of Sheather & Jones (1991), see SJ bandwidth estimation for details. Noticeably, the optimized $h$ is tripled by default. All together, this ensures that anchor genes lying closer to the expected range obtained higher weight.

Based on the assumption, standard deviation of the distance distribution ($\overline{\sigma}$) used in Equation (4) computed by Equation (5) with $\Phi^{-1}(0.75)$ being the inverse cumulative distribution function of 75th percentile.

By the end, smoothed $\beta$ and $\theta$, which is recoverd from smoothed $F$, are utilized for the next process.

Pearson residuals¶

While Chi-square test evaluates the global association between two categorical variables in a contingency table, Pearson residuals isolate the explicit contribution of each individual component (here, each gene-cell pair) to that overall dependence. The computation of specific Pearson residual ($z_{ij}$) for gene $i$ in cell $j$ is outlined in Equation (6).

First, the expected count ($\mathbb{E}[c_{ij}]$) is calculated analogously to the expected frequency in Chi-square test with the gene-specific log-scaled mean $\beta_i$ estimated in the smoothing step and the total count $C_{j}$ for each cell.

Hence count data assumably follows NB distribution, the variance ($\sigma_{ij}^2$) for residuals is derived from $\mathbb{E}[c_{ij}]$ as the mean, and smoothed overdispersion $\theta_i$.

Finally, the specific Pearson residual is measured from the deviation between the observe ($c_{ij}$, the actual count) and the expected value ($\mathbb{E}[c_{ij}]$), and scaled by standard deviation $\sigma_{ij}$ to ensure comparability.

In practice, the variance $\sigma_{ij}^2$ is bound by a lower threshold $\sigma_{\text{min}}^2$ defined in Equation (7), where the deviation is expected to be always greater than the general median, and the maximum of residual is 5.

Count correction¶

Using the stabilized Pearson residuals, the raw count matrix is corrected following Equation (8).

The gene-specific baseline mean $\mu_c(i)$ and variance $\sigma_{c}^2(i)$ used for this correction step are calculated similarly to $\mathbb{E}[c_{ij}]$ and $\sigma_{ij}^2$ in Equation (6), excepting that specifying the median of log-scaled total count of samples.

Ultimately, the corrected count value ($\mathbb{C}[c_{ij}]$) is computed by the combination of the baseline mean count across samples, and the specific expectation deviation derived from the specific residual. Finally, the corrected values are polished by rounding, and the minimum floor set at 0.

Post-processing¶

Feature selection¶

The top variable features (genes) are determined by ranking the variance of Pearson residuals across cells.

In detail, after computing the raw residuals, they are practically clipped by the range $[-\sqrt{M/30}, \sqrt{M/30}]$ (where $M$ is the total number of cells) to mitigate the distorting impact of extreme outliers. Next, genewise variance of these clipped residuals is calculated and the genes are sorted to select the top variable features (capped at 3,000 features by default).

Residualization¶

Subsequently, based on the parameter vars.to.regress, Pearson residuals are residualized against the percentage of mitochondrial genes, which reflects cell survival status. This step eliminates additional unwanted variance arising from technical errors or cell damage.

This dependency is modeled linearly, as shown in Equation (9), where the Pearson residual vector of gene $i$ ($\mathbf{z}_i$) depends on a linear model of the mitochondrial gene percentage ($p_{mt}$) and an independent residual vector ($\mathbf{r}_i$).

Finally, the resulting residuals are mean-centered within each gene.

Summary¶

To reduce technical noise without losing statistical variance, sctransform fits and regularizes a negative binomial model to stabilize variance across genes. Technical bias across cells is removed via Pearson residualization while preserving true biological variation. Furthermore, top variable features are selected, and dependence on confounding covariates is regressed out.

Notably, the corrected count matrix is reconstructed directly from the initial Pearson residuals before post-processing. Consequently, the corrected counts remain independent of additional covariate regressions, thereby preserving the baseline biological variance of the raw counts.