A Robust Manifold Graph Regularized Nonnegative Matrix Factorization Algorithm for Cancer Gene Clustering

Detecting genomes with similar expression patterns using clustering techniques plays an important role in gene expression data analysis. Non-negative matrix factorization (NMF) is an effective method for clustering the analysis of gene expression data. However, the NMF-based method is performed within the Euclidean space, and it is usually inappropriate for revealing the intrinsic geometric structure of data space. In order to overcome this shortcoming, Cai et al. proposed a novel algorithm, called graph regularized non-negative matrices factorization (GNMF). Motivated by the topological structure of the GNMF-based method, we propose improved graph regularized non-negative matrix factorization (GNMF) to facilitate the display of geometric structure of data space. Robust manifold non-negative matrix factorization (RM-GNMF) is designed for cancer gene clustering, leading to an enhancement of the GNMF-based algorithm in terms of robustness. We combine the l2,1-norm NMF with spectral clustering to conduct the wide-ranging experiments on the three known datasets. Clustering results indicate that the proposed method outperforms the previous methods, which displays the latest application of the RM-GNMF-based method in cancer gene clustering.


Introduction
With the progressive implementation of human whole-genome and microarray technologies, it is possible to simultaneously observe the expressions of numerous genes in different tissue samples. By analyzing gene expression data, genes with varying expressions in tissues and their relationships may be identified to figure out the pathogenic mechanism of cancers based on genetic changes [1]. Recently, cancer classification based on gene expression data has become a hot research topic in bioinformatics.
Due to the fact that the analysis of genome-wide expression patterns can provide unique perspectives into the structure of genetic networks, the clustering technique has been used to analyze gene expression data [2,3]. Cluster analysis is the most widespread statistical techniques for analyzing massive gene expression data. Its major task is to classify genes with similar expressions to discover groups of genes with identical features or similar biological functions, in order that people can acquire a deeper understanding about the essence of many biological phenomena such as gene functions, development, cancer, and pharmacology [4].
Currently, it has been shown that non-negative matrix factorization (NMF) [5,6] is superior to the hierarchical clustering (HC) and self-organizing map (SOM) [7] in the application of cancer samples in clustering gene expression data. Over the past few years, the NMF-based method has been used for the gene expressions of statistically analyzing data for clustering [8][9][10][11][12][13]. The main idea is to approximately factorize a non-negative data matrix into a product of two non-negative matrices, which makes sure that all elements of the matrices are non-negative. Therefore, the appearance of NMF has attracted considerable attention. Recently, various variants based on the original NMF have been developed by modifying the objective function or the constraint conditions [14][15][16]. For instance, Cai et al. proposed graph regularized non-negative matrix factorization (GNMF), giving forth to the neighboring geometric structure. It illustrates the nearest neighbor graph that preserves the neighborhood information of high-dimensional space in low-dimensional space. The GNMF reveals the intrinsic geometrical structure by incorporating a Laplacian regularization term [17], which is effective for solving clustering problems. After that, the sparse NMF [18] was proposed with sparse constraints upon the basis matrices and coefficient matrices factored by the NMF so that the sparseness may be reflected from data. The non-smooth NMF [19] can realize the global or local sparseness [20] by making basis and encoding matrices sparse simultaneously. For the sake of enhancing the robustness of the GNMF-based method in gene clustering, we propose improved robust manifold non-negative matrix factorization (RM-GNMF) by making use of the combination of l 2,1 -norm and spectral clustering with Laplacian regularization, leading to the internal geometry of data representations. It facilitates the display purposes of intrinsic geometric structure of the cancer gene data space.
This paper is organized as follows. In Section 2, we give a brief review on the NMF and the GNMF. In Section 3, we propose an improved RM-GNMF algorithm. In Section 4, we give the experimental results comparing with the previous methods. Finally, the conclusions are drawn in Section 5.

