Dimitrios M Vitsios

Fotis E Psomopoulos

Pericles A Mitkas

Christos A Ouzounis

Multi-genome Core Pathway Identification through Gene Clustering

Keywords: bioinformatics, metabolic pathways, clustering algorithm, phylogenetic analysis

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Metabolomics is the scientific study of chemical processes involving metabolites. Specifically, metabolomics is the "systematic study of the unique chemical fingerprints that are left behind specific cellular processes". The metabolome represents the collection of all reactants (enzymes, proteins or other chemical compounds) in a biological cell, tissue, organ or organism, which are the end products of cellular processes.

In biochemistry, metabolic pathways are series of chemical reactions occurring within a cell. In each pathway, a principal chemical is modified by a series of chemical reactions. Enzymes catalyze these reactions, and often require dietary minerals, vitamins, and other cofactors in order to function properly. Because of the many chemicals (a.k.a. "metabolites") that may be involved, metabolic pathways can be quite elaborate. In addition, numerous distinct pathways co-exist within a cell. This collection of pathways is called the metabolic network.

For the current study, pathway-related information was retrieved from KEGG Pathway Database (KEGG) [START_REF] Kanehisa | KEGG: kyoto encyclopedia of genes and genomes[END_REF], a collection of manually drawn pathway maps representing our knowledge on the molecular interaction and reaction networks for various processes and systems including "Metabolism", "Human Diseases" and "Cellular Processes" among others. The analysis of completely sequenced genomes can yield useful insight into the evolution and multi-level organization of organisms. With the current advances in genomics and proteomics, it has become imperative to explore its impact as reflected in the metabolic signature of each genome [START_REF] Zhang | Phylophenetic properties of metabolic pathway topologies as revealed by global analysis[END_REF]. To this end a methodology is presented, which applies a clustering algorithm to genes from different species participating in the same pathway.

Problem Description

Considering the fact that the majority of the known genomes from all three domains (Archaea, Bacteria and Eukaryota) use highly similar chemical processes in the form of metabolic pathways, it is a logical step to assume that the dissimilarities in the particular expressions in the metabolic pathways of different genomes can be viewed as the result of the evolution of a reference pathway existent at the common ancestor of those genomes. Specifically, considering that each pathway can be represented by a graph, then, starting from the graph of the reference pathway, the currently formed pathways have emerged essentially after the addition or deletion of some nodes and/or edges from the reference graph, while preserving the core body of the pathway virtually unchanged. In fact, by observing the graph representing a metabolic pathway using several genomes as a reference, several compact sets of edges (i.e. reactions) can be distinguished that seem to be highly conserved in the genomes, whereas at the same time there exist several sub-graphs that have been added or deleted, making evident the differentiation of the pathway expression between the genomes. This observation leads us to the conclusion that the transformations that have been carried out in the gene-level of the genomes through evolution due to mutations of the genome sequences have resulted to changes in the topology of the metabolic pathway while retaining the overall process.

Thus, the objective of the current study is to extract a meaningful clustering of the genes participating in a common metabolic pathway of several genomes by comparing the genes sequences and by evaluating the degree of homology between genes from different genomes. The result of the applied process is the extraction of cohesive gene groups that can give us information about the evolutionary similarity between the genomes examined or the uniqueness of some reactions in particular genomes. This knowledge can then be used in order to correlate some macroscopic differences between the genomes with the differences that arise between them in the pathway level.

Related work

Genes usually do not act individually but form functional or structure organizations, exemplified by metabolic pathways. As metabolic pathways are essential to the survival of organisms, and their evolution has been under debate for more than half a century [START_REF] Schmid | Metabolites: a helping hand for pathway evolution?[END_REF], a combined phylogenetic and phenetic analysis of pathway topology might expand the understanding of the evolutionary processes molding their form and structure.

Several groups have carried out phylogenetic analyses based on metabolic pathways, deriving phylogenetic trees from the information of individual pathways [START_REF] Forst | Evolution of metabolisms: a new method for the comparison of metabolic pathways using genomics information[END_REF][START_REF] Forst | Phylogenetic analysis of metabolic pathways[END_REF][START_REF] Heymans | Deriving phylogenetic trees from the similarity analysis of metabolic pathways[END_REF], the presence and absence of entire pathways [START_REF] Liao | Genome Comparisons Based on Profiles of Metabolic Pathway[END_REF], or the reaction content of entire pathways [START_REF] Hong | Phylogenetic analysis based on genome-scale metabolic pathway reaction content[END_REF]. These studies have provided valuable insight into the evolution of metabolism; however, as phylogenetic trees, they have generally diverged substantially from trees based on 16S rRNA, the most used molecule for phylogeny reconstruction. A common feature of phylogenetic trees based on metabolic information is that, owing to similar evolutionary pressures, organisms in similar habitats tend to be clustered together, and Aguilar et al. [START_REF] Aguila | Analysis of phenetic trees based on metabolic capabilites across the three domains of life[END_REF] therefore regarded such trees as phenetic rather than phylogenetic. Furthermore, one group showed that trees based on different subsets of metabolic networks were different [START_REF] Aguila | Analysis of phenetic trees based on metabolic capabilites across the three domains of life[END_REF], and another result also indicated a similar situation when several different pathways were used to construct trees separately [START_REF] Heymans | Deriving phylogenetic trees from the similarity analysis of metabolic pathways[END_REF].

On the other hand, phylogenetic profiles are commonly used in evolutionary studies, as they are based on sequence similarity. There are several recent approaches that either directly utilize phylogenetic profiles for functional prediction of gene clusters [START_REF] Lin | In silico prioritisation of candidate genes for prokaryotic gene function discovery: an application of phylogenetic profiles[END_REF] or combine them with other biological data sources for increased sensitivity [START_REF] Gupta | Automated Hierarchical Density Shaving: A Robust Automated Clustering and Visualization Framework for Large Biological Data Sets[END_REF]. A slightly different approach, and an intermediate step towards the work presented in this paper, is to produce a tree-like structure of gene clusters in order to reconstruct the evolutionary relationships between them [START_REF] Psomopoulos | Multi Level Clustering of Phylogenetic Profiles[END_REF]. However, it is shown that the output of large-scale reconstructions is notably more difficult to interpret biologically. Our approach extends this work by aiming to extract close gene associations from metabolic pathways through unsupervised clustering at a sequence level. This level of association can be enhanced if the phylogenetic relationship of the corresponding genomes is taken under consideration.

Algorithm overview

The algorithm accepts as input a KEGG pathway map identifier mapID, and a list of n genome identifiers that comprise the target data set. A specific clustering algorithm (MCL [START_REF] Enright | An efficient algorithm for large-scale detection of protein families[END_REF] or EM [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF]) is also selected for the current run and some parameters can be set i.e. the inflation parameter for MCL and the e-Value threshold for the construction of the homology matrix. We have determined through extensive experimentation the optimal value for the inflation parameter ('12') and the recommended value for the execution of the methodology (Table 1). Specifically, setting the inflation parameter to the optimal value stabilizes the number of clusters produced by the algorithm while at the same time yielding the optimal results (maximum values) with regards to the average intra-cluster similarity and homology.

By retrieving specific data (gene identifiers, FASTA sequences and the corresponding EC identifiers) about the k i genes (i = 1 … n) from each genome that participate in mapID, n blastable databases are constructed. In the next step, blast searches are performed, leading eventually to a matrix containing the phylogenetic profiles of all genes participating in the study. The homology between genes is determined using the default BLAST values and an e-Value threshold of 10 -5 . The data in this homology matrix P are then clustered using the MCL or EM algorithm and a custom similarity metric based on the jaccard metric, given from the following equation:

sim 𝑗𝑎𝑐 = 𝑚 11 +𝑚 00
𝑚 01 +𝑚 10 +𝑚 11 [START_REF] Kanehisa | KEGG: kyoto encyclopedia of genes and genomes[END_REF] where: m 11 represents the number of genomes where both genes (the distance of which is calculated) have a homologue, m 00 represents the number of genomes where both genes do not have a homologue and m 10 (m 01 ) represents the number of genomes where the first gene (second gene) has a homologue while the second one (first one) does not. After thorough experimentation using both the MCL and the EM clustering algorithms (data not shown), we have concluded that the clusters generated with MCL are more robust and more closely correlated to the biological aspect of the problem, while the clusters produced by EM cannot be readily interpreted. Generally, in many cases the EM's results tend to closely resemble with the respective results acquired with MCL. However, although they both generate almost the same number of clusters and have over 70% similarity between the "corresponding" clusters, they do not provide an equivalent level of granularity regarding the resulted clustering. Moreover, in all the experimental setups that were performed, EM cannot attain the distinction in separate groups of the genes belonging to a specific genome from the entire dataset, something that is achieved by MCL. On the other hand, in sharp contrast with EM, the MCL algorithm leads to a degenerate clustering (i.e. singleton cluster) when the organisms selected for a test case have a close phylogenetic relationship. Thus, in case resulting clusters form EM can be sufficiently interpreted from a biological aspect, then we can assume that EM is more advantageous than MCL in cases of limited phylogenetic diversity. For the current study however, the development of the overall methodology is based on the use of MCL as the genes' clustering algorithm.

