Mathieu Debure
email: m.debure@brgm.fr

Pascal Andreazza

Aurélien Canizarès

Sylvain Grangeon

Catherine Lerouge

Paul Mack

Benoît Madé

Patrick Simon

Emmanuel Veron

Fabienne Warmont

Marylène Vayer

Study of iron-bearing dolomite dissolution at various temperatures: Evidence for the formation of secondary nanocrystalline iron-rich phases on the dolomite surface

We investigated the dissolution of a natural Fe-containing dolomite [Ca 1.003 Mg 0.972 Fe 0.024 Mn 0.002 (CO 3 ) 2 ] under acidic conditions (pH 3-5.5) with atomic force microscopy (AFM) at 20 °C, and with batch dissolution experiments at 80 °C. Dolomite dissolution proceeded by two identified mechanisms: removal of dolomite layers through spreading and coalescence of etch pits nucleated at defect points, and stepped retreat from surface edges. The dolomite dissolution rate increased when pH decreased (from 0.410 nm s -1 at pH 3 to 0.035 nm s -1 at pH 5). Rates calculated from edge retreat (v edges ) and from etch-pit spreading rates (v sum ) were consistent; the etch-pit digging rate was almost 10 times slower than its spreading rate. Nanocrystalline secondary phases precipitated in the course of dolomite dissolution at pH 3 and 80 °C were identified as (nano)hematite, ferrihydrite and an ankerite like mineral using X-ray diffraction, transmission electron microscopy, Raman and Xray photoelectron spectrometry. In addition, Mg enrichment of the surface layer was observed at 80 °C, due to preferential release of Ca in solution. The characterizations performed at a nanocrystalline scale highlighted the role played by impurities in the dolomite dissolution/precipitation scheme and proved that two mechanisms explain the incongruent dolomite dissolution: secondary phase precipitation and preferential release of Ca over Mg.


Keywords: genome assembly, scaffolding, contig, longest simple weighted path problem, integer programming

This work focuses simultaneously on both the scaffolding and gap filling phases of de nouveau genome assembly. Given a set of contigs and their relationships-overlaps and/or remoteness in terms of distances between them-we propose an optimization-based approach for finding the genome sequence as a longest sequence that is consistent with the given contig and linkage information. Specifically, we define a graph, which we call a contig graph, that encodes information about contigs and overlaps and mate-pair distances between them, and reduce the scaffolding problem to the problem of finding a longest simple path in that graph such that as many as possible mate-pairs distances are satisfied. Since both conditions cannot generally be simultaneously satisfied, our objective function is a linear combination of the path length and penalties for distance mismatches. Unlike the shortest path problem with non-negative weights, for which efficient polynomial-time algorithms exist, the longest weighted path problem is NP-hard. We solve this problem by reformulating it as a mixed integer linear program (MILP) and develop a method that exactly solves the resulting program on genomes of up to 165 contigs and up to 6682 binary variables. An advantage of our approach is that the modeling of scaffolding as a longest path problem allows one to solve simultaneously several subtasks specific for this problem like: contig orientation and ordering, repeats, gap filling, and scaffold extension, which in other approaches are targeted as separate problems. We are not aware of previous approaches on scaffolding based on the longest path problem reduction. A drawback of the typically used strategy of constructing a set of disjoint paths, rather than a single path, is that it would require additional steps of gap filling and scaffold extension, involving additional work. Moreover, it would make impossible to find a provably optimal final solution, since, even if each separate problem is implemented optimally, their combination may not be optimal. We tested this model on a set of chloroplast and bacteria genome data and showed that it allows to assemble the complete genome as a single scaffold. Compared to the publicly available scaffolding tools that we have tested, our solution produces assemblies of significantly higher quality.

Optimisation globale appliquée à l'assemblage génomique

Résumé : Ce travail porte particulièrement sur les étapes de scaffolding et de gap-filling lors de l'assemblage de novo d'un génome. Etant donné un ensemble de contigs et leurs relations (chevauchement et/ou éloignement), nous proposons une approche basée sur l'optimisation, qui consiste à trouver le plus long assemblage possible consistant avec les distances attendues. Nous résolvons ce problème en le formulant sous la forme d'une programme mixte en nombres entiers (MILP). Notre méthode permet de résoudre exactement des instances comptant jusqu'à 6682 variables binaires pour 165 contigs.

