A robust semi-supervised NMF model for single cell RNA-seq data

Related work

Nonnegative matrix factorization (NMF) is an effective method for unsupervised clustering analysis of gene expression data. Given a nonnegative matrix X of size m × n, NMF aims to find two non-negative matrices W and H such that: ${\displaystyle X\approx WH}$

Where W ∈ R^m×k is a basis matrix, H ∈ R^k×n is a coefficient matrix. The solution to this problem can be obtained by solving the following optimization problem: ${\displaystyle {min}_{W,H}O\left(W,H\right)=mi{n}_{W,H}\frac{1}{2}\left|\right|X-WH|{|}_F^2}$

||.||_F means Frobenius norm and W and H satisfy nonnegative constrains. Since many practical problems in data mining, pattern recognition and machine learning require non-negativity constraints and appropriate low dimensional representation of original data, NMF has been successfully applied to these fields, which obtains the parts-based representation as well as enhancing the interpretability of the data (Berry et al., 2007; Chagoyen et al., 2006; Lee & Seung, 1999). In bioinformatics, it has been used to extract meaningful features that belongs to different cell types from microarray and scRNA-seq data and identify mRNA isoforms (Brunet et al., 2004; Shao & Hofer, 2017; Ye & Li, 2016).

Different variants of NMF has been put forward such as sparse NMF (SNMF) and discriminant NMF (DNMF) for microarray and RNA-seq data (Jia et al., 2015; Kim & Park, 2007). SNMF introduces a regularization term on W or H to control the degree of sparsity and generate sparser representation. DNMF incorporates Fisher’s discriminant criterion in the coefficient matrix by maximizing the distance among any samples from different classes meanwhile minimizing the dispersion between any pair of samples in the same class. The DNMF requires discriminant information to construct the objective function and has been applied in various scenes such as face recognition and facial expression recognition. Yet in most cases of scRNA-seq, we cannot know the exact class information. Moreover, none of these methods consider the technical factors including amplification, library size differences and dropouts (Buettner et al., 2015; Kharchenko, Silberstein & Scadden, 2014). In addition to the technical factors, scRNA-seq data exhibit high cell-to-cell variability in gene expression, which impeding the analysis. Overall, the technical effects and inherent variability in scRNA-seq introduce substantial noise and may corrupt the analysis of underlying biological process.

In this paper, we propose a robust semi-supervised NMF model that is robust to noises and uses marker information of cells as prior knowledge. We already have prior knowledge about the marker genes of different cells (such as MSigDB database; Liberzon et al., 2011), for instance, metastatic melanoma is a mixture of tumor cells and a variety of normal cell including T-cells, B-cells, macrophages and NK cells; thus, we can know the prior marker genes of different groups. We will show that incorporating such information into the NMF model, the clustering accuracy can be increased significantly.

Original NMF algorithm

The original NMF algorithm was introduced by Lee and Seung and the objective function is non-increasing if follows the multiplicative update rules (Lee & Seung, 2001): ${\displaystyle W_{ij}⟵W_{ij}\frac{{\left(X{H}^T\right)}_{ij}}{{\left(WH{H}^T\right)}_{ij}}{H}_{ij}⟵H_{ij}\frac{{\left(W^{T}X\right)}_{ij}}{{\left(W^{T}WH\right)}_{ij}}}$

The author proved that by repeating iteration of the update rules is guaranteed to converge to a locally optimal matrix factorization. By using the non-negative constrains, NMF can learn a low dimensional representation of the original data.

Roubst NMF algorithm

Robust NMF (rNMF) first proposed by Kong and Ding and was designed to handle outliers and noises that original NMF fail to cope with (Kong, Ding & Huang, 2011). The main difference of rNMF is using L_2,1 norm loss function instead of Frobenius norm which is defined as follow: ${\displaystyle \left|\right|X-WH|{|}_{2,1}=\sum _{j=1}^n\sqrt{\sum _{h=1}^m{\left(X-WH\right)}_{hj}^2}=\sum _{j=1}^n||x_j-W{h}_j||}$