Following that step, a clustering validation is performed to determine the similarity and consequently the significance of the clusters generated and the FASTA files for each of the clusters are retrieved from KEGG. The homology score for each pair of clusters i and j (at index i,j of the homology matrix) is the total number of homologues found between clusters i and j divided by the product of the total numbers of genes in clusters i and j. The homology between each pair of genes is determined by executing blast searches for all pairs of FASTA sequences corresponding to the genomes participating in the current test case and considering that an e-Value score lower than 10 -5 indicates the existence of homology between the examined genes. The similarity score for each pair of clusters i and j is the sum of the distances (Eq. 1) between all pairs of genes in clusters i and j divided by the product of the total numbers of genes in clusters i and j. The validation process is performed by comparing the intra-cluster similarity and homology to the respective inter-cluster values. In all cases (data not shown), it was confirmed that the ratio of intra-cluster to inter-cluster values was in the range of 10, with the notable exception of a few minor clusters were the ratio was closer to 1. Finally, a post-processing procedure is applied to the data of the clusters in order to extract additional information regarding the composition of the clusters and the distribution of similar genes across the genomes and the acquired results are used to retrieve the customized colored versions of the examined metabolic pathways by cluster and by genome. The flowchart of the algorithm is shown in Fig. 1.

Algorithm complexity and efficiency

The maximum data size processed using that algorithm involved 325 genomes and 10,327 genes, due to memory space limits, considering 4 GBs of available memory space. The corresponding execution time of the algorithm reached about 15 hours using a 3.0 GHz Quad-Core Processor. The overall complexity of the implemented methodology is shown in Fig. 2.

Software implementation

The implementation of the algorithm outlined in the paper incorporates several different software tools and libraries that are combined and interconnected on top of a Java-based application framework. Specifically, data retrieval from KEGG database is achieved through SOAP-based web-services provided by the KEGG API. The construction of the blastable databases and the consequent BLAST searches are performed using the BLAST libraries (version 2.2.25+) provided by NCBI, whereas the output of the BLAST runs is read and organized using a custom XML DOM parser. The gene clustering is performed using the MCL implementation provided by the author of the original paper [START_REF] Stijn Van Dongen | Graph Clustering by Flow Simulation[END_REF] while the EM's implementation used is the one provided by WEKA Suite [START_REF] Hall | The WEKA Data Mining Software: An Update[END_REF]. Finally, the application makes use of certain functions derived from the Matlab Bioinformatics toolbox [17] for the calculation of phylogenetic profiles using metabolic data in conjunction with information extracted from the resulted genes clustering. Those functions are accessed via Matlab generated Java components that are executed with MCR runtime engine [18]. The overall application is built on top of a common Java layer which is, besides the combination and coordination of the aforementioned heterogeneous tools and libraries, responsible for several pre/post-processing, evaluation and data mining tasks. Finally, the java tool that has been implemented is publicly available from the following url: http://olympus.ee.auth.gr/~fpsom/alignPaths-pkg_Bundle_v0.9.tar.gz

Results overview

The proposed method was applied on several different pathways and genome sets and validated through both statistical methods and literature reviews. A characteristic test case is presented here, namely the application of the method on the Glycolysis / Gluconeogenesis metabolic pathway (KEGG identifier map00010) which is a well known and extensively documented pathway.

Derived from the Greek stem glyk-, "sweet," and the word lysis, "dissolution", glycolysis is an ancient pathway employed by a host of organisms. It is the sequence of reactions that metabolizes one molecule of glucose to two molecules of pyruvate with the concomitant net production of two molecules of ATP. Glycolysis is an energy-conversion pathway in many organisms, and it is tightly controlled. Gluconeogenesis is the opposite pathway from glycolysis, generating glucose from non-carbohydrate carbon substrates such as lactate, glycerol, and glucogenic amino acids. Gluconeogenesis and glycolysis are reciprocally regulated.

The three genomes participating in the test case are Arabidopsis Thaliana (ath), Escherichia Coli K-12 MG1655 (eco) and Homo Sapiens (hsa), and were selected for sufficient phylogenetic diversity. The total number of genes in this dataset is 209, distributed across the three genomes as follows: ath: 105, eco: 39 and hsa: 65. The selected algorithm was MCL and the results are presented in detail in the following section.

Experimental results

The results of the test case are presented here along with some interesting observations. Specifically, the pathway images shown below depict the distribution of the genes both by cluster and by genome after the clustering is performed and provide useful information about the internal structure of a single pathway as a result of genomes evolution. In each pathway image (Fig. 3) the EC identifiers of the corresponding genome are highlighted according to the cluster their genes belong to. A special case is the EC identifiers highlighted in black; they contain genes from more than one cluster. However, it must be noted that each gene is assigned to a single cluster, whereas an EC identifier, corresponding to several genes, may in turn belong to different clusters.

Fig. 4 shows the distribution of the gene clusters across the pathway. It is interesting to note that Cluster 1 contains EC numbers corresponding to genes that constitute the main process of the pathway (core pathway), whereas the EC identifiers of the fourth cluster contain genes only from the human genome. The second cluster contains genes only from Homo Sapiens and Arabidopsis Thaliana, as opposed to Cluster 3 that contains genes only from Escherichia Coli and Arabidopsis Thaliana. Finally, there exist several cases where an EC identifier corresponds to genes that individually belong to different clusters. These cases are shown collectively as a fifth cluster, but are also highlighted in each cluster diagram when the EC identifier contains at list one gene of the specific cluster. For validation purposes, the produced clusters were evaluated using both the modified jaccard similarity metric (Eq. 1) and the gene homology. The average intracluster similarity and homology (≈ 1.813 and 0.182, respectively) was significantly higher than the average inter-cluster similarity and homology (≈ 0.479 and 0.016 respectively). Moreover, every similarity and homology value calculated within a cluster is higher than every corresponding value calculated between the clusters (Fig. 5). The aforementioned results are indicative to the significance and validity of the produced clusters.

Discussion

Although these are only preliminary results, some interesting observations can be made. The presented test cases were among several different experimental setups. In all cases however, the first cluster always contained EC identifiers along the main reaction chain of the pathway, leading to the tentative conclusion that it may correspond to the highly conserved genes. Moreover, by superimposing the highlighted pathway diagrams along the implied phylogenetic distance of the genomes, one may infer sub-chains of the pathway that have been transformed or

Fig. 1 .

 1 Fig. 1. Flowchart demonstrating the steps of the algorithm.

Fig. 2 .

 2 Fig. 2. Execution of the algorithm for different number of genes [a] and genomes [b]. In both cases the exponential character of the algorithm is evident.

Fig. 3 .

 3 Fig. 3. The Glycolysis/Gluconeogenesis pathway for the three genomes in the case study.

Fig. 4 .

 4 Fig. 4. The Glycolysis/Gluconeogenesis pathway for each of the four produced clusters, and the case of EC identifiers with genes from multiple clusters (highlighted in black).

Fig. 5 .

 5 Fig. 5. Inter-cluster and intra-cluster homology [a] and similarity [b] figure.

Table 1 .

 1 Statistical overview of several different runs of the methodology using MCL for various values of the inflation parameter. The reference pathway is the Glycolysis / Gluconeogenesis and the examined organisms are Escherichia Coli K-12 MG1655 (eco), Arabidopsis Thaliana (ath) και Homo Sapiens (hsa).

	MCL
	inflation
		2	6	10	12	16	20	24	30
	Results
	Number of clusters	2	3	3	4	4	4	2	1
	Average inter-cluster similarity	0.29	0.42	0.42	0.48	0.48	0.48	0.29	-
	Standard
	deviation of inter-	0	0.21	0.21	0.28	0.28	0.28	0	-
	cluster similarity
	Average intra-cluster similarity	1.92	1.39	1.39	1.81	1.81	1.81	1.92	0.77
	Standard
	deviation of intra-	1.53	0.41	0.41	0.86	0.86	0.86	1.53	0
	cluster similarity
	Average inter-cluster homology	0.035	0.02	0.02	0.016	0.016	0.016	0.035	-
	Standard
	deviation of inter-	0.035	0.017	0.017	0.014	0.014	0.014	0.035	-
	cluster homology
	Average intra-cluster homology	0.12	0.14	0.14	0.18	0.18	0.18	0.12	0.08
	Standard
	deviation of intra-	0.045	0.018	0.018	0.07	0.07	0.07	0.045	0
	cluster homology

evolved across the species. A thorough investigation of this problem, together with rigorous experimentation on several different sets of pathways / genomes may provide more information in these areas.


Vasileios G Kanas
email: vaskanas@upatras.gr

Evangelia I Zacharaki

Evangelos Dermatas

Anastasios Bezerianos

Kyriakos Sgarbas

Christos Davatzikos

Combining Outlier Detection with Random Walker for Automatic Brain Tumor Segmentation