The NMF-Based and GNMF-Based Method
The NMF-based method [5] is a linear and non-negative approximate data description for non-negative matrices. We consider an original decomposed matrix X of size m × n, where m represents data characteristics and n represents the number of samples. Based on the NMF method, the matrix X is decomposed into two non-negative matrices W ∈ R m×r and H ∈ R n×r , i.e., where r ≤ min(m, n). For a given decomposition X = W H T , sample m can be divided into r classes according to matrix H T . Each sample is placed in the highest metagene expression level in the sample, meaning that if H T ij is the largest in column j, then sample j is assigned to class i . Using the square of the Euclidean distance between X and W H T [21], we have the objective function of the NMF method According to the iterative update algorithms [6], the NMF-based method is performed on the basis of multiplicative update rules of W and H given by In order to overcome the limitation of the NMF-based method, Cai et al. [17] proposed the GNMF-based method, in which an affinity graph is generated to encode the geometrical information followed by a matrix factorization with respect to the graph structure. In contrast to the NMF-based method, it has a regular graph constraint, which preserves the advantage of the local sparse representation of the NMF-based method and preserves the similarity between the original data points after dimensionality reduction.
There are several weighting schemes, such as zero-one weighting, heat kernel weighting, and Gaussian weighting [17]. In what follows, we consider the zero-one weighting described as Based on the weight matrix Q, we obtain the objective function of the GNMF method given by where Tr(·) denotes the trace of matrices and L = D − Q. D is a diagonal matrix whose entries are column or row sums of Q with D jj = ∑ l Q jl [22]. The regularization parameter λ ≥ 0 can be used for the smoothness control of the new representation. By the iterative algorithms to minimize the objective function O GN MF , we achieve the updating rules

The RM-GNMF-Based Method for Gene Clustering
So far, we have described the NMF-based method and GNMF-based method. In what follows, we seek RM-GNMF gene clustering by making use of the combination of l 2,1 -norm and spectral clustering with the Laplacian regularization.

The l 2,1 -Norm
The l 2,1 -norm of a matrix was initially employed as a rotational invariant l 1 -norm [23], which was usually used for multi-task learning [24,25] and tensor factorization [26]. Instead of using l 2 -norm-based loss function that is sensitive to outliers, we resort to the l 2,1 -norm-based loss function and regularization [23], which is convergence-proved.
For the sake of getting over the drawbacks of the NMF-based method and enhancing the robustness of the GNMF-based method, we employ the l 2,1 -norm for the matrix factorization in the RM-GNMF-based method. For a non-negative matrix X of size m × n, the l 2,1 -norm of matrice X is defined as where data vectors are arranged in columns, and the l 2,1 -norm calculates the l 2 -norm for column vectors first. Subsequently, the matrix factorization assignment becomes

Spectral-Based Manifold Learning to Constrained GNMF
The spectral method is a classical method of analysis and algebra in mathematics. It is widely used in low dimensional representation and clustering problems of high dimensional data [27,28]. A relational matrix describing the similarity of the pair of data points is defined according to the given sample dataset, and the eigenvalues and eigenvectors of the matrices are calculated. After that, the appropriate eigenvectors are selected and the low dimensional embedding of the data is obtained. The degree matrices are defined on a given graph, such as an adjacency matrix of the graph, a Laplacian matrix, and so on [22].
Based on the spectrum of the matrices with respect to the graph, spectral theory further reveals the information contained in the graph , and establishes the connection between the discrete space and the continuous space through the techniques of geometry, analysis, and algebra. It has a wide range of applications in manifold learning. In this section, the p-nearest neighbor method can be used for establishing the relationship between each data point and its neighborhood.
For a data matrix X ∈ R m×n , we treat each column of X as a data point and each data point as a vertex, respectively. The p-nearest-neighbor graph G can be constructed with n vertices. Then the symmetric weight matrix Q ∈ R n×m is generated, where the element q ij denotes the weight of the edge joining vertices i and j and the value of q ij is given by where N p (x i ) denotes the set of p-nearest neighbors of x i . It is obvious that the matrix Q represents the affinity between the data points.
There is an assumption about manifold. Namely, if two data points x i and x j are close in the intrinsic geometric structure of the data distribution, then their presentations under a new basis will be close [29]. Therefore, we define the relationship as follows where m i and m j denote the mappings of x i and x j , respectively. The degree matrix D is a diagonal matrix given by d ii = ∑ j q ij . Obviously, d ii is the sum of all the similarities related with x i . Then, the graph Laplacian matrix is given by The graph embedding can be written as In the RM-GNMF-based method, we combine the GNMF-based method with the spectral clustering, resulting in the l 2,1 -norm constrained GNMF as follows where λ ≥ 0 is the regularization parameter. We resort to the augmented Lagrange multiplier (ALM) method to solve the above problem.
For an auxiliary variable Z = X − W H T , we rewrite the O RMGN MF in Equation (13) satisfying the constraints Z − X + W H T = 0 and H T H = I. Then, we define the augmented Lagrangian function satisfying the constraint H T H = I, where µ is the step size of update and Λ is the Lagrangian multiplier.
To optimize Equation (15), we rewrite the objective function to get the following task satisfying the constraint H T H = I.

