Thibault Broux

Tahya Bamine

Laura Simonelli

Lorenzo Stievano

François Fauth

Michel Ménétrier

Dany Carlier

Christian Masquelier

Laurence Croguennec
email: laurence.croguennec@icmcb.cnrs.fr

V IV Disproportionation Upon Sodium Extraction From Na 3 V 2 (PO 4 ) 2 F 3 Observed by Operando X-ray Absorption Spectroscopy and Solid State NMR

Among the series of polyanionic positive electrodes for Na-ion batteries having the general formula Na 3 V 2 (PO 4 ) 2 F 3-y O y (0 ≤ y ≤ 2), the composition Na 3 V 2 (PO 4 ) 2 F 3 (y = 0) has the highest theoretical energy that offers competitive electrochemical performances compared to sodium transition metal oxides. Recently the structural phase diagram from Na 3 V 2 (PO 4 ) 2 F 3 to Na 1 V 2 (PO 4 ) 2 F 3 has been thoroughly investigated by operando synchrotron X-ray diffraction revealing an unexpected structural feature for the end member composition. In fact the crystal structure of Na 1 V 2 (PO 4 ) 2 F 3 has two very different vanadium environments within each bi-octahedron that suggests a charge disproportionation of two V IV into V III and V V . This work shows an operando X-ray absorption spectroscopy at vanadium K-edge during the electrochemical extraction of Na + in order to monitor the redox processes involved in this compound. The large data set provided by this experiment has been processed by the principal component analysis combined with multivariate curve resolution. The results suggest that the bioctahedra have to be considered as the basic structural unit. The peculiar geometry of this material combined with the mixed vanadium valence, directly investigated here along the reaction, seems to allow original electronic configurations. In particular the two vanadium sites into the basic bioctahedra unit evolves from V III -V III to V III -V IV and to a final V III -V IV configuration. These observations are completed with 51 V NMR sensitive to diamagnetic V V .

Introduction

The increasing need for electrochemical energy storage, either for portable electronic devices or for larger-scale applications such as hybrid electric vehicles or static renewable energy storage systems, causes renewed interest in alternatives to Li-ion batteries in order to overcome drawbacks associated to the availability and prize of lithium resources. [START_REF] Tarascon | Is lithium the new gold?[END_REF] In this context Na-ion batteries is an emerging field owing to lower price and very large earth-abundance of sodium. [START_REF] Palomares | Naion batteries, recent advances and present challenges to become low cost energy storage systems[END_REF] In terms of energy density Na-ion batteries hardly competes with Li-ion ones due to the intrinsic properties of Na: a less negative standard reduction potential (-2.7 V vs SHE for the Na + aq /Na against -3.04 V for the Li + aq /Li one) and a higher molecular weight. Thus Na-ion technology can be addressed for targeted applications such as domestic energy storage and load leveling applications. Among several polyanionic-based electrode materials, [START_REF] Jian | Superior electrochemical performance and storage mechanism of Na 3 V 2 (PO 4 ) 3 cathode for room-temperature sodium-ion batteries[END_REF][START_REF] Sauvage | Crystal structure and electrochemical properties vs. Na + of the sodium fluorophosphate Na 1.5 VOPO 4 F 0.5[END_REF][START_REF] Shakoor | A combined first principles and experimental study on Na 3 V 2 (PO 4 ) 2 F 3 for rechargeable Na batteries[END_REF][START_REF] Liu | Preparation of NaV 1x Al x PO 4 F cathode materials for application of sodium-ion battery[END_REF][START_REF] Chihara | Cathode properties of Na 3 M 2 (PO 4 ) 2 F 3 [M= Ti, Fe, V] for sodium-ion batteries[END_REF][START_REF] Kim | Na 2 FeP 2 O 7 as a promising iron-based pyrophosphate cathode for sodium rechargeable batteries: a combined experimental and theoretical study[END_REF][START_REF] Masquelier | Polyanionic (phosphates, silicates, sulfates) frameworks as electrode materials for rechargeable Li (or Na) batteries[END_REF] one of the most promising family turns out to be the sodium-vanadium fluorophosphate compounds Na 3 V 2 (PO 4 ) 2 F 3-y O y where y can vary from 0 to 2. [START_REF] Serras | High voltage cathode materials for Na-ion batteries of general formula Na 3 V 2 O 2x (PO 4 ) 2 F 3-2x[END_REF][START_REF] Serras | Electrochemical performance of mixed valence Na 3 V 2 O 2x (PO 4 ) 2 F 3-2x /C as cathode for sodium-ion batteries[END_REF][START_REF] Park | A new high-energy cathode for a Na-ion battery with ultrahigh stability[END_REF][START_REF] Park | A Family of high-performance cathode materials for Na-ion batteries, Na 3 (VO 1-x PO 4 ) 2 F 1+2x (0≤ x≤ 1): combined first-principles and experimental study[END_REF][START_REF] Bianchini | Comprehensive investigation of the Na 3 V 2 (PO 4 ) 2 F 3 -NaV 2[END_REF][START_REF] Bianchini | Na 3 V 2 (PO 4 ) 2 F 3 revisited: a high-resolution diffraction study[END_REF][START_REF] Tsirlin | Phase separation and frustrated square lattice magnetism of Na 1.5 VOPO 4 F 0.5[END_REF][START_REF] Gover | The electrochemical insertion properties of sodium vanadium fluorophosphate, Na 3 V 2 (PO 4 ) 2 F 3[END_REF][START_REF] Song | Multifunctional dual Na 3 V 2 (PO 4 ) 2 F 3 cathode for both lithium-ion and sodium-ion batteries[END_REF][START_REF] Broux | Strong impact of the oxygen content in Na 3 V 2 (PO 4 ) 2 F 3y O y (0≤y≤0.5) on its structural and electrochemical properties[END_REF][START_REF] Ponrouch | Towards high energy density sodium ion batteries through electrolyte optimization[END_REF] This whole range of compositions from V III -rich (y = 0) to V IV -rich (y = 2), can be oxidized to V IV -rich and V V -rich compounds respectively by the deintercalation of two sodium ions. For instance in Na 3 V 2 (PO 4 ) 2 F 3 (y = 0) the extraction of 2 Na + ions has been experimentally demonstrated with two main voltage-composition plateaus at around 3.7 and 4.2 V vs Na + /Na giving a theoretical energy density of 507 Wh/kg [START_REF] Gover | The electrochemical insertion properties of sodium vanadium fluorophosphate, Na 3 V 2 (PO 4 ) 2 F 3[END_REF][START_REF] Barker | Hybrid-ion a lithium-ion cell based on a sodium insertion material[END_REF][START_REF] Barker | Li 4/3 Ti 5/3 O 4 ‖ Na 3 V 2 (PO 4 ) 2 F 3 : An example of a hybrid-ion cell Using a Non-graphitic Anode[END_REF] (128 Ah/kg at an average potential of 3.95 V), competitive with that delivered by LiFePO 4 in Li-ion batteries. As recently depicted by Bianchini et al. [START_REF] Bianchini | Na 3 V 2 (PO 4 ) 2 F 3 revisited: a high-resolution diffraction study[END_REF] the pristine Na 3 V III 2 (PO 4 ) 2 F 3 crystalizes in the Amam space group at room temperature. This structure consists of a tridimensional framework of V 2 O 8 F 3 bi-octahedra which are connected by PO 4 tetrahedra that generates large tunnels where Na + ions are mobile upon extraction/insertion reactions.

The operando structural evolution has been extensively characterized by Bianchini et al. [START_REF] Bianchini | Comprehensive investigation of the Na 3 V 2 (PO 4 ) 2 F 3 -NaV 2[END_REF] who described the phase diagram of Na 3 V III 2 (PO 4 ) 2 F 3 -Na 1 V IV 2 (PO 4 ) 2 F 3 as a complex succession of biphasic reactions and solid solution domains. One important and also unexpected feature raised by this study was the assumed charge disproportionation of 2 V IV into 1 V V and 1 V III observed at the end of the charge (i.e. after the extraction of 2 Na + ions) for the composition Na 1 V IV 2 (PO 4 ) 2 F 3 (Na 1 VPF in the following). Such a charge disproportion was for instance also observed in Li x V 2 O 5, [START_REF] Rao | Structural and electrochemical behaviour of sputtered vanadium oxide films: oxygen non-stoichiometry and lithium ion sequestration[END_REF] whereas more generally charge separation is commonly observed in inorganic vanadium-rich compounds. [START_REF] Schindler | Crystal chemical aspects of vanadium: polyhedral geometries, characteristic bond valences, and polymerization of (VO n ) polyhedra[END_REF][START_REF] Seo | Charge ordering and spin gap in NaV 2 O 5[END_REF][START_REF] Nakao | X-Ray Anomalous Scattering Study of a charge-ordered state in NaV 2 O 5[END_REF][START_REF] Ballhausen | The electronic structure of the vanadyl ion[END_REF][START_REF] Ohama | Mixed valency and charge ordering in NaV 2 O 5[END_REF] Indeed the environment stabilized for vanadium, especially in oxides, changes drastically from an octahedron to a square pyramid with the oxidation state.

The fully charged composition Na 1 V 2 (PO 4 ) 2 F 3 is described using the Cmc2 1 space group which has two distinct vanadium crystallographic sites associated to V III and V V instead of only one for the pristine material associated to V III . This hypothesis of charge disproportionation is based on bond valence sum calculations which reveal significant different values for V(1) and V(2): +3.3 and +4.7

respectively. A description of all these features is summarized in Figure 1, with a comparison between Na 3 V III 2 (PO 4 ) 2 F 3 and Na 1 V IV 2 (PO 4 ) 2 F 3 that emphasizes especially the vanadium environments of these two compositions.

