Thierry Jauniaux
email: t.jauniaux@ulg.ac.be

Frédéric Farnir

Michaël Fontaine

Jeremy Kiszka

Michael Sarlet

Freddy Coignoul

Cytochrome P450 1A1 expression in cetacean skin biopsies from the Indian Ocean

Keywords: Cytochrome P450, Biopsy, Skin, Immunohistochemistry, Mayotte

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

Abstract

The study describes cytochrome P-450 1A1 (CYPA1) expression in the skin of different cetacean species (Megaptera novaeangliae, n= 15; Stenella attenuata, n=7 and Stenella longirostris, n= 24) from the Mozambique Channel island of Mayotte. Immunohistochemical examination was performed with a monoclonal antibody against scup cytochrome CYPA1.

The sex was determined using a molecular approach consisting in the genotyping sex-specific genes. CYPA1 was detected at the junction between epidermis and blubber on dolphins only, mostly in the endothelial cells. Similar observation was obtained in the dermis of one M. novaeangliae. Immunohistochemical slides were scored to evaluate the expression of the CYPA1 and a higher expression was observed in S. longirostris, suggesting a higher exposure to pollutants for this species. The difference of expression between sexes was not significant.

Cytochrome P450 1A proteins (CYP1A1) are induced by exposure to various environmental contaminants such as polycyclic aromatic hydrocarbons (PAHs) and planar halogenated aromatic hydrocarbons (PHAHs) (e.g., non-ortho-polychlorinated biphenyls -PCBs), activating the aryl hydrocarbon receptor. These contaminants, present worldwide and in all oceans, accumulate in the marine food web. High concentrations accumulate in the tissues of marine top predators, such as the blubber of marine mammals (Aguilar et al., 2002;[START_REF] Pierce | Bioaccumulation of persistent organic pollutants in female common dolphins (Delphinus delphis) and harbour porpoises (Phocoena phocoena) from western European seas: geographical trends, causal factors and effects on reproduction and mortality[END_REF]. In addition, these stable contaminants are suspected to have adverse health effects mostly in marine top predators, especially on the immune and reproductive systems [START_REF] Hall | Individual-based model framework to assess population consequences of polychlorinated biphenyl exposure in bottlenose dolphins[END_REF][START_REF] Reijnders | Reproductive failure in common seals feeding on fish from polluted coastal waters[END_REF]. In marine mammals such as cetaceans, stable contaminants are reported to induce sensitivity to infectious diseases, such as in the harbour porpoise (Phocoena phocoena) from the North-east Atlantic with high PCBs level. Indeed, high PCB levels have been recorded in harbour porpoises dying from an infectious disease compared with individuals dying from other causes on animals stranded on the coasts of the UK, on the continental southern North Sea coast [START_REF] Jepson | Relationships between polychlorinated biphenyls and health status in harbor porpoises (Phocoena phocoena) stranded in the United Kingdom[END_REF][START_REF] Pierce | Bioaccumulation of persistent organic pollutants in female common dolphins (Delphinus delphis) and harbour porpoises (Phocoena phocoena) from western European seas: geographical trends, causal factors and effects on reproduction and mortality[END_REF] and for the German North and Baltic Seas coasts [START_REF] Beineke | Investigations of the potential influence of environmental contaminants on the thymus and spleen of harbor porpoises (Phocoena phocoena)[END_REF]. Morbillivirus epizootics in the harbour seal (Phoca vitulina) and striped dolphin (Stenella coeruleoalba) populations have been suspected to be associated with high level of PCBs (Aguilar and Borrell, 1994;[START_REF] Hall | Organochlorine levels in common seals (Phoca vituina) which were victims and survivors of the 1988 phocine distemper epizootic[END_REF].

CYP1A1 is expressed in all animal tissues including skin, liver, urine bladder and lung, and immunohistochemical investigation has be done to compare expression in tissues from beluga whales [START_REF] Wilson | Systemic effects of arctic pollutants in beluga whales indicated by CYP1A1 expression[END_REF]. In the integument, analysis of CYP1A1 expression deserves attention as a potential non-lethal technique to study the impacts of contaminants on cetaceans and such expression was described by immunohistochemistry in the integument of 17 species (six mysticetes and 11 odontocetes), from various regions of the world [START_REF] Angell | Cytochrome P450 1A1 expression in cetacean integument: implications for detecting contaminant exposure and effects[END_REF]. [START_REF] Wilson | Correlates of cytochrome P450 1A1 expression in bottlenose dolphin (Tursiops truncatus) integument biopsies[END_REF] correlated it with the PCB concentrations from bottlenose dolphins (Tursiops truncatus) resident in the Saratosa Bay. The detection of CYP1A1 has been also investigated on Mediterranean cetaceans [START_REF] Fossi | First detection of CYP1A1 and CYP2B induction in Mediterranean cetacean skin biopsies and cultured fibroblasts by Western blot analysis[END_REF]. Very little information is available on the exposure of cetaceans to contaminant and potential consequences the in Indian Ocean Sanctuary (Aguilar et al., 2002). The aim of the present study is to describe the skin histology and the cytochrome P450 1A1 expression in cutaneous biopsies of different cetacean species collected in the Mayotte waters. We hypothesised that a multi-species approach would help to understand interspecies variability of CYP1A1 expressions, with higher trophic level species having a higher expression than low trophic level species. Based on existing information in the literature and local information, we considered three species in the present study: the pantropical spotted dolphin (Stenella attenuata SA), the spinner dolphin (Stenella longirostris SL) and the humpback whale (Megaptera novaeangliae MN). The delphinids have the higher trophic level, feeding on epipelagic fishes and cephalopods [START_REF] Norris | Behavior of the Hawaiian spinner dolphin, Stenella longirostris[END_REF], while southern hemisphere humpback whales mostly feed on krill from high latitudes, being of a lower trophic level. Among the two delphinids considered in this study, the pantropical spotted dolphin has a higher trophic position than the spinner dolphin [START_REF] Gross | A preliminary study of habitat and resource partitioning among co-occurring tropical dolphins around Mayotte, southwest Indian Ocean[END_REF].

Mayotte (12 o 50'S, 45 o 10'E), the sampling area, is located in the northern Mozambique Channel between Madagascar and the African mainland, entirely surrounded by a 197 km long barrier reef. A high marine mammal diversity, 22 species including 12 delphinids [START_REF] Kiszka | Structure of a toothed cetacean community around a tropical island[END_REF] have been reported around Mayotte, the most common species being the spinner dolphin, the pantropical spotted dolphin, the Indo-Pacific bottlenose dolphin (Tursiops aduncus) and the melon-headed whale (Peponocephala electra); these are resident year-round [START_REF] Kiszka | Structure of a toothed cetacean community around a tropical island[END_REF]. The shallow and protected waters of the waters around the Comoros, including Mayotte, provide ideal habitat characteristics for wintering humpback whales during their reproductive cycle (Ersts et al., in press;[START_REF] Kiszka | Structure of a toothed cetacean community around a tropical island[END_REF].

Skin biopsies (epidermis and underlying blubber) were collected from 46 free ranging cetaceans (MN, n = 15; SA, n = 7 and SL, n = 24) from Mayotte waters, during austral winter from 2004 to 2006. The biopsies were collected using a crossbow (BARNETT Veloci-Speed® Class, 68-kg draw weight) with Finn Larsen bolts and tips (dart 20 and 40-mm long, 5-mmdiameter). The largest tips were used for humpback whales, while the smallest were used for delphinids. Animals were hit below the dorsal fin and biopsies were preserved in 90° ethanol prior to analyses. For the delphinids, only adult individuals were sampled, whereas several calves were samples in humpback whales for adult/calves comparisons. For subsequent analyses, biopsies were divided in two: one for microscopical examination, one for sex determination.

For microscopy, skin samples were embedded in paraffin wax by routine procedures and sections (3 µm) were cut and stained with haematoxylin and eosin (HE). In addition, an immunohistochemical examination with a monoclonal antibody against the scup CYP1A1 was performed on paraffin embedded sections as described by [START_REF] Angell | Cytochrome P450 1A1 expression in cetacean integument: implications for detecting contaminant exposure and effects[END_REF]. The other immunohistochemical reagents used were those included in the Enhanced Polymer One-Step Staining Procedure (Dako Envision™, Dako, Glostrup, Denmark). Test sections, on which the first layer antibody was omitted, served as negative controls.

For every microscopic slide, a quality index (Q.I.) was given from 1 (well preserved and well differentiated tissue) to 3 (altered tissue). For the immunohistochemical examination, a staining score (S.S.) was given: 0 (absence of staining); 1 (slight or focal positive staining); 2 (moderate: mild and multifocal) and 3 (intense and multiple staining). The slides were scored in three different skin levels: epidermis, junction between epidermis and blubber and finally blubber. Statistical tests were used to evaluate the effects of species and gender. Due to the non-normality of the scores, permutation statistical tests were performed: the obtained difference between the score averages for males and females were compared to the differences distribution obtained by shuffling 1000 times the scores between males and females within each of the 3 species: the obtained p-value is the rank of the observed difference in this distribution. A similar procedure was used to test for species differences [START_REF] Churgill | Empirical threshold values for quantitative trait mapping[END_REF].

PCR-based sex determination was performed using the four primers described by [START_REF] Rosel | PCR-based sex determination in Odontocete cetaceans[END_REF]. The PCR products were separated by electrophoresis on 2% agarose gels and gender was determined from the resulting banding patterns.

In total, 46 biopsies were investigated and the Q.I. was 1 for 35 samples, 2 for 10 samples and 3 for one. Examination with haematoxylin and eosin staining revealed no abnormality or lesions but four MN samples were discarded since only the superficial epidermis had been biopsied. Histologically, the junction between epidermis and blubber was rich of vascular structure. There was no microscopic difference between individuals or species except that the epidermis of MN was thicker than in the delphinids.

The immunohistochemistry (Table 1) revealed a slight positive staining (S.S.1) only observed in the epidermis of two SL, most obvious in the cytoplasm of basal cells. At the junction between epidermis and blubber, a very slight staining was observed in the cytoplasm of endothelial cells of two MN and one SA while the staining was slight on 12 SL, moderate on three and intense on two. CYP1A1 expression was stronger and most frequent in endothelial cells and rarely in the fibrocytes. Finally, in the blubber, only one MN, one SA and one SL were slightly positive. The maximum global scores were 1 for the MN, 2 for the SA and 5 for the SL, the average by species was 0.27 for MN, 0.6 for SA and 1.29 for SL.

The sex was determined on all cetaceans except on two SA. Among MN there were 12 females and 3 males (all calves); for SA, there were 2 females and 3 males and finally, there were 9 females and 15 males among SL.

For the statistics, four MN (only epidermis present) and two SA (without sexing) were not included in the analysis. Only the staining score of the junction between epidermis and blubber (endothelial cells) was considered and the permutation tests reveal that the difference was not significant among sexes (P = 0.766) but was very significant (P = 0.002) among species.

The CYP1A1 expression has been documented from North America humpback whales [START_REF] Angell | Cytochrome P450 1A1 expression in cetacean integument: implications for detecting contaminant exposure and effects[END_REF], but had never been reported in pantropical spotted dolphins and spinner dolphins. In addition, no data are available from cetaceans of the Indian Ocean. In the three species investigated, the expression of the CYP1A1 was strongest in endothelial cells at the junction between epidermis and dermis, followed by the dermis and very rarely in the epidermis. In addition, the expression was very mild in MN and in SA compared to SL. Such expression in the endothelial cells has been reported by Angel et al. (2004) and by [START_REF] Wilson | Correlates of cytochrome P450 1A1 expression in bottlenose dolphin (Tursiops truncatus) integument biopsies[END_REF] since CYP1A1 is highly inducible in the endothelium of vertebrates [START_REF] Stegeman | Cytochrome-P450ia1 Induction and Localization in Endothelium of Vertebrate (Teleost) Heart[END_REF]. Such induction in the integument endothelial cells of cetaceans may be related to the blood concentration of contaminants, agonist of the aryl hydrocarbon receptor [START_REF] Angell | Cytochrome P450 1A1 expression in cetacean integument: implications for detecting contaminant exposure and effects[END_REF]. In vitro similar inductions were reported on cetacean endothelial cells from lung and kidney after exposure to aryl hydrocarbon receptor agonists [START_REF] Garrick | Cytochrome P4501A is induced in endothelial cell lines from the kidney and lung of the bottlenose dolphin, Tursiops truncatus[END_REF].

The difference of expression between the three species investigated was statistically significant while the difference between sexes was not significant.

For [START_REF] Wilson | Correlates of cytochrome P450 1A1 expression in bottlenose dolphin (Tursiops truncatus) integument biopsies[END_REF], the CYP1A1 expression 1n bottlenose dolphins from Sarasota was neither age nor sex related, and the reproductive status did not relate to the levels of CYP1A1 in females. But, the endothelial CYP1A1 expression was positively correlated with the blubber total PCBs, although the mono-ortho PCBs concentration did not show this relationship. The contaminant concentrations appear to be stronger determinants for the expression of CYP1A1 in the skin biopsies, more than are age, sex, or reproductive status [START_REF] Wilson | Correlates of cytochrome P450 1A1 expression in bottlenose dolphin (Tursiops truncatus) integument biopsies[END_REF]. In this study, contaminant concentrations were not investigated on the blubber. As pollutants induce the expression of CYP1A1 and as SA and SL are resident in the Mayotte waters all year-round, they could have been contaminated locally through the food by the process of bio-magnification [START_REF] Pierce | Bioaccumulation of persistent organic pollutants in female common dolphins (Delphinus delphis) and harbour porpoises (Phocoena phocoena) from western European seas: geographical trends, causal factors and effects on reproduction and mortality[END_REF]. Despite the fact that the southern hemisphere has a low level of pollution, some restricted areas such as river mouths or semiclosed water masses are contaminated as a result of point source pollution (Aguilar et al., 2002). The difference in the CYP1A1 expression between the two dolphin species could be linked to a different position in the trophic web (J. Kiszka & L. Martinez, unpublished data), revealing a higher trophic level of the spinner dolphin. However, existing information on the trophic position of these two delphinid species suggest that spotted dolphins have a higher trophic level than spinner dolphins [START_REF] Gross | A preliminary study of habitat and resource partitioning among co-occurring tropical dolphins around Mayotte, southwest Indian Ocean[END_REF]. The lowest expression in humpback whales could be linked with the lowest trophic level of baleen whales in the marine food web compared to odontocetes [START_REF] Pauly | Diet composition and trophic levels of marine mammals[END_REF]. In addition, MN feed abundantly in Antarctica where the level of contamination is very low (Aguilar et al., 2002) while they starve when they are on their wintering ground in subtropical waters such as Mayotte.

This study demonstrates the expression of CYP1A1 in endothelial cells of cetaceans of the Mayotte waters. The expression could be induced by environmental contaminants that are also suspected to result in adverse health effects in cetaceans [START_REF] Hall | Individual-based model framework to assess population consequences of polychlorinated biphenyl exposure in bottlenose dolphins[END_REF][START_REF] Jepson | Relationships between polychlorinated biphenyls and health status in harbor porpoises (Phocoena phocoena) stranded in the United Kingdom[END_REF][START_REF] Pierce | Bioaccumulation of persistent organic pollutants in female common dolphins (Delphinus delphis) and harbour porpoises (Phocoena phocoena) from western European seas: geographical trends, causal factors and effects on reproduction and mortality[END_REF]. Further investigations will be helpful to identify the contaminants, their concentrations and to evaluate the potential consequences from an animal health standpoint.

Table 1 :

 1 Immunohistochemistry scoring. Means by species and by sex.

		Megaptera	Stenella	Stenella
		novaeangliae	attenuata	longirostris
		(n= 15)	(n= 7)	(n= 24)
	Female (n= 23)	0,25	1,5	1,11
	Male (n= 21)	0	0	1,4
	Not-determined (n=2)		0
	Total	0,2	0,42	1,29

ACKNOWLEDGMENTS

We acknowledge John J. Stegeman (Woods Hole Oceanographic Institution) and Harry V.

Gelboin (National Cancer Institute) who provided the monoclonal antibody against scup cytochrome P4501A. The field crew in Mayotte is thanked for its collaboration, especially Didier Fray, the personnel of the Brigade Nature, the Direction de l'Agriculture et de la Forêt and the Office National de la Chasse et de la Faune Sauvage. This research has been founded by the Collectivité Départementale de Mayotte and the French Ministry of Environment.

Biopsy sampling was conducted under scientific permit #78/DAF/2004 (September 10, 2004) and permit #032/DAF/SEF/2008 (May 16, 2008). LITERATURE CITED Aguilar, A., Borrell, A., 1994. Abnormally high polychlorinated biphenyl levels in striped dolphins (Stenella coeruleoalba) affected by the 1990-1992 Mediterranean epizootic. The Science of the Total Environment 154, 237-247. Aguilar, A., Borrell, A., Reijnders, P.J.H., 2002. Geographical and temporal variation in levels of organochlorine contaminants in marine mammals. Marine Environmental Research 53, 425-452.


Ana Arribas-Gil

Catherine Matias
email: catherine.matias@genopole.cnrs.fr

A context dependent pair hidden Markov model for statistical alignment

Keywords: Comparative genomics, Contextual alignment, DNA sequence alignment, EM algorithm, Insertion deletion model, Pair hidden Markov model, Probabilistic alignment, Sequence evolution, Statistical alignment, Stochastic expectation maximization algorithm

This article proposes a novel approach to statistical alignment of nucleotide sequences by introducing a context dependent structure on the substitution process in the underlying evolutionary model. We propose to estimate alignments and context dependent mutation rates relying on the observation of two homologous sequences. The procedure is based on a generalized pair-hidden Markov structure, where conditional on the alignment path, the nucleotide sequences follow a Markov distribution. We use a stochastic approximation expectation maximization (saem) algorithm to give accurate estimators of parameters and alignments. We provide results both on simulated data and vertebrate genomes, which are known to have a high mutation rate from CG dinucleotide. In particular, we establish that the method improves the accuracy of the alignment of a human pseudogene and its functional gene.

Introduction

Alignment of DNA sequences is the core process in comparative genomics. An alignment is a mapping of the nucleotides in one sequence onto the nucleotides in the other sequence which captures two important evolutionary processes: the substitution (of a nucleotide, by another one) and the insertion or deletion (indel, of one or several nucleotides). Thus, in this mapping, two matched nucleotides are supposed to derive from a common ancestor (they are homologous) and the mapping also allows for gaps which may be introduced into one or the other sequence. Several powerful alignment algorithms have been developed to align two or more sequences, a vast majority of them relying on dynamic programming methods that optimize an alignment score function. While score-based alignment procedures appear to be very rapid (an important point regarding to the current sizes of the databases), they have at least two major flaws: i) The choice of the score parameters requires some biological expertise and is non objective. This is particularly problematic since the resulting alignments are not robust to this choice [START_REF] Dewey | Parametric alignment of Drosophila genomes[END_REF]. ii) The highest scoring alignment might not be the most relevant from a biological point of view. One should then examine the k-highest scoring alignments (with the additional issue of selecting a value for k) but this approach often proves impractical because of the sheer number of suboptimal alignments [START_REF] Waterman | A new algorithm for best subsequence alignments with application to tRNA-rRNA comparisons[END_REF]. We refer to Chapters 2 and 4 in [START_REF] Durbin | Biological sequence analysis: probabilistic models of proteins and nucleic acids[END_REF] for an introduction to sequence alignment.

Contrarily to score-based alignment procedures, statistical (or probabilistic) alignment methods rely on an explicit modeling of the two evolutionary processes at the core of an alignment, namely the substitution and the indel processes. While there is a general agreement that substitution events are satisfactorily modeled relying on continuous time Markov processes on the nucleotide state space, very few indel models have been proposed in the literature and none of them has become a standard. The first evolutionary model for biological sequences dealing with indels was formulated by [START_REF] Thorne | An evolutionary model for maximum likelihood alignment of DNA sequences[END_REF], hereafter the TKF91 model. It provides a basis for performing alignment within a statistical framework: the optimal alignment is obtained as the evolutionary path that maximizes the likelihood of observing the sequences as being derived from a common ancestor (under this model). This has at least two major advantages over the score-based method. First, in this context, the evolutionary parameters (which play the role of an underlying scoring function) are not chosen by the user but rather estimated from the data, in a maximum likelihood framework. Second, a posterior probability is obtained on the set of alignments: the analysis does not concentrate on the highest probable path and rather provides sample alignments from this posterior distribution. Moreover, a confidence measure may be assigned to each position in the alignment. This is particularly relevant as alignment uncertainty is known to be a major flaw in comparative genomics [START_REF] Wong | Alignment uncertainty and genomic analysis[END_REF].

Statistical alignment is performed under the convenient pair-hidden Markov model, hereafter pair-HMM (see [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF][START_REF] Durbin | Biological sequence analysis: probabilistic models of proteins and nucleic acids[END_REF], enabling the use of generalizations of the expectation-maximization (em) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. Some drawbacks of the preliminary TKF91 proposal have first been improved by the same authors in what is called the TKF92 version of the model [START_REF] Thorne | Inching toward reality: an improved likelihood model of sequence evolution[END_REF]. Then, the original models have been later refined in many ways, as for instance in [START_REF] Arribas-Gil | Statistical alignment with a sequence evolution model allowing rate heterogeneity along the sequence[END_REF]; [START_REF] Knudsen | Sequence alignments and pair hidden Markov models using evolutionary history[END_REF]; [START_REF] Metzler | Statistical alignment based on fragment insertion and deletion models[END_REF]; [START_REF] Miklós | Algorithm for statistical alignment of two sequences derived from a poisson sequence length distribution[END_REF]; [START_REF] Miklós | A 'Long Indel' model for evolutionary sequence alignment[END_REF]; [START_REF] Miklós | An improved model for statistical alignment[END_REF]. We underline that an important characteristic of all those models is that the indel and the substitution processes are handled independently. An introduction to statistical alignment may be found in [START_REF] Lunter | Statistical alignment: Recent progress, new applications, and challenges[END_REF].

It is important to note that the original TKF91 model and all the later developments are primarily focused on modeling the indel process. As for the underlying substitution model, pair-HMMs procedures mostly rely on the simplest one: the one-parameter Jukes-Cantor model [START_REF] Jukes | Evolution of protein molecules[END_REF]. Surprisingly, while models of substitution processes have made important progress towards a more accurate description of reality in the past years (see for instance Felsenstein, 2003, Chapter 13), these developments have had almost no impact in the statistical alignment literature. One of the reasons for that may be the existence of important biases found in sequence alignment and resulting from incorrectly inferred indels. For instance, [START_REF] Lunter | Uncertainty in homology inferences: Assessing and improving genomic sequence alignment[END_REF] exhibit three major classes of such biases: (1) gap wander, resulting from the incorrect placement of gaps due to spurious non homologous similarity; (2) gap attraction, resulting in the joining of two closely positioned gaps into one larger gap; and (3) gap annihilation, resulting in the deletion of two indels of equal size for a typically more favorable representation as substitutions (see also [START_REF] Holmes | Dynamic programming alignment accuracy[END_REF]. This might explain a stronger focus on accurate estimation of indels positions. Another reason could be that since evolutionary parameters are to be estimated in statistical alignments procedures, their number should remain quite low.

On the contrary, score-based alignment methods offer some variety in the form of the match scoring functions, which is the counterpart of the substitution process parameters when relying on statistical alignment. Note that in this context, using different types of indels scores is hardly possible, since the dynamic programming approach limits those scores to affine (in the gap length) penalties. Nonetheless, the match scores may take different forms. One of the most recent improvements in scorebased alignment procedures concerns the development of context-based scoring schemes as in [START_REF] Gambin | Alignment with context dependent scoring function[END_REF]; [START_REF] Gambin | CTX-BLAST: context sensitive version of protein BLAST[END_REF] (see also [START_REF] Huang | A context dependent method for comparing sequences[END_REF], as a much earlier attempt). In these works, the cost of a substitution depends on the surrounding nucleotides. As a consequence, the score of an alignment depends on the order of editing operations. Another drawback of this approach lies in the issue of computing a p-value associated to the alignment significance. Indeed, in order to assess the statistical significance of an alignment, one requires the knowledge of the tail distribution of the score statistic, under the null hypothesis of non homology between the sequences. While this distribution is not even fully characterized in the simpler case of site independent scoring schemes allowing gaps, it is currently out of reach when introducing dependencies in the substitution cost.

Nonetheless, allowing for context dependence in the modeling of the evolutionary processes is important as it captures local sequence dependent fluctuations in substitution rates. It is well-known for instance that, at least in vertebrates, there is a higher level of methylated cytosine in CG dinucleotides (denoted CpG). This methylated cytosine mutates abnormally frequently to thymine which results in higher substitution rates form the pair CpG [START_REF] Bulmer | Neighboring base effects on substitution rates in pseudogenes[END_REF]. At the mutational level, context dependent substitution rates can have a significant impact on prediction accuracy [START_REF] Siepel | Computational identification of evolutionarily conserved exons[END_REF].

Note that allowing for context dependent indel rates is also a challenging issue. For instance, microsatellite expansion and contraction is known to be a context dependent process (see for instance some recent results in [START_REF] Varela | Evidence for nonindependent evolution of adjacent microsatellites in the human genome[END_REF] and suggests the use of context dependent indel rates. Such an attempt has been very recently made by Hickey and Blanchette (2011a,b), relying on tree adjoining grammars. However, introducing locally dependent indel rates in an evolutionary process implies that the substitution and the indel processes may not anymore be handled separately. This independence between the two processes is a key assumption in limiting the computational cost of pair-HMMs and removing this hypothesis is beyond the scope of the present work. Thus, we limit our approach to modeling context dependent substitution rates.

In this work, we formulate a context dependent pair hidden Markov model, so as to perform statistical alignment of two sequences allowing for context dependent substitution rates. We then propose a method to perform parameter maximum likelihood estimation as well as posterior sampling within this model. Our procedure relies on stochastic versions [START_REF] Celeux | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF][START_REF] Diebolt | Asymptotic properties of a stochastic EM algorithm for estimating mixing proportions[END_REF] of the em algorithm. A similar method has already been successfully used in [START_REF] Arribas-Gil | Statistical alignment with a sequence evolution model allowing rate heterogeneity along the sequence[END_REF] in a classical pair-HMM (with no context dependent substitution rates) in which case it gives results very similar to Gibbs sampling. Note that Holmes (2005) also suggests the use of a stochastic version of em for estimating indel rates in the multiple alignment framework, assuming the alignment is unknown. In the same spirit, [START_REF] Hobolth | A Markov chain Monte Carlo Expectation Maximization algorithm for statistical analysis of DNA sequence evolution with neighbor-dependent substitution rates[END_REF] proposes an MCMC-EM algorithm for the estimation of neighbor-dependent substitution rates from a multiple alignment and a phylogenetic tree. The advantages of stochastic versions of em algorithm is that they allow to perform maximum likelihood estimation in a hidden variables framework where the maximization is not explicit, the computation of the likelihood is computationally very expensive, and where the use of the classical em algorithm is not feasible. To our knowledge, this work is, with Hickey and Blanchette (2011a)'s approach, one of the first attempts to perform a contextual statistical alignment. This article is organized as follows. In Section 2, we first describe a context dependent pair hidden Markov model. In this framework, parameter maximum likelihood estimation combined with hidden states posterior sampling provides (a posterior probability on) contextual alignments of two sequences which is data-driven in the sense that it does not require any parameter choice. Section 3 briefly discusses asymptotic consistency results that can be obtained on the parameter estimation procedure. Then, Section 4 describes two versions of the em algorithm applied in this framework of contextual statistical alignment: the stochastic expectation maximization (sem) and the stochastic approximation expectation maximization (saem) [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] algorithms. The performance of these procedures are illustrated on synthetic data in Section 5 and a real data set is handled in Section 6.

Model

A pair-hidden Markov model is described by the distribution of a non-observed increasing path through the two-dimensional integer lattice N × N and given this hidden path, the conditional distribution of a pair of observed sequences X1:n := X1, . . . , Xn ∈ A n and Y1:m := Y1, . . . , Ym ∈ A m , with values on a finite set A (the nucleotides alphabet). The path reflects the insertion-deletion process acting on the sequences and corresponds to the bare alignment, namely the alignment without specification of the nucleotides. It is the sum of random moves from the set E = {(1, 1), (1, 0), (0, 1)} also denoted {M, IX, IY }, where M stands for 'match', IX stands for 'insert in X' and similarly IY stands for 'insert in Y '. This path represents an empty alignment of the sequences, as the move M represents an homologous site between the two sequences, the move IX represents either an insertion in the first sequence (X1:n) or a deletion in the second (Y1:m) and similarly for the move IY (see Figure 1).

A C A - T T G G - G .
We assume that the hidden path follows a Markov distribution without specifying from which insertion-deletion process it comes from. This could be for instance the TKF91 insertion-deletion process [START_REF] Thorne | An evolutionary model for maximum likelihood alignment of DNA sequences[END_REF] or any of its later variants [START_REF] Arribas-Gil | Statistical alignment with a sequence evolution model allowing rate heterogeneity along the sequence[END_REF][START_REF] Knudsen | Sequence alignments and pair hidden Markov models using evolutionary history[END_REF][START_REF] Metzler | Statistical alignment based on fragment insertion and deletion models[END_REF][START_REF] Miklós | Algorithm for statistical alignment of two sequences derived from a poisson sequence length distribution[END_REF][START_REF] Thorne | Inching toward reality: an improved likelihood model of sequence evolution[END_REF]. Given the empty alignment, we then describe the conditional distribution of the sequence letters. In a 'match' position, this distribution reflects the substitution process while in 'insert in X' or 'insert in Y' positions, this distribution reflects both the substitution process and the distribution of inserted segments.

More precisely, let {εt} t 0 be an ergodic, homogeneous and stationary Markov chain with state space E = {(1, 1), (1, 0), (0, 1)} and transition probability π (a stochastic matrix of size 3 × 3). For notational convenience, the state space E is sometimes equivalently denoted by {M, IX, IY }. In this context, we might use obvious notations as for instance πMM for the conditional probability P(εt = (1, 1)|εt-1 = (1, 1)).

Let Zt = (Nt, Mt) = t s=1 εs be the path induced by this Markov chain through the lattice N × N. In the following, a path through the lattice N × N is always assumed to be increasing, namely of the form s es for some sequence e ∈ E N . We denote by Pπ the corresponding stationary probability and by En,m the set of paths in N × N starting at (0, 0) and ending at (n, m).

For any e ∈ En,m, we also let |e| denote the length of e, which satisfies n ∨ m |e| n + m.

The hidden process generates the observed sequences in the following way. Given the hidden path {εt} t 0 , the sequences {X1:n, Y1:m} are generated according to an order one Markov chain: for any e ∈ En,m,

P(X1:n, Y1:m|ε = e) = |e| s=1 P(XN s , YM s |es, es-1, XN s-1 , YM s-1 ).
Note that the indexes Ns, Ms of the observations generated at time s are random. This is in sharp contrast with the classical hidden Markov model (HMM) and is specific to pair-HMMs (see Arribas-Gil et al., 2006, for more details).

In this work, we restrict our attention to order-one Markov conditional distributions. Note that the following developments could be generalized to higher order Markov chains at the cost of increased computational burden.

As we want to model the dependency in the substitution process only, we further constrain the model so that the dependency occurs only in successive homologous (match) positions. We thus use the following parametrization

P(XN s , YM s |es, es-1, XN s-1 , YM s-1 ) =        h(XN s , YM s |XN s-1 , YM s-1 ) if es = es-1 = (1, 1), h(XN s , YM s ) if es = (1, 1) = es-1, f (XN s ) if es = (1, 0), g(YM s ) if es = (0, 1),
where f, g are probability measures (p.m.) on A; function h is a p.m. on A 2 and h may be viewed as a stochastic matrix on A 2 . Obvious necessary conditions for the parameters to be identifiable in this model are the following

   i) ∃a, b ∈ A such that h(a, b) = f (a)g(b), ii) ∃a, b, c, d ∈ A such that h(a, b|c, d) = f (a)g(b), iii) ∃a, b, c, d ∈ A such that h(a, b|c, d) = h(a, b).
