Bioinformatics Advance Access originally published online on May 5, 2007
Bioinformatics 2007 23(12):1495-1502; doi:10.1093/bioinformatics/btm134
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sparse non-negative matrix factorizations via alternating non-negativity-constrained least squares for microarray data analysis
College of Computing, Georgia Institute of Technology, 266 Ferst Drive, Atlanta, GA 30332, USA
*To whom correspondence should be addressed.
| ABSTRACT |
|---|
|
|
|---|
Motivation: Many practical pattern recognition problems require non-negativity constraints. For example, pixels in digital images and chemical concentrations in bioinformatics are non-negative. Sparse non-negative matrix factorizations (NMFs) are useful when the degree of sparseness in the non-negative basis matrix or the non-negative coefficient matrix in an NMF needs to be controlled in approximating high-dimensional data in a lower dimensional space.
Results: In this article, we introduce a novel formulation of sparse NMF and show how the new formulation leads to a convergent sparse NMF algorithm via alternating non-negativity-constrained least squares. We apply our sparse NMF algorithm to cancer-class discovery and gene expression data analysis and offer biological analysis of the results obtained. Our experimental results illustrate that the proposed sparse NMF algorithm often achieves better clustering performance with shorter computing time compared to other existing NMF algorithms.
Availability: The software is available as supplementary material.
Contact: hskim{at}cc.gatech.edu, hpark{at}acc.gatech.edu
Supplementary information: Supplementary data are available at Bioinformatics online.
| 1 INTRODUCTION |
|---|
|
|
|---|
In many data-mining problems, dimension reduction is imperative for efficient manipulation of massive quantity of high-dimensional data. The subspace method has demonstrated its success in numerous pattern recognition tasks including efficient classification (Kim et al., 2005), clustering (Ding et al., 2002) and fast search (Berry et al., 1999). There are two general approaches for reducing dimensionality, i.e. feature extraction and feature selection. Feature extraction is transforming the existing features into a lower dimensional space, while feature selection is selecting a subset of the existing features without a transformation. For feature extraction, principal component analysis (PCA), linear discriminant analysis (LDA) and non-negative matrix factorization (NMF) have been widely used. Many practical pattern recognition problems require non-negativity constraints. For example, pixels in digital images and chemical concentrations in bioinformatics are non-negative. NMF is a useful technique in approximating these high-dimensional data.
Given a non-negative matrix A of size
, where each column of A corresponds to a data point in the m-dimensional space, and a positive integer
, NMF finds two non-negative matrices
and
so that
|
| (1) |
|
| (2) |
0 means that all elements of W and H are non-negative. Due to Since NMF may give us direct interpretation due to non-subtractive combinations of non-negative basis vectors, it has recently received much attention and it has been applied to many interesting problems including text data mining (Chagoyen et al., 2006; Lee and Seung, 1999; Pauca et al., 2004) gene expression data analysis (Brunet et al., 2004; Carmona-Saez et al., 2006; Gao and Church, 2005; Kim and Tidor, 2003; Maher et al., 2006), microarray comparative genomic hybridization (aCGH) data (Carrasco et al., 2006) and functional characterization of gene lists (Pehkonen et al., 2005). One of the interesting properties of NMF is that it often generates sparse basis vectors that allow us to discover parts-based basis vectors. NMF generated holistic basis images instead of parts-based basis images for a facial image dataset in the results presented in Li et al. (2001) and Hoyer (2004). Several approaches (Dueck et al., 2005; Hoyer, 2004; Pascual-Montano et al., 2006; Pauca et al., 2006) have been proposed to explicitly control the degree of sparseness of W and H. However, if strong sparsity constraints are imposed on the basis matrix W, some useful information for gene selection in microarray data analysis may be lost. Some other approaches (Gao and Church, 2005; Pauca et al., 2004) imposed sparsity constraints on only the coefficient matrix H, but there may be mathematical difficulties related to convergence.
In this article, we introduce a sparse NMF algorithm that can control the degree of sparseness in the coefficient matrix H via alternating non-negativity-constrained least squares, and apply it to microarray data analysis. The rest of this article is organized as follows. We give brief overviews on various NMF algorithms in Section 2.1. In Section 2.2, we introduce sparse NMF formulations and algorithms based on alternating non-negativity-constrained least squares and discuss their convergence properties. In Section 3, we describe stopping criteria, and present experimental results and biological analysis illustrating properties of the proposed sparse NMF. Summary is given in Section 4.
| 2 METHODS |
|---|
|
|
|---|
2.1 Review of NMF algorithms
Lee and Seung (1999, 2000) suggested NMF algorithms based on multiplicative update rules of W and H. The distance
|
|
q
k and 1
j
n, |
|
i
m and 1
q
k. The divergence D(A, WH) = Hoyer (2004) devised a sparse NMF algorithm based on the projected gradient descent method (SNMF/PGD) in order to constrain NMF to find solution with desired sparseness in W and H. To impose sparseness constraints on only one matrix W or H, this algorithm uses a multiplicative update rule for updating the counter matrix, which suffers from slow convergence. Puscual-Montano et al. (2006) claimed that non-smooth NMF (nsNMF) outperformed previous sparse NMF variants for their synthetic and real datasets. The nsNMF (Puscual-Montano et al., 2006) is also based on multiplicative update rules.
Pauca et al. (2006) proposed a constrained NMF (CNMF) formulation,
|
| (3) |
and ß are regularization parameters. A sparse NMF algorithm using the following least squares,
|
| (4) |
2.2 Sparse NMFs based on alternating non-negativity-constrained least squares
In order to enforce sparseness on W or H in the NMF presented in Equation (1), we introduce two formulations and the corresponding algorithms for sparse NMFs, i.e. SNMF/L for sparse W (where L denotes the sparseness imposed on the left factor) and SNMF/R for sparse H (where R denotes the sparseness imposed on the right factor). Our sparse NMF formulations that impose the sparsity on a factor of NMF utilize L1-norm minimization and the corresponding algorithms are based on alternating non-negativity constrained least squares (ANLS). Each sub-problem is solved by a fast non-negativity constrained least squares (NLS) algorithm (van Benthem and Keenan, 2004) that is improved upon the active set based NLS method. Bro and de Jong (1997) made a substantial speed improvement to Lawson and Hanson's algorithm (Lawson and Hqnson, 1974) for multiple right-hand-side cases. van Benthem and Keenan (2004) devised an algorithm that further improves the performance of NLS.
2.2.1 Formulations for Sparse NMFs
SNMF/R: To apply sparseness constraints on H, we formulate the following SNMF/R optimization problem:
|
| (5) |
|
| (6) |
|
| (7) |
SNMF/L: To impose sparseness constraints on W, we introduce the SNMF/L formulation:
|
| (8) |
2.2.2 Convergence properties of sparse NMF algorithms
We show the convergence property of the sparse NMF algorithms. Since the convergence properties of SNMF/L and SNMF/R are essentially the same, we will only discuss the case of SNMF/R in more detail. Under conditions of
,
, and
due to
, Equation (5) can rewritten as
|
| (9) |
| 3 RESULTS |
|---|
|
|
|---|
3.1 Datasets description
We used the leukemia gene expression dataset (ALLAML) (Golub et al., 1999) and the central nervous system tumors dataset (CNS) (Pomeroy et al., 2002). The ALLAML dataset contains acute lymphoblastic leukemia (ALL) that has B and T cell subtypes, and acute myelogenous leukemia (AML) that occurs more commonly in adults than in children. This gene expression dataset consists of 38 bone marrow samples (19 ALL-B, 8 ALL-T and 11 AML) with 5000 genes. The central nervous system dataset is composed of four categories of CNS tumors with 5597 genes. It consists of 34 samples representing four distinct morphologies: 10 classic medulloblastomas, 10 malignant gliomas, 10 rhabdoids and 4 normal cerebella. All datasets used by us contain only non-negative entries. All algorithms were implemented in Matlab 6.5 (MATLAB, 1992). The Matlab codes for NMF using divergence-based multiplicative update rules were obtained from Brunet et al. (2004) and modified to implement nsNMF (Puscual-Montano et al., 2006). All our experiments were performed on a P3 600 MHz machine with 512 MB memory.
3.2 Biclustering
We applied the non-negative factorization of Equation (1) to perform clustering analysis of a data matrix. The rows of a microarray data matrix A represent genes and the columns experiments. We can use the basis matrix W to divide the m genes into k gene-clusters and the coefficient matrix H to divide the n samples into k sample-clusters. Typically, gene i is assigned to gene-cluster q if the W(i, q) is the largest element in W (i, :) and sample j is assigned to sample-cluster q if the H(q, j) is the largest element in H(:, j).
3.3 Stopping criteria
NMF using divergence-based multiplicative update rules (i.e. NMF/DUR (Brunet et al., 2004) and nsNMF (Pascual-Montano et al., 2006)) in our implementation stops if
has not changed for more than 40 convergence tests (each made at 10 iterations), where
is the connectivity matrix of which entry is
if samples i and j belong to the same sample-cluster, and
if they belong to different sample-clusters.
For SNMF/L and SNMF/R, we tested convergence at every five iterations by the combined convergence criterion using the Karush-Kuhn-Tucker (KKT) optimality conditions and the convergence of positions of the largest elements in rows of W and columns of H. Our sparse NMF algorithms stop if both the positions of the largest elements in rows of W, i.e.
, and the positions of the largest elements in columns of H, i.e.
, have not changed for more than or equal to 10 convergence tests, where
is the position of the largest element in the i-th row of W and
is the positions of the largest element in the j-th column of H, and the following KKT conditions are satisfied. The KKT conditions for each objective function f(W, H) with non-negativity constraints W
0 and H
0 are
|
| (10) |
|
| (11) |
o be the KKT residual measured by the L1-vector norm,
|
| (12) |
W = #(min(W,
W)
0), and the number of the elements in H that did not converge yet, i.e.
H = #(min(H,
H)
0). We define the following normalized KKT residual:
|
| (13) |
|
| (14) |
1 is the
value in the first iteration and
is a tolerance. We used
3.4 Clustering performance comparison
To measure the performance of NMFs in clustering, we used purity and entropy. Suppose we are given l categories (true class labels), while NMF generates k clusters. Purity is given by
|
|
j
l). The larger the value of purity, the better the clustering performance. Entropy is defined as follows: |
|
The parameter
in Equation (5) is important in keeping
small. We set it to be the square of the maximal element in A. For the initialization of W in SNMF/R, the elements in the initial matrix W were randomly chosen and normalized so that the columns of the basis matrix W have unit L2-norm, i.e.
.
Tables 1 and 2 show the results of SNMF/R with various values of ß on the ALLAML dataset with k = 3 and on the CNS tumors dataset with k = 4, respectively. We compared our proposed SNMF/R with NMF based on divergence-based multiplicative update rules (Brunet et al., 2004; Lee and Seung, 2000). Average sparseness, purity and entropy were computed by running each algorithm five times with different random initializations. By increasing ß in SNMF/R, we could obtain a sparser H. SNMF/R algorithm achieved better clustering performance (higher purity, lower entropy) than NMF/DUR within certain range of ß with shorter computing time. By increasing
in SNMF/L, we could enhance the sparsity of W (results are not shown). SNMF/L can be applied to obtain parts-based basis vectors.
|
|
We compared our methods with other sparse NMF variants on the ALLAML dataset. We tested Hoyer's sparse NMF based on the projected gradient descent method by his Matlab implementation with the sparseness control parameter
= 0.5 suggested in Carmona-Saez et al. (2006) for biclustering of gene expression data due to reasonable results from numerous empirical tests. nsNMF generated average purity = 0.963 and average entropy = 0.108. The average percentages of elements in the range of
We used the model selection method proposed by Brunet et al. (2004) to determine the number of factors. We ran NMF algorithms 30 times to obtain the average connectivity matrix (i.e. consensus matrix) whose entries reflect the probability that samples i and j belong to the same cluster. To measure the dispersion of a consensus matrix C, we defined the dispersion coefficient (
) as
|
| (15) |
= 1 for a perfect consensus matrix (all entries = 0 or 1) and 0
< 1 for a scattered consensus matrix. After obtaining
k values for various k, we can determine the number of clusters from the maximal
k. Figure 1 illustrates the determination of the number of clusters in the CNS tumors dataset. SNMF/R with ß = 0.01 found perfect consensus matrices for
k
4. SNMF/R yielded finer consensus matrices (higher
k) than NMF/DUR for various k values.
|
3.5 Biological analysis
Figure 2 presents the matrices W and H obtained from SNMF/R with ß = 0.01 on the ALLAML leukemia dataset, which produced the lowest approximation error
|
A row vector of the basis matrix W indicates the contributions of a gene to the k biological pathways or processes (i.e. k columns of W). Genes can participate in more than one biological process. It is beneficial to investigate genes that have relatively large coefficient in each biological process. We selected factor-specific genes via the non-negative basis matrix
|
| (16) |
, i.e. Some chosen genes dominantly contribute to only single biological pathway or process. For instance, MB-1 gene (U05259) is most active in the first process. Transcription factor 7 (T-cell specific, HMG-box) (TCF7, X59871) is active in the second process, which is also known as T cell factor-1 (TCF-1). Some genes play a major role in the third process, for example, Interleukin 8 (IL8, M28130 [GenBank] ), DF D component of complement (adipsin) (CFD, M84526 [GenBank] ), Cystatin C (amyloid angiopathy and cerebral hemorrhage (CST3, M27891 [GenBank] ), Chemokine (C-X-C motif) ligand 2 (CXCL2, M57731), etc. Chemokine is a type of cytokines that bind to a specific cell-surface receptor and is critical to the functioning of both innate and adaptive immune responses. Total of 37 genes including MB-1, IL8, CFD and CST3 were the same genes as those found in Golub et al. (1999). Ribosomal protein S3 (RPS3, X57351) simultaneously participates in all three processes. This is reasonable since RPS3 is a housekeeping gene and ribosomal protein genes are usually over-expressed in some cancers. RPS3 encodes a ribosomal protein that is a component of the 40S subunit, where it forms a part of the domain where translation is initiated.
We used the Onto-Express (Draghici et al., 2003; Khatri et al., 2002) to investigate the enrichment of functional annotations of genes selected in each factor. Onto-Express starts by reading the input file that contains a list of GenBank accession numbers, and estimates the statistical significance of the enrichment of Gene-Ontology (GO) terms in the list with respect to a reference list. We used a list of all genes in the dataset as a reference array and hypergeometric distribution. Table 3 shows the enrichment of GO terms. We presented some significant biological processes for each factor, whose P-values were less than 0.01.
|
Figure 3 illustrates the matrices W and H obtained from SNMF/R with ß = 1.0 on the CNS tumors dataset, which produced the lowest approximation error for five runs with different random initializations of W. Only one sample (the 18th sample, Brain_MGlio_8) was incorrectly assigned to the third cluster (rhabdoids). Our gene selection method suggested total of 367 genes (cluster1: 42, cluster2: 93, cluster3: 168, cluster4: 64). To more thoroughly characterize sets of genes dominantly expressed in different factors, we used the Onto-Express. The number of genes corresponding to each GO category was compared with the number of genes expected for the GO category in the Affymetrix HuGeneFL array. Significant differences from the expected were calculated with hypergeometric distribution. Table 4 shows biological processes with a significance of P-value < 0.01. The biological processes showing significant representations in the first factor were serotonin receptor signaling pathway, feeding behavior, and central nervous system development. Serotonin (5-hydroxytryptamine, or 5-HT) is a monoamine neurotransmitter and is known to regulate human mood, emotion, sleep and appetite in the central nervous system. Two GenBank accession numbers (U49516 and M81778) for serotonin receptor signaling pathway were linked to the same gene: 5-hydroxytryptamine (serotonin) receptor 2C (HTR2C). The GO category of feeding behavior seems to be related with childhood brain tumors known as medulloblastomas. Genes involved in the second factor were (4 genes: U52155, M81886, M64752 [GenBank] , M81181 [GenBank] ) for potassium ion transport, (5 genes: X54673, M81886, M64752 [GenBank] , M19650, L32961) for synaptic transmission, (3 genes: U62801, Z19002, M93426) for central nervous system development. This second cluster contains malignant glioma that is a tumor arising from glial cells. Genes corresponding to cell adhesion, cell motility and inflammatory response were highly expressed in the third factor. Genes highly expressed in normal cerebella were (5 genes: U79245, U33632, U90065, L36069, D79998) for potassium ion transport, (6 genes: M13577, U92457, L76627, U18244, M58583, U79667) for synaptic transmission, (3 genes: U52969, M13577, U76421) for central nervous system development and (6 genes: S81944, U79245, S95936, U33632, U90065, L36069) for ion transport. Detailed description of the clusters of samples and genes selected for each factor via SNMF/R can be found in supplementary materials. We have shown that SNMF/R can be used for clustering, cancer class discovery, gene selection and biological process analysis.
|
|
| 4 CONCLUSION |
|---|
|
|
|---|
We present a novel sparse NMF algorithm via alternating non-negativity-constrained least squares. SNMF/R can be used for cancer class discovery and gene expression data analysis since it shows good biclustering performance and provides us with simple interpretation. This algorithm can be applied to many practical problems in bioinformatics and computational biology such as biomedical text mining and gene/protein microarray data analysis.
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
We would like to thank Dr Chris Ding, Dr Jean-Philippe Brunet, Dr Yuan Gao, Prof. Lars Eldén and Prof. Robert J. Plemmons for their valuable comments. In particular, we would also like to thank Prof. Chih-Jen Lin, Prof. Paul Tseng and Prof. Luigi Grippo for discussions on the convergence property. This material is based upon work supported in part by the National Science Foundation Grants ACI-0305543 and CCF-0621889. Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.
Conflict of Interest: none declared.
| FOOTNOTES |
|---|
Associate Editor: David Rocke
Received on November 2, 2006; revised on February 19, 2007; accepted on April 1, 2007
| REFERENCES |
|---|
|
|
|---|
Berry MW, et al. Matrices, vector spaces, and information retrieval. SIAM Rev (1999) 41:335–362.[CrossRef]
Berry MW, et al. Algorithms and applications for approximate nonnegative matrix factorization. Comput. Stat. Data Anal (2006) (to appear).
Bro R, de Jong S. A fast non-negativity-constrained least squares algorithm. J. Chemometrics (1997) 11:393–401.[CrossRef]
Brunet JP, et al. Metagenes and molecular pattern discovery using matrix factorization. Proc. Natl Acad. Sci. USA (2004) 101:4164–4169.
Carmona-Saez P, et al. Biclustering of gene expression data by non-smooth non-negative matrix factorization. BMC Bioinformatics (2006) 7:78.[CrossRef][Medline]
Carrasco DR, et al. High-resolution genomic profiles define distinct clinico-pathogenetic subgroups of multiple myeloma patients. Cancer Cell (2006) 9:313–325.[CrossRef][Web of Science][Medline]
Chagoyen M, et al. Discovering semantic features in the literature: a foundation for building functional associations. BMC Bioinformatics (2006) 7:41.[CrossRef][Medline]
Ding C, et al. Adaptive dimension reduction for clustering high dimensional data. (2002) In Proceedings of the 2nd IEEE International Conference on Data Mining: Maebashi, Japan.
Draghici S, et al. Onto-tools, the toolkit of the modern biologist: Onto-Express, Onto-Compare, Onto-Design and Onto-Translate. Nucleic Acids Res (2003) 31:3775–3781.
Dueck D, et al. Multi-way clustering of microarray data using probabilistic sparse matrix factorization. Bioinformatics (2005) 21(Suppl. 1):i144–i151.[Abstract]
Gao Y, Church G. Improving molecular cancer class discovery through sparse non-negative matrix factorization. Bioinformatics (2005) 21:3970–3975.
Golub TR, et al. Molecular classification of cancer: class discovery and class prediction by gene expression monitoring. Science (1999) 286:531–537.
Gonzales EF, Zhang Y. Accelerating the Lee-Seung algorithm for non-negative matrix factorization. Technical report (2005) Department of Computational and Applied Mathematics, Rice University.
Grippo L, Sciandrone M. On the convergence of the block nonlinear Gauss-Seidel method under convex constraints. Operations Res. Lett (2000) 26:127–136.[CrossRef]
Hoyer PO. Non-negative matrix factorization with sparseness constraints. J. Machine Learning Res (2004) 5:1457–1469.
Khatri P, et al. Profiling gene expression using onto-express. Genomics (2002) 79:266–270.[CrossRef][Web of Science][Medline]
Kim H, et al. Dimension reduction in text classification with support vector machines. J. Machine Learning Res (2005) 6:37–53.
Kim PM, Tidor B. Subsystem identification through dimensionality reduction of large-scale gene expression data. Genome Res (2003) 13:1706–1718.
Lawson CL, Hanson RJ. Solving Least Squares Problems (1974) NJ: Prentice-Hall, Englewood Cliffs.
Lee DD, Seung HS. Learning the parts of objects by non-negative matrix factorization. Nature (1999) 401:788–791.[CrossRef][Medline]
Lee DD, Seung HS. Algorithms for non-negative matrix factorization. (2000) Proceedings of Neural Information Processing Systems. 556–562.
Li SZ, et al. Learning spatially localized parts-based representations. (2001) In Proceedings IEEE Conference on Computer Vision and Pattern Recognition (CVPR). 207–212.
Maher EA, et al. Marked genomic differences characterize primary and secondary glioblastoma subtypes and identify two distinct molecular and clinical secondary glioblastoma entities. Cancer Res (2006) 66:11502–11513.
MATLAB. User's Guide. (1992) The MathWorks, Inc. Natick, MA 01760.
Pascual-Montano A, et al. Nonsmooth nonnegative matrix factorization (nsNMF). IEEE, Trans. Pattern Anal. Machine Intell (2006) 28:403–415.[CrossRef]
Pauca VP, et al. Text mining using non-negative matrix factorizations. (2004) Proceedings SIAM International Conference on Data Mining (SDM'04).
Pauca VP, et al. Nonnegative matrix factorization for spectral data analysis. Linear Algebra and Applications (2006) (to appear).
Pehkonen P, et al. Theme discovery from gene lists for identification and viewing of multiple functional groups. BMC Bioinformatics (2005) 6:162.[CrossRef][Medline]
Pomeroy SL, et al. Prediction of central nervous system embryonal tumour outcome based on gene expression. Nature (2002) 415:436–442.[CrossRef][Medline]
Tibshirani R. Regression shrinkage and selection via LASSO. J. Roy. Statist. Soc. B (1996) 58:267–288.
van Benthem MH, Keenan MR. Fast algorithm for the solution of large-scale non-negativity-constrained least squares problems. J. Chemometrics (2004) 18:441–450.[CrossRef]
This article has been cited by other articles:
![]() |
D. Greene, G. Cagney, N. Krogan, and P. Cunningham Ensemble non-negative matrix factorization methods for clustering protein-protein interactions Bioinformatics, August 1, 2008; 24(15): 1722 - 1728. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||








