Generalized Clustergrams for Overlapping Biclusters

1 Introduction

Biological processes are extremely complex, showing a hierarchical organization at various levels of granularity. On the other hand, many clustering methods require the number of clusters to be given as input. But in real gene expression data [Eisen et al., 1998] one cannot unequivocally determine a well-defined number of clusters, as coarser-grained clusters may exhibit progressively finer-grained structure. For example, Figure 1 shows the evolution of the error of nonnegative decompositions (biclusterings) as the number of clusters is increased in the case of two datasets: a synthetic dataset with 5 biclusters (Figure 1a) and a large colon cancer dataset (Figure 1b).

Note that in the synthetic dataset, a well-defined number of clusters (nc = 5) can be discerned, while the colon cancer dataset shows a steep drop in error until nc = 3, followed by a series of progressively smaller ones. The latter could represent either overfitting or finer-grained substructure and only a more in-depth biological analysis can settle the issue. Anyhow, it is not enough from a biological point of view to determine a unique number of clusters and to analyze the resulting clusters at that fixed dimensionality. Instead, we need to perform some sort of multi-scale analysis of the clusters and their evolution as the dimensionality (nc) of the clustering is varied.

For this, it is essential to be able to determine the relationships between the clusters generated at different nc 1 , which is possible only if we can guarantee the stability of the clusters and if we have a mathematically sound method of comparing biclusters. (We deal with biclusters [Cheng and Church, 2000] since in the case of gene expression data, as in many other domains, objects/genes tend to be correlated only for certain subsets of samples, i.e. specific biological contexts.)

We achieve biclustering stability using the metaclustering approach of [Badea, 2005; Badea and Tilivea, 2007], which is based on a positive tensor factorization (PTF) of the biclusters obtained in various repeated clustering runs. Bicluster comparisons are "built into" this approach in an elegant manner.

These two key features allow us to construct a series of decompositions at varying nc and to trace the evolution of the individual clusters as nc changes. This generalizes hierarchical clustering dendrograms to a more complicated setting where items that have been grouped (at a certain nc) can be separated later on at a lower dimensionality.

Hierarchical clustering is probably the most frequently used clustering method in the domain of gene expression data analysis [Eisen et al., 1998]. This is due not only to its simplicity, but also to the very intuitive nature of the clustering dendrograms it produces, which graphically depict the evolution of clusters at various dimensionalities. Since choosing a unique well-defined number of clusters is usually problematic in the case of gene expression data, being able to trace the dimensional evolution of clusters is essential in this domain.

1 For example, we may want to analyze in more detail the cluster merges and splits as nc is varied.

However, extensive experimental evaluations [Thalamuthu et al., 2006] have demonstrated that hierarchical clustering shows mediocre performance especially for noisy data, with many "unrelated" (or, so called "scattered") genes. Thus, it might seem that performance has been sacrificed for better visualization and that an inherent trade-off exists between the two. In the following we show that we can trace the dimensional evolution of clusters in a more complicated setting involving biclustering based on PTF by developing a generalization of dendrograms for the case of overlapping biclusters.

a. synthetic data (5 biclusters)

b. colon adenocarcinoma

2 Stable meta-clustering with PTF

As already mentioned in the Introduction, the ability of clustering genes and samples simultaneously ('biclustering') is essential for analyzing gene expression data, where genes tend to be co-expressed only for certain subsets of samples, corresponding to specific biological contexts. In the following, we use the biclustering approach based on positive tensor factorizations (PTF) of [Badea, 2005; Badea and Tilivea, 2007], which we briefly review here. The increased stability of the meta-clustering approach is essential for being able to follow the dimensional evolution of biclusters.

2.1 Biclustering using Nonnegative Matrix Factorizations with offset

An elegant method of biclustering consists in factorizing the gene expression matrix X as a product of an ns nc (samples clusters) matrix A and an nc ng (clusters genes) matrix S 2

subject to additional nonnegativity constraints:

which express the obvious fact that expression levels and cluster membership degrees cannot be negative.

Factorization (1) differs from the standard Nonnegatuive Matrix Factorization (NMF) [Lee and Seung, 1999; 2000] by the additional "gene offset" So, whose main role consists in absorbing the constant expression levels of genes, thereby making the cluster samples Scg "cleaner".

The factorization (1-2) can be regarded more formally as a constrained optimization problem:

subject to the nonnegativity constraints (2). This problem can be solved using an iterative algorithm with the following multiplicative update rules (which can be easily derived using the method of Lee and Seung [Lee and Seung, 2000]):