(1)

The whole parameter set is then given by

Θ = {θ = (π, f, g, h, h); θ satisfies (1)}.
Statistical alignment of two sequences consists in maximizing with respect to θ the following criterion wn,m(θ) [START_REF] Durbin | Biological sequence analysis: probabilistic models of proteins and nucleic acids[END_REF]. This criterion plays the role of a log-likelihood in the pair-HMM. (For a discussion on the quantities playing the role of likelihoods in pair-HMMs, see [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF]. For any integers n, m 1 and any observations X1:n and Y1:m, let wn,m(θ) := log Q θ (X1:n, Y1:m) := log P θ (∃s 1, Zs = (n, m); X1:n, Y1:m).

(

The criterion wn,m is more explicitly defined as

wn,m(θ) = log n+m s=n∨m e∈En,m,|e|=s Pπ(ε1:s = e1:s) × s k=1 f (Xn k ) 1{e k =(1,0)} g(Ym k ) 1{e k =(0,1)} × h(Xn k , Ym k ) 1{e k =(1,1),e k-1 =(1,1)} h(Xn k , Ym k |Xn k-1 , Ym k-1 ) 1{e k =e k-1 =(1,1)} ,
where (n k , m k ) = k t=1 et. Now, we define the parameter estimator as θn,m = arg max θ∈Θ wn,m(θ).

In the classical pair-HMM (namely without accounting for dependencies in the substitution process), consistency results for the above estimator θn,m were obtained in [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF] and are extended here to the context dependent case in Section 3. Then in Section 4, we develop an saem algorithm to compute θn,m and provide an a posteriori probability distribution over the set of alignments. Indeed, a main issue is to obtain an accurate alignment of the sequences. Score-based alignment methods heavily rely on the Viterbi algorithm for this purpose (see [START_REF] Durbin | Biological sequence analysis: probabilistic models of proteins and nucleic acids[END_REF], which provides the most probable a posteriori alignment. However, optimal alignments often look different from typical ones, and providing a unique alignment as the result of the estimation procedure may not be very informative. For this reason, it may be more interesting to provide the probability distribution of hidden states over each pair of nucleotides from the observed sequences (for the parameter value θn,m). This can be done under the pair-HMM framework and gives us a reliability measure of any alignment of the two sequences (see [START_REF] Arribas-Gil | Statistical alignment with a sequence evolution model allowing rate heterogeneity along the sequence[END_REF].

Consistency results for parameter estimation

In this part, we give generalizations of results first obtained in [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF], to the context dependent pair-HMM.

As the proofs follow the same lines as in this reference, we shall omit them. We first introduce some notations and formulate some assumptions that will be needed. For any δ > 0, we let

Θ δ = {θ ∈ Θ; ∀k, θ k δ} and Θ0 = {θ ∈ Θ; ∀k, θ k > 0}.
The true parameter value will be denoted by θ0 and is assumed to belong to the set Θ0. Probabilities and expectations under this parameter value are denoted by P0 and E0, respectively. We now introduce a notation for the marginal distributions of We then let Θmarg be a subset of parameters in Θ0 satisfying some assumptions on the marginals of h and h.

Θmarg = θ ∈ Θ0; ∀b, c ∈ A, hX = f, hY = g, hX (•|b, c) = f (•), hY (•|b, c) = g(•) .
Note that for instance, the parametrization used in Section 5 satisfies θ ∈ Θmarg. Finally, we introduce the process wt defined in a similar way as in (2) by

wt(θ) := log Q θ (X1:N t , Y1:M t ).
As previously noted, the length t of an alignment of sequences with respective sizes n and m satisfies n ∨ m t n + m. Thus, asymptotic results for t → +∞ will imply equivalent ones for n, m → +∞. In other words, consistency results obtained when t → +∞ are valid for long enough observed sequences, even if one does not know the length t of the true underlying alignment. We now state the main theorem.

Theorem 1: For any θ ∈ Θ0, we have i) Renormalized log-likelihood t -1 wt(θ) converges P0-almost surely and in L1, as t tends to infinity, to

w(θ) = lim t→∞ 1 t E0 (log Q θ (X1:N t , Y1:M t )) = sup t 1 t E0 (log Q θ (X1:N t , Y1:M t )) .
ii) Moreover, w(θ0) w(θ) and strict inequality is valid as long as at least one of these conditions is satisfied

A. θ ∈ Θ0 and ∀λ > 0, E θ (ε1) = λE0(ε1), B. θ0, θ ∈ Θmarg and either f = f0 or g = g0.
iii) The family of functions {t -1 wt(θ)} t 1 is uniformly equicontinuous on the set Θ δ , for any δ > 0.

The three statements in Theorem 1 are the key ingredients to establish the consistency of θn,m, following Wald's classical proof [START_REF] Wald | Note on the consistency of the maximum likelihood estimate[END_REF] of maximum likelihood estimators consistency (see also van der Vaart, 1998, Chapter 5). The most subtle part is point ii), that states some sufficient conditions under which the maximum of the limiting function w is attained only at the true parameter value θ0. For instance, condition A means that the main directions of the paths {Zt} t 1 generated under the two parameter values θ and θ0 are different. As a consequence, these two parameters may be easily distinguished from each other relying on the distributions P θ and P0. Note that these two different cases (condition A or B satisfied) are certainly not exhaustive and we refer to [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF] for simulation results indicating that point ii) might be true under other scenarios. Note in particular that for evolutionary models assuming that insertions and deletions happen at the same rate, E θ (ǫ1) is proportional to (1, 1) for any parameter value θ. In this case, condition A) is never satisfied and we are not able to prove that h = h0 or h = h0 implies strict inequality w(θ0) > w(θ). This is due to the dependency structure on the sequences, that make it difficult to link the difference w(θ0) -w(θ) with a Kullback-Leibler divergence between densities h and h0 or between h and h0 (see the proof of Theorem 2 in [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF], for more details).

In general, one is not interested in the full (context) pair-HMM parameters but rather in a sub-parametrization induced by the choice of a specific evolutionary model (see for instance the sub-parametrization proposed in Section 5). Let β → θ(β) be a continuous parametrization from some set B to Θ. For any δ > 0, let B δ = θ -1 (Θ δ ). We assume that β0 := θ -1 (θ0) belongs to B δ for some δ > 0. We then let βn,m := arg max

β∈B δ wn,m(θ(β)).
Besides from giving results on the frequentist estimator βn,m, we consider Bayesian estimates of the parameters. Let ν be a prior probability measure on the set B δ . Markov chain Monte Carlo (MCMC) algorithms approximate the random distribution νn,m, interpreted as the posterior measure given observations X1:n and Y1:m and defined by

νn,m(dβ) = Q θ(β) (X1:n, Y1:m)ν(dβ) B δ Q θ(β ′ ) (X1:n, Y1:m)ν(dβ ′ ) .
With these definitions at hand, we can now state the following corollary.

Corollary 1: If the maximizer of β → w(θ(β)) over the set B δ is reduced to the singleton {β0}, then we have the following results i) βn,m converges P0-almost-surely to β0, ii) If ν puts some weight on a neighborhood of β0, then the sequence of posterior measures νn,m converges in distribution, P0-almost surely, to the Dirac mass at β0.

