cdev: a ground-truth based measure to evaluate RNA-seq normalization performance

Simulations

Simulation of RNA-seq counts was performed in the same manner as done by Law et al. (2014) and later by Evans, Hardin & Stoebel (2017). Specifically, baseline relative abundance was generated for n genes in m samples of two conditions, according to an empirical distribution obtained from real data (Law et al., 2014). In m/2 samples of condition 2, a k percentage of these values were then multiplied by the fold-change f, creating k × n DE genes. From these baselines, counts were generated by sampling from the corresponding negative binomial distribution. There are two types of simulations used in this study. In the homogeneous fold-change scenario, a single fold-change was applied for all DE genes. For a more realistic scenario, i.e., heterogeneous fold-change, random values of fold-change were sampled from the uniform distribution Unif(1, 10) for up-regulated genes, or Unif(−10,  − 1) for down-regulated genes, such that 50% of the DE genes are up-regulated, and 50% are down-regulated. This choice is arbitrary, since the distribution of fold-change in real data sets comes in all shapes and sizes, and even shifted by normalization.

A summary of different simulation scenarios were given in Table 1.

Compilation of RNA-seq data

RNA-seq data were collected from the Gene Expression Omnibus (GEO) database based on the following criteria: having ERCC spike-ins added at a constant concentration across a sufficiently large number of samples. The list of studies where these assays were originally performed is available in Table 2.

All RNA-seq assays from said studies were downloaded from the Sequence Read Archive (SRA) in FASTQ format and re-processed by a single workflow. In this workflow, the raw reads were aligned to the corresponding transcriptome using STAR v2.7 (Dobin et al., 2013), and quantified using RSEM v1.3 (Li & Dewey, 2011). Prior to alignment and quantification, each species-specific transcriptome was extended by a set of 92 ERCC spike-in sequences (Jiang et al., 2011) to create the software-specific indices. Metadata were extracted and compiled manually from all available information in GEO, SRA, the sequence identifiers in the FASTQ files, and original publication when necessary. The biological samples were then separated into individual data sets, such that every sample in the same set was added with the same amount of ERCC spike-ins being processed in the same experimental procedure (Table 3). The data sets from each project were then wrapped in R data packages (Table 2) to facilitate further analysis. All data packages are publicly available on GitHub.

Ground-truth normalization

An expression profiling experiment of n genes on m samples typically results in an n × m matrix Y, the element y_ij of which captures the measurement of gene i in sample j. In RNA-seq, Y is the read count matrix output by the quantification process. To enable between-sample comparison, the read counts in Y requires proper normalization. The normalized expression matrix X is obtained by scaling Y against a set of m factors ν = {ν₁, ν₂, …, ν_m} such that every gene in sample j is scaled by the same factor ν_j. Let R be the set of reference genes whose expression levels are stable across all samples. If this set is identified, the scaling factor ν_j can be computed for each sample as follows

${\displaystyle {\nu }_j=\frac{\sum _{i\in R}^{}{y}_{ij}}{{\left(\prod _{k=1}^m\sum _{i\in R}^{}{y}_{ik}\right)}^{1/m}}.}$

Note that any multiplication of ν by a scalar results in a valid normalization. To obtain expression levels at a comparable magnitude with the raw counts, we chose to scale these factors by their geometric mean. This step facilitates the visual comparison of expression matrices but does not critically affect the outcome of differential expression analysis (Fig. S8).

In simulations, the reference set R is composed of all non-DE genes, i.e., those with fold-change set to one. In real data sets, R contains the ERCC spike-ins detected at sufficiently high number of read counts, that is, at least t reads in every sample. The threshold t ranges from 10 to 100 depending on the data sets, in order to eliminate unreliable measurements while retaining enough references to compute normalization factors. Manual examination of each data set was documented in the R notebook published as part of the source code.

Condition-number based deviation (cdev)