In order to directly investigate this proposed disproportionation reaction and thus getting a better insight on the redox processes involved during the sodium extraction from Na 3 V 2 (PO 4 ) 2 F 3 an operando vanadium K-edge X-ray absorption spectroscopy (XAS) study has been performed during the cycling of a Na//Na 3 V 2 (PO 4 ) 2 F 3 battery. The collected XAS data are directly sensitive to the vanadium oxidation state and to the local structure around the vanadium centers. The whole dataset has been analyzed using chemometrics through principal component analysis (PCA) and multivariate curve resolution (MCR). Due to the peculiar environment of the vanadium in this material and its complex behavior during oxidation, the use of linear combinations of well-chosen XAS spectra, either standards or identified compositions upon sodium extraction, turned out to be not possible. The combined use of PCA and MCR offers a blind-source separation method to process a large data set of this operando XAS experiment by expressing the data based on their similarities and differences. These methods mainly used for following in situ catalysis reaction. [START_REF] Voronov | Multivariate curve resolution applied to in situ X-ray absorption spectroscopy data: An efficient tool for data processing and analysis[END_REF][START_REF] Rochet | Influence of the preparation conditions of oxidic NiMo/Al 2 O 3 catalysts on the sulfidation ability: a quick-XAS and raman spectroscopic study[END_REF][START_REF] Cassinelli | Timeresolved XAS/MS/Raman monitoring of mutual copper self-reduction and ethanol dehydrogenation reactions[END_REF][START_REF] Cassinelli | Multivariate curve resolution analysis applied to time-resolved synchrotron X-ray Absorption Spectroscopy monitoring of the activation of copper alumina catalyst[END_REF] Note that this method has also been used to analyze lithium insertion in graphite and lithium extraction in LiFePO 4 via neutron powder diffraction [START_REF] Rodriguez | In situ analysis of LiFePO 4 batteries: signal extraction by multivariate analysis[END_REF] which constitutes a relevant approach to investigate the operando behavior of batteries upon cycling. 51 V MAS NMR (magic-angle spinning nuclear magnetic resonance) was performed on materials obtained at the end of the charge, close to the composition Na 1 V IV 2 (PO 4 ) 2 F 3 , in order to confirm the formation of the diamagnetic ion V V (t 2g 0 ) which is in that case the probed nucleus. By the reported results of the combined use of XAS and NMR we were then able to follow the evolution of the local electronic and structural properties as a function of Na content in operando conditions and to reveal the coexistence of V III and V V at the end of the charge.

Experimental part

The Na 3 V 2 (PO 4 ) 2 F 3 powder sample was synthetized according to the method already described in details elsewhere. [START_REF] Bianchini | Na 3 V 2 (PO 4 ) 2 F 3 revisited: a high-resolution diffraction study[END_REF] From the direct comparison (given in supplementary information, Figure S1) with a previous XAS study performed on Na 3 V 2 (PO 4 ) 2 F 3-y O y (0 ≤ y ≤ 0.5) [START_REF] Broux | Strong impact of the oxygen content in Na 3 V 2 (PO 4 ) 2 F 3y O y (0≤y≤0.5) on its structural and electrochemical properties[END_REF] the amount of V IV defect (y) was estimated to be 0.03 in the as-prepared material. Considering the small content of defects it was neglected in the discussion of the overall mechanism.

Electrodes were prepared by mixing Na 3 V 2 (PO 4 ) 2 F 3 with Carbon Super P (Csp) as conductive additive and polytetrafluoroethylene (PTFE) as the binder (60:30:10 wt%). The mixture was ground in a mortar until a homogeneous film was obtained. Small pellets of 8 mm diameter and ∼300 μm thickness were used in an in situ cell already described in details elsewhere. [START_REF] Leriche | An electrochemical cell for operando study of lithium batteries using synchrotron radiation[END_REF] The electrochemically active part of the cell is obtained from the stacking of the positive electrode, glass fiber separators (Whatman) wetted with a 1 M solution of NaPF 6 (Strem Chemical; 99%) in a mixture of ethylene carbonate and dimethyl carbonate (EC:DMC = 1:1) and sodium metal as the counter electrode. This stacking is inserted between two beryllium windows acting as X-ray-transparent current collectors. The internal pressure within the cell is applied through a spring to ensure an optimal electrical contact. The comparison of the electrochemistry obtained in this in-situ cell is compared in Figure S2 (given in supplementary information) to that obtained in classical coin cells, excellent reproducibility was obtained.

Vanadium K-edge XAS measurements were collected in the fluorescence mode at the CLAESS beamline [START_REF] Simonelli | CLAESS: The hard x-ray absorption beamline of the ALBA CELLS synchrotron[END_REF] of the ALBA Synchrotron (Barcelona, Spain) using a single channel silicon drift detector.

The synchrotron radiation of a wiggler source was monochromatized by means of a Si(111) double crystal monochromator. Rh-coated collimating and toroidal mirrors were used to optimize the energy resolution and to reject higher harmonics and to control the beam size at the sample position (0.7 x 1.5 mm 2 ). The calibration of the monochromator was performed by measuring the spectrum of a vanadium metal foil in transmission mode. The absolute energy reproducibility of the measured spectra was ± 0.1 eV. The incoming energy resolution around the vanadium K-edge can be estimated below 0. First of all, the operando XANES spectra were normalized to the atomic absorption estimated by a linear fit to the data in the EXAFS (Extended X-ray Absorption Fine Structure) region after a pre-edge linear background subtraction. To increase the signal-to-noise ratio the scans were merged in bunches of 5 which give a total of 40 spectra representing the complete charge. Then the data were globally analysed using a statistical tool named PCA using the computer program Matlab. 36 PCA is a chemometric factor analysis tool generally used to discover the minimal particular structures in multivariate spectral data sets. Here it allows to determine the number of independent component contributing to the whole series of collected spectra during the electrochemical cycling. The number of principal components determined was then used as the basis for Multivariate Curve Resolution-Alternating Least Squares (MCR-ALS) analysis. [START_REF] Jaumot | MCR-ALS GUI 2.0: New features and applications[END_REF][START_REF] Jaumot | A graphical user-friendly interface for MCR-ALS: a new tool for multivariate curve resolution in MATLAB[END_REF] This algorithm allows the stepwise reconstruction of the spectral components which is necessary for interpreting the whole set of operando spectra.

The same electrochemical setup, using classical cells and without binder, was used to prepare the deintercalated compositions near Na 1 V IV 2 (PO 4 ) 2 F 3 . The cells were charged in galvanostatic mode at a C/20 rate up to 4.3 V vs Na + /Na, and maintained at this voltage for several hours (C rate is corresponding to the exchange of one electron per hour per mole of compound). The materials as obtained were subsequently washed in DMC to remove any traces of the electrolyte salt and then dried in vacuum. The so-obtained powders were stored under inert atmosphere and subsequently characterized by capillary synchrotron XRD and 51 V MAS NMR. The ex situ characterization of the phase pure Na 1 V 2 (PO 4 ) 2 F 3 is difficult due to its lack of stability because of its high equilibrium potential. It is revealed by the presence, as a secondary phase, of the more reduced composition Na 1.3 V 2 (PO 4 ) 2 F 3 already identified in a previous work. [START_REF] Bianchini | Comprehensive investigation of the Na 3 V 2 (PO 4 ) 2 F 3 -NaV 2[END_REF] Its presence as a secondary phase is not an issue since it does not prevent the detection of the possible V V signature in the major phase Na 1 VPF.

High resolution Synchrotron powder X-ray diffraction (SXRD) data were collected using the MSPD diffractometer at ALBA (Barcelona, Spain). The powders were placed in a 0.5 mm diameter capillary and data recorded in Debye-Scherrer geometry with a wavelength of 0.9540 Å in the 2θ angular range of 1 -70 ° with a 0.005 ° step and an accumulation time of 6 minutes. 51 V MAS NMR spectra were recorded on a Bruker Avance III 500 MHz spectrometer (131.5 MHz resonance frequency for 51 V) using a standard Bruker 2.5 mm MAS probe with a 30 kHz typical spinning frequency. The recycle time D 0 = 0.2 s is long enough to avoid T 1 saturation effects. The baseline distortions resulting from the spectrometer dead time (5-10 µs) were removed computationally using a polynomial baseline correction routine. α-VOPO 4 with a single V site that exhibits a signal located at -765 ppm [START_REF] Siegel | 51 V magic angle spinning NMR in VOPO 4 phases[END_REF] was used as secondary reference.

Results and discussion

Figure 2 displays the electrochemical data obtained during the operando XAS experiment at the rate of C/10 up to 4.4 V vs. Na + /Na as a function of the number of normalized collected spectra. The soobtained electrochemical curve is identical to those already reported in previous studies. [START_REF] Bianchini | Comprehensive investigation of the Na 3 V 2 (PO 4 ) 2 F 3 -NaV 2[END_REF][START_REF] Bianchini | Na 3 V 2 (PO 4 ) 2 F 3 revisited: a high-resolution diffraction study[END_REF] Figure 2 displays also a series of selected XAS spectra obtained during the charge between Na 3 V III 2 (PO 4 ) 2 F 3 (Na 3 VPF -purple line) and Na 1 V IV 2 (PO 4 ) 2 F 3 (Na 1 VPF -brown line). The inset focuses on the pre-edge region.

Vanadium K-edge XANES exhibits several features as already largely testified by experimental and theoretical studies. [START_REF] Wong | K-edge absorption spectra of selected vanadium compounds[END_REF][START_REF] Poumellec | Experimental and theoretical studies of dipole and quadrupole contributions to the vanadium K-edge XANES for VOPO 4 ⋅2H 2 O xerogel[END_REF][START_REF] Wu | Quadrupolar transitions and medium-range-order effects in metal K-edge x-ray absorption spectra of 3d transitionmetal compounds[END_REF] The main absorption edge, at around 5485 eV, is related to the valence state of the transition metal and corresponds to the strong contribution of the 1s → 4p dipole allowed transition. The weak pre-edge contributions around 5470 eV, arise from 1s → 3d forbidden quadrupolar transitions that become allowed by orbital 3d-p hybridizations. [START_REF] Wong | K-edge absorption spectra of selected vanadium compounds[END_REF] Both energy positions and intensity distributions of the pre-peak signals are linked to the electronic structure and the local structural environment of the probed element.