Computation of Z
For the given W and H, we can solve Z in Equation (15) by using the update of Z related to the following issue We need the following Lemma to solve Z in Equation (15). Please see the Appendix A1 for a detailed proof. Lemma 1. Given a matrix W = [w 1 , · · · , w n ] ∈ R m×n and a positive scalar λ, Z * is the optimal solution of and the i-th column of Z * is given by

Computation of W and H
For the given other parameters, we solve the optimal W. The update of W amounts to solve Let X − Z + Λ µ = M. The problem in Equation (20) can be rewritten as If the partial derivative of W is set to be 0, we obtain Then, we derive the optimal H. Taking W = MH, the update problem of H can be expressed as satisfying the constraint H T H = I. We have Therefore, the optimal H r+1 can be achieved by counting eigenvectors of H r+1 = (h 1 , · · · , h k ).

Updating of Λ and µ
The update standard of Λ and µ can be described as follows where p is the nearest neighbor graph parameter. The detailed process of the RM-GNMF-based method is listed in Algorithm 1.

Input:
The dataset X = [x 1 , x 2 , · · · , x n ] ∈ R m×n , a predefined number of clusters k, parameters µ, λ, the nearest neighbor graph parameter p, maximum iteration number t max . Initialization:

Repeat
Fix other parameters, and then update Z by formula (17); Fix other parameters, and then update W by: where U and V are the left and right singular values of the SVD decomposition; Fix other parameters, and then update Λ and µ by formulas (26)(27); t = t + 1; Until t ≥ t max . Output: matrix W ∈ R m×k , matrix H ∈ R k×n .

Results and Discussion
In this section, we evaluate the performance of the proposed method on the three gene expression datasets. We compare the RM-GNMF-based method with the NMF-based method [6], the l 2,1 -NMF-based method [23], the LNMF-based method [20], and the GNMF-based method [17].

Datasets
In order to evaluate the performance of proposed RM-GNMF algorithm, the clustering experiment was conducted on several gene expressions datasets of cancer patients. Three classical genetic datasets were used in the experiment, including leukemia [1], colon, and GLI_85 [30]. These gene expression datasets are downloaded from: http://featureselection.asu.edu/datasets.php.The colon cancer datasets consist of the gene expression profiles of 2000 genes for 62 tissue samples among which 40 are colon cancer tissues and 22 are normal tissues. The leukemia datasets consist of 7129 genes and 72 samples (47 ALL and 25 AML).
A brief description of experimental datasets is described in Table 1. More detailed information on these datasets can be found in the relevant references, and these datasets are available for download from the reference website.

Evaluation Metrics
For the sake of evaluating the clustering results, we use the clustering accuracy and normalized mutual information (NMI) to demonstrate the performance of the proposed algorithm.
Clustering accuracy can be calculated as where c i is the cluster label of x i , and l i is the true class label i-th sample, n denotes the total number of samples, and δ(map(c i ), l i ) is a delta function. If map(c i ) = l i , we obtain δ(map(c i ), l i ) = 1, where map(c i ) is the mapping function that maps the cluster label c i into the actual label l i . Otherwise, we have δ(map(c i ), l i ) = 0. We can find the best mapping by the Kuhn-Munkres method [31]. NMI can be described as where n i is the size of the i-th cluster andn j is the size of the j-th class, n i,j is the number of data between the intersections, and N denotes the number of clusters. We perform 100 experiments under each target feature dimension, taking the mean of the accurate and NMI values as the experimental results.