Keywords: brain neoplasms, segmentation, random walks, outlier detection, k-means, magnetic resonance imaging

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Diagnosis and characterization of brain neoplasms appears of greatest importance for remedial management. The procedure of diagnosis consists of two main steps, the extraction of the pathological region in the brain (tumor segmentation), and the assessment of malignancy of brain neoplasms (tumor classification). These two steps are essential for treatment planning and future assessment. Despite various efforts and promising results in automatic brain tumor segmentation using MRI data, accurate segmentation and characterization of abnormalities is still a demanding and difficult task due to the large variability in shape, location and image intensity of brain neoplasms. Moreover, some neoplasms also distort the surrounding healthy structures or appear together with edema that changes intensity properties of the nearby region making the extraction of neoplastic area a very complex assignment. There have been considerable efforts to develop automated computer algorithms for locating neoplasms in brain using MRI data. A review of pattern recognition methods for automatic brain tumor segmentation is presented in [START_REF] Corso | Efficient Multilevel Brain Tumor Segmentation With Integrated Bayesian Model Classification[END_REF]. Other methods are based on learning a statistical distribution intensity model of the abnormal tissue types using manually labeled samples [START_REF] Verma | Multiparametric Tissue Characterization of Brain Neoplasms and Their Recurrence Using Pattern Classification of MR Images[END_REF]. Such supervised classification methods are subject to inter-but also intra-rater variability 1 . In other relevant study [START_REF] Menze | A generative model for brain tumor segmentation in multi-modal images[END_REF], an atlas of the tumoral region has been extracted in order to augment a probabilistic map of healthy tissue priors aiming at identifying tumor boundaries. They incorporate prior knowledge about the normal brain tissue, which may not be accessible for all situations. In [START_REF] Popuri | 3D variational brain tumor segmentation using Dirichlet priors on a clustered feature set[END_REF] it is propose a supervised variational method that incorporates appearance priors to better disambiguate the tumor from the surrounding deformed brain tissue. Yet, these appearance priors are not suitable for all patients. Nan Zhang et. al [START_REF] Zhang | Kernel feature selection to fuse multi-spectral MRI images for brain tumor segmentation[END_REF] integrate Support Vector Machine classification with a selection of the features in a kernel space to learn the tumor profile from the first MRI examination. Then, they proceed to refine the tumor boundaries using a region growing technique to follow up the brain tumor evolution. Though, they used a supervised learning process at the beginning of the patients' treatment.

Neoplastic tissue classification generally requires information of several MR protocols and contrasts, as T1, T1 contrast-enhanced T2, FLAIR (FL) or diffusion tensor imaging (DTI), yet not all sequences are acquired in clinical routine, especially during regular follow-up. The focus of this work is the development of an automatic method for brain tumor segmentation combining multiple MR sequences acquired in general clinical practice. The MR sequences used in this study were the following: axial 3D T1-weighted (T1), Fluid Attenuated Inversion Recovery (FLAIR), axial 3D contrast enhanced T1 contrast-enhanced (T1CE) images. We have chosen these modalities due to the complementary information they incorporate. In T1CE, the vital neoplastic tissue appears hyperintense and in FL datasets, the perifocal edema is hyperintense.

The proposed framework applies a hybrid methodology incorporating unsupervised clustering in the first step and the random walker algorithm [START_REF] Grady | Random Walks for Image Segmentation[END_REF] in the second step. The main contribution of this work is the automatic generation of seed points using a succession of procedures, which exploit the information provided by multi-modal MRI datasets. Such an approach is applicable to enhancing masses with different malignancy grades. It aims to segment normal and abnormal tissue including the nonnecrotic enhancing neoplastic tissue (ET) and the peritumoral edema with possible inclusion of neoplastic infiltration (EDM). Its application to several datasets with various tumor shapes, intensities and locations shows that it can automatically detect and segment enhancing brain neoplasms and edematous regions with a good quality.

The rest of a paper is organized as follows. In the next section, the methods incorporated on our proposed framework are explained. In section 3, the data are described and experimental results are shown. Finally the last section is devoted to some discussion and concluding remarks.

Materials and Methods

In this paper, we focus on the segmentation process of the different abnormal tissue types (neoplastic and edematous) and healthy brain tissue, i.e. white matter (WM), gray matter (GM), cerebrospinal fluid (CSF). The main purpose is to develop an automated tool to avoid the manual tracing of all tissue types, which is extremely time consuming when performed in 3D data and also highly dependent on the expertise of operators. Specifically, we develop a method based on random walker algorithm. This algorithm treats images as a graph, where user-defined seeds are required. Our contribution focuses on automatically finding seed points and segmenting intracranial neoplastic regions by combining several MRI sequences. This unsupervised approach incorporates imaging data which are acquired in a routine clinical protocol, such as multi-parametric conventional MRI. During our experiments we concluded that T1CE, FL and T1 images provide us with adequate information appropriate to our analysis, as inputs to random walker. Our approach could be separated in the following core steps (Fig. 1), which will be discussed in detail during this paper:

o Data acquisition and preprocessing.

o Automatically calculating seeds based on unsupervised clustering and outlier detection.

o Segmentation of all issue types using the random walker algorithm. Particularly, based on the hypothesis that neoplastic and edematous tissue profile deviates from the healthy tissue profile, we first applied a non-parametric distancebased outlier detection methodology to approximately segment tumor and edema. We verified our hypothesis by plotting the T1CE and FL image intensities for abnormal and normal tissue based on the expert defined ROIs, as illustrated in Fig. 2. Generally it can be noticed that healthy tissues and abnormalities overlap with each other. This makes thresholding techniques inaccurate for the segmentation process. However, it can be observed that some of the pathological tissue (for both neoplastic and edematous regions) deviates significantly from the healthy tissue. The aim is to automatically detect these voxels and use them later as seeds to the random walker algorithm in order to better determine the separation boundary between normal and abnormal tissue based on local image properties. More details on each methodological component are given next.

Data description

The data and ROIs are acquired from a previous study [START_REF] Zacharaki | Classification of brain tumor type and grade using MRI texture and shape in a machine learning scheme[END_REF] and described also here for completeness. Twenty-six patients (age 45-83 years) with a diagnosis of brain neoplasm were selected out of 97 with the criterion of the presence of tumor enhancement. All patients underwent biopsy or surgical resection of the neoplasm with histopathological diagnosis. The total of 26 brain masses were graded based on WHO criteria as metastasis (14), meningiomas grade I (3), gliomas grade II (2), gliomas grade III (2), and glioblastomas grade IV [START_REF] Zhang | Kernel feature selection to fuse multi-spectral MRI images for brain tumor segmentation[END_REF].

The images were preprocessed following a number of steps including noise reduction, bias-field correction, co-registration of all sequences (T1, T2, T1CE, FL), skull stripping and histogram matching to make the intensities comparable across patients. For further information about the preprocessing step see [START_REF] Zacharaki | Classification of brain tumor type and grade using MRI texture and shape in a machine learning scheme[END_REF].

Extraction of healthy tissue clusters

The k-means clustering is a method of cluster analysis which aims to partition n observations into k clusters in which each observation belongs to the cluster with the nearest mean. The healthy tissue priors are calculated by segmenting the patient's T1 image into WM, GM and CSF/necrosis using the k-means algorithm. Several distance metrics were tested, however we used the Euclidean distance metric because it provided more accurate segmentation results for healthy tissue. During our analysis we combined tumor necrosis with the CSF cluster since the intensity profiles are similar and because necrosis is a non-vital tissue and consequently cannot be considered as malignancy. As a note, fuzzy c-means algorithm [START_REF] Bezdek | Pattern Recognition with Fuzzy Objective Function Algorithms[END_REF] was also tested without producing more accurate results.

Outlier detection

After calculating the cluster centers of healthy tissue, abnormalities are detected as samples that are distant to all three cluster centers (for WM, GM, CSF). The Mahalanobis distance (MD) is used as distance metric, and is calculated for voxel j as follows:

(1) where is the feature representation of voxel j, is the center location of cluster i and the inverse intensity variance of the data in class i and it expresses the probability of each voxel to be an outlier. We implemented univariate analysis (p=1), e.g. T1CE and FL as two distinct feature vectors for ET and EDM were used, respectively. By calculating MD for every voxel and then thresholding the obtained MD maps for ET and EDM, we detect the neoplastic tissue and edematous tissue, as well as intensity artifacts.

The selected threshold value is not selected the same for all subjects because images might differ in contrast and brightness levels even if preprocessing is applied to factor out this variance. The threshold is automatically determined by examining the histogram of the MD maps for each subject and by selecting the 4% largest MD values for ET and EDM.

False Positive reduction

At this point, we have noticed that the thresholding technique introduces some false positives that must be eliminated. For example, the thresholding technique when applied on T1CE image detected pixels belonging to cerebellum. In order to reduce false positives we constraint the number of possible tumor clusters to equal a predefined parameter (N f ) describing the number of tumor foci (e.g. 1 or 2). The retained tumor clusters were the N f largest (regarding total volume) spatially connected components. Connected-component analysis was performed by application of the union-find algorithm [START_REF] Sedgewick | Algorithms in C[END_REF].