Firstly, the overall 2.5 eV shift of the main absorption edge towards higher energy from Na 3 V III 2 (PO 4 ) 2 F 3 (around 5484.5 eV) to Na 1 V IV 2 (PO 4 ) 2 F 3 (around 5487 eV) is consistent with an increase of the vanadium average oxidation state. Pre-edge lines are complex to identify and attribute precisely since there are several contributions all along the electrochemical oxidation of NVPF. It can be stated however that their evolution is consistent with the increase of vanadium oxidation state since the contributions at higher energy are growing and the overall integrated intensity is increasing. In fact, an increase of the contribution at higher energy can be ascribed to a progressive modification/distortion of the local symmetry [START_REF] Wong | K-edge absorption spectra of selected vanadium compounds[END_REF] that leads to an increasing 3d-4p orbital mixing, that means a larger overlap of the 3d vanadium orbitals with the 2p oxygen orbitals. This mixing enhances the theoretically forbidden 1s → 3d transition, inducing an increase in the pre-edge intensity. [START_REF] Salem | Energy shift and structure of the K-absorption edge of vanadium in some vanadium compounds[END_REF] During sodium deintercalation, the vanadium oxidation state is expected to increase linearly, and the absorption edge is supposed to evolve in the same way. However the energy position of the main absorption edge shows a nonlinear evolution, in a way strongly correlated to the electrochemical data.

In fact, a rapid evolution is observed between spectra 17 and 20, together with the sudden increase of potential in the voltage curve (Figure 3). This particular evolution suggests that the whole set of data cannot be described as a simple linear combination of the initial and final spectra, as it could have been expected for a process characterized by the gradual vanadium oxidation from the average V III to V IV oxidation state. In fact more complex processes seem to occur during the redox reaction suggesting unexpected local electronic rearrangements. From this finding and the impossibility to do relevant comparisons with known references due to the specific environment of vanadium in the Na 3 V 2 (PO 4 ) 2 F 3 framework (i.e. in V 2 O 8 F 3 bi-octahedra displayed in Figure 1) data were analyzed using principal component analysis (PCA) combined with multivariate curve resolution (MCR). In order to avoid an over-interpretation of the experimental data and to extract the maximum amount of useful information from the operando XANES spectra PCA was performed on the whole dataset. In this method each spectrum is considered as an n-dimensions vector (n corresponding to the number of measured absorption points). All measured spectra form a matrix, which can be reduced to its lowest dimensionality by the use of orthogonal factor space: this factor analysis is the principal component analysis (i.e. PCA).

The first principal component has the largest initial variance, i.e., alone it takes into account as much as possible the whole set of data. The following principal components are built orthogonal to the first one, and allow the reproduction of the experimental spectra with the highest possible variance via adapted linear combinations of them. An orthogonal basis of uncorrelated components whose the dimension is less than or equal to the number of original spectra is thus obtained. In the different linear combinations which are built to reconstruct each experimental spectrum, the multiplication factors of the principal components are commonly called scores.

The variance plot obtained from PCA analysis given in Figure 4 shows rather clearly that about 99.9 % of the variance of the XANES spectra can be described using only three principal components, the residual part being supposedly due to experimental noise. A visualisation of the evolution of the respective scores comforts this hypothesis, with both shape and contribution of the first three components varying gradually through the different processes (as expected in a gradual evolution of the studied material during oxidation) and those of the following ones changing in a totally erratic way, thus probably representing mainly experimental noise.

Even though the so-obtained principal components are only orthogonal mathematical functions and not real XANES spectra, their determination reveals the number of independent spectral components that can be used to fit the whole set of experimental data. Several methodologies are then available to reconstruct the corresponding real XANES components as well as their evolution such as the MCR-ALS (Multivariate Curve Resolution -Alternating Least Squares) analysis tool. [START_REF] Jaumot | MCR-ALS GUI 2.0: New features and applications[END_REF][START_REF] Jaumot | A graphical user-friendly interface for MCR-ALS: a new tool for multivariate curve resolution in MATLAB[END_REF] A detailed description of this method from a theoretical point of view is given by Tauler et al. [START_REF] Tauler | Multivariate curve resolution applied to second order data[END_REF][START_REF] De Juan | Chemometrics applied to unravel multicomponent processes and mixtures: revisiting latest trends in multivariate resolution[END_REF] who was the first to propose this method for the analysis of in situ spectroscopic data. On the other hand, the intrinsic limits of the application of this method are discussed in ref. [START_REF] Ruckebusch | Multivariate curve resolution: a review of advanced and tailored applications and challenges[END_REF] In this work, the MCR-ALS analysis was performed by applying the following constraints for both spectral components and their respective concentrations, in order to reproduce a physically and chemically meaningful behavior: (i) non-negativity of the concentrations of the components, (ii) unimodality for one of the three components (only a single highest value of concentration) and (iii) closure (sum of the 3 components always equal to 100% of the intensity). MCR-ALS method converged after 25 iterations with a variance of 99.998 %.

As a result, the three XANES spectra shown in Figure 5 were obtained together with the evolution of their concentration along the electrochemical reaction path. According to the evolution of their relative concentrations, the three XANES spectra can be described as follows: the first component (C1) corresponds to the pristine electrode material Na 3 V III 2 (PO 4 ) 2 F 3 and the third component (C3) to the desodiated material Na 1 V IV 2 (PO 4 ) 2 F 3 obtained at the charged state. However the second component (C2) seems to be an intermediate component since it cannot be identified as a given composition.

Concerning the evolution of their relative intensities, on the first voltage plateau it appears that C1 is continuously decreasing while C2 is increasing. On the other hand once the voltage jump is reached C3 starts growing while both C1 and C2 decrease gradually.

The XANES spectra extracted from the MCR-ALS analysis are in good agreement with the observations made considering the quantitative evolutions of the components: C1 and C3 are, as expected, almost identical to the first (number 1) and the last (number 40) operando measured XANES spectra, respectively. Since C2 never reaches 100 % of the relative quantity it cannot be identified with a measured spectrum. Comparing the absorption energy of the main edge jump obtained for each component, 5484.5 eV for C1, 5485.8 eV for C2 and 5486.9 eV for C3, it appears that it evolves linearly as a function of the oxidation state of vanadium, as expected, if the component C2 represents a composition with an average oxidation state between V III and V IV , very close to V +3. [START_REF] Shakoor | A combined first principles and experimental study on Na 3 V 2 (PO 4 ) 2 F 3 for rechargeable Na batteries[END_REF] .

Besides, during the first voltage plateau only V III and V IV are expected to be present, so C2 can reasonably be ascribed to a mixture of about 50 % of V III and 50 % of V IV . Assuming these hypotheses, i.e. an average oxidation state of +3 for C1, +3.5 for C2 and +4 for C3, the overall average oxidation state which can be calculated from the relative quantities evolution is compared in Figure 6 to the actual oxidation state expected from the electrochemical experiment. It shows a rather good agreement despite a small deviation from the theoretical oxidation state given by the straight blue line.

It is worth noticing that the evolution of the three components is much more simple compared to the complex phase diagram previously reported by Bianchini et al. [START_REF] Bianchini | Comprehensive investigation of the Na 3 V 2 (PO 4 ) 2 F 3 -NaV 2[END_REF] that revealed complex structures showing for some of them commensurate modulations and thus complex charge and/or ionic rearrangements upon sodium extraction. Considering the peculiar framework whose basic units are the V 2 O 8 F 3 bi-octahedra linked to each other by phosphates (see Figure 1), and that each vanadium has another vanadium in its local environment (V-V interatomic distance within V 2 O 8 F 3 bi-octahedron:

3.96 Å), and only one (as the next V-V distance is: 4.73 Å), it is reasonable to assume that the components obtained here which are associated to local environments seen by XAS reflect the local pair arrangement. The three components can thus be expressed by the local pair arrangements V III -V III for C1, V III -V IV for C2 and V IV -V IV for C3 (or V III -V V in case of charge disproportionation) within vanadium bi-octahedron.

The pre-edge show four main contributions observed at 5466, 5468, 5469.5 and 5471 eV with different relative intensities along sodium deintercalation. As no fully relevant standards were found, it was difficult to give a full interpretation for each of the contribution. Nevertheless the increase from C1 to C3 of the whole pre-edge integrated intensity and the shift towards higher energies of the average energy position are in agreement with the decrease of the number of d electrons in the vanadium electronic configuration and to the increasing distortion of the octahedral coordination of the vanadium centres, thus with an increasing probability of the 1s to 3d transitions.

In Figure 7 we report the measured absorption pre-peak, the corresponding MCR calculated one and their corresponding difference. Considering the residues of the difference between the measured and the calculated data, they are subtle but still observable within the pre-edge. They could be attributed to slight modifications of the components along the reaction that are too small to be taken into account by the PCA as well by the MCR analysis. This could be responsible for the small deviation from the actual oxidation state mentioned before and shown in Figure 6, this limited "discrepancy" highlighting the accuracy of this mathematical analysis.

Nevertheless at this point of the study there is no clear proof of the presence of V V in the material.

Comparing the spectrum obtained previously by Serras et al. [START_REF] Serras | Electrochemical Na extraction/insertion of Na 3 V 2 O 2x (PO 4 ) 2 F 3-2x[END_REF] at the end of the charge for the compound Na 1 V 2 (PO 4 ) 2 F 1.4 O 1.6 (showing an average oxidation state of +4.5 and thus a non-negligible amount of V V ) gives only indirect but already interesting information. In this case the presence of the V V is mainly associated to the growth of a peak around 5471 eV in the pre-edge region. In our case, as depicted in Figures 2 and5 a similar phenomenon is observed as a peak at 5471 eV rises concomitant with the appearance of C3 which might indicate that the charge disproportionation is occurring at the first voltage jump from the extraction of the second Na + ion. It is interesting to notice, that there is a delay in the apparition of the peak at 5471 eV in the spectra calculated using the MCR components compared to the experimental ones. This observation could explain the slightly lower average oxidation state of vanadium found by MCR compared to the one expected from electrochemistry.

In order to ensure the presence of V V ions in the fully oxidized compound 51 V MAS NMR was used. Indeed, NMR signals have been observed by high resolution MAS NMR at room temperature only for compounds contained diamagnetic V V ions as observed for α-VOPO 4 that was used as secondary reference. [START_REF] Jordan | Crystal structure of α-VPO 5[END_REF] Firstly, in order to assess the quality of the compound measured by 51 V MAS NMR, a SXRD pattern was preliminary collected in a sealed capillary to prevent contact with ambient air. SXRD data of the near fully charged material were analyzed by Rietveld refinement using Cmc2 1 and I4/mmm space groups to describe Na 1 VPF and Na 1.3 VPF respectively as displayed in Figure 8. This material is composed of 63 % of Na 1 V IV 2 (PO 4 ) 2 F 3 and 27 % of Na 1.3 V 2 (PO 4 ) 2 F 3 and its corresponding 51 V NMR MAS spectra are compared to those of Na 3 V 2 (PO 4 ) 2 F 3, and α-VOPO 4 . In agreement with literature α-VOPO 4 exhibits a single narrow signal (-765 ppm) and a set of spinning side bands. [START_REF] Siegel | 51 V magic angle spinning NMR in VOPO 4 phases[END_REF] The absence of significant quadrupolar second-order lineshape indicates that the quadrupolar interaction is moderate.