L'un des avantages de notre modélisation est qu'elle permet de résoudre simultanément différentes sous taches spécifiques à l'assemblage complet d'un génome : déterminer l'orientation de chaque contig, ordonnancer les répétitions, compléter les scaffolds et les étendre en un assemblage unique etc. . . Ce qui constitue d'ordinaire un ensemble de problèmes à traiter séquentiellement. (Nécessitant alors des solutions spécifiques pour chacun d'entre eux). A notre connaissance, aucune autre approche ne propose de réduire le problème de l'assemblage génomique à celui de la recherche du plus long chemin.

Mots-clés : Pas de motclef 1 Introduction

Modern Next-Generation Sequencing (NGS) techniques are not able to output the whole genome sequence in one large sequence, but instead output billions of short DNA sequences, called reads. These reads are extremely redundant, erroneous, and, consequently, unusable practically. Genome assembly is the challenging computational task, consisting in reconstructing the full genome sequence from these fragmented raw data. This is a complex procedure, usually composed of three main steps: (1) generation of contigs, which are long consensual DNA sequences issued from the overlapping of the reads; (2) orientation and ordering of these contigs, called the scaffolding; (3) gap-filling and contig extension. A notable weakness of the contemporary approaches for de-novo assembly is to consider the above process as a set of independent tasks and not be able to propose a global optimal solution.

The first step generates a list of contigs (assembled ungapped sequences) that represent the "easily assembled regions" of the genome. Building contigs is currently supported by methods using a specific data structure called de-Bruijn graph [?].

Whereas the main challenge in the first step is in the sheer size of the raw data, in the second step, scaffolding, the data is of moderate size, but the problem remains largely open because of its NP-hard complexity [?]). The goal here is to provide a reliable order and orientation of the contigs in order to link them together into scaffolds-a sequence of contigs that overlap or are separated by gaps of a given length. These gaps are generated during sequencing based on pairedend or mate pair reads [?], and have distance information associated with them. Specifically, they can be represented as couples of fragments separated by a known distance (called insert size) and provide information about distance between contigs.

While the ultimate goal of the genome assembly is to generate a complete genome, the scaffolding phase usually produces a set of multiple scaffolds that, in addition, may contain inside them regions that have not been completely predicted. Further steps, such as gap-filling and contig extension are generally required to complete the genome.

This paper focuses on the scaffolding phase. Given a set of contigs and their relationshipsoverlaps and/or remoteness in terms of distances between them (insert sizes)-we propose an optimization-based approach for finding the genome sequence as the longest sequence that is consistent with the given contigs and linkage information. Specifically, we define a graph, which we call contig graph, that encodes information about contigs and distances between them, and reduce the scaffolding problem to finding a longest simple path in that graph such that as many as possible mate-pairs distances are satisfied (we call hereafter such path just a longest path). Since both conditions cannot generally be simultaneously satisfied, our objective function is a linear combination of them.

Unlike the shortest path problem with non-negative weights, for which efficient polynomialtime algorithms exist, the longest weighted path problem is NP-hard [?], which means that no polynomial time solution is likely to exist for it. We solve this problem by reformulating it as a mixed integer linear program (MILP) and develop a method that exactly solves the resulting program on genomes of up to 165 contigs and up to 6682 binary variables. We analyze the performance of the algorithm on several chloroplast genomes and compare it to other scaffolding algorithms.

Most previous work on scaffolding is heuristics based, e.g., SSPACE [?], GRASS [?], and BESST [?]. Such algorithms may find in some cases good solutions, but their accuracies cannot be guaranteed or predicted. In contrast, our method always finds a longest path in the contig graph. There is no guarantee that the genome sequence corresponds to a longest path, but our experiments show that that is the case in many instances or, if not, there is a very small difference between the two. Exact algorithms for the scaffolding problem are presented in [?],