Final seeds acquisition

Once the seeds regions have been segmented, we selected as seeds voxels with higher likelihood to be correct detection. Those are voxels in the interior of the regions since the uncertainty is higher towards the boundaries due to the diffusive nature of many high grade tumors and due to partial volume effects. Particularly, we applied skeletonization operation on the labeled image, because the skeleton retains spatial and topological properties of the object while containing all the information necessary to reconstruct the object. Assuming the boundaries of the region are unclear, the skeletonization technique retains seeds belonging to the specific region with higher certainty. Furthermore, we implemented the same morphological operations to the segmented healthy tissue (section 2.2) in order to obtain seeds for the healthy tissues (WM, GM and CSF).

Random Walker algorithm (Theoretical background)

We propose a multi-parametric framework for automated brain tumor segmentation by extracting intensity attributes and defining a precise notion for a graph. For a given image, we can use the random walk model or circuit model to solve the segmentation problem [START_REF] Grady | Random Walks for Image Segmentation[END_REF]. The random walk method consists of computing the probability that a "random walker", starting at an arbitrary node, first reaches a pre-assigned node with a particular label. The starting node is then given the label with the highest probability. Assuming that the image is a given graph, each pixel in the graph represents a vertex of the circuit, which connects with 4 neighborhood pixels (if the algorithm is applied on a 2D slice) with weights. The total weight between pixel i and pixel j is then defined as [START_REF] Verma | Multiparametric Tissue Characterization of Brain Neoplasms and Their Recurrence Using Pattern Classification of MR Images[END_REF] where h i and h j are the intensities of pixels i, j respectively and λ is a constant penalizing dissimilarity. Low edge weights indicate high probabilities of region margin evidence between two neighboring pixels and avoid a random walker crossing these boundaries. We used 3-dimensional vectors for h i,j corresponding to T1CE, T1 and FL intensities. Equation ( 2) is a typical Gaussian weighting function and during our experiments the free parameter λ was assigned to its default value (λ=90). Results showed that the segmentation was not sensitive to the choice of λ. The steps of random walker algorithm are the following:

(i) Construct an undirected weighted graph representing the image. Each pixel corresponds to vertices of the graph and image intensities are mapped to edge weights (equation ( 2)).

(ii) Obtain automatically a set of marked nodes with K labels. In general, seeds are produced interactively by the user. However, we generate the seeds automatically, as described in the previous section. Moreover, in our case the number of labels are K=6 (ET, EDM, GM, WM, CSF, background of the image).

(iii) Solve a linear system of equations (see [START_REF] Grady | Random Walks for Image Segmentation[END_REF] for more details) to obtain probabilities for the unmarked nodes.

(iv) Assign to each node the label corresponding to maximum probability obtaining a final segmentation of the image.

Method evaluation

The results of our automated method are compared with manually segmented ROIs. Manual segmentation was performed by experts by tracing hyperintensities in FL images and labeling them as EDM (including possible tumor infiltration). Hyperintensities in T1CE were marked as ET. The automatic segmentation of edematous and neoplastic regions is assessed using the Dice coefficient as measure of spatial overlap. The Dice coefficient (DSC) over two segmentations, L1 and L2, is the volume of intersection of the two segmentations divided by the mean of the two segmentation volumes:

(3)

3 Experimental Results

Segmentation of enhancing tumor, edematous region and healthy tissue

Fig. 3 and Fig. 4 illustrate some examples of enhancing tumor and edema segmentation by the proposed method and Fig. 5 demonstrates the quantitative segmentation assessment of our method for both ET and EDM tissues. Overall good agreement with manual segmentation is observed. In general, the results shown in Fig. 3 and Fig. 4 are encouraging and prove the efficiency of the method to segment particular tumor tissues. It's noteworthy that this method can utilize the characteristics of the enhancing neoplastic regions and take more general and complete pathologic information into account if different weighting functions or features are chosen. A limitation of random walker algorithm is observed in situations, where the number of unmarked regions exceeds the number of seed/label groups. In such cases ambiguous segmentation results could occur. The probability inside an unseeded region could be defined as the average of its neighboring regions, weighted by the contrast between the re-gions and the level of mutual surface area. Therefore, an unmarked region that shares an equal surface area with two labeled regions will be assigned to the region for which it has a lower contrast. On the other hand, if an unmarked region has the same contrast with two labeled regions, it will be assigned to the seeded region with which it shares a greater boundary. However, if the contrast and the overlapping surface area of an ambiguous region with two seeded regions are equal, the ambiguous region would be divided in half with respect to the two labels [START_REF] Grady | Random Walks for Image Segmentation[END_REF]. The above hypothesis becomes clear, if we examine patient 10, where necrotic tissue is surrounded by the ET region. Such a case is plausible since we have not taken into account necrosis as a separate class; thus the calculation of DSC is affected.

In addition, observing either T1CE or FL sequences, we notice that in most cases the ROIs are only approximately drawn by the experts and do not follow the exact boundaries visible in the MR images. This might be due to either insufficiency of the incorporated modalities to accurately reflect boundaries or due to incorporation of additional knowledge by the expert during manual segmentation. For example experts might better appreciate the biophysical tissue changes. Recent approaches try to predict these changes by modeling tumor growth and infiltration instead of only looking at image intensity characteristics [START_REF] Gooya | Deformable Registration of Tumor Invaded Images Using EM Algorithm and Glioma Growth Modeling[END_REF].

Conclusion

The segmentation of brain neoplasms is a complex task. Our system attempted to segment the enhancing tumor area, the region that appears bright on T1CE sequences, and edematous area, the region that appears bright on FL sequences. Note that it is not sufficient to simply threshold T1CE images by brightness because other tissues can have the same range of intensities. Here, we present an unsupervised hybrid method, combining outlier detection with random walker. The proposed method includes two steps: 1) The automatic generation of seeds and 2) implementation of multiparametric random walker algorithm. Firstly, abnormal regions were extracted using an outlier detection process. Secondly, reduction of false positives was implemented using con-nected components analysis resulting in the construction of ET and EDM components. The seeds for abnormal and normal tissue were obtained by retaining voxels with high certainty level. Finally, multilabel multiparametric random walker was applied and the segmentation results were validated against manual segmentation. The results showed that our method can produce very satisfactory segmentation outcomes while being fully automated. In this work, we used the typical Gaussian weighting function and investigated only intensity attributes (equation ( 2)). The analysis in this paper concerns only enhancing tumors with possible edema; Consequently, in future, we intend to extract more subtle attributes using more advanced methods, such as wavelet transform or shape features to exploit the information provided by MRI sequences.

Fig. 1 .

 1 Fig. 1. Proposed framework for brain enhancing tumor (ET) and edema (EDM) segmentation

Fig. 2 .

 2 Fig. 2. T1CE and FL intensity distribution for healthy tissue (WM, GM, CSF), enhancing tumor (ET) and non-enhancing (NET) or edematous tissue (EDM)

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Four patients diagnosed with brain neoplasms. From top to bottom: T1CE sequence, FL sequence, ROIs traced by expert (yellow: EDM, light blue: ET),

The same rater selecting different samples in two different time points

Acknowledgments

This research was supported by a Marie Curie International Reintegration Grant within the 7th European Community Framework Programme and an NIH grant R01 NS042645.


Evangelia I Zacharaki

Angeliki Skoura

Li An

Desmond J Smith

Vasileios Megalooikonomou

Using an Atlas-Based Approach in the Analysis of Gene Expression Maps Obtained by Voxelation

Keywords: gene expression maps, brain mapping, voxelation, gene function, gene ontology

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The mammalian brain is a complex organ exhibiting a rich variety of gene expression patterns across a broad range of cell types. Expression of genes is manifested by the production of RNA transcripts within cells and recent advances in the quantitative detection of mRNA and proteins on a genomic scale permit localization of gene products onto maps of the brain. Concerning that about 40% of the proteins encoded in eukaryotic genomes are proteins of unknown function, a challenging issue in the field of analysis of gene expression maps (GEMs) is their association with gene function information in order to reveal functional characterization of unannotated genes [START_REF] Horan | Annotating genes of known and unknown function by large-scale coexpression analysis[END_REF]. Preliminary results concerning the amounts of transcripts detected within different tissues or the same tissues under different states (e.g., physiological or disease), which are visualized using GEMs, shed light on the etiology and pathology of neurological diseases and may lead to the discovery of biomarkers and of proteins responsible for diseases [START_REF] Zhang | Region-specific protein abundance changes in the brain of MPTP-induced Parkinson's disease mouse model[END_REF].

The methodologies for spatial mapping of proteins and transcripts in the brain are various, offering compelling information. A cheap and fast method is voxelation [START_REF] Chin | A genome-scale map of expression for a mouse brain section obtained using voxelation[END_REF]. Voxelation allows acquisition of both transcript and protein mapping data in parallel simplifying co-registration of multiple genes, however it offers gene expression maps of intermediate resolution. According to this approach, the brain is divided into spatially registered voxels (cubes) and using microarrays or mass spectroscopy spatial images with quantitative information on transcripts or proteins are reconstructed.