Definition. Assuming that a ground-truth expression matrix, denoted A, is known, the quality of the normalized expression matrix X can then be measured by the condition-number based deviation (cdev) between X and A, denoted cdev(X, A). To compute cdev(X, A), we first define the transformation from X to A. Let B be the matrix that transforms X to A, i.e., XB = A. As we are interested in recovering the pattern of differential expression rather than the absolute concentration of each RNA species, any multiplication of the ground truth A by a scalar is considered a valid normalization. Equivalently, X is considered valid if B = αI, where α is any scalar. Hence the quality of X can be measured by how much B deviates from the identity form. Such deviation is quantified by the condition number of B, defined as κ(B) =  ∥ B ∥ ⋅ ∥ B⁻¹ ∥. Choosing ∥• ∥ to be ℓ₂ norm, the condition number can be calculated using the ratio between the largest and smallest singular values, that is, $\kappa \left(\mathbf{B}\right)=\frac{{\sigma }_{max}\left(\mathbf{B}\right)}{{\sigma }_{min}\left(\mathbf{B}\right)}$. To compute B, one may recognize that B = (X^TX)⁻¹X^TA, in which (X^TX)⁻¹X^T is the Moore–Penrose pseudoinverse of X, denoted as X^†. Hence the definition of cdev can be written as follows

${\displaystyle cdev\left(\mathbf{X},\mathbf{A}\right)\equiv \kappa \left({\mathbf{X}}^{†}\mathbf{A}\right)=\frac{{\sigma }_{max}\left({\mathbf{X}}^{†}\mathbf{A}\right)}{{\sigma }_{min}\left({\mathbf{X}}^{†}\mathbf{A}\right)}.}$ The pseudoinverse X^† can be computed using the singular value decomposition X = UΣV^T, that is, X^† = VΣ⁻¹U (Alter & Golub, 2004; Alter et al., 2004).

Existence and uniqueness. The existence of cdev(X, A) depends on that of the pseudoinverse X^† and the singular value decomposition of X^†A, both of which are guaranteed to exist. Hence cdev(X, A) always exists for any expression matrices X and A. As the singular values of X^†A are unique, so is cdev(X, A).

Symmetry. One should note that, cdev is not symmetric, i.e., cdev(X, A) ≠ cdev(A, X) for generally given matrices X and A. However, when the transformation from X to A involves only column-scaling, we have cdev(X, A) = cdev(A, X). To see why, we first notice that the transformation matrix B, by which XB = A, is now diagonal,

${\displaystyle \mathbf{B}=\left[\begin{array}{cccc}\hfill b_{11}\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill b_{22}\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill b_{mm}\hfill \end{array}\right].}$ The transformation from A to X is then C = B⁻¹. By definition of cdev,

${\displaystyle cdev\left(\mathbf{A},\mathbf{X}\right)=\kappa \left(\mathbf{C}\right)=\kappa \left({\mathbf{B}}^{-1}\right)=\frac{1/{\sigma }_{min}\left(\mathbf{B}\right)}{1/{\sigma }_{max}\left(\mathbf{B}\right)}=\kappa \left(\mathbf{B}\right)=cdev\left(\mathbf{X},\mathbf{A}\right).}$ Consequently, when used to measure the changes in expression matrices made by global scaling, cdev(X, A) is a symmetric measure.

Interpretation. cdev ranges from 1 to infinity, with a smaller cdev indicating more similar expression patterns. The definition of cdev is closely related to that of a valid normalization. Specifically, any multiplication X′ of the matrix X by a scalar results in cdev(X, X′) = 1.

Because any non-empty subset of references result in a valid normalization, the result of scaling Y by the full set of references R, denoted X_R, should be similar to that of scaling by the subset S ⊂ R, denoted X_S. Hence, cdev(X_S, X_R) should be close to 1. This implication is helpful in validating the behavior of cdev on real data sets (See the Results section).