No signals were observed for 51 V in Na 3 V III 2 (PO 4 ) 2 F 3 , as expected. The NMR spectrum shows however, some residual spinning sidebands of [START_REF] Rao | Structural and electrochemical behaviour of sputtered vanadium oxide films: oxygen non-stoichiometry and lithium ion sequestration[END_REF] Na MAS NMR signal, as the sodium ( 23 Na) and vanadium ( 51 V) resonance frequencies are very close 132.3 and 131.5 MHz respectively. For the charged compound composed of 63 % of Na 1 V 2 (PO 4 ) 2 F 3 and 27 % of Na 1.3 V 2 (PO 4 ) 2 F 3 , the spectrum exhibits clearly signals located in the [-650, -850] ppm region and their corresponding spinning sidebands. These signals are assigned to V 5+ ions present in the charged compound. The signal located around -680 ppm may be assigned to V 5+ in Na 1 V 2 (PO 4 ) 2 F 3 which is the major phase present and the others to V 5+ ions in Na 1.3 V 2 (PO 4 ) 2 F 3 that may exhibit different environments versus Na + .

The presence of V V in the near charged state of the battery indicates that V V has to be considered in the redox processes. As previously proposed by Bianchini et al., [START_REF] Bianchini | Comprehensive investigation of the Na 3 V 2 (PO 4 ) 2 F 3 -NaV 2[END_REF] charge disproportionation really occurs and this information is contained in the third component that can be ascribed to Na 1 V III V V (PO 4 ) 2 F 3 . The C3 quantity starts growing concomitantly with the voltage jump indicating that a partial charge disproportionation occurs from the average composition Na 2 V III V IV (PO 4 ) 2 F 3 .

In summary a possible explanation can be given if V 2 O 8 F 3 bi-octahedra within this peculiar framework is considered to be one entity. During the first voltage plateau it is easier to oxidize V III -V III rather than a V III -V IV entity whence the formation of V III -V IV pairs up to reach Na 2 V III V IV (PO 4 ) 2 F 3 (i.e. average oxidation state +3.5). Once the Na 2 V III V IV (PO 4 ) 2 F 3 configuration is reached further sodium extraction (i.e. oxidation) leads to the formation of V III -V V pairs. The appearance of the charge disproportionation at this point of the charge is likely due to the fact that the structure does not allow the formation of V IV -V IV pair in V 2 (PO 4 ) 2 F 3 bi-octahedra.

Conclusion

Based on structural observations it has been proposed by Bianchini et al. that a charge disproportionation of 2V IV into V III and V V occurs in Na 1 V IV 2 (PO 4 ) 2 F 3 . [START_REF] Bianchini | Comprehensive investigation of the Na 3 V 2 (PO 4 ) 2 F 3 -NaV 2[END_REF] In this study 51 V MAS NMR has confirmed the presence of the diamagnetic ions V V in the near fully deintercalated materials supporting that the charge disproportionation actually occurs. Furthermore, operando vanadium Kedge XANES study has shown that this phenomena does not occur only at the end of the charge, but immediately after the extraction of one Na + ion i.e. after reaching the composition Na 2 V III V IV (PO 4 ) 2 F 3 .

One way to understand this feature is to consider the bi-octahedra, and not the octahedra, as the basic structural unit. This peculiar geometry combined with mixed valence seems to allow electronic exchanges that induce original configurations. The characterization of these electronic exchanges is crucial in order to establish the groundwork required to fully understand the sodium extraction/insertion mechanism combined with electronic properties in this compound and might be extended to other related compounds having this bi-octahedron geometry with the composition V 2 (PO 4 ) 2 F 3-x O x

Figure captions

Figure 1: Cristal structures of Na 3 V 2 (PO 4 ) 2 F 3 and Na 1 V 2 (PO 4 ) 2 F 3 with PO 4 in green, VO 4 F 2 octahedra in blue and Na atoms in red. In Na 1 V 2 (PO 4 ) 2 F 3 significant different vanadium environments are highlighted with light and dark blue corresponding to regular and distorted environments respectively.

  4 eV. The incoming and, in the case of calibration, outcoming flux have been measured by customized ionization chambers. During the electrochemical cycle, 200 individual spectra were recorded in continuous mode (around 4 min per spectrum) in the XANES (X-ray Absorption Near-Edge Structure) region of vanadium K-edge (5465 eV) with equidistant energy steps of ΔE = 0.2 eV.

Figure 2 :

 2 Figure 2: Top: voltage as a function of the spectrum number and corresponding to the galvanostatic charge obtained at C/10 during operando XAS data collection. Bottom: operando vanadium K-edge XANES spectra collected during the charge. The inset is focused on the pre-edge region.

Figure 3 :

 3 Figure 3: Left: voltage (black line) as a function of the spectrum number and corresponding to the galvanostatic charge obtained at C/10 during operando XAS data collection with the associated 1 st derivative (blue line). Right: 2D projection of operando vanadium K-edge XANES spectra collected during the charge and emphasizing discontinuities at voltage jumps.

Figure 4 :

 4 Figure 4: Variance plot of PCA analysis. Evolution of the scores of the 6 first components.

Figure 5 :

 5 Figure 5: Top: three components extracted from MCR analysis. The inset is focused on the pre-edge region. Bottom: evolution of the relative quantities of the three principal components compared with the experimental electrochemical data (black solid line).

Figure 6 :

 6 Figure 6: Experimental electrochemical curve (black solid line) associated to the expected oxidation state (blue solid line). The blue dotted line represents the oxidation state calculated from the relative quantities of the three components associated with their observed oxidation state: C1 = +3; C2 = +3.5 and C3 = +4.

Figure 7 :

 7 Figure 7: 2D projection of operando vanadium K-edge XANES spectra of the pre-edge -Top: measured data -Middle: reconstructed data from MCR analysis -Bottom: difference between the measured data and the reconstructed data.

Figure 8 :

 8 Figure 8: Top: Rietveld refinement of the electrochemically oxidized compound using Cmc2 1 and I4/mmm space groups for Na 1 VPF and Na 1.35 VPF respectively. Bottom: 51 V MAS NMR spectra NMR of the electrochemically oxidized compound (red line), of the pristine material Na 3 VPF (purple line) and of V V PO 4 O (green line).

Figure 1 Na

 1 Figure 1

Figure 3 Please this figure has to be a two columns figure Na 3

 33 Figure 3

Acknowledgments

The authors thank Matteo Bianchini (ILL/LRCS/ICMCB) for his technical assistance and fruitful discussions, also C. Marini and W. Olszewski at CLAESS beamline at ALBA (Barcelona, Spain) for their help while preparing the XAS experiment and MSPD beamline at ALBA for Synchrotron X-ray diffraction, and G. Gachot (LRCS) for fruitful discussions on the electrolyte. The authors also acknowledge RS2E for the funding of TB's postdoctoral fellowship. This project has received funding from Région Nouvelle Aquitaine, the French National Research Agency (STORE-EX Labex Project ANR-10-LABX-76-01 and SODIUM Descartes project ANR-13-RESC-0001-02), and the European Union's Horizon 2020 research and innovation program under grant agreement No 646433-NAIADES.

Supporting Information

Comparison of XAS spectra for Na 3 V 2 (PO 4 ) 2 F 3-y O y (0  y  0.5) and electrochemical results obtained in in-situ cells and in classical coin cells


Kerian Thuillier

AlPha: A Mixed Integer Linear Programming Approach for Genome Haplotyping

Keywords: metagenomic, genome assembly, milp, multicommidity-flow

Background: Nowadays, biological samples often contain dna of several species or individuals. When sequenced, we loose the origin of each fragments of dna and their position on the molecule. The problem of associating a dna fragment to a version of the genome and to find its position is called the haplotyping problem. Results: In this report, we present two new Mixed Integer Linear Programming (milp) methods based on multicommodity-flow to solve this problem. Unlike previous approaches, the resolution avoid the use of heuristics by using global optimisation.

Conclusion:

We do not currently known the quality of our results. However, proposed a comparison between the solving time complexity of our two milp models. We also proposed an experimental protocol to test the solution quality.

Introduction

A sequencer is a machine able to read short fragments of dna called reads. A read is a short sequence of characters representing nucleotides. As only short fragments of the genome are read by the sequencer, we loose the read position on the genome. Reads can also be paired, in this case we know that both reads come from the same molecule. Standard genome assembly consists in ordering all the reads to reconstruct the whole genome. However, such assembly reconstructs the genome by consensus. When an assembler has the choice between several reads to complete the assembled sequence, it adds the most-read read. Thus, standard genome assembly could not reconstruct more than one genome from a set of reads.

A sample contains often more than one individual or comes from a polyploid organism, an organism which has several versions of its genome. For example, human beings are diploids, they have two versions of the human genome: one version inherited from the father and the other inherited from the mother. The different versions of one genome are called haplotypes. During the sequencing process, all the haplotypes are sequenced together. Reconstructing all the haplotypes from a set of reads is called haplotype-aware genome assembly. Recently proposed algorithmic solutions solve this problem [START_REF] Baaijens | Overlap graph-based generation of haplotigs for diploids and polyploids[END_REF][START_REF] Baaijens | Strain-aware assembly of genomes from mixed samples using variation graphs[END_REF]. However, these solutions are error-prone and dedicated to small bacterial genomes.

In this report, we propose two new approaches for haplotype-aware genome assembly based on Mixed Integer Linear Programming (milp). The first one relies on fact graph (see section 3.1), while the second one relies on alleles graph (see section 4.1). These graphes come from the output of the DiscoSNP [START_REF] Peterlongo | Discosnp++: de novo detection of small variants from raw unassembled read set(s)[END_REF] algorithm. DiscoSNP is a de novo1 variant detection tool. A variant is a single nucleotide difference between several genomes or haplotypes.

First, we will introduce both the genome assembly problem and the haplotyping problem (Section 1). Then, we will present our first approach relying on fact graph (Section 3) and its limits. Afterward, we will introduce our second approach based on allele graph (Section 4). The advantage and the limits of this approach are discussed in this section. Finally, we will present our result validation protocol (Section 5) and develop our work through a discussion (Section 6).