The analysis of GEMs involves the application of feature extraction techniques combined with data mining methodologies such as clustering, classification and similarity search. Gene information from other aspects, such as Gene Ontology, is usually employed to validate biological hypothesis or to strengthen the fidelity of research outcomes. For example, aiming at the identification of unannotated genes An et al. [START_REF] An | Analysis of multiplex gene expression maps obtained by voxelation[END_REF] analyzed GEMs by extracting wavelet features and by using a multiple clustering technique. The authors confirmed the hypothesis that a subset of genes with similar expression maps display function similarity, where the identification of function similarity was based on Gene Ontology. Regarding the analysis of GEMs for the discrimination between normal and disease, Brown et al. [START_REF] Brown | Multiplex three-dimensional brain gene expression mapping in a mouse model of Parkinson's disease[END_REF] investigated the expression differences between normal and Parkinson's disease (PD) mice brains using voxel expression maps of 9000 genes acquired by voxelation. The analysis was based on two genexgene matrices (normal and PD matrix) whose elements represented the cross correlation of corresponding GEMs. Gene expression in normal and PD brains revealed significant global expression differences when averaged across voxels for known genes. Moreover, the singular Value Decomposition method was applied to the two matrices and global shifts of gene expression between normal and PD brains were indicated. Concerning spatial differentiation through a GEM produced by voxelation, a clustering analysis of gene expression patterns from mouse voxelation data revealed four distinct groups of genes corresponding to different mapping patterns [START_REF] Park | Dual spatial maps of transcript and protein abundance in the mouse brain[END_REF].

In this paper, we extend previous work [START_REF] An | Analysis of multiplex gene expression maps obtained by voxelation[END_REF] exploring gene expressions differences with regard to brain anatomy. Our aim is the identification of genes whose expressions display similar anatomical distribution regarding specific brain regions such as white matter, gray matter and the hippocampal region. We also want to investigate if the gene clusters with similar expression patterns have also similar gene function. Finally, we examine if we can extract more informative clusters by down-weighting inconsistent measurements, such as in voxels with high partial volume effect. Our investigation concludes that clusters of genes with similar localized expression patterns display functional similarity indicating that our work has the potential of creating comprehensive atlases of gene and protein expression in the mammalian brain.

Methods

The association of gene expression in the brain anatomy with functional activity can provide a better understanding of the role of the gene's products. In this study we are investigating the hypothesis that genes with similar expression maps have similar gene functions. For this purpose GEM's similarity is calculated as in previous work [START_REF] An | Analysis of multiplex gene expression maps obtained by voxelation[END_REF] based on the expression patterns acquired with the voxelation technique. Voxelation data however have much lower resolution than single cell resolution data, thus suffer from partial volume effect, meaning that the acquired expression values represent an average over the gene expression of all cells in each voxel. This limitation becomes especially prominent in regions where different tissue types mix, whereas in homogenous regions where similar expression patterns are expected, averaging does not alter significantly the gene's expression profile. Thus, in this study we examine whether partial volume effect and unreliable measurements can affect GEMs similarity and therefore alter the relationship between gene expression and function. The idea is that measurements on untrustworthy locations should have less effect on the calculation of the overall similarity between genes. Such regions include the background and ventricles and also voxels with high partial volume effect. Next we describe the construction of four spatial maps: three of them represent spatial distribution of different tissues and one map reflects the confidence on measurements. We furthermore explain how these spatial maps are used in the calculation of similarity between GEMs.

Brain partitioning and confidence map

Let Ω be the set of genes and 𝐱 𝑖 ∈ 𝑅 𝑛 , 𝑖 ∈ Ω, be the expression profile of gene i,

where 𝑛 is the number of voxels of the particular slice (at the level of the striatum) of mouse brain we consider. We explore the brain's anatomical morphology by mapping a mouse brain atlas [START_REF]C57BL/6J Coronal[END_REF] on the space of the GEM as illustrated in Fig. 1 (left). Then the registered atlas image is partitioned into three regions: (i) gray matter (GM) in the cerebral cortex and anterior cingulated area, (ii) white matter (WM) including the striatum and caudoputamen and (iii) hippocampal region (HR) including the nucleus accumbens, substantia innominata, diagonal band nucleus and medial septal nucleus and excluding the lateral septal nucleus. The three brain segments are visualized in Fig. 1 (in the midde) and are used to construct spatial maps by assigning the value of 1 to a voxel if it belongs to the corresponding brain segment, or 0 if it doesn't belong. Voxels on region boundaries are assigned a value equal to the partial volume in each tissue compartment. The acquired spatial maps are denoted as w GM , w WM , w HR ∈ 𝑅 𝑛 for GM, WM and HR, respectively. The amount of partial volume for each brain voxel 𝑗 is then calculated by the following measure of fuzziness:

𝑤 𝑃𝑉 (𝑗) = �1 -�𝑤 𝐺𝑀 2 (𝑗) + 𝑤 𝑊𝑀 2 (𝑗) + 𝑤 𝐻𝑅 2 (𝑗)� (1)
It is easy to see that the more equally distributed is the tissue to the three compartments, the higher is w PV (𝑗). The uncertainty map is calculated by averaging the amount of partial volume and the volume outside brain tissue (background or ventricular regions). A confidence map, 𝐰 𝐂 ∈ 𝑅 𝑛 , indicating the certainty of each voxel value, is then defined as the negative of the uncertainty map as illustrated in Fig. 1 on

𝑤 𝐶 (𝑗) = 1 - 𝑤 𝑃𝑉 (𝑗)+ 1-�𝑤 𝐺𝑀 (𝑗)+𝑤 𝑊𝑀 (𝑗)+𝑤 𝐻𝑅 (𝑗)� 2 (2)
the right and defined in equation 2.

Next, due to the inherent bilateral symmetry of the mouse brain and lack of "handedness" or speech-centers in mice, we decrease the amount of data by retaining the voxels of only one brain hemisphere. Similarly the GEMs, left and right hemispheres are averaged and only one hemisphere is retained, thus decreasing noise. For all GEMs and spatial maps the number of voxels is therefore reduced to 𝑛 = 42.

Definition of gene similarity of expression and function

The gene expression maps similarity between two genes, 𝐱 1 and 𝐱 2 , is defined as the squared weighted Euclidean distance function, D, formalized below:

D(𝐱 1 , 𝐱 2 ) = ∑ 𝑤(𝑗)�𝑥 1 (𝑗)-𝑥 2 (𝑗)� 2 𝑛 𝑗=1 ∑ 𝑤(𝑗) 𝑛 𝑗=1 𝑛 (3)
The weight vector w is used to emphasize dissimilarity on selective spatial locations. By incorporating such a weight vector, we can test two hypotheses. First, we can investigate whether by down-weighting the measurements on locations with high uncertainty we can form a more informative similarity measure that is not affected by partial volume and artifacts due to ventricles. This hypothesis is tested by using the confidence map w C as weight vector. The second hypothesis is that the gene function might correlate with gene expression in specific anatomic locations. Thus we investigate whether genes with similar expression in some anatomic locations have similar gene functions. For this purpose we use the three spatial maps w GM , w WM , w HR as weight vectors in equation ( 3) and investigate each region independently.

The gene function similarity is calculated using Lin's method [START_REF] Lin | An Information-Theoretic Definition of Similarity[END_REF] to evaluate function distance in Gene Ontology structure. The similarity values are obtained within each of the three categories of Gene Ontology (GO version: January 2009), and are based on frequencies from the Mouse Genome Informatics (MGI) annotation dataset (MGI version: 01/31/2009). The three categories of gene function refer to "Cellular Component", "Molecular Function" and "Biological Process". The function similarity between two clusters of genes is calculated as average pairwise function distance, as explained in [START_REF] An | Analysis of multiplex gene expression maps obtained by voxelation[END_REF]. The significance of the function distance (p-value) is calculated as a percentile in respect to average pairwise function distances. Due to the huge number of possible gene combinations, groups of genes are randomly selected and the corresponding average distances are calculated. The p-value then indicates how small is the respective function similarity with respect to the average function distances. More details are provided in [START_REF] An | Analysis of multiplex gene expression maps obtained by voxelation[END_REF].

Clustering analysis

For each spatial map, clustering analysis was performed by two sets of experiments.

In the first set we selected prototype genes with diverse expression patterns as queries and detected genes with similar expression maps. Then we calculated the average function distance in each cluster of similar genes. Similarity was assessed by the pvalue; smaller p-values result in smaller clusters. In the second set of experiments we attempted to find clusters that have both similar GEMs and similar gene functions. First clusters of GEMs were determined by the k-means algorithm using the weighted Euclidean distance function (Eq. 3). According to this criterion, the clusters consisted of genes with similar expression in the investigated region of interests (all 4 spatial maps were tested). Then only the clusters with significant expression maps similarity and average function similarity were retained, whereas the rest of the clusters were further split into an increasing number of smaller clusters until they reached the significance threshold (p-value = 0.05) for both gene expression and function. Thus the parameter K in the k-means algorithm (representing number of clusters) was not predefined, but calculated in a hierarchical fashion [START_REF] An | Analysis of multiplex gene expression maps obtained by voxelation[END_REF]. The first set of experiments will help us investigate whether the average function distance for each group of genes is reduced when specific spatial maps are used. The results of the second set of experiments will be used to extract the common gene expression patterns for each significant cluster and examine whether these patterns are related with specific anatomical locations. Moreover, connectivity relations might be revealed, if distinctive expression patterns will be identified in locations different from the applied spatial maps.