In this formulation, the error of each data point is ||x_j − Wh_j|| rather than squared, thus the errors caused by outliers and noises do not dominate the objective function compared with the squared one. Different from the revision on the objective function, Zhang proposed a method which introduce an error matrix S ∈ ℝ^n×m to capture the errors and handle the extreme data points (Zhang et al., 2011): ${\displaystyle {min}_{W,H,S}\left|\right|X-WH-S|{|}_F^{2}s.t.W\ge 0,H\ge 0,|\left|S\right|{|}_0\le \upsilon }$

Where ν is the parameter that specifies the maximum number of nonzero elements in S.

Robust semi-supervised NMF algorithm

Motivated by robust NMF formulation, we adopt similar method to deal with the noises in scRNA-seq data under the assumption that the data matrix may be corrupted by noises and the noises are sparse (we use rNMF to represent this model below). In addition, we incorporate the prior information into a graph regularization term to maintain the intrinsic geometrical and discriminating structure of the data in the parts-based representation: ${\displaystyle {min}_{W,H,S}\left|\right|X-WH-S|{|}_F^2+\alpha |\left|S\right|{|}_1+\beta Tr\left(HL{H}^T\right)s.t.W\ge 0,H\ge 0}$

Where Tr(.) denotes the trace of a matrix and α, β is the regularization parameters. Since the 𝔩₀ − norm is difficult to solve, we replace the 𝔩₀ − norm constraint with a 𝔩₁ − norm regularizer, which is a common technique for sparse solution. The second graph regularization term maintains the local geometrical structure of the original data matrix X in the low dimensional representation, i.e., the matrix H. Briefly, if two data points x_i and x_j are close to each other in matrix X, then h_i and h_j, points in the low dimensional representation matrix H, are also close to each other. This property is usually referred to as local invariance assumption (Belkin & Niyogi, 2002; Cai, Wang & He, 2009) and has been applied in the development of dimensionality reduction algorithms and semi-supervised learning algorithms (Belkin & Niyogi, 2002; Zhu & Lafferty, 2005).

Chung et al. demonstrate that the local geometrical structure can be modeled by a nearest neighbor graph on a scatter of data points (Chung & Graham, 1997). For each data point x_i in X, we need to find its neighbors and put edges between x_i and its neighbors. In our scRNA-seq problem, we suppose that all the cells are the nodes in a graph, but different group of cells have different weighting. We use marker genes of cell group to construct the weight matrix Q, which can be defined by various ways. In this paper, we chose the heat kernel weighting and the weighting between x_i and x_j are: ${\displaystyle Q_{ij}=e^{-\frac{\left|\right|x_i-x_j|{|}^2}{\sigma }}}$

where σ is parameter that controls the weighting and we set σ = 1 in the experiment. Note that only the rows (genes in rows and samples in columns for matrix X) that we select as cell markers of X will be used to compute the weighting. As shown in Fig. 1, the expression matrix X includes five groups and each group has its own highly-expressed genes (yellow rectangular block). The rows which are the markers of different groups are selected to construct the weight matrix.

We use Euclidean distance to measure the dissimilarity in the low dimensional representation of X: ${\displaystyle D\left(h_i,h_j\right)=\left|\right|h_i-h_j|{|}^2}$

then the local geometrical structure preserving criterion can be written as: ${\displaystyle \Re =\sum _{\mathit{i},\mathit{j}=1}^{\mathit{n}}\left|\right|h_i-h_j|{|}^2{Q}_{ij}=\sum _{\mathit{i}=1}^{\mathit{n}}{h}_i^T{h}_i{D}_{ii}-\sum _{\mathit{i},\mathit{j}=1}^{\mathit{n}}{h}_i^T{h}_j{Q}_{ij}=Tr\left(HD{H}^T\right)-Tr\left(HQ{H}^T\right)=Tr\left(HL{H}^T\right)}$

Where D is a diagonal matrix and D_ii = ∑_jQ_ij. L is called graph Laplacian and L = D − Q. We can see that by minimizing ℜ, if data point x_i and x_j are close (Q_ij is big), h_i and h_j will also be close to each other, so the distance relation between points in original data matrix X can be preserved in low dimensional matrix H. This technique was first designed by Cai et al. and introduced to solve the embedded structure problem, which assume that the data is usually sampled from a low dimensional manifold embedded in a high dimensional ambient space (Cai et al., 2011). We show that by encoding the geometrical information in the NMF model, the clustering accuracy can be greatly improved.

