James Fishbaugh

Marcel Prastawa

Guido Gerig

Stanley Durrleman

GEODESIC REGRESSION OF IMAGE AND SHAPE DATA FOR IMPROVED MODELING OF 4D TRAJECTORIES

teaching and research institutions in France or abroad, or from public or private research centers.

Geodesic regression of image and shape data for improved modeling of 4D trajectories

INTRODUCTION

Analysis of longitudinal data incorporating both spatial and temporal information is essential for various clinical tasks such as predicting patient outcome and measuring efficacy of different therapeutic strategies. A crucial tool for longitudinal analysis is regression of observed data, which enables interpolation to generate continuous evolution models as well as extrapolation to predict future observations. Regression models are also necessary for conducting population studies comparing the change trajectories of different subjects.

In medical imaging, it is important to consider image data in anatomical context, which motivates regression on image and shape data in different combinations (a multi-object complex). A variety of regression schemes have been proposed on images or shapes, although available methods do not handle them jointly. For example, the extension of kernel regression for image data [START_REF] Davis | Population shape regression from random design data[END_REF] or piecewise linear regression for time series of images [START_REF] Khan | Representation of timevarying shapes in the large deformation diffeomorphic framework[END_REF] and shapes [START_REF] Durrleman | Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape data[END_REF]. Combining intensity and geometric information has been explored for registration [START_REF] Cachier | Iconic feature based nonrigid registration: the pasha algorithm[END_REF].

To conduct statistical analysis on 4D (3D space + time) data, it is particularly useful to consider compact generative regression models which have low number of parameters situated at only one chosen time point. Geodesic regression is such a model and is fully characterized by baseline images and shapes (the intercept) and the tangent vector defining the geodesic at the baseline objects (the slope). Geodesic regression frameworks for images [START_REF] Niethammer | Geodesic regression for image time-series[END_REF][START_REF] Singh | A vector momenta formulation of diffeomorphisms for improved geodesic regression and atlas construction[END_REF] and for shapes [START_REF] Fishbaugh | Geodesic shape regression in the framework of currents[END_REF] have been proposed using the LDDMM setting. However, no notion of how to combine images and shape data is provided.

We propose a novel geodesic regression framework that leverages image and shape data together to estimate a single deformation of the ambient space. We use the currents representation for geometric data that allows flexible representation of a wide variety of shape objects such as point sets, curves, or surface meshes, without the need for point correspondence between shapes. Compared to image regression alone, shape data provides anatomical information that constrains the regression, especially in cases where images have low contrast, by placing larger weights on regions with anatomical importance. Compared to shape regression alone, image information provides data in areas where segmentations are not available, as well as providing context to regions surrounding anatomical objects. Our framework uses the control-point parameterization of geodesic flows introduced in [START_REF] Durrleman | Optimal data-driven sparse parameterization of diffeomorphisms for population analysis[END_REF], which makes the parameterization of the deformation independent from the data. This allows us to keep a reasonable dimension of the parameterization, which is determined by the complexity of anatomical changes in time, and not the sampling of the data. We therefore combine image and shape data without introducing a complexity overhead.

METHODOLOGY

We perform regression on observed time-series data of images I i and shapes X i , each acquired at time-point t i . Shape data may consist of a mix of point sets, curves, or surface meshes where all vertices are concatenated into one vector X i .

We use the control point formulation of [START_REF] Durrleman | Optimal data-driven sparse parameterization of diffeomorphisms for population analysis[END_REF] to generate geodesic flows of diffeomorphisms. Let S 0 = {c 0,k , α 0,k } be a set of momentum vectors α 0,k attached to control points c 0,k distributed in the image domain. Geodesic flows are computed by evolving control points and momenta by integrating the following Hamiltonian equations over the time interval of interest:

             ċk (t) = N i=1 K(c k (t), c i (t))α i (t) αk (t) = - N i=1 α k (t) T α i (t)∇ 1 K(c k (t), c i (t)) (1)
with initial conditions c k (0) = c 0,k and α k (0) = α 0,k (assuming starting time-point to be 0), and K is a Gaussian kernel with variance σ 2

V which controls the spatial scale of deformation. For simplicity, we write these equations as Ṡ(t) = F (S(t)) with S(0) = S 0 . The convolution of the momenta defines the following time-varying velocity field: v(t, x) = N i=1 K(x, c i (t))α i (t) for any point x in the domain. The velocity is used to deform the domain: a particle at point x at time 0 moves to φ(t, x) at later time t, where φ(t, x) follows the integral curve of ∂φ(t,x) ∂t = v(t, φ(t, x)) starting with φ(0, x) = x. In this formulation, the velocity of the particle is given by the field v(t, .) at its current location. It has been shown in [START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF] that for all t, φ(t, .) is a 3D diffeomorphism.

Following this model, the vertices of a given baseline shape complex concatenated into a vector X 0 move at time t to X(t) = φ(t, X 0 ), which satisfies the ordinary differential equation (ODE): Ẋ(t) = v(t, X(t)) with X(0) = X 0 . To make explicit the dependency of the equation of motion on S(t), we write it as: Ẋ(t) = G(X(t), S(t)).

A given baseline image I 0 is also deformed by the flow Fig. 1. Conceptual overview of geodesic regression on multiobject complexes containing both image and shape data. The framework estimates parameters at t = 0 which consist of the baseline image I 0 and shape X 0 along with the deformation model parameterized by control points c 0 and initial momenta α 0 such that overall distance between the deformed objects and the observations are minimal.

of diffeomorphisms and its trajectory is given as

I(t) = I 0 • φ(t, .) -1 .
The inverse flow satisfies the equation ∂φ(t,.) -1 ∂t = -dφ(t, .) -1 v(t, .). For the sake of simplicity, we denote Y (t, .) = φ(t, .) -1 , a L 2 function that maps the point x to its position at time t under the inverse flow φ -1 (t, x). This maps satisfies Ẏ (t, .) = -dY (t, .)v(t, .) = H(Y (t, .), S(t)), where we make explicit the dependency on S(t). At time t, the intensity of the warped baseline image at voxel position x is given by I(t, x) = I 0 (Y (t, x)) using 3D interpolation.

A conceptual overview of our framework is shown in Fig. 1 where regression is performed by minimizing the overall distance between the observations and the deformed baseline objects (shapes and/or images). Let d(X(t i ), X i ) be a metric between the deformed baseline shape complex X 0 at time t i and the data shape complex X i . This metric may be a weighted sum over each component of the shape complex of the currents metric between sets of curves or surface meshes. This term essentially depends on X(t i ) and is denoted A(X(t i )). Similarly, we have a metric d(I(t i ), I i ) denoted as B(Y (t i , .)) that is the sum of squared differences between the deformed baseline image I 0 • Y (t i , .) and the observed image I i .

The geodesic regression problem amounts to finding the deformation parameters S 0 and baseline anatomical configuration (I 0 , X 0 ) such that the following criterion is minimized:

E(S 0 , I 0 , X 0 ) = ti λ It i A(X(t i )) + λ St i B(Y (t i , .)) + L(S 0 ) (2) subject to      Ṡ(t) = F (S(t)) S(0) = S 0 Ẋ(t) = G(X(t), S(t)) X(0) = X 0 Ẏ (t, .) = H(Y (t, .), S(t)) Y (0, .) = Id (3)
where the regularizer L(S 0 ) = N i,j=1 α T 0,i K(c 0,i , c 0,j )α 0,j is the squared norm of initial velocity and weights on image and shape matching λ It i and λ St i .

As shown in the supplemental material (www.cs.utah. edu/ ˜jfishbau/docs/isbi2014_eqns.pdf), the gradient is computed by integrating 3 linear ODEs with source terms from final time-point T f back to time-point 0:

∇ S0 E = ξ(0) + ∇ S0 L ∇ X0 E = η(0) 1 2 ∇ I0 E = tj Splat Y (tj ,.) (I 0 • Y (t j , .) -I i ) with              η(t) = -∂ 1 G(t) T η(t) - ti ∇ X(ti) A δ(t -t i ) θ(t) = -∂ 1 H(t) † θ(t) - ti ∇ Y (ti,.) B δ(t -t i ) ξ(t) = -∂ 2 G(t) T η(t) -∂ 2 H(t) † θ(t) -d S(t) F T ξ(t) with final conditions η(T f ) = θ(T f ) = ξ(T f ) = 0.
The vector η is same size as X 0 , which brings back to time t = 0 the gradients of the data matching terms, and is used to update the position of the vertices of the baseline shape complex. Similarly θ is of the same size as Y (0, .) (an image of vectors in practice) which integrates the successive gradients of the image matching terms that acts as jumps in the differential equation. Finally, ξ is a variable of the same size as S 0 which is used at time t = 0 to update the deformation parameters (the position of the control points and their momentum vectors). The gradient with respect to the baseline image involves the splatting of the current residual images at positions Y (t i , .) as done in [START_REF] Durrleman | Optimal data-driven sparse parameterization of diffeomorphisms for population analysis[END_REF].

