Cluster analysis of gene expression dynamics

Autoregressive Models.

Let S = {y_j1, … y_jt, … y_jn} be a stationary time series of continuous values. The series follows an autoregressive model of order p, say ar(p), if the value of the series at time t > p is a linear function of the values observed in the previous p steps. More formally, we can describe the model in matrix form as

1

where y_j is the vector (y_j(p+1), … , y_jn)^T, X_j is the (n − p) × q regression matrix whose tth row is (1, y_j(t−1), … , y_j(t−p)), for t > p, and q = p + 1. The elements of the vector β_j = {β_j0, β_j1, … , β_jp} are the autoregressive coefficients, and ɛ_j = (ɛ_j(p+1), … , ɛ_jn)^T is a vector of uncorrelated errors that we assume normally distributed, with expected value E(ɛ_jt) = 0 and variance V(ɛ_jt) = σ, for any t. The value p is the autoregressive order and specifies that, at each time point t, y_jt is independent of the past history given the previous p steps. The time series is stationary if it is invariant by temporal translations. Formally, stationarity requires that the coefficients β_j are such that the roots of the polynomial f(u) = 1 − ∑ β_jhu^h have moduli greater than unity. The model in Eq. 1 represents the evolution of the process around its mean μ, which is related to the β_j coefficients by the equation μ = β_j0/(1 − ∑ β_jh). In particular, μ is well defined as long as ∑ β_j ≠ 1. When the autoregressive order p = 0, the series S becomes a sample of independent observations from a normal distribution with mean μ = β_j0 and variance σ. Note that the model in Eq. 1 is a special case of a state-space model (12).

Probabilistic Scoring Metric.

We describe a set of c clusters of time series as a statistical model M_c, consisting of c autoregressive models with coefficients β_k and variance σ. Each cluster C_k groups m_k time series that are jointly modeled as

where the vector y_k and the matrix X_k are defined by stacking the m_k vectors y_kj and regression matrices X_kj, one for each time series, as follows

Note that we now label the vectors y_j assigned to the same cluster C_k with the double subscript kj, and k denotes the cluster membership. The vector ɛ_k is the vector of uncorrelated errors with zero expected value and constant variance σ. Given a set of possible clustering models, the task is to rank them according to their posterior probability. The posterior probability of each clustering model M_c is

where P(M_c) is the prior probability of M_c, y consists of the data {y_k}, and the quantity f(y|M_c) is the marginal likelihood. The marginal likelihood f(y|M_c) is the solution of the integral

where θ is the vector of parameters specifying the clustering model M_c, and f(θ|M_c) is its prior density. Assuming uniform prior distributions on the model parameters and independence of the time series conditional on the cluster membership, f(y|M_c) can be computed as

2

where n_k is the dimension of the vector y_k, and rss_k= y(I_n − X_k(XX_k)⁻¹X)y_k is the residual sum of squares in cluster C_k. When all clustering models are a priori equally likely, the posterior probability P(M_c|y) is proportional to the marginal likelihood f(y|M_c), which becomes our probabilistic scoring metric.

Agglomerative Bayesian Clustering.

The Bayesian approach to the clustering task is to choose the model Mc with maximum posterior probability. As the number of clustering models grows exponentially with the number of time series, we use an agglomerative, finite-horizon search strategy that iteratively merges time series into clusters. The procedure starts by assuming that each of the m observed time series is generated by a different process. Thus, the initial model M_m consists of m clusters, one for each time series, with score f(y|M_m). The next step is the computation of the marginal likelihood of the m(m − 1) models in which two of the m series are merged into one cluster. The model M_m−1 with maximal marginal likelihood is chosen and, if f(y|M_m) ≥ f(y|M_m−1), no merging is accepted and the procedure stops. If f(y|M_m) < f(y|M_m−1), the merging is accepted, a cluster C_k merging the two time series is created, and the procedure is repeated on the new set of m − 1 clusters, consisting of the remaining m − 2 time series and the cluster C_k.

Heuristic Search.

Although the agglomerative strategy makes the search process feasible, the computational effort can still be extremely demanding when the number m of time series is large. To reduce this effort further, we use a heuristic strategy based on a measure of similarity between time series. The intuition behind this strategy is that the merging of two similar time series has better chances of increasing the marginal likelihood of the model. The heuristic search starts by computing the m(m − 1) pairwise similarity measures of the m time series and selects the model M_m−1 in which the two closest time series are merged into one cluster. If f(y|M_m−1) > f(y|M_m), the merging is accepted, the two time series are merged into a single cluster. The average profile of this cluster is computed by averaging the two observed time series, and the procedure is repeated on the new set of m − 1time series, containing the new cluster profile. If this merging is rejected, the procedure is repeated on pairs of time series with decreasing degree of similarity until an acceptable merging is found. If no acceptable merging is found, the procedure stops. Note that the decision of merging two clusters is actually made on the basis of the posterior probability of the model and that the similarity measure is used only to improve efficiency and limit the risk of falling into local maxima.

