Photo by Sven Mieke on Unsplash

Normalizing single cell RNA sequencing data — Pitfalls and Recommendations

The goals of a single cell RNA sequencing (scRNA-seq) project are often Identification of subpopulations and Differential Gene Expression Analysis. To avoid the ‘curse of dimensionality’, Highly Variable Genes (HVGs) are used for cluster analysis. Several studies have shown that selection of HVGs is sensitive to the choice of method used for normalization of raw count matrices.

Why Normalization?

Raw read counts cannot be directly used to compare gene expression between cells, as they are confounded by technical and ‘uninteresting’ biological variations. There are QC steps and other methods available to filter and regress uninteresting biological variations. While PCR amplification bias is often taken care by use of Unique Molecular Identifiers (UMIs), normalization is required to remove effects of other technical variations like differences in sequencing depth, cell lysis and reverse transcription efficiency.

The primary goal of a normalization process is to remove the influence of technical effects, with preservation of true biological heterogeneity. In a well normalized dataset, the variance of a gene should be independent of gene abundance and sequencing depth of a cell. The ‘truly’ differentially expressed genes should exhibit high variance across cell types, while housekeeping genes should exhibit low variance.

Normalization is thus, a critical pre- processing step that greatly impacts the downstream applications of scRNA analysis. Unfortunately, scRNA datasets are typically normalized using methods inherited from bulk RNA sequencing methods, which as we will see shortly are inappropriate due to the nature of technical variations and inherent complexity of these datasets. In this blog post we will see the limitations of global scaling methods in context to scRNA-seq analysis. We will also discuss the potential of recently introduced SCNorm and SCTransform normalization methods that have been specifically tailored for single cell analysis.

Global Scaling methods

Traditionally, raw expression counts across cells were normalized for sequencing depth with RPKM (Reads per Kilobase Million), FPKM (Fragments per Kilobase Million) or TPM (Transcripts per Million) methods. To understand how they work, please watch this video. While these methods work well for within-sample normalization, they have been widely recognized as unsuitable for differential expression analysis in-between samples. For example, consider a case where expression of two genes A and B are compared between two conditions -control and treated. Gene A is expressed at the same level in both conditions, while Gene B has a 2-fold higher expression in treated cells. TPM normalization converts absolute expression into relative expression and therefore, one might conclude that Gene A is differentially expressed, although this effect is only a consequence of its comparison with Gene B.

TPM normalization is unsuitable for differential expression analysis.

Alternative approaches were developed for between-sample normalizations; TMM (trimmed mean of M-values) and DESeq being most popular. Both strategies follow the same motivation: to bring cell-specific measures onto a common scale by standardizing a quantity of interest across cells, while assuming that most genes are not differentially expressed. The global scale factors calculated by these methods based on non-differentially expressed genes are applied to samples (one for each sample). All genes in cells of a sample are then scaled using the same factor. These methods demonstrate excellent performance in bulk RNA-seq, but they are compromised in the single-cell setting because of an abundance of zero-expression values. Also, these methods assume that the underlying RNA content is constant for all cells in a sample, and that a single scaling factor can be applied for all genes. Global scaling factor based normalization strategies therfore underperform for a typical zero inflated and highly heterogenous scRNA dataset.

Gene group based methods

To address the inherent problems with the global scaling approach, two interesting normalization methods have recently been introduced -SCnorm (2017) and SCTransform (Seurat package v3, 2019).

SCnorm

SCnorm is an R package available on Bioconductor. For every gene, SCnorm estimates the dependence of gene expression on sequencing depth by a quantile regression. Genes with similar dependence are then grouped, and a second quantile regression is used to estimate scale factors for each group. Every gene group is finally adjusted for sequencing depth using the group specific scale factor to produce normalized estimates of expression.

A: Raw Counts vs sequence depth, B: Global Scale Factor normalized vs sequence depth, C:SCnorm count vs sequence depth for 3 genes in a single cell dataset, edited from Bacher et al.

The above figure shows count-depth relationships for three genes from a single cell dataset. Panel A is unnormalized or raw expression counts. As is evident, a global scale factor based method (panel B) does a poor job of normalization as compared to SCnorm in panel C.

SCTransform

SCTransform is an R package available with Seurat v3. The method models UMI counts using a regularized negative binomial model to remove variation due to sequencing depth. Briefly, the method first constructs a generalized linear model (GLM) for each gene using sequencing depth as an independent variable and UMI count as response or dependent variable. The parameter estimates are then regularized (or adjusted) based on gene expression. A second round of negative binomial regression is applied using the regularized parameters. The output of this model (residuals) is the normalized expression levels for each gene.

The take-home messge here is that instead of using a constant factor to normalize all genes, both SCnorm and SCTransform methods learn gene- group specific factors. The factors individually address low, moderately and highly expressed genes, removing the influence of technical variation with preservation of true biological heterogeneity.

Summary

The choice of normalization method impacts the selection of Highly Variable Genes and thereby all downstream analyses of scRNA data. A direct application of bulk RNA normalization methods to scRNA datasets is inappropriate. It is worth to update the analysis pipeline by choosing SCNorm or SCTransform normalization methods and make the most of recent advances in technology.