RESULTS AND DISCUSSION

Pediatric Brain Development: We explore the impact of joint image and shape regression in modeling pediatric brain development. The data consists of T1W images of the same healthy child observed at 6, 12, and 25 months of age. Regression on images alone is difficult in this case due to the very low contrast in the 6 month old image. Despite the low contrast, tissue segmentations can still be reliably and consistently estimated [START_REF] Vardhan | Quantifying regional growth patterns through longitudinal analysis of distances between multimodal MR intensity distributions[END_REF]. We estimate a geodesic model using only T1W images and a model jointly on images and white matter surfaces to emphasize the development of the tissue interface. We initialize 120 control points on a regular grid with the deformation kernel σ V = 20 mm. Finally, due to limited contrast at 6 months, we estimate the baseline at 25 months and follow the evolution backwards in time.

The results of geodesic regression are shown for several snapshots in time in Fig. 2. The model estimated using only images mostly captures the scale change, but does not cap-Fig. 2. Images and deformations estimated by geodesic regression using images alone (top) and jointly on images and white matter surfaces (bottom). Regression jointly on image and shape results in a more realistic evolution which captures detailed changes in brain tissue in addition to the increase in brain size. In both cases, geodesic regression was estimated backwards in time. Volume is measured continuously from the modeled shape trajectories, not fitted to discrete volume measurements. The model estimated on images alone fails to capture the volume loss. Evolution of caudates for the image only model is not estimated, but instead we shoot the baseline caudate shapes along the geodesic estimated from images alone. Note: measurements extracted continuously from non-linearly deforming shapes can produce either linear or non-linear trends with no prior assumption of linearity. ture much deformation in the interior of the brain. The model estimated jointly on image and shape captures more detailed development as white matter stretches and expands. Neurodegeneration in Huntington's Disease: Next, we investigate the application of joint image and shape regression to Huntington's disease (HD) where accurate 4D models are needed to measure the effectiveness of therapies or drug treatments. In HD, degeneration of the caudate has been shown to be significant [START_REF] Eh | Rate of caudate atrophy in presymptomatic and symptomatic stages of huntington's disease[END_REF]. Here we explore T1W image data from a single patient diagnosed with HD scanned at 58, 59, and 60 years of age. Sub-cortical structures are segmented, manually verified, and cleaned. Models are estimated using only T1W images as well as T1W images plus caudate surfaces. Control points are initialized on a regular grid with 10 mm spacing with kernel σ V = 10 mm.

The trajectory of caudate volume extracted after regression is shown in Fig. 3. The model estimated from images alone fails to capture the volume loss observed in both caudates, and rather, shows an increase in right caudate volume. By incorporating caudate shape data in model estimation, we Conclusions: We presented a novel geodesic regression framework that jointly considers image and shape information in the LDDMM framework, where dense diffeomorphisms are built using a control point formulation. This formulation decouples deformation parameters from input object parameters (e.g., voxels, surface points) providing greater flexibility and consistency in mapping different object types across time. Our regression model seamlessly handles images and multiobject complexes consisting of curves, and/or surfaces in different combinations. Experiments show that our framework effectively combines image and shape information to estimate a single deformation of the ambient space, resulting in improved modeling of 4D trajectories.

Fig. 3 .

 3 Fig. 3. Caudate volume extracted continuously after regression compared to observed caudate volumes (circles and x's).Volume is measured continuously from the modeled shape trajectories, not fitted to discrete volume measurements. The model estimated on images alone fails to capture the volume loss. Evolution of caudates for the image only model is not estimated, but instead we shoot the baseline caudate shapes along the geodesic estimated from images alone. Note: measurements extracted continuously from non-linearly deforming shapes can produce either linear or non-linear trends with no prior assumption of linearity.

Fig. 4 .

 4 Fig.4. Top) Evolution estimated on images alone. Evolution of caudates are not estimated, but instead we shoot the baseline caudate shapes along the estimated geodesic. Bottom) Evolution estimated jointly using images and caudate shapes. Regression on images alone results in a slight expansion of ventricles, but does not capture the shrinking of caudates. Our method is able to capture both the expansion of ventricles and the shrinking of caudates.

* Supported by grants: RO1 HD055741 (ACE, project IBIS), U54 EB005149 (NA-MIC), and U01 NS082086 (HD) and the Utah Science Technology and Research (USTAR) initiative at the University of Utah.


Figure 1: Sequential acquisition of genetic alterations in tumor evolution.

