Model order selection for bio-molecular data clustering

Data perturbations using randomized maps with bounded distortions

A major requirement for clustering algorithms is the reproducibility of their solutions with other data sets drawn from the same source; this is particularly true with bio-molecular data, where the robustness of the solutions is of paramount importance in bio-medical applications. From this standpoint the reliability of a clustering solution is tied to its stability: we may consider reliable a cluster if it is stable, that is if it is maintained across multiple data sets drawn from the same source. In real cases, however, we may dispose only of limited data, and hence we need to introduce multiple " small" perturbations into the original data to simulate multiple "similar" samples from the same underlying unknown distribution. By applying appropriate indices based on similarity measures between clusterings we can then estimate the stability and hence the reliability of the clustering solutions.

We propose to perturb the original data using random projections μ : ℝ ^d → ℝ ^d' from high d-dimensional spaces to lower d'-dimensional subspaces. A related approach is presented in [17], where the authors proposed to perturb the data randomly choosing a subset of the original features (random subspace projection [18]); the authors did not propose any principled method to choose the dimension of the projected subspace, but a key problem consists in finding a d' such that for every pair of data p, q ∈ ℝ ^d , the distances between the projections μ(p) and μ(q) are approximately preserved with high probability. A natural measure of the approximation is the distortion dist _μ :

If dist _μ (p, q) = 1, the distances are preserved; if 1 - ε ≤ dist _μ (p, q) ≤ 1 + ε, we say that an ε-distortion level is introduced.

In [15] we experimentally showed that random subspace projections used in [17] may introduce large distortions into gene expression data, thus introducing bias into stability indices based on this kind of random projections. For these reasons we propose to apply randomized maps with guaranteed low distortions, according to the Johnson-Lindenstrauss (JL) lemma [19], that we restate in the following way: Given a d-dimensional data set D = {p₁, p₂,...,p _n } ⊂ ℝ ^d and a distortion level ε, randomly choosing a d'-dimensional subspace S ⊂ ℝ ^d , with d' = c logn/ε², where c is a suitable constant, with high probability (say ≥ 0.95) the random projection μ :ℝ ^d → S verifies 1 - ε ≤ dist _μ (p _i , p _j ) ≤ 1 + ε for all p _i ≠ p _j .

In practice, using randomized maps that obey the JL lemma, we may perturb the data introducing only bounded distortions, approximately preserving the metric structure of the original data [15]. Note that the dimension of the projected subspace depends only on the cardinality of the original data and the desired ε-distortion, and not from the dimension d of the original space.

The embedding exhibited in [19] consists in projections from ℝ ^d in random d'-dimensional subspaces. Similar results may be obtained by using simpler maps [20, 21], represented through random d' × d matrices $R=1/\sqrt{d^{\prime }}(r_{ij})$, where r _ij are random variables such that:

E[r _ij ] = 0, Var[r _ij ] = 1

Using the above randomized maps (with the exception of RS projections), the JL lemma guarantees that, with high probability, the "compressed" examples of the data set represented by the matrix D _R = RD have approximately the same distance (up to a ε-distortion level) of the corresponding examples in the original space, represented by the columns of the matrix D, as long as d' ≥ c logn/ε².

We propose a general MOSRAM (Model Order Selection by RAndomized Maps) algorithmic scheme, that implements the above ideas about random projection with bounded distortions to generate a set of similarity indices of clusterings obtained by pairs of randomly projected data. The main difference with respect to the method proposed in [13] is that by our approach we perturb the original data using a randomized mapping μ : ℝ ^d → ℝ ^d' :

Stability indices based on the distribution of the similarity measures

Using the similarity measures obtained through the MOSRAM algorithm, we may compute stability indices to assess the reliability of clustering solutions.

More precisely, let $\mathcal{C}$ be a clustering algorithm, ρ a random perturbation procedure (e.g. a resampling or a random projection) and sim a suitable similarity measure between two clusterings (e.g. the Fowlkes and Mallows similarity).

We may define the random variable S _k , 0 ≤ S _k ≤ 1:

