A Douzal Chouakria

A Diallo

F Giroud

Adaptive clustering of time series

Keywords: Microarray technology, Gene expression data, proximity measure, classification, clustering

HAL is

Introduction

DNA microarray technology allows to monitor simultaneously the expression levels of thousands of genes during important biological processes and across collections of related experiments. Clustering and classification techniques have proved to be helpful to understand gene function, gene regulation, and cellular processes. Though most cells in our bodies contain the same genes, not all of them intervene in each cell: genes are turned on, or expressed when needed. Such specific genes define the molecular pattern related to a specific function of a cell and in most cases appear as organized in molecular regulation networks. To know how cells achieve such specialization, scientists need to identify which genes each type of cell expresses. Microarray technology now allows us to look at many genes at once and determine which are expressed in a specific cell type, i.e. which transcriptome (set of all mRNA or "transcripts") is characteristic of its particular function [START_REF] Eisen | DNA arrays for analysis of gene expression[END_REF]). Researchers are using this powerful technology to learn which genes are turned on or off in diseased versus healthy human tissues for example. The genes that are expressed differently in the two tissues may be involved in causing the disease. In this paper we will be interested in the dynamic progression of cell division cycle through the four distinct phases: G 1 , S, G 2 and M phases. The expression levels of a set of studied genes are then observed at a specific instants of time during cell division cycle. The identification of the set of genes highly characterizing each cell cycle phase is generally based on a prior set of reference genes. Each studied gene is assigned a cell cycle phase by it's peak similarity to the reference genes. This classical approach suffer of two limitations. On the one hand, the most widely used proximity measures between gene expression profiles are based on the closeness of the values regardless to the similarity with respect to (w.r.t.) the genes expression behavior. On the other hand, many different ill-founded sets of reference genes are proposed in the literature, and biologists are not agree about those of genes best characterizing the observed cell cycle phases. Our aim is twice. We propose a new dissimilarity index for gene expression profiles to include both proximity measures w.r.t. expressed values and w.r.t. genes expression behavior. An adaptive unsupervised classification, based on the proposed dissimilarity index, is then performed to identify the cell cycle phases of the studied genes. Finally and assessed by a biological knowledge, we propose a new well-justified set of reference genes. The paper is organized as follows: the next section gives the definition and the properties of the new dissimilarity index. Section 3 presents the human HeLa cell line application, and gives the principal of the proposed adaptive unsupervised classification for cell cycle genes identification. Section 4 performs a comparative analysis and discuss the main obtained results.

Proximity measure between genes expression profile

For clustering or classifying a set of gene expression profiles evolving over time, the commonly used proximity measures are the euclidean distance. Let g 1 = (u 1 , ..., u p ) and g 2 = (v 1 , ..., v p ) be the expression levels of two genes g 1 , g 2 observed at the instant of times (t 1 , ..., t p ). The Euclidean distance δ E between g 1 and g 2 is defined as:

δ E (g 1 , g 2 ) = p i=1 (u i -v i ) 2 1 2
. It stems directly from the above definition that the closeness between two gene expression profiles depends on the closeness of the values regardless to the gene expression behavior. Our aim is to propose a dissimilarity index including both gene expression behavior and values proximity measures. A necessary prior step to the design of such a dissimilarity is to define what we mean about similar gene expression behaviors, and specify the main characteristics that the dissimilarity would measure. We distinguish at least two important characteristics of temporal data. On the one hand, the temporal data where only occurring events, and not their instants of time, are determinant for the proximity evaluation. For instance, in voice processing domain only the occurring syllables are used to identify words; the flow rate being specific to each person. On the other hand, the temporal data where both occurring events and their instants of time are determinant, for instance, ECG, delay response to a treatment, etc. The gene expression data lie within the scope of the latter case.

Behavior proximity measures