1 Haplotype-aware genome assembly

Genome assembly Problem

Modern sequencers allow reading a genome by slicing it in fragments. Those fragments are called reads. A read is a sequence of about 200 characters: A, C, G or T. Note that sequencer are error-prone, a read could contain errors. A sequencer output is a list of billions of reads. The output often does not contain any information about the read positions. Some modern sequencers can provide more information as the distance on the genome between some reads. An example of sequencer output is provided Figure 1.

Genome assembly problem is a scheduling problem. Each read is ordered to reconstruct the genome. Read position is determined by its neighbours. A read is positioned between two reads with which it overlaps. An overlap is a read's subsequence which is prefix for a read and suffix for another read. Thus, the Genome assembly problem could formally be defined as below.

GENOME ASSEMBLY PROBLEM INPUT:

R a set of reads OUTPUT:

S the longuest sequence of reads such that: 1. each element of S is in R; 2. each read of R is at most once in S; 3. two successives reads of S are overlapping. The number of unused element of R is minimised.

The genome assembly problem is considered as an np-complet problem due to the method used to solve it: hamiltonian path or longuest path [START_REF] Medvedev | Computability of models for sequence assembly[END_REF].

Genome assemblers such as SPAdes [START_REF] Bankevich | Spades: A new genome assembly algorithm and its applications to single-cell sequencing[END_REF] are commonly used to reconstruct genomes. However, these algorithms are not error-proof and have some restrictions. First, genomes contain a lot of repeats (nucleotide sequences repeated several times in the genome). The handling of repeats is an unresolved problem and an active research topic. Secondly, standard assemblers can only reconstruct one genome from a set of reads provided by one individual.

Haplotyping Problem

Humans are diploid, each chromosome is present twice. Those chromosomes are almost identical but have slight variations (about 1% of divergence). For humans, one version of these chromosomes was inherited from the mother while the second was inherited from the father. Chromosomes inherited from the mother form a haplotype as those inherited from the father. In the same way, a cup of sea water contains billions of micro-organisms. Each organism is a haplotype. When these elements are sequenced, reads from all haplotypes are merged. Thus, standard genome assembly strategies do not work. For example, if the human genome is sequenced and assembled with standard strategies, only one haplotype will be produced. This haplotype will be a matchup of both real haplotypes. Even more absurd, if the water cup is sequenced and assembled, the reconstructed genome will be a chimera.

That's where the haplotyping problem takes place. Let R be a set of reads obtained by sequencing k haplotypes {h 1 ... h k }. The number of haplotypes k is unknown. The haplotyping problem consists in finding k, reallocating each element of R to its associated haplotype and assembling them. In the end, each read r ∈ R has a position in a haplotype or is ignored. Two versions of this problem exists: first one is the version described above, while second one has k fixed or minimised.

HAPLOTYPING PROBLEM INPUT:

R a set of reads (k if k is fixed) OUTPUT:

k sequences of reads {S_1, ..., S_k} such that: 1. for each sequence S_i, each element of S_i is in R; 2. a read is at most in one sequence S_i; 3. for each S_i, two successives reads of S_i are overlapping. The number of unused element of R is minimised.

The haplotyping problem is a well-known np-complete problem. The haplotyping problem is the same as the genome assembly problem if there is only one haplotype. As this last problem is np-complet, we deduce that haplotyping problem is too.

In recent papers, haplotype-aware assemblers have been proposed [START_REF] Baaijens | Overlap graph-based generation of haplotigs for diploids and polyploids[END_REF][START_REF] Baaijens | Strain-aware assembly of genomes from mixed samples using variation graphs[END_REF]. However, these approaches are error-prone and can only process small bacterial genomes.

DiscoSNP: Variant Phasing

As different haplotypes provided by a single species are almost identical, it is more convenient to simplify the problem. Instead of working with reads, we could work with variants -Single Nucleotide Polymorphism (snp). A snp is a variation of only one nucleotide between two haplotypes. Working with snp instead of reads is only removing information that is known to be shared between haplotypes. It allows to remove long non discriminative genomic regions.

DiscoSNP

DiscoSNP [START_REF] Peterlongo | Discosnp++: de novo detection of small variants from raw unassembled read set(s)[END_REF] is a reference-free variant phasing algorithm. From a set of reads, Dis-coSNP computes the variants and phases them. The different versions of a variant are called alleles. For instance, the pair sequences accaagttt and accatgttt determine the variant a versus t (central position), while a is one allele of this variant and t the other allele. Two alleles are phased if and only if there is a read overlapping both alleles. Phased alleles come from the same haplotype. A group of phased alleles is called a fact. Facts can be paired. Two facts are paired if they are overlapped by two paired reads (i.e. information provided by sequencers meaning that both reads come from the same haplotype). Figure 2 illustrates the phasing principle.

Outputs

DiscoSNP returns three files that contains the fact information:

• one text file containing all the raw facts and raw paired facts with their coverage, i.e. the number of reads mapping them. Each line of this file is a fact or a paired facts. A fact is described as below:

1h_0;2l_-30; => 5 Here, 1h is the allele h of the variant 1 while 2l is the allele l of the variant 2. The value after the underscore represent the distance of the allele to the previous one: the allele 2l start 30 nucleotide before the end of the allele 1h. This fact has been seen by 5 reads. A paired fact is described as below: 1h_0;2l_-30; 5h_0;6h_10;7l_-30; => 2 This line represents two paired facts: 1h_0; 2l_ -30; and 5h_0; 6h_10; 7l_ -30;. Due to sequencing constraints, the second fact must be reversed to get the real paired facts. The real paired facts is thus 1h_0; 2l_ -30; and -7l_0; -6h_ -30; -5h_10 (the reverse of 5h_0; 6h_10; 7l_ -30;) which has been seen by 2 reads.

• two fasta files2 containing alleles informations. These files contain the nucleotide sequences associated with each allele.

Post-process

Before being processed, DiscoSNP output is simplified. There are two phased of simplification. First one consists in filtering the uncoherent alleles. Second one consists in spliting the facts into working zone. A working zone is a cycle-free connected components where erroneous data were removed.

3 AlPha-Fact: milp approach for genome haplotyping

Our first model is based on the interaction between the facts. These interactions are modelled as a fact graph. This model is implemented according to the Mixed Integer Linear Programming (milp) paradigm. This paradigm relies upon declarative programming; problems are written in the form of linear objective function under linear constraints [START_REF] Bradley | Applied Mathematical Programming, chapter 9-Integer Programming[END_REF]. Used variables could be either integers or real values.

Fact Graph

From DiscoSNP output, we can build a fact graph. A fact graph is a graph modelling interactions between facts. There are four interactions between facts:

Overlaps There is an overlaps between two facts f 1 and f 2 iff the last phased alleles of f 1 are the same as the first phased alleles of f 2 . For example, let f 1 = 1h; 2l; 3h;, f 2 = 2l; 3h; 4l; and f 3 = 2l; 3l; 4h then there is an ovelaps of size 2 between f 1 and f 2 , and no overlaps between f 1 and f 3 .

Successive Two facts f 1 and f 2 are successives iff the last variant phased by f 1 are known as preceding the first variant phased by f 2 . For example, let f 1 = 1h; 2l; 3h;, f 2 = 4h; 5l; and f 3 = 3l; 4l; then there is no overlaps between those three facts. However, we know from f 3 that variant 3 comes before variant 4. Thus, we conclude that f 1 and f 2 are successives.

Incompatible Two facts f 1 and f 2 are incompatibles iff they phased different alleles of the same variant. For example, let f 1 = 1h; 2l; and f 2 = 1l; 2l; then f 1 and f 2 are incompatibles: f 1 phases 1h while f 2 phases 1l which are two alleles of variant 1.

Paired Two facts f 1 and f 2 are paired iff they have been seen as paired in the DiscoSNP output.

Let formally define a fact graph G = (V, E, Ext, A, weight) where:

V: set of facts representing the graph vertices.

E: set of edges, there are four different types of edges corresponding to the four interactions between facts. First one is the overlaps edge (f 1 , f 2 , overlaps) representing the overlap of two facts f 1 and f 2 . Second one is the successive edge (f 1 , f 2 , successive) modelling the successive relation between the facts f 1 and f 2 . Third one is the incompatible edge (f 1 , f 2 , incompatible). Last one is the linked edge (f 1 , f 2 , links) meaning that f 1 and f 2 are paired. All these edges are directed except incompatible edges.

Ext: Ext is a subset of vertices called extremities. Any v in Ext is such that either its set of predecessors, or its set of successors is an empty set (according to the overlaps and successive edges).

A: the set of alleles phased with each fact such that A[v] is the set of alleles associated with the fact v ∈ V .

weight: weight associated with each overlaps and links edge. For an overlaps edge, weight[u, v, overlaps] represents the overlap size, i.e. the number of overlaped variants. For a links edge, weight[u, v, links] represents the abundance of the edge, i.e. the number of paired reads associated to the edge.

We will denote by ab[al] the abundance of an allele al ∈ AL where AL is the set of all the alleles. The allele abundance ab[al] is the number of reads containing al.

Note that some facts can be striclty included in others. For example, we can have both following facts in the graph: f 1 = 1h; 2l; 3h; 4l; and f 2 = 2l; 3h;. These facts could not be removed since they can overlaps with other facts. Removing them could delete real path. Let f 3 = 1l; 2l; and f 4 = 3h; 4l; 5l;. There are overlaps between (f 3 , f 2 ), (f 2 , f 4 ) and (f 3 , f 4 ). If f 2 is removed, the relation between f 3 and f 4 is definitively lost.

The figure 7 is an example of 4 haplotypes of 5 variants. Each variant has two alleles. A full example of fact graph over the haplotypes of figure 7 is available in figures 8 and 9.

milp model

Multicommodity-flow based approach