Algorithms

The approach we use here to maximize criterion (2) and provide a posterior distribution on the set of alignments, given two observed sequences relies on a stochastic version [START_REF] Celeux | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF][START_REF] Diebolt | Asymptotic properties of a stochastic EM algorithm for estimating mixing proportions[END_REF] of em algorithm [START_REF] Baum | A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains[END_REF][START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. This approach has already been successfully used in [START_REF] Arribas-Gil | Statistical alignment with a sequence evolution model allowing rate heterogeneity along the sequence[END_REF] in a classical pair-HMM (meaning with no context dependent substitution rates) in which case it gives results very similar to Gibbs sampling. It allows to perform maximum likelihood estimation in a framework in which the maximization is not explicit, the computation of the likelihood is computationally very expensive, and in which the use of em algorithm is not feasible (see Appendix for details).

Forward and backward equations

We describe below the forward and backward probabilities used in the later procedures. The formulas are easily obtained by generalizing those from the classical pair-HMM [START_REF] Arribas-Gil | Estimation dans des modèles à variables cachées : alignement de séquences biologiques et modèles d'évolution[END_REF][START_REF] Durbin | Biological sequence analysis: probabilistic models of proteins and nucleic acids[END_REF].

The forward equations are obtained as follows. For any value u ∈ E , let E u n,m = {e ∈ En,m, e |e| = u} be the set of paths ending at (n, m) with last step being u. Then the forward probabilities are defined for any i n and j m, as α u (i, j) = P θ (∃s 1, Zs = (i, j), εs = u, X1:i, Y1:j) = e∈E u n,m P θ (ε = e, X1:i, Y1:j).

These forward probabilities are computed recursively in the following way. We first set the initial values

∀u ∈ E , α u (0, -1) = α u (-1, 0) = 0 ∀u ∈ {IX , IY }, α u (0, 0) = 0 and α M (0, 0) = 1.
Then for i = 0, . . . , n and j = 0, . . . , m except (i, j) = (0, 0), we recursively compute α M (i, j) = P θ (∃s 1, εs = (1, 1), Zs = (i, j), X1:i, Y1:j)

= u∈E P θ (∃s 1, εs = (1, 1), Zs = (i, j), εs-1 = u, X1:i, Y1:j) = h(Xi, Yj|Xi-1, Yj-1)πMM α M (i -1, j -1) + h(Xi, Yj)πI X M α I X (i -1, j -1) + h(Xi, Yj)πI Y M α I Y (i -1, j -1), (3)
and in the same way

α I X (i, j) = f (Xi) u∈E πuI X α u (i -1, j), α I Y (i, j) = g(Yj) u∈E πuI Y α u (i, j -1).
Note that the log-likelihood of the observed sequences is then simply obtained from the forward probabilities by taking wn,m(θ) = log P θ (∃s 1, Zs = (n, m), X1:n, Y1:m) = log(α

I X (n, m) + α I Y (n, m) + α M (n, m)).
Note that the forward equations could thus be used to compute the log-likelihood of two observed sequences and a numerical optimization of this quantity could give maximum likelihood estimates. However, such a strategy fails short as soon as the dimension of the parameter space is not small. Indeed, the computational cost of calculating the log-likelihood for a given parameter value is O(nm), and thus the computational cost of its maximization over a grid is O(p k nm), where k is the number of parameters and p the number of values considered for each parameter.

We now describe the backward probabilities and their recursive computation. For any value u ∈ E , the backward probabilities are defined as β u (i, j) = P θ (Xi+1:n, Yj+1:m|∃s 1, Zs = (i, j), εs = u, Xi, Yj).

We initialize these values by taking, for any u ∈ E , β u (n, m) = 1 and β u (n, m + 1) = β u (n + 1, m) = 0, and recursively compute, for i = n, . . . , 1 and j = m, . . . , 1, except (n, m), the quantities

β u (i, j) = πuM β M (i + 1, j + 1){h(Xi+1, Yj+1)1 u∈{I X ,I Y } + h(Xi+1, Yj+1|Xi, Yj)1u=M } + πuI X β I X (i + 1, j)f (Xi+1) + πuI Y β I Y (i, j + 1)g(Yj+1). (4)
Note that contrarily to the case of HMM, the forward-backward probabilities do not give rise to the classical em strategy, as we do not obtain from these equations the conditional expectation of the complete log-likelihood, given the two observed sequences (see Appendix for further details).

sem and saem algorithms

The forward-backward probabilities may be used to simulate the hidden path, conditional on the observed sequences. Thus, em algorithm may be replaced by sem [START_REF] Celeux | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF] or saem [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] procedures. Indeed, let us explain the basic idea of these stochastic versions of em algorithm. Denoting by Lnm the random value s 1 such that Zs = (n, m) (the first and only hitting time for the point (n, m), which is not necessarily finite), the complete log-likelihood is log P θ (X1:n, Y1:m, Lnm, ε1:L nm ) (see Appendix for details) and for a current parameter value θ ′ , its conditional expectation writes Q θ ′ (θ) := E θ ′ (log P θ (X1:n, Y1:m, Lnm, ε1:L nm )|X1:n, Y1:m). Now, the idea of these stochastic approximations is to replace the computation of Q θ ′ (θ) by an approximation obtained from the simulation (under parameter value θ ′ ) of a number of hidden paths {ε1:s} in En,m, with possibly different lengths. These paths satisfy the property that for each obtained value of s, the paths with length s have the same occurrence probability as P θ ′ (ε|X1:n, Y1:m, Lnm = s). Note that the length of each simulated path is random. Then the maximization with respect to θ of the complete log-likelihood is done on the basis of the simulated alignments. As a consequence, this maximization is performed via a simple counting of specific events and does not require numerical procedures. In particular, iteration r of sem algorithm writes -Simulation step: generate one hidden path e r of some random length s with the same distribution as P θ (r-1) (ε1:s|X1:n, Y1:m, Lnm = s).

