Introduction

Top Abstract Introduction Mathematical Framework:... Biological Data Analysis:... Biological Data Analysis: ... Conclusions References

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 genome-wide expression data. SVD is also known as Karhunen-Loève expansion in pattern recognition (11) and as principal-component 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 data-driven, 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 genome-wide 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 genome-wide expression data, in which both the mathematical variables and operations may be assigned biological meaning.

	Mathematical Framework: Singular Value Decomposition

Top Abstract Introduction Mathematical Framework:... Biological Data Analysis:... Biological Data Analysis: ... Conclusions References

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 genome-wide expression levels in M different samplesi.e., under M different experimental conditions. Let the matrix ê, of size N-genes × M-arrays, 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 genome-wide relative expression measured by the mth array.

SVD (10) is then linear transformation of the expression data from the N-genes × M-arrays 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 L-eigengenes × L-eigenarrays, which satisfies k||l  _lkl  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 ^T define the N-genes × L-eigenarrays and the L-eigengenes × M-arrays basis sets, respectively. The vector in the lth row of the matrix ^T, _l|  l|^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 genome-wide expression in the lth eigenarray. The eigengenes and eigenarrays are orthonormal superpositions of the genes and arrays, such that the transformation matrices û and 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 data-driven, except in degenerate subspaces.

SVD Calculation. According to Eqs. 1 and 4, the M-arrays × M-arrays symmetric correlation matrix â = ê^Tê = ^2T is represented in the L-eigengenes × L-eigengenes space by the diagonal matrix ². The N-genes × N-genes correlation matrix  = êê^T = û²û^T is represented in the L-eigenarrays × L-eigenarrays space also by ², 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 and onto ê to obtain û = ê¹.

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)  (|_l, ... , |_m)}, and eigenarrays {(|_l, ... , |_m)  (|_l, ... , |_m)}, for all orthogonal , ^T = Î. We choose a unique rotation 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 y-axis, vs. that of |_l (or |_l) along the x-axis. 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¹  (or r_m  a_m|a_m¹ ). The angular distance of each gene (or array) from the x-axis 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: Elutriation-Synchronized Cell Cycle

Top Abstract Introduction Mathematical Framework:... Biological Data Analysis:... Biological Data Analysis: ... Conclusions References