Implementation. Although cdev(X, A) can be computed for any X and A of the same shape using the definition above, it is much more trivial and efficient to compute in practice, when X is transformed to (or from) A by global scaling. In particular, the singular values of diagonal matrix B are now its elements. These elements also define the normalization vector v which can be re-constructed by dividing an arbitrary row X[i, ] of X by the corresponding row A[i, ] of A, as long as the feature i is not zero in any sample.

Normalization methods

As an illustration for the use of cdev, we compared six common normalization methods, including Total Count (TC), Upper Quartile (UQ) (Bullard et al., 2010), trimmed mean of M-values (TMM) (Anders & Huber, 2010), DESeq (Love, Huber & Anders, 2014), DEGES/TMM (Differentially Expressed Gene Elimination Strategy) (Kadota, Nishiyama & Shimizu, 2012), and PoissonSeq (Li et al., 2012). These methods were chosen due to the simplicity of computation and/or availability of R implementation which allows quick integration to the current analysis.

Each of these methods produce an m-vector v, the element v_j of which is used to scale sample j. While some methods produce close-to-1 scaling factors, some others require further normalization of these factors, for example, by the geometric mean of the elements v_j, to keep the magnitude comparable. As the magnitude of v does not affect the validity of its resulting normalization, a simplified description of v resulted by each method is provided below to highlight the differences between them.

In TC, the scaling factor v_j is proportional to the total count of all genes measured in sample j, i.e., ${\displaystyle v_j^{\mathrm{TC}}\propto \sum _{i=1}^n{y}_{ij}.}$

In UQ, the factor is computed by the upper quartile instead of the summation over gene counts, ${\displaystyle v_j^{UQ}\propto {\text{Upper-Quartile}}_{i\in G^{\ast }}{y}_{ij},}$

in which G^∗ is the set of genes with reads in at least one sample.

In TMM, the scaling factor is computed from the ratio of a gene count in sample j and the reference sample, arbitrarily selected as the first one, ${\displaystyle {\nu }_j^{\text{TMM}}\propto \frac{1}{\left|G^{\text{TMM}}\right|}\sum _{i\in G^{\text{TMM}}}\frac{y_{ij}}{y_{i1}},}$

in which G^TMM denotes the set of genes retained by TMM criteria.

In DESeq, scaling factors are computed from the medians of the ratio between a count in sample j and that in the pseudo-reference sample, which is computed by taking the geometric mean of each gene over the samples, ${\displaystyle {\nu }_j^{\text{DESeq}}\propto {\mathrm{median}}_i\left[\frac{y_{ij}}{{\left(\prod _{k=1}^m{y}_{ik}\right)}^{1/m}}\right].}$

DEGES is an iterative strategy in which the count matrix is scaled by a normalization method (such as TMM, DESeq, etc.), DE genes are identified by a DE caller (such as edgeR, DESeq, etc.), and the normalization is repeated with DE genes now removed (Kadota, Nishiyama & Shimizu, 2012; Sun et al., 2013). There are multiple choices for the normalizer and the DE caller, as well as the number of iterations. In this work, we chose TMM as the first-step normalizer, and edgeR as the second-step DE caller. With one iteration, the workflow is TMM-edgeR-TMM, which we denoted as DEGES/TMM. The scaling factor at iteration t is computed using the same method as that of TMM, while restricting the set of genes to that of G^t, ${\displaystyle {\nu }_j^{\text{DEGES/TMM}}\left(t\right)\propto \frac{1}{\left|G^{\text{t}}\right|}\sum _{i\in G^{\text{t}}}\frac{y_{ij}}{y_{i1}},}$

in which G^t is the set of genes identified as non-DE at iteration t.

PoissonSeq is also an iterative protocol which starts with fitting a Poisson log linear model under the null hypothesis that no gene is DE, then removes those with low ranking for goodness-of-fit, and fits again under the new hypothesis until convergence. The scaling factor is computed in the similar manner to that of TC, ${\displaystyle {\nu }_j^{\text{PoissonSeq}}\propto \frac{\sum _{i\in G^t}{y}_{ij}}{\sum _{k=1}^m\sum _{i\in G^t}{y}_{ik}},}$