Single genome assembly milp models are based on the longest path problem. The longest path problem is a well-known NP-hard problem. This strategy consists in spreading a flow on an elementary path such that this flow is maximal [START_REF] Andonov | Complete Assembly of Circular and Chloroplast Genomes Based on Global Optimization[END_REF]. According to this, we choose to base our approach on the multicommodity-flow principle [START_REF] Ahuja | Network flows: theory, algorithms, and applications[END_REF]. Multicommodityflow lets spread k flows in a graph, where k is the number of commodities and is fixed by the user. The k flows are optimised simultaneously. Let k be the number of haplotypes belonging to a sample and G = (V, E, Ext, A, weight) be a fact graph. Let's denote by Φ i al a real variable representing the quantity of the haplotype i passing through the allele al (i.e. how many times the haplotype i has been detected traversing the allele al). Multicommodity-flow principle defines that the sum of each commodity's flow is lower-bounded by the abundances of alleles [eq. 1].

∀al ∈ AL, ab[al] ≤ k i=1 Φ i al (1)
To simplify the next constraints, we will give them from a commodity point of view.

Elementary path with abundance

In our context, a commodity can be seen as an elementary path representing a haplotype. The value of the commodity is the haplotype abundance. The haplotype abundance is the number of times this haplotype has been sequenced.

Elementary path

An elementary path is a path that does not pass through the same vertex twice. To model such path, we use gat: a technique described in [START_REF] François | Global optimization for scaffolding and completing genome assemblies[END_REF]. Let's define three binary variables i v , s v and t v denoting the state of the node v in the path, respectively intermediary node, source node and terminal node. gat constraints are the following. A vertex could only have one state [eq. 2].

∀v ∈ V, s v + i v + t v ≤ 1 (2)
There is only one starting vertex [eq. 3] and only one terminal vertex on a path [eq. 4].

v∈Ext s v = 1 ( 3
) v∈Ext t v = 1 (4)
A fact v ∈ V could be a source or a target only if it belongs to Ext [eq. 5].

v∈V \Ext

(t v + s v ) = 0 (5)
Now that vertex states are defined, we must add constraints to link them in an elementary path. Let's denote by x (u,v) a binary variable representing the activation of an edge (u, v, overlaps) ∈ E. An edge (u, v, overlaps) is activated if it is used in the elementary path, i.e. u is a source or intermediary vertex and v is a terminal or intermediary vertex. Finding an elementary path can be perceived as spreading a flow worth 1. Thus, a law of conservation is necessary: a vertex has a single input iff it is terminal or intermediary [eq. 6], and a single output iff it is source or intermediary [eq. 7].

∀v ∈ V, (u,v,overlaps)∈E x (u,v) = i v + t v (6) ∀v ∈ V, (v,w,overlaps)∈E x (v,w) = s v + i v (7)
All these constraints define a path. To define an elementary path, we must remove loops.

To remove loops, we use the Miller-Tucker-Zemlin (mtz) constraint [START_REF] Miller | Integer programming formulation of the traveling salesman problems[END_REF]. mtz constraint consists in labelling each used vertex with a real number. Let's call mtz v this value where v ∈ V . The idea is that each vertex's value is strictly greater than its predecessor's value [eq. 8].

∀(u, v, overlaps) ∈ E, mtz v -mtz u ≥ x (u,v) -|V | × (1 -x (u,v) ) (8)

Path abundance

Each haplotype has an abundance. This abundance can be computed by spreading a flow on the elementary path representing the haplotype. The model must be completed with new constraints. Let f (u,v) be the flow passing through the edge (u, v, overlaps) ∈ E. Let's denote by F F an upper bound for the flow. A flow passes through an edge (u, v, overlaps) ∈ E iff the edge is activated (x (u,v) = 1) [eq. 9].

∀(u, v, overlaps) ∈ E, x (u,v) ≤ f (u,v) ≤ x (u,v) × F F (9)
According to the conservation law, for internal vertices the quantity of the input flow equals the quantity of the exit flow [eq. [START_REF] Laso-Jadart | metaVaR: introducing metavariant species models for reference-free metagenomic-based population genomics[END_REF][START_REF] Medvedev | Computability of models for sequence assembly[END_REF]. No constraint is applied for source and terminal vertices, they are not subject to the conservation law (no inputs or outputs).

∀v ∈ V, (u,v,overlaps)∈E f (u,v) -F F × (1 -i v ) ≤ (v,w,overlaps)∈E f (v,w) (10)
∀v ∈ V,

(v,w,overlaps)∈E f (v,w) ≤ (u,v,overlaps)∈E f (u,v) + F F × (1 -i v ) (11)
Thus, we can formulate the abundance of the path. Let f be a real variable representing the abundance of the path. The path abundance is equal to the outgoing flow of the path's source [eq. 12].

∀v ∈ V, (v,w,overlaps)∈E f (v,w) -F F ×(1-s v ) ≤ f ≤ (v,w,overlaps)∈E f (v,w) +F F ×(1-s v ) (12)

Haplotype specifications

Sequencers and DiscoSNP can provide us complementary information about the haplotypes. This information is useful and needs to be used to prevent the creation of chimeras.

Incompatible edges

Firstly, an incompatible edge (u, v, incompatible) ∈ E represents an incompatible combination of facts: u and v do not belong to the same haplotype [eq. 13].

∀(u, v, incompatible) ∈ E, (s u + i u + t u ) + (s v + i v + t v ) ≤ 1 (13)

Links edges

Secondly, a link edge (u, v, links) ∈ E is a data provided by sequencers. Alleles phased by the paired facts have been seen on the same molecule. That means that u and v belong to the same haplotype. As data provided by sequencers is error-prone, we could not force these links to be satisfied. Let's denote by l (u,v) a binary variable representing the links validation. A link edge is satisfied iff u and v belong to the path [eq. 14] and u is an ancestor of v [eq. 15].

∀(u, v, links) ∈ E, 2 × l (u,v) ≤ (s u + i u + t u ) + (s v + i v + t v ) (14)
∀(u, v, links) ∈ E, l (u,v) -(1 -l (u,v) ) × |V | ≤ mtz v -mtz u (15)

Variants

Each fact phases a set of alleles. The use of alleles' data allows to specify the haplotype. Let a i al be a variable such that a i al = 1 iff the haplotype i passes through the allele al,

otherwise a i al = 0 [eq. 16-17]. ∀i ∈ [1; k], ∀v ∈ V, ∀al ∈ A[v], (s i v + i i v + t i v ) ≤ a i al (16) ∀i ∈ [1; k], ∀al ∈ AL, a i al ≤ v∈V al∈A[v] (s i v + i i v + t i v ) (17)
Moreover, all the alleles must be used by at least one haplotype [eq. 18].

∀al ∈ Al, 1 ≤ i∈[1;k] a i al (18)
Let define Φ i al formulated in section 3.

∀i ∈ [1; k], ∀al ∈ AL, f i -F F × (1 -a i al ) ≤ Φ i al ≤ f i + F F × (1 -a i al ) (19) ∀i ∈ [1; k], ∀al ∈ AL, -F F × a i al ≤ Φ i al ≤ F F × a i al (20)
Note that f i is the abundance of the i th path. This real variable is the extended version of f (see section 3.2.2) for the different commodities.

Objective Function

dna sequencers are error-prone: we are almost sure of what we have seen, but it is possible to have missed some information. Thus, all the haplotype could not be of the same length. The length of the haplotypes should be as close as possible to each other. Let us introduce the real variable δ l modelling this difference. This variable is defined according to the following constraints [eq. 21].

∀i 1 ∈ [1; k], i 2 ∈ [1; k], al∈AL (a i 1 al -a i 2 al ) ≤ δ l (21)
The allele abundance is the number of times that we are sure to have seen the allele. It is possible to have missed some allele occurrences during the sequencing process. Thus, the abundance of haplotypes should be as close as possible to the abundance of alleles. Let δ v al be the difference between the allele's abundance and the flow passing through the allele on each haplotype, where al ∈ AL. This variable is defined by the constraints below [eq. 22].

∀al ∈ AL, (

i∈[1;k] Φ i al ) -ab[al] ≤ δ v al ( 22
)
Note that according to equation 1, δ v al ≥ 0. The objective function of our model is a multi-objective one. It must minimise the difference between alleles abundance and paths flow, and the difference between the length of each haplotype. It must also maximise the number of validated links edges, since it is a powerful information given by the sequencer (see section 3.2.3). To model this function, we must introduce three real variables fixed by the user: α 1 , α 2 , α 3 with α i ≥ 0. These variables are coefficients used to weight the three objectives. Our objective function H k for a given k is as below [eq. 23].

max H k = α 1 × 1≤i≤k (u,v,links)∈E l i (u,v) × weight (u,v) -α 2 × al∈AL δ v al -α 3 × δ l (23)
The solution of the fact graph presented in figure 8 is available in figure 9. The solution satisfied all the constraints presented in this report and the three links edges. It is the solution find by the model.

Results

The model was tested on a simulated dataset. This dataset was obtained by simulating mutations on sars-cov-2 genome (Wuhan strain). This data set is composed of 3 761 facts, 13 847 overlaps and successives, and 36 385 links (see figure 14).

Despite the lack of results showing the assembly quality, the model shows a scalability problem. Small instances of about 100s facts take more than 3 hours to be solved, while we have not been able to process bigger instances in less than 6 hours. For example, we have not been able to process the connected component shown in figure 15.

AlPha-Allele: milp approach for genome haplotyping

To solve the scalability problem of previous approach, we proposed an other model relying on allele graph instead of fact graph. The idea behind this model is that allele graphs are smaller and more regular than fact graphs. The model implementation is still made according to the milp paradigm.

Allele Graph

From DiscoSNP output, we can build an allele graph. An allele graph is a graph modelling the alleles: each vertex models an allele and each edge models the allele order. Let formally define an allele graph G = (V, A, E, F, P ) where:

V: set of variants, it represents the graph columns.

A: set of alleles (v, a) ∈ A, such that v ∈ V is a variant and a ∈ {h, l, -} is one allele of variant v. The abundance of the allele (v, a) will be denoted by ab((v, a)).

E: set of edges, ((v 1 , a 1 ), (v 2 , a 2 )) ∈ E with (v 1 , a 1 ) and (v 2 , a 2 ) two alleles. There is a link between two alleles iff there is at least one fact containing both alleles.

F: set of facts, a fact f is an ordered sequence of alleles. It models the raw facts returned by DiscoSNP. It is a set of tuple, such that (n, (v, al)) ∈ F where n is the fact id and (v, al) ∈ A an allele.

P: set of paired facts, (f 1 , f 2 ) ∈ P such that f 1 and f 2 are two facts. It is a set of couple of facts id. Figure 10 is an example of allele graph. Each vertex is an allele, each graph column is a variant. This graph is constructed with the same input (figure 7) as the fact graph in figure 8. Figure 11 show the original haplotypes on this allele graph.

Properties

Unlike fact graphs, allele graphs have few interesting properties:

1. allele graphs are directed acyclic graph. It is impossible to pass through the same allele twice;

2. allele graphs are n-partite graphs. It is possible to split the vertex in n = |V | parts, such that two vertices of the same parts are not linked. Thus, no elementary path can pass through the same variant twice.

3. allele graphs are almost regular: almost each variant has at most one successor and at most one predecessor.

Linearisation

As previously explained, allele graphs are almost regular. There are few irregularities in the graphs. There are three types of irregularities which are illustrated in figure 3:

Bubbles: a variant has at least two variant successors. It create two concurrent branches in the graph (figure 3-a).

Shortcuts: a variant has at least two variant successors, one of these successors is a descendant of the others variant successors (figure 3-b).

Deadends: one allele of a variant does not have any successors while the others alleles have some (figure 3-c).

Let denote by branches two sequences of variants such that they have the same predecessor and successor, and that there are no edge between element of the two branches.

In order to remove these unwanted structures, we developed a linearisation method. It consists in ordering the variants such that:

• Let v i be a variant and v i+1 one of its successors, thus id(v i ) < id(v i+1 ).

• Let b 1 and b 2 be two concurrent branches, u the variant preceding the two branches and v the variant succeding the two branches. Thus, we must have ∀w ∈ b 1 , id(u) < id(w) < id(v) and ∀w ∈ b 2 , id(u) < id(w) < id(v). • Two distinct variants u, v could not have the same id, id(u) = id(v).

It is impossible to use classical ordering algorithm. Such algorithms could not properly order the variants of concurrent branches. To order these branches, let introduce the notion of dummy alleles. A dummy allele is a fake allele used to model skipped variant.

To order the variant, we apply the following algorithm:

1. Finding the smallest concurrent branches b i and b j . It can be considered as finding the smallest cycle basis of the graph.

2. Ordering the two branches. We choose to order the branches according to the name of the first element of each branch. Let denote by b 1 and b 2 both branches such that b 1 is the fist one.

3. Retrieving the u the branches predecessor and v the branches successor. An example of this algorithm applied to a bubble can be seen figure 4. This algorithm was applied on the allele graph shown in figure 10, the result is shown in figure 12. The graph with the original haplotype (i.e. the solution) is shown in figure 13.

milp Model

Like milp model explained in section 3.2, this new model is based on multicommiditycommodity flow.

Thereafter, we will denote by G = (V, A, E, F, P ) the linearised allele graph and by k the number of haplotypes. The set of linear constraints is divided in 4 categories: haplotype constraints, flow constraints, fact constraints and the objective function.

Haplotype Constraints

Firstly, let's define the constraints defining haplotypes. Let a i (v,al) be a binary variable equaled to 1 iff the allele al of variant v belongs to the i th path. As the allele graph is linear, each haplotype must pass through each variant exactly once [eq. 24].

∀i ∈ [1; k], ∀v ∈ V, (v,al)∈A a i (v,al) = 1 (24)
Let last(V ) be the last variant of the graph, the only variant without any successor, and next(v) for v ∈ V the successor of variant v. Two activated alleles of successive variants are linked by an edge in the allele graph [eq. 25].

∀i ∈ [1; k], ∀v ∈ V \ { last(V )}, ∀(v, al) ∈ A, a i (v,al) ≤ w=next(v) ((v,al),(w,bl))∈E a i (w,bl) (25)

Flow Constraints

Let's now define the abundance of an haplotype, i.e. the flow passing through the commodities. Let f i (v,al) be a real variable modelling the flow of the allele al of variant v on the i th commodity. Let F F be an upperbound of the flow. The flow on inactivated alleles is null [eq. 26].

∀i ∈ [1; k], ∀(v, al) ∈ A, a i (v,al) ≤ f i (v,al) ≤ F F × a i (v,al) (26)
The flow on a commodity is constant between two successive variants [eq. 27].

∀i ∈ [1; k], ∀v ∈ V \ {last(V )}, (v,al)∈A f i (v,al) = w=next(v) (w,bl) f i (w,bl) (27) 4.2

.3 Fact Constraints

A fact, i.e. a tuple of alleles, is satisfied by an haplotype iff all the alleles of the fact belongs to the haplotype. Let q i n be a binary variable modelling the satisfaction of fact n. Fact satisfaction can be modelled by the following linear equation [eq. 28].

∀i ∈ [1; k], ∀(n, (v, al)) ∈ F, q i n ≤ a i (v,al) (28)
Thus, let's denote by q n a binary variable equalled to 1 iff the fact n is globally satisfied. A fact is said globally satisfied iff at least one haplotype satisfied it [eq. 29].

∀(n, _) ∈ F, q n ≤ i∈[1;k] q i n (29)
A paired fact (n 1 , n 2 ) ∈ P is satisfied by an haplotype i iff n 1 and n 2 are satisfied by i. Let p i (n,m) be a binary variable modelling the satisfaction of the paired fact (n, m). Paired facts satisfaction can be defined as below [eq. 30].

∀i ∈ [1; k], ∀(n, m) ∈ P, p i (n,m) ≤ q i n p i (n,m) ≤ q i m (30)
Thus, let's denote by q (n,m) a binary variable equalled to 1 iff the paired fact (n, m) is globally satisfied. A paired fact is said globally satisfied iff at least one haplotype satisfied it [eq. 31].

∀(n, m) ∈ P, p (n,m) ≤ i∈[1;k] p i (n,m) (31)

Objective Function

Multi-objective function are harder to solve and to set up. To stay with a mono-objective function, we choose to force the number of satisfied facts and paired facts. Let α ∈ [0, 1] be the percentage of facts to satisfied and β ∈ [0, 1] the percentage of paired facts to satisfied. Following equations force a certain percentage of facts [eq. 32] and paired facts [eq. 33] satisfaction.

α × |F | ≤ α × n∈F q i n (32) β × |P | ≤ β × (n,m)∈P p i (n,m) (33)
Thus, the only objective is to minimise the gap between the sum of haplotype abundances on each allele and the real allele abundances. Let A r be the set of real alleles, i.e. a) is a dummy allele}. Let define a real value modelling the biggest gap. This variable can defined thanks to following equation [eq. 34].