The validation of clustering is performed by the commonly used ratio of intercluster distance (D inter ) to intra-cluster distance (D intra ). The intra-cluster distance is defined as the average distance of each point to its cluster centroid,

𝐷 𝑖𝑛𝑡𝑟𝑎 = 1 𝑁 ∑ ∑ �𝑥 𝑗 -𝜇 𝑖 � 2 𝑥 𝑗 ∈𝑆 𝑖 𝑘 𝑖=1
, where N is the total number of data points, 𝑆 𝑖 , 𝑖 = 1,2, …, 𝑘, 𝑘 is the number of clusters and 𝜇 𝑖 is the centroid of the cluster 𝑆 𝑖 ). The inter-cluster distance is the minimum of the distances between each pair of cluster centroids,

𝐷 𝑖𝑛𝑡𝑒𝑟 = min (�𝜇 𝑖 -𝜇 𝑗 � 2 , 𝑖 = 1,2, … , 𝑘 -1, 𝑗 = 𝑖 + 1, … , 𝑘).
Fig. 2. Gene expression maps of the selected prototype genes 3 Results

Clustering based on prototype genes

Fig. 2 shows the gene expression maps for 6 genes selected as queries. PPP1r1b is strongly expressed in striatum, Ndn and HSLOH11 are expressed in hypothalamus, Serpinb1a is weakly expressed in striatum, Nfix is expressed in a gradient pattern in cortex and Pbx3 is expressed in striatum and adjacent ventral structures. For each prototype gene, we detected increasing number of genes (7, 15, 23,..., 78) with similar expression maps based on the use of different spatial maps and we calculated the average function distance in the group. The function similarity was considered with respect to the three function categories, Cellular Component, Molecular Function and Biological Process. Some indicative results are shown for the prototype genes PPP1r1b and Nfix in Tables 1 and2, correspondingly. We highlight p-values that are smaller than 0.05.

A gene example demonstrating the importance of the use of the confidence map w C is presented in Table 1. It can be seen the p-values regarding Biological Process are much smaller when applying the confidence map in the computation of the similarity of GEMs leading to the conclusion that the expression map of the gene Pbx3 is affected by partial volume effects and artifacts due to background and ventricles. The use of the confidence map which reflects the down-weighting of unreliable regions improves the similarity rate between the query gene Pbx3 and retrieved similar genes, confirming the hypothesis that a more informative similarity measure is obtained when incorporating weights of significance in the calculation of GEM's similarity.

Regarding the hypothesis that genes expression in specific anatomic regions might be correlated with similar gene functions, we calculated the gene expression map similarity using each one of the three spatial maps w GM , w WM , w HR and using no mask. In most cases the use of a spatial mask provides better function similarity (i.e. smaller p-value) results compared to the absence of any mask. Among the three spatial maps, each gene is associated mainly with one of them; for example PPP1r1b displays higher function similarity when focusing on WM region, whereas Nfix displays higher results when focusing on GM region. Furthermore, the comparison with the function similarity results when no spatial map is used, revealed that the Table 1. Comparison of similarity results of the gene Pbx3 with respect to the use of the confidence map w C .

Number of similar genes

Using the confidence map 𝐰 𝐂 Without using a confidence map [START_REF] Lin | An Information-Theoretic Definition of Similarity[END_REF] for the gene PPP1r1b.

The clusters of similar GEMs are created using the proposed spatially-oriented methodology focusing on WM region (left) and without using a spatial map for comparison (right).

proposed spatially-oriented approach of GEMs achieves not only higher function similarity (e.g. regarding the biological process for PPP1r1b) but also revealed new function categories that are related with this gene (e.g. cellular component and biological process for PPP1r1b). These comparative results regarding PPP1r1b gene are illustrated in Table 2. Obviously, focusing on the anatomic region of WM the retrieved genes display much higher function similarity indicating that specific anatomic regions of the brain play important role in identification of gene function.

Number of similar genes

Using the spatial map 𝐰 𝐖𝐌 Without using a spatial map

Clustering based on all genes

The GEMs of each significant cluster obtained by the hierarchical k-means algorithm for different function categories are averaged and some of them are illustrated in Fig. 3. The average maps, each of them representing one cluster, are shown only for the first 6 significant clusters for each one of the three function categories when the proposed confidence map was used (right column of Fig. 3) and when no spatial map was used in the clustering process (left column of Fig. 3). We also searched for significant clusters with low p-value of functions distance in any one of the three function categories; the corresponding results are shown in the last row of 3. The cardinality of each cluster is presented above the corresponding average GEM. Comparing the average GEMs in Fig. 3, we remark that many significant clusters are present regardless the use of the confidence map; a characteristic example is the cluster 2 in biological process without any spatial map and the cluster 4 in biological process with the use of confidence map. Although red-like pixels of GEMS indicate strong positive expression while blue-like pixels indicate strong negative expression, both of them are useful as encode important information for the analysis of gene expression. Furthermore, the location of such a pixel plays an important role; for example a red pixel in the region of ventricles is not informative. As it can be seen in Fig. 3, the clusters when using the confidence map contain more informative pixels, i.e. pixels which both display strong gene expression and are localized in meaningful regions within the brain. Table 3 shows the clustering validity score (D inter / D intra ), with the highest validity score for each function category highlighted. The results indicate that when function similarity for each of the three function categories is sought independently, the best clustering is achieved when the confidence map is used. When function similarity for any of the three function categories is aimed, the best clustering is achieved when the gene expression in GM is considered.

Discussion and Conclusions

Clustering analysis on voxelation data showed that the group of genes that was identified as similar to a target gene shares very similar gene functions in at least one gene function category. Moreover for some genes when a spatial map was used in the calculation of GEMs similarity, the average function similarity was increased or a new function category was revealed. By clustering GEMs of genes with known and unknown function together, the proposed approach has the potential to be used in predicting unknown gene functions.

Fig. 1 .

 1 Fig. 1. Mouse brain partitioning. Left: Brain atlas with superimposed voxelation grid (red voxels indicate background). Middle: Brain tissue maps based on atlas [7] at bregma=0 (green: GM, blue: WM, red: HR). Right: Confidence map w C (the darker, the less certainty).

Table 2 .

 2 Average Function Distance (p-values) using Lin's method

		Cellular	Molecular	Biological	Cellular	Molecular	Biological
		Component	Function	Process	Component	Function	Process
	7	1.00	0.00	0.00	1.00	0.00	0.60
	15	1.00	0.00	0.00	1.00	0.08	0.38
	23	1.00	0.00	0.00	1.00	0.00	0.00
	31	0.88	0.00	0.00	0.98	0.00	0.04
	39	0.24	0.00	0.00	0.86	0.00	0.00
	47	0.16	0.00	0.02	0.59	0.00	0.06
	55	0.02	0.00	0.00	0.30	0.00	0.82
	63	0.76	0.00	0.32	0.11	0.00	0.73
	70	0.59	0.00	0.92	0.02	0.00	0.94
	78	0.04	0.00	0.97	0.07	0.00	0.98

Acknowledgments

This research has been co-financed by the European Union (European Social Fund -ESF) and Greek national funds through the Operational Program "Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF) -Research Funding Program: Τhales. Investing in knowledge society through the European Social Fund.

Confidence Map


Antonis Lambrou
email: a.lambrou@cs.rhu.ac.uk

Harris Papadopoulos
email: h.papadopoulos@frederick.ac.cy

Ilia Nouretdinov
email: i.nouretdinov@cs.rhu.ac.uk

Alexander Gammerman
email: a.gammerman@cs.rhu.ac.uk

Reliable probability estimates based on Support Vector Machines for large multiclass datasets

Keywords: Support Vector Machine, well calibrated probabilities, multiclass, Inductive Venn Predictor, Machine Learning

Venn Predictors (VPs) are machine learning algorithms that can provide well calibrated multiprobability outputs for their predictions. The only drawback of Venn Predictors is their computational inefficiency, especially in the case of large datasets. In this work, we propose an Inductive Venn Predictor (IVP) which overcomes the computational inefficiency problem of the original Venn Prediction framework. Each VP is defined by a taxonomy which separates the data into categories. We develop an IVP with a taxonomy derived from a multiclass Support Vector Machine (SVM), and we compare our method with other probabilistic methods for SVMs, namely Platt's method, SVM Binning, and SVM with Isotonic Regression. We show that these methods do not always provide well calibrated outputs, while our IVP will always guarantee this property under the i.i.d. assumption.

Introduction