Inria but the focus of that work is on finding structural properties of the contig graph that will make the optimization problem of polynomial complexity. In [?], integer linear programming is used to model the scaffolding problem, with an objective to maximize the number of links that are satisfied. In contrast, our objective is to maximize the length of the resulting scaffold. Moreover, we aim at producing a single path or cycle, rather than a set of paths and cycles. While our focus is on accuracy, [?] focuses on efficiency, and indeed their algorithm, being a kind of heuristics, is faster than ours.

The contributions of this study are as follows:

• Our modeling of the scaffolding problem as a longest path problem allows to solve simultaneously the set of subtasks specific for this problem like: contigs orientation and ordering, repeats, gap filling and scaffold extension, which in other approaches are separate phases.

• The scaffolding problem is reduced to finding a longest paths in a particular graph. In addition, these paths need to satisfy a set of distances between couples of vertices along these paths. We are not aware of previous approaches on scaffolding based on the longest path problem. The disadvantage of having a set of disjoint paths, rather than a single path, is that it would require additional steps of gap filling and scaffold extension, involving additional work. Moreover, it would make it impossible to produce an optimal final solution, since, even if each separate problem is implemented optimally, their combination may not be optimal.

• We formulate the above problem as a mixed integer linear program (MILP) with several interesting properties like: cycles elimination constraints, using binary variables for the edges of the graph only. Vertices are modeled with real variables, but we prove that the integrality of these variables follows from other constraints. Moreover, the commonly used approach for solving the longest weighted simple path when the initial (source) and final (target) vertices are unknown consists in artificially adding these two vertices and 2|V | edges in the graph G = (V, E) [?]. This increases by 2|V | the number of binary variables, which is a drawback when the density of the graph is small (as is the case of scaffolding graph).

In contrast, our modeling does not require such a graph transformation and requires fewer binary variables.

• We tested this model on a set of chloroplast and bacteria genome data and showed that it allows to assemble the complete genome as a single scaffold. None of the publicly available scaffolding tools that we have tested targets single scaffolds (this is corroborated by the obtained numerical results).

• Our numerical experiments indicate that the relaxation of the mixed integer model is tight and produces upper bounds of excellent quality. This suggests a promising direction of research towards the scalability of our approach.

2 Modeling the scaffolding problem

Graph Modeling

We model the problem of scaffolding as path finding in a directed graph G = (V, E) that we call a contig graph, where both vertices V and edges E are weighted. The set of vertices V is generated based on the set C of the contigs according the following rules: the contig i is represented by at least two vertices v i and v i (forward/inverse orientation respectively). If the contig i is repeated k i times, it generates 2k i vertices. The values of k i are computed in the contigs-generation phase. Denote N = i∈C k i , therefore |V | = 2N . The edges are generated following given patterns-a set of known overlaps/distances between the contigs. Any edge is given in the graph G in its forward/inverse orientation. We denote by e ij the edge joining vertices v i and v j and the inverse of edge e ij is e j i . For any i, the weight w i on a vertex v i corresponds to the length of the contig i, while the weight l ij on the edge e ij corresponds to the value of the overlap/distance between contigs i and j. The problem then is to find a path in the graph G such that the total length (the sum over the traversed vertices and edges) is maximized, while a set of additional constraints are also satisfied:

• For any i, either vertex v i or v i is visited (participates in the path).

• The orientations of the nodes does not contradict the constraints imposed by mate-pairs. This is at least partially enforced by the construction of G.

To any edge e ∈ E we associate a variable x e . Its value is set to 1, if the corresponding edge participates in the assembled genome sequence (the associated path in our case), otherwise its value is set to 0. There are two kinds of edges: edges corresponding to overlaps between contigs, denote them by O (from overlaps), and edges associated with mate-pairs relationships, denote them by L (from links). We therefore have E = L∪O. Let l e be the length of the edge e = (u, v). We have l e < 0 and |l e | < min {w(u), w(v)}, ∀e ∈ O, and l e > 0 ∀e ∈ L. Let w v be the length of the contig corresponding to vertex v and denote W = v∈V w v .

Let

A + (v) ⊂ E (resp. A -(v) ⊂ E
) denote the subset of arcs in E leaving (resp. entering) node v.

Mixed Integer Linear Programming Formulation