Several measures can be used to assess the similarity of two time series, both model-free, such as Euclidean distance, correlation and lag-correlation, and model-based, such as Kullback–Leiber distance. Model-free distances are calculated on the raw data. Because the method uses these similarity measures as heuristic tools rather than scoring metrics, we can actually assess the efficiency of each of these measures to drive the search process toward the model with maximum posterior probability. In this respect, the Euclidean distance of two time series S_i = {y_i1, … , y_1n} and S_j = {y_j1, … , y_jn}, computed as

performs best on the short time series of our data set. This finding is consistent with the results of ref. 9, claiming a better overall performance of Euclidean distance in standard hierarchical clustering of gene expression profiles.

Validation.

Standard statistical diagnostics are used as independent assessment measures of the cluster model found by the heuristic search. Once the procedure terminates, the coefficients β_k of the ar(p) model associated with each cluster C_k are estimated as β̂_k = (XX_k)⁻¹Xy_k, while σ̂ =rss_k/(n_k − q) is the estimate of the within-cluster variance σ. The parameter estimates can be used to compute the fitted values for the series in each cluster as ŷ_kj = X_kjβ̂_k, from which we compute the standardized residuals r_kj = (y_kj − ŷ_kj)/σ̂_k. If ar(p) models provide an accurate approximation of the processes generating the time series, the standardized residuals should behave like a random sample from a standard normal distribution. A normal probability plot or the residuals histogram per cluster are used to assess normality. Departures from normality cast doubt on the autoregressive assumption, so that some data transformation, such as a logarithmic transformation, may be needed. Plots of fitted vs. observed values and of fitted values vs. standardized residuals in each cluster provide further diagnostics. To choose the best autoregressive order, we repeat the clustering for p = 0, 1, … , w, for some preset w—by using the same p for every clustering model—and compute a goodness of fit score defined as

where c is the number of clusters, n_k is the size of the vector y_k in C_k, q = p + 1, where p is the autoregressive order, and rss_k is the residual sum of squares of cluster C_k. This score is derived by averaging the log-scores cumulated by the series assigned to each clusters. The derivation of this score is detailed in the technical report available from the Supporting Appendixes, which are published as supporting information on the PNAS web site, www.pnas.org. The resulting score trades off model complexity—measured by the quantity cq + ∑_kn_klog(n_k − q)—with lack of fit—measured by the quantity ∑n_klog(rss_k), and it generalizes the well known Akaike information criterion goodness of fit criterion of ref. 26 to a set of autoregressive models. We then choose the clustering model with the autoregressive order p that maximizes this goodness of fit score.

Display.

As in ref. 6, a colored map is created by displaying the rows of the original data table according to the pairwise merging of the heuristic search. In our case, a set of binary trees (dendrogram), one for each cluster, is appended to the colored map. Each branching node is labeled with the ratio between the marginal likelihood of the merging accepted at that node and the marginal likelihood of the model without this merging. This ratio measures how many times the model accepting the merging is more likely than the model refusing it.

Statistical Analysis.

Euclidean distance gave the best results: it systematically returned clustering models with higher posterior probability than correlation-based distances. The number of clusters found for p = 0, 1, 2, 3 varied between 4 (p = 0, 1) and 3 (p = 2, 3). To choose a clustering model among these four, we used the goodness of fit score described in Methods. The scores for the four models were, for increasing p, 10130.78, 13187.15, 11980.38, and 11031.12, and the model with autoregressive order p = 1 was therefore selected. This model merges the 517 gene time series into four clusters of 3, 216, 293, and 5 time series, with estimates of the autoregressive coefficients and within-cluster variance β̂₁₀ = 0.518; β̂₁₁ = 0.708; σ̂ = 0.606 in cluster 1, β̂₂₀ = 0.136; β̂₂₁ = 0.776; σ̂ = 0.166 in cluster 2, β̂₃₀ = −0.132; β̂₃₁ = 0.722; σ̂ = 0.091 in cluster 3; and β̂₄₀ = −0.661; β̂₄₁ = 0.328; σ̂ = 0.207 in cluster 4. In this model, merging any of these clusters decreases the posterior probability of the clustering model of at least 10.05 times, a strong evidence in favor of their separation (27). Colored maps and dendrogram of the four clusters are displayed in Fig. 2.

