SC-JNMF: single-cell clustering integrating multiple quantification methods based on joint non-negative matrix factorization

The outline of SC-JNMF

We proposed SC-JNMF, a method that extracts latent factors from different gene expression profiles at the same time using joint-NMF, which can express matrix data as the product of lower-dimensional matrices, one of which is shared. SC-JNMF extracted the latent factors in different gene expression profiles using a similar approach to NMF and used them for cell clustering and gene analysis.

The outline of our method is showed in Fig. 1. We created two different gene expression profiles using different quantification methods from the same RNA-seq library. Then, we extracted the common factor using joint-NMF and extended it to perform multiple NMFs in parallel with two different basis factors W₁, W₂ (derived from the different methods) and shared factors H (derived from the original RNA-seq library). Finally, we performed cell clustering using these extracted factors and any appropriate clustering methods, such as hierarchical clustering. Although the conventional NMF clustering methods that perform clustering directly with the factorized matrix strongly depend on the rank, our method performs robust clustering by using hierarchical clustering.

Quantification and normalization

As the first step of SC-JNMF, we quantified gene expressions with different quantification pipelines, and obtained a set of gene expression profiles. Parameters and references could be set as recommended by each quantification pipeline.

For preprocessing before performing joint matrix factorization, we applied a gene filter to the input data, similar to the SC3 method (Kiselev et al., 2017), and then log₂(x + 1) transformed the data. Finally, we normalized the data so that the sum of each gene L¹ norm was 1 because it should be compared the relative expression values between cells for each gene to perform the cell clustering. This L₁ normalization is possible to consider smaller expression values than using the L₂ norm and to keep the values non-negative.

Joint-NMF

We considered the given data as a non-negative matrix $\mathit{D}\in {\mathbb{R}}_{+}^{N\times M}$. Non-negative matrix factorization found the basis matrix $\mathit{W}\in {\mathbb{R}}_{+}^{N\times k}$ and coefficient matrix $\mathit{H}\in {\mathbb{R}}_{+}^{k\times M}$, where all the elements were non-negative, such that these matrix products approximated the input data matrix. (1)${\displaystyle \mathit{D}\approx \mathit{W}\mathit{H}}$ NMF minimized the distance between matrix D and matrix WH. Here, we considered a Euclidean norm distance. Thus, the objective function that NMF minimizes was as follows: (2)${\displaystyle \underset{\mathit{W},\mathit{H}\ge 0}{min}L:={\parallel \mathit{D}-\mathit{W}\mathit{H}\parallel }_{\mathcal{F}}^2}$ Here, we considered a simultaneous matrix factorization of two matrices using NMF. Given the two input data as non-negative matrices ${\mathit{D}}_1\in {\mathbb{R}}_{+}^{N_1\times M}$ and ${\mathit{D}}_2\in {\mathbb{R}}_{+}^{N_2\times M}$, joint-NMF found the basis matrices ${\mathit{W}}_1\in {\mathbb{R}}_{+}^{N_1\times k},{\mathit{W}}_2\in {\mathbb{R}}_{+}^{N_2\times k}$ corresponding to each approximated input matrix and the common coefficient matrix $\mathit{H}\in {\mathbb{R}}_{+}^{k\times M}$, minimizing the distances between the given matrices and the approximated matrices. In our proposed method, we added the L1 norm constraint to this objective function for H so that the factorized shared matrix was sparse. Thus, SC-JNMF found the matrix W₁, W₂, H that minimized the following objective function: (3)${\displaystyle \underset{{\mathit{W}}_1,{\mathit{W}}_2,\mathit{H},{\lambda }_n\ge 0}{min}L:={\parallel {\mathit{D}}_1-{\mathit{W}}_1\mathit{H}\parallel }_{\mathcal{F}}^2+{\lambda }_1{\parallel {\mathit{D}}_2-{\mathit{W}}_2\mathit{H}\parallel }_{\mathcal{F}}^2+{\lambda }_2\sum _k\parallel {\mathit{H}}^{k,\ast }{\parallel }_1}$ where, λ₁ is a balance parameter for the losses of reconstruction matrices and λ₂ is a parameter for row vector sparsity regularization of shared factorized matrix.

