A permutation test for determining significance of clusters with applications to spatial and gene expression data

Measuring Clustering Quality

We use the pooled within-cluster variance matrix to measure the goodness of clustering [21]. This is similar to the idea behind hypothesis testing in Multivariate Analysis of Variance (MANOVA), which involves decomposing the total variance into between-cluster and within-cluster components and comparing their relative sizes [22]. We compute the between, within, and total variance matrices, respectively, for g groups with n_l elements in group l (l = 1, …, g) as follows:

Here x_lj is the element j in group l, x̄_l is the average of the elements in group l, and x̄ is the average over all the elements. We say that the difference between groups exists if the within-cluster variance is sufficiently small in some metric when compared to the between-cluster variance. There are several formal statistics to assess this equality of multivariate means, they include Wilks’ lambda statistic, λ = |W|=|B + W| = 1=|BW⁻¹ + I|; the Lawley-Hotelling Trace, tr(BW⁻¹); the Pillai Trace, tr(B(B+W) ⁻¹); and Roy’s largest root, the maximum eigenvalue of BW⁻¹. If one assumes normality, then as sample sizes increase, these tests become equivalent, and p-values can be obtained from the appropriate asymptotic distributions. Wilks’ lambda is the most common statistic, as it is equivalent to the likelihood ratio test under the assumption of normality. Advantages and disadvantages of various statistics are described, for example, in [23].

The major disadvantage of using these asymptotic results, however, is that they require the assumption of normality. While this may be valid in some cases, there is usually no a priori reason to assume it to be true. In particular, it does not seem to apply to microarray data, which provides the original motivation for this work. Some have noted that the gene expression levels might be close to a log-normal distribution (marginally), but this assumption is difficult to justify in many cases [24]. Without normality, asymptotic distributions do not seem to hold readily, and thus testing hypotheses is not feasible. This is of course even more difficult for small samples, and we cannot standardize the statistic to obtain a p-value.

In general, B and W offer an attractive way of measuring goodness of clusters and they should be used together to compute a single measure. In the permutation test we propose below, however, we note that T = W + B always remains constant across all permutations. Therefore, it is sufficient to consider only how W varies for our statistic. A slightly different measure of cluster quality is the Mann-Whitney U-test [25, 26] and it appears, in our experience, to give comparable results.

Using the singular values

We now discuss some choices regarding the use of W. In the univariate case, the group sum-of-squares is used as a natural criterion in the analysis of variance (ANOVA). The simplest generalization of this sum-of-squares to the multivariate case is the trace, $\text{tr}(\mathbf{W})={\sum }_{i=1}^n{W}_{ii}$, where the within-group sum-of-squares is computed over all the variables. In the context of hierarchical clustering, this idea is related to that of the Ward method [27], in which the pair of clusters that minimizes the sum-of-squares is combined among all possible pairs at each step of the algorithm.

While the trace criterion is often used, there are some problems associated with it. One problem is that it implicitly tries to produce spherical clusters in p dimensions. Thus, when this is used as a clustering measure, this can result in an incorrect set of clusters. For example, it will fail to identify ellipsoidal clusters that are very close to each other [1]. Another problem with the trace is that it is scale-dependent. When the raw data are standardized differently, different clusters will result by this criterion.

Some prefer the determinant, |W|. However, it still makes the implicit assumption that all the clusters have the same shape, but this requirement is less stringent than that of the trace. The determinant criterion is scale-independent and is not affected by normalization of the data. For more discussion of the use of the determinant instead of the trace, see [23]. There are other advantages to using |W|, such that highly correlated components are not given excessive weight [28].

While the determinant is preferred in general, a problem occurs when the number of components is greater than the number of samples (p > n). This is the case when clustering the samples in microarray data, which has a distinguishing characteristic that the number of genes typically far exceeds the number of samples. There are several ways of screening genes to assess their potential relevance and filter them to obtain a smaller set of genes, but in most case p is still much larger than n and, as a result, W is singular.

A solution in this case is to compute the quantity ${\sum }_{i=1}^r{\sigma }_i$, where r is the rank of W (rank(W) = r < min(n, p)) and σ_i are the (non-zero) singular values of W. Recall that by singular value decomposition, W_n×p = U_n×n D_n×p V_p×p, where U and V^T are orthonormal (UU^T = U^TU = I_n×n and VV^T = V^TV = I_p×p), and the entries of the diagonal matrix D are the singular values σ₁ ≥ σ₂ ≥ … ≥ σ_r ≥ 0 = ··· = 0.

