Alain Audouard

Enric Canadell

Jean-Yves Fortin

Tatyana G Prokhorova

Vladimir N Laukhin

David Vignolles

Eduard B Yagubskii

De Haas-van Alphen oscillations in the compensated organic metal alpha-'pseudo-kappa'

Keywords: 71, 18, +y Fermi surface: calculations and measurements; effective mass, g factor -71, 20, Rv Polymers and organic compounds

Introduction

Provided that no phase transition occurs as the temperature is lowered, Fermi surface (FS) of two-dimensional (2D) organic metals is generally rather simple and achieve model systems for quantum oscillations physics. Indeed, in the numerous cases where the compound possesses two carriers (generally holes) per unit cell, the FS originates from a single orbit with an area equal to that of the first Brillouin zone [1]. In an extended zone scheme, these orbits overlap either along (i) one or (ii) two directions yielding in magnetic field either (i) the model linear chain of orbits coupled by magnetic breakdown (MB) proposed by Pippard in the early sixties [2,[START_REF] Shoenberg | Magnetic Oscillations in Metals[END_REF] for which all the orbits are of the same (hole) type or (ii) a set of compensated electron-and hole-type orbits (i.e. the sum of the hole-type orbits cross section is equal to that of the electron-type orbits).

The former case (i) has been widely studied: it is known that quantum oscillations spectra are strongly affected by field-induced chemical potential oscillations yielding many frequency combinations, the field and temperature dependence of which cannot be accounted for by the Lifshitz-Kosevich (LK) formula [START_REF] Wosnitza | Fermi Surfaces of Low-Dimensional Organic Metals and Superconductors[END_REF]5,6]. In contrast, to our best knowledge, the second case (ii) has only received little attention up to now even though MB between hole-and electron-type orbits is relevant, for instance, for recently studied oxide superconductors [7,8,9]. Numerical resolution of the grand potential equation for the case (ii) reveals a e-mail: alain.audouard@lncmi.cnrs.fr b e-mail: fortin@ijl.nancy-universite.fr that chemical potential oscillations are strongly damped for compensated orbits, even in the case of 2D metals [10,11]. As a consequence, contrary to case (i), the LK formula is predicted to account for the field and temperature dependence of the oscillations amplitude in case (ii). However, depending on the MB probability and taking into account the opposite sign of electron and hole orbits, a given and Solv is an organic solvent) [13,14].

As reported in Fig. 1(a), the unit cell of α-'pseudo-κ'- [13] contains two different donor planes with different packing. One of them, with a 'pseudo-κ' structure is insulating while the other with an α-type structure is metallic. As a consequence, spacing between conducting layers is as large as 3.64 nm, ensuring negligibly small interlayer transfer integral, hence avoiding effects due to FS corrugation [15]. According to band structure calculations, the FS of the metallic plane, displayed in Fig. 1(b), is composed of two electron-type and one hole-type compensated orbits (i.e. the hole orbit area is twice the electron orbits area). In the extended zone scheme, this set of three orbits is isolated from the other sets. In other words, we are not dealing with a network of coupled orbits. In contrast, each of the orbits within a set is liable to be connected to the other by MB. This feature is shared by the FS studied in Ref. [10], which is composed of one hole and one electron compensated orbit.

(ET) 4 H 3 O[Fe(C 2 O 4 ) 3 ]•(C 6 H 4 Br 2 )
The aim of this paper is, in the first step, to provide calculation of the Fourier amplitudes of de Haas-van Alphen (dHvA) oscillations spectra for the FS of Fig. 1(b). In the second step, magnetic torque oscillations of α-'pseudo-κ'-

(ET) 4 H 3 O[Fe(C 2 O 4 ) 3 ]•(C 6 H 4 Br 2 )
, yielding dHvA spectra, are studied in magnetic fields up to 55 T. Whereas magnetoresistance data of this organic metal, measured in magnetic fields below 16 T reveal only one frequency attributed to the hole orbit [13], the two frequencies cor-responding to electron and hole orbits area are observed.

It is shown that the field and temperature dependence of the Fourier amplitudes are in agreement with the reported calculations.

Experimental

Crystals were synthesized by electrocrystallization technique as reported by Zorina et al. [13]. Two crystals denoted as crystal #1 and #2 hereafter were studied. Their approximate dimensions were 0.12 × 0.1 × 0.04 mm 3 .

Magnetic torque was measured with a commercial piezoresistive microcantilever, in pulsed magnetic fields of up to 55 T with a pulse decay duration of 0.32 s. Variations of the cantilever piezoresistance were measured in the temperature range from 1.4 K to 4.2 K with a Wheatstone bridge with an ac excitation at a frequency of 63 kHz.

The angle between the normal to the conducting plane and the magnetic field direction was θ = 15 • and 32 • for crystal #1 while θ was varied from 15 • to 71 • thanks to a rotating sample holder for crystal #2.

Calculations of the Fourier amplitudes

In the framework of the Lifshitz-Kosevich and Falicov-Stachowiak models [START_REF] Shoenberg | Magnetic Oscillations in Metals[END_REF]16], the oscillatory part of the magnetization M osc for a set of 2D orbits η can be written as

M osc = - η p F η π A pη sin(2πp F η B cos θ ) ( 1
)
where θ is the angle between the normal to the conducting plane and the field direction. The index η stands for all the closed orbits allowed by the FS topology, including MB orbits, with fundamental frequencies F η . They do not include orbits which can be associated with an harmonic of a simpler trajectory, the latter being accounted for by the index pF η (with p > 1) where p is the harmonic order. The amplitude A pη of the Fourier component with frequency pF η depend on parameters such as temperature, magnetic field, effective masses (m η ), Dingle temperatures (T Dη ), MB field (B 0 ), and effective Landé factors (g * η ). These amplitudes can be expressed as

A pη = (-1) psη R MB pη R pη /(pm η )
. MB damping factor is given by

R MB pη = C pη (ip 0 ) ptη q pbη 0
, where C pη is the symmetry factor of the orbit pη, t η and b η are the number of tunnelings and reflections, respectively, encountered by a quasi-particle during its path and 2s η is the number of turning points around the orbit η. The tunneling and reflection probabilities at a MB junction are given, in agreement with the Chambers approximation, by p 2 0 = exp(-B 0 /B cos θ) and

q 2 0 = 1 -p 2 0
, respectively [START_REF] Chambers | Proc. Phys. Soc[END_REF]. The damping factor R pη can be written as the product of thermal, Dingle and spin

damping factors (R pη = R T pη R D pη R s pη ) which are given by R T pη = X pη / sinh(X pη ), R D pη = exp(-u 0 T Dη pm η /B cos θ),
and R s pη = cos(pπg * η m η /2 cos θ), respectively, where

X pη = u 0 T pm η /B cos θ and u 0 = 2π 2 k 2 B m e /e .
In the case of the FS of Fig. with S h = -2S e , all Fourier components with frequencies [START_REF] Shoenberg | Magnetic Oscillations in Metals[END_REF]16], where n e and n h are the (positive) numbers of electron and hole orbits involved, contribute. Therefore, the dominant frequencies F e , 2F e = F h and 3F e we will consider in the following, arise from infinite orbits combinations with n e = 2n h ± 1, n e = 2n h ± 2 and n e = 2n h ± 3, respectively. To obtain the amplitude of these Fourier components, the Fermi surface is modeled in Fig. 2 by a linear finite chain of e and h orbits.

F nee+n h h = |n e F e -n h F h | = |n e -2n h |F e ,

and effective masses m nee+n

h h = n h m h + n e m e
To count all the allowed paths from one of the arbitrary starting and ending points α, β, β, γ in Fig. 2, we consider paths of length n, where n is an integer corresponding to the number of elementary steps from one point to another that form and define the complete closed trajectory.

Turning points are located at the arrows in Fig. 2 and going through one of these is tantamount to adding a factor i in the wave-function amplitude. We introduce the vector amplitude |α, β, β, γ corresponding to the state of the quasi-particle. Any closed path is therefore defined by its length n and its starting point and state s 0 among the set (α, β, β, γ). Initially the state |s 0 is filled with zeroes except for the component corresponding to s 0 , for example |α = |1, 0, 0, 0 . We then define each step factor of the trajectory by moving from one point among the set (α, β, β, γ) to the next one, by following the direction imposed by the field, with the adequate changes in the phase and amplitude. For example, as shown in Fig. 2,

β = -ip 0 α exp(iS e /2 + iS h /2) + iq 0 β exp(iS h /2). Intro-
ducing parameters x = exp(iS e /2) and y = exp(iS h /2), with y = 1/x 2 , we can write a transfer matrix for all the elementary steps on the Fermi surface:

T (x, y) =            -q 0 x 2 -ip 0 xy 0 0 0 0 iq 0 y -p 0 x -p 0 x iq 0 y 0 0 0 0 -ip 0 xy -q 0 x 2           
Since we consider closed paths and conserved current, the output vector is the same as the input vector |s 0 .

Then the number of all possible paths of length n and starting from a point s 0 is given by the number s 0 |T (x, y) n |s 0 .

The generalized damping factor A pη in Eq. ( 1) is computed for all possible orbit combinations contributing to the same frequency pF η . The combinatorial factor (-1) psη R MB pη in A pη , combining the different ways of drawing the orbits pη on the FS, is precisely defined as the coefficient of x 2ne y 2n h , as well as x -2ne y -2n h , in the following polynomial function

A(x, y) = 1 0 dz z s0 n≥0 W s0 s 0 |T n ( √ zx, zy)|s 0 (2)
weighted by damping factors for each orbit, and corresponding to integers (n e , n h ). W s0 is the weight of the point s 0 in the quasi-particle path. We will take S e (corresponding to the frequency F e ) has a weight proportional to 10p 4 0 q 0 /5 = 2p 4 0 q 0 (the factor 5 = n e + 2n h is the number of cyclic permutations of the same orbit

W α = W γ = 1 and W β = W β = 1/
ee √ he √ h = e √ he √ he = • • • ). It is different from the orbit ee √ h √
he for example, which has a weight proportional to -2p 3 0 q 2 0 . We find that the first three amplitudes can be written as