Support Vector Machines (SVMs) [START_REF] Vapnik | The nature of statistical learning theory[END_REF] are widely used in the field of Machine Learning for classification or regression analysis. To date, several efforts have been made in order to map the unthresholded SVM outputs into probability estimates. Some of these methods are Platt's method [START_REF] Platt | Probabilistic outputs for Support Vector Machines and comparisons to regularized likelihood methods[END_REF], SVM binning [START_REF] Drish | Obtaining calibrated probability estimates from Support Vector Machines[END_REF], and SVM with Isotonic Regression [START_REF] Zadrozny | Transforming classifier scores into accurate multiclass probability estimates[END_REF]. Nevertheless, there is no guarantee provided that the probability estimates produced by these methods will always be well calibrated. In fact as our experiments show, they can become quite misleading.

In this work, we develop a Venn Predictor (VP) based on the SVM classifier in order to produce probability estimates that are always guaranteed to be well calibrated. Venn Prediction is a novel machine learning framework that can be combined with conventional classifiers for producing well calibrated multiprobability predictions under the i.i.d. assumption. In [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF], the Venn Prediction framework is described thoroughly and a proof of the validity of its probabilities is given.

In order to overcome the computational inefficiency problem of the original Venn Prediction approach, which renders it not suitable for application to large datasets, we propose an Inductive Venn Predictor (IVP) based on the idea of Inductive Conformal Prediction. As it was shown in many studies, see e.g. [START_REF] Papadopoulos | Inductive Conformal Prediction: Theory and application to Neural Networks[END_REF][START_REF] Papadopoulos | Inductive Confidence Machines for Regression[END_REF][START_REF] Papadopoulos | Qualified predictions for large data sets in the case of pattern recognition[END_REF], Inductive Conformal Predictors are as computationally efficient as the conventional algorithms they are based on. The same is true for the proposed IVP, which is based on SVMs for multiclass tasks. We experiment on two multiclassification datasets, the Car Evaluation [START_REF] Bohanec | Knowledge acquisition and explanation for multi-attribute decision making[END_REF] and the Wine Quality [START_REF] Cortez | Modeling wine preferences by data mining from physicochemical properties[END_REF] datasets, which are freely available at the University of California, Irvine (UCI) machine learning repository [START_REF] Frank | UCI machine learning repository[END_REF]. We compare our method with Platt's method, SVM Binning, and SVM with Isotonic Regression. We demonstrate that these methods do not always provide well calibrated results, while our method can always guarantee this property under the i.i.d. assumption.

The rest of the paper is structured as follows. In section 2, we outline related work that has been conducted for estimating probabilities. In section 3, we describe the Venn Prediction framework, propose the Inductive version of the framework, and explain the taxonomy we used with multiclass SVM. In section 4, we detail our experimental settings and the obtained results. Finally, in section 5, we give our conclusions and future plans.

Related Work

In this section, we provide related work that has been conducted on methods that convert the unthresholded output f (x i ) of the SVM decision rule into a probability estimate. Hereon, f (x i ) will also be referred as the SVM score of the example x i . We examine Platt's method [START_REF] Platt | Probabilistic outputs for Support Vector Machines and comparisons to regularized likelihood methods[END_REF], SVM binning [START_REF] Drish | Obtaining calibrated probability estimates from Support Vector Machines[END_REF], and SVM with Isotonic Regression [START_REF] Zadrozny | Transforming classifier scores into accurate multiclass probability estimates[END_REF]. Moreover, we describe the approach we have followed for extending the binary SVM into multiclass SVM.

Platt's method

Platt introduced a method in [START_REF] Platt | Probabilistic outputs for Support Vector Machines and comparisons to regularized likelihood methods[END_REF] to estimate posterior probabilities based on the decision function f by fitting a sigmoid:

P (Y j = 1|f (x i )) = 1 1 + exp(Af (x i ) + B) , (1)
where Y j ∈ {-1, 1}. The best parameters A and B are determined so that they minimise the negative log-likelihood of the training data. Platt uses a Levenberg-Marquardt (LM) optimisation algorithm to solve this. As indicated in [START_REF] Platt | Probabilistic outputs for Support Vector Machines and comparisons to regularized likelihood methods[END_REF], any method for optimisation can be used. In this work, we use an improved implementation of Platt's method which uses Newton's method with backtracking for optimisation. Further details of this approach are described in [START_REF] Hsuan-Tien Lin | A note on Platt's probabilistic outputs for Support Vector Machines[END_REF].

SVM Binning

The SVM binning method [START_REF] Drish | Obtaining calibrated probability estimates from Support Vector Machines[END_REF] sorts the training examples according to their SVM scores, and then divides them into b equal sized sets, or bins, each having an upper and lower bound. Given a test example x i , it is placed in a bin according to its SVM score. The corresponding probability P (Y j = 1|x i ) is the fraction of positive training examples that fall within that bin.

There is no imposed lower or upper bound on SVM scores. Therefore, when using this method it is possible for some scores from the test examples to fall below or above the low and high scores, respectively, of the training examples. If this happens the corresponding probability P (Y j = 1|x i ) is that of the nearest bin to the score of x i .

Isotonic Regression

Isotonic regression has been used in order to map the SVM scores into probability estimates in [START_REF] Zadrozny | Transforming classifier scores into accurate multiclass probability estimates[END_REF]. An isotonic function has a monotonically increasing trend. If the scores of the SVM are ranked correctly, we can assume that the probability P (Y j = 1|x i ) will be increasing as the SVM scores increase. Therefore, we can use isotonic regression to map SVM scores into probability estimates. The most common algorithm used for isotonic regression is the Pair-Adjacent-Violators (PAV) algorithm.

The algorithm learns the probability estimate g(x i ) for each ranked example x i . First, we set g(x i ) = 1 if x i is a positive example, and g(x i ) = 0 otherwise. If g is already isotonic the function has been learned. Otherwise, there must be an example where g(x i-1 ) > g(x i ). The two examples x i-1 and x i are called pair-adjacent violators, because they violate the isotonic assumption. The values of g(x i-1 ) and g(x i ) are then replaced by their average, such that their values no longer violate the isotonic assumption. This process is repeated until an isotonic set of values is obtained. In the end, we have a list of probability estimates together with the adjacent SVM scores of the training examples. When a new example arrives, we assign the mapped probability estimate based on the score that x i has obtained from the SVM decision rule. Normally, there will be intervals of scores with the same probability estimates. Since there are no imposed boundaries on the SVM scores, the lowest interval begins from -∞ and the highest interval ends at +∞.

Multiclass SVM

The original SVM works only for binary classification problems. In this work we apply the one-against-all procedure [14] to extend the SVM for multiclass tasks. In one-against-all, we train a binary SVM classifier for each class using as positives the examples that belong to that class, and as negatives all other examples. We then convert the SVM scores of each classifier into probability estimates based on the methods described in the previous subsections, and then we combine the binary probability estimates to obtain multiclass probabilities.

The probabilities are combined by finding the probability P (Y j = 1|x i ) of each class j = 1, . . . , c and then by normalizing the probabilities of all classes to 1. The largest probability is then used to classify the example.

Venn Prediction

Venn Prediction has been introduced in [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF] where the interested reader can find a more detailed description of the framework. Since then VPs have been developed based on k-Nearest Neighbours [START_REF] Dashevskiy | Predictions with confidence in applications[END_REF], Nearest Centroid [START_REF] Dashevskiy | Reliable probabilistic classification and its application to internet traffic[END_REF] and Neural Networks [START_REF] Papadopoulos | Reliable probabilistic prediction for medical decision support[END_REF]. Furthermore, a VP based on a binary SVM has been developed in [START_REF] Zhou | A comparison of Venn Machine with Platt's method in probabilistic outputs[END_REF], and has been compared with Platt's method in the batch setting.

Typically, we have a training set1 of the form {z 1 , . . . , z n-1 }, where each z i ∈ Z is a pair (x i , y i ) consisting of the object x i and its classification y i . For a new object x n , we intend to estimate its probability of belonging to each class Y j ∈ {Y 1 , . . . , Y c }. The Venn Prediction framework assigns each one of the possible classifications Y j to x n and divides all examples {(x 1 , y 1 ), . . . , (x n , Y j )} into a number of categories based on a taxonomy. A taxonomy is a sequence A n , n = 1, . . . , N of finite measurable partitions of the space Z (n) × Z, where Z (n) is the set of all multisets of elements of Z of length n. We will write A n ({z 1 , . . . , z n }, z i ) for the category of the partition A n that contains ({z 1 , . . . , z n }, z i ). Every taxonomy A 1 , A 2 , . . . , A N defines a different VP. In the next subsection, we define a taxonomy based on the output of the SVM.

After partitioning the examples into categories using a taxonomy, the empirical probability of each classification Y k in the category τ new that contains (x n , Y j ) will be