where e is a column vector of 1 of size equal to the number of samples and is a regularization parameter (a very small positive number).

2 Xsg represents the gene expression level of gene g in sample s, Scg the membership degree of gene g in cluster c and Asc the mean expression level of cluster (biological process) c in sample s.

The algorithm initializes A, S and So with random entries, so that (slightly) different solutions may be obtained in different runs. (This is due to the non-convex nature of the optimization problem (3), which in general has many different local minima.)

We can view the different solutions obtained by the generalized NMFO algorithm as overfitted solutions, whose consensus we'll need to construct.

We have observed experimentally that adding offsets to standard NMF leads to significant improvements in the quality of the recovered clusters.

More precisely, the genes with little variation are reconstructed by the standard NMF algorithm from combinations of clusters, while NMFO uses the additional degrees of freedom So to produce null cluster membership degrees Scg for these genes. Moreover, NMFO recovers with much more accuracy than standard NMF the original sample clusters, the standard NMF algorithm being confused by the cluster overlaps. This improvement in recovery of the original clusters is very important in our application, where we aim at a correct sub-classification of samples.

2.2 Meta-clustering with PTF

Unfortunately, virtually all clustering methods that are currently used for gene expression data analysis tend to produce highly unstable clusters, especially when clustering genes. (The instability manifests itself either w.r.t. the initialization of the algorithm, as in the case of k-means and NMF, or w.r.t. small perturbations of the dataset in the case of deterministic algorithms, such as hierarchical clustering.)

A frequently used method to obtain more stable clusters consists in building a consensus of several individual clusterings constructed from different NMFO initializations.

More precisely, starting with a number of NMFO runs

a consensus biclustering is constructed using a Positive Tensor Factorization (PTF) [Welling and Weber, 2001] of the biclusters 3 , which simultaneously determines the bicluster correpondence and the consensus biclustering ( , ) [Badea, 2005; Badea and Tilivea, 2007]:

where s are samples, g, genes, c clusters and k metaclusters (or "consensus clusters"). 4 and represent the consensus of A (i) and S (i) respectively. More precisely, the columns k of and the corresponding rows k of make up a base set of bicluster prototypes k k out of which all biclusters of all individual runs can be recomposed, while encodes the (bi)cluster-metacluster correspondence. The factorization (6) can be computed using the multiplicative update rules from [Badea, 2005; Badea and Tilivea, 2007]:

3 A tensor factorization is needed instead of a matrix factorization since biclusters are matrices.

4 To simplify the notation, the indices i and c were merged into a single index (ic).

where ' ' and ' ' represent element-wise multiplication and division of matrices, while ' ' is ordinary matrix multiplication. After convergence of the PTF update rules, the rows of are normalized to unit norm to make the gene clusters directly comparable to each other, whereas the columns of are normalized such that (r is the c i k ic r , ) (

number of runs). Then, NMFO initialized with , ( ) , 0 is run 5 to produce the final factorization X A S + e So.

The nonnegativity constraints of PTF meta-clustering are essential both for allowing the interpretation of k k. as consensus biclusters, as well as for obtaining sparse factorizations. In practice, the rows of the correspondence matrix tend to contain typically one or only very few significant entries.

3 PTF for tracing the dimensional evolution of biclusters