Figure 2: Measuring alterations in a cancer genome. A. Whole genome karyotyping image (SKY) of breast cancer cell line MDA-MB-468 (http://www.path.cam.ac.uk/pawefish/); B. High resolution arrays (BAC, CGH) representing relative copy number variation profile; C. High resolution Affymetrix SNP chip 6.0 array profiles of a primary breast tumor: allelic imbalance (AI), copy number variation (CNV) and recognized absolute copy number (CN) profiles; D. Circus plot showing chromosomal translocations and copy number variation obtained by NGS for a primary breast tumor (figure from (Natrajan et al., 2012)).

Introduction

Abnormal genetic content is observed in many tumors and is considered as one of the hallmarks of cancer [START_REF] Hanahan | Hallmarks of cancer: the next generation[END_REF]. Evolution of cancer is thought to be tightly connected with evolution of cancer genome [START_REF] Podlaha | Evolution of the cancer genome[END_REF]. Studying the cancer genome is important as it may unravel the key genomic events or some particular genomic features, which could shed light on tumor biology and have clinical implications. When considering the genome of an advanced tumor, we observe the "end" point of its evolution often showing numerous acquired rearrangements (Figures 1 &2). Exact genetic pathways and chronology of events acquisition remain unclear and have only started to be unraveled using whole genome sequencing [START_REF] Stratton | The cancer genome[END_REF]. Reconstructed history of several tumor genomes supported the hypothesis of sequential acquisition of genomic events (with possibly variable density in time) and clonal development [START_REF] Nik-Zainal | The life history of 21 breast cancers[END_REF]. Although there is still debate whether genomic alterations are causes or consequences of cancer, numerous genomic features have already been linked to initiation and progression in many types of cancer [START_REF] Tran | Cancer genomics: technology, discovery, and translation[END_REF]. Some of these genomic features were introduced into clinical practice largely contributing to the diagnostics and treatment choices [START_REF] Stuart | Linking somatic genetic alterations in cancer to therapeutics[END_REF].

In the last decades, several high throughput techniques have been developed to measure genetic alterations in cancer. Spectral Karyotyping (SKY) or other mitotic chromosome imaging gives a general view on chromosomal content in a cancer cell, but the resolution is low, and fine description of genetic alterations based on these images is largely impossible (Figure 2A). Array comparative genomic hybridization (CGH) has been widely in use for the last 20 years showing relative pattern of genetic alterations at rather high resolution (Figure 2B). Single nucleotide polymorphism-based platforms (SNP-arrays) have progressively replaced CGH in tumor characterization, as they allowed estimation of absolute copy numbers and allelic contents (Figure 2C). Next Generation Sequencing (NGS) gives the most complete information about cancer genome, including point mutations and chromosomal translocations at base pair resolution (Figure 2D).

In this chapter we consider basic hypothesis, problem statements and technological and computational solutions for analysis of copy number alterations in tumor genomes. We provide a data mining technique (based on the GAP method described in [START_REF] Popova | Genome Alteration Print (GAP): a tool to visualize and mine complex cancer genomic profiles obtained by SNP arrays[END_REF]) which allows extraction of absolute copy numbers and allelic contents from the whole genome copy number variation and allelic imbalance profiles obtained by SNP arrays or NGS.

Basic Hypothesis in Exploratory Analysis of Cancer Genome

What events are essential for tumor initiation and evolution and how they are manifested in the tumor genome? In the simplistic view, there are two major contributors to cancer development: first, tumor promoting events, i.e. "switching on" oncogenes (growth factors, etc); second, elimination of antitumor protection, i.e. "switching off" tumor suppressor genes (genome guardians, checkpoints, etc) [START_REF] Grander | How do mutated oncogenes and tumor suppressor genes cause cancer?[END_REF]. Major genomic mechanisms of oncogene activation are gain-of-function mutations, chromosome translocations (which, for example, place the gene under an activating promoter) and amplification of chromosome region containing the gene. Major genomic mechanisms resulting in inactivation of tumor suppressor gene are mutations and deletions of chromosome region containing the gene. The Knudson's two hit hypothesis for recessive tumor suppressor gene implies tumor initiation when two copies of the gene are inactivated (i.e. in the normal diploid human genome, where autosomic genes exist in two copies) [START_REF] Knudson | Mutation and cancer: statistical study of retinoblastoma[END_REF]. Double deletion of a locus or combination of mutation and loss of normal copy (resulting in so called loss of heterozygosity [LOH]) are the genomic indicators of possible tumor suppressor genes. The association between recurrent deletions and tumor suppressors as well as amplifications and oncogenes has been confirmed by the analysis of more than 1300 tumor genomes [START_REF] Beroukhim | The landscape of somatic copy-number alteration across human cancers[END_REF].

Besides local genomic indications for lost and acquired genes the whole pattern of alterations in a tumor genome, such as chromosomal content (ploidy), number of chromosomal breaks, complexity of intraand inter-chromosomal rearrangements, may evidence disruption of some pathways, such as DNA repair or chromosome maintenance. For example, BRCA1 -/-breast tumors display more rearranged genomes compared to non-BRCA1 tumors; alternatively, non-BRCA1 tumors with highly rearranged genomes could be mutated in another gene from the DNA repair pathway, such as BRCA2 [START_REF] Popova | Ploidy and large-scale genomic instability consistently identify basal-like breast carcinomas with BRCA1/2 inactivation[END_REF]. In this regard analysis of tumor genome profiles could contribute to the tumor classification problem.

A number of studies have shown that although each tumor genome is unique, some alterations observed within the type of tumor/tissue are highly recurrent. For example, more than 90% of triple negative breast tumors carry a mutation in the TP53 gene and display LOH on the chromosome 17 [START_REF] Manie | High frequency of TP53 mutation in BRCA1 and sporadic basal-like carcinomas but not in BRCA1 luminal breast tumors[END_REF]; almost all clear cell renal cell carcinomas have deletion in the 3p chromosome arm [START_REF] Hagenkord | Clinical genomics of renal epithelial tumors[END_REF]; majority of low grade estrogen receptor positive ductal breast carcinomas have 16q deletion and half of them display 16p gain [START_REF] Natrajan | Loss of 16q in high grade breast cancer is associated with estrogen receptor status: Evidence for progression in tumors with a luminal phenotype?[END_REF]; numerous subtypes of cancers have amplification of 8q; etc. Some of these recurrent alterations were found to target certain tumor related genes as the minimal shared region of alteration clearly identifies a gene locus. Other recurrent alterations comprise the whole chromosome arm. In this case the list of possible functionally related targets is large and uncertain. On the other hand, it was shown that some genomic regions "prefer" to be gained or lost independent of the type of tumor [START_REF] Beroukhim | The landscape of somatic copy-number alteration across human cancers[END_REF]. The current understanding associates some of recurrent alterations with chromosome structure or function (for example, fragile sites, long repeats favoring chromosome recombination, formation of di-centric chromosomes leading to breaks and fusion cycles, interference of transcription and replication in large genes, etc.) [START_REF] Bignell | Signatures of mutation and selection in the cancer genome[END_REF][START_REF] Cassidy | Genome instability mechanisms and the structure of cancer genomes[END_REF]. Genomic alterations "accompanying" causative (driver) mutations or arising as a consequence of functioning due to extensive growth were referred as passenger events [START_REF] Stratton | The cancer genome[END_REF]. Thus, the main challenge in tumor genome interpretation is distinguishing alterations associated with the driver and the passenger events.

To conclude, for unraveling cancer genome complexity we need to annotate alterations in tumor genome in order to find recurrent amplifications, deletions, gains, losses, and LOH, which possibly target cancer related genes. We need to be able to characterize tumor ploidy and the level of rearrangements to describe cancer genomes as completely as possible.

Measuring Genetic Alterations by SNP Based Technique

The normal human genome contains two copies of each autosome (chromosomes 1-22) and two copies of the sex chromosome X in females (XX) or one copy of each of the sex chromosomes X and Y (XY) in males. These two alleles have paternal and maternal origins and are 99.9% identical in DNA sequence. A difference between two alleles is referred as a genetic variation and includes, in particular, Single Nucleotide Polymorphisms (SNPs) (Figure 3). We say that a SNP is homozygous if paternal and maternal alleles have identical nucleotides; the SNP is heterozygous if paternal and maternal alleles have different nucleotides. Status of SNP in individual genome is defined as homozygous (two identical nucleotides) or heterozygous (distinct nucleotides). Single nucleotide variation is designated as a SNP, if its minor allele frequency exceeds 1% in at least one population; otherwise it is called variant or mutation.

Homozygous SNPs are called non-informative because they do not allow the two alleles to be distinguished. At the moment, near 30 million of SNPs are annotated in the dbSNP database (to be annotated as a SNP, variant frequency at the specific locus should exceed 1% in, at least, one population [START_REF] Sherry | dbSNP: the NCBI database of genetic variation[END_REF]).

The cancer genome is characterized by abnormal genomic content, meaning, the number of chromosomes and their structure are different from the normal state. Variation in copy numbers is a result of losses of one of the two alleles (or both alleles, so called homozygous deletion, rarely comprising a large genomic region), gains of up to several copies of one or both alleles, and copy neutral LOH (two identical alleles, also called uniparental disomy) [START_REF] Albertson | Chromosome aberrations in solid tumors[END_REF]. Which allele is lost or gained in a cancer genome is essentially unknown unless the parents of the patient are also genotyped (that is not a common practice in tumor genomic studies). That is why alleles are designated arbitrarily as A and B, and characterization of allelic content in tumor genome consists of indicating copy number (CN) and major allele (MA) counts (for example, CN = 3 and MA = 2 correspond to genotypes AAB and ABB, and could be described as "two identical alleles out of three"). Allelic content is characterized by B Allele Frequency (BAF) or Allelic Difference (AD), which represent the ratio of the major allele counts and the copy number or the difference between major and minor allele (MI) counts (MI = CN -MA), respectively:

BAF = MA/CN = n c B /CN, (1)
AD = MA-MI = n c B -n c A , ( 2
)
where CN is copy number, n c B ,n c A are numbers of B and A alleles in cancer genome (historically, B allele counts are associated with the major allele counts). Plotting genomic states in two dimensions, namely, allelic content (x axis) and copy number (y axis) gives visual representation of all possible allelic contents and copy numbers, which could be observed in a tumor genome (Figure 4). Monoallelic genomic states (B, BB, BBB, etc., designated by red color in the table and red rectangles in Figure 4) represent genomic states with LOH.

Measuring genetic alterations in a cancer genome with SNP-based technology consists of measuring allelic contents in numerous SNP loci along the genome. If the technology provided perfect measurement of the number of A and B alleles in each SNP locus, each genomic region would be characterized by the pair of values (CN, MA) corresponding to one of the genomic states shown in Figure 4. However, there is no such perfect technology available yet. Whole genome profiles of genetic alterations provided by SNParrays (or NGS) represent relative copy number variation (CNV) and allelic imbalance (AI) profiles (Figures 2C &5, see Appendix A for SNP array platforms description). Allelic imbalance profile characterizes allelic content by BAF or AD, depending on the normalization applied to the measured profile. These profiles are affected by technological noise and experimental variation (biological sample preparation); furthermore, the measured tumor sample often has significant admixture of the normal genome (so called normal contamination, as the single cell technology is not yet a common place, the measured tumor sample is usually a mixture of tumor and normal stromal cells). Thus, annotation of measured relative variation profiles with absolute copy numbers and major allele counts represents a problem statement for further data mining.

Genome Alteration Print (GAP) for Mining Genetic Alterations

Here we describe a data mining technique developed for SNP arrays, which allows extraction of absolute copy numbers and allelic contents from the whole genome CNV and AI profiles. The method is based on the structure denoted by Genome Alteration Print (GAP) [START_REF] Popova | Genome Alteration Print (GAP): a tool to visualize and mine complex cancer genomic profiles obtained by SNP arrays[END_REF]. GAP is a two dimensional representation of segmented CNV and AI profiles of a measured tumor sample and characterizes the spectrum of rearrangements presented in the tumor (Figure 5, right panels). In order to build the GAP of a tumor, CNV and AI profiles are segmented; the values are smoothed within the segments by median (CNV) and mode (AI); the segments are plotted on the AI × CNV plane as the circles of the radius proportional to the segment size (see Appendix B for details).

GAP patterns of the breast tumor sample measured on the Illumina and Affymetrix platforms (Figure 5, right panels) resemble the patterns of copy numbers and allelic contents represented by BAF and AD respectively (Figure 4). Tumor GAPs differ from the copy number / allelic content plots by (1) uneven distribution of CN levels along the vertical axis (due to the original CNV profiles represented in the log-scale: Log R Ratio and Log2Ratio in Illumina and Affymetrix platforms, respectively); and (2) a shift on the horizontal axis toward 0.5 or 0 of BAF or AD, respectively (due to the normal contamination present in the tumor).

In order to account for normal contamination in copy number / allelic content template we considered the sample containing proportion p of normal cells and (1-p) of tumor cells. Adding proportion p of normal heterozygous (AB) signal to the copy numbers and allelic contents we easily can obtain the model patterns accounting for the normal contamination (Figure 6):

CN p = (1 -p) ⋅CN + 2p, (3)
BAF p = (1 -p) ⋅ n c B + pn n B (1 -p) ⋅CN + 2p , ( 4
)
AD p = (n c B -n c A ) ⋅ (1 -p) + (n n B -n n A ) ⋅ p, ( 5
)
where

c A = CN -n c B ; n n B = 1, if n c B < CN; n n B = 1,2, if n c B = CN).
These two patterns represent the model templates for mining genetic alteration measured by SNP arrays. It is worth noting that LOH are now represented by two sets of genomic states (designated by red quadrangles in the Figure 6) and the distance between the two LOH sets characterizes the level of normal contamination. Because of normal contamination SNPs with acquired homozygosity had addition of normal heterozygous signal, which shifted their AI towards heterozygous states (for example, (1 -p)⋅ BB + p⋅AB). SNPs homozygous in the germline did not change their position ((1 -p)⋅ BB +p⋅BB = BB) (compare Figures 4 &6). Thus, moderate level of normal contamination could help distinguish segments with germline and acquired homozygosity. However, increasing the normal contamination proportion leads to the pattern shrinkage towards 2 copies, in extremity (when p is close to 1) converging to a normal genome. Acquired and germline homozygosities are not distinguishable in the case of pure tumor sample or cell line.