We applied a multiplicative update algorithm to optimize the objective function same as conventional NMF. By applying Jensen’s inequalities to the first and second terms of the objective function, the function to be minimized could be rewritten as follows: (4)${\displaystyle \underset{{\mathit{W}}_1,{\mathit{W}}_2,\mathit{H},{\lambda }_n\ge 0}{min}L:=\sum _{i,j}\left(|{\mathit{D}}_1^{i,j}{|}^2-2{\mathit{D}}_1^{i,j}\sum _k{\mathit{W}}_1^{i,k}{\mathit{H}}^{k,j}\right.\left.+\sum _k\frac{{\left({\mathit{W}}_1^{i,k}\right)}^2{\left({\mathit{H}}^{k,j}\right)}^2}{c_1^{i,j,k}}\right)+{\lambda }_1\sum _{i,j}\left(|{\mathit{D}}_2^{i,j}{|}^2-2{\mathit{D}}_2^{i,j}\sum _k{\mathit{W}}_2^{i,k}{\mathit{H}}^{k,j}\right.\left.+\sum _k\frac{{\left({\mathit{W}}_2^{i,k}\right)}^2{\left({\mathit{H}}^{k,j}\right)}^2}{c_2^{i,j,k}}\right)+{\lambda }_2\parallel \mathit{H}{\parallel }_1}$ where, (5)${\displaystyle c_1^{i,j,k}=\frac{{\mathit{W}}_1^{i,k}{\mathit{H}}^{k,j}}{\sum _{k^{\prime }}{\mathit{W}}_1^{i,k^{\prime }}{\mathit{H}}^{k^{\prime },j}}}$ (6)${\displaystyle c_2^{i,j,k}=\frac{{\mathit{W}}_2^{i,k}{\mathit{H}}^{k,j}}{\sum _{k^{\prime }}{\mathit{W}}_2^{i,k^{\prime }}{\mathit{H}}^{k^{\prime },j}}}$ We found each element of W₁, W₂, and H that minimized the objective function by performing partial differentiation. (7)${\displaystyle \frac{\partial L}{\partial W_1^{i,k}}=\sum _j\left(-2D_1^{i,j}{H}^{k,j}+\frac{2W_1^{i,k}{\left(H^{k,j}\right)}^2}{c_1^{i,j,k}}\right)}$ (8)${\displaystyle \frac{\partial L}{\partial W_2^{i,k}}={\lambda }_1\sum _j\left(-2D_2^{i,j}{H}^{k,j}+\frac{2W_2^{i,k}{\left(H^{k,j}\right)}^2}{c_2^{i,j,k}}\right)}$ (9)${\displaystyle \frac{\partial L}{\partial H^{k,j}}=\sum _i\left(-2D_1^{i,j}{W}_1^{i,k}+\frac{2{\left(W_1^{i,k}\right)}^2{H}^{k,j}}{c_1^{i,j,k}}\right)+{\lambda }_1\sum _{i^{\prime }}\left(-2D_2^{i^{\prime },j}{W}_2^{i^{\prime },k}+\frac{2{\left(W_2^{i^{\prime },k}\right)}^2{H}^{k,j}}{c_2^{i^{\prime },j,k}}\right)+{\lambda }_2}$

The objective function are minimized when these are 0. Thus, the variable updates became: (10)${\displaystyle {\mathit{W}}_1={\mathit{W}}_1\frac{{\mathit{D}}_1{\mathit{H}}^{\top }}{{\left[\mathit{H}{\left[{\mathit{W}}_1\mathit{H}\right]}^{\top }\right]}^{\top }}}$ (11)${\displaystyle {\mathit{W}}_2={\mathit{W}}_2\frac{{\lambda }_1{\mathit{D}}_2{\mathit{H}}^{\top }}{{\lambda }_1{\left[\mathit{H}{\left[{\mathit{W}}_2\mathit{H}\right]}^{\top }\right]}^{\top }}}$ (12)${\displaystyle \mathit{H}=\mathit{H}\frac{{\left[{\mathit{D}}_1^{\top }{\mathit{W}}_1\right]}^{\top }+{\lambda }_1{\left[{\mathit{D}}_2^{\top }{\mathit{W}}_2\right]}^{\top }-{\lambda }_2/2}{{\mathit{W}}_1^{\top }{\mathit{W}}_1\mathit{H}+{\lambda }_1{\mathit{W}}_2^{\top }{\mathit{W}}_2\mathit{H}}}$

Applications

Dataset

We used six different datasets including RNA-seq library and quantified gene expression levels measured using each method in previous studies (Table 1). Only Monaco dataset is a bulk RNA-seq, the others are scRNA-seq libraries. Because Monaco dataset is a set of sorted cell types, it is suitable for the evaluation of clustering accuracy. CellBench is a set of datasets designed to benchmark various single-cell data analysis. In this study, we used the dataset named “sc_10x_5cl”, single cells from the mixture of five cell lines. In advance, we removed any cell types that were unclear in previous studies.

To perform our method SC-JNMF and compare the accuracy of the clustering methods, We alternatively quantified gene expression values in each cell from the RNA-seq library using Salmon(v1.0.0) with GENCODE references and annotations (Mouse Release 21/ Human Release 32). Only CellBench dataset was quantified by Using Salmon with Alevin (Srivastava et al., 2019), kallisto with BUStools and STAR with Cell Ranger to process barcodes in short reads.

In addition, for the first 4 datasets, we also prepared other gene expression profiles using kallisto(v0.46.2), STAR(v2.7.3), to evaluate our method for the combinations of quantification methods and the sample size.

Settings of clustering methods