in which G^t is the set of genes remained at iteration t.

Differential expression analysis

Differential expression (DE) analysis was used as a downstream indicator of normalization quality. Among a great number of methods developed for calling DE genes, we chose limma in consideration of its speed and robust performance relative to other methods (Soneson & Delorenzi, 2013). The method was implemented in the R package limma (Law et al., 2014; Ritchie et al., 2015).

Prior to DE calling, all the raw counts were increased by a small count (averaging 1), log-transformed, and normalized, using the function edgeR::cpm() from the edgeR package (Robinson, McCarthy & Smyth, 2010). The resulted log-cpm matrix was then passed to limma functions, following the limma-trend protocol (Law et al., 2014).

To verify that the magnitude of scaling factors does not affect the outcome of DE analysis regardless of the specific choice of DE caller, we performed DE gene detection using three methods limma-trend, DESeq (Love, Huber & Anders, 2014) and edgeR (Robinson, McCarthy & Smyth, 2010), as implemented in the R packages limma, DESeq2, and edgeR, respectively.

All the analyses were done in R computing environment (R Core Team, 2021).

Real data sets as experimental ground-truth

We presented here a diverse compilation of RNA-seq data sets with ERCC spike-ins, covering five different species: Homo sapiens, Rattus norvegicus, Mus musculus, Xenopus tropicalis, and Drosophila melanogaster. This collection is by far the most diverse set of publicly available RNA-seq assays with ERCC spike-ins.

Some earlier observations were reproduced on these data sets. First, the read counts of ERCC spike-ins are highly correlated with one another, especially those at high and moderate levels. This phenomenon was consistently observed in all data sets (Fig. S2), corroborating what was predicted theoretically and observed empirically on a mouse data set (Tran et al., 2020). Second, differences in library preparation protocol led to uneven efficiency among the spike-in sequences. Such protocol-dependent biases were observed earlier by Qing et al. (2013) (on data sets not available publicly) and by Owens et al. (2016) (on frog development time series re-processed into xtro1m and xtro2 in this work). Earlier studies such as Qing et al. (2013) and Risso et al. (2014a) also noted the fluctuation of reads mapped to ERCC spike-ins across all samples. Although we did find the same behavior, it should be noted that such fluctuation is an expected behavior of un-normalized read counts, and does not invalidate the use of these sequences as experimental references (Fig. S1).

Within each data set, the samples were spiked with the same concentration of ERCC spike-in RNAs and prepared with the same experimental protocol. Consequently, spike-ins are the true references in each data set. After normalizing by these references (see Methods section for the exact formula), the levels of ERCC spike-ins become well equalized (last columns in Fig. S3).

cdev is indicative of both upstream and downstream performance

In order for cdev to be a useful metric, it must be indicative of normalization quality which may be controlled by the quality of the reference set. It was shown that including differential genes in this set results in an invalid normalization (Tran et al., 2020). Using simulated data, it was clear that cdev is linearly correlated with the reference set quality, that is, inversely proportional to the fraction of differential genes included in the scaling set (Fig. 2, left column). This is true in both simulation scenarios, whether the fold-change is fixed for all DE genes (homogeneous fold-change), or varied over a continuous range (heterogeneous fold-change).

In line with conventional assessment, we expect cdev to be predictive of DE detection performance, for the improvement of normalization quality should lead to that of DE detection.

On the same simulated data as mentioned above, cdev is highly indicative of DE calling success which is summarized by the area under the receiver-operating characteristic (ROC) curve (AUC) (Fig. 2). When fold change is fixed for all DE genes, AUC as a function of cdev has a reverse sigmoid shape. The performance of DE calling drops significantly when cdev reaches a high enough threshold (Figs. 2A–2F). When fold change of DE genes spans a continuous range, this relation gets closer to a linear form (Figs. 2G–2J).