We define the similarity w.r.t gene expression behavior by considering two features. On the one hand, the strength of the monotonicity and, on the other hand, the closeness of the growth rates. Without loss of generality, assume that g 1 and g 2 values lie in [0, D]. g 1 and g 2 are similar w.r.t. behavior if at any observed period [t i , t i+1 ] they increase or decrease simultaneously (monotonicity), with a growth rate (closeness of growth rates). One can quantify this similarity concept by considering the classical Pearson correlation coefficient, however this correlation leads to an overestimation when taking into account the temporal dependency between measurements. For more discussion about alternative approaches see [START_REF] Douzal Chouakria | Adaptive dissimilarity index for measuring time series proximity[END_REF].

To overcome this problem we propose the following temporal correlation coefficient:

cort(g 1 , g 2 ) = p-1 i=1 (u (i+1) -u i )(v (i+1) -v i ) p-1 i=1 (u (i+1) -u i ) 2 p-1 i=1 (v (i+1) -v i ) 2
where cort(g 1 , g 2 ) belongs to the interval [-1, 1]. The value cort(g 1 , g 2 ) = 1 signifies that in any observed period [t i , t i+1 ], the genes g 1 and g 2 increase or decrease simultaneously with the same growth rate (similar behavior). The value cort(g 1 , g 2 ) = -1 means that in any observed period [t i , t i+1 ] where g 1 increases, g 2 decreases and vice-versa with a same growth rate (in value; opposite behavior). Finally, cort(g 1 , g 2 ) = 0 expresses that there is no monotonicity between g 1 and g 2 and their growth rates are stochastically linearly independent (different behaviors). For more details about temporal correlation see Chouakria [START_REF] Douzal | Compression Technique Preserving Correlations of a Multivariate Temporal Sequence[END_REF]. Now we present the new dissimilarity index based on the temporal correlation coefficient as a behavior proximity measure.

Dissimilarity index for gene expression profiles

The aim is to provide a new dissimilarity index which would cover both the euclidean distance for the proximity w.r.t. values and the temporal correlation for the proximity w.r.t. behavior. Let us first describe the main specifications of the new dissimilarity index. The dissimilarity index should modulate the proximity w.r.t. values according to the proximity w.r.t. behavior. For the same proximity w.r.t. values, the dissimilarity measure should be dependent on the proximity w.r.t. behavior, and for the same proximity w.r.t. behavior, the dissimilarity should depend on the proximity w.r.t. values. The resulting dissimilarity measure should also allow to adjust the weight contribution between both quantities. The modulating function will increase conventional measure when the temporal correlation decreases from 0 to -1. The resultant dissimilarity should approach the conventional measure if the temporal correlation is zero. The modulating function should decrease the conventional measure when the temporal correlation increases from 0 to +1. According to the specifications above, we propose a dissimilarity index D based on an automatic adaptive tuning function defined as follows:

D(g 1 , g 2 ) = f (cort(g 1 , g 2 )) • δ E (g 1 , g 2 )
where f (x) is an exponential adaptive tuning function:

f (x) = 2 1 + exp(k x)
, k ≥ 0 with f (0) = 1. Figure 1 shows the adaptive tuning effect for several values of k ≥ 0. In the case FIGURE 1. The adaptive tuning effect of genes with different behavior (i.e. with cort near 0), f (x) is near 1 whatever the value of the weight k and D is approximately equal to δ E . However, in the case of genes with an opposite or similar behavior (i.e. with |cort| near to 1), the parameter k modulates the contributions of the proximity w.r. 3 Adaptive unsupervised classification for identifying of cell cycle regulated genes

HeLa Cell line Data description