All the GAPs obtained for a large set of SNP-arrays measuring breast, prostate, lung and some other tumors were found resembling the model patterns of copy numbers and allelic contents defined above. However, the dynamics of change in measured CNV profiles was observed to have high degree of variation from sample to sample. To account for this experimental variation we introduced the coefficient of CNV contraction q, which is supposed to adjust the difference in scales between model CN and measured CNV profiles (Figure 6).

The common way of annotating a SNP array profile consisted thus in finding the model template which fits best to the measured GAP. Finding the model template means defining three parameters, namely, proportion of normal contamination (p), coefficient of CNV contraction (q), and position of 2 copies on the CNV scale (C), such that the model template maximally corresponds to the measured GAP. GAP plane is annotated with CNs and MAs from the model template and genomic segments are annotated accordingly (Figure 7).

A quality of fitness of the measured GAP to the model template is assessed by the genome coverage in terms of number of SNPs that are explained by the model. However, fitting regular structures might not result in one unique best solution. To avoid ambiguity, we filter redundant templates and choose interpretation with the lowest possible CN set and minimal possible normal contamination. All parameters of the model template were set up and extensively tested on a large cohort of tumor genomes and under manual control of recognition quality (visual representation of tumor GAPs and model templates, as shown in Figure 7, allows judging on correct or incorrect recognition). Current implementation of the GAP method provides near 80% rate of correctly found model templates. Moreover, complicated cases could be annotated manually by introducing three parameters, which (according to a user) provide a proper correspondence between tumor GAP and the model template. More details and algorithm of recognition are presented in Appendix C; R scripts and full details of the application are available at (http: //bioinfo-out.curie.fr/projects/snp_gap/).

Validation of the GAP Method

Our group used the GAP method to process numerous SNP array profiles mainly of breast tumors. As a result of the recognition procedure we obtained the level of normal contamination (p), and the profile of segmental copy numbers and allelic contents (segmental genotypes) for each tumor (one example is shown in Figure 2C). With the current implementation recognition of absolute copy number ranged from 0 to 8 copies (all segments with the copy number variation exceeding 8-copy level were ascribed 8-copy status). Thus, 22 possible segmental genotypes were discriminated (copy number / major allele count): B (1/1); BB (2/2) and AB (2/1); BBB (3/3) and ABB (3/2); BBBB (4/4), ABBB (4/3) and AABB (4/2); etc.

The way we model the proportion of normal contamination (p) for BAF evaluation has been already described in a number of publications and confirmed by the dilution series and computer simulations [START_REF] Nancarrow | SiDCoN: a tool to aid scoring of DNA copy number changes in SNP chip data[END_REF]Staaf et al., 2008a). To verify our evaluation of the normal contamination and AD we considered the dilution series of the lung cancer cell line (H-1395) available on Affymetrix SNP 6.0 plat- form (GEO GSE29172) [START_REF] Rasmussen | Allele-specific copy number analysis of tumor samples with aneuploidy and tumor heterogeneity[END_REF]. For the proportions 0, 0.3, 0.5, and 0.7 of the normal cells admixture, we obtained the model parameter 2p to be 0.1, 0.44, 0.6, and 0.76, respectively (see Figure 6 B). As far as we can conclude from this small series, 2p could be considered as a reasonable upper estimation of the normal contamination in Affymetrix SNP 6.0. It is worth noting that (i) the model parameter p significantly underestimated normal contamination; (ii) 70% contaminated tumor displayed almost flat merely recognizable profile.

In order to validate our SNP array interpretation in the total set of copy numbers, we compared DNA contents inferred from recognized copy number profile with those provided by flow cytometry (FCM). DNA index provided by FCM characterizes DNA content of tumor genome relative to normal diploid genome, for which DNA index is set to 1. DNA indexes from SNP arrays were estimated by averaging segmental copy numbers divided by 2. DNA indexes inferred from SNP-arrays are in close correspondence with actual tumor DNA indexes measured by FCM analysis for 20 out of 23 tested samples of breast carcinomas (Figure 8A).