A r = A \ {(v, al) ∈ A | (v,
∀(v, al) ∈ A r , -≤ i∈[1;k] f i (v,al) -ab((v, al)) ≤ (34)
The objective function can be defined as: minimise

Results

This model was tested on the same dataset as previously (see section 3.3). Despite the lack of results showing the assembly quality, we see an improvement of computation time. Let's see the impact of linearisation on computation time. The computation time for a normal graph and a linearised graph is shown figure 5. We can see that the linearisation process drastically improve the computation time. The number of variants has an impact on the computation time as we can see it in figure 6. This figure compares computation times of linearised and normal graphs according to the number of variants. It shows that normal graph with less than 50 variants can be processed in less than 10 800 seconds. In contrast, linearised graph with 200 variants can be processed in less than 10 800 seconds.

Results validation

At this stage, we do not have any idea on the model assembly quality. The model is still under development but can already prefectly solve toy examples like the one shown in figure 8. Real instances are more complex than the toy example shown in figure 8. An example of a full Covid-19 instance is shown figure 14 with a zoom on one of its connected components in figure 15. This instance is about 500 times bigger than our toy example. The problem is solved separately on each connected component. Thus, we reduce the problem to instances about 20 times bigger than the toy example. To validate the model results on real data sets, we have established an experimental protocol. The objective of this project is to use the reconstructed haplotypes in population genomic. In particular, model results will be used to compute the nucleotide diversity. Nucleotide diversity is a metric used to measure the degree of polymorphism within a population (i.e. the population diversity). Thus, we can validate our results by comparing the nucleotide diversity of our solution with the real nucleotide diversity. Note that this kind of validation could be used only when the original haplotypes are known. We formulate the following results validation protocol.

1. We take the genome of α species. Let [S 1 ... S α ] be these genomes. From each genome S i , we compute β haplotypes [h i,1 ... h i,β ] by mutating the original genome according to a mutation rate ν.

2. Let H = {h i,j | 1 ≤ i ≤ α, 1 ≤ j ≤ β} be the set of all the haplotypes generated for each species. We sample with replacement n haplotypes from H. This last step is repeated γ times. Let [I 1 ... I γ ] be these samples.

3. The original nucleotide diversity value is computed for each sample. Let π i j be the real nucleotide diversity for the i th sample and the j th species.

4. From each sample I i , we simulate illumina reads with SimSeq [START_REF] Benidt | SimSeq: a nonparametric approach to simulation of RNA-sequence datasets[END_REF]. Abundances are generated for each haplotype during the simulation process.

5. Simulated reads are filtered by species thanks to MetavaR [START_REF] Laso-Jadart | metaVaR: introducing metavariant species models for reference-free metagenomic-based population genomics[END_REF]. MetavaR is a de novo algorithm used to cluster reads by species. For each reads, we have a probability of belonging to a cluster. Each cluster is then sent to DiscoSNP [START_REF] Peterlongo | Discosnp++: de novo detection of small variants from raw unassembled read set(s)[END_REF] Figure 6: Solving time in seconds according to the number of variants. Blue line is the solving time for linearised graph. Red line is the solving time for the normal graph. Black line is the time limits, i.e. the process is kill after 10 800 seconds. The time limit is reached for normal graph with at least 50 variants, while it is not reached for linearised graph of 200 variants.

to create fact graphs. Notice that DiscoSNP filters reads to reduce the number of erroneous facts.

6. Finally, we use our model to solve each instance. The nucleotide diversity πi j is computed for each species and compared to the original one π i j . Another solution would be to compare the reconstructed haplotypes to the original ones. However, it is more complicated to have perfectly reconstructed haplotypes. We do not know yet if the reconstructed haplotypes have variant permutations (i.e. the haplotypes are correct within a few variant permutations). These permutations should not impact nucleotide diversity. Thus, we can first study our results with respect to nucleotide diversity before focusing on the quality of sequence reconstruction.

Discussion

In this report, we propose a de novo approach for haplotype-aware genome assembly. In contrast to similar approaches [START_REF] Baaijens | Overlap graph-based generation of haplotigs for diploids and polyploids[END_REF][START_REF] Baaijens | Strain-aware assembly of genomes from mixed samples using variation graphs[END_REF], our model aims to reconstruct all haplotypes thanks to global optimisation. We do not use heuristics to reconstruct haplotypes from the input data. Currently, the number k of haplotypes to reconstruct is fixed by users. We aim to automate this according to the following heuristic: k is the minimum number of haplotypes needed to solve the problem. If k is too small, the model will be infeasible. Thus, finding the minimum k could be done by first trying to solve the model with k = 1.

We then increase k until the model is feasible. Finding that a model is infeasible is really quick and is usually determined during the presolve process. The presolve process is a process during which the problem is simplified before being sent to the solver.

The approach introduced in this report is based on the multi-commodity flow strategy. Almost each model's variable is multiplied k times (where k is the number of commodities). Thus, the computation time will increase as a function of k. In order to reduce the impact of k on the model, we have thought about another (not implemented) approach. This approach relies on extracting a set of paths from the linearised allele graph. Let's see how to compute this set of paths: Remember that each fact and paired fact has a weight.

1. Extract one of the paths that satisfies the most facts and paired facts.

2. Reduce the weight of each satisfied fact and paired fact by one.

3. Remove all the facts and paired facts with a null weight.

4. Restart at step 1 until all facts have been removed. Thus, each extracted paths is a candidate for our assembled haplotypes. In order to filter these paths and select the assembled haplotypes, we proposed to solve a simple linear problem: finding coefficients for each path such that they minimise the difference between haplotypes weights and alleles weights.

Conclusion

In this report, we introduce two milp approaches for the Haplotyping problem. From what we know, we present the first mathematical formulations of the Haplotyping problem. Unlike already existing approaches [START_REF] Baaijens | Overlap graph-based generation of haplotigs for diploids and polyploids[END_REF][START_REF] Baaijens | Strain-aware assembly of genomes from mixed samples using variation graphs[END_REF], ours is not based on heuristics and solve the problem with global optimisation. To realise it, we start by presenting the biological context: genome assembly problem (Section 1.1) and haplotyping problem (Section 1.2). Then, we present both our milp model based on multicommodity-flow [START_REF] Ahuja | Network flows: theory, algorithms, and applications[END_REF] (Sections 3 and 4). This model aims to reconstruct highly reliable haplotypes fragments that will be used in population genomic. Finally, we finish with the presentation of our results validation protocol (Section 5) and a discussion (Section 6).

We currently only have results on toy examples and we do not have any on real data. However, we provide computation times for several real instances. This data shows that the computation time is better for the alleles approach with linearised graph. Result quality has not been studied. However, we proposed a full experimental protocol that will be used to validate the model results.

We introduce a promising new milp approach for haplotyping and the first milp formulation of this problem. However, further tests need to be done and the experimental protocol must be implemented. Moreover, this stage the number of haplotypes to reconstruct is fixed by users. We aim to automate this task. In the discussion (Section 6), we also proposed another (not implemented yet) approach allowing to divide by k the number of model's variables (where k is the number of haplotypes to reconstruct). The figure on the left is a visual representation of the 4 haplotypes. Each colored line represents a haplotype. Straigh lines mean that the haplotype passes through the variant, dotted lines mean that the haplotype does not pass through the variant. For example, the haplotype H1 (blue line) passes through the following alleles: 1h;2h;3l;4l;6h. It does not pass through the variant 5.

In this example, haplotypes pass either through the variant 4 or through the variant 5. On the top right is a summary table of the haplotypes and their abundances. On the bottom right is a summary of the allele abundances. The abundance of an allele is the sum of the abundance of the haplotypes passing through it. The abundance of a variant is the sum of its allele abundances. This instance is about 20 times bigger than the toy example presented in figure 8.

Figure 1 :

 1 Figure 1: Example of sequencer output. It contains 100 millions of reads. The lines starting with '>' are comments provided by the sequencers, here it only contains the identifier of the read (the line below). The read is the sequence of 'A', 'C', 'G' and 'T'.

Figure 2 :

 2 Figure 2: Illustration of the phasing problem: facts and paired facts construction. In the top figure, 'T'-'A' and 'C'-'G' are two variants such that 'T' and 'A' are the two alleles of the first variant, while C and G are the two alleles of the second variant. Red reads phase 'T' and 'C' into a fact. Purple reads phase 'A' and 'G' into another fact. In the bottom figure, there is a new variant 'C'-'A' (with 'C' and 'A' being its alleles). The dot edges between the reads represent paired reads. Red reads phase 'T' and 'C', and their paired reads phase 'C'. The facts formed by these paired reads are paired.

2 . 1 .

 21 It is a real variable representing the flow passing through the allele al ∈ AL in the haplotype i ∈ [1; k]. The flow passing through an allele al is equal to the flow of the path i if this path passes through al [eq. 19-20].

Figure 3 :

 3 Figure 3: a) Bubble example. Purple nodes are the bubble extremities. The left one has two variant successors which create two branches (blue one and red one). These branches merge into the right purple node. b) Shortcut example. Left purple node has two variant successors. One of these successors is the right purple node which is also a successor of the blue branch. c) Deadend example. Purple node does not have any successors while others nodes of the column have some.