We compared the accuracy of cell clustering using our proposed method to that obtained using other major unsupervised clustering methods, including LSNMF (Lin, 2007; Zitnik & Zupan, 2012), SC3 (Kiselev et al., 2017), and Seurat (Satija et al., 2015).

For our proposed method, each run consisted of ten runs and extracted an ARI score of the top combined rank in terms of sparseness and loss (the top ranking of sparseness is maximum and the top ranking of loss is minimum). The rank was determined according to the results of our experiment described in the Supplemental Information and Fig. S1 as follows: Treutlein, k = 5; Pollen, k = 8; Segerstolpe, k = 25; Xin, k = 16; Monaco, k = 8; CellBench, k = 5. In addition, we determined the regularization parameter λ₁ = |geneset1|/|geneset2|. Hierarchical clustering could not determine the number of clusters; therefore, we set the same number reported in previous studies in advance.

In the LSNMF method, similar to our joint-NMF, classification by hierarchical clustering was performed using a matrix of factors of the cells, and the rank and the number of clusters were set to the same values as ours. SC3 method was performed with the default parameters. For LSNMF and SC3, each run was performed ten times, considering random initial values. For Seurat, we plotted only the highest ARI score in the runs as the resolution parameter in the FindClusters function was increased from 0 to 1 in 0.1 steps. The parameters of the FindNeighbors function were set to the default values (dims = 1:10, k.param = 20).

Accuracy of cell clustering

Figure 2 shows the ARI score of the original (quantified previously) and alternative (quantified using Salmon) data for each clustering method. We ran a two-way experiment using λ₂ = 1 and λ₂ = 10.

In the Pollen, Xin, and CellBench datasets, our proposed method performed accurate cell clustering. In contrast, in the Treutlein, Segerstolpe, and Monaco datasets, the ARI score of the proposed method was higher than that of LSNMF, but it was not the highest score. Additionally, the resulting clusters of our method were stable compared to the conventional NMF method (except for the CellBench dataset).

Comparison of quantification methods in SC-JNMF

SC-JNMF performed cell clustering using one or two gene expression profiles. In this study, we compared the combinations of quantification methods (kallisto, STAR, Salmon) by performing cell clustering and calculating ARI using SC-JNMF. To evaluate the impact of “joint” NMF, we also performed factorization and clustering using a single gene expression profile with setting to the parameter λ₁ = 0. We ran SC-JNMF repeatedly for ten trials with random initial values and the parameter λ₂ = 10.

Figure 3 shows the ARI of the cell clustering using SC-JNMF in each quantification method and their combination. As a result, for each dataset, the accuracies of kallisto & STAR and Salmon & STAR combination outperformed using the single quantification method. The accuracy of Salmon & kallisto combination was lower than that of other combinations.

Additionally, we counted the frequency of correlation coefficients for each cell and each gene in each dataset and combination (Fig. 4). In Salmon & kallisto combination, the correlation coefficients were inclined toward a higher value than the other combinations. These results indicated that SC-JNMF had a lower accuracy when using similar gene expression profiles (e.g., Salmon & kallisto combination). Therefore, two gene expression profiles with different properties, having a lower correlation of cells and genes, were suitable for SC-JNMF.

Size effect for SC-JNMF and NMF

To further understand the benefit of NMF-based methods, the effect of sample size (the number of cells in the RNA-seq library) against clustering accuracy was evaluated. We randomly subsampled the 25%, 50%, and 75% cells from gene expression profiles, and ran SC-JNMF (λ₂ = 10) and LSNMF. Each trial was performed ten times.

Figure 5 shows the ARI of the cell clustering using LSNMF and SC-JNMF for each dataset and subsample rate. As a result, for each dataset, the accuracies of both LSNMF and SC-JNMF were increased depends on the subsample rate. This result indicated that the methods based on NMF benefit more from a dataset including a larger number of cells.

Gene analysis using factorized matrix

We found marker genes of the Xin dataset using factorized matrices. The Xin dataset contained data on 1492 single cells and four classes (alpha, beta, delta, and PP) in the pancreas. We showed the results in Fig. 6. The factor of the coefficient matrix showed some characteristic patterns in each cell cluster (e.g., Factor 2 and Factor 5 showed higher values in delta and PP cells) (Fig. 6A). Next, we calculated the correlation of factors between the basis matrices in the common genes (Fig. 6B). The factors in the basis matrices showed a similar tendency. We also showed the loadings of marked factors for delta cells in the coefficient matrix in Fig. 6C. Almost all genes showed similar factor loadings between base1 and base2; however, we confirmed that some of the genes only observed in either gene expression profile also had high values. We also showed the marker genes of delta cells detected using Scanpy Wolf, Angerer & Theis (2018) in the scatter plot. The factor loadings of these genes tended to be higher than those of others, regardless of whether the gene was observed in both gene expression profiles or not. All factors of Xin dataset were shown in Fig. S2.