A e = - 2q 0 m e R e - 2p 2 0 q 0 m e + m h R e+h + 2[p 4 0 q 0 -p 2 0 q 3 0 ] 3m e + m h R 3e+h + • • • , (3)
A 2e = q 2 0 2m e R 2e - q 2 0 m h R h + 3p 4 0 q 2 0 -2p 2 0 q 4 0 2(m e + m h ) R 2e+2h + 2p 2 0 q 4 0 -3p 4 0 q 2 0 4m e + m h R 4e+h + • • • ,(4) A 3e = - 2q 3 0 9m e R 3e + 2p 2 0 q 3 0 m e + 2m h R e+2h + • • • . (5)
where R pη = R T pη R D pη R s pη . As an example, the first term -(2q 3 0 /9m e )R 3e of A 3e in Eq. 5, comes from the third harmonics of the e orbit. Indeed, according to Eq. 1, this factor is given by F e (-2q 3 0 )R 3e /(3m e ) which can be rewritten as 3F e (-2q 3 0 )R 3e /(9m e ). Magnetization can then be expanded as

M osc = - F e π A e sin 2π F e B - 2F e π A 2e sin 2π 2F e B - 3F e π A 3e sin 2π 3F e B + • • • (6)
Examples of Fourier spectra deduced from Eq. 6 are reported in Fig. 3 increase which may lead to errors in effective mass determination [11]. Contribution of the orbit 2e is substantial in any case, which indicates that this orbit needs to be considered, even in the case of large MB field, scattering rate and temperature, for correct data analysis. In contrast, complex orbits such as 2e + 2h have significant contribution for clean compounds with moderate MB field at low temperature, only. Similar conclusions can be derived regarding the component with frequency F e for which the contribution of e + h is always significant in the range explored, except for large MB fields, while 3e + h is negligible at high temperature and large scattering rate. Of course, all the considered amplitudes also depend on the field-and temperature-independent spin damping factor R s pη through the product g * η m η .

Results and discusion

This section is devoted to de Haas-van Alphen oscilla- Long dashed, short dashed, dotted and dash-dotted lines are the components linked to the h, 2e, 2e + 2h and 4e + h orbits, respectively. The magnetic breakdown damping factors for the orbits 2e + 2h and 4e + h are given by R 2e+2h = 3p 4 0 q 2 0 -2p 2 0 q 4 0 and R 4e+h = 2p 2 0 q 4 0 -3p 4 0 q 2 0 , where q0 and p0 are the reflection and magnetic breakdown probability, respectively.

tions of α-'pseudo-κ'-(ET) 4 H 3 O[Fe(C 2 O 4 ) 3 ] •(C 6 H 4 Br 2 ),
the FS of which (see Fig. the hole and electron orbits cross section area S h = 8.9 ± 0.2 % and S e = 4.4 ± 0.3 %, respectively, of the FBZ area. These data are in agreement with magnetoresistance data (S h = 8.9 % of the FBZ area) and band structure calculations at room temperature (S h = 7.6 % of the FBZ area and S e =S h /2) [13]. Compared to the Fourier component with the frequency F h , that with the frequency F e , which was not observed in Ref. [13], has a relatively small amplitude at θ = 15 • and 32 • , and can only be reliably studied at θ = 49 • .

It must be noticed that the studied compound contains magnetic Fe 3+ ions. As reported in the case of λ-(BETS) 2 FeCl 4 (where BETS stands for bis-ethylenedithiotetraselenafulvalene) [19,20], these magnetic ions induce an exchange field leading to angle-dependent splitting of the oscillation frequency which may alter the oscillations amplitude. Nevertheless, such a splitting is not observed neither in the considered compound nor in the compound with the same composition and β" structure in fields of up to 55 T [21,22,23]. This could be due to the lower, by a factor of two, Fe 3+ concentration in the present case hence to a reduced exchange coupling constant.

A naive analysis of the temperature dependence of A τ h , relevant to F h , assuming that only one orbit contributes (i.e. through Eq. 1) yields an effective mass m h (θ)×cos(θ) of 0.96 ± 0.05, 0.93 ± 0.02 and 0.69 ± 0.07 at θ = 15 • , 32 • and 49 • , respectively. Within the same hypothesis, the magnetoresistance data at θ = 0 • of Ref. [13] yield m h = 1.11 ± 0.04 [START_REF] The | mass value m h = 1.26[END_REF]. In other words, the product m h (θ) × cos(θ) would monotonously decreases as θ increases. Hence, at variance with the angle dependence of the frequency, the cosine law (m h (θ) = m h (0)/ cos θ) consistent with a 2D FS, would not be followed for the effective mass within this assumption. This result strongly suggests that Eq. 1 is unable to account for the temperature dependence of the amplitude. Besides, it can be checked that the field dependence of the amplitude cannot be accounted for by Eq. 1 as well. As discussed in the preceding section and reported hereafter, other orbits (such as 2e, 2e + 2h, etc.)

with effective masses different from each other enter the oscillatory spectra.

According to Eqs. 6, oscillatory torque amplitudes A τ e and A τ h , of the Fourier components with frequencies F e and F h = 2F e observed in Fig. 4, are given by:

A τ e = τ 0 tan(θ)B F e π A e , A τ h = τ 0 tan(θ)B 2F e π A 2e , (7)
respectively, where τ 0 is a prefactor depending on the cantilever stiffness, crystal geometry, etc. and θ is the angle between the field direction and the normal to the conducting plane. Amplitudes A e and A 2e in Eqs. 7 are given by Eqs. 3 and 4, respectively. We will limit ourselves to small n h and n e values. Namely, in addition to the basic electron and hole orbits, only the second harmonic of the electron orbit (2e) and the MB orbits composed of one hole and one electron orbit (e + h), 3 electron and one hole orbits (3e + h), 2 electron and 2 hole orbits (2e + 2h)

and 4 electron and one hole orbits (4e + h) are taken into account. In short, the contribution of the MB orbits and harmonics composed of more than 5 individual orbits are neglected. It must be kept in mind that spin damping fac-tors may influence the sign of a given contribution or, in other words, induce a π dephasing, as observed on either side of a spin zero angle [5]. .27 [27]).

Composites orbits, e+h and 2e have significant contributions to A τ e and A τ h , respectively, as observed in Fig. 5(c).

Nevertheless, higher order terms (linked to 2e + 2h and 3e + h) are small, due to rather large Dingle temperatures. Finally, as reported in Fig. 5c for the data at 49 • , contributions of e and 3e + h to A e have opposite signs.

This result which also holds at 15 • and 32 • , explains why this amplitude is small compared to A τ h .

Summary and conclusion

Fourier

  component can involve combinations of an infinite set of orbits, the contribution of which must be included in the Fourier amplitudes calculation. From the experimental side, only few quasi-2D compensated metals have been synthesized yet. Such FS are observed in e.g. (BEDO-TTF) 2 ReO 4 •H 2 O [12] (where BEDO-TTF stands for the bis-ethylenedioxy-tetrathiafulvalene molecule) and, more recently synthesized, (ET) 4 H 3 O[Fe(C 2 O 4 ) 3 ]• Solv (where ET stands for the bis-ethylenedithio-tetrathiafulvalene molecule

Fig. 1 .

 1 Fig. 1. (color on line) (a) Crystal structure of α-'pseudo-κ'-(ET)4H3O[Fe(C2O4)3]•(C6H4Br2). α-type and 'pseudo-κ'-type layers are conducting and insulating, respectively. (b) Fermi surface relevant to the α-type layers, according to Zorina et al.[13]. Blue and green areas mark electron-and hole-type orbits, respectively. Ellipses in dashed lines mark the orbits from which the Fermi surface originates. Labels Γ , X and Y refer to the (0, 0), (a * /2, 0) and (0, b * /2) points of the first Brillouin zone, displayed as a rectangle in solid line.

  1(b), which is modeled by Fig. 2, and owing to the orbits compensation predicted by band structure calculations, each electron (e) and hole (h) orbit contributes to the dHvA oscillations spectrum with the frequency F e and F h = 2F e , respectively. Besides, MB orbits composed of several individual orbits (η = n e e + n h h) are also liable to contribute to the spectrum. Indeed, all contributing trajectories for a given frequency are accounted for by considering the amplitude and phase variation exp(iS η ) = exp[2iπhA η /eB)] of the wave-function, where A η is the area of the trajectory in the Brillouin zone, which is directly identified to the frequency F η = |A η |h/e. Taking into account the opposite sign of electron and hole surface area, S nee+n h h = n e S e + n h S h

Fig. 2 .

 2 Fig. 2. Compensated Fermi surface with classical representation of the quasi-particle orbits. Phase variation of the wavefunction is given by exp(iSe) around the electronic orbit e and exp(iS h ) around the hole orbit h, where S e(h) is the electron (hole) surface area. Either magnetic breakdown occurs with the probability amplitude ip0 between two orbits or reflection is allowed with the probability amplitude q0 = (1p 2 0 ) 1/2 .

2 .

 2 Indeed points β and β are equivalent in the sense that they define the same set of paths and belong to the same orbit. Finally, the integration over z performed so that to remove the number of cyclic permutations n e + 2n h of the same fundamental orbit by adding a compensatory factor 1/(n e + 2n h ), and to take into account the harmonics coefficients with a weighting factor 1/p. Each path can be simply decomposed using e (or x) for a portion of trajectory around the electron orbits, and √ h (or y) around the hole orbit, and z is added for each e or √ h encountered on the trajectory. For example the orbit ee √ he √ h (corresponding to n e = 3, n h = 1) with effective mass 3m e + m h and phase

  (a) for various temperatures and MB fields. Landé factors (g * e = 2, g * h = 2) and effective masses (m e = 1, m h = 1), which otherwise are close to those deduced from the data reported in the next section, are chosen so that the absolute value of the spin damping factors is equal to 1 (|R s pη |=1) in order to avoid any spurious effect due to spin-zero phenomenon. First, A 3e is always small compared to A e and A 2e . Corresponding field-dependent amplitudes A 2e are given in Figs. 3(b), (c), (d). Despite the effective mass m η of a given η orbit increases as the number of individual orbits involved increases, a clear contribution of the orbits 2e, 2e + 2h and 4e + h to the amplitude is observed in Fig. 3(b). Their relative contributions decrease as the Dingle temperatures (not shown), the temperature (see Fig. 3(c)) and the MB field (see Fig. 3(d))