4 . 6 .

 46 Removing all the edges between u and the first element of b 2 , and the last element of b 1 and v. 5. Adding a dummy allele to each variant of b 1 and b 2 . Adding edge between the dummy alleles of two consecutives variants of each branch. Adding edges between u and the first dummy alleles of b 1 , the last dummy allele of b 1 and all the first non-dummy alleles of b 2 , all the last non-dummy alleles of b 1 and the first dummy allele of b 2 , and the last dummy allele of b 2 and v.

Figure 4 :

 4 Figure 4: Example of linearisation on a bubble. The two blue blocks are the concurrent branches Each block is a variant and each circle is an allele. Yellow circles are dummy alleles.

Figure 5 :

 5 Figure 5: Solving times in seconds for each instance. Blue stack is the computation time of the linearised instance. Red stack is the computation time of the normal instance.

Figure 7 :

 7 Figure 7: Example of 4 haplotypes of 6 variants.The figure on the left is a visual representation of the 4 haplotypes. Each colored line represents a haplotype. Straigh lines mean that the haplotype passes through the variant, dotted lines mean that the haplotype does not pass through the variant. For example, the haplotype H1 (blue line) passes through the following alleles: 1h;2h;3l;4l;6h. It does not pass through the variant 5.In this example, haplotypes pass either through the variant 4 or through the variant 5. On the top right is a summary table of the haplotypes and their abundances. On the bottom right is a summary of the allele abundances. The abundance of an allele is the sum of the abundance of the haplotypes passing through it. The abundance of a variant is the sum of its allele abundances.

Figure 8 :

 8 Figure 8: Example of fact graph obtains from the haplotypes shown in figure 7. On the left is the set of alleles associated with each fact. On the right is the fact graph. The nodes are the facts and blue nodes are the extremities. Black arrows are the overlaps edges and green dotted lines are the links edges. The value over each edge is its weight. The incompatibles edges are not shown on the figure due to visibility purposes. There are no successive in this example.

Figure 9 :

 9 Figure 9: Fact graph presented figure 8 with the reconstructed haplotypes of figure 7. The four links edges are satisfied: (1, 12) is satisfied by H4, (1, 4) by H1, (6, 7) by H2 and (8, 9) by H3.

Figure 10 :

 10 Figure 10: Example of allele graph obtains from the haplotypes shown in figure 7. On the left is the set of alleles associated with each fact. On the right is the allele graph. Each node is an allele (blues nodes are extremities alleles). All the nodes of each column belong to the same variant. Black arrows represent the allele succession according to the facts. Notice that a path could either passes by the variant 4 or by the variant 5. On the bottom-left is the set of paired facts.

Figure 11 :

 11 Figure 11: Allele graph presented figure 10 with the reconstructed haplotypes of figure 7. The four links edges are satisfied: (1, 12) is satisfied by H4, (1, 4) by H1, (6, 7) by H2 and (8, 9) by H3.

Figure 12 :

 12 Figure 12: Example of simplified allele graph obtains from allele graph shown in figure 10. On the left is the set of alleles associated with each fact. On the right is the allele graph. Each node is an allele (blues nodes are extremities alleles). All the nodes of each column belong to the same variant. Black arrows represent the allele succession according to the facts. Unlike the allele graph shown figure 10, each path passes through each variant. The alleles 4-of variant 4 and 5-of variant 5 are ghost alleles. Ghost alleles are used to allow a path to skip unwanted variants while still passing through each variant. On the bottom-left is the set of paired facts.

Figure 13 :

 13 Figure 13: Simplified allele graph presented figure 12 with the reconstructed haplotypes of figure 7. The four links edges are satisfied: (1, 12) is satisfied by H4, (1, 4) by H1, (6, 7) by H2 and (8, 9) by H3. Unlike the solution shown figure 11, each haplotype passes through each variant.

Figure 14 :

 14 Figure 14: Example of a fact graph of a real data set: 10 haplotypes of sars-cov-2.This fact graph contains 3 761 facts (including 435 extremities) divided into 25 connected components. The facts are the blue squares, while the green squares are the extremities. These facts phase 1 860 alleles. There are 13 257 overlaps edges (blue edges), 590 successive edges (yellow edges), 32 295 incompatible edges (not shown for visibility purposes) and 36 385 links edges (not shown for visibility purposes). This instance is about 500 times bigger than the toy example instance presented in figure8.

Figure 15 :

 15 Figure 15: Zoom on one connected component of the sars-cov-2 fact graph. This connected component is an average one. It contains 210 facts (including 18 extremities). The facts are the blue squares, and the extremities the green squares. There are 677 overlaps edges (blue edges), 39 successive edges (yellow edges), 1 447 incompatible edges (not shown for visibility purposes) and 1 521 links edges (not shown for visibility purposes).This instance is about 20 times bigger than the toy example presented in figure8.

It does not use any reference genome to perform the assembly.

The fasta format is used to stock nucleotide sequences

Acknowledgement

I would like to thank Pierre Peterlongo and Rumen Andonov for giving me this project and for helping me to realise it. I would like to also thank Amin Madoui who helped me to understand the biological base needed for this work. Finally, I thank the whole GenScale team for giving me the opportunity to present it to roadef2020.