-Maximization step: θ (r) = arg max θ∈Θ log P θ (X1:n, Y1:m, ε = e r ).

And iteration r of saem algorithm writes -Simulation step: generate m(r) hidden paths e r (j), j = 1, . . . , m(r) with resulting random lengths sr,j , each one having the same distribution as P θ (r-1) (ε1:s r,j |X1:n, Y1:m, Lnm = sr,j).

-Stochastic approximation step: update

Qr(θ) = Qr-1(θ) + γr 1 m(r) m(r) j=1 log P θ (X1:n, Y1:m, ε = e r (j)) -Qr-1(θ) ,
where {γr} r 1 is a decreasing sequence of positive step size and Q0(θ) = 1/m(0) m(0) j=1 log P θ (X1:n, Y1:m, ε = e 0 (j)) . -Maximization step: θ (r) = arg max θ∈Θ Qr(θ).

The simulation step is common to both algorithms and consists in drawing paths of some unknown length s with same distribution as P θ (r-1) (ε1:s|X1:n, Y1:m, Lnm = s). This is performed via the backwards sampling based on the forward probabilities as we now explain. We may write P θ (r-1) (ε1:s|X1:n, Y1:m, Lnm = s) = P θ (r-1) (εs|X1:n, Y1:m, Lnm = s)