Fig. 3 .

 3 Fig. 3. (color on line) (a) Fourier spectra of the oscillatory magnetization calculated at a mean field of 40 T from Eq. 6 with me = 1, m h = 1, g * e = 2, g * h = 2, TDe = 1 K, T Dh = 1 K. Labels b, c, d correspond to T = 1.5 K and B0 = 2 T, T = 5 K and B0 = 2 T, T = 1.5 K and B0 = 100 T, respectively. (b), (c), (d) Field dependence of the amplitude of the Fourier component with frequency F h = 2Fe corresponding to Labels b, c, d, respectively, of (a). Solid lines are calculated from Eq. 4.

  1(b)) is relevant to the above calculations. Note that this FS can be considered as resulting from the hybridization of a series of ellipses centered at the Γ point. The area of these ellipses, which is related to the total number of holes per unit cell in the HOMO bands, is twice the area of the cross section of the FBZ because the repeat unit of the layer contains eight ET molecules with an average charge of +1/2, i.e. a total of four electrons per repeat unit. It is interesting to point out how this FS differs from that of other α-type ET salts, such as the α-(ET) 2 [MHg(SCN) 4 ] family. As pointed out by Mori et al. [18], the closed and open portions of the

Fig. 4 .

 4 Fig. 4. Fourier analysis in the field range 25-55 T of the magnetic torque data reported in the insert for various directions of the magnetic field with respect to the normal to the conducting plane (angle θ).

Fig. 5 .

 5 Fig. 5. ln(|A τ |) vs 1/B at various temperatures for (a) A τ h (solid circles) and A τ e (solid squares) at θ = 49 • , (b) A τ h at θ = 15 • (solid squares) and 32 • (solid circles) and (c) field dependence of A τ h (solid circles) and A τ e (solid squares) at θ = 49 • at 1.5 K. Solid lines are fits of Eqs. 7, obtained with me = 0.93, m h = 0.88, g * e = 2.27, g * h = 2.06, TDe = 5 K, T Dh = 4 K and B0 = 2 T. Long dashed, short dashed, dotted and dash-dotted lines in (c) are the contributions of the various components entering the fittings (see text and Fig. 3).

Fig. 6 .

 6 Fig. 6. Angle dependence of the amplitude |A τ h | at T = 1.5 K and B = 40 T. Solid line is the best fit of Eq.7 to the data, obtained with the same set of parameters as in Fig. 5. Long dashed, short dashed and dotted lines in the insert are the contributions of h, 2e and 2e + 2h, respectively.