Given a stable set of biclusters (A (i) , S (i) , So (i) ) generated for progressively larger numbers of clusters } ,..., 2 { max c c n i n :

we aim at determining the relationships between the biclusters at nc = i and those at larger nc (e.g. nc = i + 1). (For example, bicluster c1 at nc = i may be very similar to bicluster c2 at nc = i + 1. Or, bicluster c3 at nc = i may result from merging biclusters c4 and c5 at nc = i + 1.)

To achieve this, we start with a "reference" factorization that is fine-grained enough to approximate any lower-dimensional factorization (for example with ) , , ( ) ( ) ( ) ( ref ref ref So S A c c n n max 6 ).

We then express the biclusters (A (i) , S (i) , So (i) ) obtained at various dimensionalities in terms of the fine-grained biclusters of the reference factorization as follows:

g i i g So 0 ) ( 00 ) ( i i r ) ( 00 5 0 is obtained from the 1-dimensional NMF decomposition with the normalization .

max c n 6 We can estimate using the error curve as in Figure 1.

where , , for c � 0 and for k � 0. 1 ) ( 0 i s A 1 ) ( 0 ref s A 0 ) ( 0 i c 0 ) ( 0 i k

More precisely, (9) can be expressed as an optimization problem

subject to the constraints . 0 ) (i ck

Note that (10) is similar to the positive tensor factorization (6), except that in (10) only are free variables. Moreover, since in (10) the factorizations (i) and the reference factorization have different dimensionalities, the consensus tensor may contain several significant entries , ,… on a given row c (corresponding to a cluster c formed by superposing (merging) reference clusters k1, k2, …). ) (i ck ) (i ck ) ( 1 i ck ) ( 2 i ck

It can be shown that the decomposition (10) is computable using the following multiplicative update rule:

where we have merged indices i and c as (ic):

The cluster correspondence matrix generalizes hierarchical clustering dendrograms and will be called in the following "generalized clustergram". For example, Figure 2.c shows for the clustering dendrogram from Figure 2.a. can thus be used to trace the evolution of individual clusters as nc is varied. In the case of hierarchical clustering, once two items have been grouped (at a given nc), they remain grouped for all smaller nc. The clustergram from Figure 2.b and the generalized clustergram (correspondence matrix) from Figure 2.c can be used to describe more complex cluster evolutions such as that shown in Figure 3, which presents the evolution of PTF decompositions for nc ranging from 2 to 8 using a 15-dimensional reference factorization of a synthetic dataset with 5 partially overlapping biclusters (for reasons of layout, Fig. 3 presents the transpose of ).

For visualization purposes, the rows c of k ic ) ( were permuted such that the cluster c from decomposition i matches cluster c of decomposition i - 1 for all i and c = 1, …, i - 1. This was achieved in a greedy manner as follows:

for i = 2 … max for c = 1, …, i - 1 c n let max k ) 1 ( max arg i ck k and s.t. c' is not assigned ) ( ' ' max max max arg i k c c c assign row c of decomposition i to position c max assign remaining row of decomposition i to position i The tests of the algorithm on syn thetic datasets showed that it is quite effective at tracing the incremental evolution of the individual biclusters as the number of clusters is increased, despite the variability of the individual NMF clustering runs. (As far as we know, this is the first algorithm to achieve this.) For example, Figure 3 shows that the clusters for nc=5 are essentially preserved when going to nc=6, except for cluster 3:5 (i.e. cluster 3 for nc=5), which splits into two overlapping clusters 2:6 and 5:6.

4 Multi-scale analysis of a colon cancer dataset

We have applied our approach to a large gene expression dataset of colon adenocarcinoma. Our preliminary investigations of this dataset have shown that the disease has several distinct subclasses, but no unique well-defined dimensionality could be determined (as discussed in the Introduction - see also Figure 1.b).

The most frequent colon cancer type, sporadic colon adenocarcinoma, is very heterogeneous and its best current classification based on the presence or absence of microsatellite instabilities (MSI-L, MSI-H and MSS) [Jass et al., 1999] is far from ideal from the point of view of gene expression. A more refined analysis of the biclusters generated at different dimensionalities is therefore needed in order to be able to determine the biologically relevant subclasses of this disease.

To obtain a more accurate subclassification based on gene expression profiles, we have applied our approach to a large dataset (204 samples) containing 182 colon adenocarcinoma samples from the expO database [expO] and 22 control ("normal") samples from [Hong et al., 2007]. (All of these had been measured on Affymetrix U133 Plus 2.0 chips.)

The combined raw scanning data was preprocessed with the RMA normalization and summarization algorithm [Irizarry et al., 2003]. (The logarithmic form of the gene expression matrix was subsequently used, since gene expression values are approximately log-normally distributed.) After eliminating the probe-sets (genes) with relatively low expression as well as those with a nearly constant expression value 7 , we were left with 3666 probe-sets. Finally, the Euclidean norms of the expression levels for the individual genes were normalized to 1 to disallow genes with higher absolute expression values to overshadow the other genes in the factorization.

To estimate the number of clusters (nc), we compared the dimensional evolution of the decomposition error F F rel X So e S A X / of PTF meta-clustering of the real dataset with that corresponding to a randomized dataset 8 (similar to [Kim and Tidor, 2003], see Figure 1.b).

7 Only genes with an average expression value over 100 and with a standard deviation above 150 were retained.

Figure 4 depicts the evolution of PTF decompositions of the colon adenocarcinoma dataset for nc ranging from 2 to 8 using a 15-dimensional reference factorization. Note that as we go from nc = 2 to nc = 3, a new cluster appears, involving fine-grained subclusters (rows) k = 7,14,6,2. This cluster is largely conserved at nc = 4, where another new cluster (involving k = 4,8) is introduced. At nc = 5, the latter suffers a complex transition: it essentially splits in two, but each of the two subclusters picks up additional contributions, so that they end up as clusters 5:5 (corresponding to k = 8,7,9) and respectively 2:5 (k = 4,15). The fact that additional subclusters are involved at nc = 5 as compared to nc = 4 suggests that there are at least 5 biologically significant clusters.

At nc = 6, the main change w.r.t. nc = 5 consists in the appearance of cluster 4:6 (k = 8,2,10). 9 It may be interesting to note that subcluster k = 10 is formed virtually only by probesets associated to the XIST ('X inactive-specific transcript') gene, which is the major effector of X chromosome inactivation (normally expressed only in females by inactivated X chromosomes). Actually, cluster 4:6 is generated, together with 5:6, from 5:5. At nc = 7, the two clusters 4:6 and 5:6 are merged and split again to 2:7 (which is almost identical to 5:5) and 7:7 (which contains subclusters k = 2,10). Therefore, the transition from nc = 5 to nc = 7 amounts to adding two new clusters 4:7 and 7:7.

For a biological interpretation of the biclusters, we have used several high-throughput studies of microsatelliteinstability in colon cancer (e.g. [Banerjea et al., 2004]). Clusters 1:5, 1:6 and 3:7 correspond to a very well defined microsatellite stable subtype (MSS as defined in [Jass et al., 1999]). More specifically, the keratin 23 gene, which we find specifically over-expressed in these clusters was previously known to be specific to MSS colon cancer [Birkenkamp et al., 2007]. Some of the most important genes involved in this bicluster are the following transcription factors: ASCL2, DACH1, FOXQ1, EREG, TNRC9, all previously related to colon cancer and which presumably control the various biological processed involved in microsatellite stable colon carcinoma.

Clusters 5:5, 5:6, 2:7 correspond to the microsatellite instability-high (MSI-H) subtype, which very interestingly over-expresses the developmental homeobox genes HOXC6, PRRX1, HOP. Moreover, the genes overexpressed in this cluster are typically overexpressed in MSI-H tumors with BRAF (rather than KRAS) mutations (using [Kim et al., 2006]).

8 The randomized dataset was obtained by randomly permuting for each gene its expression levels in the various samples. The original distribution of the gene expression levels is thereby preserved.

9 Note that in our method, unlike in the case hierarchical clustering, subclusters that have been merged at a given dimension (for instance k=7,8,9 at nc=5) can later on be separated if the this improves the overall clustering (e.g. at nc=6, k=8 is separated from the other two subclusters and merged with k=2,10).

The normal colon samples are grouped in clusters 4:5, 3:6, 5:7, whose genes are down-regulated in colon cancer. For instance, the down-regulation of carcinoembryonic antigen CEACAM7 is an early event in colorectal tumorigenesis [Thompson et al., 1997], while the critical cell cycle gene CDKN2B is frequently inactivated in colon cancer [Ishiguro et al., 2006].

Our approach to tracing the dimensional evolution of biclusters has proved a useful tool for analyzing the complex subclassification of colon adenocarcinoma that seems to emerge from such gene expression studies.

5 Conclusions

Many real-life data mining datasets and most gene expression datasets exhibit complex substructures at various levels and thus do not have unique well-defined numbers of clusters. Therefore, it is essential to be able to trace the evolution of gene expression biclusters as the number of dimensions of the clustering is varied.

Hierarchical clustering dendrograms are visually very informative in this regard, but have certain essential drawbacks related to their rather mediocre performance for noisy data, as well as to their unidimensional nature (for gene expression data, biclustering methods are more appropriate, since genes tend to be co-expressed only for certain subsets of samples, in certain specific biological contexts). Also, hierarchical clustering produces tree-like cluster merging structures, which may not reflect the dimensional evolution of real-life (biological) processes.

In this paper, we present an original approach that enables tracing the evolution of biclusters as the clustering dimensionality is varied. The method heavily relies on the stability of a positive tensor factorization-based metaclustering of NMF decompositions. We also show that the cluster correspondence matrix of the decomposition of the factorizations produced at various nc (w.r.t. a common reference factorization) can play the role of a generalization of hierarchical clustering dendrograms.

Finally, we applied our approach to a large colon adenocarcinoma dataset, which shows a complex, non-tree-like dimensional evolution of subclusters. Our method has already proven to be useful in the difficult task of determining the biologically relevant subclasses of colon adenocarcinoma.