We also validated our recognition procedure by comparing inferred number of chromosomes with known karyotypes or SKY data of 25 breast cancer cell lines (http://www.lgcstandards-atcc.org/, http://www.path.cam.ac.uk/~pawefish/). The number of chromosomes was estimated by the sum of the copy numbers detected at the pericentric regions. The status of the pericentric region of each chromosome arm was defined by the corresponding juxta-centromeric segment when the latter contained 500 SNPs or more (Affymetrix SNP 6.0). When not measurable, missing values were substituted by the modal copy number of the considered chromosome arm (3.4±2.2 out of 41 chromosome arms per genome were substituted in the series).

Error rate was less than 2 chromosomes per sample (1.7±2.3) in the considered set of breast cancer cell lines (Table 1). To conclude, copy number recognition by the GAP method results in correct predictions in the majority of tested cases.

Other Approaches to Mine SNP Arrays and General Problems

Numerous studies have contributed to better understanding cancer genomics and in particular mining SNP array measurements [START_REF] Assie | SNP arrays in heterogeneous tissue: highly accurate collection of both germline and somatic genetic information from unpaired single tumor samples[END_REF][START_REF] Attiyeh | Genomic copy number determination in cancer cells from single nucleotide polymorphism microarrays based on quantitative genotyping corrected for aneuploidy[END_REF][START_REF] Gardina | Ploidy status and copy number aberrations in primary glioblastomas defined by integrated analysis of allelic ratios, signal ratios and loss of heterozygosity using 500K SNP Mapping Arrays[END_REF][START_REF] Lamy | A Hidden Markov Model to estimate population mixture and allelic copy-numbers in cancers using Affymetrix SNP arrays[END_REF]. However, it is the copy number / allelic content structure (Figure 4 & 6) standing behind the SNP-based genomic profiling (with linear or non-linear transformation depending on the physical properties of a measuring device and a way of data normalization) that should govern any copy number recognition procedure. Simplicity of copy number / allelic content pattern implies various possible approaches for copy number recognition. Existing integrated approaches includes Hidden Markov Models (PICNIC [START_REF] Greenman | PICNIC: an algorithm to predict absolute allelic copy number variation with microarray cancer data[END_REF], GPHMM [START_REF] Li | GPHMM: an integrated hidden Markov model for identification of copy number alteration and loss of heterozygosity in complex tumor samples using whole genome SNP arrays[END_REF]) and pattern recognition strategies (GAP [START_REF] Popova | Genome Alteration Print (GAP): a tool to visualize and mine complex cancer genomic profiles obtained by SNP arrays[END_REF], ASCAT (Van Loo et al., 2010), TAPS [START_REF] Rasmussen | Allele-specific copy number analysis of tumor samples with aneuploidy and tumor heterogeneity[END_REF]), Absolut [START_REF] Carter | Absolute quantification of somatic DNA alterations in human cancer[END_REF].

Quality of a SNP array mining procedure should be estimated by correspondence between reported tumor characteristics (such as CN and MA profiles, number of breakpoints, percent of tumor cells) and the actual tumor states in a large set of samples and under various conditions. Problems that could affect recognition potential of all the methods mentioned above could be subdivided into three groups: (i) technical problems associated with quality of measured SNP array profile and physical properties of a measuring device; (ii) tumor heterogeneity; (iii) ambiguity in interpretation of SNP array profiles.

The first class of problems associated with quality of measured SNP array profile includes various effects of unknown origin resulted in non-specific variation in CNV profile, such as high percent of outlying values; low-frequency fluctuations, so called waves in CNV profile; high-frequency saw-tooth waves, not associated with CN change; etc. AI profiles appeared to be more stable and not affected by these variations. Thus, data mining techniques more relying on AI profiles could be more successful and give more reliable estimation of tumor rearrangements in these cases. However, low quality profiles are more exceptional than regular cases in the latest generation of SNP arrays.

The problems associated with physical properties of measuring devices also include the signal saturation effects, which hamper accurate detection of high copy number levels.

The second class of problems is intra-tumor heterogeneity, which include normal contamination and subclonal structure [START_REF] Shipitsin | Molecular definition of breast tumor heterogeneity[END_REF]. Although presumably arisen from a single cell (monoclonal proliferation), cancer progression leads to the sub-populations bearing different genomic alterations (subclones) coexisting in most tumor samples [START_REF] Navin | Inferring tumor progression from genomic heterogeneity[END_REF]. Variation observed by SNP-array in a tumor genome reflects genomic alterations shared by all tumor cells and subclonal events shared by only subpopulation of tumor. CN and MA status of an alteration specific to sub-clones is generally indefinable from SNP arrays as the measured signal reflects the sum of unknown subclonal signals in unknown proportions. Intermediate CNV levels due to subclonal events could be interpreted as actual CN level and overall pattern of copy number set could be thus overestimated.

Another type of problems in mining SNP arrays is associated with classes of tumor genomes not distinguishable by SNP array profiles. This follows from the relative nature of measured profile and regular structure of allelic content / copy number pattern. For example, there is no way to distinguish the genomes displaying AB and B states versus AABB and BB ones, if there are no other events (such as ABB) in the latter case. However, this is probably a rare situation (as we can conclude at least for breast cancer genomes).

Table 1: Number of chromosomes in the cell lines according to GAP estimation, ATCC description and SKY images. Error was calculated as minimal absolute difference between GAP estimation and ATCC and/or SKY numbers.

Obviously, all aforementioned methods performed well on the high quality SNP array profiles with moderate level of normal contamination, low complexity and a good contrast. The differences between the methods start to arise when a tumor sample or its measured profile have compromised quality. PICNIC has problems with normal contamination as it was designed for cell lines, GPHMM tolerated up to 90% of normal contamination well but might be sensitive to sub-clones, GAP usually overestimates the number of breakpoints (an additional smoothing procedure is needed) and looses recognition quality when normal contamination exceeds 70%.

Accurate evaluation of performance of several SNP array mining approaches with respect to the normal contamination and profile complexity (using simulated data) has been done recently by an independent group [START_REF] Mosen-Ansorena | Comparison of methods to detect copy number alterations in cancer using simulated and real genotyping data[END_REF]. The GAP method showed the best performance among all approaches tested. We consider that the major advantage of the GAP method is its visualization based strategy which allowed setting up an automatic recognition procedure and fitting a number of hidden parameters based on a large set of tumor genomes and under clear-cut manual control.

Next Generation Sequencing (NGS)

The recently developed Next generation sequencing (NGS) technology -through whole-genome, wholeexome and whole-transcriptome approaches -provides the most complete view on the cancer genome and has a potential to uncover all spectrum of somatic variation, deduce tumor history and sub-clonal structure [START_REF] Greenman | Estimation of rearrangement phylogeny for cancer genomes[END_REF][START_REF] Meyerson | Advances in understanding cancer genomes through second-generation sequencing[END_REF]. Extracting absolute copy numbers belongs to the "must have" package for future NGS processing pipelines providing cancer genome characterization from point mutations to large-scale chromosomal rearrangements. Here, we want to outline current methodological and technological issues related to the large-scale absolute copy number and allelic content recognition based on NGS data.

NGS technology consists of shearing tumor genomic DNA into small fragments (200-4000bp), amplifying, and reading fragments using Illumina, SOLiD, PacBio or IonTorrent Analyzers [START_REF] Metzker | Sequencing technologies -the next generation[END_REF]. The single-end or paired-end reads are mapped back to the reference genome, providing coverage of each genomic position by a number of short fragments (reads). Window-averaged number of reads along the genome characterizes CNV profile; reads covering known annotated SNP positions could be used to obtain BAF profile [START_REF] Nielsen | Genotype and SNP calling from next-generation sequencing data[END_REF]. As soon as CNV and BAF profiles are obtained, the same strategy of tumor genome annotation that is described for SNP arrays could be applied [START_REF] Boeva | Control-FREEC: a tool for assessing copy number and allelic content using next-generation sequencing data[END_REF]. Due to potentially linear relationship between the number of reads and actual DNA copy number, segmented CNV and BAF profiles should follow the pattern shown in Figure 6A, where coefficient q is no more necessary and dynamics of copy number layers itself could provide the level of normal contamination (see equation ( 3)) [START_REF] Boeva | Control-free calling of copy number alterations in deep-sequencing data using GC-content normalization[END_REF][START_REF] Gusnanto | Correcting for cancer genome size and tumour cell content enables better estimation of copy number alterations from next-generation sequence data[END_REF].

However, at the moment a number of difficulties exist in extracting tractable CNV and BAF profiles from NGS data. Firstly, the number of reads along the genome depends heavily on regional GC content and read mappability, which result in large unspecific variations in CNV profile. Secondly, the reference allele usually has more reads due to asymmetric mapping procedure and BAF calculation results in a very noisy profile. Several normalization strategies using a matched normal sample or correcting for GC-content and/or mappability have been suggested [START_REF] Boeva | Control-free calling of copy number alterations in deep-sequencing data using GC-content normalization[END_REF]. In the case of targeted sequencing (for example, whole exome sequencing), in addition to mappability and GC-content, a bias resulting from uneven capture in different targeted regions should also be corrected [START_REF] Li | CONTRA: copy number analysis for targeted resequencing[END_REF][START_REF] Sathirapongsasuti | Exome sequencing-based copy-number variation and loss of heterozygosity detection: ExomeCNV[END_REF]. Mapping reads on the genome indexed with SNPs and exclusion of low quality reads before BAF calculation could improve genotype prediction.

Alternatively, huge numbers of reads provided by NGS could help annotating alterations by sequencing-specific approaches. In the case of paired-end sequencing, predicted alterations can be supported and clarified by the fragments encompassing the junction, which could be detected using abnormal reads mapping [START_REF] Medvedev | Detecting copy number variation with mated short reads[END_REF]. Moreover, "split-read" mapping (i.e. partial mapping of reads encompassing a junction) can provide base-pair resolution of alterations [START_REF] Wang | CREST maps somatic structural variation in cancer genomes with base-pair resolution[END_REF].

To conclude, accurate evaluation of copy numbers and allelic contents directly from the read counts is still not entirely implemented into powerful bioinformatics tools; mostly because of the short time frame since NGS data are available for the scientific community. Bioinformatics developments and NGS technology itself advance at a staggering rate that will undoubtedly result in new computational solutions in the future.

Extracting biologically relevant information from a recognized tumor genome profile represents an important issue. Here we address the detection of tumor ploidy and annotation of gains and losses with the example of breast tumor genomes.

According to the genomic content of a tumor cell, one could distinguish diploid, tetraploid, neardiploid, over-diploid, near-tetraploid or simply aneuploid (highly differing from diploid or tetraploid) status. We have suggested a way to ascribe tumor ploidy status based on the actual genomic content distribution in a large cohort of breast tumors. Inferred DNA index in the cohort of 250 breast tumors showed a bimodal distribution (Figure 8B) similar to those demonstrated for the genomes of various types of cancers [START_REF] Storchova | The consequences of tetraploidy and aneuploidy[END_REF]. This observation supports the hypothesis of frequent duplication of the whole genome during cancer progression, which explains bimodality [START_REF] Carter | Absolute quantification of somatic DNA alterations in human cancer[END_REF][START_REF] Storchova | The consequences of tetraploidy and aneuploidy[END_REF]. Two modes of the DNA index in the breast cancer cohort were found to be 1 and 1.7 units with the cut-off distinguishing two modes approximately equal to 1.3 (1 corresponds to the normal diploid genome). Thus, tumor genomes with DNA index less than 1.3 were considered to have a ploidy of two and are called "near-diploid genomes"; tumor genomes with DNA index equal or higher than 1.3 were considered to have a ploidy of four and are called "near-tetraploid genomes". Distribution of genomic states in near-diploid and near-tetraploid genomes was generally consistent with attributed ploidy. Near-diploid genomes had genomic DNA mainly presented in 1, 2, 3 copies (B, AB, ABB segmental genotypes) while near-tetraploid genomes had the well represented 2, 3, 4, 5-copy layers (BB, ABB, AABB, AABBB genotypes). Moreover, a significant bias observed in distribution of genomic states (such as high frequency of BB compared to AB states) in the small series of near-tetraploid basal-like breast tumors evidences the late (after a number of alteration events) whole genome duplication [START_REF] Popova | Genome Alteration Print (GAP): a tool to visualize and mine complex cancer genomic profiles obtained by SNP arrays[END_REF]. The same conclusion for other types of tumors has been published recently [START_REF] Carter | Absolute quantification of somatic DNA alterations in human cancer[END_REF].

Relative genomic alterations such as gains, losses, amplifications should be called according to the ploidy of a tumor sample. We suggested the following procedure: for near-diploid tumors copy loss, gain, amplification is called for the segments with ≤ 1, ≥ 3 and ≥ 4 copies, respectively. For near-tetraploid tumors copy loss, gain, amplification is called for the segments with ≤ 2, ≥ 6 and ≥ 8 copies, respectively. Thus, genomic segments in either of 3, 4, 5-copy states are considered to be not altered in a near-tetraploid genome. Indeed, according to the hypothesis of the late whole genome duplication, 3 and 5 copies represent "secondary" alteration events. Moreover, small relative changes in copy numbers are unlikely to strongly affect gene expression. An additional evidence to support our approach is the similar rate of gains and losses obtained for near-diploid and near-tetraploid basal-like breast tumors. With this approach the most of known recurrent alterations in total set of basal-like breast tumors displayed 80-100% recurrence rate (as compared to 50-60% reported based on the array CGH).

It is worth noting that genetic pathways in other types of tumors could be different from that described for breast cancers. For example, neuroblastoma is known to frequently display genomic content close to the near-triploid state (DNA index of 1.5), implying a high rate of whole chromosome gains or some other mechanisms of genome transformation [START_REF] Brodeur | Neuroblastoma: biological insights into a clinical enigma[END_REF]. For this reason, DNA index cut-off for ploidy attribution (as well as annotation of gain and losses) obtained for the breast tumors might not be suitable for the other types of tumors. However, we believe that the way of annotating genetic alteration should reflect (as far as possible) tumor genome evolution.

In this chapter we present an approach for mining genetic alterations measured by the SNP-based techniques. The main advantages of the GAP method are its visual and contextual simplicity and natural interpretation. Moreover, results of automatic recognition could be easily monitored manually, which increases the confidence in evaluating important cases. Systematic investigation, description and comparative analysis of cancer genomes by SNP arrays and NGS are currently ongoing. Clinical trials have been designed to translate some of the genomic findings into clinic to improve cancer diagnostics and treatment. We hope to see encouraging results in the near future.

Figure 3 :

 3 Figure 3: Single Nucleotide Polymorphism (SNPs) in individual genome and in population.Status of SNP in individual genome is defined as homozygous (two identical nucleotides) or heterozygous (distinct nucleotides). Single nucleotide variation is designated as a SNP, if its minor allele frequency exceeds 1% in at least one population; otherwise it is called variant or mutation.

Figure 4 :

 4 Figure 4: All possible allelic contents for DNA copy number from 1 copy to 5 copies. A. Pattern of genomic states if allelic content is characterized by B allele frequency (BAF) (Illumina); B. Pattern of genomic states if allelic content is characterized by Allelic difference (AD) (Affymetrix SNP 6.0, CytoScanHD). C. Evaluation of BAF and AD from CN and MA. CN: copy number; MA: major allele counts; homozygous genomic states are shown in red; position of each pair (CN, MA) is represented on the graphs by the small squares (MA: x axis; CN: y axis); red rectangles indicate homozygous genomic states (LOH); only CN ≤ 5 are shown, however, the structure could be easily extended to the higher CN levels.

Figure 5 :

 5 Figure 5: The whole genome profiles of genomic rearrangements for one primary breast tumor measured on two platforms and their GAP plots. A. Illumina 300K SNP array profiles (left) and corresponding GAP plot (right): Copy Number Variation (CNV) profile is represented by log R ratio (LRR), Allelic Imbalances (AI) are represented by B allele frequency (BAF); B. Affymetrix SNP 6.0 profiles (left) and corresponding GAP plot (right): CNV profile is represented by Log 2 Ratio (L2R), AI are represented by Allelic Differences (AD). GAP plots display combined sideview projections of segmented CNV and AI profiles; each region of genome is represented by two symmetrical circles; position of a circle is defined by median of CNV (y axis) and (lower) mode of AI (x axis); the size of a circle is proportional to the segment size; rectangles indicate the copy number / allelic content patterns.

Figure 6 :

 6 Figure 6: Model templates and three parameters modifying the model (indicated by the blue arrows (p,q) and stars (C)) to fit to a tumor GAP. A. Template for Illumina platform. B. Template for Affymetrix platform. p is normal contamination proportion; q is coefficient of CNV contraction; the stars correspond to 2 copy centering position on CNV scale (parameter C). Homozygous states are designated by red quadrants.

Figure 7 :

 7 Figure 7: GAP plots and results of automatic recognition procedure of copy number and allelic content in three primary breast tumors. The best fitting model template is shown by colored rectangles, designated by the ratio: CN/MA. A. and B. GAPs and recognition templates for the breast tumor sample shown in Figure 5 on Illumina (A) and Affymetrix (B) platforms; C. and D. Examples of model fitting for two over-diploid breast tumors on Illumina (C) and Affymetrix (D) platforms.

Figure 8 :

 8 Figure 8: DNA indexes for the breast tumor genomes. A. Validation of copy number attribution based on the GAP method: each point represents DNA index of a primary breast tumor measured by SNP array (x-axis) and by flow cytometry (FCM) (y-axis). B. Distribution of DNA indexes obtained from the SNP arrays of a big cohort of breast tumors.

Acknowledgements

The authors would like to thank Dr. Amaury Dumont, Jaydutt Bhalshankar and reviewers for productive comments and accurate proofreading of the manuscript.

The latest evaluation of the SNP-arrays on the Illumina platform is essentially similar to all previous versions (only the number of SNPs was increased and significant number of CN probes were included;

CN probes are intended to measure germline CNV s, which we are not considering here). All conclusions presented here for Illumina platform equally concern all versions of Illumina SNP arrays. In contrast, the latest versions of Affymetrix platform, Affymetrix SNP 6.0 and CytoScan HD, producing higher quality data, compared with the previous version of SNP-arrays (Affymetrix 100K, 250K, 500K) were subjected to specific normalization. All conclusions presented here were evaluated for Affymetrix SNP 6.0 platform and successfully tested on CytoScan HD with some minor modifications. Normalization of previous versions of Affymetrix (except SNP chip 5.0) could be obtained by Aroma package http://www.aroma-project. org and should be treated within the Illumina framework.

For the Illumina platform CNV profile is represented by the Log R ratios (LRR), which are the logtransformed ratios of experimental and normal reference SNP intensities, centered at zero for each sample. Allelic imbalances are represented by BAF, which are the normalized proportions of the B allele signal. Normalization procedure designed by Illumina produces not symmetrical BAF profiles, thus, we recommend to symmetrize BAF profile by applying the tQN algorithm for quantile normalization (Staaf et al., 2008a).

For the Affymetrix platform CNV profile is represented by the Log 2 ratios (L2R), which are the logtransformed ratios of experimental and normal reference SNP intensities. Allelic imbalances are represented by the AD, which represents the difference of A signal and B signal each standardized with respect to their median values in the reference.

Appendix B CNV and AI Profiles Segmentation

The circular binary segmentation (CBS) algorithm (DNAcopy package, Bioconductor) [START_REF] Venkatraman | A faster circular binary segmentation algorithm for the analysis of array CGH data[END_REF] was used for segmentation of all profiles presented here. There exist numerous methods of optimal CNV profile segmentation, and any of them could be used instead of CBS. However, the GAP method tolerated "false positive" breakpoints in either profile well, while "false negative" breakpoints could be misleading for pattern recognition. Thus, segmentation method should be as sensitive as possible but could tolerate less specificity.

BAF and AD profiles are supposed to be symmetric, thus both profiles were transformed in order to reduce redundancy: BAF = 0.5 + abs(BAF -0.5) and AD = abs(AD).

BAF and AD profiles represent two mode distributions in any genomic segment with heterozygous genotype (in tumor or germline): lower mode corresponds to heterozygous SNPs (informative) and higher mode corresponds to homozygous SNPs (non-informative). Because segmentation of a two mode profile is less trivial than that of a one mode profile, we attempted to reduce BAF and AD profiles to one mode. In the case of BAF, non-informative homozygous SNPs do not depend on copy number and could be filtered out based on the threshold (BAF > 0.97, Illumina 300K, as suggested, for example, in (Staaf et al., 2008a)).

In the case of AD, non-informative homozygous SNPs depend on copy number and simple threshold is not applicable. However, due to the fact that AD values for homozygous genomic states equal to the corresponding copy numbers, we applied AD filtering based on the L2R profile. To filter out homozygous mode, we defined a profile, representing a joint neighborhood of 2 L2R . Those points in AD profile which fall into the joint neighborhood of 2 L2R were filtered out.

Filtered AI (BAF or AD) profiles were segmented by CBS algorithm. CNV and AI breakpoints were united, providing segmentation of the genomic profile into alteration units with presumably stable copy number and allelic content state.

Appendix C Algorithm of Automatic Recognition of the GAP Pattern

1. Consider GAP of a tumor sample to be a set of genomic segments (alteration units, obtained after segmentation of CNV and AI profiles) each characterized by the three values: the median of CNV , the (lower) mode of AI, and the length in SNP counts:

2. For each p (proportion of normal contamination), q (coefficient of experimental variation), and model centering C (corresponding to the position of 2 copies on the CNV axis) we define a model GAP

where CN = 1,...,5; MA = n c B = [CN/2],...,CN; n n B = 1, if n c B < CN ; n n B = 1,2 , if n c B = CN ; approximation for L2R values was found experimentally, by considering distributions in a large set of Affymetrix SNP 6.0 arrays. Each (CNV CN ,AI CN/MA ) pair is included into the GAP T with a certain rectangle neighborhood, providing a model recognition template (similar to ones in Figure 6).

3. We define a grid in the parameter space: p ∈ [0;0.7](2p ∈ [0;0.7] for GAP T Affymetrix ), q ∈ [0.1;1] , C ∈ [-1;1], with the step = 0.05. 4. For each set of parameters (p,q,C) we calculate goodness-of-fit criterion: K p,q,C = ∑ i,(CNV i ,AI i )∈GAP T L i , (CNV i ,AI i ) ∈ GAP T if (CNV i ,AI i ) falls into any model rectangle defined by (CNV CN ,AI CN/MA ) pair from GAP T . 5. We order parameter sets by the criterion value and choose the interpretation with maximal goodnessof-fit: (p * ,q * ,C * ) = arg max (K p,q,C ), after filtering out "redundant" interpretations (in particular, such that superimposition of GAP T and tumor GAP results in empty acquired homozygosity states, etc). This procedure results in approximately 80% rate of correctly recognized models.

6. We annotate the segments from the tumor GAP by the closest CN and MA from the model template (for tumor annotation we enlarge the model GAP T to 8 copies).

7. We smooth resulting CN/MA profiles and exclude redundant breakpoints. R scripts are available at http://bioinfo-out.curie.fr/projects/snp_gap. Due to regular updating, some details or constants in implemented code may be slightly different from those described above.


Ivan Obradović
email: ivano@afrodita.rcub.bg.ac.yu

Cvetana Krstev
email: cvetana@matf.bg.ac.yu

Gordana Pavlović-Lažetić
email: gordana@matf.bg.ac.yu

Duško Vitas
email: vitas@matf.bg.ac.yu

Corpus Based Validation of WordNet Using Frequency Parameters

In this paper we define a set of frequency parameters to be used in synset validation based on corpora. These parameters indicate the coverage of the corpus by wordnet literals, the importance of one sense of a literal in comparison to the others, as well as the importance of one literal in a synset in comparison to other literals in the same synset. The obtained results can be used in synset refinement, as well as in information retrieval tasks.

Introduction

The main goal of BalkaNet, the Balkan wordnet project (BWN) is the development of a multilingual database with wordnets for a Bulgarian, Czech, Greek, Romanian, Serbian and Turkish [START_REF] Stamou | BALKANET: A Multilingual Semantic Network for Balkan Languages[END_REF]. In its initial phase, Balkanet followed the approach similar to that of EuroWordNet (EWN) developing monolingual wordnets interconnected through an interlingual index (ILI). The development of databases started with a translation of a common set of concepts named Base Concepts in EWN, using the Princeton WordNet (PWN) as the source.

The Serbian WordNet (SWN) has been developed according to this common approach. In the absence of both an explanatory dictionary and an English/Serbian dictionary in electronic form, the translation of English synsets from PWN was done manually, while preserving the PWN semantic structure. The fact that a Serbian dictionary of synonyms does not exist even in paper form made this task even more difficult. In order to establish a relation with the reference six volume explanatory Serbian dictionary of Matica Srpska (RMS), the senses attributed in SWN to literals, or words representing sysnset lemmas in general correspond to the ones in this dictionary. Since the RMS dictionary was published in 1971, new senses had to be attributed in SWN to some of the existing literals but also new literals had to be added. Another reason for refinement of senses defined by RMS is due to the fact that concepts, and hence literal senses in PWN are far more fine grained that the ones in RMS.

The conditions under which SWN has been developed brought up the question of validation of Serbian synsets on corpora. The idea to semantically tag corpora with senses from WordNet has first been realized within the SemCor project [START_REF] Miller | Using a semantic concordance for sense identification[END_REF]. The use of monolingual and multilingual corpora for synset validation leading to the introduction of new literals or removal of existing ones from a synset has already been tackled in [START_REF] Krstev | Corpora Issues in Validation of Serbian Wordnet[END_REF][START_REF] Obradović | Application of Intex in Refinement and Validation of Serbian Wordnet[END_REF]. Further refinement of this approach is presented in this paper. In order to establish more precise criteria for synset validation a set of numerical parameters related to literal-sense pair frequency in corpora has been developed.

Frequency Parameters

In order to evaluate the quality of a synset in terms of the comprehensiveness and adequacy of literals used for the lexicalization of a particular concept on one hand, and to establish an ordering among literals within a synset which may be used in information retrieval tasks, on the other, we define a set of indices as numerical measures of relevance of particular literals to synsets they are used in.

Let S be the finite set of all synsets within a wordnet: S = {S i |S i is a synset describing a specific concept, i = 1, 2, . . . , N S }; N S is the total number of synsets within a wordnet. Let L be the finite set of all literals used as lexicalizations of one or more concepts: L = {L k |L k is a literal used in at least one synset of the wordnet, k = 1, 2, . . . , N L }; N L the total number of different literals used within the wordnet. When a literal L k ∈ L is used as a lexicalization of a specific concept described by the synset S i , it is used in a specific sense (a sense tag is attached to the literal). Omitting the index k of the literal we shall mark all literal-sense pairs within a nonempty synset S i ∈ S in a sequence as L S i j ( j = 1, 2, . . . , n i ), where n i ≥ 1 is the total number of literals within the nonempty synset S i .

We shall define the indices for literals within the wordnet, with the aim to determine the relevance of a particular literal to a synset it is used in. In order to determine these indices for a literal a search is performed on a corpus and all occurrences of the selected literal as well as its inflectional forms are identified within a context of a predefined length. We shall first denote the total number of occurrences of a literal L k within the corpus, regardless of its sense, as L C k . The next step is a time-consuming one since it requires manual identification of the sense in which the literal has been used in every particular concordance line identified in the corpus. When this task is completed then the number of occurrences of a literal within the corpus in each specific sense is established. For the senses covered by the wordnet, the appropriate synset S i the literal belongs to can then be identified. We then proceed taking into account only these senses, and denote the number of times the literal L k has been used for lexicalization of a concept described by the synset S i as L S C i j . The sum of these numbers obtained for all possible senses of a literal covered by the wordnet yields L W N k , namely, the number of cases when a literal has been used within the corpus as a lexicalization of a concept represented in the wordnet. It is clear that

L W N k ≤ L C
k , and that the target of each wordnet should be that for all literals, ideally,

L W N k = L C
k holds. This would mean that all possible sense usages of a literal identified within the corpus have been covered by wordnet synsets.

If we want to express the relevance of a particular literal L k to a particular synset S i within a corpus, then we should compare the number of occurrences of this literal in the corpus denoting the concept represented by the synset S i , that is L S C i j , to the total number of occurrences of this literal within the corpus, namely L C k . Thus we define the overall synset relevance index of a literal as the ratio of the number of times this literal has been used in a specific sense and the total number of occurrences of this literal in the corpus, namely:

I C ik = L S C i j /L C k
where the literal from L S i j equals the literal L k . The index range is 0 < I C ik ≤ 1, where I C ik = 1 means that the literal L k is used in one and only one sense, and that is to lexicalize the concept described by the synset S i .

Since the wordnet coverage of the senses of a literal does not always have to be complete, we define the wordnet synset relevance index as the relevance of a particular literal L k to a particular synset S i within a more restricted part of the corpus, that is, the part already covered by the wordnet. This index is defined as the ratio of the number of times this literal has been used in a specific sense and the total number of occurrences of a literal within the corpus denoting concepts represented in the wordnet (L W N k ), namely:

I W N ik = L S C i j /L W N k ,
where the literal from L S i j is the literal L k . As is the case with I C ik , the index range is

0 < I W N ik ≤ 1
I W NC k = L W N k /L C k .
The index ranges between 0 and 1, and in case of full coverage is equal to 1.

All previous indices evaluated the relevance of a literal to a synset regardless of possible other literals within that synset. In order to compare the relevance of a literal within a synset in comparison to other literals denoting the same concept we define the local synset relevance index of the literal L k as the ratio of the number of occurrences of this literal in the corpus denoting the concept represented by the synset S i , that is L S C i j , and S C i , the number of occurrences of all literals denoting this concept (i.e. belonging to synset S i ):

I L ik = L S C i j /S C i , S C i = n i j =1 L S C i j .
It should be noted that the range of the index is 0 < I L ik ≤ 1 where I L ik = 1, holds when either the synset has only one literal, or other literals from that synset have not appeared in the corpus.

The Validation Procedure

In order to test this approach a subset of literal strings, that we called main strings has been chosen among those nouns and verbs that have the most senses in Serbian wordnet. Next, a subcorpus has been compiled consisting of contemporary newspaper texts comprising 1.7MW. Concordances were produced for all the inflectional forms of these nouns and verbs.

In the next step all the synsets in which the main strings appear have been identified, as well as literal strings, that we called supporting strings, that occur beside them in these synsets. For these supporting strings concordances have also been produced. The main and supporting literal strings form the "lexical sample" as defined by the SENSEVAL project [START_REF] Kilgarriff | English SENSEVAL: Report and Results[END_REF]. The produced concordances (around 10.000) have than been manually analyzed by lexicographers. In the first step the concordance lines containing the homograph forms have been rejected. In the remaining lines the senses have been identified according to the RMS dictionary and SWN, and marked using the same sense labels.

On the basis of the obtained results tables have been produced and the indices introduced in the section 2 calculated. These data for the noun lice and the verb proizvesti are given in Tables 1 and2. For each of the main strings only the senses that are present in SWN are represented. The frequency of occurrence of the these senses in the corpus is given in column

L S C i j . The row L W N k
represents the frequency of the occurrence of all the senses of a string that are covered by SWN, while the row other represents the frequency of the occurrence of those senses that are not yet covered. L C k is the sum of these two data, and represents the total frequency of the occurrence of the main string, while the index I W NC k represents their ratio. Among 12 main strings that have been analyzed, three had the value of this index 1, which means that for these strings all the senses identified in RMS dictionary (and perhaps some more) have been included in SWN. For all analyzed literals this index ranges from 0.246 to 1.

Table 2. The frequency indices for the lemma proizvesti obtained on newspaper corpus

Synset proizvesti

L S C i j

prouzrokovati:1 potaknuti:2x iznedriti:1 proizvoditi:3 napraviti:1a

S C i I C ik I W N ik I L ik produce:3,. . . 1a
I W NC k 0.985 I C ik I C ik I C ik I L ik I L ik
1.000 1.000 0.000 0.930 0.132

I L ik I L ik I L ik I L ik I L ik
0.816 0.026 0.000 0.570 0.113

The parameter S C i gives the overall occurrence of the synset, that is all its literals, in the corpus. The indices I C ik , I W N ik , and I L ik in the upper part of the table refer to the main string, while the same indices in the lower part refer to the appropriate supporting strings. The first one is the ratio L S C i j /L C k : for instance, for the sense 1a of the main string lice, this index is 0.063, which means that this sense represents 6.3% of all the occurrences of this string in corpus. The second index is the ratio L S C i j /L W N k . For the same sense of the string lice its value is 0.085 meaning that it covers 8.5% of all the occurrences that represent senses from SWN. Finally, the third index is the ratio L S C i j /S C i . For the sense 2a of the string lice the value of this index is 0.026, meaning that of all occurrences of this synset, 2.6% were represented by this particular literal.

If for some string the value of its index I L ik is close to 0 it can indicate that it has been misplaced in the synset, especially in the cases when both indices L C k and S C i are considerably greater than 0. For instance, that is the case for the string napraviti:1a (Table 2). The string napraviti has a considerable frequency on corpus (L C k = 159), and the synset to which the literal string napraviti:1a belongs also has a considerable frequency (S C i = 186).

However, its local synset relevance index is relatively low (I L ik = 0.113), and the synonymy of the literal string napraviti:1a with the main string proizvesti should be reconsidered.

The calculated indices enable the ordering of the literal strings in a synset. This can be useful for information retrieval (IR) tasks that are seen as one of the most interesting applications of BWN. Especially, strings that have a low value of I L and a high value of I C and which are not necessarily misplaced in a synset, can be neglected in IR tasks, thus reducing the recall but improving the precision. In order to test the impact of the nature of the corpus to index values the validation procedure was performed on a small literary corpus of 0.5 MW for aselected number of literals. The results obtained show that the index values can be largely affected by the nature of the corpus. Thus, for example, the values of both I C ik and I W N ik have dramatically changed for senses 1a and 2a of the noun lice (Table 3). This does not come as too much of a surprise since meaning 2a ("A part of a person that is used to refer to a person") is more used in newspaper texts whereas the meaning 1a ("The front of the human head. . . ") in literature.

The changes seem to be far less dramatic for the indices I L ik , but in order to draw some final conclusions the literals should be tested on a larger corpus.

Conclusion

The applied procedure confirmed the importance of the validation of synsets on a corpus. The adequacy of placement of each literal and its sense in a synset can not be fully assessed without analyzing its appearances in the concordance lines. The frequency indices can serve as useful numerical indicators in this assessment procedure. However, to get a fair estimate of a literal in terms of these parameters, the procedure needs to be applied on a large and balanced corpus. To that end automatic or/and semi-automatic procedures need to be developed in order to alleviate the time-consuming task of manual concordance analysis.

  , where I W N ik = 1 means that the literal L k is used in one and only one sense. Since L W N k ≤ L C k , then I W N ik ≥ I C ik . As, ideally, L W N k = L C k should hold for every literal, in an ideal case I W N ik = I C ik should also be true. In order to evaluate how close a particular literal L k is to the ideal case, namely when all its possible senses are covered by the wordnet, we should compare the number of occurrences of a literal within the corpus denoting concepts represented in the wordnet L W N k to the total number occurrences of the literal within the corpus L C k . We therefore define the wordnet coverage index of a literal L k , namely

Table 1 .

 1 The frequency parameters for the lemma lice obtained on newspaper corpus

	Synset	lice L S C i j uloga:1a lik:3 strana:1b S C i	I C ik I W N ik	I L ik
	face, human face	1a	33	*	*	* 33 0.063 0.085 1.000
	face:6	2a 353	*	*	* 353 0.675 0.912 1.000
	character:4, role:2,	2b	1	34	3	* 38 0.002 0.003 0.026
	theatrical role:1,. . .,
	face:14	3	0	*	*	*	0 0.000 0.000 0.000
	side:5,	5a	0	*	*	5	5 0.000 0.000 0.000
		L W N k	387
		other 136					0.260	*	*
		L C k I W NC k	523 0.740	298 I C ik	20 I C ik	861 I C ik
				0.114 0.150	0.006
				I L ik	I L ik	I L ik
				0.895 0.079	1.000

Table 3 .

 3 The frequency parameters for the lemma lice obtained on literary corpus

	Synset	lice L S C i j uloga:1a lik:3 strana:1b S C i	I C ik I W N ik	I L ik
	face, human face	1a 380	*	*	* 380 0.936 0.979 1.000
	face:6	2a	3	*	*	*	3 0.007 0.008 1.000
	character:4, role:2,	2b	3	6	1	* 10 0.007 0.008 0.300
	face:14	3	0	*	*	*	0 0.000 0.000 0.000
	side:5,	5a	2	*	*	4	6 0.005 0.005 0.333
		L W N k	388
		other	18					0.044	*	*
		L C k I W NC k	406 0.956	22 I C ik	25 I C ik	287 I C ik
				0.273 0.040	0.014
				I L ik	I L ik	I L ik
				0.600 0.100	0.667