cdev is agnostic to the magnitude of scaling factors, as are DE gene calling procedures. This behavior is demonstrated empirically on simulated data sets. DE performance is almost identical on expression matrices scaled by the same-direction vectors spanning several orders of magnitude. Not only limma, but also count-based methods such as DESeq and edgeR, were not affected by the magnitude of scaling factors (Fig. S8).

Simplistic simulation enables us to tear apart the effects of different data characteristics on cdev. By varying one parameter at a time, we found that the range of cdev is not affected by the number of genes, nor the fraction of DE genes, but determined by the level of difference among the samples, represented by the fold change between two simulated conditions (Figs. 2A, 2C, 2E). Fold-change also dictates how tolerant DE detection can be toward bad normalization. Specifically, when fold change is as high as 8, a normalized expression matrix with cdev of 3 can still result in an AUC larger than 0.5. At smaller fold change such as 2, this level of deviation from the ground truth is sufficient to drop the AUC of DE detection to zero. This behavior makes sense, for it is easier to call DE when the fold change is large.

Since the behavior of cdev on real data sets could not be validated with the same metrics as that on simulation data, we examined if it behaves as theoretically predicted. Specifically, cdev(A′, A) should be close to 1 with A being normalized by the true reference set R, and A′ by any non-empty subset S of R. Meanwhile, cdev(A_rand, A) where A_rand is resulted from a random set of genes should be larger. These expectations were observed in all the real data sets (Fig. 3). In highly similar samples, such as those from brain regions of mice with identical genetic background (data sets mpreg1a, mpreg1b, mpreg2a, mpreg2b), cdev(A_rand, A) is larger than cdev(A′, A), but never exceeds 10. In highly different samples, such as those from 16 fly strains (dmelAC, dmelAV), random normalization results in very large deviation from the ground truth, with cdev(A_rand, A) mostly around 100. Data sets with average heterogeneity, such as those of rat tissues (rbm1, rbm2) or sarcoma samples (sarcoma), exhibit cdev(A_rand, A) in the middle range.

These observations illustrate the value range of cdev in practice, at the same time implies the importance of normalization with respect to the heterogeneity of transcriptomic data. Specifically, when the expression profiles are highly similar, wrong normalization is not as detrimental as it would be on heterogeneous ones. Based on the behavior of cdev on real data sets (Fig. 3), we may recommend acceptable values of cdev below which DE performance is better than random. Specifically, on extremely diverse expression profiles (dmelAC, dmelAV), cdev up to 10 may be acceptable for DE analysis. On moderately diverse expression profiles (rbm1, rbm2, rnor1, rnor2, sarcoma), cdev below 5 may be necessary to obtain good enough DE calling. On closely related expression profiles (mpreg1b, mpreg2b, yfv), cdev should be as low as 2 in order to provide a good enough DE results. As cdev is a continuous measure of normalization quality, these values should not be treated as hard threshold, but rather a guideline to understand the normalization outcome. The intent of cdev is not to predict downstream analysis performance, but rather to enable the measurement of normalization quality independent from downstream analyses, as would be argued in the Discussion section.

cdev as a robust alternative to fold-change MSE

Among a few ground-truth based measures of normalization quality, fold-change MSE is an intuitive metric. As the name implies, a fold-change MSE closer to zero indicates a better normalization. In Maza et al. (2013), fold-change MSE was defined for all genes, while in Evans, Hardin & Stoebel (2017), it was narrowed down to the true reference (i.e., non-DE) genes. The latter version requires less information from the ground truth and works equally well as the former, hence chosen to compare with cdev. When defined for true references, the expected fold-changes are all one.