p Yj (Y k ) = |{(x * , y * ) ∈ τ new : y * = Y k }| |τ new | . ( 2
)
This is a probability distribution for the class of x n . So after assigning all possible classifications to x n we get a set of probability distributions P n = {p Yj : Y j ∈ {Y 1 , . . . , Y c }} that compose the multi-probability prediction of the VP. As proved in [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF], these are automatically well calibrated, regardless of the taxonomy used. The maximum and minimum probabilities obtained for each label Y k amongst all distributions {p Yj : Y j ∈ {Y 1 , . . . , Y c }}, define the interval for the probability of the new example belonging to Y k . We denote these probabilities as U (Y k ) and L(Y k ), respectively. The VP outputs the prediction ŷn = Y k best , where

k best = arg max k=1,...,c p(k), (3)
and p(k) is the mean of the probabilities obtained for label Y k amongst all probability distributions. The probability interval for this prediction is

[L(Y k ), U (Y k )].

Inductive Venn Prediction

The transductive nature of the original Venn Prediction framework is computationally inefficient, since it requires training the underlying algorithm for every possible class of each new test example. To address this problem we follow the idea of the Inductive Conformal Prediction, and propose an efficient Inductive Venn Predictor (IVP). Our approach splits the available training examples into two parts, the proper training set and the calibration set. We then use the proper training set to train the underlying algorithm and the calibration set to calculate the set of probability distributions for each new example. Specifically, on each step of the algorithm in the online mode, we make a Venn Prediction analogue to a step of the Inductive Conformal Prediction ( [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF], p.98). For each number of available training examples n-1, we select q ≤ n-1 examples to form the training set for the SVM classifier and use the remaining examples as the calibration set. For the taxonomy the training examples z 1 , . . . , z q are considered as fixed parameters. The original taxonomy function A is transformed to another taxonomy A such that A n-q ({z q+1 , . . . , z n }, z i ) = A q+1 ({z 1 , . . . , z q }, z i ), for i = q+1, . . . , n. Although slightly different VPs are applied on different steps, we will see that the validity of the outputs is not affected in practice.

SVM Venn Predictor

We define a taxonomy based on the output of the multiclass SVM. As explained in section 3, the validity of a VP is guaranteed under the i.i.d. assumption, regardless of the taxonomy used. For instance, a taxonomy that puts all examples in one single category would still give a valid predictor. Nevertheless, the performance of each VP is highly affected by the information provided from the categories defined in a taxonomy.

In this work, our taxonomy is simply based on the largest SVM score of the multiclass SVM. Therefore, each example is categorized according to the SVM classification. This taxonomy will give c categories and it is the simplest taxonomy we may define using the output of the SVM. If the SVM is good at classifying examples, then each category should contain sufficient information for the VP to perform well in terms of accuracy.

Experiments and results

In order to show the validity of the probability estimates of our method, we conduct experiments in the on-line mode. Initially all examples are test examples and they are added to the training set one by one after a prediction for each one is made. We calculate the cumulative average accuracy of the predictor, and the cumulative average probability. The cumulative average accuracy is calculated as the total accuracy of all tested examples, divided by the total number of tested examples. In the same way we calculate the cumulative average probability. If the methods provide well calibrated probability estimates, the cumulative average accuracy should be near the cumulative average probability. We test all algorithms described in this paper: Platt's method; SVM Binning; SVM with Isotonic Regression (SVM-IR); and our SVM IVP. In our experiments, we did not try to improve the accuracy of these methods, instead we focused our work on testing the validity of the probability estimates. The underlying SVM algorithm that we have used works with the RBF kernel. We test each algorithm two times, one with a RBF parameter set to an optimal value, and another with a RBF parameter set to the optimal value divided by 10 (we do this in order to test the difference in the results when the predictors do not perform so well). The optimal value for each experiment was chosen based on offline tests (10-fold cross validation) that have been conducted with a standard SVM predictor. The standard SVM predictor was tested with the RBF parameter ranges of [0.1, 1] with steps of 0.1, and [START_REF] Bohanec | Knowledge acquisition and explanation for multi-attribute decision making[END_REF][START_REF] Drish | Obtaining calibrated probability estimates from Support Vector Machines[END_REF] with steps of 1. The number of bins for the SVM Binning method was set to b = 10. In our experiments with the IVP we have set q = 0.7(n -1) . In the next two subsections, we describe our results on two multi-classification datasets.

Car evaluation dataset

The Car Evaluation dataset was derived from hierarchical decision model [START_REF] Bohanec | Knowledge acquisition and explanation for multi-attribute decision making[END_REF] and is available at [START_REF] Frank | UCI machine learning repository[END_REF]. The dataset contains 1728 examples with 6 features for each example. There are 4 classes for this dataset which describe the car acceptability based on features that describe the price, technology, and comfort of a car. In Figure 1, we show the results of the four methods on the Car Evaluation dataset. The best RBF parameter for this dataset is 0.2. For the first three methods we plot the cumulative average probability for the output classifications along with their cumulative accuracy, while for the proposed approach we plot the upper and lower cumulative probability for the output classifications along with their cumulative accuracy. One would expect the curves in each plot to be relatively near if the probabilities produced by the corresponding method were well calibrated. However this is true only for the IVP in both experiments and for Platt's method only with the optimal RBF parameter. When the RBF parameter is 0.2 the accuracy is around 90% for all methods which is the expected accuracy on this dataset. In contrast, when we set the RBF parameter to 0.02 the accuracy is reduced to around 70% (which is near the percentage of the first class), while the probability estimates are near 100% for all methods except the IVP. As shown in the last row of the graphs, the IVP probability estimates are automatically lowered to around 68%, which is near the actual accuracy.

To confirm our observations from the graphs we calculated the 2-sided pvalues of obtaining a total accuracy with the observed deviation from the expected accuracy given the probabilities produced by each method. In the case of the IVP we used the mean of the upper and lower probabilities as the probability of each prediction being correct. The p-values obtained for the outputs of the Platts's method with the RBF parameter set to 0.2 and the IVP with both parameter values were above 0.15. However, the p-values in all other cases were below 10 -50 . This shows that the probabilistic outputs produced by the three methods can be far from being well calibrated. Even for Platt's method, a wrong selection of the RBF parameter leads to misleading outputs. This does not happen with Venn Prediction which produces well calibrated outputs regardless of the underlying algorithm or the taxonomy used.

Red Wine quality dataset

The Red Wine quality dataset contains 1599 examples of physicochemical features of red variants of the "Vinho Verde" wine [START_REF] Cortez | Modeling wine preferences by data mining from physicochemical properties[END_REF]. This dataset can be used as a regression or a classification problem. Each example has a quality score from 1 to 10. In this work, we have used the scores as 10 different classes from 1 to 10. This dataset is particularly difficult and requires some pre-proccessing to remove redundant features, or even reduce the number of classes. For instance, some classes have very few or even no examples in the training set. In our experiments, we have intentionally left the dataset to its original state in order to demonstrate the reliability of our probability estimates on difficult problems where the underlying algorithm may not be able to fit the data very well. In Figure 2, we show the online results of the four methods on the Wine quality dataset. The best RBF parameter on this dataset is 0.6. From the results, we can see that Platt's method, SVM-Binning, and SVM-IR did not give reliable probability estimates (due to the difficulty of the task), whereas the IVP has automatically lowered the probability estimates and has given well calibrated results in both cases. The 2-sided p-values for the IVP were above 0.3, whereas for all other methods with both RBF parameters the p-values were below 10 -50 .

Conclusion

In this work, we have examined existing methods that convert SVM scores into probability estimates. We have shown that there is no trust in the probability estimates produced by these methods, especially when the underlying algorithm is not well configured or when the dataset is difficult. For the purpose of overcoming this limitation, we have developed an Inductive Venn Predictor based on SVMs, which guarantees (under the i.i.d. assumption) that the probability estimates will be well calibrated, regardless of the configuration of the algorithm or the difficulty of the task. Additionally, the proposed IVP overcomes the computational inefficiency problem which renders the original Venn Prediction framework unsuitable in the case of large datasets. Our future aim is to improve the performance of our IVP in terms of accuracy by introducing better taxonomies that can be derived from the unthresholded scores of SVMs. Furthermore, we wish to investigate whether the one-against-all procedure is one of the causes for the non calibrated probability estimates of the existing methods, and we would like to compare our IVP with other multiclass procedures.

  on Car Evaluation dataset with RBF parameter 0.02. on Car Evaluation dataset with RBF parameter 0.2. low bound upper bound accuracy

Fig. 1 .

 1 Fig. 1. Online experiments of all four methods on the Car evaluation dataset, RBF parameter is 0.02 on the left column and 0.2 on the right column.

  method on Red Wine dataset with RBF parameter 0.06. on Wine quality dataset with RBF parameter 0.06. on Wine Quality dataset with RBF parameter 0.6. low bound upper bound accuracy

Fig. 2 .

 2 Fig. 2. Online experiments of all four methods on the Wine quality dataset, RBF parameter is 0.06 on the left column and 0.6 on the right column.

The training set is in fact a multiset, as it can contain some examples more than once.

Acknowledgments. This work was supported by the Cyprus Research Promotion Foundation through research contract TPE/ORIZO/0609(BIE)/24 ("Development of New Venn Prediction Methods for Osteoporosis Risk Assessment").