S _k = sim($\mathcal{C}$(D₁, k), $\mathcal{C}$(D₂, k))     (2)

where D₁ = ρ⁽¹⁾(D) and D₂ = ρ⁽²⁾(D) are obtained through random and independent perturbations of the data set D; the intuitive idea is that if S _k is concentrated close to 1, the corresponding clustering is stable with respect to a given controlled perturbation and hence it is reliable.

Let f _k (s) be the density function of S _k and F _k (s) its cumulative distribution function. A parameter of concentration implicitly used in [13] is the integral g(k) of the cumulative distribution:

Note that if S _k is centered in 1, g(k) is close to 0, and hence it can be used as a measure of stability. Moreover, the following facts show that g(k) is strictly related to both the expectation E[S _k ] and the variance Var[S _k ] of the random variable S _k :

A χ²-based test for the assessment of the significance of the solutions

In this section we propose a method for automatically finding the " optimal" number of clusters and to detect significant and possibly multi-level structures simultaneously present in the data. First of all, let us consider the vector (ξ₂, ξ₃,...,ξ_H+1) (eq. 6) computed by using the output of the MOSRAM algorithm. We may perform a sorting of this vector:

where p is the permutation index such that ξ_p(1)≥ ξ_p(2)≥ ... ≥ ξ_(H). Roughly speaking, this ordering represents the "most reliable" p(1)-clustering down to the least reliable p(H)-clustering; exploiting this we would establish which are the significant clusterings (if any) discovered in the data.

To this end, for each k ∈ $\mathcal{K}$ = {2, 3,...,H + 1}, let us consider the random variable S _k defined in eq. 2, whose expectation is our proposed stability index. For all k and for a fixed threshold t ^o ∈ [0, 1] consider the Bernoulli random variable B _k = I(S _k > t ^o ), where I is the indicator function: I(P) = 1 if P is True, I(P) = 0 if P is False. The sum $X_k={\displaystyle {\sum }_{j=1}^m{B}_k^j}$ of i.i.d. copies of B _k is distributed according to a binomial distribution with parameters m and θ _k = Prob(I(S _k > t ^o )).

If we hypothesize that all the binomial populations are independently drawn from the same distribution (i.e. θ _k = θ, for all k ∈ $\mathcal{K}$), for sufficiently large values of m the random variables $\frac{X_k-m{\theta }_k}{\sqrt{m{\theta }_k(1-{\theta }_k)}}$ are independent and approximately normally distributed. Consider now the random variable:

This variable is known to be distributed as a χ² with |$\mathcal{K}$| - 1 degrees of freedom, informally because the constraint $\widehat{\theta }$ between the random variables X _k , k ∈ $\mathcal{K}$ introduces a dependence between them, thus leading to a loss of one degree of freedom. By estimating the variance mθ(1 - θ), with the statistic m $\widehat{\theta }$(1 - $\widehat{\theta }$), we conclude that the following statistic

is approximately distributed according to ${\chi }_{\left|\mathcal{K}\right|-1}^2$ (see, e.g. [24] chapter 12, or [25] chapter 30 for more details).

A realization x _k of the random variable X _k (and the corresponding realization y of Y) can be computed by using the output of the MOSRAM algorithm:

Using y, we can test the following alternative hypotheses:

If $y\ge {\chi }_{\alpha ,\left|\mathcal{K}\right|-1}^2$ we may reject the null hypothesis at α significance level, that is we may conclude that with probability 1 - α the considered proportions are different, and hence that at least one k-clustering significantly differs from the others.

The output of the proposed procedure is the set of the remaining (top sorted) k-clusterings that correspond to the set of the estimate stable number of clusters (at α significance level). Equivalently, following the sorting of ξ, we may compute the p-value (probability of committing an error if we reject the null hypothesis) for all the ordered groups of clusterings from the p(1)...p(H) to the the p(1), p(2) group, each time removing the bottom ranked clustering from the ξ vector. Note that if the set of the remaining top-ranked clusterings contains more than one clustering, we may find multiple structures simultaneously present in the data (at α significance level).

Detection of multiple levels of structure in synthetic data