s-1 k=1 P θ (r-1) (ε k |ε k+1:s , X1:n, Y1:m, Lnm = s) = P θ (r-1) (εs|X1:n, Y1:m, Zs = (n, m)) s-1 k=1 P θ (r-1) (ε k |ε k+1 , Z k+1 , X1:N k+1 , Y1:M k+1 ). (5)
The last equality comes from the following: Conditional on ε k+1 , the random variable ε k does not anymore depend on the observations after time k + 2. If we moreover condition on the variables ε k+1:s and Lnm = s, then we know which point (i, j) on the lattice corresponds to Z k+1 = (n, m) -s t=k+2 εt = (i, j). The values of the observed sequence after time k + 2 correspond exactly to Xi+1:n, Yj+1:m. Then conditional on {ε k+1:s , X1:n, Y1:m, Lnm = s}, the variable ε k only depends on ε k+1 (according to the Markov property) and on the values Z k+1 , X1:N k+1 , Y1:M k+1 .

So, given ε k+1 , the position Z k+1 = (i, j) on the lattice and the observed values X1:i, Y1:j , we want to sample ε k from P θ (r-1) (ε k |ε k+1 , X1:i, Y1:j, Z k+1 = (i, j)), for k = s, s -1, . . . , 1. Up to a renormalizing constant, each step corresponds to sampling from P θ (r-1) (ε k , ε k+1 |X1:i, Y1:j, Z k+1 = (i, j)) or equivalently from πε k ,ε k+1 α ε k ((i, j) -ε k+1 ). Then the backwards sampling writes -e r s is sampled to be v with probability