Spellman et al. (3) monitored genome-wide 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 30-min intervals. The elutriation dataset we analyze (see supplemental data and Mathematica notebook at www.pnas.org and at http://genome-www.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 |₁, 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 time-invariant additive constants due to uncontrolled experimental variables may be superimposed on the data. We infer that |₁ represents experimental additive constants superimposed on a steady gene expression state, and assume that |₁ 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 |₃ fits a normalized sine function of period T, sin(2t/T). We infer that |₃ 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 2t/T. However, while |₂ shows decreasing expression on transition from t = 0 to 30 min, |₄ shows increasing expression. We infer that |₂ and |₄ 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|, 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 steady-state levels of expression for these gene and array as calculated by SVD. Therefore, e_C,nm² 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²) 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 |_1LV, which captures more than 80% of the overall information in this dataset (Fig. 9b), describes a weak initial transient increase superimposed on a time-invariant scale of expression variance. The initial transient increase in the scale of expression variance may be a response to the elutriation. The time-invariant 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 time-invariant multiplicative constants due to uncontrolled experimental variables may be superimposed on the data. We filter out |_1LV, removing the steady scale of expression variance, ê_LV  ê_CLV = ê_LV  1,LV|_1LV LV1|.

View larger version (35K): [in this window] [in a new window]   Fig. 1.   Normalized elutriation eigengenes. (a) Raster display of NT, the expression of 14 eigengenes in 14 arrays. (b) Bar chart of the fractions of eigenexpression, showing that |1N and |2N capture about 20% of the overall normalized expression each, and a high entropy d = 0.88. (c) Line-joined graphs of the expression levels of |1N (red) and |2N (blue) in the 14 arrays fit dashed graphs of normalized sine (red) and cosine (blue) of period T = 390 min and phase  = 2/13, respectively.

Data Sorting. Consider the normalized expression of the 14 elutriation arrays {|a_m} in the subspace spanned by |_1N and |_2N, 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₁₁. This suggests that |_1N and |_2N 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 |_2N to |_1N and back to |_2N, 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₁ stage through G₁, S, S/G₂, and G₂/M and back to M/G₁.

View larger version (40K): [in this window] [in a new window]   Fig. 2.   Normalized elutriation expression in the subspace associated with the cell cycle. (a) Array correlation with |1N along the y-axis vs. that with |2N along the x-axis, color-coded according to the classification of the arrays into the five cell cycle stages, M/G1 (yellow), G1 (green), S (blue), S/G2 (red), and G2/M (orange). The dashed unit and half-unit circles outline 100% and 25% of overall normalized array expression in the |1N and |2N subspace. (b) Correlation of each gene with |1N vs. that with |2N, for 784 cell cycle regulated genes, color-coded according to the classification by Spellman et al. (3).

View larger version (77K): [in this window] [in a new window]   Fig. 3.   Genes sorted by relative correlation with |1N and |2N of normalized elutriation. (a) Normalized elutriation expression of the sorted 5,981 genes in the 14 arrays, showing traveling wave of expression. (b) Eigenarrays expression; the expression of |1N and |2N, the eigenarrays corresponding to |1N and |2N, displays the sorting. (c) Expression levels of |1N (red) and |2N (green) fit normalized sine and cosine functions of period Z  N  1 = 5,980 and phase   2/13 (blue), respectively.

	Biological Data Analysis:  Factor-Synchronized Cell Cycle and CLB2 and CLN3 Overactivations

Top Abstract Introduction Mathematical Framework:... Biological Data Analysis:... Biological Data Analysis: ... Conclusions References

Spellman et al. (3) also monitored genome-wide 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 7-min intervals for 119 min. They also measured, in two independent experiments, mRNA levels of yeast strain cultures with overactivated CLB2, which encodes a G₂/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₁/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 |_1RN and |_2RN fit normalized sine and cosine functions of two 66-min periods during the cell cycle, from t = 7 to 119 min, and initial phase   /4, respectively (Fig. 4c). While |_2RN describes steady-state expression in the CLB2- and CLN3-overactive arrays, |_1RN describes underexpression in the CLB2-overactive arrays and overexpression in the CLN3-overactive arrays.

View larger version (43K): [in this window] [in a new window]   Fig. 4.   Rotated normalized  factor, CLB2, and CLN3 eigengenes. (a) Raster display of RNT, where |1RN = 21|1N, |2RN = 1|2N, and |3RN = 2|3N. (b) |1RN, |2RN and |3RN capture 20% of the overall normalized expression each. (c) Expression levels of |1RN (red) and |2RN (blue) fit dashed graphs of normalized sine (red) and cosine (blue) of period T/2 = 66 min and phase /4, respectively, and |3RN (green) fits dashed graph of normalized sine of period T = 112 min and phase /8, from t = 7 to t = 119 min during the cell cycle.

Upon sorting the 4,579 genes in the subspace spanned by |_1RN and |_2RN (Fig. 5b), |_1RN is correlated with genes that peak late in the cell cycle stage G₁ and early in S, among them CLN3, and we associate |_1RN with the cell cycle expression oscillations that start at the transition from G₁ to S and are dependent on CLN3, which encodes a G₁/S cyclin. Also, |_1RN is anticorrelated with genes that peak late in G₂/M and early in M/G₁, among them CLB2, and therefore we associate |_1RN with the oscillations that start at the transition from G₂/M to M/G₁ and are dependent on CLB2, which encodes a G₂/M cyclin. Similarly, |_2RN is correlated with genes that peak late in M/G₁ and early in G₁, anticorrelated with genes that peak late in S and early in S/G₂, and uncorrelated with CLB2 and CLN3. We, therefore, associate |_2RN with the oscillations that start at the transition from M/G₁ to G₁ (and appear to be CLB2- and CLN3-independent), and |_2RN with the oscillations that start at the transition from S to S/G₂ (and appear to be CLB2- and CLN3-independent).

View larger version (41K): [in this window] [in a new window]   Fig. 5.   Rotated normalized  factor, CLB2, and CLN3 expression in the subspace associated with the cell cycle. (a) Array correlation with |1RN along the y-axis vs. that with |2RN along the x-axis, color-coded according to the classification of the arrays into the five cell cycle stages, M/G1 (yellow), G1 (green), S (blue), S/G2 (red), and G2/M (orange). The dashed unit and half-unit circles outline 100% and 25% of overall normalized array expression in the |1RN and |2RN subspace. (b) Correlation of each gene with |1RN vs. that with |2RN, for 638 cell cycle regulated genes, color-coded according to the classification by Spellman et al. (3).

Upon sorting the 22 arrays in the subspace spanned by |_1RN and |_2RN (Fig. 5a), |_1RN is correlated with the arrays |a₁₃ and |a₁₄, as well as with |a₂₁ and |a₂₂, which measure the CLN3-overactive samples. We therefore associate |_1RN with the cell cycle cellular state of transition from G₁ to S, which is simulated by CLN3 overactivation. Also, |_1RN is anticorrelated with the arrays |a₉ and |a₁₀, as well as with |a₁₉ and |a₂₀, which measure the CLB2-overactive samples. We associate |_1RN with the cellular transition from G₂/M to M/G₁, which is simulated by CLB2 overactivation. Similarly, |_2RN appears to be correlated with |a₂, |a₃, |a₁₁, and |a₁₂, anticorrelated with |a₆, |a₇, |a₁₆, and |a₁₇, and uncorrelated with |a₁₉, |a₂₀, |a₂₁, or |a₂₂. We therefore associate |_2RN with the cellular transition from M/G₁ to G₁ (which appears to be CLB2- and CLN3-independent), and |_2RN with the cellular transition from S to S/G₂ (which also appears to be CLB2- and CLN3-independent).

With all 4,579 genes sorted the gene variations of |_1RN and |_2RN 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 CLN3-overactive arrays approximately fits standing waves, constant across the arrays and varying sinusoidally across genes only, which appear similar to |_1RN and |_1RN, respectively (Fig. 6a).

View larger version (86K): [in this window] [in a new window]   Fig. 6.   Genes sorted by relative correlation with |1RN and |2RN of rotated normalized  factor, CLB2, and CLN3. (a) Normalized expression of the sorted 4,579 genes in the 22 arrays, showing traveling wave of expression from t = 0 to 119 min during the cell cycle and standing waves of expression in the CLB2- and CLN3-overactive arrays. (b) Eigenarrays expression; the expression of |1RN and |2RN, the eigenarrays corresponding to |1RN and |2RN, displays the sorting. (c) Expression levels of |1RN (red) and |2RN (green) fit normalized sine and cosine functions of period Z  N  1 = 4,578 and phase /8 (blue), respectively.

	Conclusions

Top Abstract Introduction Mathematical Framework:... Biological Data Analysis:... Biological Data Analysis: ... Conclusions