To show the ability of our method to discover multiple structures simultaneously present in the data, we propose an experiment with a 1000-dimensional synthetic multivariate gaussian data set (sample1) with relatively low cardinality (60 examples), characterized by a two-level hierarchical structure, highlighted by the projection of the data into the two main principal components (Figure 1): indeed a two-level structure, with respectively 2 and 6 clusters is self-evident in the data.

Two clusterings (using the Prediction Around Medoid algorithm) are detected at 10^-4 significance level by applying our MOSRAM algorithm and the proposed χ²-based statistical procedure. In particular we performed 100 pairs of Bernoulli projections with a distortion bounded to 1.2 (ε = 0.2), yielding to random projections from 1000 to 479-dimensional subspaces. Indeed Table 1 reports the sorted means of the stability measures together with their variance and the corresponding p-values computed according to the proposed χ²-based statistical test, showing that 2 and 6-clusterings are the best scored, as well as the most significant k-clusterings discovered in the data. This situation is depicted in Figure 2, where the histograms of the similarity measures for k = 2 and k = 6 clusters are tightly concentrated near 1, showing that these clusterings are very stable, while for other values of k the similarity measures are spread across multiple values. Note that the p-values of Table 1 (as well as the p-values of Table 2 and 3) refer to the probability of committing an error if we reject the null hypothesis. The clusterings with p ≥ α are considered equally reliable (in this case the null hypothesis cannot be rejected), while the clusterings for which p <α are considered less reliable (at α significance level).

k	mean	variance	p-value
2	1.0000	0.0000	1.0000
6	1.0000	0.0000	1.0000
7	0.9217	0.0016	0.0000
8	0.8711	0.0033	0.0000
9	0.8132	0.0042	0.0000
5	0.8090	0.0104	0.0000
3	0.8072	0.0157	0.0000
10	0.7715	0.0056	0.0000
4	0.7642	0.0158	0.0000

k	mean	variance	p-value
2	0.8285	0.0077	1.0000
3	0.8060	0.0124	0.7328
4	0.6589	0.0060	2.3279e-06
5	0.6012	0.0073	9.5199e-11
6	0.5424	0.0057	6.3282e-15
7	0.5160	0.0062	0.0000
8	0.4865	0.0050	0.0000
9	0.4819	0.0060	0.0000
10	0.4744	0.0049	0.0000

k	mean	variance	p-value
2	0.9566	0.0028	1.0000
3	0.7900	0.0149	0.0000
4	0.6963	0.0128	0.0000
5	0.6387	0.0075	0.0000
6	0.6135	0.0082	0.0000
7	0.6129	0.0079	0.0000
9	0.5864	0.0063	0.0000
8	0.5792	0.0079	0.0000
10	0.5744	0.0058	0.0000

Experiments with DNA microarray data

To show the effectiveness of our methods with gene expression data we applied MOSRAM and the proposed statistical test to Leukemia [29] and Lymphoma [1] samples. These data sets have been analyzed with other model order selection algorithms previously proposed [10, 13, 30–32]: at the end of this section we compare the results obtained with the cited methods with our proposed MOSRAM algorithm.

Leukemia

This well known data set [29] is composed by 72 leukemia samples analyzed with oligonucleotide Affymetrix microarrays. The Leukemia data set is composed by a group of 25 acute myeloid leukemia (AML) samples and another group of 47 acute lymphoblastic leukemia (ALL) samples, that can be subdivided into 38 B-Cell and 9 T-Cell subgroups, resulting in a two-level hierarchical structure. We applied the same pre-processing steps performed by the authors of the Leukemia study [29], obtaining 3571 genes from the original 7129 gene expression values. We further selected the 100 genes with the highest variance across samples, since low variance genes are unlikely to be informative for the purpose of clustering [10, 31]. We analyzed both the 3571-dimensional data and the data restricted to the 100 genes with highest variance, using respectively Bernoulli projections with ε ∈ {0.1, 0.2, 0.3, 0.4} and projections to 80-dimensional subspaces. In both cases the k-means clustering algorithm has been applied.

