Singular value decomposition for genome-wide expression data processing and modeling

SVD Calculation.

According to Eqs. 1 and 4, the M-arrays × M-arrays symmetric correlation matrix â = ê^Tê = v̂ɛ̂²v̂^T 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 v̂ and ɛ̂ onto ê to obtain û = êv̂ɛ̂⁻¹.

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̂^TR̂ = Î. 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 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).

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(2πt/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 2πt/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 = ê − ɛ₁|α₁〉 〈γ₁|, 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 |γ₁〉_LV, 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 |γ₁〉_LV, removing the steady scale of expression variance, ê_LV → ê_CLV = ê_LV − ɛ_1,LV|α₁〉_LV LV〈γ₁|.

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 steady-state 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, |γ₁〉_N and |γ₂〉_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 |γ₁〉_N and |γ₂〉_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 |γ₁〉_N and |γ₂〉_N represent cell cycle expression oscillations, and assume that the corresponding eigenarrays |α₁〉_N and |α₂〉_N represent the corresponding cell cycle cellular states. Upon sorting of the genes (and arrays) according to |γ₁〉_N and |γ₂〉_N (and |α₁〉_N and |α₂〉_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₁. The decay to zero in the time variation of |γ₂〉_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 |α₁〉_N and |α₂〉_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₁₁〉. This suggests that |α₁〉_N and |α₂〉_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 |α₂〉_N to |α₁〉_N and back to |α₂〉_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₁ stage through G₁, S, S/G₂, and G₂/M and back to M/G₁.

Because |α₁〉_N is correlated with the arrays |a₄〉, |a₅〉, |a₆〉, and |a₇〉 and is anticorrelated with |a₁₃〉 and |a₁₄〉, we associate |α₁〉_N with the cell cycle cellular state of transition from G₁ to S, and −|α₁〉_N with the transition from G₂/M to M/G₁. Similarly, |α₂〉_N is correlated with |a₂〉 and |a₃〉, and therefore we associate |α₂〉_N with the transition from M/G₁ to G₁. Also, |α₂〉_N is anticorrelated with |a₈〉 and |a₁₀〉, and therefore we associate −|α₂〉_N with the transition from S to S/G₂. With these associations the phase of |a₁〉, φ₁ = −θ ≈ −2π/13, corresponds to the 30-min delay between the start of the experiment and that of the cell cycle stage G₁, which is also present in the inferred cell cycle expression oscillations |γ₁〉_N and |γ₂〉_N.

Consider also the expression of the 5,981 genes {|g_n〉} in the subspace spanned by |γ₁〉_N and |γ₂〉_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 |γ₂〉_N to that of |γ₁〉_N and back to |γ₂〉_N, starting at φ₁ ≈ −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 α factor-synchronized cell cycle expression there is much better agreement between the two classifications (Fig. 5b).

With all 5,981 genes sorted, the gene variations of |α₁〉_N and |α₂〉_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).