α v (n, m) α M (n, m) + α I X (n, m) + α I Y (n, m) ,
-For k < s and Z k+1 = (i, j), i > 1, j > 1, if e r k+1 = M , then e r k is sampled to be M with probability α M (i -1, j -1)πMM h(Xi, Yj|Xi-1, Yj-1) α M (i, j) ,

and to be u = IX or IY with probability α u (i -1, j -1)πuM h(Xi, Yj) α M (i, j) .

Note that these latter probabilities are correctly renormalized according to Equation (3). -For k < s and Z k+1 = (i, j), i > 1, j > 1, if e r k+1 = IX then e r k is sampled to be u with probability

α u (i -1, j)πuI X f (Xi) α I X (i, j) , (6)
-For k < s and Z k+1 = (i, j), i > 1, j > 1, if e r k+1 = IY then e r k is sampled to be u with probability

α u (i, j -1)πuI Y g(Yj) α I Y (i, j) , (7)
-For k < s and Z k+1 = (1, j), j > 1, if e r k+1 = M or IX then all the values e r l , 1 l k are fixed to IY , otherwise e r k is sampled to be u with probability given by ( 7).

-For k < s and Z k+1 = (i, 1), i > 1, if e r k+1 = M or IY then all the values e r l , 1 l k are fixed to IX , otherwise e r k is sampled to be u with probability given by ( 6).

According to Equation ( 5), we thus have sampled paths e r satisfying the following property.

Proposition 1: The paths e ∈ En,m sampled from the above scheme with parameter value θ occur with probability P θ (ε = e|X1:n, Y1:m, Lnm = |e|).

We now consider the complete log-likelihood obtained from the observed sequences and one simulated path. We have a,b,c,d) log h(a, b|c, d),

log P θ (X1:n, Y1:m, ε = e r ) = u∈E 1 e r 1 =u log πu + |e r | k=2 u,v∈E 2 1 e r k-1 =u,e r k =v log πuv + |e r | k=1 a∈A [1 e r k =I X ,Xn k =a log f (a) + 1 e r k =I Y ,Ym k =a log g(a)] + |e r | k=1 a,b∈A 2 1 e r k =M,Xn k =a,Ym k =b log h(a, b) + |e r | k=2 a,b,c,d∈A 4 1 e r k =e r k-1 =M,(Xn k ,Ym k ,Xn k-1 ,Ym k-1 )=(
where we recall that (n k , m k ) = k i=1 e r i (for simplicity, we omit the subscript r). Then in the maximization step of sem, the components of θ are updated to the values: ∀u, v ∈ E , a, b, c, d ∈ A,

π (r) u = Ner (u) v∈E Ner (v) , π (r) u,v = Ner (u, v) w∈E Ner (u, w) , f (r) (a) = Ner (a|IX) a ′ ∈A Ner (a ′ |IX , ) , g (r) (a) = Ner (a|IY ) a ′ ∈A Ner (a ′ |IY , ) , h (r) (a, b) = Ner (a, b|M ) a ′ ,b ′ ∈A Ner (a ′ , b ′ |M ) , h(r) (a, b|c, d) = Ner (a, b|M M ; c, d) a ′ ,b ′ ∈A Ner (a ′ , b ′ |M M ; c, d) , (8)
where we use obvious notations for the counts Ne in sequence e. Namely for any u, v ∈ E , and any a, b, c, d ∈ A, we let

Ne(u) = k 1{e k = u}, Ne(u, v) = 1{e k-1 = u, e k = v}
; c, d) = 1{e k = e k-1 = M, (XN k , YM k , XN k-1 , YM k-1 ) = (a, b, c, d)}.
In the maximization step of saem we have to take into account the m(r) hidden paths and the stochastic approximation of Qr(θ). Then the values of θ are updated as in ( 8) by replacing the counts Ner (•) by their counterparts Ñr(•), where

Ñr = Ñr-1 + γr 1 m(r) m(r)
j=1 N e r (j) -Ñr-1 . From a theoretical point of view, the convergence properties of the two algorithms are different. In saem, the sequence {θ r } r 1 converges, under general conditions, to a local maximum of the log-likelihood. It is important to note that the stochastic perturbation introduced with respect to the original em algorithm allows avoiding saddle points, which are possible attractive stationary points of the sequence generated by em [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF]. In sem, the sequence {θ r } r 1 does not converge pointwise, but it is an homogeneous Markov chain which is ergodic under appropriate conditions, see for instance [START_REF] Diebolt | A stochastic EM algorithm for approximating the maximum likelihood estimate[END_REF].

In this context, maximum likelihood estimators are sensitive to overfitting if there is insufficient data. Indeed, if an event has never been observed in the sequences, the corresponding estimator is not well defined, as both the numerator and denominator in (8) are zero. Pseudo-counts (Durbin et al., 1998, Chapter 3) or other smoothing techniques such as the one proposed in [START_REF] Kneser | Improved backing-off for m-gram language modeling[END_REF] may be used to solve this problem. However, when using saem algorithm to perform maximum likelihood estimation, this problem is minimized and no smoothing strategies are required. Indeed, since at iteration r we

Table 1

Parameter values in the two data sets.

Data set 1 α = 0.4 β = 0.2 γ = 0.06 λ = 0.04 Data set 2 α = 0.5 β = 0.15 γ = 0.05 λ = 0.02 generate m(r) hidden paths, the counts Ñr (computed as an average over those alignments) are rarely equal to 0, even during the first iterations when γr is typically set to 1. However, if we want to completely avoid this possibility, we may replace event counts by pseudo-counts in the case of unobserved events, during these first iterations. Note also that in practice we do not generate whole alignments at each iteration step, but only sub-alignments of a certain given length within the current alignment. This strategy combines an MCMC procedure within saem and is known to preserve convergence properties of saem algorithm (see [START_REF] Arribas-Gil | Statistical alignment with a sequence evolution model allowing rate heterogeneity along the sequence[END_REF], for further details). In terms of event counts, this strategy implies that a current low estimate of an event occurrence probability will have some impact on next iteration estimate. Nonetheless, at the last iterations of the algorithm, when typically γr < 1, current counts are averaged with counts from previous iterations, which provides a smoothing procedure by itself.

Simulations