Fig. 5

 5 Fig. 5 displays field dependence at various temperatures of either A τ e or A τ h for θ = 15 • and 32 • (crystal #1) and 49 • (crystal #2) while Fig. 6 displays the angle dependence of A τ h at T=1.5 K and B = 40 T for crystal #2. Solid lines in these figures are the best fits of Eqs. 7 to the Fourier amplitudes: the same set of parameters stands for all the data, namely m e = 0.93 ± 0.04, m h = 0.88 ± 0.04, g * e = 2.27 ± 0.12, g * h = 2.06 ± 0.09. Dingle temperatures and MB field, which jointly govern the field dependence of the amplitude are obtained with a large uncertainty.Dingle temperatures are rather large, of the order of few K, whereas B 0 is in the range between 0 and 4 T which is rather small, in agreement with the FS of Fig.1. Owing to the small size of e and h orbits, compared to that of

Fourier spectra of

  de Haas-van Alphen oscillations of compensated two-dimensional metals with Fermi surface composed of one hole and two electron components have been considered. The two main Fourier components observed have frequencies F e and F h = 2F e , corresponding to the electron and hole orbits area. Nevertheless, it is demonstrated that, taking into account the opposite sign of electron and hole orbits, a given Fourier component involves combination of several orbits, the contribution of which must be included in the calculations. Such FS, which is a textbook case, is observed in the strongly 2D organic metal α-'pseudo-κ'-(ET) 4 H 3 O[Fe(C 2 O 4 ) 3 ]•(C 6 H 4 Br 2 ). Magnetic torque oscillations of this compound have been studied up to 55 T for various directions of the magnetic field with respect to the conducting plane. It is demonstrated that data analysis performed assuming that only single electron and single hole orbits contribute to Fourier components with frequency F e and 2F e , respectively, cannot account for the data. In other words, additional orbits generated by tunneling and reflection at magnetic breakdown junctions must be taken into account. Calculations are in good quantitative agreement with the data.

This work has been supported by EuroMagNET II under the EU Contract No. 228043, and MINECO-Spain (Projects FIS2012-37549-C05-05 and CSD 2007-00041).


Catherine Belleannée
email: catherine.belleannee@irisa.fr

Olivier Sallou
email: olivier.sallou@irisa.fr

Jacques Nicolas
email: jacques.nicolas@irisa.fr

Logol: Expressive Pattern Matching in sequences Application to Ribosomal Frameshift Modeling

Most of the current practice of pattern matching tools is oriented towards finding efficient ways to compare sequences. This is useful but insufficient: as the knowledge and understanding of some functional or structural aspects of living systems improve, analysts in molecular biology progressively shift from mere classification tasks to modeling tasks. People need to be able to express global sequence architectures and check various hypotheses on the way their sequences are structured. It appears necessary to offer generic tools for this task, allowing to build more expressive models of biological sequence families, on the basis of their content and structure. This article introduces Logol, a new application designed to achieve pattern matching in possibly large sequences with customized biological patterns. Logol consists in both a language for describing patterns, and the associated parser for effective pattern search in sequences (RNA, DNA or protein) with such patterns. The Logol language, based on an high level grammatical formalism, allows to express flexible patterns (with mispairings and indels) composed of both sequential elements (such as motifs) and structural elements (such as repeats or pseudoknots). Its expressive power is presented through an application using the main components of the language : the identification of -1 programmed ribosomal frameshifting (PRF) events in messenger RNA sequences. Logol allows the design of sophisticated patterns, and their search in large nucleic or amino acid sequences. It is available on the GenOuest bioinformatics platform at http://logol.genouest.org. The core application is a command-line application, available for different operating systems. The Logol suite also includes interfaces, e.g. an interface for graphically drawing the pattern.

Background

During last decade, a number of pattern matching tools have been proposed, and some of them are used extensively and helpful. Depending on systems under study, modeling needs may vary from looking for all exact occurrences of a given string in a protein bank, to looking for approximated occurrences of a given transposon in a full genome, or locating pseudoknots in a RNA sequence. This section proposes a quick overview of the existing software diversity, which shows there exists still room for a new pattern matching tool, both flexible (high "genericity") and with the capacity to represent complex structures (high expressive power).

General Purpose Pattern Matching

Some tools have been designed for the analysis of several types of sequences (DNA, RNA, proteins) with an generic expressiveness, i.e. without targeting the recognition of a particular motif family. Among these general tools, two tendencies can be observed, efficiency-oriented and expressiveness-oriented software.

One of the most advanced software from the point of view of efficiency is the Vmatch suite (http://www.vmatch.de) that offers a large variety of search facilities in very large sequences. It is based on a careful implementation of enhanced suffix trees for the computation of a sequence index that provides a fast access to every substring in that sequence. If the search for a motif contains some rare substrings, this technique is particularly efficient. The software Vmatch is the core search engine used in a number of more specialized tools working on specific sequence structures (e.g. "tandem-repeats" or LTR retrotransposons). Another highly generic tool is Biogrep [START_REF] Jensen | Biogrep: A multi-threaded pattern matcher for large pattern sets[END_REF], designed with the objective of quickly recognizing a large set of simple motifs (typically more than 100) in biological sequence banks. Biogrep allows queries in the POSIX language, a standard format of extended regular expressions, and can look for patterns in parallel on a set of processors.

The other approach for the analysis of biological sequences is more concerned with modeling the peculiarities of biological objects in the most relevant and expressive way. A major contribution in this respect is the work of D. Searls who laid the foundations for the research in this domain. He was the first to supervise developments allowing users to design biological grammars and to apply them for the large scale analysis of their genomic sequences [START_REF] Searls | A syntactic pattern recognition system for DNA sequences[END_REF][START_REF] Dong | Gene structure prediction by linguistic methods[END_REF][START_REF] Searls | String variable grammar: A logic grammar formalism for the biological language of DNA[END_REF]. D. Searls has introduced a very practical object in algebraic grammars, the string variable, which allows to elegantly express the notion of copy (either direct or reverse). He has implemented the resulting logical formalism, called SVG -for String Variable Grammars-, in the (no longer available) Genlang tool [START_REF] Dong | Gene structure prediction by linguistic methods[END_REF]. The direct copy (e.g. X...X) allows to search for two occurrences of a same unknown string, using optionally some indication on the string size. The reverse copy (e.g. X... ∼X) introduces in addition the notion of reverse complement and allows thus the representation of biological palindromes like stem-loops (Stem, Loop, ∼Stem) or pseudo-knots (Stem1, Loop1,Stem2, Loop2, ∼Stem1, Loop3,∼Stem2). Genlang, Stan [START_REF] Nicolas | Suffix-tree analyser (stan): looking for nucleotidic and peptidic patterns in chromosomes[END_REF] (developed in our research team), Patscan [START_REF] Dsouza | Searching for patterns in genomic data[END_REF] and Patsearch [START_REF] Pesole | Patsearch: a pattern matcher software that finds functional elements in nucleotide and protein sequences and assesses their statistical significance[END_REF] are all tools belonging to this family. Thanks to string variables and other additional components, these languages offer the possibility to combine easily in a single model informations on the sequence and on the structure of a molecule.

Dedicated Pattern Matching

It is not possible to provide here an exhaustive review of the profusion of specific tools that have been made available to bioanalysts. Some are specific to a sequence family and others to a particular motif type. A famous one dedicated to proteins is ScanProsite [START_REF] De Castro | Scanprosite: detection of prosite signature matches and prorule-associated functional and structural residues in proteins[END_REF], where motifs are built upon regular expressions that are searched either by a query in a precomputed database or with the algorithm ps-scan [START_REF] Gattiker | Scanprosite: a reference implementation of a prosite scanning tool[END_REF]. A number of tools are dedicated to RNA sequences, in response to the increasing needs of structure exploration in the complex RNA world boosted by the recent importance of non coding RNA studies. For instance, RNAmotif [START_REF] Macke | Rnamotif, an rna secondary structure definition and search algorithm[END_REF], RNAbob [START_REF] Eddy | Rnabob: a program to search for rna secondary structure motifs in sequence databases[END_REF], Hypasearch [START_REF] Graf | HyPaLib: a Database of RNAs and RNA Structural Elements defined by Hybrid Patterns[END_REF][START_REF] Strothmann | The syntax and semantics of a language for describing complex patterns in biological sequences[END_REF] and Palingol [START_REF] Billoud | Palingol: a declarative programming language to describe nucleic acids' secondary structures and to scan sequence database[END_REF] have been designed for the description of patterns as a succession of stems and loops, usually offering the possibility of choosing either a standard Watson-Crick pairing (A-U, G-C) or a pairing including Wobble (A-U, G-C, G-U). A more recent tool in this category, Structator [START_REF] Meyer | Structator: fast index-based search for rna sequencestructure patterns[END_REF], significantly improves the parsing time by making use of an index structure that is suited for the analysis of palindromic structures, the affix arrays. Patterns may also contain some sequence information on words that have to be present in particular places of the stems or the loops. RNAmotif is probablly the most popular in this category.

Logol Language

In this landscape, we designed Logol, a new general grammatical pattern matching tool, in order to greatly enhance the range of admissible patterns.

Basics: a Grammatical Model with Constraints

• String variable grammars: With the objective of being a general and expressive language allowing a natural expression of composite patterns, the Logol language has been designed on the basis of String Variable Grammars (SVG) introduced by D.Searls [START_REF] Searls | A syntactic pattern recognition system for DNA sequences[END_REF][START_REF] Dong | Gene structure prediction by linguistic methods[END_REF][START_REF] Searls | String variable grammar: A logic grammar formalism for the biological language of DNA[END_REF]. As already mentioned in the previous section, while it is easy to express motifs and gaps by regular expressions (e.g. PROSITE [START_REF] De Castro | Scanprosite: detection of prosite signature matches and prorule-associated functional and structural residues in proteins[END_REF]), SVG allow to express structures beyond the capabilities of regular languages such as palindromes (e.g. in stem-loops and pseudo-knots) and repeats (duplicated substrings), which are recorded by string variables. In fact, SVG and Logol grammars lay even beyond the possibilities of context-free grammars (XML-like), in a class that A. Joshi called "mildly context sensitive languages" [START_REF] Joshi | The convergence of midly contextsensitive grammars[END_REF]. Starting from the sound basis of SVG grammars, the Logol language proposes several extensions -most notably by adopting a constraint approach-with the goal to allow the expression of realistic biological motifs. The rest of the section introduces its main constituents.

• First steps in Logol: Let us first present a very simple Logol grammar: mod1()==*> SEQ1 mod1()==>"aaa" The first line is the top level instruction (called 'rule'). The rule is identified by the constant '==*> SEQ1' and it triggers the parsing of a sequence for a particular grammatical model (here, mod1). The second line provides the model (i.e. pattern) definition itself, here the string "aaa". It triggers the search for all occurrences of "aaa" in a genomic sequence.

The next grammar describes a slightly more interesting pattern made of two distant copies of a same string. The size of the string (in range [START_REF] Dsouza | Searching for patterns in genomic data[END_REF][START_REF] Forest | Modélisation et détection automatique de sites de décalage de cadre en -1 dans les génomes eucaryotes[END_REF]) and the distance between the two copies (in range [START_REF] Abouelhoda | Replacing suffix trees with enhanced suffix arrays[END_REF][START_REF] Graf | HyPaLib: a Database of RNAs and RNA Structural Elements defined by Hybrid Patterns[END_REF]) are bounded but the content of the string itself is kept free: mod2()==*> SEQ1 mod2() ==> X1:# [START_REF] Dsouza | Searching for patterns in genomic data[END_REF][START_REF] Forest | Modélisation et détection automatique de sites de décalage de cadre en -1 dans les génomes eucaryotes[END_REF], .*:{# [START_REF] Abouelhoda | Replacing suffix trees with enhanced suffix arrays[END_REF][START_REF] Graf | HyPaLib: a Database of RNAs and RNA Structural Elements defined by Hybrid Patterns[END_REF]}, X1 The model mod2() reads as follows: X1 denotes a string variable; any string made of 5 to 8 letters can be an instance of X1. '.*' denotes a space ('gap'). It is constrained to have a size between 1 and 10 characters. A second occurrence of X1 is waited for after the gap. For example, acuggcccgacuggcacuggc is an instance of this pattern on the input sequence uucagacuggcccgacuggcacuggccac, with X1 = acuggc.

• All matches: Logol returns all the instances of a model. For instance, when launched on the sequence cagaaaacgccgaaacuggc with the model "aaa", it returns the three possible matches: in position 3, 4 and 12.

• A powerful feature: instance saving. The Logol language supports an alternative way to express the pattern mod2. That is: mod2() ==> X1:{# [START_REF] Dsouza | Searching for patterns in genomic data[END_REF][START_REF] Forest | Modélisation et détection automatique de sites de décalage de cadre en -1 dans les génomes eucaryotes[END_REF],_IX1}, .*:{#[1,10]}, ?IX1. In this case, the string corresponding to the occurrence value of X1 is saved (using '_') in a new variable (named here IX1). After a gap of length 1 to 10, the same string IX1 is required again and called back using '?'. This complicated version of mod2() is only shown for the purpose of introducing the notion of instance saving that will be fully used in the next paragraphs. Actually, the various instances of a variable are not necessarily exact copies and this explicit naming process (here _IX1) makes it possible to distinguish one instance from another. Furthermore, such a mechanism allows to save some instance in any part of a model and refer to it elsewhere in the model.

Constraints

Logol modeling is based on a constraint approach. The various constraint types applicable to a model element may be split into two categories, string constraints and structure constraints. String constraints delimit the start (@), the end (@@), the content (?) and the length (#) of admissible strings. Structure constraints include cost constraints ($ for mismatch count, $$ for indel count) and composition constraints (%). These two categories of constraints are written in two separated sets, as in the following model: mod1() ==> X1:{# [START_REF] Eddy | Rnabob: a program to search for rna secondary structure motifs in sequence databases[END_REF][START_REF] Firth | Computational resources for studying recoding[END_REF],@ [START_REF] De Castro | Scanprosite: detection of prosite signature matches and prorule-associated functional and structural residues in proteins[END_REF][START_REF] Jensen | Biogrep: A multi-threaded pattern matcher for large pattern sets[END_REF]}:{% "a":50}. Here mod1 looks for a string whose size is in range [START_REF] Eddy | Rnabob: a program to search for rna secondary structure motifs in sequence databases[END_REF][START_REF] Firth | Computational resources for studying recoding[END_REF], which starts at a position in range [START_REF] De Castro | Scanprosite: detection of prosite signature matches and prorule-associated functional and structural residues in proteins[END_REF][START_REF] Jensen | Biogrep: A multi-threaded pattern matcher for large pattern sets[END_REF] and contains at least 50% of a. For example, the instances of this pattern in the sequence ccaaaacgtacgtttttttcccccc are aaacgt, aaacgta and aacgta (positions start at 0 in a sequence).

• Non exact copies: mismatch and indel cost constraints ($ and $$) Genomic sequences evolve through a duplication process prone to errors or mutations. Elementary variations (on one position) between a model and its instances are taken into account through two dedicated cost counters: the counter of mismatches (i.e. substitutions) and the counter of indels. A mismatch cost constraint is defined by a $[m,n] expression, where m and n are integers. This constraint allows from m to n substitutions. For example, aaaa, acaa and aagt are all instances of the pattern "aaaa":{$[0,2]}. A mismatch constraint can also take the form of a rate: p$[m,n]. Here, m (resp. n) designates the minimum (resp. maximum) allowed percentage of substitutions. For example, aaaa, acaa and aagt are all instances of the pattern "aaaa":{p$[0,50]}. Indels are defined similarly, by setting the indel cost constraints $$[m,n] and p$$[m,n]. Thus, "aaaa":{$$[0,1]} is accepting, among others, the strings aaaa (no indel), aaa (one deletion) or aaaca (one insertion). Here is a new example to further illustrate the concept of instance saving:

X1:{#[5,8],_I1},.*:{#[1,7]}, ?I1:{_I2}:{$[1,1]}, .*:{#[1,7]}, ?I2:{$[1,1]}
This model allows to look for 3 instances of a same string successively deriving from each other (e.g I1 =aaaaa, I2 = aaaca and I3 = agaca). The second pattern, ?I1: {_I2}:{$[1,1]}, reads as follow: the expected string must be similar to the previous I1 string (aaaaa here), apart from 1 mismatch ($ [START_REF] Abouelhoda | Replacing suffix trees with enhanced suffix arrays[END_REF][START_REF] Abouelhoda | Replacing suffix trees with enhanced suffix arrays[END_REF]). The matched string (aaaca) is saved in I2 ({_I2}) for further use ({?I2}). This individualization of instances allows to adjust fine notions of sequence evolution.

• Letter frequencies: composition constraints (%): Some properties like hydrophobic regions in proteins or GC content in RNA correspond to statistical expectations on a particular segment composition rather than the search of a well-defined element. Logol proposes the expression of composition constraints that check the relative frequency of given letters in a sequence. Thus X1:{#[2,43]}:{% "gc":65} describes a segment of length 2 to 43 characters with a GC rate of at least 65%.

Operators

• Negation: Also called negative content constraint, negation can be used in order to exclude some values in a motif. It is denoted by the exclamation mark symbol, !. Thus ("aaa"| "ttt"), !"ga":{#[2,2]} refers to a string made of 5 characters, the first three being 3 a or 3 t, and the next two being anything but the word ga.

• Morphism: A morphism is a function that applies a transformation to a string by substituting letters or substrings. It can be used in direct (+) or reverse (-) direction. Each user can define its own morphisms, but some are already defined. For instance, "wc" transforms a RNA sequence into its complement sequence, applying the Watson-Crick pairing (A-U, G-C). Thus, the pattern +"wc" "acuggc" represents the string "ugaccg" and -"wc" "acuggc" represents the string "gccagu". The morphism -"wc" produces the reverse complement of a string and can be used to describe biological palindromes such as stem-loops. The next example provides a pattern for the recognition of stem-loops whose stem length varies between 5 and 11 and loop size between 1 and 9. Moreover, the Watson-Crick pairing is not required to be perfect: up to 2 substitutions and 1 indel are allowed. STEM1:{# [START_REF] Dsouza | Searching for patterns in genomic data[END_REF][START_REF] Jensen | Biogrep: A multi-threaded pattern matcher for large pattern sets[END_REF],_IS1}, .*:{#[1,9]}, -"wc" ?IS1 :{$[0,2],$$[0,1]} In this description, the content of STEM1 (first strand of the stem) is saved in IS1, (_IS1). The second stem strand is then defined as the exact reverse complement of the previous content (that is -"wc" ?IS1), except for 2 mismatch and 1 indel.

• Repeats: Tandem repeats are frequent genomic structures made of directly adjacent copies of a same entity that may contain only a few letters (microsatellites) or be longer and reach size over 100 nucleic acids (minisatellites). In Logol, such structures can be handled by applying a special constructor, repeat, which manages the characteristics of series of occurrences. Its standard format is: repeat(<entity>,<distance>)+<occurrence number>. For instance, repeat("acgt",[0,3])+ [START_REF] Firth | Computational resources for studying recoding[END_REF]38] states that substring acgt is repeated from 7 to 38 times, using a spacing of at most 3 characters between 2 repeats.

• Views and scope of constraints: Constraints (on the content, the size...) can be set on various parts of a model. They can be imposed to elementary entities like strings or variables as it has been shown previously, or to a set of entities that have themselves individually their own constraints. If the set represents contiguous elements, it is called a view. In Logol syntax, a view is delimited by parentheses. In the following example, the model considers strings built from the concatenation of the instances of 3 variables X1, X2 and X3, each one having length up to 10 characters. A supplementary constraint on the view made of the whole string (X1,X2,X3) requires that its total length is bounded between 8 and 20 characters.

(X1:{#[1,10]}, X2:{#[1,10]}, X3:{#[1,10]} ) : {#[8,20]}
It is also possible to set some constraints on a collection of non-contiguous elements (for instance on the two segments that form the stem of a stem-loop in a RNA structure). Such constraints are set in this case in a specific global module, the control panel. The following example details a stem-loop structure made of two stem elements that have to contain globally at least 30% of C. • Multiple analyses: The coexistence of alternative structures in a same region is certainly amongst the important features of biological sequences. Gene overlapping for instance has been found in all kingdoms of life, including viruses and higher eukaryotes. Logol allows to model such situations by stating alternative models in the grammar top rule (==*> SEQ1). Then, a sequence is accepted only if it contains an instance of each possible alternative. For instance the grammar: mod1().mod2()==*> SEQ1 mod1() ==> "yvcpfdgcnk" mod2() ==> "nklkshil" accepts the sequences containing both the strings yvcpf dgcnk and nklkshil, independently of their positions, being overlapping or not. Parameter-passing is possible between alternative models, e.g. to settle the respective positions of alternative elements. mod1(SAVE1).mod2(SAVE1)==*> SEQ1 mod1(SAVE1) ==> "aata":{_SAVE1},X1:{#[30,30]}:{% "gc":60} mod2(SAVE1) ==> "gggcaa":{@[@SAVE1 -20,@SAVE1 + 20]} The above model is looking for an instance of string aata that is both followed by a GC-reach area and contains a neighboring occurrence of string gggcaa. Indeed, @[@SAVE1 -20,@SAVE1 + 20] constrains the gggcaa string to be located 20 nt before or after the aata string.

Logol Implementation

Input /output Specifications

The Logol software is in charge of matching a Logol pattern against one or more (DNA, RNA or protein) sequences, in order to point out all the pattern occurrences within the sequences.

To this end, it needs two main inputs: a Fasta file with the sequences to be analyzed, and a textual file with the grammatical rules of the pattern. The tool accepts also a configuration file, setting some parsing parameters. This allows, among other, to limit the scope of the search by limiting the maximum number of matches, to choose the indexing tool (Vmatch or Cassiopee, see below), or to detect and filter irrelevant match variants.

The application outputs a zip archive containing one XML file per input sequence. Output files record all the details of the matches, including the matching rule, the location of the match, its length, and the number of substitutions and indel. It also keeps the match information in a tree hierarchy. A repeat, for example, will be decomposed in an array of matches. Thus, it is possible to analyze the result of the global match, but also the details of any element of the grammar. In addition to the XML output, with match details, the analyzer can also output the results in Fasta or GFF3 format. Those formats ease exploiting results within workflows using other tools.

The core application to launch a Logol analysis is a command-line application, available for different operating systems. However, there exists also user interfaces, for more comfort. Among them, the model designer let the user draw the model graphically, and the web application converts it into a Logol grammar.

Sequence Analysis

Pattern matching is performed in two stages. At first the Logol pattern is deciphered by the grammar analyzer, then it is applied on the input sequences by the sequence matcher.

• Parser: The grammar analyzer is a Java program. Its role is to decode the grammar to generate a script used by the sequence matcher, and to launch the calls to the sequence matcher. The generated script is a Prolog file which uses a dedicated library containing Prolog predicates for each kind of grammar element (spacers, repeats, ...). The grammar follows a DSL (Domain Specific Language) analyzed by the Antlr library (http://www.antlr.org/). The Prolog programming language has been chosen for its flexibility and conciseness in expressing parsers, due to its built-in ability for backtracking on partial solutions and natural handling of non determinism. However, the implementation could have been achieved with any other language.

To generate the script, several parsing runs are achieved. The first parsing stage gathers information on each element (expected position range, minimum and maximum size, number of allowed errors, ...).

A second parsing stage tries to solve cases where a variable is used but will be instanciated later in the model, e.g in "acgt",?X1,"cgta":{ X1 }. Indeed, though the grammar itself does not require a left to right reading, we use in practice a left to right parsing of the sequence. To manage such cases, the tool finds the variables in this specific situation and applies a dedicated search technique. It also uses information gathered in the first step to add as many constraints as possible on the variable (length, content, ...) in order to reduce the search space.

The last step, using information from previous stages, generates the Prolog script that will be used by the sequence matcher.

Once the script is generated, the analyzer tries to split the input sequence in smaller parts. Indeed, if grammar analysis or input parameters show match length to be smaller than a fixed integer N, then sequences can be cut in several parts (according to configuration, but at least 2N long). This is used to parallelize the search (multi-threading or using a DRMAA compliant cluster). The analyzer triggers sequence matcher runs on each sequence part and merges the results. In case of multiple input sequences, each sequence analysis will also be parallelized.

If multi-thread is used, the program will limit the number of parallel analysis according to the configuration. If DRMAA is used, the program will also try to use multi-thread on the remote node if sequences can be cut in smaller parts.

• Sequence matcher: The sequence matcher is a Prolog compiled script that loads the script generated by the grammar analyzer. It has been tested with Sicstus Prolog (sicstus.sics.se) and SWI-Prolog (swi-prolog.org). It scans the input sequence with the input script rules, trying to match each rule one by one. When a complete rule is matched, it records the match.

For each rule element, the matcher takes a chunk of the sequence and tries to apply the rule on the chunk. If it matches, it goes to the end location of the match and tries to apply the next rule. The matcher records all the details of the match in an XML format. The matcher will optionally apply a filter to delete redundant matches at different levels since it is possible to get two matches at the same location that only differ by their parsing structures.

In case of spacers in a model, the matcher calls an external program using indexing sequence techniques to directly look for positions of the following words. Two possibilities are offered by Logol to perform indexing: VMatch or Cassiopee.

VMatch [START_REF] Abouelhoda | Replacing suffix trees with enhanced suffix arrays[END_REF] is a suffix array search tool that supports substitution and indel search. An index is created at startup and the matcher calls the VMatch program to search for a pattern in sequences. VMatch is not open source, but is free for academics. The tool is not delivered with the Logol software suite and needs a manual installation. It is efficient for large sequences.

Cassiopee (https://github.com/osallou/cassiopee) is a Ruby tool, developed for Logol, though it can be used independently. It scans the sequence to match a pattern with error support. This tool has been developed to provide a complete open source solution, but it is not as efficient as VMatch for large sequences. Then, VMatch usage is the recommended choice when performance is crucial. The tool selection is made in the configuration file so that it can be adapted for each analysis.

Illustration: Modeling -1 Ribosomal Frameshifting

RNA recoding is a fundamental biological mechanism that cells use to expand the number of proteins assembled from a single DNA code. There are several types of RNA editing modifying the standard translation of a messenger RNA by the ribosome, which seem largely directed by the 3D conformation adopted by the RNA molecule. The -1 programmed ribosomal frameshifting (PRF) is a recoding event which occurs when the ribosome is moving rearward exactly 1 nt on a 'slippery site', X XXY YYZ, where X,Y and Z are nucleotides. The ribosome reads the first X nucleotide two times in this case. Indeed, while standard translation processes codons ABX XXY YYZ CDE ..., the codons processed by PRF are ABX XXX YYY ZCD .... The typical structure promoting a -1 frameshifting event, which we call in the following "PRF pattern", is sketched in Figure 1. It is made sequentially of a start codon, a number of codons, a heptameric slippery site XXXYYYZ placed in -1 phase, a few nucleotide spacer, and a characteristic stable secondary structure. The secondary structure is the obstacle stopping the ribosome during the heptamer translation and triggering a movement one nucleotide backwards that causes reading frameshifting. The secondary structure is usually made of a "H-type pseudo-knot" including two nested stem-loops. A number of tools exists for the detection of putative sites where a -1 frameshifting event might occur [START_REF] Firth | Computational resources for studying recoding[END_REF], but this detection process remains an active research topic since the PRF pattern is not universal (the characteristic features of the heptamer, the spacer, and the secondary structure depend on the organism) and the detection of pseudo-knots is a difficult issue. Many methods proceed by successive filtering steps like KnotInFrame [START_REF] Theis | Knotinframe: prediction of -1 ribosomal frameshift events[END_REF], one of the most advanced tool in this category. KnotInFrame initially detects all heptamers XXXYYYZ, then looks for potential pseudo-knots downstream of this motif, using a dedicated RNA folding procedure.

PRF Logol Model

The complexity of the PRF pattern makes it a good candidate to investigate the expressivity of Logol. In order to elaborate the corresponding model, one needs to use a number of Logol features like multi-analysis, negative content constraints, repeated motifs or search for biological palindromic structures. We present here the most prominent aspects of the model.

• Multi-analysis of two overlapping ORFs: Among the mandatory structural features for the occurrence of a -1 frameshifting event, some are concerning the reading frame setting. The standard translation occurs in a sequence with an open reading frame: a start codon (AUG) followed by a number of non stop codons (triplets), and a stop codon (UGA,UAG ou UAA) terminating the translation. All these codons are in 0 phase. The alternative translation, in case of -1 frameshifting, starts on the same start codon but moves backward one nucleotide on the slippery site, leading to proceed further on triplets in -1 phase until a stop codon is reached, also in -1 phase.

In order to possibly generate a -1 frameshifting event, a RNA sequence should thus contain both an open reading frame in 0 phase (a start followed by a sufficient numbers of codons ending by a stop codon) and at a constrained distance from the start, a series of codons ending by a stop in -1 phase. This check is triggered in Logol by a multiple analysis recognizing alternative patterns, "ORF" and "ORFminus", on the same string with a parameter-passing between the two models in order to share the common start. The ORF model thus contains a repeat that accepts up to 300 non stop codons (a value consistent with the literature), that is: repeat(notstop(),[0,0])+[0,300], where notstop stands for a model built from a view, that accepts a string of length 3 that is not a stop. This is achieved by a negative content constraint on the view.

• Slippery site and spacer: The PRF pattern describes 3 segments: the slippery site, the spacer and the secondary structure. The Logol model for the slippery site respects the consensus: it is an heptameric motif in the form XXXYYYZ, which must be positioned in -1 phase, where X is any nucleotide repeated 3 times, Y is the base A or U repeated 3 times, and Z differs from base G, according to the Logol pattern mod3 below. mod3()==> mod4(), (("aaa")|("uuu")), ! "g":{#[1,1]} mod4()==> (("aaa")|("ccc")|("uuu")|("ggg")) The spacer is straightforwardly described by a gap element (of size less than 10 in this case).

• Pseudo-knots: The most efficient secondary structure for -1 frameshifts is the pseudo-knot of type H (two interwined stem-loops, cf 1), even if it is not the sole existing structure. It is thus the structure that has been modeled here. We first provide a simplified Logol grammar for pseudo-knot structures. In this grammar, STEM15 refers to the first strand (in the 5' direction) of the first stem and -"wc" ?IS15 refers to its 2 nd strand (in the 3' direction), which is its reverse complement up to 4 mismatches. STEM25 and -"wc" ?IS25 refer to the two elements of the2 nd stem. The gaps refer to the loop elements between stems. Our validation process on real data (next paragraph) resulted in a significant refinement of this model (cf Figure 2). The final model makes use of a great variety of Logol language elements. Among new elements, the model integrates the count of the GC ratio in stems, the separate treatment of nucleotides at the end of stems in order to forbid mismatches at these positions, or the possibility of non-canonical Wobble pairing (G-U), called wcw here, at particular stem positions [START_REF] Rocheteau | Recherche d'éléments structurés dans les génomes par modèles logiques[END_REF]. An excerpt from the Logol grammar dedicated to the first stem follows, the whole model being presented in figure 2: • A first model validation: In order to test and refine our -1 frameshift model, we have elaborated a sequence test set [START_REF] Rocheteau | Recherche d'éléments structurés dans les génomes par modèles logiques[END_REF] around sequences known to produce -1 frameshift events. Thirty proven sequences ("validated -1 frameshift") from the reference base Recode2 (recode.genetics.utah.edu) have been completed by random sequences obtained by shuffling 100 times each reference sequence, using Shuffleseq (emboss.bioinformatics.nl/cgi-bin/emboss/). This procedure keeps the sequence lengths and nucleotidic ratios of the reference. The final set thus contains 30 "positive" sequences (in which one expects to find the PRF pattern at the right position) and 300 "negative" sequences (in which one expects a minimum number of PRF pattern occurrences). This validation work [START_REF] Rocheteau | Recherche d'éléments structurés dans les génomes par modèles logiques[END_REF] led us to perform comparisons between the Logol prediction of pseudo-knots and those made by "DotKnot", a pseudo-knot prediction software using RNAfold to compute the probability of each folding structure (http://dotknot.csse.uwa.edu.au). It brought about some significant changes in our initial model: parameter tuning has concerned wooble pairing, GC ratio and the mismatch rate in the stems. Ultimately, the Logol model finds about 100 matches per sequence on the Recode2 reference set, among them being in most cases the desired match. It is possible to compute a posteriori a quality score for each stem and sort accordingly the matches [START_REF] Rocheteau | Recherche d'éléments structurés dans les génomes par modèles logiques[END_REF]. This procedure leads to the right frameshift area prediction in 20 cases over the 30. The score is based on a pairing cost function proposed by J.P. Forest [START_REF] Forest | Modélisation et détection automatique de sites de décalage de cadre en -1 dans les génomes eucaryotes[END_REF] for the stem: {GCpairing=+3, AUpairing=+2, GUpairing=+1, mismatch=-2}.

• Runtime To give an idea of the performances, parsing the largest reference sequence (30Kb) with the final Logol model takes 1'30s (on a PC Intel X5550, 144Go RAM), while the KnotinFrame answer is immediate on such a sequence. The complete analysis of the Bacillus subtilis genomic sequence (str168NC 000964.3 , 4.2 Mbp) produces 7000 matches in 2 hours. KnotinFrame web site does not accept such a large sequence.

Conclusion

The Logol pattern matching tool has been conceived to allow the modeling and search of realistic structures in biological sequences. It has been designed to be expressive but also evolutive, in order to ease the introduction of new features. The fact that the language has proved fairly well suited to model the complex pattern of ribosomal frameshift, whereas it was not designed for this task, seems an encouraging sign on the genericity of the language elements. The tool is operational and available on the GenOuest bioinformatics platform, under CeCILL license. Although efforts have been made to offer a wide access to Logol functionalities through either command-line or graphical interface inputs, we welcome any user feed-back to increase its ergonomic features and actual range of applicability.

  controls:{ % "c"[mod1.ISTEM1,mod1.ISTEM2]>=30 } mod1()==> STEM1:{#[2,18],_ISTEM1},.*:{#[1,10]}, -"wc" ?ISTEM1:{_ISTEM2} mod1()==*>SEQ1

Fig. 1 .

 1 Fig. 1. "PRF pattern"= typical structure promoting a -1 frameshifting event.

  STEM15:{#[4,16],_IS15},.*:{#[1,5]}, STEM25:{#[3,8],_IS25},.*:{#[0,4]}, -"wc" ?IS15 :{$[0,4]},.*:{#[4,40]},-"wc" ?IS25 :{$[0,2]}

  // 50% of GC pairing in Stem1 => 25% of C in [Stem1.5' + Stem1.3'] controls:{ % "c"[mod5.IA5,mod5.IS15,mod5.IZ5,mod5.IZ3,mod5.S13,mod5.IA3]>=25} mod5()==> (A5:{#[1,1],_IA5},S15:{#[2,14],_IS15},Z5:{#[1,1],_IZ5}):{%"gc":50}, LOOP1:{#[1,5]}, ... stuff deleted .... -"wcw" ?IZ5:{_IZ3}, -"wc" ?IS15:{_S13}:{p$[0,34]}, -"wcw" ?IA5:{_IA3}, ...

Fig. 2 .

 2 Fig. 2. Global overview of the final Logol model for the PRF pseudo-knot


Pierre-Anthony Grenier
email: pierre-anthony.grenier@ensicaen.fr

Luc Brun
email: luc.brun@ensicaen.fr

Didier Villemin
email: didier.villemin@ensicaen.fr

Incorporating molecule's stereisomerism within the machine learning framework

An important field of chemoinformatics consists in the prediction of molecule's properties, and within this field, graph kernels constitute a powerful framework thanks to their ability to combine a natural encoding of molecules by graphs, with classical statistical tools. Unfortunately some molecules encoded by a same graph and differing only by the three dimensional orientation of their atoms in space have different properties. Such molecules are called stereoisomers. These latter properties can not be predicted by usual graph methods which do not encode stereoisomerism. In this paper we propose to encode the stereoisomerism property of each atom of a molecule by a local subgraph. A kernel between bags of such subgraphs provides a similarity measure incorporating stereoisomerism properties. We then propose two extensions of this kernel incorporating in each sub graph information about its surroundings.

Introduction

A molecular graph is a graph G = (V, E, µ, ν), where each node v ∈ V encodes an atom and each edge e ∈ E a bond between two atoms. The labelling functions µ and ν associate to each vertex and each edge a label encoding respectively the nature of the atom (carbon, oxygen,. . . ) and the type of the bond (single, double, triple or aromatic). However, those graphs have a limitation: they do not encode the spatial configuration of atoms. Some molecules, called stereoisomers, are associated to a same molecular graph but differ by the relative positioning of their atoms.

Most of stereoisomers are characterized by the three dimensional orientation of the direct neighbors of a single atom or two connected atoms. We can consider for example, a carbon atom, with four neighbors, each of them located on a summit of a tetrahedron. If we permute two of the atoms, we obtain a different spatial configuration and hence an alternative stereoisomer (Figure 1(a)). An atom is called a stereocenter if a permutation of two atoms belonging to its neighborhood produces a different stereoisomer. We should stress here that, to a large extend, stereoisomerism is independent of a particular embedding of a molecule. Indeed, in Figure 1(a), any particular embedding keeping the same relative positioning of atoms H, Cl, Br and F according to the central carbon atom C, would correspond to a same stereoisomer. In the same way, two connected atoms form a stereocenter if a permutation of the positions of two atoms belonging to the union of their neighborhoods produces a different stereoisomer (Figure 1(b)). According to chemical experts [START_REF] Jacques | Enantiomers, racemates, and resolutions[END_REF], within molecules currently used in chemistry, 98% of stereocenters correspond either to carbons with four neighbors, called asymmetric carbon (Figure 1(a)) or to couples of two carbons adjacent through a double bond (Figure 1(b)). We thus restrict the present paper to such cases.

Graph kernels [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF][START_REF] Gaüzère | Two New Graphs Kernels in Chemoinformatics[END_REF], provide a measure of similarity between graphs. Under the assumption that a kernel k is symmetric and definite positive, the value k(G, G ′ ), where G and G ′ encode two graphs, corresponds to a scalar product between two vectors Ψ (G) and Ψ (G ′ ) in an Hilbert space. This latter property allows us to combine graph kernels with usual machine learning methods such as SVM or kernel ridge regression by using the well known kernel trick, which consists in replacing the scalar product between Ψ (G) and Ψ (G ′ ) by k(G, G ′ ) in these algorithms.

Up to now, only few methods have attempted to incorporate stereoisomerism within the graph kernel framework. Brown et al. [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF] have proposed to incorporate this information through an extension of the tree-pattern kernel [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF]. One drawback of this method is that, patterns which encode stereo information, and patterns which do not, are combined without any weighting in the final kernel value. So for a property only related to stereoisomerism, patterns that do not encode stereo information may be assimilated to noise which deteriorates the prediction. Grenier et al. [START_REF] Grenier | Treelet kernel incorporating chiral information[END_REF] have introduced the minimal subtree which characterizes a stereocenter within an acyclic molecule. They also proposed a kernel based on this minimal subtree, which takes into account stereoisomerism. This kernel is however restricted to acyclic graphs.

Based on [START_REF] Grenier | Treelet kernel incorporating chiral information[END_REF], we present in Section 2 an encoding of molecules distinguishing stereoisomers. Section 3 present the construction of a subgraph, which allows to characterizes locally a stereocenter. Then in Section 4, we use this subgraph to propose new graph kernels valid for cyclic as well as acyclic molecules, thus overcoming the main limitation of [START_REF] Grenier | Treelet kernel incorporating chiral information[END_REF]. We finally present in Section 5 results obtained using those kernels and compare these results with state of the art methods.

Ordered Graph and Stereo Vertices

The spatial configuration of the neighbors of each atom may be encoded through an ordering of its neighborhood. For example, considering the left part of Figure 1(a), and looking at the central carbon from the hydrogen atom (H), the sequence of remaining neighbors of the carbon: Cl, Br and F may be considered as lying on a plane and are encountered clockwise. Thus, this spatial configuration is encoded by the sequence H, Cl, Br, F and the sequence H, Br, Cl, F encodes the second configuration.

In order to encode this information in a graph, we introduce the notion of ordered graph. An ordered graph G = (V, E, µ, ν, ord) is a molecular graph G m = (V, E, µ, ν) together with a function ord : V → V * which maps each vertex to an ordered list of its neighbors. Two ordered graphs G and

G ′ are isomorphic (G ≃ o G ′ ) if there exists an isomorphism f between their respective molecular graphs G m and G ′ m such that ord ′ (f (v)) = (f (v 1 ) . . . f (v n )) with ord(v) = (v 1 . . . v n ) (where N (v) = {v 1 , . . . , v n } denotes the neighborhood of v).
In this case f is called an ordered isomorphism between G and G ′ .

However, different ordered graphs may encode a same molecule. In the example of the left part of Figure 1(a), if we look to the central carbon from a different neighbor, we can obtain a different sequence, for example F, Br, Cl, H, that represents the same configuration but now considered from the atom F. We thus have to define an equivalence relationship between ordered graphs, such that two ordered graphs are equivalent if they represent a same configuration.

To do so, we introduce the notion of re-ordering function σ, which associates to each vertex v ∈ V of degree n a permutation σ(v) on {1, . . . , n}, which allows to re-order its neighborhood. The graph with re-ordered neighborhoods σ(G) is obtained by mapping for each vertex v its order ord(v) = v 1 . . . .v n onto the sequence v σ(v)(1) . . . .v σ(v)(n) where σ(v) is the permutation applied on v.

In order to define a permutation σ(v) for each vertex of a graph, we first introduce the notion of potential asymmetric carbon which corresponds to a carbon with four neighbors. Such a vertex corresponds to a stereocenter if one permutation of two of its neighbors provides a different stereoisomer (Section 1). Permutations associated to a potential asymmetric carbon correspond to all even permutations of its four neighbors [START_REF] Petitjean | Chirality in metric spaces[END_REF]. For a double bond between two carbons, permutations associated to each carbon of the double bound must have a same parity. Finally, for any vertex which does not correspond to a potential asymmetric carbon nor to a carbon of a double bond, we do not search to characterize its spatial configuration. So these vertices are associated to all possible permutations of their neighbors.

The set of re-ordering functions, transforming an ordered graph into another one representing the same configuration is called a valid family of re-ordering functions Σ [START_REF] Grenier | Incorporating stereo information within the graph kernel framework[END_REF]. We say that it exists an equivalent ordered isomorphism f between G and

G ′ according to Σ if it exists σ ∈ Σ such that f is an ordered isomorphism between σ(G) and G ′ (σ(G) ≃ o G ′ ).
The equivalent order relationship defines an equivalence relationship [START_REF] Grenier | Incorporating stereo information within the graph kernel framework[END_REF] and two different stereoisomers are encoded by non equivalent ordered graphs. We denote by IsomEqOrd(G, G ′ ) the set of equivalent ordered isomorphism between G and G ′ .

Potential asymmetric carbons, and double bonds between carbons, are not necessarily stereocenters. For example if the label of vertex Br of Figure 1(a) is replaced by Cl, both left and right molecules of Figure 1(a) would be identical. In the same way, if the label of the vertex F in Figure 1(b) is replaced by Cl, the left and right molecules of this figure also become identical. For those cases, any permutation in the ordered list of the carbons would lead to an equivalent ordered graph. We thus define a stereo vertex as a vertex for which any permutation of two of its neighbors produces a non-equivalent ordered graph:

Definition 1 (Stereo vertex). Let G = (V, E, µ, ν, ord) be an ordered graph. A vertex v ∈ V is called a stereo vertex iff: ∀(i, j) ∈ {1, . . . , |N (v)|} 2 , i = j, ∄f ∈ IsomEqOrd(G, τ v i,j (G)) with f (v) = v.
(1) where τ v i,j (G) corresponds to an ordered graph deduced from G by permuting nodes of index i and j in ord(v).

Minimal Stereo SubGraph

Definition 1 is based on the whole graph G to test if a vertex v is a stereo vertex. However, given a stereo vertex s, one can observe that on some configurations, the removal of some vertices far from s should not change its stereo property. In order to obtain a more local characterization of a stereo vertex, we should thus determine a vertex induced subgraph H of G, including s, large enough to characterize the stereo property of s (i.e. ∀(i, j) ∈ {1, . . . , |N (s)|} 2 , i = j, ∄f ∈ IsomEqOrd(H, τ s i,j (H)) with f (s) = s), but sufficiently small to encode only the relevant information characterizing the stereo vertex s. Such a subgraph is called a minimal stereo subgraph of s.

We now present an heuristic, used to compute a minimal stereo subgraph of a stereo vertex. We focus our attention on asymmetric carbons. Let H be a subgraph of G containing a stereo vertex s corresponding to an asymmetric carbon. We say that the stereo property of s is not captured by H if (Definition 1):

∃(i, j) ∈ {1, . . . , |N (s)|} 2 , i = j, ∃f ∈ IsomEqOrd(H, τ s i,j (H)) with f (s) = s (2)
To define a minimal stereo subgraph of s, we consider a finite sequence (H k s ) n k=1 of vertex induced subgraphs of G. The first element of this sequence H 1 s is the smaller vertex induced subgraph for which we can test (2) :

V (H 1 s ) = {s} ∪ N (s)
where V (H 1 s ) and N (s) denote respectively the set of vertices of H 1 s and the set of neighbors of s in G.

If the current vertex induced subgraph H k s does not capture the stereo property of s, we know by [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF], that it exists some isomorphisms f of equivalent ordered graphs between H k s and τ s i,j (H k s ) with i = j and f (s) = s. Let us consider such an isomorphism f . By definition of equivalent ordered isomorphism, it exists σ ∈ Σ such that f is an ordered isomorphism between H k s and σ τ s i,j (H k s ) . By definition of ordered isomorphisms, and since f (s) = s, we have:

∀l ∈ {1, . . . , |N (s)|}, f (v l ) = v σ(s)•τ s i,j (l) . with ord(s) = v 1 , . . . , v n .
As σ(s) is an even permutation, σ(s) • τ s i,j is an odd one. Hence it exists l in {1, . . . , |N (s)|} such that l = σ(s) • τ s i,j (l) and we have f

(v l ) = v l and f (2) (v l ) = f (v l ).
In other words, any equivalent ordered isomorphism corresponding to equation (2) maps at least two vertices in the neighborhood of s in H k s onto a different vertex in the same neighborhood. Let us denote by E k f the set of vertices of H k s connected to s by a path whose all vertices are mapped onto other vertices by f :

E k f = {v ∈ V (H k s ) | ∃c = (v 0 , . . . , v q ) ∈ H k s with v 0 = s and v q = v s.t. ∀r ∈ {1, . . . , q}, f (v r ) = v r } (3)
For any equivalent ordered isomorphism f satisfying (2), the set E k f is not empty since it contains at least 2 vertices. A vertex v belongs to E k f if neither its label nor its neighborhood in H k s allows to differentiate it from f (v). The basic idea of our algorithm consists in enforcing constraints on each v ∈ E k f at iteration k + 1 by adding to H k s the neighborhood in G of all vertices belonging to E k f . This last set is denoted by N (E k f ). The set of vertices of the vertex induced subgraph H k+1 s is thus defined by:

V (H k+1 s ) = V (H k s ) ∪ f ∈F k s N (E k f ) (4)
where F k s denotes all equivalent ordered isomorphisms satisfying (2). Since f ∈ F k s implies that E k f is not empty, adding iteratively constraints on the existence of vertices in E k f removes f from F k s . The algorithm stops when the set F k s becomes empty. Note that such a condition must be satisfied since s is a stereocenter and hence the whole molecule does not satisfy [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF].

The intermediate vertex induced subgraphs found by our algorithm are illustrated in Fig. 2. Note that at iteration 2, it exists an equivalent ordered isomorphism f ∈ F 2 C mapping the path CCO (bottom right of the figure) onto the same path located on the top right part of Fig 2 . In this case E 2 f contains the three carbons of these two paths and both oxygen atoms. The oxygen atoms belong to E 2 f since their neighborhoods in H 2 C does not allow to differentiate them (Fig. 2). At iteration 3, the neighborhood in G of these oxygen atoms are added to H 3 C , hence adding N and Br which allow to differentiate both paths and thus removes the equivalent ordered isomorphism f from F 3 C .

Stereo Kernel

Given an ordered G, we can associate a minimal stereo subgraph to each of its stereocenter. A same stereo subgraph may be present more than once in a given molecule, we thus need to associate a unique code to each such subgraph in order to enumerate efficiently the eventual multiple occurences of a stereo subgraph within a molecule. To do so, we use [START_REF] Wipke | Stereochemically unique naming algorithm[END_REF], which associates to each molecule a unique code which allows to test the existence of an equivalent ordered isomorphism between two stereo subgraphs, unlike [START_REF] Bonnici | A subgraph isomorphism algorithm and its application to biochemical data[END_REF] which allows to find efficiently all isomorphisms between two graphs. We can thus compute the set of minimal stereo subgraphs H(G) together with the spectrum spec(G) = (f H (G)) H∈H(G) which encodes the frequency f H (G) of each H ∈ H(G). The set H(G) and the spectrum spec provide a characterisation of each stereo center of G and hence describe the stereoisomerism of G.

The comparison of the spectrum of two ordered graphs, is then used to define a kernel between two molecules taking into account the stereoisomerism:

k(G, G ′ ) = H∈H(G)∩H(G ′ ) K(f H (G), f H (G ′ )). ( 5
)
where K denotes a kernel between real values (e.g. Gaussian, intersection or polynomial). The choice of a particular kernel, together with its parameters is performed through cross-validation.

Augmented labels

Cycles are important sub-parts of molecules, and thus two atoms with identical label could have different influence if one of them is included in a cycle. Thus it can be useful, during the computation of a minimal stereo subgraph, to consider two atoms with a same label, but not included in a same number of cycles, as different.

To do so, we first compute the set of relevant cycles of each ordered graph. Relevant cycles are defined as cycles of a graph which can not be deduced from shorter cycles [START_REF] Vismara | Union of all the minimum cycle bases of a graph[END_REF]. We can associate to each vertex v, the number n v of relevant cycles to which it belongs. Then, for an ordered graph G, this information is added to the label of each of its vertex v (µ A (v) = µ(v).n v ) to obtain a new ordered graph G A . The method described in Section 3 is then applied on the ordered graph G A . We thus obtain a different set of minimal stereo subgraphs H(G A ), composed of smaller stereo subgraphs where nodes encode a more global information. As in Section 4.1 we define from this set of subgraphs a spectrum encoding the frequency of each subgraph, and define a kernel between graphs by comparing those spectrum:

k(G, G ′ ) = H∈H(G A )∩H(G ′ A ) K(f H (G), f H (G ′ )). ( 6
)

Expanded Subgraphs

Equations 5 and 6 allow to compare two molecules through the distribution of their stereo subgraphs. However those kernels are based on a binary similarity measure between configurations: either two stereo subgraphs are different, and thus the configurations encoded by those subgraphs are dissimilar, or the subgraphs are identical and thus the configurations are similar. By adding information about the adjacency relationships between these stereo subgraphs and the remaining part of the molecule, we can obtain a finer measure of similarity between configurations around stereocenters.

To take into account the adjacency relationships between a stereo subgraph H s and its surrounding, we consider larger vertex induced subgraphs than H s . Let H be a subgraphs of G, the neighborhood N (H) of H is the set of vertices of G -H which are neighbors of a vertex of H:

N (H) = {v ∈ V (G) -V (H) | ∃(u, v) ∈ E s.t u ∈ V (H)}
The set of vertex induced subgraphs obtained by adding k of its neighbors to H s can be used to construct a kernel between graphs. We can also consider subgraphs where vertex located farther from the stereo subgraph than its direct neighborhood are added. However we have to limit the number of vertices we add, in order to keep a local information. Moreover the number of subgraphs increases quickly with the number of added vertices. Indeed, C k N subgraphs can be constructed by adding k vertices of N (H s ) to H s , with N = |N (H s )|. We thus, have to determine a number of vertex to add, which is large enough to characterize the adjacency relationships between a stereo subgraph and the remaining part of a molecule, but also sufficiently small to keep a local information. In our experiment, we have considered subgraphs obtained by adding up to three different neighbors of H s and those obtained by adding one neighbor v of H s , and one neighbor of v not included in the neighborhood of H s .

For each minimal stereo subgraph, we have a set of subgraphs which encodes its adjacency relationships with other parts of the molecule. As in section 4.1, we associate to those subgraphs a unique code. By adding those subgraphs to the set of minimal stereo subgraphs, we obtain a new set of subgraphs H E (G), for which we can associate a spectrum which encodes the frequency of each subgraphs H ∈ H E (G). We thus define a kernel between ordered graphs, which takes into account stereoisomerism, and the adjacency relationships of each stereo subgraphs with its surrounding:

k(G, G ′ ) = H∈H E (G)∩H E (G ′ ) K(f H (G), f H (G ′ )). ( 7
)

Experiments

Our first experiment is based on a dataset composed of all the stereoisomers of the perindoprilate [START_REF] Castillo-Garit | Atombased stochastic and non-stochastic 3d-chiral bilinear indices and their applications to central chirality codification[END_REF]. As this molecule has 5 stereocenters, the dataset is composed of 2 5 = 32 molecules. In this dataset, we try to predict if a molecule inhibit the angiotensin-converting enzyme (ACE). The dataset is split into a training set of 23 molecules, and a test set of 9 molecules, as in [START_REF] Castillo-Garit | Atombased stochastic and non-stochastic 3d-chiral bilinear indices and their applications to central chirality codification[END_REF].

Table 1 shows the results obtained by our kernels and the adaptation of the tree pattern kernel to stereoisomerism [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF]. All these kernels are combined with the standard SVM method [START_REF] Cortes | Support-vector networks[END_REF] to classify molecules. As all molecules in the dataset are stereoisomers of each other, methods which do not include stereoisomerism information [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF][START_REF] Gaüzère | Two New Graphs Kernels in Chemoinformatics[END_REF] are unable to differentiate any molecule of this dataset and are consequently unable to predict the considered property. Moreover, information not related to stereoisomerism included in kernel [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF] consists of the same patterns for all molecules. This leads to add a constant shift to all values of the kernel and hence does not deteriorate the prediction for this dataset. Two stereocenters of the molecules of the dataset have a same minimal stereo subgraph, however one of them contains vertices belonging to a cycle. The stereo kernel (Section 4.1) is not able to distinguish these stereocenters, and hence misclassified molecules containing these stereocenters. By using augmented labels (Section 4.2), these two stereocenters are distinguished, and this distinction allows us to reach a prediction accuracy of 100%. The expanded subgraph (Section 4.3) may also help to differentiate the two stereocenters, but for this dataset, one molecule of the testset remains misclassified due to the same stereocenters which are not sufficiently discriminated by this kernel. [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF] 100 100

The second dataset is a dataset of synthetic vitamin D derivatives, used in [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF]. This dataset is composed of 69 molecules containing cycles, with an average of 9 stereocenters per molecule. This dataset is associated to a regression problem, which consists in predicting the biological activities of each molecules. Each kernel is test by using it with the standard SVM regression method [START_REF] Drucker | Support vector regression machines[END_REF].

After normalizing the values of the dataset, the standard deviation of the biological activities is equal to 0.258. To choose the different parameters and estimate the performance of each kernel on this dataset we use a nested crossvalidation. The outer cross-validation is a leave-one-out procedure, used to compute an error for each molecule of the dataset. For each fold, we use another leave-one-out procedure on the remaining molecules, to compute a validation error. Parameters which provide the lowest root mean squared error on the validation are selected. We obtain for each molecule an error, and report in Table 2, the mean of this distribution of errors together with the confidence interval at 95% of this distribution.

Greatest errors in Table 2 are obtained by methods [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF][START_REF] Gaüzère | Two New Graphs Kernels in Chemoinformatics[END_REF] which do not include stereo information. The adaptation of the tree pattern kernel to stereoisomerism [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF] improves the results over the two previous methods hence showing the insight of adding stereoisomerism information. Our kernel with no extensions obtain results not as good as [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF]. For this experiment the modification of label to incorporate information about cycles, decrease our results. However, the addition of information about relationships between minimal stereo subgraphs, and remaining part of the molecules, allow us to obtain better results than [START_REF] Brown | Compound analysis via graph kernels incorporating chirality[END_REF]. In this case the best results are obtained by considering only subgraphs including one neighbor of H s . Our methods use a subgraph isomorphism algorithm, but the minimal stereo subgraph are small and thus we have small execution times as we can see in Table 2.

Conclusion

The study of stereoisomers constitutes an important subfield of chemistry and thus a major challenge in chemoinformatics. We have proposed in this paper, a graph kernel based on an explicit enumeration of all the stereo subgraphs of a molecule. Each stereo subgraph is associated to a stereo vertex and encodes the part of the graph which provides the stereo property to this vertex. Based on the notion of stereo subgraphs we propose to describe a molecule by its bag of stereo subgraphs. The similarity between two molecules is then encoded through a graph kernel based on the similarity of both bags. Moreover we propose two extensions of this kernel. One extension consists in adding to the labels of the graphs information about cycles of molecules. The second one consists in considering larger subgraphs encoding relationships between each stereo subgraph and the remaining part of the molecule. Experiments related to stereoisomerism properties demonstrate the relevance of our approach and of its extensions.

  (a) Two different spatial configurations of the neighbors of a carbon (b) Two different spatial configurations of two carbons linked by a double bond.

Fig. 1 .

 1 Fig. 1. Two types of stereocenters.

Fig. 2 . 3 k=1 4

 234 Fig. 2. An asymmetric carbon and its associated sequence (H k C ) 3 k=1

Table 1 .

 1 Classification of the ACE inhibitory activity of perindopirilates stereoisomers

	Accuracy Accuracy

Table 2 .

 2 Prediction of the biological activity of synthetic vitamin D derivatives.

	Mean Confidence Gram's matrices

Acknowledgements

The authors wish to thank the association CRIHAN for their computing resources.