This paper focuses on a specific biological events occurring during cell proliferation, this process insuring the multiplication of cells, which is drastically aberrant in cancer cells. The cell cycle, or cell division cycle, is the series of events between one cell division and the next one. The cell cycle consists of progression along four distinct phases: G 1 , S (DNA synthesis or DNA replication), G 2 and M phases. A molecular surveillance system monitors the cell's progress through the cell cycle and checkpoints ensure that a cell divides only when it has completed all of the molecular prerequisites for producing healthy daughter cells. These restriction points mark the inter-phases transitions, G 1 /S is the first such transition. The genome-wide program of gene expression during the cell division cycle aims to determine the genes well expressed during studied cell cycles [START_REF] Spellman | Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization[END_REF], [START_REF] Oliva | The cell cycle-regulated genes of Schizosaccharomyces pombe[END_REF] and [START_REF] Cho | Transcriptional regulation and function during the human cell cycle[END_REF]). Looking at the transcriptome of proliferating synchronized cells leads to the construction of gene expression profiles along time, i.e. during cell cycle progress. This application is concerned with the analysis of experimental transcriptomic data from the human Hela cell line published in [START_REF] Whitfield | Identification of Genes Periodically Expressed in the Human Cell Cycle and Their Expression in Tumors Molecular[END_REF] (http://genome-www.stanford.edu/Human-CellCycle/Hela/). Our study will focus on the 1099 genes, recorded in the third experimentation of the HeLa application. Genes are described by their expressing levels, during the cell-cycle progression, along 48 instants of times after cell synchronization.

Conventional identification of cell cycle genes

Let's illustrate the approach proposed by [START_REF] Whitfield | Identification of Genes Periodically Expressed in the Human Cell Cycle and Their Expression in Tumors Molecular[END_REF] to identify the cell cycle genes of the HeLa application. Authors consider a set of 20 reference genes characterizing the following 5 cell cycle phases and transitions: S, G 1 /S, G 2 , G 2 /M , and M/G 1 . The set of 20 genes is composed of 5 classes of 4 reference genes per phase (table 2). Figure 4 gives, for each cell cycle phase, the expression profiles of the 4 reference genes. Authors argued the selection of the 20 reference genes by their peaks expression in each cell cycle phase. Each of the 1099 studied genes is then assigned a cell cycle phase of the most similar 4 reference genes. The used similarity is based on the expression values regardless to the gene expression behavior. If we observe in detail the 20 gene expression profiles given in figure 2 we find some contradictions. First, the reference genes CDC2, CCNF, CCNA2 characterizing the G 2 phase don't peak at G 2 but at G 2 /M . Similarly, the G 2 /M reference genes BUB1 and PLK peak at M/G 1 instead at G 2 /M . These observations are supported by the annotations of Genecards database (http://www.genecards.org/) and KEGG molecular pathway database (http://www.genome.ad.jp/kegg/kegg2.html).

M <-> M <-> M <-> <---> S <--> G2 <---> G1 <------> S <-> G2 <---> G1 <------> S <--> G2 <------> G1 <-> S CCNE1 E2F1
M <-> M <-> M <-> <---> S <--> G2 <---> G1 <------> S <-> G2 <---> G1 <------> S <--> G2 <------> G1 <-> S RFC4
M <-> M <-> M <-> <---> S <--> G2 <---> G1 <------> S <-> G2 <---> G1 <------> S <--> G2 <------> G1 <-> S CDC2 TOP2A
M <-> M <-> M <-> <---> S <--> G2 <---> G1 <------> S <-> G2 <---> G1 <------> S <--> G2 <------> G1 <-> S STK15 BUB1
M <-> M <-> M <-> <---> S <--> G2 <---> G1 <------> S <-> G2 <---> G1 <------> S <--> G2 <------> G1 <-> S PTTG1 RAD21 VEGFC CDKN3 M/G1
FIGURE 2. Gene expression profiles for the 20 reference genes whose expression peaks in each phase of the cell cycle : G1/S, S, G2, G2/M and M/G1. The double arrowed lines delimit the time duration for each cell cycle phase : G1, S, G2 and M.

Adaptive unsupervised classification for identifying genes cell cycle phases

Our purpose is first to identify the cell cycle phases of the studied genes; then to determine the set of genes well characterizing each cell cycle phase. For this, we propose to use the classical partitioning around medoids (PAM) method through an adaptive approach based on the new dissimilarity index. The idea of PAM method is that in order to classify objects into nb clusters, it selects nb objects called representative objects, the clusters are then obtained by assigning each remaining object to the nearest representative object. The clue is that the representative objects must be selected so that they minimize the average dissimilarity to all the other objects of the same cluster. In this work, PAM method is preferred to the familiar k-means approach for mainly two reasons. First, it is more robust with respect to outliers, which are numerous in genes expression data due to measurement errors. Secondly, it allows a more detailed analysis of the partition by providing clustering characteristics and a graphical display (a so-called silhouette plot). In particular, the extraction of the set of genes well characterizing each cluster (i.e. each cell cycle phase) is based on the silhouette width of each classified object indicating whether an object is well classified, misclassified or lies on the boundary of a cluster. The quality of a partition is estimated by the average of the silhouette width for all the classified objects. For more details about PAM see [START_REF] Kaufman | Finding Groups in Data. An Introduction to Cluster Analysis[END_REF]. The main idea of the adaptive approach is to perform, for several values of k (k=0,...,6 per a lag of 0.01), a PAM (number of clusters=5, for the 5 cell cycle phases in interest) method on the whole 1099 genes based on the dissimilarity index D k to look for the best value of k maximizing the average silhouette width of the obtained partition. Let P k * (k*=5.9) be such a partition and figure 3 the associated silhouette plot. In the literature, there are 43 genes (about 10 genes per phase) identified as involving in the cell cycle division process [START_REF] Whitfield | Identification of Genes Periodically Expressed in the Human Cell Cycle and Their Expression in Tumors Molecular[END_REF]); therefor we extract from each cluster of P k * a FIGURE 3. Silhouette width of P k * =5.9

