New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Singular value decomposition for genomewide expression data processing and modeling

Contributed by David Botstein
Abstract
We describe the use of singular value decomposition in transforming genomewide expression data from genes × arrays space to reduced diagonalized “eigengenes” × “eigenarrays” space, where the eigengenes (or eigenarrays) are unique orthonormal superpositions of the genes (or arrays). Normalizing the data by filtering out the eigengenes (and eigenarrays) that are inferred to represent noise or experimental artifacts enables meaningful comparison of the expression of different genes across different arrays in different experiments. Sorting the data according to the eigengenes and eigenarrays gives a global picture of the dynamics of gene expression, in which individual genes and arrays appear to be classified into groups of similar regulation and function, or similar cellular state and biological phenotype, respectively. After normalization and sorting, the significant eigengenes and eigenarrays can be associated with observed genomewide effects of regulators, or with measured samples, in which these regulators are overactive or underactive, respectively.
DNA microarray technology (1, 2) and genome sequencing have advanced to the point that it is now possible to monitor gene expression levels on a genomic scale (3). These new data promise to enhance fundamental understanding of life on the molecular level, from regulation of gene expression and gene function to cellular mechanisms, and may prove useful in medical diagnosis, treatment, and drug design. Analysis of these new data requires mathematical tools that are adaptable to the large quantities of data, while reducing the complexity of the data to make them comprehensible. Analysis so far has been limited to identification of genes and arrays with similar expression patterns by using clustering methods (4–9).
We describe the use of singular value decomposition (SVD) (10) in analyzing genomewide expression data. SVD is also known as Karhunen–Loève expansion in pattern recognition (11) and as principalcomponent analysis in statistics (12). SVD is a linear transformation of the expression data from the genes × arrays space to the reduced “eigengenes” × “eigenarrays” space. In this space the data are diagonalized, such that each eigengene is expressed only in the corresponding eigenarray, with the corresponding “eigenexpression” level indicating their relative significance. The eigengenes and eigenarrays are unique, and therefore also datadriven, orthonormal superpositions of the genes and arrays, respectively.
We show that several significant eigengenes and the corresponding eigenarrays capture most of the expression information. Normalizing the data by filtering out the eigengenes (and the corresponding eigenarrays) that are inferred to represent noise or experimental artifacts enables meaningful comparison of the expression of different genes across different arrays in different experiments. Such normalization may improve any further analysis of the expression data. Sorting the data according to the correlations of the genes (and arrays) with eigengenes (and eigenarrays) gives a global picture of the dynamics of gene expression, in which individual genes and arrays appear to be classified into groups of similar regulation and function, or similar cellular state and biological phenotype, respectively. These groups of genes (or arrays) are not defined by overall similarity in expression, but only by similarity in the expression of any chosen subset of eigengenes (or eigenarrays). Upon comparing two or more similar experiments, with a regulator being overactive or underactive in one but normally expressed in the others, the expression pattern of one of the significant eigengenes may be correlated with the expression patterns of this regulator and its targets. This eigengene, therefore, can be associated with the observed genomewide effect of the regulator. The expression pattern of the corresponding eigenarray is correlated with the expression patterns observed in samples in which the regulator is overactive or underactive. This eigenarray, therefore, can be associated with these samples.
We conclude that SVD provides a useful mathematical framework for processing and modeling genomewide expression data, in which both the mathematical variables and operations may be assigned biological meaning.
Mathematical Framework: Singular Value Decomposition
The relative expression levels of N genes of a model organism, which may constitute almost the entire genome of this organism, in a single sample, are probed simultaneously by a single microarray. A series of M arrays, which are almost identical physically, probe the genomewide expression levels in M different samples—i.e., under M different experimental conditions. Let the matrix ê, of size Ngenes × Marrays, tabulate the full expression data. Each element of ê satisfies 〈nêm〉 ≡ e_{nm} for all 1 ≤ n ≤ N and 1 ≤ m ≤ M, where e_{nm} is the relative expression level of the nth gene in the mth sample as measured by the mth array.§ The vector in the nth row of the matrix ê, 〈g_{n} ≡ 〈nê, lists the relative expression of the nth gene across the different samples which correspond to the different arrays. The vector in the mth column of the matrix ê, a_{m}〉 ≡ êm〉, lists the genomewide relative expression measured by the mth array.
SVD (10) is then linear transformation of the expression data from the Ngenes × Marrays space to the reduced L“eigenarrays” × L“eigengenes” space, where L = min{M, N} (see Fig. 7 in supplemental material at www.pnas.org), 1 In this space the data are represented by the diagonal nonnegative matrix ê, of size Leigengenes × Leigenarrays, which satisfies 〈kɛ̂l〉 ≡ ɛ_{l}δ_{kl} ≥ 0 for all 1 ≤ k,l ≤ L, such that the lth eigengene is expressed only in the corresponding lth eigenarray, with the corresponding “eigenexpression” level ɛ_{l}. Therefore, the expression of each eigengene (or eigenarray) is decoupled from that of all other eigengenes (or eigenarrays). The “fraction of eigenexpression,” 2 indicates the relative significance of the lth eigengene and eigenarray in terms of the fraction of the overall expression that they capture. Assume also that the eigenexpression levels are arranged in decreasing order of significance, such that ɛ_{1} ≥ ɛ_{2} ≥ … ≥ ɛ_{L} ≥ 0. “Shannon entropy” of a dataset, 3 measures the complexity of the data from the distribution of the overall expression between the different eigengenes (and eigenarrays), where d = 0 corresponds to an ordered and redundant dataset in which all expression is captured by a single eigengene (and eigenarray), and d = 1 corresponds to a disordered and random dataset where all eigengenes (and eigenarrays) are equally expressed.
The transformation matrices û and v̂^{T} define the Ngenes × Leigenarrays and the Leigengenes × Marrays basis sets, respectively. The vector in the lth row of the matrix v̂^{T}, 〈γ_{l} ≡ 〈lv̂^{T}, lists the expression of the lth eigengene across the different arrays. The vector in the lth column of the matrix û, α_{l}〉 ≡ ûl〉, lists the genomewide expression in the lth eigenarray. The eigengenes and eigenarrays are orthonormal superpositions of the genes and arrays, such that the transformation matrices û and v̂ are both orthogonal 4 where Î is the identity matrix. Therefore, the expression of each eigengene (or eigenarray) is not only decoupled but also decorrelated from that of all other eigengenes (or eigenarrays). The eigengenes and eigenarrays are unique, except in degenerate subspaces, defined by subsets of equal eigenexpression levels, and except for a phase factor of ±1, such that each eigengene (or eigenarray) captures both parallel and antiparallel gene (or array) expression patterns. Therefore, SVD is datadriven, except in degenerate subspaces.
SVD Calculation.
According to Eqs. 1 and 4, the Marrays × Marrays symmetric correlation matrix â = ê^{T}ê = v̂ɛ̂^{2}v̂^{T} is represented in the Leigengenes × Leigengenes space by the diagonal matrix ɛ̂^{2}. The Ngenes × Ngenes correlation matrix ĝ = êê^{T} = ûɛ̂^{2}û^{T} is represented in the Leigenarrays × Leigenarrays space also by ɛ̂^{2}, where for L = min{M, N} = M, ĝ has a null subspace of at least N − M null eigenvalues. We, therefore, calculate the SVD of a dataset ê, with M ≪ N, by diagonalizing â, and then projecting the resulting v̂ and ɛ̂ onto ê to obtain û = êv̂ɛ̂^{−1}.
Pattern Inference.
The decorrelation of the eigengenes (and eigenarrays) suggests the possibility that some of the eigengenes (and the corresponding eigenarrays) represent independent regulatory programs or processes (and corresponding cellular states). We infer that an eigengene γ_{l}〉 represents a regulatory program or process from its expression pattern across all arrays, when this pattern is biologically interpretable. This inference may be supported by a corresponding coherent biological theme reflected in the functions of the genes, whose expression patterns correlate or anticorrelate with the pattern of this eigengene. With this we assume that the corresponding eigenarray α_{l}〉 (which lists the amplitude of this eigengene pattern in the expression of each gene g_{n}〉 relative to all other genes 〈nα_{l}〉 = 〈g_{n}γ_{l}〉/ɛ_{l}) represents the cellular state which corresponds to this process. We infer that the eigenarray α_{l}〉 represents a cellular state from the arrays whose expression patterns correlate or anticorrelate with the pattern of this eigenarray. Upon sorting of the genes, this inference may be supported by the expression pattern of this eigenarray across all genes, when this pattern is biologically interpretable.
Data Normalization.
The decoupling of the eigengenes and eigenarrays allows filtering the data without eliminating genes or arrays from the dataset. We filter any of the eigengenes γ_{l}〉 (and the corresponding eigenarray α_{l}〉) ê → ê − ɛ_{l}α_{l}〉 〈γ_{l}, by substituting zero for the eigenexpression level ɛ_{l} = 0 in the diagonal matrix ɛ̂ and reconstructing the data according to Eq. 1. We normalize the data by filtering out those eigengenes (and eigenarrays) that are inferred to represent noise or experimental artifacts.
Degenerate Subspace Rotation.
The uniqueness of the eigengenes and eigenarrays does not hold in a degenerate subspace, defined by equal eigenexpression levels. We approximate significant similar eigenexpression levels ɛ_{l} ≈ ɛ_{l+1} ≈ … ≈ ɛ_{m} with ɛ_{l} = … = ɛ_{m} = . Therefore, Eqs. 1–4 remain valid upon rotation of the corresponding eigengenes {(γ_{l}〉, … , γ_{m}〉) → R̂(γ_{l}〉, … , γ_{m}〉)}, and eigenarrays {(α_{l}〉, … , α_{m}〉) → R̂(α_{l}〉, … , α_{m}〉)}, for all orthogonal R̂, R̂^{T}R̂ = Î. We choose a unique rotation R̂ by subjecting the rotated eigengenes to m − l constraints, such that these constrained eigengenes may be advantageous in interpreting and presenting the expression data.
Data Sorting.
Inferring that eigengenes (and eigenarrays) represent independent processes (and cellular states) allows sorting the data by similarity in the expression of any chosen subset of these eigengenes (and eigenarrays), rather than by overall similarity in expression. Given two eigengenes γ_{k}〉 and γ_{l}〉 (or eigenarrays α_{k}〉 and α_{l}〉), we plot the correlation of γ_{k}〉 with each gene g_{n}〉, 〈γ_{k}g_{n}〉/〈g_{n}g_{n}〉 (or α_{k}〉 with each array a_{m}〉) along the yaxis, vs. that of γ_{l}〉 (or α_{l}〉) along the xaxis. In this plot, the distance of each gene (or array) from the origin is its amplitude of expression in the subspace spanned by γ_{k}〉 and γ_{l}〉 (or α_{k}〉 and α_{l}〉), relative to its overall expression r_{n} ≡ 〈g_{n}g_{n}〉^{−1} (or r_{m} ≡ 〈a_{m}a_{m}〉^{−1} ). The angular distance of each gene (or array) from the xaxis is its phase in the transition from the expression pattern γ_{l}〉 to γ_{k}〉 and back to γ_{l}〉 (or α_{l}〉 to α_{k}〉 and back to α_{l}〉) tan φ_{n} ≡ 〈γ_{k}g_{n}〉/〈γ_{l}g_{n}〉, (or tan φ_{m} ≡ 〈α_{k}a_{n}〉/〈α_{l}a_{m}〉). We sort the genes (or arrays) according to φ_{n} (or φ_{m}).
Biological Data Analysis: ElutriationSynchronized Cell Cycle
Spellman et al. (3) monitored genomewide mRNA levels, for 6,108 ORFs of the budding yeast Saccharomyces cerevisiae simultaneously, over approximately one cell cycle period, T ≈ 390 min, in a yeast culture synchronized by elutriation, relative to a reference mRNA from an asynchronous yeast culture, at 30min intervals. The elutriation dataset we analyze (see supplemental data and Mathematica notebook at www.pnas.org and at http://genomewww.stanford.edu/SVD/) tabulates the measured ratios of gene expression levels for the N = 5,981 genes, 784 of which were classified by Spellman et al. as cell cycle regulated, with no missing data in the M = 14 arrays.
Pattern Inference.
Consider the 14 eigengenes of the elutriation dataset. The first and most significant eigengene γ_{1}〉, which describes time invariant relative expression during the cell cycle (Fig. 8a at www.pnas.org), captures more than 90% of the overall relative expression in this experiment (Fig. 8b). The entropy of the dataset, therefore, is low d = 0.14 ≪ 1. This suggests that the underlying processes are manifested by weak perturbations of a steady state of expression. This also suggests that timeinvariant additive constants due to uncontrolled experimental variables may be superimposed on the data. We infer that γ_{1}〉 represents experimental additive constants superimposed on a steady gene expression state, and assume that α_{1}〉 represents the corresponding steady cellular state. The second, third, and fourth eigengenes, which show oscillations during the cell cycle (Fig. 8c), capture about 3%, 1%, and 0.5% of the overall relative expression, respectively. The time variation of γ_{3}〉 fits a normalized sine function of period T, sin(2πt/T). We infer that γ_{3}〉 represents expression oscillation, which is consistent with gene expression oscillations during a cell cycle. The time variations of the second and fourth eigengenes fit a cosine function of period T with the amplitude of a normalized cosine with this period, cos 2πt/T. However, while γ_{2}〉 shows decreasing expression on transition from t = 0 to 30 min, γ_{4}〉 shows increasing expression. We infer that γ_{2}〉 and γ_{4}〉 represent initial transient increase and decrease in expression in response to the elutriation, respectively, superimposed on expression oscillation during the cell cycle.
Data Normalization.
We filter out the first eigengene and eigenarray of the elutriation dataset, ê → ê_{C} = ê − ɛ_{1}α_{1}〉 〈γ_{1}, removing the steady state of expression. Each of the elements of the dataset ê_{C}, 〈nê_{C}m〉 ≡ e_{C,nm}, is the difference of the measured expression of the nth gene in the mth array from the steadystate levels of expression for these gene and array as calculated by SVD. Therefore, e_{C,nm}^{2} is the variance in the measured expression of the nth gene in the mth array. Let ê_{LV} tabulate the natural logarithm of the variances in elutriation expression, such that each element of ê_{LV} satisfies 〈nê_{LV}m〉 ≡ log(e_{C,nm}^{2}) for all 1 ≤ n ≤ N and 1 ≤ m ≤ M, and consider the eigengenes of ê_{LV} (Fig. 9a in supplemental material at www.pnas.org). The first eigengene γ_{1}〉_{LV}, which captures more than 80% of the overall information in this dataset (Fig. 9b), describes a weak initial transient increase superimposed on a timeinvariant scale of expression variance. The initial transient increase in the scale of expression variance may be a response to the elutriation. The timeinvariant scale of expression variance suggests that a steady scale of experimental as well as biological uncertainty is associated with the expression data. This also suggests that timeinvariant multiplicative constants due to uncontrolled experimental variables may be superimposed on the data. We filter out γ_{1}〉_{LV}, removing the steady scale of expression variance, ê_{LV} → ê_{CLV} = ê_{LV} − ɛ_{1,LV}α_{1}〉_{LV LV}〈γ_{1}.
The normalized elutriation dataset ê_{N}, where each of its elements satisfies 〈nê_{N}m〉 ≡ sign(e_{C,nm}), tabulates for each gene and array expression patterns that are approximately centered at the steadystate expression level (i.e., of approximately zero arithmetic means), with variances which are approximately normalized by the steady scale of expression variance (i.e., of approximately unit geometric means). The first and second eigengenes, γ_{1}〉_{N} and γ_{2}〉_{N}, of ê_{N} (Fig. 1a), which are of similar significance, capture together more than 40% of the overall normalized expression (Fig. 1b). The time variations of γ_{1}〉_{N} and γ_{2}〉_{N} fit normalized sine and cosine functions of period T and initial phase θ ≈ 2π/13, sin(2πt/T − θ) and cos(2πt/T − θ), respectively (Fig. 1c). We infer that γ_{1}〉_{N} and γ_{2}〉_{N} represent cell cycle expression oscillations, and assume that the corresponding eigenarrays α_{1}〉_{N} and α_{2}〉_{N} represent the corresponding cell cycle cellular states. Upon sorting of the genes (and arrays) according to γ_{1}〉_{N} and γ_{2}〉_{N} (and α_{1}〉_{N} and α_{2}〉_{N}), the initial phase θ ≈ 2π/13 can be interpreted as a delay of 30 min between the start of the experiment and that of the cell cycle stage G_{1}. The decay to zero in the time variation of γ_{2}〉_{N} at t = 360 and 390 min can be interpreted as dephasing in time of the initially synchronized yeast culture.
Data Sorting.
Consider the normalized expression of the 14 elutriation arrays {a_{m}〉} in the subspace spanned by α_{1}〉_{N} and α_{2}〉_{N}, which is assumed to approximately represent all cell cycle cellular states (Fig. 2a). All arrays have at least 25% of their normalized expression in this subspace, with their distances from the origin satisfying 0.5 ≤ r_{m} < 1, except for the eleventh array a_{11}〉. This suggests that α_{1}〉_{N} and α_{2}〉_{N} are sufficient to approximate the elutriation array expression. The sorting of the arrays according to their phases {φ_{m}}, which describes the transition from the expression pattern α_{2}〉_{N} to α_{1}〉_{N} and back to α_{2}〉_{N}, gives an array order which is similar to that of the cell cycle time points measured by the arrays, an order that describes the progress of the cell cycle expression from the M/G_{1} stage through G_{1}, S, S/G_{2}, and G_{2}/M and back to M/G_{1}.
Because α_{1}〉_{N} is correlated with the arrays a_{4}〉, a_{5}〉, a_{6}〉, and a_{7}〉 and is anticorrelated with a_{13}〉 and a_{14}〉, we associate α_{1}〉_{N} with the cell cycle cellular state of transition from G_{1} to S, and −α_{1}〉_{N} with the transition from G_{2}/M to M/G_{1}. Similarly, α_{2}〉_{N} is correlated with a_{2}〉 and a_{3}〉, and therefore we associate α_{2}〉_{N} with the transition from M/G_{1} to G_{1}. Also, α_{2}〉_{N} is anticorrelated with a_{8}〉 and a_{10}〉, and therefore we associate −α_{2}〉_{N} with the transition from S to S/G_{2}. With these associations the phase of a_{1}〉, φ_{1} = −θ ≈ −2π/13, corresponds to the 30min delay between the start of the experiment and that of the cell cycle stage G_{1}, which is also present in the inferred cell cycle expression oscillations γ_{1}〉_{N} and γ_{2}〉_{N}.
Consider also the expression of the 5,981 genes {g_{n}〉} in the subspace spanned by γ_{1}〉_{N} and γ_{2}〉_{N}, which is inferred to approximately represent all cell cycle expression oscillations (Fig. 10 in supplemental material at www.pnas.org). One may expect that genes that have almost all of their normalized expression in this subspace with r_{n} ≈ 1 are cell cycle regulated, and that genes that have almost no expression in this subspace with r_{n} ≈ 0, are not regulated by the cell cycle at all. Indeed, of the 784 genes that were classified by Spellman et al. (3) as cell cycle regulated, 641 have more than 25% of their normalized expression in this subspace (Fig. 2b). We sort all 5,981 genes according to their phases {φ_{n}}, to describe the transition from the expression pattern γ_{2}〉_{N} to that of γ_{1}〉_{N} and back to γ_{2}〉_{N}, starting at φ_{1} ≈ −2π/13. One may expect this to order the genes according to the stages in the cell cycle in which their expression patterns peak. However, for the 784 cell cycle regulated genes this sorting gives a classification of the genes into the five cell cycle stages, which is somewhat different than the classification by Spellman et al. This may be due to the poor quality of the elutriation expression data, as synchronization by elutriation was not very effective in this experiment. For the α factorsynchronized cell cycle expression there is much better agreement between the two classifications (Fig. 5b).
With all 5,981 genes sorted, the gene variations of α_{1}〉_{N} and α_{2}〉_{N} fit normalized sine and cosine functions of period Z ≡ N − 1 = 5,980 and initial phase θ ≈ 2π/13, − sin(2πz/Z − θ) and cos(2πz/Z − θ), respectively, where z ≡ n − 1 (Fig. 3 b and c). The sorted and normalized elutriation expression fit approximately a traveling wave of expression, varying sinusoidally across both genes and arrays, such that the expression of the nth gene in the mth array satisfies 〈nê_{N}m〉 ∝ −2 cos[2π(t/T − z/Z)]/ (Fig. 3a).
Biological Data Analysis: α FactorSynchronized Cell Cycle and CLB2 and CLN3 Overactivations
Spellman et al. (3) also monitored genomewide mRNA levels, for 6,108 yeast ORFs simultaneously, over approximately two cell cycle periods, in a yeast culture synchronized by α factor, relative to a reference mRNA from an asynchronous yeast culture, at 7min intervals for 119 min. They also measured, in two independent experiments, mRNA levels of yeast strain cultures with overactivated CLB2, which encodes a G_{2}/M cyclin, both at t = 40 min relative to their levels at the start of overactivation at t = 0. Two additional independent experiments measured mRNA levels of strain cultures with overactivated CLN3, which encodes a G_{1}/S cyclin, at t = 30 and 40 min relative to their levels at the start of overactivation at t = 0. The dataset for the α factor, CLB2, and CLN3 experiments we analyze (see supplemental data and Mathematica notebook at www.pnas.org) tabulates the ratios of gene expression levels for the N = 4,579 genes, 638 of which were classified by Spellman et al. as cell cycle regulated, with no missing data in the M = 22 arrays.
After data normalization and degenerate subspace rotation (see Appendix in supplemental material at www.pnas.org), the time variations of γ_{1}〉_{RN} and γ_{2}〉_{RN} fit normalized sine and cosine functions of two 66min periods during the cell cycle, from t = 7 to 119 min, and initial phase θ ≈ π/4, respectively (Fig. 4c). While γ_{2}〉_{RN} describes steadystate expression in the CLB2 and CLN3overactive arrays, γ_{1}〉_{RN} describes underexpression in the CLB2overactive arrays and overexpression in the CLN3overactive arrays.
Upon sorting the 4,579 genes in the subspace spanned by γ_{1}〉_{RN} and γ_{2}〉_{RN} (Fig. 5b), γ_{1}〉_{RN} is correlated with genes that peak late in the cell cycle stage G_{1} and early in S, among them CLN3, and we associate γ_{1}〉_{RN} with the cell cycle expression oscillations that start at the transition from G_{1} to S and are dependent on CLN3, which encodes a G_{1}/S cyclin. Also, γ_{1}〉_{RN} is anticorrelated with genes that peak late in G_{2}/M and early in M/G_{1}, among them CLB2, and therefore we associate −γ_{1}〉_{RN} with the oscillations that start at the transition from G_{2}/M to M/G_{1} and are dependent on CLB2, which encodes a G_{2}/M cyclin. Similarly, γ_{2}〉_{RN} is correlated with genes that peak late in M/G_{1} and early in G_{1}, anticorrelated with genes that peak late in S and early in S/G_{2}, and uncorrelated with CLB2 and CLN3. We, therefore, associate γ_{2}〉_{RN} with the oscillations that start at the transition from M/G_{1} to G_{1} (and appear to be CLB2 and CLN3independent), and −γ_{2}〉_{RN} with the oscillations that start at the transition from S to S/G_{2} (and appear to be CLB2 and CLN3independent).
Upon sorting the 22 arrays in the subspace spanned by α_{1}〉_{RN} and α_{2}〉_{RN} (Fig. 5a), α_{1}〉_{RN} is correlated with the arrays a_{13}〉 and a_{14}〉, as well as with a_{21}〉 and a_{22}〉, which measure the CLN3overactive samples. We therefore associate α_{1}〉_{RN} with the cell cycle cellular state of transition from G_{1} to S, which is simulated by CLN3 overactivation. Also, α_{1}〉_{RN} is anticorrelated with the arrays a_{9}〉 and a_{10}〉, as well as with a_{19}〉 and a_{20}〉, which measure the CLB2overactive samples. We associate −α_{1}〉_{RN} with the cellular transition from G_{2}/M to M/G_{1}, which is simulated by CLB2 overactivation. Similarly, α_{2}〉_{RN} appears to be correlated with a_{2}〉, a_{3}〉, a_{11}〉, and a_{12}〉, anticorrelated with a_{6}〉, a_{7}〉, a_{16}〉, and a_{17}〉, and uncorrelated with a_{19}〉, a_{20}〉, a_{21}〉, or a_{22}〉. We therefore associate α_{2}〉_{RN} with the cellular transition from M/G_{1} to G_{1} (which appears to be CLB2 and CLN3independent), and −α_{2}〉_{RN} with the cellular transition from S to S/G_{2} (which also appears to be CLB2 and CLN3independent).
With all 4,579 genes sorted the gene variations of α_{1}〉_{RN} and α_{2}〉_{RN} fit normalized sine and cosine functions of period Z ≡ N − 1 = 4,578 and initial phase π/8, respectively (Fig. 6 b and c). The normalized and sorted cell cycle expression approximately fits a traveling wave, varying sinusoidally across both genes and arrays. The normalized and sorted expression in the CLB2 and CLN3overactive arrays approximately fits standing waves, constant across the arrays and varying sinusoidally across genes only, which appear similar to −α_{1}〉_{RN} and α_{1}〉_{RN}, respectively (Fig. 6a).
Conclusions
We have shown that SVD provides a useful mathematical framework for processing and modeling genomewide expression data, in which both the mathematical variables and operations may be assigned biological meaning.
Acknowledgments
We thank S. Kim for insightful discussions, G. Sherlock for technical assistance and careful reading, and J. Doyle and P. Green for thoughtful reviews of this manuscript. This work was supported by a grant from the National Cancer Institute (National Institutes of Health, CA77097). O.A. is an Alfred P. Sloan and U.S. Department of Energy Postdoctoral Fellow in Computational Molecular Biology, and a National Human Genome Research Institute Individual Mentored Research Scientist Development Awardee in Genomic Research and Analysis (National Institutes of Health, 1 K01 HG0003801). P.O.B. is an Associate Investigator of the Howard Hughes Medical Institute.
Footnotes
Abbreviation
 SVD,
 singular value decomposition
 Accepted June 15, 2000.
 Copyright © 2000, The National Academy of Sciences
References
 ↵
 ↵
 Schena M,
 Shalon D,
 Davis R W,
 Brown P O
 ↵
 Spellman P T,
 Sherlock G,
 Zhang M Q,
 Iyer V R,
 Anders K,
 Eisen M B,
 Brown P O,
 Botstein D,
 Futcher B
 ↵

 Eisen M B,
 Spellman P T,
 Brown P O,
 Botstein D

 Alon U,
 Barkai N,
 Notterman D A,
 Gish K,
 Ybarra S,
 Mack D,
 Levine A J

 Tamayo P,
 Slonim D,
 Mesirov J,
 Zhu Q,
 Kitareewan S,
 Dmitrovsky E,
 Lander E S,
 Golub T R
 ↵
 Brown M P S,
 Grundy W N,
 Lin D,
 Cristianini N,
 Sugnet C W,
 Furey T S,
 Ares M Jr,
 Haussler D
 ↵
 Golub G H,
 Van Loan C F
 ↵
 Mallat S G
 ↵
 Anderson T W