For these simulations, we consider nucleotide sequences with alphabet A = {A, C, G, T } and a constrained parameter set, motivated by the choice of an underlying evolutionary model on these sequences. As for the generation of the hidden Markov chain, we consider the TKF91 indel model with equal value for insertion and deletion rate λ (see Arribas-Gil, 2007, for more details). The transition matrix of the hidden path is thus given by

π = 1 1 + λ     e -λ 1 -e -λ λ λe -λ 1 -e -λ λ 1 + λ - λ 1 -e -λ e -λ 1 -e -λ λ     ,
where the order of the states in the matrix is M , IX and IY . Then, we describe the conditional distribution of the sequences. For any x ∈ A, we set f (x) = g(x) = µx, where µ is going to be the stationary distribution of the substitution process, according to the parametrization of h and h considered below. We set the distribution h according to a simple substitution process with equal substitution rate γ and different nucleotides frequencies (modified Jukes Cantor). Let

h(x, y) = µx(1 -e -γ )µy if x = y µx[µx(1 -e -γ ) + e -γ ] if x = y.
The parameter h accounting for the context in the substitution process differs from h only in a C C match context. In this way, we want to take into account a possibly higher substitution rate from the pair CpG to CpA. More precisely, we let

P(XN s = x, YM s = y|es = es-1 = M, XN s-1 = x ′ , YM s-1 = y ′ ) = hC (x, y) if x ′ = y ′ = C, h (x, y) otherwise.
In a C C match context, we use a model initially introduced in two slightly different forms by [START_REF] Felsenstein | Distance methods for inferring phylogenies: A justification[END_REF] and [START_REF] Hasegawa | Dating of the human-ape splitting by a molecular clock of mitochondrial DNA[END_REF]. This model considers different substitution rates for transitions (i.e. substitutions within the chemical class of purines R = {A, G} or pyrimidines Y = {C, T }) and for transversions (substitutions modifying the chemical class) as well as different nucleotide frequencies µ. Denoting by x the other nucleotide in the same chemical class as x (namely x = x and either {x,

x} ∈ R or {x, x} ∈ Y), we get hC (x, y) =      µxe -(α+β) + µxe -β (1 -e -α ) µx µx+µx + µx(1 -e -β )µx if y = x, µxe -β (1 -e -α ) µx µx +µx + µx(1 -e -β )µx if y = x, µx(1 -e -β )µy otherwise.
We refer to [START_REF] Felsenstein | Inferring phylogenies[END_REF] for more details on substitution processes.

According to this parametrization we have simulated two sets of data with an alignment length of 2000 base pairs (bp) and parameter values described in Table 1. In these data sets, the substitution rate γ is larger than the indel rate λ as expected in biological sequences, and the transition and transversion substitution rates in a C C match context (α and β) are larger than the regular substitution rate γ. In both sets of data we have set the stationary distribution of the substitution process to µA = 0.225 µC = 0.275 µG = 0.275 µT = 0.225, in order to obtain a GC content (proportion of Gs or Cs) larger than 50%. For each one of the data sets we have conducted two estimation procedures. On the one hand, we have estimated the whole parameter (π, f, g, h, h) via the saem algorithm as explained in Section 4.2. On the other hand, we have combined the saem algorithm with numerical optimization to estimate only the evolutionary parameters λ, γ, α and β. The first procedure has the advantages of being robust against misspecification of the underlying evolutionary model and relying on a explicit maximization step based on the counts of events. The second procedure is computationally more expensive, however, it is more parsimonious and it provides a straightforward evolutionary interpretation of the parameter values. In both cases we have performed 150 iterations of the saem algorithm with the parameter γr of the stochastic approximation set to 1 for r = 1, . . . , 100 and to 1/(r -100) for r = 101, . . . , 150. The number of simulated hidden paths is m(r) = 5, r = 1, . . . , 20 and m(r) = 10, r = 21, . . . , 150. We used the same initial values of the parameters for every simulation run over the two simulated data sets. These initial values, estimates and standard deviations are given in Tables 2,3 (first procedure) and 4 (second procedure). In Figure 2 we present the distribution of estimates obtained with the second estimation approach for the two simulated data sets. We can observe that both procedures provide unbiased and very accurate estimates. In Figure 3 we show the convergence of the saem estimates on a single simulation run from the second data set. We use a logarithmic scale for the x-axis to highlight the very fast convergence of the parameters to their true values.

Application to real data

In this section we illustrate our method through an alignment of the human alpha-globin pseudogene (HBPA1) with its functional counterpart, the human alpha-globin gene (HBA1). This example is inspired by the work of [START_REF] Bulmer | Neighboring base effects on substitution rates in pseudogenes[END_REF], who studied the neighbouring base effects on substitution rates in a series of vertebrate pseudogenes, including the human alphaglobin, concluding that there is an increase in the frequency of substitutions from the dinucleotide CpG. We have extracted the sequences of the gene and the pseudogene, who are located in the human alpha-globin gene cluster on chromosome 16, from the UCSC Genome Browser database [START_REF] Fujita | The UCSC Genome Browser database: update 2011[END_REF]. We have considered the whole sequences (coding and non-coding sequence in the case of the functional gene) for the alignment, obtaining sequences lengths of 842 bp for HBA1 and 812 bp for HBPA1. Because of the presence of introns and exons in the HBA1 sequence, it is natural to consider a model allowing for two different substitution behaviors along the sequence. That is why we have used a model for which the state space of the hidden Markov chain is {M1, M2, IX, IY }, where Mi, i = 1, 2, stand for two different match states, with different associated substitution processes hi, hi, i = 1, 2. This kind of model may also handle fragments insertion and deletion in contrast to single nucleotide indels (see [START_REF] Arribas-Gil | Statistical alignment with a sequence evolution model allowing rate heterogeneity along the sequence[END_REF], for more details). As in the simulation studies of the previous section, h is the substitution matrix in a C C match context, and h is the substitution matrix in any other case. In this way, we want to take into account a possibly higher transition rate from nucleotide G occurring in the dinucleotide CpG. In order to avoid any possible misspecification of the underlying indel and substitution processes, we have decided not to parametrize any of those, conducting an estimation approach equivalent to the first procedure described in the previous section. The only difference here lies in the dimension of the model parameter, since this parameter is now (π, f, g, h1, h2, h1, h2), where π is a 4 × 4 transition probability matrix, and hi, hi, i = 1, 2 are four different 4 × 4 stochastic vectors, yielding a total number of 12 + 3 + 3 + 4 × 15 = 78 free parameters. We have run saem algorithm on the sequences, performing 500 iterations with γr set to 1 for r = 1, . . . , 400 and to 1/(r -400) for r = 401, . . . , 500. The number of simulated hidden paths is m(r) = 5, r = 1, . . . , 20 and m(r) = 10, r = 21, . . . , 500.

In Table 5, we present the estimated substitution probabilities from nucleotide G. As expected, we found an increase in the substitution occurrence for G in a C C match context. We present in Figure 4 the posterior probability distribution over the set of possible alignment columns, match, insertion and deletion (the two match states have been merged) for every pair of aligned positions in a consensus alignment obtained at the convergence stage (last iterations) of saem algorithm.

In order to compare our method to a traditional pair-HMM alignment, we have also run saem algorithm, with the same settings as before, on a simpler model without substitution context-dependence. In this model we consider again two different substitution regimes, so the parameter vector is now (π, f, g, h1, h2), yielding a total number of 12 + 3 + 3 + 2 × 15 = 48 free parameters. We present in Figure 5 the posterior probability distribution over the set of possible alignment columns, match, insertion or deletion (the two match states have been merged) for every pair of aligned positions in a consensus alignment obtained at the convergence stage (last iterations) of saem algorithm. This posterior probability is quite similar to the one obtained with the context-dependent model. In order to further compare the two models, we have calculated a BIC-type criterion penalizing the maximum likelihood by the number of parameters. Specifically, we have computed BIC = -2L + k log(n), where L is the value of the log-likelihood at θML, n is the sample size and k is the total number of free parameters. We have taken n = 842, the length of the longest sequence, since the length of the alignment is unknown, but is equivalent to those of the observed sequences. Indexing by 1 the context-dependent model and by 2 the non context-dependent one, we get BIC1 = -2 • (-7128.50) + 78 log(842) = 14782.39 BIC2 = -2 • (-7337.40) + 48 log(842) = 14998.12, which means that the context-dependent model better fits this data set.

Discussion

We have proposed a model for statistical sequence alignment accounting for context-dependence in the substitution process. This model extends the pair hidden Markov model in the sense that, conditional on the alignment, the presence of particular nucleotides at the different sites of a sequence are no longer independent. We have adapted the traditional dynamic programing algorithms to the new framework, and we have proposed an efficient estimation method based on saem algorithm. We have obtained asymptotic results for maximum likelihood estimation on the new model, and through a simulation study, we have shown the accuracy of our algorithms on finite sample estimation. Finally, we have illustrated our methods with the alignment of the human alpha-globin pseudogene and its functional counterpart. We have compared the new model with a classical pair- HMM through a model selection criterion, concluding that taking into account the context-dependence on the substitution process improves the fitting of the data.

Figure 1 .

 1 Figure 1. Graphical representation of an alignment between two sequences X = AAT G and Y = CT GG. The displayed alignment is