Comparative analysis and discussion

Let's first discuss the PAM obtained results. The optimal partition P k * maximizing the average silhouette is obtained for k*=5.9. This value means that the 5 main patterns of gene expression profiles are distinctive essentially through their behaviors (table 1). Figure 3, reveals an average silhouette (sw) of 0.33 which indicates that the clustering structure is no better than reasonable. However, if limited to the 50 kernel genes, the average silhouette coefficient is about 0.67, which means that the kernel sets are well separated from each other. Figure 3 indicates that the cluster 2 (G 1 /S) possesses the largest sw of 0.47, which means that this cluster is well separated from the other clusters, whereas the cluster 5 (G 1 ) possesses a rather very narrow sw of 0.08, which means that the cluster G 1 is not very clearly separated from the other clusters.

According to the biological knowledge, note that the well studied genes CCNE1, CCNA2, and CCNB1 known as mitotic cyclins, respectively classed into the clusters G 1 /S, G 2 /M and M/G 1 , appear in the expected biological temporal order during the cell division cycle (G 1 , S, G 2 and M ). As a support to the obtained G 1 /S cluster note that E2F1, a transcription factor known as a key regulator of cell cycle progression involved in the control of cell cycle progression from G 1 to S, exhibits the largest silhouette (sw=0.832). The genes CCNE1 (sw=0.786) and MCM6 (sw=0.812) known as activated and induced respectively by E2F1, are also classified into the G 1 /S cluster. It's shown experimentally that the CCNA2 gene promotes G 2 /M transition.

Our proposed approach classified well CCNA2 into the cluster G 2 /M (sw=0.458), whereas it is selected as a G 2 Whitfield reference gene. The gene UBE2C (sw=0.779) belonging to the G 2 /M cluster is well assessed by the biological knowledge: it represents an enzyme of the ubiquitin pathway regulating destruction of mitotic cyclins near the end of mitosis (G 2 /M transition). Finally, note all Whitfield reference genes labeled as G 2 phase are classified into the cluster G 2 /M and except STK15, all the reference genes labeled as G 2 /M are classified into the cluster M/G 1 , which clearly corroborate the contradictions discussed in the paragraph 3.2.

Conclusion

This paper focuses on an alternative adaptive clustering for the identification of the cell cycle genes. We propose a new dissimilarity index for gene expression profiles to include both proximity measures w.r.t. expressed values and w.r.t. genes expression behavior. An adaptive unsupervised classification, based on the proposed dissimilarity index, is then performed to identify the cell cycle phases of the studied genes. Finally, we propose a new well-justified set of reference genes, assessed by a published biological knowledge.