We associate a binary variable for any edge of the graph, i.e. ∀e ∈ O : x e ∈ {0, 1} and ∀e ∈ L : g e ∈ {0, 1}.

(

Furthermore, to any vertex v ∈ V we associate three variables, i v , s v , and t v , which stand respectively for intermediate, source, and target for some path, and satisfy

0 ≤ i v ≤ 1, 0 ≤ s v ≤ 1, 0 ≤ t v ≤ 1. (2)
All three variables are set to zero when the associated vertex v participates in none of the paths. Otherwise, v can be either a source/initial (noted by

s v = 1, t v = 0, i v = 0), or a target/final (t v = 1, s v = 0, i v = 0)
, or an intermediate vertex, in which case the equalities i v = 1, t v = 0 and s v = 0 hold. Moreover, each vertex (or its inverse) can be visited at most once, i.e.

∀(v, v ) :

i v + i v + s v + s v + t v + t v ≤ 1. (3)
The four possibles states for a vertex v (to belong to none of the paths, or otherwise, to be a source, a target, or an intermediate vertex in some path) are provided by the following two constraints

s v + i v = e∈A + (v) x e , t v + i v = e∈A -(v) x e . (4)
Finally, only one sequence (a single path) is searched for

v∈V s v = 1 and v∈V t v = 1. (5)
Inria Theorem 1. The real variables i v , s v , t v , ∀v ∈ V take binary values.

We introduce a continuous variable f e ∈ R + to express the quantity of the flow circulating along the arc e ∈ E ∀e ∈ E : 0 ≤ f e ≤ W x e .

For e ∈ O, the value of x e is set to 1, if the arc e carries some flow and 0, otherwise. In other words, no flow can use the arc e when x e = 0.

We use the flows f e in the following constraints, ∀v ∈ V ,

e∈A -(v) f e - e∈A + (v) f e ≥ (i v + t v )(w v + e∈A -(v) l e x e ) -W s v , W s v ≤ e∈A + (v) f e . ( 7
)
The purpose of the last two constraints is manifold. When a vertex v is a source (s v = 1), (??) generates and outputs from it an initial flow of sufficiently big value (W is enough in our case). When v is an intermediate vertex (i v = 1), constraint (??) forces the flow to decrease by at least l (u,v) + w v units when it moves from vertex u to its adjacent vertex v. The value of the flow thus is decreasing and this feature forbids cycles in the context of (??). When v is a final vertex, (??) is simply a valid inequality for the input flow.

We furthermore observe that because of (??), the constraint (??) can be written as follows

∀v ∈ V : e∈A -(v) f e - e∈A + (v) f e ≥ (i v + t v )w v + e∈A -(v) l e x e -W s v . (8)
The constraint (??) is linear and we keep it in our model instead of (??). Furthermore, binary variables g e are associated with links. For (s, t) ∈ L, the value of g (s,t) is set to 1 only if both vertices s and t belong to the selected path and the length of the considered path between them is in the given interval [L (s,t) , L (s,t) ]. Constraints related to links are :

g (s,t) ≤ s s + i s + t s and g (s,t) ≤ s t + i t + t t (9)
∀(s, t) ∈ L :

e∈A + (s) f e - e∈A -(t) f e ≥ L (s,t) g (s,t) -M (1 -g (s,t) ) (10)
∀(s, t) ∈ L :

e∈A + (s) f e - e∈A -(t) f e ≤ L (s,t) g (s,t) ) + M (1 -g (s,t) ), (11)
where M is some big constant. We search for a long path in the graph and such that as much as possible mate-paired distances are satisfied. The objective hence is :

max e∈O x e l e + v∈V w v (i v + s v + t v ) + p e∈L g e ( 12
)
where p is a parameter to be chosen as appropriate (currently p = 1).

Computational results

Here we present the results obtained on a small set of chloroplast and bacteria genomes given in Table ? parameter unitig instead of contig (a unitig is a special kind of a high-confidence contig). Based on these unitigs, the overlaps between them, as well as the mate-pair distances, we generated a graph as explained in Section ??.