The symmetry of the standardized residuals in Fig. 1, together with the lack of any significant patterns in the scatter plot of the fitted values vs. the standardized residuals and the closeness of fitted and observed values suggests that ar(1) models provide a good approximation of the processes generating these time series. This impression is reinforced further by the averages of the fitted time series in each cluster, shown in Fig. 1, which follow closely their respective cluster average profiles.

Comparison with Correlation-Based Clustering.

The most evident difference between the model in Fig. 2 and the model obtained by the original analysis is the number of clusters: our method detects four distinct clusters, characterized by the autoregressive models described above, while hierarchical clustering merges all 517 genes into a single cluster. Across the four clusters, both average profiles and averages of the fitted time series appear to capture different dynamics. Iyer et al. (8) identify, by visual inspection, eight subgroups of genes—labeled A, B, C, … , I, J—from eight large contiguous patches of color. With the exception of a few genes, our cluster 2 merges the subgroups of time series labeled as D, E, F, G, H, I, and J, and cluster 3 merges the majority of time series assigned to subgroups A, B, and C. Interestingly enough, the members of subgroups A, B, and C differ, on average, by one single value. Similarly, groups D and G differ by a single value, as well as F, H, J, and I.

Our cluster 1 collects the temporal patterns of three genes—IL-8, prostaglandin-endoperoxide synthase 2, and IL-6 (IFN-β2). These time series were assigned by ref. 8 to the subgroups F, I, and J, respectively. Cluster 4 collects the time series of five genes—receptor tyrosine kinase-like orphan receptor, TRABID protein, death-associated protein kinase, DKFZP586G1122 protein, and transcription termination factor-like protein. Three of these time series were assigned by ref. 8 to the A and B subgroups. These two smaller clusters—clusters 1 and 4—are particularly noteworthy because they illustrate how our method automatically identifies islands of particular expression profiles. The first of these two clusters merges cytokines involved in the processes of the inflammatory response and chemotaxis and the signal transduction and cell–cell signaling underlying these processes. The cluster includes IL-8, IL-6, and prostaglandin-endoperoxide synthase 2, which catalyzes the rate-limiting step in the formation of inflammatory prostaglandins. The second small cluster includes genes that are known to be involved in the cell-death/apoptosis processes. This cluster includes kinases and several transcription factors reported to be involved in these processes. The cluster includes receptor tyrosine kinase-like orphan receptor 2, TRAF-binding protein domain, and Death-associated protein kinase. The cluster also includes the transcription termination factor-like protein, which plays a central role in the control of rRNA and mRNA synthesis in mammalian mitochondria (28), and DKFZP586G1122 protein, which has unknown function but has strong homology with murine zinc finger protein Hzf expressed in hematopoiesis.

The number of clusters found by our algorithm is directly inferred from the data, which also provide evidence in favor of a temporal dependency of the observations: the goodness of fit score of the ar(0) clustering model, where the observations are assumed to be marginally independent, is lower than the goodness of fit score of the ar(1) clustering model, which assumes that each observation depends on its immediate predecessor. The allocation of the 20 repeated genes in the data set seems to support our claim that identifying subgroups of genes by visual inspection may overfit the data: with the exception of the two repeats of the DKFZP566O1646 protein, our model assigns each group of repeated genes to the same cluster, whereas four of the repeated genes are assigned to different subgroups in ref. 8. Details are shown in Table 1. The risks of overfitting by visual inspection can be easily appreciated by looking at the color patterns in Fig. 2. As the dendrogram is built, genes with highly similar temporal profiles are merged first, thus producing subtrees with similar patterns of colors. However, according to our analysis, the data do not provide enough evidence to conclude that such subtrees contain time series generated by different processes and they are therefore merged into a single cluster.

An example of this phenomenon is shown in detail by Fig. 3, which enlarges part of the dendrogram in Fig. 2. The subtree on the top half of the figure merges 29 time series that appear to be more homogenous to visual inspection, and the large Bayes factors, in log scale, which shows that at each step of the iterative procedure, merging the time series determines a model which is more likely than the model determined by not merging them. Similarly, the bottom half subtree merges 16 time series that appear to be more similar. The Bayes factors attached to the terminal node of the picture are exp(33), meaning that the model in which the two subtrees are merged together is exp(33) times more likely than the model in which these subtrees are taken as two separate clusters.