Optimization algorithm for robust semi-supervised NMF

In this section, we derive the iterative multiplicative update rules for the robust semi-supervised NMF based on the coordinate descent method (Sha et al., 2007). The objective function is not convex for W, H and S together, so it is unrealistic to expect an algorithm to find the global minima. It is natural to optimize W, H and S separately since the objective function is convex while holding the other two variables as constant.

The objective function of rssNMF can be rewritten as: ${\displaystyle \mathfrak{O}=\left|\right|X-WH-S|{|}_F^2+\alpha |\left|S\right|{|}_1+\beta Tr\left(HL{H}^T\right)=Tr\left(\left(X-WH-S\right){\left(X-WH-S\right)}^T\right)+\beta Tr\left(HL{H}^T\right)+\alpha |\left|S\right|{|}_1=Tr\left(X{X}^T\right)-2Tr\left(X{H}^T{W}^T\right)-2Tr\left(X{S}^T\right)+2Tr\left(WH{S}^T\right)+Tr\left(WH{H}^T{W}^T\right)+Tr\left(S{S}^T\right)+\beta Tr\left(HL{H}^T\right)+\alpha |\left|S\right|{|}_1}$

After introducing the Lagrange multiplier ΨandΦ for the constrains, the Lagrange function 𝔏 is stated as ${\displaystyle \mathfrak{L}=Tr\left(X{X}^T\right)-2Tr\left(X{H}^T{W}^T\right)-2Tr\left(X{S}^T\right)+2Tr\left(WH{S}^T\right)+Tr\left(WH{H}^T{W}^T\right)+Tr\left(S{S}^T\right)+\beta Tr\left(HL{H}^T\right)+\alpha |\left|S\right|{|}_1+Tr\left(\Psi W\right)+Tr\left(\Phi H\right)}$

and the partial derivatives of 𝔏 with respect to WandH are: ${\displaystyle \frac{\partial \mathfrak{L}}{\partial \mathrm{W}}=-2X{H}^T+2S{H}^T+2WH{H}^T+\Psi \frac{\partial \mathfrak{L}}{\partial \mathrm{H}}=-2W^{T}X+2W^{T}S+2W^{T}WH+2\beta HD-2\beta HQ+\Phi }$

By KKT conditions ψ_hiw_hi = 0andϕ_ijh_ij = 0, we can obtain: ${\displaystyle -{\left(X{H}^T\right)}_{ij}{w}_{ij}+{\left(S{H}^T\right)}_{ij}{w}_{ij}+{\left(WH{H}^T\right)}_{ij}{w}_{ij}=0-{\left(W^{T}X\right)}_{ij}{h}_{ij}+{\left(W^{T}S\right)}_{ij}{h}_{ij}+{\left(W^{T}WH\right)}_{ij}{h}_{ij}+{\left(\beta HD\right)}_{ij}{h}_{ij}-{\left(\beta HQ\right)}_{ij}{h}_{ij}=0}$

Then the updating rules for W and H are: (1)${\displaystyle W_{ij}⟵W_{ij}\frac{{\left(X{H}^T\right)}_{ij}}{{\left(WH{H}^T+S{H}^T\right)}_{ij}}}$ (2)${\displaystyle H_{ij}⟵H_{ij}\frac{{\left(W^{T}X+\beta HQ\right)}_{ij}}{{\left(W^{T}WH+W^{T}S+\beta HD\right)}_{ij}}}$

For fix W and H, we update S via the soft-thresholding operator (Hale, Yin & Zhang, 2008) and the optimization problem for S is: ${\displaystyle {min}_S\left|\right|X-WH-S|{|}_F^2+\alpha |\left|S\right|{|}_1}$