M <-> M <-> M <-> <---> S <--> G2 <---> G1 <------> S <-> G2 <---> G1 <------> S <--> G2 <------> G1

  )=2/(1 +exp(k x))

FIGURE 4 .

 4 FIGURE 4. Kernel Gene expression profiles during cell division cycle

TABLE 1 .

 1 t. values and w.r.t. behavior to the dissimilarity index D. As k increases, the contribution of the proximity w.r.t. behavior 1 -2/(1 + exp(k |cort|)) increases, whereas the contribution of the proximity w.r.t. values 2/(1 + exp(k |cort|)) decreases. For instance, for k = 0, the proximity w.r.t. behavior contributes at 0% to D whereas the proximity w.r.t. values contributes at 100% to D (the value of D is totally determined by δ E ). For k = 2, the proximity w.r.t. behavior contributes at 76.2% to D whereas the proximity w.r.t. values contributes at 23.8% to D (23.8% of the value of D is determined by δ E and the remaining 76.2% by cort). Table1summarizes, in the case of similar or opposite behavior (|cort|=1), the contributions of both, proximity w.r.t. behavior and w.r.t. values, to the dissimilarity index D. Let us add two remarks: First, if k = 0 the proposed dissimilarity index D is identical to δ E ; hence D could be considered as an extension of δ E to both, behavior and value proximity measures. The second point is that if Contribution of the proximity w.r.t. behavior and w.r.t. values to D as a function of k δ E approaches 0 (i.e., the genes expression are close w.r.t. values), then cort approaches 1 (i.e., the genes expression are similar w.r.t. behavior) and D approaches 0. Finally, we can check easily that D verifies the identity and the symmetry properties of a distance, but not the triangular inequality.

		Proximity w.r.t. behavior	Proximity w.r.t. values
		Contribution (%)	Contribution (%)
	k = 0	0	100
	k = 1	46.2	53.7
	k = 2	76.2	23.8
	k = 3	90.5	9.4
	k ≥ 5	; 100	; 0

TABLE 2 .

 2 The 20 reference genes.

		1.5
		1.0
	log(base2) expression ratio	-0.5 0.0 0.5
		-1.0
		-1.5
		0 0	2	4	6	8	10 10	12	14	16	18	20 20	22	24	26	28	30 30	32	34	36	38	40 40	42	44	46
													time(hrs)

kernel set of the 10 well-classified genes maximizing the silhouette width. Table 3 gives for each cluster the set of kernel genes (Gene Type = K). We indicate for each kernel gene it's name, it's Whitfield assignment phase, the number of the neighbor cluster and it's silhouette width (sw) indicating if it is well-classified (sw close to 1) or misclassified (sw close to -1). We indicate also the set of Whitfield reference genes (table 2) belonging to each cluster (Gene Type = R). Figure 4, visualizes for each cluster the expression profiles of the associated kernel genes. The observation of the progression of the kernel genes during the cell division cycle reveals that: kernel genes of the cluster 1 peak clearly at S phase, kernel genes of cluster 2 peak at G 1 /S, kernel genes of cluster 3 peak at G 2 /M phase, kernel genes of cluster 4 peak at M/G 1 and finally, kernel genes of the cluster 5 peak at G 1 phase. Let's remark that due to the asynchronization of cells, it's more reliable that interpretations will be limited to the earlier cell cycles. According to that, each cluster is assigned the cell cycle phase of its kernel set genes (indicated in the last column of the table 3). Finally, each of the 1099 studied genes is then assigned the cell cycle phase of the cluster it belongs to.

TABLE 3. The 50 Kernel genes (Gene Type = K) well-characterized the cell cycle phases: S, G1/S, G2/M , M/G1 and G1, with the classification of the 20 Whitfield reference genes (Gene Type = R) through the 5 obtained clusters.