Our results were obtained on an Intel(R) Xeon(R) CPU E5-2670 v2 @ 2.50GHz with 20 cores, 64 GB of RAM, and using Gurobi solver for solving the MILP models. We compared our results with the results obtained by three of the most recent scaffolding tools-SSPACE [?], BESST [?], and Scaffmatch [?]. In order to evaluate the quality of the produced scaffolds, we applied the QUAST tool [?]. The results are shown on Table ??. We observe that our tool GST (from Genscale Scaffolding Tool) is the only one that consistently assembles the complete genome (an unique scaffold in #scaffolds column) with more than 98% (and in four cases at least 99.9%) correctly predicted genome fraction and zero misassembles.

Our next computational experiments focussed on comparing various relaxations and other related formulations for the MILP model described in the previous section. Let us denote in the sequel by BR (Basic Real) the model defined by the linear constraints (??-??) and (??-??) and objective function (??). Let BB (from Basic Binary) denote the same model except that constraint (??) is substituted by its binary variant, i.e. ∀v ∈ V : i v ∈ {0, 1} and s v ∈ {0, 1} and t v ∈ {0, 1}.

According to Theorem ??, the models BB and BR are equivalent in quality of the results. We are interested here in their running time behavior. We study as well the linear programming relaxation of BR (denoted by BR LP ) where the binary constraints (??) are substituted by ∀e ∈ O : 0 ≤ x e ≤ 1 and ∀e ∈ L : 0 ≤ g e ≤ 1.

Furthermore, let us omit from the objective function the last term that is associated to matepairs, as well all constraints in BR model that relate to to mate-pairs distances (i.e. constraints (??, (??), (??)). We call this model LP (Longest Path) since it simply targets finding the longest path in the contig graph. Any solution of this model can be extended to a solution of BR model by a completion g e = 0 ∀e ∈ L. Its optimal value yields a lower bound for the main model BR.

By experimentally analyzing these models (see Table ?? in the Appendix for the data) we observe that, as expected, the model BR significantly outperforms BB in running time. The upper bounds computed by the linear relaxation BR LP are extremely close to the exact values computed by BR model. Furthermore, the quality of the lower bound found by the longest path approach LP is also very good. Interestingly, we observed for the given instances that this value is close to the genome's size. We presume that the LP model can be used for predicting the length of the genome when it is unknown.

Conclusion

We developed and tested algorithms for scaffolding based on a version of the longest path problem and MILP representation. Our algorithms significantly outperform three of the best known scaffolding algorithms with respect to the quality of the scaffolds. Regardless of that, we consider the current results as a work in progress. The biggest challenge is to extend the method to much bigger genomes. We plan to use some additional ideas and careful implementation to increase the scalability without sacrificing the accuracy of the results. Table 3: Performance of the basic model, some of its relaxations and related formulations. The symbol * indicates that the execution has been stopped by time limit.

B Computational results -model comparison

Inria

Table 1 :

 1 ?. Synthetic sequencing reads have been generated for these instances applying ART simulator[?]. For the read assembly step required to produce contigs we applied minia [?] with Scaffolding datasets.

	Datasets	Total	#unitigs	#nodes #edges	#mate-pairs
	Acinetobacter	3 598 621	165	676	8344	4430
	Wolbachia	1 080 084	100	452	7552	2972
	Aethionema Cordifolium	154 167	83	166	898	600
	Atropa belladonna	156 687	18	36	114	46
	Angiopteris Evecta	153 901	16	32	144	74
	Acorus Calamus	153 821	15	30	134	26

Table 2 :

 2 Performance of different solvers on the datasets from Table??. GST is our tool.

	RR n°9050

RR n°9050

Appendix

A Proof of Theorem ??

Proof. Let us analyse the four possibles cases deduced from (??). Denote S + = e∈A + (v) x e and S -= e∈A -(v) x e . i) S + = 0 and S -= 0.

In this case it follows from (??) that

The above is equivalent to

iii) S + = 1 and S -= 0.

The above is equivalent to s v + i v = 1 and t v + i v = 0, which leads to s v -t v = 1. Hence, from (??), we have t v = 0 and therefore s v = 1 and i v = 0. iv) S + = 0 and S -= 1.

This case is analogous to iii) and we have s

Contents

RR n°9050

RESEARCH CENTRE RENNES -BRETAGNE ATLANTIQUE