Parameter Selection
The RM-GNMF-based method involves two essential parameters, i.e., the regularization parameter λ and the regularity coefficient µ determining the penalty for infeasibility.
We set the parameters λ and µ in the range of λ ∈ {0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1} and µ ∈ {10 −1 , 10 −2 , 10 −3 , 10 −4 , 10 −5 , 10 −6 , 10 −7 } . We use the cross-validation method to get the best parameter values λ = 0.05 and µ = 10 −3 . In order to intuitively to analyze the influence of parameters λ and µ of the RM-GNMF-based method on the accuracy of clustering, Figure 1 shows the variation on clustering accuracy when the two parameters are modified. The three subgraphs in Figure 1 correspond to three gene expression datasets respectively. As can be viewed in Figure 1, the parameter µ = 10 −3 can get higher ACC. With the change of regular parameter λ, the change of ACC is relatively flat, and the clustering accuracy is higher when the value of λ is smaller. Therefore, we set λ = 0.05, µ = 10 −3 in the follow-up experiments.

Clustering Results
In Table 2, we demonstrate the clustering results on the colon, GLI_85 and leukemia datasets, respectively. Reported is the mean of clustering results from 100 runs of different NMF methods together. Table 2. Clustering results on different datasets. NMF: non-negative matrix factorization; GNMF: graph regularized non-negative matrices factorization; RM-GNMF: robust manifold non-negative matrix factorization; NMI: normalized mutual information. It can be found that the RM-GNMF-based method outperforms the original NMF-based method, while the RM-GNMF-based method achieves the best performance compared with the other three datasets. The clustering accuracies of the RM-GNMF-based method are 66.13%, 75.29%, and 65.28% for the colon, GLI_85, and leukemia datasets, respectively.

Colon
Our tests on several gene expression profiling datasets of cancer patients consistently indicate that the RM-GNMF-based method achieves significant improvements in comparison with the NMF-based method, the l 2,1 -NMF-based method, the LNMF-based method, and the GNMF-based method, in terms of cancer prediction accuracy.
As shown in Figure 2, the RM-GNMF-based method always gives birth to better clustering results than other NMF-based method using the three original datasets. To demonstrate the robustness of our approach to data changes, we add uniform noise onto the three gene expression datasets. A disturbed matrix Y noise is generated by adding independent uniform noise, defined as follows: where Y is the original matrix, r is a random number generated by a uniform distribution on the interval [0, max], and max is the maximum expression of Y . The experimental results with noise added are shown in Figure 3. It can be seen that the clustering result of RM-GNMF algorithm is still stable with the addition of noise, which shows that RM-GNMF algorithm is robust. In order to verify the results obtained from the algorithms in the experiments, we import the clustering result of the comparison methods into STAC web platform to perform the statistical test (http://tec.citius.usc.es/stac/). We selected the Friedman test of non-parametric multiple groups; the significance level is 0.05. The analysis results obtained are presented in Tables 3 and 4.  From the above test results it can be concluded that H0 is rejected. Hence, we believe that the clustering results of five algorithms are significantly different.

Conclusions
We have proposed the RM-GNMF-based method with the l 2,1 -norm and spectral-based manifold learning. This algorithm is suitable for cancer gene expression data clustering with an elegant geometric structure. Our tests on several gene expression profiling datasets of cancer patients consistently indicate that the RM-GNMF-based method achieves significant improvements in comparison with the NMF-based method, the l 2,1 -NMF-based method, the LNMF-based method, and the GNMF-based method, in terms of cancer prediction accuracy and robustness.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Lemma A1. With the given matrix W = [w 1 , · · · , w n ] ∈ R m×n and the positive scalar λ, Z * is the optimal solution of and the i-th column of Z * can be calculated as Proof. The objection function in Equation (A1) is equivalent to the following equation which can be solved in a decoupled manner After taking derivative with respect to z i , we get where r is a subgradient vector and ||r|| 2 ≤ 1. For z i = 0, we get where λ ≥ ||w i ||. For z i = 0, we get Combining Equation (A6) with Equation (A7), we obtain where α = ||z i || 2 ||z i || 2 +λ > 0. Plugging z i in Equation (A8) to Equation (A7), we solve α, which is substituted back into Equation (A8). After performing the above-mentioned steps, we obtain