Figure 3 summarizes the results using gene expression levels of the genes with highest variance. Table 2 reports the sorted means of the stability measures together with their variance and the corresponding p-values computed according to the proposed χ²-based statistical test. By these results 2 and 3 clusters are correctly predicted at α = 10^-5 significance level. Indeed the empirical means of the stability measures (eq. 6) for 2 and 3 clusters are quite similar and the corresponding lines (black and red) of the empirical cumulative distributions (ecdfs) cross several times, while the other ecdfs are clearly apart from them (Figure 3). The same results are approximately obtained also using the 3571-dimensional data with random projections to 428-dimensional subspaces (ε = 0.2), but 2 and 3 clusters are predicted at α = 10^-13 significance level. Similar results are also achieved with ε = 0.1 and ε = 0.3, while with ε = 0.4 the results are less reliable due to the relatively large distortion induced (data not shown). Also using the PAM [27] and hierarchical clustering algorithms with the Ward method [33] we obtained a two-level structure with 2 and 3 clusters at α = 10^-5 significance level.

Lymphoma

Three different lymphoid malignancies are represented in the Lymphoma gene expression data set [1]: Diffuse Large B-Cell Lymphoma (DLBCL), Follicular Lymphoma (FL) and Chronic Lymphocytic Leukemia (CLL). The gene expression measurements are obtained with a cDNA microarray specialized for genes related to lymphoid diseases, the Lymphochip, which provides expression levels for 4026 genes [34]. The 62 available samples are subdivided in 42 DLBCL, 11 CLL and 9 FL. We performed pre-processing of the data according to [1], replacing missing values with 0 and then normalizing the data to zero mean and unit variance across genes. As a final step, according to [10], we further selected the 200 genes with highest variance across samples, obtaining a resulting data set with 62 samples and 200 genes. As in the previous experiment, we processed both the high-dimensional original data and the data with the reduced set of high-variance genes, using respectively Bernoulli projections with ε ∈ {0.1, 0.2, 0.3, 0.4} and projections to 160-dimensional subspaces. The k-means clustering algorithm has been applied.

The results with highest variance genes are summarized in Figure 4 and Table 3. The statistical test identifies as significant only the 2-clustering. Indeed, looking at the ecdf of the stability index values (left), the 2-clustering (black) is clearly separated from the others. The 3 (red) and 4-clustering (green) graphs, are quite distinct from the others, as shown also by the corresponding empirical mean of the stability index values (Figure 4), but they are also clearly separated from the 2-clustering curve. Accordingly, our proposed χ²-based test found a significant difference between the 2-clustering and all the others. Similar results are obtained also with hierarchical clustering and PAM algorithms. Using all the 4026 genes and Bernoulli random projections to 413-dimensional (ε = 0.2) subspaces with the Ward's hierarchical clustering algorithm our method finds as significant the 3-clustering as well as the 2-clustering. In this case also similar results are obtained with ε = 0.1 and ε = 0.3. It is worth noting that the subdivision of Lymphoma samples in 3 classes (DLBCL, CLL and FL) have been defined on histopathological and morphological basis and it has been shown that this classification does not correspond to the bio-molecular characteristics and to clinical outcome classes of non Hodgkin lymphomas. In particular studies based on the gene expression signatures of the DLBCL patients [1] and on their supervised analysis [35], showed the existence of two subclasses of DLBCLs. Moreover Shipp et al. [36] highlighted that FL patients frequently evolve over time and acquire the clinical features of DLBCLs, and Lange et al. [10] found that a 3-clustering solution groups together FL, CLL and a subgroup of DLBCLs, while another subgroup of DLBCLs sets up another cluster, even if the overall stability of the clustering is lower with respect to the 2-clustering solution. The relationships between FL and subgroups of DLBCL patients are confirmed also by recent studies on the individual stability of the clusters in DLBCL and FL patients [15]. These considerations show also that the stability analysis of patients clusters in DNA microarray analysis are only the first step to discover significant subclasses of pathologies at bio-molecular level, while another necessary step is represented by the bio-medical validation.