h

  and h. For any b, c ∈ A, let hX (•) := a∈A h(•, a), hY (•) := a∈A h(a, •), hX (•|b, c) := a∈A h(•, a|b, c), hY (•|b, c) := a∈A h(a, •|b, c).

  , Ne(a|IX ) = 1{e k = IX , XN k = a} and similarly for Ne(a|IY ), Ne(a, b|M ) and also Ne(a, b|M M

Figure 2 .

 2 Figure 2. Distribution of saem estimates of the reduced parameter vector (α, β, γ, λ) over 100 simulation runs for the two simulated data sets (Data set 1 on top, Data set 2 on bottom). The real values of the parameters are displayed in dotted line.

Figure 3 .

 3 Figure3. saem estimates of the reduced parameter vector (α, β, γ, λ) over iterations in a single simulation run from the second data set. A logarithmic scale is used for the x-axis. The true parameter values are displayed in dotted line.

  G, namely φi(G, •)/ a∈A φi(G, a), i = 1, 2, φ = h, h. The probability of transition for G increases from 0.08 in the general context to 0.88 in the C C match context in the first substitution regime, and from 0.02 to 0.15 in the second.

Figure 4 .

 4 Figure 4. Posterior probability distribution of the alignment states at the maximum likelihood parameter estimate with the substitution context-dependent model.

Figure 5 .

 5 Figure 5. Posterior probability distribution of the alignment states at the maximum likelihood parameter estimate without taking into account possible substitution context-dependence.

Table 2

 2 Mean values and standard deviations of estimates over 100 simulation runs for Data Set 1. The estimation has been performed without taking into account the specific parametrization used to generate the data.

		Initial value True value	Estimate (sd)
	πMM	0.85	0.9238	0.9184 (0.0101)
	πMI X	0.075	0.0377	0.0407 (0.0074)
	πMI Y	0.075	0.0385	0.0410 (0.0077)
	πI X M	0.85	0.9424	0.8807 (0.0998)
	πI X I X	0.075	0.0385	0.0387 (0.0308)
	πI X I Y	0.075	0.0191	0.0805 (0.1034)
	πI Y M	0.85	0.9238	0.8788 (0.0982)
	πI Y I X	0.075	0.0377	0.0821 (0.0987)
	πI Y I Y	0.075	0.0385	0.0391 (0.0305)
	h(A, A)	0.0625	0.2148	0.2155 (0.0109)
	h(A, C)	0.0625	0.0036	0.0031 (0.0019)
	h(A, G)	0.0625	0.0036	0.0032 (0.0019)
	h(A, T )	0.0625	0.0029	0.0027 (0.0018)
	h(C, A)	0.0625	0.0036	0.0034 (0.0021)
	h(C, C)	0.0625	0.2634	0.2628 (0.0130)
	h(C, G)	0.0625	0.0044	0.0044 (0.0024)
	h(C, T )	0.0625	0.0036	0.0034 (0.0021)
	h(G, A)	0.0625	0.0036	0.0030 (0.0017)
	h(G, C)	0.0625	0.0044	0.0041 (0.0023)
	h(G, G)	0.0625	0.2634	0.2653 (0.0120)
	h(G, T )	0.0625	0.0036	0.0032 (0.0018)
	h(T, A)	0.0625	0.0029	0.0026 (0.0016)
	h(T, C)	0.0625	0.0036	0.0031 (0.0017)
	h(T, G)	0.0625	0.0036	0.0035 (0.0021)
	h(T, T )	0.0625	0.2148	0.2166 (0.0113)
	hC (A, A)	0.0625	0.1600	0.1626 (0.0169)
	hC (A, C)	0.0625	0.0112	0.0108 (0.0053)
	hC (A, G)	0.0625	0.0446	0.0454 (0.0111)
	hC (A, T )	0.0625	0.0092	0.0091 (0.0051)
	hC (C, A)	0.0625	0.0112	0.0102 (0.0056)
	hC (C, C)	0.0625	0.2055	0.2072 (0.0193)
	hC (C, G)	0.0625	0.0137	0.0126 (0.0056)
	hC (C, T )	0.0625	0.0446	0.0429 (0.0108)
	hC (G, A)	0.0625	0.0446	0.0440 (0.0098)
	hC (G, C)	0.0625	0.0137	0.0140 (0.0061)
	hC (G, G)	0.0625	0.2055	0.2050 (0.0174)
	hC (G, T )	0.0625	0.0112	0.0110 (0.0053)
	hC (T, A)	0.0625	0.0092	0.0086 (0.0047)
	hC (T, C)	0.0625	0.0446	0.0456 (0.0106)
	hC (T, G)	0.0625	0.0112	0.0111 (0.0062)
	hC (T, T )	0.0625	0.1600	0.1598 (0.0197)

Table 3

 3 Mean values and standard deviations of estimates over 100 simulation runs for Data set 2. The estimation has been performed without taking into account the specific parametrization used to generate the data.

		Initial value True value	Estimate (sd)
	πMM	0.85	0.9610	0.9559 (0.0071)
	πMI X	0.075	0.0194	0.0221 (0.0051)
	πMI Y	0.075	0.0196	0.0219 (0.0055)
	πI X M	0.85	0.9706	0.8477 (0.1262)
	πI X I X	0.075	0.0196	0.0199 (0.0221)
	πI X I Y	0.075	0.0098	0.1323 (0.1289)
	πI Y M	0.85	0.9610	0.8407 (0.1321)
	πI Y I X	0.075	0.0194	0.1367 (0.1379)
	πI Y I Y	0.075	0.0196	0.0226 (0.0252)
	h(A, A)	0.0625	0.2165	0.2168 (0.0114)
	h(A, C)	0.0625	0.0030	0.0027 (0.0015)
	h(A, G)	0.0625	0.0030	0.0027 (0.0014)
	h(A, T )	0.0625	0.0025	0.0024 (0.0015)
	h(C, A)	0.0625	0.0030	0.0026 (0.0017)
	h(C, C)	0.0625	0.2653	0.2673 (0.0122)
	h(C, G)	0.0625	0.0037	0.0033 (0.0020)
	h(C, T )	0.0625	0.0030	0.0027 (0.0016)
	h(G, A)	0.0625	0.0030	0.0023 (0.0015)
	h(G, C)	0.0625	0.0037	0.0033 (0.0020)
	h(G, G)	0.0625	0.2653	0.2643 (0.0096)
	h(G, T )	0.0625	0.0030	0.0023 (0.0013)
	h(T, A)	0.0625	0.0025	0.0024 (0.0015)
	h(T, C)	0.0625	0.0030	0.0025 (0.0016)
	h(T, G)	0.0625	0.0030	0.0027 (0.0015)
	h(T, T )	0.0625	0.2165	0.2196 (0.0102)
	hC (A, A)	0.0625	0.1588	0.1595 (0.0159)
	hC (A, C)	0.0625	0.0086	0.0069 (0.0038)
	hC (A, G)	0.0625	0.0505	0.0504 (0.0096)
	hC (A, T )	0.0625	0.0071	0.0067 (0.0042)
	hC (C, A)	0.0625	0.0086	0.0087 (0.0045)
	hC (C, C)	0.0625	0.2053	0.2053 (0.0199)
	hC (C, G)	0.0625	0.0105	0.0103 (0.0046)
	hC (C, T )	0.0625	0.0505	0.0504 (0.0101)
	hC (G, A)	0.0625	0.0505	0.0515 (0.0102)
	hC (G, C)	0.0625	0.0105	0.0097 (0.0052)
	hC (G, G)	0.0625	0.2053	0.2036 (0.0191)
	hC (G, T )	0.0625	0.0086	0.0078 (0.0043)
	hC (T, A)	0.0625	0.0071	0.0078 (0.0041)
	hC (T, C)	0.0625	0.0505	0.0515 (0.0098)
	hC (T, G)	0.0625	0.0086	0.0084 (0.0040)
	hC (T, T )	0.0625	0.1588	0.1613 (0.0166)

Table 4

 4 Mean values and standard deviations of estimates over 100 simulation runs for the two data sets. The estimation has been performed by numerical optimization on the reduced parameter vector (α, β, γ, λ).

			Data set 1
		Initial value True value	Estimate (sd)
	α	0.8	0.4	0.400 (0.0750)
	β	0.25	0.2	0.1985 (0.0372)
	γ	0.1	0.06	0.0602 (0.0087)
	λ	0.08	0.04	0.0402 (0.0040)
			Data set 2
		Initial value True value	Estimate (sd)
	α	0.8	0.5	0.5091 (0.0772)
	β	0.25	0.15	0.1493 (0.0296)
	γ	0.1	0.05	0.0499 (0.0079)
	λ	0.08	0.02	0.0196 (0.0024)

Table 5

 5 Relative substitution probabilities for

Appendix

The need for a stochastic approximation Let us explain why em algorithm may not be applied in pair-HMM. Denoting by Lnm the random value s 1 such that Zs = (n, m) (the first and only hitting time for the point (n, m), which is not necessarily finite), the complete log-likelihood writes log P θ (X1:n, Y1:m, Lnm, ε1:L nm ) = n+m s=n∨m 1L nm=s u∈E