So when the matrix is 1) singular because p > n; or 2) p < n after filtering but still nearly singular due to highly co-expressed genes, we can use the product of the singular values in place of the determinant. In general, deciding on the right threshold for counting a singular value as zero is not simple; but our threshold based on the product of the largest singular value and the machine precision appears to work well. For computing the singular values, fast and stable algorithms are available. If the number of genes is very large or the computation is limited by the memory size, we can resort to the trace, which is not affected by the problem of a singular matrix. In our applications below, the trace and the product of singular values give similar results, with the trace generally giving slightly larger p-values.

Implementation of the the permutation approach in this setting is simple. Following Pitman [29], in spirit, at each node of the dendrogram we compare the value of the statistics to the values of the statistic we would have obtained if the group memberships had been randomly assigned. In each random sample, the same data points are used and the number of elements in each of the two groups remains the same, but different sets of elements are assigned to each group. We perform this membership permutation many times, and compute a p-value by counting the proportion of times that the same or a more extreme value of the statistic is obtained in the random samples. For small n, we can do a complete enumeration; in general, we can only sample from the permutation distribution.

An illustrative example

Figure 1(a) shows a two-dimensional sample dataset with a clearly defined cluster structure. Two sets of ten data points, numbered {1{10} and {11{20}, were generated from two multivariate normal distributions with equal identity covariance matrices but with means 4 units apart in each dimension. In Figure 1(b), we performed a hierarchical clustering using Euclidean distance and average linkage. The p-values shown are from the permutation test based on |W|, using 1000 samples at each node. One can see the presence of two distinct clusters, as evidenced by the low p-value of <0.001 at the top division. Following the left branch, the p-value is 0.010, which suggests that the left cluster that consists of {6,1,5,3,4} is distinct from the right cluster containing {9,8,10,2,7}. There is no clearly defined structure past that point if we are to use the typical p = 0.05 value. Note that without these numbers, it is difficult to determine which nodes may be more significant just from the dendrogram. For example, the node joining {12,18} and {16,20,14,11,13} is similar in appearance to the one joining {6,1,5,3,4} and {9,8,10,2,7}. But the first node has a larger p-value associated with it than the second one. These and other similar observations are consistent with the structure of the data shown in Figure 1(a). Note also that for low dimensional data, such as this example, one may resort to examining the original data. But, this option is generally not available, especially for high dimensional data.

In Figure 2(a), we list the p-values computed on the same data using the tr(W) criterion. The values are mostly similar to those obtained previously, except near the bottom of the dendrogram where the number of objects within each group is too small to give meaningful results. One reason for the dissimilarity is the different assumptions underlying each criterion. For example, because the trace tends to be smaller for spherical structures, the p-values may be inflated slightly if there are adjacent elongated structures correctly captured by the hierarchical clustering. In that case, some cluster-membership permutations may result in more spherical clusters with a smaller trace even though the clustering is incorrect.

After identifying significant nodes, one can proceed to combine the non-significant branches, as is done in [26]. Using a cut-off p-value of 0.05, we revise the dendrogram shown in Figure 1(b) and display it in Figure 2(b). There are only two nodes with p-values less than 0.05 and so there are only three clusters left. Looking at the data in Figure 1(a), we see that the new dendrogram is a good summary of the data. The cut-off p-value can be chosen to reflect the desired level of detail in the dendrogram. Since the repeated tests on subsets of the same data are not independent, we can use a “correction” as explained in Westfall and Young[30]. Since we are more intent on displaying the properties of the dendrogram trimming algorithm, we chose to use the simpler Bonferroni correction idea as a guide in selecting an initial cut-off value, e.g., 0.05/(n − 1), since there are n − 1 nodes.

Spatial data

In the analysis of spatial data, one is often interested in determining if there is any clustering of cases. This clustering can represent any number of interesting aspects, such as excess cases of disease [32], point sources of exposure [33] or the existence of health disparities within a community [34]. Consider the spatial data collected by Alt and Vach [31] who map a medieval burial site in Neresheim, Baden-Wurttemberg, Germany. The authors address an anthropologic question regarding burial patterns; specifically to determine if members of the same family tend to be buried in close proximity as a familial unit. To answer this question, the authors examined skeletal remains of 152 grave sites for evidence of a known hereditary genetic feature: missing or reduced wisdom teeth. The locations of the 30 grave sites with dental feature (referred as cases) are presented in Figure 3(a). In this example, if the proposed hypothesis is correct, and familial units are buried together, then we would expect to see several clusters distributed throughout the grave site. In fact, each (extended) family might be considered its own cluster.

We perform hierarchical clustering on the spatial locations of the cases using Euclidean distances and average linkage. Figure 3(b) shows the resulting dendrogram, from which it is nearly impossible to determine which groupings are important. We calculate p-values for each node from the permutation test based on |W|, using 1000 simulations at each node, which are also shown in Figure 3(b). These p-values suggest 4 distinct clusters, where the node significance is less than Bonferoni corrected p-value, p = 0.001724 among the cases. These clusters are composed of the the cases numbered {15,23,24,30}, {4,5,7,25,26,27,28,29}, {1,2,3,6,8,21}, and {9,10,11,12,13,14,16,17,18,19,20,22}.

Once we have calculated the p-values for the nodes, we once again, combine the non-significant branches to simplify the dendrogram. We use a Bonferoni p-value cutoff of 0.05/(n − 1) = 0.00172, where n is the number of nodes, to determine significance, and we display the resulting dendrogram in Figure 4(a).

We then turn our attention to examining the locations of the non-affected graves, the controls. If the hypothesis of familial burials is correct, we would expect to see many small clusters, similar in size, to the clusters of the cases. We perform the same analysis on the controls that we applied to the cases and the resultant dendrogram is presented in Figure 4(b). As expected, we see evidence of many small clusters within the spatial distribution of the non-affected graves.

Application to microarray data

We apply our technique to three microarray datasets. The first two are for clustering patient samples and the third is for clustering genes. The advantage of the technique becomes more apparent as the size of the dataset grows larger, but we concentrate on clustering samples rather than genes for illustration due to space constraints. We first consider RNA data studied in [35], which was used to demonstrate how cigarette smoking status effects the the human airway epithelial cell transcriptomethe based on gene expression levels. Among the 72 patients, there are 34 patients identified as ‘current smokers’ and 23 identified as ‘never smokers’. The remaining 15 subjects are identified as ‘former smokers’, which we do not include in our classification, for the sake of clarity. Those patients numbered {1–34} are ‘current smokers’ and those numbered {35–57} are ‘never smokers’. From the more than 8000 genes, we select 848 genes, highly correlated to smoking status, for classification purposes. There are still more dimensions than objects, however, and we use the trace criterion for the permutation test. The result is plotted in Figure 5(a). After some normalizations to the data, the hierarchical clustering gives a good result, with clearly defined clustering between never smokers and current smokers.

In Figure 5(b), we use the cut-off p-value of .05/(n − 1) = 0.000893 (Bonferroni correction) to eliminate the nodes that are not significant. This pruning of the dendrogram provides a much simpler and more informative view of the data.

In Figure 6(a), we apply the same technique to a dataset on embryonal tumors of the central nervous system [36]. This dataset consists of several types of heterogeneous tumors about which little is known, including medulloblastomas, malignant gliomas, renal/extrarenal rhabdoid tumors, supratenotorial primitive neuroectodermal tumors, as well as normal cerebella. A main conclusion of the paper is molecular profiles from those of other embryonal CNS tumors[36]. For normalization of the data, we have tried to follow the steps described in the supplemental information of the article, with the additional step of using variance filter to reduce the number of genes to 100. Different numbers of genes gave similar results.

There are more disease categories here than in the epithelial data, but the significance-adjusted dendrogram appears to capture the important part of the dendrogram very well. As was the case in [36], PNET samples are dispersed among different clusters; NCER (normal) samples are separated from the rest at the top branch; most other samples are clustered with those of the same type. The modified dendrogram presents a clear and concise description for easier interpretation.

For both the epithelial and the CNS tumor data, using the trace criterion appears to inflate the p-values slightly, causing more branching in the tree. In the second case, for example, it results in an additional subdivision of the MD branch. In both cases, the Bonferroni corrected p-value of 0.05 appears to give a clear result. For Figures 5(a) and 6, we used 10000 permutations at each node.

We can also apply the permutation test with genes as objects and samples (experiments) as variables, although the size of the dendrogram makes it harder to examine the details in that case. In Figure 7, we cluster the genes using the compendium data of S. cerevisiae experiments [17]. We eliminate those genes with too many missing values and carry out the missing value imputation using k nearest neighbors for the rest. Because the dimension of W is equal to the number of experiments, the computational cost for clustering genes is low; however, we selected only the 500 genes with the largest variations due to the space available for visualization. For clustering genes, the W matrix is not singular, unless the experiments are nearly identical. Figure 7(a) and 7(b) show the clustering result before and after the modification by the permutation test, again using 10000 samples and p-value cut-off of 0.05/(n −1). As before, the simpler structure facilitates the examination and interpretation of the co-expressed genes. Further work would involve annotation of the main clusters identified using Gene Ontology and other annotation databases.