Theorem 1: Define the soft-thresholding operator as below: ${\displaystyle T_{\upsilon }\left(z\right)=\left\{\begin{array}{c}z-\upsilon ,ifz>\upsilon \hfill \\ z+\upsilon ,ifz<-\upsilon \hfill \\ 0,otherwise\hfill \end{array}\right.}$

Where z ∈ ℝ and υ > 0 and this operator can also be applied to vectors or matrices by element-wise operation. For the following 𝔩₁ − norm problem: ${\displaystyle {min}_v\frac{1}{2}\left|\right|x-v|{|}_F^2+\alpha |\left|v\right|{|}_1}$

the unique solution of v is given by T_α(x) (Hale, Yin & Zhang, 2008). Similarly, we can get the update rule for S: (3)${\displaystyle \mathrm{S}⟵T_{\frac{\alpha }{2}}\left(X-WH\right)}$

Based on above analysis, our algorithm for solving rssNMF is presented in Algorithm 1.

The proof of convergence analysis of our algorithm essentially follows the idea of graph NMF and ensures that the objective function of rssNMF is nonincreasing under the updating rules in Eqs. (1), (2) and (3) and could converge to a stationary point (Cai et al., 2010). Since the idea is similar to the method in the paper (Cai et al., 2010) so we did not prove it here.

Dataset

Cell markers, usually referring to surface molecules on cell membrane, are different in different kinds of cells. Surface molecules that appear only on a particular type of cell are called cell markers and often used for cell type identification. For example, cancer cell specific markers often appear on the surface of cancer cells and are utilized as targets for anticancer drugs. Cell markers have been a hotspot in molecular biology study and researchers have accumulated lots of marker data for different kinds of cells, such as MSigDB (Liberzon et al., 2011). Here, we extend the concept of cell markers. We use ”cell marker genes” to represent those genes that can distinguish different types of cell, which are usually only expressed in specific type of cells or relatively highly expressed, such as differentially expressed genes. In single-cell RNA-sequencing, we already know where the sequencing samples come from and then we can have a modest degree of prior knowledge. That is, we have known what kinds of cells are in the experimental samples although there are still many kinds of cells remaining unknown. For instance, when the cancer tissue was sequenced, the samples contained a large number of normal cells in addition to cancer cells, such as lymphocytes, myeloid populations and cancer-associated fibroblasts. Taking melanoma as an example, melanoma is often known as malignant melanoma, which is a kind of cancer that stems from the pigment-containing cells known as melanocytes. Melanoma is a mixture of tumor cells and a variety of normal cells including T-cells, B-cells, macrophages and NK cells. This feature has been found in many studies, thus we can know the prior marker genes of different groups of cells in advance. We will show that incorporating such information into the NMF model, the clustering accuracy can be increased significantly.

Twelve scRNA-seq datasets are used in the experiment. All the datasets are normalized expression level (FPKM, RPKM or CPM) or counts and number of cells range from 56 to 3005 (Table 1). Almost allthe datasets except dataset “Zeisel” have true labels since they are collected from different time points during embryonic development and we provide the accession number in Table 1. The labels of dataset “Zeisel” comes from computation and can be seen as a silver dataset.

Compared algorithms

To benchmark rssNMF, six clustering algorithms are used: NMF, rNMF, ssNMF, KMeans, Hierarchical Clustering and SC3.

NMF: F-norm formulation is adopted to cluster scRNA-seq data.

rNMF: similar to the proposed rssNMF but without the graph regularization term: ${\displaystyle {min}_{W,H,S}\left|\right|X-WH-S|{|}_F^2+\alpha |\left|S\right|{|}_{1}s.t.W\ge 0,H\ge 0}$

ssNMF: ssNMF removes the matrix S which copes with noises and outliers but keep the graph regularization term: ${\displaystyle {min}_{W,H}\left|\right|X-WH|{|}_F^2+\beta Tr\left(HL{H}^T\right)s.t.W\ge 0,H\ge 0}$

KMeans: a canonical distance-based iterative algorithm. Euclidean distance is used in the experiment.

Hierarchical clustering: we use a division-based algorithm which initially starts with all observations in a single cluster and divide samples until each cluster only contain one observation. Euclidean distance is used in the experiment and “ward.D2” method was used in R hclust() function.

SC3: consensus clustering algorithm for single-cell RNA-seq data, which is a benchmark in comparison to various clustering methods in scRNA-seq clustering analysis (Duò, Robinson & Soneson, 2018; Kiselev et al., 2017).

After factorization, we use KMeans to cluster matrix H and obtain the clustering results. The other common method is categorizing sample j directly by the largest coefficient in column h_j.

To evaluate the clustering performance of robust semi-supervised NMF, we use adjusted Rand index (ARI, Supplementary File) which ranges from −1 to +1. Since cell labels are available so ARI can be calculated to measure the similarity between two data clustering results.

Prior genes selection and parameter setting

For all the datasets, we use differentially expressed genes as marker genes. For all the clusters in each dataset, we use R package “Deseq2” to identify the differentially expressed genes and select 20 differentially expressed genes with largest variance as cell marker genes for each group to construct weighting matrix. We set the same parameter value for rNMF, ssNMF and rssNMF and results are obtained by running the six methods 30 times on each dataset. For the factoring matrices W and H, we use standard normal distribution for random initialization.

Experimental results

Table 2 shows the clustering results on twelve datasets and we can see that rssNMF performs better than the most methods on almost all datasets. Here, we consider a basic circumstance that we only know one group of cells in the experiment and shows the results that one group of genes is added. All these datasets can be downloaded on NCBI Gene Expression Omnibus (GEO). The final results are evaluated by taking the average value of ARI over 30 runs. Compared with SC3, rssNMF gets better performance over seven datasets. We notice that the overall ARI on these seven datasets is relatively low (smaller than 0.4) which indicates that the cells hard to discriminate. Under such circumstance, the prior information can help improve the clustering accuracy than purely unsupervised clustering algorithm.

Our rssNMF model has two critical parameters: regularization parameter α and β. Figure 2 and Fig. S1 shows the how the average performance of rssNMF varies with the parameters α and β. As we can see in the Fig. 2, the performance of rssNMF is quite stable with respect to the parameter α on all four datasets. If α∕2 > max_ij(X − WH)_ij, all the elements in S will be zero, thus when the parameter α is large enough, rNMF will be equivalent to the original NMF.

The parameter β controls the weights of the graph regularization term and we test different weights range from 2 to 20,000 (Fig. S1). The performances are stable when β is smaller than 200 in the four datasets. When β is too large, clustering accuracy decreases a lot for the graph regularization term dominates the whole formulation and covers the main information.

To investigate how the prior information influences the clustering accuracy, we test the performance of rssNMF with respect to different number of group information. For each group in one dataset, 20 marker genes are selected as prior information to compute the weighting matrix. Figure S2 shows how the ARI changes as the group marker information increases (20 markers genes of one group, 40 markers gene for two groups and so on). We can see that as the number of group information increases, the performance fluctuates but clustering accuracy is higher than only adding one group information.

To evaluate the clustering stability, we calculate the consensus matrix for NMF, rNMF and rssNMF over 30 runs. A consensus matrix M is a n × n matrix that stores, for each pair of samples, the proportion of clustering runs in which two samples are clustered together. A consensus matrix can be obtained by taking the average over the connectivity matrices for all runs. The connectivity matrix C is also an n × n matrix and defined based on a single run: ${\displaystyle C_{ij}=\left\{\begin{array}{c}1,i\mathrm{and}j\text{are in the same clulster}\hfill \\ 0,i\mathrm{and}j\text{are not in the same cluster}\hfill \end{array}\right.}$

So, if sample i and j are in the same cluster for 30 runs, then M_ij will equal to 1, corresponding to red (0 corresponds to blue). Figure 3 present the factorizing matrices W, H and consensus matrix taking the average over 30 runs and the consensus matrix from NMF, rNMF and rssNMF on the dataset Yan. The color bar represents different group of cells or genes. More specifically, we assign genes to different cell groups according to the gene score, which is described in section Biological Significance Analysis section. The values in each row of W an column of H have been normalized between 0 and 1. We can see that all of the three methods have a stable clustering result on most samples and only small part of samples are clustered into multiple clusters in dataset Yan. Generally, rssNMF has the highest clustering accuracy and better stability than the other two methods.

Biological significance analysis

In this section, we choose dataset Yan for biological significance analysis and dataset Yan are human embryonic stem cells including zygote, 2-cell, 4-cell, 8-cell, morula and blastocyst. In Fig. 3, each column of matrix W defines a metagene (biological procesess or pathways) or a metasample (cell cluster) then entry w_ij can be regarded as the coefficient of gene I in metagene or metasample j (Brunet et al., 2004). Correspondingly, a column vector h_j represents the expression contribution of sample j to k biological processes, and H_ij can be seen as the contribution of metagene or metasample i in sample j. The class of samples is denoted by column annotation bar. We can see that for some samples (columns of H), it is hard to visually determine which cluster it belongs to since the depth of color is close in different clusters. Another obvious feature is that compared with NMF and rNMF, the pattern of H in rssNMF is clearer and more concentrated, which makes the classification characteristics of samples more explicit. We suppose that the addition of prior information makes the differences between different clusters more obvious, thus achieving better classification results.

For the basis matrix W, the ith row of W represents the contribution of gene i to all metagenes so the expression level of a specific gene in one cell is determined by the linear combination of its contribution to all biological processes. From Fig. S3 we can see that a gene can participate in multiple biological processes, but it is more reasonable to focus on the process with the largest value. We define a gene score to determine cluster-specific genes and the assignment of cluster-specific genes is shown in row annotation bar of matrix W. The gene score for the ith gene is: ${\displaystyle Gene\text{\_}score\left(i\right)=1+\frac{1}{log2\left(k\right)}\sum _{j=1}^{k}p\left(i,j\right)log2\left(p\left(i,j\right)\right),}$

where p(i, Ω) is the probability that the ith gene contributes to cluster Ω, i.e., $p\left(i,\Omega \right)=W_{i\Omega }∕{\sum }_{j=1}^k{W}_{ij}$. The gene score is a real value and ranges from [0, 1]. The higher gene score is, the more cluster-specific the corresponding gene. Simply we can think that a gene is highly-expressed if the gene has high Gene_score in a specific cluster over other clusters.

As shown in the A, D and G which are the basis matrix W from NMF, rNMF and rssNMF in the Fig. 3, we assign all the genes to a specific cluster using gene score by different colors as shown in the left color bar. For the basis matrix W, the rows indicate genes and the columns indicate clusters andthe depth of the color in the matrix indicates the size of the value in W. Deep color represents a larger value and we can assign the cluster genes visually. After classifying each gene, we need to validate the biological significance of feature genes detected in different cluster. Table S1 present the results of enrichment analysis for cluster-specific genes.

We use the KEGG database to investigate the functional genes selected for all clusters. KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies. Through functional enrichment analysis of genes in each cell cluster, we can understand what biological processes are involved in the cell cluster and further infer the role of the cluster in the biological microenvironment.

Embryonic stem cells are highly undifferentiated cells. It has developmental omnipotence and can differentiate all tissues and organs of adult animals, including germ cells. The research and utilization of ES cells is one of the core topics in the field of bioengineering. The enrichment results (Table S1) show that embryonic stem cells are in an active metabolic and proliferative state, such as Ribosome pathway in cluster 3, 5, 6, 7 and 8. Ribosomes are organelles that are responsible for making proteins and are widespread. The manufacture of ribosomes requires hundreds of cytokines that are not found in mature particles. In the absence of these factors, ribosome production will stagnate. Once ribosome production is stagnant, cell growth is terminated even under optimal growth conditions. Besides, various metabolic pathways such as TCA cycle, oxidative Phosphorylation and cholesterol metabolism also occur in multiple clusters, indicating that these clusters are in an active proliferation state. In addition, we can infer the differentiation path of stem cells from the results of enrichment. Cluster 6 contains two distinctive pathways, Cardiac muscle contraction and Retrograde endocannabinoid signaling. Both Cardiac muscle contraction and Retrograde endocannabinoid signaling are involved in nerve signal transduction and are usually active in the nervous system and cardiac conduction system. Endogenous cannabinoids (endocannabinoids) serve as retrograde messengers at synapses in various regions of the brain. Therefore, it can be inferred that cells in cluster 6 are differentiating toward the ectoderm and will form epidermis and nervous system in the future. Thermogenesis pathway appears in cluster 5, 6 and 7 and this pathway usually activates in brown adipose tissues. Brown adipose tissues originate from mesenchymal stem cells which differentiate from ectoderm. This is consistent with our view that cluster 7 belongs to ectoderm cell groups.