In order to compare the two metrics on simulation data, we normalized an expression matrix by random sets of genes, computed both cdev and fold-change MSE with respect to the corresponding ground truth. To ensure these randomly normalized matrices cover a wide enough range of quality, we varied the fraction of differential genes allowed in the set of scaling genes. On this simulated data, the two metrics are highly correlated (Fig. 4).

On real data, the relationship becomes more noisy (Fig. 5). On a closer look, the “noise” is actually caused by multiple slopes of the linear relation between cdev and fold-change MSE. This complicated behavior happens due to the ambiguity of fold-change MSE by which many improper normalizations are able to score zero fold-change error (Fig. 6E). The problem is exaggerated on multi-factorial experiments, such as that of the lymphoma data set. In this experiment, different cell lines originating from the parental cells were treated with THZ1 (Zhao et al., 2019), resulting in samples of varied cell lines and THZ1 concentrations. When the data set is composed of samples from a single cell line, such as I10 or L20, the relationship between cdev and fold-change MSE is cleanly linear. When the samples include those from the two cell lines, the two linear functions with different slopes aggregate and make the correlation seem noisy (Fig. 5). Discordance between the two metrics is caused by large cdev and small MSE. One such example was visualized in Fig. 6. In this case, the random set of genes could not equalize the true references, nor themselves, resulting in a poorly normalized expression matrix. This faulty result was able to score a low fold-change MSE, but not a low cdev. Projections on the first two principal components of the ERCC spike-ins showed that random references could not remove the variation between two organs (heart and kidney), while the true ones properly removed such variation. Although there may seem to be a few instances with “high” fold-change MSE and “low” cdev, they are in fact part of a cdev-MSE function with a lower slope, thus do not contribute to the discordance between the two metrics. Inspection of such cases confirm that these random normalizations are clearly worse than the true one (Fig. S7).

Application of cdev to compare common normalization methods

cdev is specifically devised to measure performance by global scaling normalization. To illustrate its use, we compare six common normalizers. Four of these are one-step methods in which the scaling factors are computed using a closed form formula. The other two (DEGES/TMM and PoissonSeq) are iterative procedures in which the set of non-differentially expressed features are refined until convergence. In all the following evaluation, the ground truth is computed by scaling against the set of known references. For simulated data, this set includes the genes with fold-change set to 1 by the simulator, such that the number of DE (and non-DE) genes satisfies the required percentage of differential expression. For real data, it includes selected ERCC spike-ins which have sufficiently high number of reads. In both simulation and real data, the true references were able to equalize themselves across all samples (Figs. 7 and 8). All normalizers perform worse when the input data becomes more challenging, i.e., contains a larger fraction of DE genes. At the highest DE fraction simulated, none of the normalizers could stabilize the expression levels of non-DE (i.e., reference) genes (Fig. 7D). As the quality of normalization decreases, so is the performance of DE detection. The false positive rate (FPR) and AUC of DE calling similarly decreases with increasing DE percentage, demonstrating a linear relationship with cdev (Fig. 7E). Although false negative rate (FNR), i.e., miss rate, was used in previous studies (Maza et al., 2013), we found it not informative for DE detection performance (Fig. S5). Apparently, one may achieve zero FNR by calling all genes differential. Meanwhile, incorrect normalization is more likely to increase than reduce the variation, hence increasing the chance of false discovery, but not that of missing the true differential genes. Based on the high concordance between cdev and fold-change MSE on simulation data (Fig. 4), one may expect that similar results would be observed with fold-change MSE. Our analysis indeed confirmed that (Fig. S6).

Real data sets are more challenging, for none of the normalizers could stabilize the true references properly (Fig. S3). Figure 8 shows such results on the Rat Bodymap data sets as an example. In this case, the cdev between the raw count matrices and their corresponding ground truths are even smaller than those produced by the normalizers, suggesting that normalization may induce more variation, driving the expression profiles to deviate further from the true patterns. These results imply that on such data sets, no normalization may be better than doing with the wrong assumption.