these simple structure spaces, that a global optimum is computable in polynomial time for any dataset [7]. Due to this latter property, it is necessary to introduce variance in the dataset for each component learning task, in order to converge to diverse simple models, allowing different significant features of the final model to be represented in the population of simple models. This variance is obtained by: 1) perturbing the original dataset through sampling with replacement; 2) forbidding the use of some edges in the spanning arborescence, the edge blacklist being randomly generated for each component.

The paper is organized as follows. We first provide the background material through the description of required theory and previous works in section 2. Then, we describe our mixture algorithm in section 3. We finally validate in section 4 our algorithm with promising results on experiments over the challenging biological network inference application from time series data from the popular DREAM D8C1 recent challenge [8], before discussing the algorithm and its perspectives in section 5.

## 2 Related Works

Network inference has been studied for decades in bioinformatics context and many solutions can be grouped into the following families.

Feature selection methods break the global network inference problem into a set of *N* subproblems, *N* being the number of features in the studied biological system. The objective is to find the best neighbors of each feature in the biological model to learn, using various methods. A main strategy is to learn and rank, for each feature *j*, a set of (possibly linear) regression coefficients describing the function of all other features to *j*. Recent examples of this strategy include: the HLICORN method [9,10], where linear coefficients are computed between each possible gene and previously discovered coregulator sets; the GENIE3 [11] algorithm, where a supervised classification tree is built for each node, seen as the target, and potential parents are ranked according to their decreasing order of entropy gain between the tree's layers; BORUTA [12], a wrapper type of method around random forest classification; TIGRESS [13], which uses a *perturb and combine* strategy as this paper does. Note that this latter strategy is still significantly different from our proposal. First, TIGRESS reduces the learning problem into an independent learning problem for each node, while our proposal constrains the learned models globally, allowing here to enforce global sparsity and more generally allowing rich extensions to the method (e.g. global prior addition). Furthermore, the models to combine in this work are obtained by Lasso regression and forward selection strategy, without any quality guarantee, unlike our strategy. Finally, the approach only introduces variance in the process via data perturbation, while we propose both data and edge space sampling for each component, the latter showing significant experimental impact (cf. Section 4).

Ordinary Differential Equations (ODE) methods aim at modelling a dynamic biological system as a set of differential equations involving its features (genes, proteins, etc.). In practice, these methods can find highly accurate parameters for given equations, and have even been extended for learning the equations themselves [14], but their computational cost make them only suitable for small biological systems, which is not our target.

Pairwise score approaches aim at finding networks optimizing a sum of feature to feature scores measuring the level of correlation or dependency between them. A famous algorithm in this category is the ARACNE [15] method which computes a superset of the Chow-Liu spanning tree induced by the learning dataset, by first computing a complete graph with all pairwise mutual informations (MI), and then iteratively flagging edges for removal using a MI inequality for every triangle in the graph. Such posterior optimization, called *transitive reduction*, is indeed often very important in algorithms of this family since distinguishing between direct or indirect dependencies is not necessarily possible with the used measures.

(Dynamic) Bayesian network (DBN) learning algorithms [16] aim at finding the best factorization of the joint distribution involving the features at different consecutive time steps, as a product of conditional distributions (one per feature). Most DBN learning algorithms targetting problems of any size try to find a local optimum in the model space, through the use of a heuristic, going from one model in the space to another through iterative local perturbation of a best candidate obtained so far (best-first strategy) [17] or a set of best candidates (beam search) [18]. These algorithms often add some extra-mechanics to proceed further the first local optima, as for example random restarts or tabu search [19].

Most methods in the above mentioned families consider the input dataset as is to make network inference. In the context of small datasets, the data fragmentation issue together with the possibility of noise presence in
the dataset can misconduct the learning algorithms. In order to abstract from a single model learned from such dataset, ensemble learning strategies have been designed to implement a "wisdom of crowds" principle in a machine learning context. Historically used in a fully supervised context to average classification predictions and turn them into final majority decisions, these principles have been transposed for probability density estimation and structure learning contexts, a good example of such strategy for probabilistic models having been theorized under the Bayesian Model Averaging (BMA) [4] framework.

In BMA, the objective is to find an expected model, defined by a set of *expected features* it must satisfy, such as the dependency between its variables, the underlying graph edges or paths and so forth, assuming those features are relevant in the model space, by combining different *component* models. In practice, considering the whole space is intractable and one must rely on a subset to approximate the result, such as using Markov Chain Monte Carlo methods [20] or bagging [21].

Typical BMA methods consider the whole model space for learning which leads to the question of which heuristic to use and with which extra mechanics. This problem has been recently partially tackled in [5] but choices still remain. A recent work by Schnitzler et al. [22] considers combining more simple component models instead, namely spanning trees, and proposes extensions of the Chow-Liu algorithm [23] which show good results in their non-dynamic density estimation context. The algorithm proposed in this paper is close but differs in two main ways: it focuses at outputing an expected structure as in traditional BMA for structure learning, and thus require to make choices between the different component properties, while their work outputs a result in the form of a linear combination of each component's density probability function; it focuses on time series datasets instead of i.i.d. propositional data. Note that a preliminary adaptation of Schnitzler work for structure learning has been proposed in [24], but authors use the component models as an intermediate product to refine potential parent candidates of a greedy algorithm, while we propose an algorithm which directly infers a model from the components.

Among BMA methods, it is important to mention the seminal work in [4] about BMA structure learning through the computation of expected features using bagging on the Bayesian networks space. However, in this work, learned expected features are used as constraints for heuristic learning, which does not include component edges themselves, and whose result does not lead to significant improvement compared to the unconstrained learning counterpart. Finally, in this paper, while most work focuses on combining unconstrained Bayesian networks, few indications are also given on the results obtained by merging spanning trees, which apparently gave worse results than with unconstrained Bayesian networks. However, this can be explained by the low introduction of variance, unlike our proposal, which can lead to poor space exploration because of the components global optimality. Also, chosen experiment datasets can not be considered as being high dimensional ones and the variables/individuals ratio is less subject to data fragmentation, which can explain that more biased models perform worse. We will demonstrate in section 4 that combining arborescences in an adequate way can actually give very good results in this context.

# 3 Combining Spanning Arborescences for Network Inference

In this section, we propose a general algorithm for network inference from the combination of spanning arborescences. Several application oriented details are willingly not given, such as the specific scoring function used and the value of hyperparameters. These informations are detailed in Section 4 for the problem at hand.

## 3.1 Data representation

Let us consider a matrix representation of a dataset *D* consisting of *n* ordered time stamps (*D* rows) over *N* variables (*D* columns). Each column *j* is a sequence describing an observed variable over *n* time steps ⟨*D*<sub>*ij*</sub>⟩<sub>*i*∈{1,...,n}</sub> and each row *i* describes the state of a system at time step *i* over the *N* observed variables. Our goal is to find a model of the system of interest in terms of dependencies between the observed variables at different time steps of the system. In this paper, we assume that this system is a Markov process, i.e. that each time step state only depends on the immediate previous step state, and that the transition from each state to the next state is driven by the same underlying model.

The first step for the proposed algorithm to work is to transform the $n \times N$ dataset into a $(n-1) \times 2N$ dataset $D^t$ describing 2 consecutive time slices of the system. The transformation consists in concatenating
every pair of consecutive time steps from *D* into a "dynamic" example in *D*<sup>*t*</sup>, i.e., we have for each row *D*<sub>*i*</sub><sup>*t*</sup>:

$$D_i^t = [D_{i,1} \dots D_{i,N} D_{i+1,1} \dots D_{i+1,N}].$$

The given assumptions and the consequent representation allow transforming the ordered structure learning problem from *D* into a simpler i.i.d. structure learning problem in the dynamic examples space of *D*<sup>*t*</sup>.

## 3.2 Learning component models

From a dynamic dataset *D*<sup>*t*</sup>, the first learning step of the proposed algorithm is to compute a set of component models, i.e. simple models which will be combined in the second part of the algorithm.

Considering a number *m* of simple models to learn, we first compute *m* local perturbations of *D*<sup>*t*</sup>, denoted by {*D*<sup>*t*[*k*]</sup>}<sub>*1*≤*k*≤*m*</sub>, by sampling from *D*<sup>*t*</sup> with replacement (bagging strategy [21]). Then, for each *D*<sup>*t*[*k*]</sup>, a directed graph *G*<sup>*k*</sup> = (*V*, *E*<sup>*k*</sup>), with *V* having one vertex for each of the *N* original variables, is built by first randomly choosing α·*N*·(*N* − 1)/2 undirected edges and then computing the two directed scores *s*(*A* → *B*) and *s*(*B* → *A*) for each of them. Finally, each graph *G*<sup>*k*</sup> is searched for its optimum spanning arborescence *A*<sup>*k*</sup> with respect to the score *s*, using the Edmonds algorithm [6].

Even if the built graphs only have one node per original variable (as opposed to two nodes, one for timestep *t* and one for timestep *t* + 1), the semantics of an arc *X* → *Y* measures the influence of *X* at time *t* over *Y* at time *t* + 1. This semantics is taken into account during scores' computation. Indeed, a score computation *s*(*X* → *Y*) is actually a score involving *D*<sub>*x*</sub><sup>*t*[*k*]</sup> and *D*<sub>(*y*+*N*)</sub><sup>*t*[*k*]</sup> columns of *D*<sup>*t*</sup>, where *x* and *y* are the indices of variables *X* and *Y* in *D*. Concerning the choice of score itself, many asymmetrical scores can be used. A simple one with conditionals interpretation is the conditional entropy *H*(*X*|*Y*). Bayesian scores, like Bayesian Dirichlet variants [7] can also be used with the advantage of being able to add prior information, relatively to each edge, to the learning task.

The edge sampling step before spanning arborescence computation is important in order to counterbalance the determinism of the Edmonds algorithm due to its global optimality. The choice of α at this step is critical and will be discussed in Section 4. On one hand, it must be low enough to avoid that the optimality of the spanning algorithm restricts the component models diversity too much before the combination step. On the other hand, it should also be high enough so that the spanning algorithm still discriminates between different models and the resulting components are not just the consequences of random sampling. The choice of sampling undirected edges instead of directly sampling directed ones is of major importance to this purpose, since it allows to obtain strongly connected *G*<sup>*k*</sup> graphs, to guarantee the existence of a spanning arborescence for a wide range of α.

To conclude, the first learning step ensures each component to be sparse, and globally optimal in the considered arborescences space, with respect to *s*. Edmonds algorithm, like Kruskal one in the undirected graphs space, allows to get this optimality for each directed component in polynomial time.

## 3.3 Computing the composite model

Once the *m* component models have been learned, the second learning step aims at combining them into a composite model. In the Bayesian Model Averaging framework, this step is achieved by computing a set of *expected features* {*E*(*f*<sub>*i*</sub>)*}$_{*i*}$ for the composite model, each *E*(*f*<sub>*i*</sub>) being inferred from each component features set {*f*<sub>*i*</sub><sup>*k*</sup>}<sub>*1*≤*k*≤*m*</sub>. In this paper, the feature space consists in the set of all possible edges in *V*<sup>2</sup>, and an expected edge score is computed by counting how often that edge was present in the arborescence *A*<sup>*k*</sup>, considering it was present in the initial weighted graph *G*<sup>*k*</sup>. Formally, we have for all (*A*, *B*) ∈ *V*<sup>2</sup>:

$$\mathbf{E}(f_{A \to B}) \approx \frac{|\{k | (A, B) \in \text{edges}(A^k)\}|}{\alpha}.$$

More complex features could be considered, such as paths instead of edges or ancestor/descendant relationships, as in [4] (although the authors do not combine them in a single model). We leave these problems for future work since it would require more complex combination rules, requiring transitive reduction techniques [25], a difficult problem in the case where the input graph has cycles or is weighted.

The computation of all edges' expected scores in the composite model directly provides a ranking for those edges. A combined model is finally built from such ranking by choosing the *k*-best edges or all edges whose
score exceeds a given threshold. The ranking itself can be used to compare the learning results with an optimal model through an AUROC evaluation.

The Algorithm 1 summarizes the proposed approach.

<caption>**Algorithm 1 The learning algorithm**</caption>

**Require:**
* $D^t$: a dynamic 2 slices of time learning dataset,
* $m$: a number of component models to learn;
* $s$: an edge directed weighting score;
* $\alpha$: a density for graphs setup pre-spanning arborescence;
* $\sigma$: an edge weight threshold for final edges keep decision;
**Ensure:** a structural model from $t$ to $t+1$ of the system
**for** $1 \le k \le m$ **do**
    * $D^{t[k]} := \text{sample_with_replacement_from}(D^t)$
    * $G^k := \text{build_strongly_connected_graph}(D^{t[k]}, s, \alpha)$
    * $A^k := \text{edmonds_spanning_arborescence}(G^k)$
    * $F^k := \text{edges}(A^k)$
**end for**
* $\forall (A, B) \in V^2: \mathbf{E}(f_{A \to B}) := |\{k | (A, B) \in F^k\}| / \alpha$
* $\text{return } G = \text{choose_top_edges}(\{\mathbf{E}(f_{A \to B})\}, \sigma)$

## 3.4 Complexity

Following the decomposition of the Algorithm 1, the time complexity can be expressed as the sum of two terms: one for the components computation, and another for the combination step. The components computation complexity is $m \cdot (s + g + e)$, where $s$ (resp. $g$, $e$) is the complexity of sampling (resp. connected graph construction and Edmonds algorithm). The complexity of the sampling step is negligible here, but the construction of the graph $G^k$ is in $O(\alpha N(N-1)) \approx O(N^2)$, as is the Edmonds algorithm computation with the Tarjan optimization for dense graphs [26] ($O(N^2 \log N)$ for sparse ones). Thus the components computation part is in $O(mN^2)$.

The combination part is trickier since it consists in a succession of joins between the component edgelsists for further counting. Depending on the join algorithm used, this part can become the bottleneck of the overall learning approach. Indeed, a simple nested join has a time complexity in $O(PQ)$ where $P$ and $Q$ are the number of rows in each table to join together. In our algorithm, this leads to a complexity in $O(N^{2m})$. However, it is possible to considerably improve this step using better strategies, such as hash joins, running in $O(P+Q)$, thus leading to a linear complexity in our settings. Note that in a purely sequential algorithm, it is not really necessary to compute joins since component edges can be counted just after arborescence computation. However, since this method is highly parallelizable due to the independence of components learning and of combinations order, it is preferable to consider this solution since the parallelization gain overcomes the joins cost in practice.

Overall, the proposed approach is the sum of a quadratic and a linear step (in a parallel configuration), and thus is of quadratic complexity.

# 4 DREAM 8 (HPN-DREAM) SC1B Network Inference Challenge Results

In this section, we validate our method and compare it with other algorithms through a dataset of the recent HPN-DREAM 8 Breast Cancer closed challenge [8].

## 4.1 Challenge and evaluation method description

The DREAM 8 SC1B subchallenge learning objective is to find the network of a synthetic biological model built using state of the art methods and biological knowledge. Simulation of this model led to the production of several time series involving 20 biological features. The data used to perform the validation of our algorithm was firstly pre-processed as described in section 3 to produce a single two-time slice dataset, the resulting data containing 80 $t$ to $t+1$ examples over 40 temporal features.

The evaluation of learning results for this task is achieved by an official tool, the DREAMTools python package [27], through the computation of an AUROC score against the golden standard. In addition to com-
<figure><img src="image_16.png" /><figcaption>**Fig. 1.** (Top) mean and sds of AUC computed by DREAMTools over 50 computations as a function of the number of combined models, as well as the samples and edge ratio used for components learning. (Bottom) mean and sds of AUC computed by DREAMTools over 50 computations as a function of the edge and sample ratios. Only convergence values for increasing *m* are plotted.</figcaption></figure>

puting scores the same way from one algorithm to another, this package also provides the expected ranking an algorithm would have reached if the challenge were still open, using all final results from the more than 100 official submissions, which allows for a cheap comparison with many algorithms of all families described in section 2.

In order to quantify the impact of several parameters on our algorithm learning quality, we have tested the method with different parametrization of the number of combined modes *m*, the ratio of samples *n* contained in each data perturbation, and the ratio of edges *α* present in each graph before each component learning. We used BDeu [7] gain as edge score, the difference between the BDeu score of the *A* → *B* local structure and the no edge one. Namely for an edge *A* → *B*: $BDeu Intr parents(B) = \{A\}) - BDeu Intr parents(B) = \emptyset).$

## 4.2 Results

Results for many parametrizations, given in Figure 1 (Top), show different clear trends. Firstly, we can see that for small edge ratios, the obtained AUC seems to monotonically increase with the number of combined models, until reaching plateaus. For bigger ratios, the trend seems to be mostly observable, but the higher the
sampling ratio, the lower the minimum edge ratio needs to be to show this trend. Additionally, we can observe
that the convergence AUC value tends to increase whenever any of the edge or the sample ratio decrease, which
is clearer in Figure 1(Bottom). These results seem to indicate that focusing on smaller parts of the available
information for each component, while aggregating a higher number of them for final consensus, seem to give
the best results, which confirms the requirement for components diversity in order to give a good consensus.
Extra experiments done and not displayed here show that smaller edge and sample ratios break the observed
trends. For edge ratio, this is predictable, since the minimum value displayed of 0.05 corresponds to an average
number of 2 neighbors per node (considering we add the reverse edges for each sampled edge, to ensure we
can compute a spanning arborescence), which is the minimal number of neighbors required for the algorithm
to make a choice. Lower values actually lead the spanning arborescence algorithm to just select most available
edges in the graphs it is given. For samples ratio, it seems to indicate a limit from which the dataset is too small
to capture faithful enough information.

Concerning the expected ranking for the different results, our approach is very promising since it reached the 3rd position for the best mean AUC obtained over the different parametrizations, outperforming GENIE3, ARACNE, all heuristic oriented Bayesian network methods, as well as all linear and most non-linear regression methods, all ODE and all ensemble learning solutions.

## 4.3 Discussion

At the moment, the gap with the best performance is of 0.045. A particularity of the considered DREAM subchallenge is that 3 out of the 20 biological features are actually fake nodes, supposed to have no correlation with the others. This shows a limitation of our approach in its current form: learning spanning trees means that every node will get one parent per component, even if there is no true correlation. Note that this problem is not necessarily easy to solve, since there is also a tendency for such spurious correlations to be non-uniformly distributed. Indeed, the optimization of the spanning arborescence score encourages to keep the apparently more correlated pairs of nodes, so the ones with the most biased noise are chosen. Since sample and edge samplings are uniformly done, there is a high probability for a restricted number of parents to appear in each component for a fake node. In practice, this means that a simple pruning of the components is not enough. Future work will address this issue.

# 5 Conclusion

In this paper, we have presented a network inference learning algorithm based on the combination of mul-
tiple spanning arborescences learned over multiple perturbation of the original dataset, with enforced diversity
through edge sampling, showing promising results in practice on a recent DREAM challenge.

Experiments have more particularly shown that combining more models together with more diversity, in-
volving a decrease of both sample and edge ratios in the currently defined parameter space, leads to better
convergence values. It is encouraged to use this strategy in a quite extreme way, since best performances ob-
tained in the experiments are achieved by situations where both ratios are very low. The only warning would be
to still allow the Edmonds algorithm to have choice, in order not to make the components completely random.
We have also seen in section 4 the impact of fake nodes on the results, and the difficulty of identifying them
whenever the spanning arborescence assign most nodes a parent. This issue has a significant impact on the
current results since removing edges involving fake nodes would lead to a top position of the approach.

Future works will address the limitations of the current algorithm, such as its sensitivity to fake nodes.
More advanced extensions will also be investigated, such as the introduction of priors, really important in
biological contexts, modular capabilities, which is becoming a standard in recent methods to abstract from a
model complexity, and different component combination rules to preserve extra properties in the consensus
model, such as paths or path lengths.

## References

[1] Michael Hecker, Sandro Lambeck, Susanne Toepfer, Eugene van Someren, and Reinhard Guthke. Gene regulatory network inference: Data integration in dynamic models—a review. *Biosystems*, 96(1):86–103, 2009.

[2] Karen Sachs, Omar Perez, Dana Pe'er, Douglas A. Lauffenburger, and Garry P. Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. *Science*, 308(5721):523–529, 2005.
[3] Ludwig Geistlinger, Gergely Csaba, Simon Dirmeier, Robert Küffner, and Ralf Zimmer. A comprehensive gene regulatory network for the diauxic shift in saccharomyces cerevisiae. *Nucleic acids research*, page gkt631, 2013.

[4] Nir Friedman, Moises Goldszmidt, and Abraham Wyner. Data analysis with bayesian networks: A bootstrap approach. In *Proceedings of the Fifteenth conference on Uncertainty in artificial intelligence*, pages 196–205. Morgan Kaufmann Publishers Inc., 1999.

[5] Bradley M Broom, Kim-Anh Do, and Devika Subramanian. Model averaging strategies for structure learning in bayesian networks with limited data. *BMC bioinformatics*, 13(13):S10, 2012.

[6] Jack Edmonds. Optimum branchings. *Mathematics and the Decision Sciences, Part*, 1:335–345, 1968.

[7] David Heckerman, Dan Geiger, and David M Chickering. Learning bayesian networks: The combination of knowledge and statistical data. *Machine learning*, 20(3):197–243, 1995.

[8] Steven M Hill, Laura M Heiser, Thomas Cokelaer, Michael Unger, Nicole K Nesser, Daniel E Carlin, Yang Zhang, Artem Sokolov, Evan O Paull, Chris K Wong, et al. Inferring causal molecular networks: empirical assessment through a community-based effort. *Nature methods*, 13(4):310–318, 2016.

[9] I Chebil, Rémy Nicolle, G Santini, Céline Rouveirol, and Mohamed Elati. Hybrid method inference for the construction of cooperative regulatory network in human. *IEEE transactions on nanobiology*, 13(2):97–103, 2014.

[10] Rémy Nicolle, François Radvanyi, and Mohamed Elati. Coregnet: reconstruction and integrated analysis of co-regulatory networks. *Bioinformatics*, page btv305, 2015.

[11] Va Anh Huynh-Thu, Alexandre Irrthum, Louis Wehenkel, and Pierre Geurts. Inferring regulatory networks from expression data using tree-based methods. *PLOS ONE*, 5(9):1–10, 09 2010.

[12] Miron B Kursa and Witold R Rudnicki. Feature selection with the boruta package. *Journal of Statistical Software*, 36(11):1–13, 2010.

[13] Anne-Claire Haury, Fantine Mordelet, Paola Vera-Licona, and Jean-Philippe Vert. Tigress: Trustful inference of gene regulation using stability selection. *BMC Systems Biology*, 6(1):145, 2012.

[14] Alex Greenfield, Christoph Hafemeister, and Richard Bonneau. Robust data-driven incorporation of prior knowledge into the inference of dynamic regulatory networks. *Bioinformatics*, 29(8):1060–1067, 2013.

[15] Adam A Margolin, Ilya Nemenman, Katia Basso, Chris Wiggins, Gustavo Stolovitzky, Riccardo D Favera, and Andrea Califano. Aracne: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context. *BMC bioinformatics*, 7(Suppl 1):S7, 2006.

[16] Min Zou and Suzanne D Conzen. A new dynamic bayesian network (dbn) approach for identifying gene regulatory networks from time course microarray data. *Bioinformatics*, 21(1):71–79, 2005.

[17] Judea Pearl. Heuristics: intelligent search strategies for computer problem solving. 1984.

[18] Peter Norvig. *Paradigms of artificial intelligence programming: case studies in Common LISP*. Morgan Kaufmann, 1992.

[19] Jan Karel Lenstra. *Local search in combinatorial optimization*. Princeton University Press, 2003.

[20] Walter R Gilks. Markov chain monte carlo. *Encyclopedia of Biostatistics*, 2005.

[21] Leo Breiman. Bagging predictors. *Machine learning*, 24(2):123–140, 1996.

[22] François Schnitzler, Sourour Ammar, Philippe Leray, Pierre Geurts, and Louis Wehenkel. *Efficiently Approximating Markov Tree Bagging for High-Dimensional Density Estimation*, pages 113–128. Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.

[23] C Chow and C Liu. Approximating discrete probability distributions with dependence trees. *IEEE transactions on Information Theory*, 14(3):462–467, 1968.

[24] Sourour Ammar and Philippe Leray. *Mixture of Markov Trees for Bayesian Network Structure Learning with Small Datasets in High Dimensional Space*, pages 229–238. Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.

[25] Andrea Pinna, Sandra Heise, Robert J Flassig, Alberto De La Fuente, and Steffen Klamt. Reconstruction of large-scale regulatory networks based on perturbation graphs and transitive reduction: improved methods and their evaluation. *BMC systems biology*, 7(1):73, 2013.

[26] Robert Endre Tarjan. Finding optimum branchings. *Networks*, 7(1):25–35, 1977.

[27] Thomas Cokelaer, Mukesh Bansal, Christopher Bare, Erhan Bilal, Brian M Bot, Elias Chaibub Neto, Federica Eduati, Alberto de la Fuente, Mehmet Gönen, Steven M Hill, et al. Dreamtools: a python package for scoring collaborative challenges. *F1000Research*, 4, 2015.
Long-term Tracking of Budding Yeast Cells with CellStar

Cristian VERSARI¹, Szymon STOMA², Kirill BATMANOV¹, Artémis LLAMOSI³,⁴, Filip MROZ⁵, Adam

KACZMAREK⁵, Matt DEYELL³, Cédric LHOUSSAINE¹, Pascal HERSEN³ and Gregory BATT⁴

¹ CRISTAL - Univ. Lille 1, bât M3 ext. Av. Carl Gauss, 59655, Villeneuve d'Ascq, France

² ScopeM - ETH Zurich, Wolfgang-Pauli-Str. 14, 8093, Zurich, Switzerland

³ Lab. MSC, bât. Condorcet 10 rue A. Domon et L. Duquet, 75205, Paris cedex 13, France

⁴ Inria Saclay, Bât A. Turing 1 rue H. D'Estienne d'Orves, 91120, Palaiseau, France

⁵ Inst. of Computer Science - Univ. of Wroclaw, Joliot-Curie 15, 50-383, Wroclaw, Poland

Corresponding Author: cristian.versari@univ-lille1.fr

Abstract

Observing cellular processes at the single cell level is often necessary to understand how cells respond to endogenous and environmental changes. Used in combination with fluorescence reporter techniques, flow cytometry and time-lapse microscopy are arguably the two most widely employed quantitative single-cell observation approaches. The former provides great statistical details on the diversity of the studied cell population, whereas the latter provides longitudinal information on single cells: individual cells can be tracked in time. This is a decisive advantage to investigate a number of important biological problems, including chronological aging, epigenetic heritability, and dynamic features such as cell-cycle and circadian oscillations in non-synchronized cell populations.

However, the capability to extract single cell traces from microscopy images in a fully-automated manner is a necessary prerequisite to obtain conclusions that are valid and biologically relevant in long lasting experiments. Incorrect assignments (e.g. two cells exchanged at some time point) can possibly hide interesting features, or worse, create spurious information. Although such incorrect assignments are expected to be relatively rare at each time point, a simple analysis shows that the number of correct traces decreases rapidly with the duration of the experiment: even for yeast cells that have relatively regular shapes, no solution has been proposed that reaches the high quality required for long-term experiments for segmentation and tracking based on brightfield images.

In this demo we present CellStar [1], a tool chain designed to achieve good performance in long-term experiments. The key features are the use of a new variant of parametrized active rays for segmentation, a neighbourhood-preserving criterion for tracking, and the use of an iterative approach that incrementally improves segmentation and tracking quality. A graphical user interface enables manual corrections of segmentation and tracking errors and their use for the automated correction of other, related errors and for parameter learning.

In [1] we created a benchmark dataset with manually analysed images and compared CellStar with six other tools, showing its high performance, notably in long-term tracking. As a community effort, we set up a website, the Yeast Image Toolkit, with the benchmark and the Evaluation Platform to gather this and additional information provided by others.

References

[1] Cristian Versari, Szymon Stoma, Kirill Batmanov, Artémis Llamosi, Filip Mroz, Adam Kaczmarek, Matt Deyell, Cédric Lhoussaine, Pascal Hersen and Gregory Batt. Long-term tracking of budding yeast cells in brightfield microscopy: CellStar and the Evaluation Platform. *Journal of The Royal Society Interface*, 14(127), 2016.
JOBIM 2017 Highlight

Listeriomics: A Multi-Omics Interactive Web Platform for Systems Biology of the
Model Pathogen *Listeria*

Christophe Bécavin¹,²,³,⁴, Mikael Koutero¹,²,³, Nicolas Tchitchek⁵, Franck Cerutti⁶, Pierre Lechat⁴, Nicolas Maillet⁴, Claire Hoede⁶, Hélène Chiapello⁶, Christine Gaspin⁶, and Pascale Cossart¹,²,³

¹ Institut Pasteur, Unité des Interactions Bactéries-Cellules, Département de Biologie Cellulaire et Infection, F-75015 Paris, France

² INSERM, U604, F-75015 Paris, France

³ INRA, USC2020, F-75015 Paris, France

⁴ Institut Pasteur - Bioinformatics and Biostatistics Hub - C3BI, USR 3756 IP CNRS - Paris, France

⁵ CEA, Division of Immuno-Virology, Institute of Emerging Diseases and Innovative Therapies, IDMIT Center, Fontenay-aux-Roses 92265, France

⁶ MIAT, Université de Toulouse, INRA, Castanet-Tolosan, France

Corresponding Author: christophe.becavin@pasteur.fr

**Paper Reference:** Bécavin et al. (2016) Listeriomics: An Interactive Web Platform for Systems Biology of Listeria. mSystems. http://dx.doi.org/10.1128/mSystems.00186-16

**Abstract:** As for many model organisms, the amount of Listeria omics data produced has recently increased exponentially. There are now >80 published complete Listeria genomes, around 350 different transcriptomic data sets, and 25 proteomic data sets available. The analysis of these data sets through a systems biology approach and the generation of tools for biologists to browse these various data are a challenge for bioinformaticians. We have developed a web-based platform, named Listeriomics, that integrates different tools for omics data analyses, i.e., (i) an interactive genome viewer to display gene expression arrays, tiling arrays, and sequencing data sets along with proteomics and genomics data sets; (ii) an expression and protein atlas that connects every gene, small RNA, antisense RNA, or protein with the most relevant omics data; (iii) a specific tool for exploring protein conservation through the Listeria phylogenomic tree; and (iv) a coexpression network tool for the discovery of potential new regulations. Our platform integrates all the complete Listeria species genomes, transcriptomes, and proteomes published to date. This website allows navigation among all these data sets with enriched metadata in a user-friendly format and can be used as a central database for systems biology analysis.

**Link:** http://listeriomics.pasteur.fr/

**Keywords** Listeria, genomics, transcriptomics, proteomics, systems biology

**Summary**

Listeria monocytogenes is a foodborne pathogen responsible for foodborne infections with a mortality rate of 25%. This pathogen is responsible for gastroenteritis, sepsis, and meningitis and can cross three host barriers, the intestinal, placental, and blood-brain barriers [1]. Over the past three decades *Listeria* has become a model organism for host-pathogen interactions, leading to critical discoveries in a broad range of fields [2], including virulence-factor regulation, cell biology, and bacterial pathophysiology. To study these mechanisms, several genomics, transcriptomics, and proteomics data sets have been produced.

The analysis of these data sets through a systems biology approach and the generation of tools for biologists to browse these various data are a challenge for bioinformaticians. We have developed a web-based platform, named Listeriomics (http://listeriomics.pasteur.fr/), that integrates different tools (see Figure 1) for omics data analyses, i.e., (i) an interactive genome viewer to display gene expression arrays, tiling arrays, and
sequencing data sets along with proteomics and genomics data sets; (ii) an *expression and protein atlas* that connects every gene, small RNA, antisense RNA, or protein with the most relevant omics data; (iii) a specific tool for exploring *protein conservation through the Listeria phylogenomic tree*; and (iv) a *coexpression network tool* for the discovery of potential new regulations.

To our knowledge, none of the referent databases dedicated to model organisms such as *E. coli* or *B. subtilis* integrates as many data sets and visualization tools as the Listeriomics resource does. User experience and feedback from our collaborators using the Listeriomics interface [3] for the past 5 years were driving forces in organizing and improving the way to access data and tools. Our main purpose was to design an *easy-to-use website with a dynamic interface for biologists* wanting to access the different heterogeneous “omics” data sets available for Listeria. As such our website *should interest the JOBIM community as it shows an example of extensive multi-omics data integration for model organism*.

The website is developed with an in-house platform which we named BACNET. It is written in Java language, with BioJava and Eclipse RAP APIs. We designed the platform so bioinformaticians wishing to *create their own multi-omics website for another organism can do so with few efforts*. We believed the description of the BACNET platform will also interest the JOBIM community.

<figure><img src="image_17.png" /><figcaption>Figure 1: Overview of the different tools available in Listeriomics website</figcaption></figure>

## References

[1] P. Cossart, “Illuminating the landscape of host-pathogen interactions with the bacterium Listeria monocytogenes.,” *Proc. Natl. Acad. Sci. U. S. A.*, Nov. 2011.

[2] A. Lebreton, F. Stavru, S. Brisse, and P. Cossart, “1926–2016: 90 Years of listeriology,” *Microbes Infect.*, vol. 18, no. 12, pp. 711–723, Dec. 2016.

[3] F. Impens, N. Rolhion, L. Radoshevich, C. Bécavin, M. Duval, J. Mellin, F. García del Portillo, M. G. Pucciarelli, A. H. Williams, and P. Cossart, “N-terminomics identifies Prli42 as a membrane miniprotein conserved in Firmicutes and critical for stressosome activation in Listeria monocytogenes,” *Nat. Microbiol.*, vol. 2, no. February, p. 17005, 2017.
# AskOmics, a web tool to integrate and query biological data using semantic web technologies

Xavier GARNIER<fnref n="1" />,<fnref n="2" />, Anthony BRETAUDEAU<fnref n="1" />,<fnref n="3" />, Olivier FILANGI<fnref n="1" />,<fnref n="3" />, Fabrice LEGEAI<fnref n="1" />,<fnref n="4" />, Anne SIEGEL<fnref n="2" /> and Olivier DAMERON<fnref n="2" />

Corresponding author: xavier.garnier@irisa.fr

Research programs involving genetics, genomics and epigenetics are quickly growing. Lot of experiments producing large amount of data are now feasible in the frame of a laboratory. As well, the tools analyzing the data generated by these experiments are often available. Results of these experiments can be stored and structured into files and loaded into databases, and lots of these databases are public. However, these bases are usually implemented according to schemes and techniques that do not allow their interoperability in an easy manner.

The technologies from the Semantic Web, especially RDF and SPARQL are one of the key elements for combining databases, which has led to the emergence of linked data. It is based on triples (subject, predicate and objects) describing the relationships between elements stored into interoperable triplestores allowing distributed querying. Because of its flexibility, versatility and ontology-awareness, numerous biological databases, such as UniProt, PubChem or ChEBI and Reactome at EBI, give access to their data via a SPARQL endpoint.

We present AskOmics, a new software, that uses the Semantic web technologies, which helps to integrate multiple format of data and query them through a user-friendly interface. AskOmics is a free and open-source software (AGPL licence) available on GitHub (https://github.com/askomics/askomics).

AskOmics supports both intuitive data integration and querying while shielding a non-expert user from most of the technical difficulties underlying the web semantic technologies. Because large and heterogeneous biological datasets are often difficult to integrate, AskOmics users can provide simple tabulation-separated files (TSV), that are transformed automatically into RDF triples, then stored into a triplestore. Finally, for data querying, AskOmics provides a visually intuitive interface to obtain a comprehensive view of the biological study.

During data integration, user provides input files in common formats (currently TSV and GFF) to be converted into RDF triples. AskOmics generates triples corresponding to the data (the content), and also triples which describe the data (the abstraction). Triples are loaded in a triplestore in order to persist data and optimize queries.

The query interface is composed of a dynamic graph at the left and a right view for filtering attributes. On the graph, each node represents an entity. Entities are linked between them with arrows. Attributes of the selected entities are displayed on the right view. To build the graph, AskOmics query the abstraction. Users build their queries by starting from a node of interest and sequentially select its neighbors and filter on attributes, creating a path on the abstraction. This path is converted into a SPARQL query and sent to the triplestore. Finally, the results are displayed as a table and can be downloaded as a TSV file.

AskOmics has been applied successfully to the analysis of large scale datasets on the aphid embryogenesis, on the variability of Brassicaceae in response to clubroot disease, and to the analysis of biological pathways.

<footnote n="1">¹ Institut de Génétique, Environnement et Protection des Plantes (IGEPP) - Institut national de la recherche agronomique (INRA) : UMR1349, Agrocampus Ouest - Agrocampus Ouest, UMR1349 IGEPP, F-35 042 Rennes, France</footnote>

<footnote n="2">² DYLISS (INRIA - IRISA) - INRIA, Université de Rennes 1, CNRS : UMR6074 - Campus de Beaulieu, F-35 042 Rennes Cedex, France</footnote>

<footnote n="3">³ Plateforme bioinformatique GenOuest - Université de Rennes 1, Biogenouest - France</footnote>

<footnote n="4">⁴ GENSCALE (INRIA - IRISA) - Université de Rennes 1, CNRS : UMR6074, INRIA - Campus de Beaulieu, F-35 042 Rennes Cedex, France</footnote>
# Regulatory and signaling network assembly through Linked Open Data

Marie LEFEBVRE<fnref n="1" />, Jérémie BOURDON<fnref n="2" />, Carito GUZIOLOWSKI<fnref n="2" /> and Alban GAIGNARD<fnref n="1" />

Corresponding author: marie.lefebvre@univ-nantes.fr, alban.gaignard@univ-nantes.fr

## 1 Introduction and problem statement

Nowadays, huge efforts address the organization of biological knowledge through linked open databases. These databases can be automatically queried to reconstruct a large variety of biological networks such as regulatory or signaling networks. Assembling networks still implies manual operations due to (i) source-specific identification of biological entities, (ii) source-specific semantics for entity-entity relationships, (iii) proliferating heterogeneous life-science databases with redundant information and (iv) the difficulty of recovering the logical flow of a biological pathway due to the bidirectionality of chemical reactions. Homogenization of biological networks is therefore costly and error-prone. Existing tools such as the ReactomeFIPlugIn of Cytoscape 3.0[1] or STRING[2] allow to link entities to each other or to identify an entity's membership to a single pathway. Nevertheless, they are still limited in the global modeling aspects (logical rules inferred from the knowledge representation). Here, we present a framework to automate the assembly of regulatory and signaling networks in the context of tumor cells modeling.

## 2 Approach

Our framework is based on Semantic Web technologies. It addresses (i) the uniform identification of multi-source biological entities, (ii) the description of labeled directed graphs through RDF, and (iii) the use of BioPAX [3] as a semantic reference. We consider a list of target gene names or IDs as entry points. The first step consists in retrieving transcription factors (TFs) controlling these target genes. Then, the second step consists in considering the TFs as new entry points for the reconstruction algorithm. The full regulatory network is finally assembled by iteratively applying the second step until no new TFs can be found.

## 3 Demonstration

To assemble networks our algorithm queries PathwayCommons[4] through its SPARQL endpoint and retrieves a graph of TFs associated to target genes. We developed a web tool that displays the biological network assembly step by step, allowing users to interact with the reconstruction process and to visually shape the network. Through this web tool, it is also possible to launch a command line tool (Java) to address larger scale input gene lists. These tools have been deployed on the BiRD Cloud infrastructure. From a list of 1800 targets genes, we were able to assemble in less than 3 minutes a graph of 1474 nodes and 12303 edges.

## 4 Discussion and Conclusion

As future works, we aim at integrating drug-target informations (e.g. KEGGdrug, DrugBank) through SPARQL federated queries to get insights on (i) tumor cells growth and (ii) drug response on patient and cell lines gene expression data. Our tool is freely available at https://github.com/symmetric-group/bionets-demo

### Acknowledgements

This work was supported by the BiRD bioinformatics facility, the SymETRIC project and the GRIOTE project.

### References

<footnote n="1">¹ Nantes Academic Hospital, CHU de Nantes, France</footnote>

<footnote n="2">² LS2N - UMR 6004, University of Nantes, Ecole Centrale de Nantes, France</footnote>

<footnote>[1] W., Guanming et al. A human functional protein interaction network and its application to cancer data analysis. *Genome Biology*, 11(5):R53, 2010.</footnote>

<footnote>[2] D., Szklarczyk et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. *Nucleic Acids Res.*, 43:D447–D452, 2015.</footnote>

<footnote>[3] E., Demir et al. The biopax community standard for pathway data sharing. *Nature Biotechnology*, 28:935–942, 2010.</footnote>

<footnote>[4] E.G., Cerami et al. Pathway commons, a web resource for biological pathway data. *Nucleic Acids Res.*, 39:D685–90, 2011.</footnote>
Scientific Workflows for Computational Reproducibility in the Life Sciences:
Status, Challenges and Opportunities

Sarah COHEN-BOULAKIA<sup>1,2,3</sup>, Khalid BELHAJJAME<sup>4</sup>, Olivier COLLIN<sup>5</sup>, Jérôme CHOPARD<sup>6</sup>, Christine

FROIDEVAUX<sup>1</sup>, Alban GAIGNARD<sup>7</sup>, Konrad HINSEN<sup>8</sup>, Pierre LARMANDE<sup>9</sup>, Yvan LE BRAS<sup>10</sup>, Frédéric LEMOINE<sup>11</sup>, Fabien MAREUIL<sup>12,13</sup>, Hervé MÉNAGER<sup>12,13</sup>, Christophe PRADAL<sup>14,2</sup> and Christophe

BLANCHET<sup>15</sup>

<sup>1</sup> LRI, Univ. Paris-Sud, CNRS UMR 8623, Univ. Paris-Saclay

<sup>2</sup> Inria, VirtualPlants, Montpellier

<sup>3</sup> Inria, Zenith, Montpellier

<sup>4</sup> Lamsade, Univ. Paris-Dauphine, CNRS UMR 7243, Paris

<sup>5</sup> IRISA, Rennes, France

<sup>6</sup> INRA, UMR729, MISTEA, Montpellier

<sup>7</sup> Nantes Academic Hospital, CHU de Nantes

<sup>8</sup> Centre de Biophysique Moléculaire, CNRS UPR4301, Orléans

<sup>9</sup> IRD, DIADE, F-34394 Montpellier

<sup>10</sup> EnginesOn / INRIA, Rennes

<sup>11</sup> Institut Pasteur, Unité Bioinformatique Evolutive, C3BI, USR 3756 IP CNRS, Paris

<sup>12</sup> Institut Pasteur, Hub Bioinformatique et Biostatistique, C3BI, USR 3756 IP CNRS, Paris

<sup>13</sup> Institut Pasteur, Centre d'Informatique pour la Biologie, DSI, Paris

<sup>14</sup> CIRAD, UMR AGAP, Montpellier

<sup>15</sup> Institut Français de Bioinformatique, IFB-core, CNRS, UMS 3601, Gif-sur-Yvette

Corresponding author: cohen@lri.fr

**Reference paper:** Cohen-Boulakia et al. (2017) Scientific Workflows for Computational Reproducibility in the Life Sciences: Status, Challenges and Opportunities. *Future Generation Computer Systems*.
http://dx.doi.org/10.1016/j.future.2017.01.012

**Abstract** With the development of new experimental technologies, an avalanche of data has to be computationally analyzed for scientific advancements and discoveries to emerge. Faced with the complexity of analysis pipelines, the large number of computational tools, and the enormous amount of data to manage, there is compelling evidence that many if not most scientific discoveries will not stand the test of time: increasing the reproducibility of computed results is of paramount importance. The objective we set out in this paper is to place scientific workflows in the context of reproducibility: We define several levels of reproducibility; we characterize and define the criteria that need to be catered for by reproducibility-friendly scientific workflow systems; we use such criteria to place several representative and widely used workflow systems and companion tools within such a framework; we discuss the remaining challenges posed by reproducible scientific workflows in the life sciences. Our study was guided by three use cases from the French community, involving in silico experiments.

**Keywords** Reproducibility; Scientific Workflows; Provenance; Packaging environments.

Novel technologies have led to the generation of very large volumes of data at an unprecedented rate. This is particularly true for the life sciences, where, for instance, innovations in Next Generation Sequencing (NGS) have led to a revolution in genome sequencing. Current instruments can sequence several hundreds of human genomes in one week whereas more than ten years have been necessary for the first human genome. Many laboratories have thus acquired NGS machines, resulting in an avalanche of data which has to be further analyzed using a series of tools and programs for new scientific knowledge and discoveries to emanate. The same kind of situation occurs in completely different biological domains, such as plant phenotyping which aims at understanding the complexity of interactions between plants and environments in order to accelerate the discovery of new genes and traits thus optimize the use of genetic diversity under different environments. Here, thousands of plants are grown in controlled environments, capturing a lot of information and generating huge amounts of raw data to be stored and then analyzed by very complex computational analysis pipelines for scientific advancements and discoveries to emerge.

Faced with the complexity of analysis pipelines designed, the number of computational tools available and the amount of data to manage, there is compelling evidence that the large majority of scientific discoveries will not stand the test of time: increasing reproducibility of results is of paramount importance.
Over the recent years, many authors have drawn attention to the rise of purely computational experiments which are not reproducible (see in particular [1]). Major reproducibility issues have been highlighted in a very large number of cases: while [2] has shown that even when very specific tools were used, textual description of the methodology followed was not sufficient to repeat experiments, [3] has focused on top impact factor papers and shown that insufficient data were made available by the authors to make experiments reproducible, despite the data publication policies recently put in place by publishers.

Scientific communities in different domains have started to act in an attempt to address this problem. Pres-
tigious conferences (such as two major conferences from the database community, namely, VLDB<sup>15</sup> and SIG-
MOD<sup>16</sup>) and journals such as PNAS<sup>17</sup>, Biostatistics, Nature, and Science, to name only a few, encourage or
require published results to be accompanied by all the information necessary to reproduce them. However,
making their results reproducible remains a very difficult and extremely time-consuming task for most authors.
In the meantime, considerable efforts have been put into the development of scientific workflow management
systems. They aim at supporting scientists in developing, running, and monitoring chains of data analysis pro-
grams. A variety of systems (Galaxy, OpenAlea,...) have reached a level of maturity that allows them to be
used by scientists for their bioinformatics experiments, including analysis of NGS or plant phenotyping data.
By capturing the exact methodology followed by scientists (in terms of experimental steps associated with tools
used) scientific workflows play a major role in the reproducibility of experiments.

The propose of this paper is thus to better understand the core problematic of reproducibility in the specific context of scientific workflow systems. We aim to provide answers to the following key points: How can we define the different levels of reproducibility that can be achieved when a workflow is used to implement an *in silico* experiment? Which are the criteria of scientific workflow systems that make them *reproducibility-friendly*? What is concretely offered by the scientific workflow systems in use in the life science community to deal with reproducibility? Which are the open problems to be tackled in computer science (in algorithmics, systems, knowledge representation etc.) which may have huge impact on the problems of reproducing experiments when using scientific workflow systems?

Accordingly, we make the following five contributions: We present three real use cases involving in silico experiments, and elicit concrete reproducibility issues that they raise. We define several kinds of reproducibility that can be reached when scientific workflows are used to perform experiments. We characterize and define the criteria that need to be catered for by *reproducibility-friendly* scientific workflow systems. Using the framework of the criteria identified, we place several representative and widely used workflow systems and companion tools within such a framework. We go on to discuss the challenges posed by reproducible scientific workflows in the life sciences and describe the remaining opportunities of research in several areas of computer science which may address them.

This paper is the result of a large collaborative work between several members of the French Bioinformatics community. Ongoing work includes organizing hackathons to concretely test the ability of the various workflow systems to deal with reproducibility.

Acknowledgements

The authors acknowledge the support of GDR CNRS MaDICS, programme CPER Région Bretagne "CeSGO", and programme Région Pays de la Loire "Connect Talent" (SyMeTRIC). We acknowledge funding by the call "Infrastructures in Biology and Health" in the framework of the French "Investments for the Future" (ANR-11-INBS-0012 and ANR-11-INBS-0013).

References

[1] Victoria Stodden, Peixuan Guo, and Zhaokun Ma. Toward reproducible computational research: an empirical analysis of data and code policy adoption by journals. PloS one, 8(6):e67111, 2013.

[2] Anton Nekrutenko and James Taylor. Next-generation sequencing data interpretation: enhancing reproducibility and accessibility. *Nature Reviews Genetics*, 13(9):667–672, 2012.

[3] Alawi A Alsheikh-Ali, Waqas Qureshi, Mouaz H Al-Mallah, and John PA Ioannidis. Public availability of published research data in high-impact journals. *PloS one*, 6(9):e24357, 2011.

<footnote n="15">15. International conference on Very Large Data Bases</footnote>

<footnote n="16">16. ACM's Special Interest Group on Management Of Data.</footnote>

<footnote n="17">17. http://www.pnas.org/site/authors/format.xhtml</footnote>
# Sequanix: A Dynamic Graphical Interface for Snakemake Workflows

Dimitri DESVILLECHABROL¹, Rachel LEGENDRE¹,², Christiane BOUCHIER¹, Sean KENNEDY¹ and Thomas COKELAER¹,²

¹ Institut Pasteur - Biomics Pole - CITECH - Paris, France

² Institut Pasteur - Bioinformatics and Biostatistics Hub - C3BI, USR 3756 IP CNRS - Paris, France

Corresponding author: dimitri.desvillechabrol@pasteur.fr, thomas.cokelaer@pasteur.fr

Sequencing platforms produce large and heterogeneous data sets. For example, the Biomics pole (Institut Pasteur, CITECH) provides sequencing data related to transcriptomics, genomics and metagenomics. Before delivering results, platforms need to perform quality controls. Nevertheless, these platforms may also provide dedicated analysis (e.g. variant caller, de-novo).

There are many sequencing pipelines available. We provide some of them in the Sequana project. The particularity of the pipelines available in Sequana is that they are based on the Snakemake framework. Snakemake pipelines are written in Python with a rule-based syntax; a configuration file is also required (YAML format) [1]. The Snakemake pipelines targets an audience of developers since they require the pipelines to be run on cluster with command line arguments; configuration file also need to be edited.

Sequana project is made of (i) a Python library, (ii) a set of Snakemake pipelines and (iii) standalones [2]. In order to expose the Snakemake pipelines to all kind of users, we develop Sequanix, a standalone Graphical User Interface (GUI) based on Python and PyQt. Sequanix can launch Snakemake pipelines available in Sequana without the need for command line interface. Indeed, the GUI can be used on a cluster (with a display). So the knowledge of editor such as VIM is not necessary anymore. Therefore, Sequanix fills the gap between Snakemake developers and their end-users.

Pipelines provided by Sequana are directly available in a dropdown box. When a pipeline is chosen, the embedded configuration file is automatically loaded as a form. This form can be edited without the need for a text editor. It also have convenient features: if a parameter name ends in _file or _browser, a dedicated widget (file or directory browser) is shown instead of a simple line edit; if comments are available in the configuration file, we interpret and show them as tooltips in the GUI. Once a pipeline and a working directory are set, the project can be saved, the workflow visualised as a directed acyclic graph (DAG) and finally the pipeline can be executed. The interface displays the Snakemake standard output and the progress.

Sequanix exposes all Sequana pipelines (snakemake-based) within a graphical interface. Yet, many Snakemake pipelines are developed by a large community of developers especially in the NGS field. So, we extended the GUI so that any external Snakemake pipeline can also be imported and executed through the interface making Sequanix a generic tool for any Snakemake pipeline.

## References

[1] Johannes Köster and Sven Rahmann. Snakemake - a scalable bioinformatics workflow engine. *Bioinformatics*, 28(19):2520-2522, 2012.

[2] Sequana: a set of flexible genomic pipelines for processing and reporting NGS analysis (http://github.com/sequana/sequana). Documentation available on http://sequana.readthedocs.io
# Biosphère : un portail web haut niveau pour une utilisation bioinformatique des clouds

Bryan BRANCOTTE¹, Mohamed BEDRI¹, Jonathan LORENZO¹, Sandrine PERRIN¹, Frédéric SÉNÉ¹, Awa SEPOU NGAÎLO¹, Christophe BLANCHET¹, Jean-François GIBRAT¹

¹ CNRS, UMS 3601 ; Institut Français de Bioinformatique, IFB-core, Avenue de la Terrasse, F-91190 Gif-sur-Yvette, France

Corresponding Author: bryan.brancotte@france-bioinformatique.fr

Le cloud, par ses promesses de machines à grande mémoire et clusters de calcul conséquents, se pose comme un sérieux candidat pour répondre aux problématiques tant scientifiques que techniques amenées par le déluge de données en bioinformatique. Afin de s'appuyer sur le cloud, il faut rendre son utilisation simple, intuitive, et conserver un haut niveau de sécurité. C'est pour répondre à ce défi technique, et dans le cadre du projet CYCLONE, que nous avons construit le portail web Biosphère (https://biosphere.france-bioinformatique.fr).

Le projet CYCLONE (http://www.cyclone-project.eu) est un projet Horizon-2020 « Actions d'innovation » financé par la Commission Européenne (projet 644925). Il vise à répondre aux challenges posés par des applications, dont certaines en bioinformatique, en intégrant et améliorant des solutions open-source pour la gestion de clouds. Les objectifs sont une gestion unifiée de différents clouds, le déploiement et le maintien de plateformes de traitement de données complexes dans une architecture multi-cloud, une réactivité et une élasticité dans l'utilisation des ressources proposées. Dans le cadre de CYCLONE, l'Institut Français de Bioinformatique (IFB) est en charge de la définition des besoins et cas d'utilisation en bioinformatique, la proposition de machines virtuelles prédéfinies (appliances) y répondant, la formation des utilisateurs à l'utilisation de ces appliances, la création de nouvelles appliances, et la proposition d'un portail à destination des utilisateurs.

Biosphère est un portail web à haut niveau d'abstraction qui permet à l'utilisateur de lancer les appliances "en un clic" sur différents clouds sans connaissance technique particulière, ou de configurer leur déploiement. Ce portail lui permet aussi de surveiller et gérer les déploiements, et la consommation des ressources. Le portail Biosphère contient un catalogue d'appliances bioinformatiques (RAINBio[1]), et permet d'explorer ces appliances et outils d'après les termes de l'ontologie EDAM et une recherche textuelle.

Afin d'assurer un authentification simple, de confiance, et utilisable autant dans le portail que dans les appliances, nous utilisons le gestionnaire d'identités développé dans le cadre du projet CYCLONE. Ce gestionnaire s'appuie sur la fédération d'identités académique européenne eduGAIN.

Une fois l'appliance choisie, il est possible de la lancer en un clic, et de configurer son déploiement. Lorsque l'appliance est basée sur une unique machine virtuelle (VM), il est possible de choisir la mémoire et le nombre de CPU; et dans le cas d'un cluster, le nombre de noeuds de calcul. Le portail Biosphère permet aussi de choisir le cloud utilisé, le cas échéant de lancer l'appliance sur plusieurs clouds, et de modifier ensuite dynamiquement le nombre de noeuds de calcul.

Dans les déploiements, l'authentification s'appuie sur la fédération d'identités eduGain, tant pour les connexions en terminal distant, qu'en bureau distant. L'isolation réseau d'un déploiement, par exemple pourun cluster de calcul, est assurée par le composant CNSMO, développé dans le cadre de CYCLONE, fournissant un réseau privé virtuel (VPN), basé sur le logiciel communautaire OpenVPN (exécuté dans un conteneur Docker).

Le portail Biosphère permet une gestion simplifiée des ressources utilisées. Afin de rendre le cloud accessible à tous, le portail Biosphère s'appuie sur les développements du projet CYCLONE. Il propose ainsi un haut niveau de sécurité réalisé par une authentification reposant sur la fédération d'identités académique européenne eduGAIN, et un placement des appliances au sein de réseaux isolés et sécurisés. Le portail permet aussi une allocation dynamique des ressources de calcul, et des déploiements multi-cloud. L'ensemble de ces fonctionnalités étant proposées dans une interface se voulant intuitive.

## References

[1] GROSJEAN Marie, HÉRIVEAU Claudia, JULLIEN Renaud, COLLIN Olivier, GIBRAT Jean-François, and BLANCHET Christophe. A RAINBio over the Life Sciences Cloud. Post-118 JOBIM 2015
# Biodjango, an open framework for bioinformatics publishing

Ennys GHEYOUCHE, Stéphane TÉLETCHÉA

¹ UFIP, UMR 6286 CNRS Université de Nantes, 2 rue de la Houssinière, 44322 Nantes cedex 3

Corresponding Author: stephane.teletchea@univ-nantes.fr

Most of recent bioinformatics methods are available to the scientific community through web-based portals. Without technical knowledge, users can rapidly evaluate or use the method presented. This simplicity is often a tradeoff between the focused service provided and the integration with external resources. These web servers also lack of consistency for the presentation of services, for the organization of data and results.

With the advent of specialized editions of journals, some guidelines start to be set up for web services. Users nowadays expect to see a description of the method, a contact page, a simple form for job submission and a results page with sample examples. These web services are still mostly independent of external resources.

In this work, we present Biodjango, an integration of the Django framework for the presentation of web services for bioinformatics methods. Within this extensive framework, we provide integrated methods for linking external biological data to the analysis performed. On top of biopython or biodjango-specific methods, it is possible to handle Uniprot entries for protein annotation, Gene Ontology vocabulary, NCBI data and Bibliography management. For each biodjango application, a reference template is provided with examples so users building their own web service can rapidly adapt these applications for their needs. The application is made for modularity so users can pick only sub-parts of biodjango as required. To simplify the management of job submission, progression and display of results, biodjango offers their management with a simplified scheduler-like mechanism. To ease user adoption of biodjango, an extensive documentation is provided for a rapid set up of a biodjango-derived web service, this documentation will be updated regularly from the remarks and demands of the community.

We expect the biodjango project to accelerate web deployment of web services for the scientific community, so bioinformaticians will be able to dedicate more time for the development of innovative methods.
# New Generation Phylogeny.fr: Refactoring Phylogeny.fr for Innovative Phylogenetic Services

Damien CORREIA<fnref n="1" />, Vincent LEFORT<fnref n="2" />, Olivia DOPPELT-AZEROUAL<fnref n="3" />, Fabien MAREUIL<fnref n="3" />,

Sarah COHEN-BOULAKIA<fnref n="4" /> and Olivier GASCUEL<fnref n="2" />,<fnref n="3" />

Corresponding Author: damien.correia@pasteur.fr

## 1 Context

Phylogenetic analyses aim at reconstructing the evolutionary history of biological objects from molecules to species, and populations. Phylogeny plays a major role in highly diverse domains such as predicting biological functions or measuring the biodiversity of ecosystems. As a consequence, a plethora of approaches have been designed in various communities, resulting in a large number of programs available. Faced with the increasing need to perform phylogenetic analysis, and the difficulty for scientists to determine which program to use at each step of an analysis and how to combine the use of programs together, we designed Phylogeny.fr [1] almost ten years ago.

While Phylogeny.fr runs 50,000 data analysis per month, we are now faced with two major difficulties: (i) the series of programs used need to be updated and (ii) phylogeny.fr is used in contexts where its technology is reaching its limits (e.g., when simultaneously used by hundreds of students or when the server is used through batch scripts). Refactoring Phylogeny.fr is thus of paramount importance.

In the meantime, new solutions have emerged to help users manage their analyses: scientific workflow systems such as Galaxy [2] have reached a level of maturity making them particularly suitable for complex and large-scale analyses.

## 2 Sketch of the demonstration

In this demonstration, we introduce the first release of "NGPhylogeny.fr" (NGPhylogeny.pasteur.fr), developed under Affero GPL v2 licence within a python Web framework (Django), in which we have refactored phylogeny.fr by designing a scalable environment, an easy-to-use web interface and a series of modular Galaxy workflows to perform a large variety of phylogenetic analyses. All programs have been updated or replaced while some others have been added (such as *Noisy* to trim the alignment sites). Default parameter settings have also been revised. Our demonstration will be based on real datasets. We will show (i) how "NGPhylogeny.fr" can be used in a functional genomics context to quickly analyze large sets of protein superfamilies, (ii) how in depth studies can be launched and (iii) how "NGPhylogeny.fr" can be installed on a wide variety of configurations.

**Acknowledgements:** This project is supported by ANR-11-INBS-0013 - IFB.

## References

[1] Phylogeny.fr: robust phylogenetic analysis for the non-specialist. Dereeper A, Guignon V, Blanc G, Audic S, Buffet S, Chevenet F, Dufayard JF, Guindon S, Lefort V, Lescot M, Claverie JM, Gascuel O. *Nucleic Acids Res.* 2008, 36(Web Server issue):W465-9

[2] Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Goecks J, Nekrutenko A, Taylor J; Galaxy Team. *Genome Biol.* 2010;11(8):R86.

<footnote n="1">¹ Institut Français de Bioinformatique, IFB-Core - UMS 3601, Gif-sur-Yvette, France</footnote>

<footnote n="2">² Laboratoire d'Informatique, de Robotique et de Microélectronique de Montpellier, LIRMM - UMR 5506 CNRS et Université de Montpellier, France</footnote>

<footnote n="3">³ Unité de Bioinformatique Evolutive et Hub Bioinformatique et Biostatistique, C3BI - USR 3756 Institut Pasteur et CNRS, Paris, France</footnote>

<footnote n="4">⁴ Laboratoire de Recherche en Informatique, LRI - CNRS UMR 8623 et Université Paris-Saclay, Orsay, France</footnote>
Conférence invitée

Franca FRATERNALI

<figure><img src="image_18.png" /></figure>

Randall Division of Cellular and Molecular Biophysics, King's college, London,
United Kingdom

Unraveling the Good and the Bad in Protein Networks: Functional
versus Dysfunctional Interactions

In the last years, protein interactome comparisons have highlighted conserved modules that might represent common functional cores of ancestral origin. However, recent analyses of protein-protein interaction networks (PPINs) have led to a debate about the influence of the experimental method on the quality and biological relevance of these interaction data. It is crucial to know to what extent discrepancies between the networks of different species reflect sampling biases of the respective experimental methods, as opposed to topological features due to biological functionality. This requires new, precise and practical mathematical tools to quantify and compare the topological structures of networks at high resolution. To this end, we have studied the relationship between structured random graph ensembles and real biological signaling networks, focusing on the number of short loops in networks, which represent complexes in PPINs. By combination of a method for graph dynamics and an algorithm for loop counting, we estimated the relative importance of loops in biological networks compared to random graphs. We found that loops are a predominant feature of PPINs, suggesting that enrichment of their occurrence has a key functional role.

Nevertheless, one must keep in mind that not all the interactions between proteins result
in a functional role that benefits the cell. One example is protein aggregation, resulting in
neurotoxic assemblies that lead ultimately to cell death. We investigate in detail the case of
interactions between fragments of the Prion protein (PrP) constituted by only the helices
H2 and H3 of the entire protein. We have investigated the molecular mechanisms of the
self-assembly process in solution by Molecular Dynamics. Our simulations show that this
process occurs by assembly of small modules of four monomers that precede the creation of
a base of six to eight H2H3 monomers; starting from this base, other H2H3 units attach to
it in various configurations, assembling short filaments.

1) Chung SS, Pandini A, Annibale A, Coolen AC, Thomas NS, Fraternali F. Bridging topo-logical and functional information in protein interaction networks by short loops profiling. Sci Rep. 2015; 5:8540.

2) Chakroun N, Fornili A, Prigent S, Kleinjung J, Dreiss CA, Rezaei H, Fraternali F. De-crypting Prion Protein Conversion into a β-Rich Conformer by Molecular Dynamics. J Chem Theory Comput. 2013; 5:2455-2465.

3) Chakroun N, Prigent S, Dreiss CA, Noinville S, Chapuis C, Fraternali F, Rezaei H. The oligomerization properties of prion protein are restricted to the H2H3 domain. FASEB J. 2010; 23:222-31.
# Use of cross-docking simulations for identification of protein-protein interactions sites: the case of proteins with multiple binding sites

Nathalie LAGARDE¹, Lydie VAMPARYS¹, Benoist LAURENT¹, Alessandra CARBONE²,³ and Sophie SACQUIN-MORA¹

¹ Laboratoire de Biochimie Théorique, CNRS UPR9080, Institut de Biologie Physico-Chimique, 13 rue Pierre et Marie Curie, 75005, Paris, France

² Laboratoire de Biologie Computationnelle et Quantitative, CNRS UMR7238, UPMC Univ-Paris 6, Sorbonne Université, 15 rue de l'Ecole de Médecine, 75006, Paris, France

³ Institut Universitaire de France, 75005, Paris, France

Corresponding Author: lagarde@ibpc.fr, sacquin@ibpc.fr

## 1 Introduction

Understanding protein-protein interactions (PPI) is essential to decipher the mechanism of numerous biological functions. Some *in silico* methods can be used to investigate PPI. In particular, cross-docking simulations of large datasets of proteins can be used to predict interface residues [1-3]. In this study, we discussed the ability of the cross-docking method to detect the multiple binding sites on protein surfaces.

## 2 Methods

358 proteins extracted from 138 unique PDB structures were used for this study. The MAXDo algorithm [1] was used with a rigid-body docking approach and a reduced protein representation, a coarse-grain protein model developed by Zacharias [4]. Binding site predictions resulting from evolutionary sequence analysis produced with JET [5] were used to restrict the initial search space. For each surface residue, its Protein Interface Propensity (PIP) was computed and used to predict binding sites on the protein surface.

## 3 Results

For a large number of proteins, alternative interfaces different from the reference experimental interfaces were predicted. However, about 70 % of these interfaces were not false positives but correspond to interfaces with other partners (other chain of the same PDB not included in the database, nucleic acid molecule or homo-/hetero-dimerization interface). We compared the use of two different scoring schemes accounting for multiple binding sites, for evaluating the binding sites prediction. The first score was obtained by comparing the predicted interface with one single global reference experimental interface generated by concatening all the existing experimental interfaces. In the second score, the predicted interface was compared to each experimental interface separately, and only the interface associated with the best performance was kept.

## 4 Conclusion

Using cross-docking simulations on a large dataset of proteins, accurate binding sites predictions could be realized, including proteins which present multiple binding sites.

## References

[1] Sophie Sacquin-Mora, Alessandra Carbone, and Richard Lavery. Identification of Protein Interaction Partners and Protein-Protein Interaction Sites. *Journal of Molecular Biology*, (382):1276-1289, 2008.

[2] Anne Lopes, Sophie Sacquin-Mora, Viktoriya Dimitrova, Elodie Laine, Yann Ponty and Alessandra Carbone. Protein-Protein Interactions in Crowded Environment: An Analysis via Cross-Docking Simulations and Evolutionary Information. *PLoS Computational Biology*, (9/12):e1003369,1-18, 2013.

[3] Lydie Vamparys, Benoist Laurent, Alessandra Carbone and Sophie Sacquin-Mora. Great interactions: How binding incorrect partners can teach us about protein recognition and function. *Proteins*, (84/10):1408-1421, 2016.

[4] Martin Zacharias. Protein-protein docking with a reduced protein model accounting for side-chain flexibility. *Protein Science*, (12):1271-1282, 2003.

[5] Stefan Engelen, Ladislas A Trojan, Sophie Sacquin-Mora, Richard Lavery and Alessandra Carbone. Joint Evolutionary Trees: A Large-Scale Method To Predict Protein Interfaces Based on Sequence Sampling. *PLoS Computational Biology*, (5/1):1-17, 2009.
*In silico* developments for the study of glycosylation applied to extracellular matrix proteins

Camille BESANÇON¹, Alexandre SÉBASTIEN BLAISE¹, Manuel DAUCHEZ¹,², Jessica JONQUET¹, Nicolas BELLOY¹,², Stéphanie BAUD¹,²

¹ UMCR CNRS/URCA N° 7369, Université de Reims Champagne-Ardenne, UFR Sciences Exactes et Naturelles, Moulin de la Housse, 51687 Reims Cedex 2, France

² Plateau de Modélisation Moléculaire Multiéchelle, Université de Reims Champagne-Ardenne, UFR Sciences Exactes et Naturelles, Moulin de la Housse, 51687 Reims Cedex 2, France

Corresponding Author: camille.besancon@etudiant.univ-reims.fr

N-glycosylations play an important role in protein functions and some alterations of glycosylation such as sialic-acid hydrolysis can alter protein activity. Some alterations are also known as pathological markers. Considering their biological importance, it is essential to study glycosylated proteins. However, *in vitro* study of N-glycans can be very difficult because of the structural diversity and the many reactive groups of the glycan chains. Molecular dynamics can be a useful tool to overcome this problem and give access to conformational informations through exhaustive sampling. At the beginning of this project, we lacked a tool, enabling to realize in a single workflow, from the building of a glycan structure to the analysis of the trajectories. To adress this issue, we want to use these studies to propose a complete tool for the building, the simulation and the analysis of glycosylated proteins that will be implemented in GROMACS.

The present work is based on a previous study from our laboratory: we investigated the structural influence of sialic-acid on N-glycan combining *in vitro* and *in silico* approaches. A new method to evaluate the conformational landscape has been developed by projecting 3D vectors on a geometric plan. We aligned the glycan core on the Z axis in an orthogonal space before projecting the position of the glycan branches on the XY plan. With this method, along with the clustering of representative structures and the measurement of dihedral angles values of the glycosidic linkage, it was possible to show the importance of sialic acids and their role in the structural flexibility of N-glycans [1].

To investigate the structural influence of glycosylation on protein structure, we focus on the glycosylated Insulin Receptor (IR). It was previously shown that the hydrolysis of sialic acids on the IR by neuraminidase-1 sialidase induced insulin resistance and disrupted cell glucose uptake [2]. Simultaneously, we will work on new methods to specifically analyze protein / glycan complexes. We're aiming at developing a method for the measurement of the protein surface covered by the glycan in order to investigate the structural role of the sugar chains on protein surfaces.

At this time, using the online tool CHARMM-GUI [3] to generate the structures, we have created a first glycan library that will be implemented in our GROMACS module. We first concentrated on biologically relevant N-glycans (pathological markers identified by literature search [4,5]). We then started molecular dynamics simulations on some of these glycans. Our aim is to further our analysis and our results on the matter of N-glycan structure and flexibility. With this data set, we will be able to start implementing the GROMACS module and test our analysis methods. From this work, we aim at providing a complete tool paired with specific analysis methods thus enabling the building, the simulation of glycosylated proteins.

References

[1] Guillot, A., Dauchez, M., Belloy, N., Jonquet, J., Duca, L., Romier, B., Maurice, P., Debelle, L., Martiny, L., Durlach, V., et al. Impact of sialic acids on the molecular dynamic of bi-antennary and tri-antennary glycans. *Sci. Rep.* **6**, 35666, 2016.

[2] Hayes, G. R. & Lockwood, D. H. The role of cell surface sialic acid in insulin receptor function and insulin action. *J. Biol. Chem.* **261**, 2791–2798 (1986).

[3] Jo, S., Song, K.C., Desaire, H., MacKerell, A.D., and Im, W. (2011). Glycan Reader: Automated Sugar Identification and Simulation Preparation for Carbohydrates and Glycoproteins. *J. Comput. Chem.* **32**, 3135–3141.

[4] Kawar, Z.S., Haslam, S.M., Morris, H.R., Dell, A., and Cummings, R.D. (2005). Novel poly-GalNAcbeta1-4GlcNAc (LacdiNAc) and fucosylated poly-LacdiNAc N-glycans from mammalian cells expressing beta1,4 Nacetylgalactosaminyltransferase and alpha1,3-fucosyltransferase. *J. Biol. Chem.* **280**, 12810–12819.

[5] Lee, L.Y., Thaysen-Andersen, M., Baker, M.S., Packer, N.H., Hancock, W.S., and Fanayan, S. (2014). Comprehensive Nglycome profiling of cultured human epithelial breast cells identifies unique secretome N-glycosylation signatures enabling tumorigenic subtype classification. *J. Proteome Res.* **13**, 4783–4795.
# HG-CoLoR: Hybrid Graph for the error Correction of Long Reads

Pierre MORISSE, Thierry LECROQ and Arnaud LEFEBVRE

Normandie-Université, UNIROUEN, LITIS EA4108, 76821, Mt-St-Aignan, France

Corresponding author: pierre.morisse2@univ-rouen.fr

**Abstract** The recent rise of long read sequencing technologies allows the solving of assembly problems for large and complex genomes that were, until then, unsolvable with the use of short read sequencing technologies alone. Despite the fact that they can reach lengths of tens of kbps, these long reads are very noisy, and can reach an error rate as high as 30%, involving mandatory error correction before using them to efficiently solve assembly problems. However, as the vast majority of these errors are insertions and deletions, classical error correction tools developed for short reads, which mainly focus on substitution errors, are not effective for correcting long reads. Therefore, several new methods specifically designed for long read error correction have recently been developed. In particular, NaS, instead of directly correcting the long reads, proposes to use them as templates in order to produce assemblies of related accurate short reads, and use them as corrections. Following this idea, we introduce HG-CoLoR, a new tool for the production of such corrected long reads, that gets rid of the need to align the short reads against each other, which is the bottleneck from NaS. Indeed, HG-CoLoR focuses on a seed-and-extend approach based on a hybrid graph built from the short reads. Our experiments show that, while producing comparable results both in terms of length and accuracy of the corrected long reads, HG-CoLoR is several times faster than NaS, and also yields better assembly results than other state-of-the-art long read hybrid error correction methods. HG-CoLoR is available from https://github.com/pierre-morisse/HG-CoLoR.

**Keywords** NGS, long reads, correction, assembly

## 1 Introduction

Since a few years, long read sequencing technologies are being developed, and allow the solving of assembly problems for large and complex genomes that were, until then, unsolvable with the use of short reads sequencing technologies alone. The two major actors of these long read sequencing technologies are Pacific Biosciences and Oxford Nanopore, which, with the release of the MinION device, allows a low-cost and easy long read sequencing.

However, even though long reads can reach lengths of tens of kbps, they also reach a very high error rate of around 15% for Pacific Biosciences, and up to 30% for Oxford Nanopore, the vast majority of these errors being insertions and deletions. Due to this high error rate, correcting these long reads before using them to efficiently solve assembly problems is mandatory. Many methods are available for short read correction, but these methods are not applicable to long reads, on the one hand because of their much higher error rate, and on the other hand, because most of the error correction tools for short reads focus on substitution errors, the dominant error type in Illumina data, whereas insertions and deletions are more common in long reads.

Recently, several methods for long read correction have been developed. These methods can be divided into two main categories: either the long reads are selfcorrected by aligning them against each other (HGAP [1], PBcR [2]), or either a hybrid strategy is adopted, in which the long reads are corrected with the help of accurate short reads (LSC [3], proovread [4], CoLoRMap [5]). de Bruijn graph [6] based methods, where the long reads are mapped on the graph, and erroneous regions corrected by traversing its paths, also started to develop recently, in the hybrid case (LoRDEC [7], Jabba [8]), as well as in the non-hybrid case (LoRMA [9]).

NaS [10], instead of locally correcting the long reads, uses them as templates to produce corrected long reads from assemblies of related accurate short reads. The short reads are mapped both on these templates, and against each other, in order to associate a subset of related short reads to each template. Each subset in then assembled, and the obtained contig is used as the correction of the associated template.

In this paper, we introduce HG-CoLoR, a new long read hybrid error correction method that combines both the main idea from NaS to produce corrected long reads from assemblies of accurate short reads, and the use
of a graph, in order to get rid of the time consuming step of aligning all the short reads against each other. HG-CoLoR indeed focuses on a seed-and-extend approach where the seeds, which are short reads that align correctly on the long reads, are used as anchor points on a graph that is traversed in order to link them together and to produce the corrected long reads. This graph is actually a hybrid structure between a de Bruijn graph and an overlap graph [12]. It is defined from the short reads' *k*-mers, and allows to compute perfect overlaps of variable length between these *k*-mers. It is not explicitly built, but the use of PgSA [11] allows to simulate its traversal.

Our experiments show that, while producing comparable results both in terms of length and accuracy of the corrected long reads, HG-CoLoR is several times faster than NaS, and also yields better assembly results than other state-of-the-art long read hybrid error correction methods.

For the sake of understanding, we first give an overview of NaS, and describe our hybrid graph and the way PgSA allows its traversal, before introducing HG-CoLoR.

## 2 NaS Overview

NaS is a hybrid method for the error correction of long reads that, unlike other methods, uses long reads as templates rather than directly correcting them. Short reads are mapped both on these templates and against each other, in order to gather different subsets of short reads, each subset related to one given template. Each subset is then assembled, and the produced contig is used as the correction of the related template. More precisely, a corrected long read is produced from a template as follows.

First, the short reads are aligned on the template using BLAT [13] in fast mode, or LAST [14] in sensitive mode, in order to find seeds, which are short reads that align correctly on the template. Then, once these seeds have been found, all the short reads are aligned against each other, and similar reads, which are reads that share a certain number of non-overlapping *k*-mers with the seeds, are recruited with the help of Commet [15]. Finally, the obtained subset of short reads is assembled using Newbler (unpublished), and a contig is produced, and used as the correction of the initial template

The reads recruitment is the most crucial step of the method, as it allows to retrieve short reads corresponding to low quality regions of the template. However, this step is also the bottleneck of the whole NaS pipeline, as it is responsible for 70% of the total runtime on average.

NaS is able to generate corrected long reads up to 60 kbps, that align entirely on the reference genome and that span repetitive regions. On average, the accuracy of the corrected long reads produced by NaS reaches 99.75%, without any significant length drop compared to the input long reads. Moreover, these corrected long reads also yield highly contiguous assembly results, and thus show that focusing on the production of corrected long reads from assemblies of accurate short reads, rather than on the local correction of the long reads, is an interesting alternative to classical long read hybrid error correction.

## 3 Hybrid graph

As previously mentioned, the graph used by HG-CoLoR is a hybrid structure between a de Bruijn graph and an overlap graph. This hybrid graph is not explicitly built, but can be traversed with the help of PgSA, which is a data structure that allows the indexing of a set of reads of constant length, in order to answer different queries, for a given string *f*. For place sake, we do not detail how the index is built, the complete list of queries, nor how they are processed. For more details, one can refer to [11]. We simply mention that PgSA supports querying for variable lengths of *f* without recomputing the index, and that one of the queries returns the positions of all the occurrences of *f* in the different reads of the set.

This way, using PgSA to index a set of reads, and looping over the aforementioned query, allows to compute perfect overlaps of variable length between the reads, thus simulating the use of an overlap graph. In the same fashion, indexing the *k*-mers from a set of reads, and looping over the aforementioned query, fixing the length of the queries strings as *k* − 1, allows to compute perfect overlaps of length *k* − 1 between the *k*-mers, thus simulating the use of a de Bruijn graph. However, indexing the *k*-mers from a set of reads, and looping over the aforementioned query, of course, also allows to compute perfect overlaps of variable length between the different *k*-mers, thus simulating the use of a hybrid structure between a de Bruijn graph and an overlap graph. To the best of our knowledge, this is the first time such a structure is mentioned. For better understanding, an example of a simple graph is given in Figure 1.
<figure><img src="image_19.png" /><figcaption>**Fig. 1.** A simple example of our graph, when fixing the length of the *k*-mers to 6, computing overlaps of minimum length 3, and building from the three following reads: AAGCTTAC, CTTACGTA, GTATACTG. Numbers on the edges of the graph represent the overlap length between the *k*-mers.</figcaption></figure>

# 4 HG-CoLoR description

HG-CoLoR, like NaS, aims to use erroneous long reads as templates, and to produce corrected long reads from assemblies of short reads related to these templates. However, its main objective is to get rid of the time consuming step of reads recruiting, that requires the mapping of all the short reads against each other. To do so, it focuses on a seed-and-extend approach where the seeds are found in the same way as NaS, and where the *k*-mers from the short reads, and their reverse-complements, are indexed with PgSA, to allow the traversal of the previously described graph. This graph is traversed, in order to extend and link together the seeds, used as anchor points, by directly assembling the short reads' *k*-mers during the traversal. HG-CoLoR's workflow is summarized in Figure 2, and its four main steps are described below.

<figure><img src="image_20.png" /><figcaption>**Fig. 2.** HG-CoLoR's workflow. First, the short reads are corrected in order to get rid of as much sequencing errors as possible. Then, all the *k*-mers from the corrected short reads, and their reverse-complements, are obtained with Jellyfish, and indexed with PgSA, to allow the traversal of the graph. The corrected short reads are aligned on the long reads with BLASR to find seeds, and each long read is then considered as a template, and processed independently. For a given template, the graph is traversed in order to extend and link together the associated seeds, used as anchor points. Then, the tips of the sequence obtained after the seeds linking step are extended in both directions by traversing the graph, to reach the initial template's borders. Finally, the corrected long read is output.</figcaption></figure>

## 4.1 Short reads correction and indexing

Even though short reads are very accurate prior to any correction, as HG-CoLoR seeks to arrange their *k*-mers into a graph structure, and traverse it to extend and link the seeds together, it needs to get rid of as much sequencing errors as it can in this data. Thus, prior to any other step, the short reads are corrected with the help of QuorUM [16], which is able to provide a good raise of the accuracy in very little time. Then, the *k*-mers from the corrected short reads, and their reverse-complements, are extracted with Jellyfish [17], and indexed with PgSA, in order to allow the traversal of the graph during the following steps.

## 4.2 Seeds retrieving and merging

Like with NaS, the seeds are found by mapping the corrected short reads on the long reads, used as templates. This is done with the help of BLASR [18], an alignment tool specifically designed to align long reads dominated by insertion and deletion errors. Then, each template is processed independently, and two phases of analyze and merging are applied to the associated seeds. First, if the mapping positions of a given couple of
seeds imply that they overlap on the template over a sufficient length, their assumed overlapping sequences are compared, and the two seeds are merged accordingly. If the mapping positions indicate that the two seeds do overlap on the template, but not over a sufficient length, or if the assumed overlapping sequences do not coincide, only the seed with the best alignment score is kept. Then, once all the seeds with overlapping mapping positions have been merged or filtered out, sequence overlaps between consecutive seeds are computed. As in the previous step, if a given seed overlaps the following one over a sufficient length, the two seeds are merged.

## 4.3 Seeds linking

Once the seeds have been found and merged for all of the templates, HG-CoLoR once again processes each template independently and attempts to link the related seeds together by considering them as couples, and traversing the graph. The rightmost *k*-mer of the left seed (source) and the leftmost *k*-mer of the right seed (destination) are used as anchor points, and the source is extended with perfectly overlapping *k*-mers from the corrected short reads, found by following the paths of the graph, until the destination is reached. When facing branching paths, every possible path is explored with the use of backtracking, to find the one that will allow correct linking of the source to destination. Of course, HG-CoLoR explores these different paths in decreasing order of the overlaps lengths, which means that edges representing longer overlaps are always explored before those representing shorter ones. It also only explores edges that represent overlaps that are longer than a defined minimum length. Moreover, as short reads from a different region of the reference genome can align on the template and can be used as seeds, thus leading to impossible linkings, a threshold on the maximum number of backtracks is set, to avoid useless important runtime and intensive computation.

If this threshold is reached, and no path has been found to link the source to the destination, the current linking iteration is given up. When such a situation occurs, two different cases have to be taken into account. In the first case, if no seeds have been linked so far, the current source is simply ignored, and a new linking iteration is computed for the next couple of seeds. In the second case, if seeds have already been linked previously, the source remains the same, the destination seed that could not be reached is ignored, and the destination is defined as the next seed for the next linking iteration. An illustration of these two different cases is given in Figure 3.

However, in the second case, as this process of skipping a seed in the middle of the template can provoke an important number of failed linking attempts, if seeds from a wrong region are present in great proportion on the template, a threshold on the maximum number of seeds that can be skipped is set. Once this threshold is reached, if the sequence obtained from the previously linked seeds could not be extended to reach one of the remaining seeds, HG-CoLoR attempts to produce a fragmented corrected long read: the part corresponding to the seeds linked so far is output, and the graph is traversed again, in order to try to link the remaining seeds together, independently of the previous part.

<figure><img src="image_21.png" /><figcaption>**Fig. 3.** Illustration of the different cases of seed skips. Hatched lines represent the templates, standard segments represent the seeds, and bold segments represent the sequences obtained from the previously linked seeds. First case (left): No seeds have been linked so far, the current source seed is simply ignored, and both the source and the destination are moved to the next couple of seeds. Second case (right): Seeds have already been linked previously, the source remains the same, the destination seed that could not be reached is ignored, and the destination is defined as the next seed.</figcaption></figure>

## 4.4 Tips extension

Finally, it is obvious that the seeds do not always map right at the beginning and until the end of the templates. Thus, in order to get as close as possible to the original templates' lengths, once all the seeds of a given template have been linked, HG-CoLoR keeps on traversing the graph and extending the tips of the produced corrected long read, on the left of the leftmost seed, and on the right of the rightmost seed, until they reach the template's borders, or a branching path. Indeed, in the case of tips extension, when facing a
branching path, HG-CoLoR has no clue as to which path to chose and continue the extension with, nor any anchor points, unlike when it attempts to link two seeds together. Therefore, backtracking is useless and the extension is simply stopped when such a situation occurs. In the case of fragmented corrected long reads, as HG-CoLoR can not properly rely on the template's borders, every fragment is extended until a branching path is reached.

# 5 Results and discussion

We compare the quality of our corrected long reads with those produced by NaS, and also with those produced by two other state-of-the-art hybrid error correction methods, namely CoLoRMap and Jabba (which is more recent than LoRDEC). We compare the results both in terms of alignment identity of the corrected long reads, and in terms of quality of the assemblies that could be generated from these reads.

## 5.1 Parameters

We ran multiple rounds of correction with HG-CoLoR on the different datasets to experiment with the parameters, and find the combination that would produce the best results. Thereby, we found that a *k*-mer size of 64 for the graph construction yielded the best compromise between accuracy, genome coverage, and average length of the output corrected long reads. The minimum overlap length to allow the merging of two seeds during the second step was set to 63, accordingly to the *k*-mer size chosen for the graph construction. The minimum overlap length allowed to explore an edge during the graph traversal was set to 59, as decreasing it more yielded unsatisfying results, and increasing it would make our graph closer to an actual de Brujin graph than to the hybrid graph it's supposed to be. The maximum number of backtracks was set to 1,125, as decreasing it more drastically impacted the quality of the produced corrected long reads, and increasing it, even to very large values, barely yielded better results, but greatly increased the runtime. For the same reason, the maximum number of seed skips was set to 5. For the mapping of the short reads on the long reads, BLASR was used with default parameters except for bestn, that was set to 30 instead of 10. Yet again, increasing this parameter to larger values only impacted the runtime, and did not improve the correction results enough to be interesting, while decreasing it induced a drop of the number of output corrected long reads. Finally, GNU Parallel [19] was used to allow HG-CoLoR to run on multiple processes. CoLoRMap was run with default parameters. Following the authors' recommendations, before running Jabba, the short reads were corrected with Karect [20], and the de Brujin graph was constructed and corrected with Brownie (unpublished). As choosing the same value as the one used for HG-CoLoR led to worse results, a *k*-mer size of 75 was chosen for the graph construction, as recommended by the authors. All tools were run with 16 processes.

## 5.2 Datasets

As we mainly seek to compare our results with NaS, we use the same data to allow a better comparison. This data is composed of both long Oxford Nanopore reads and short Illumina reads for three different genomes: *Acinetobacter baylyi*, *Escherichia coli*, and *Saccharomyces cerevisae*. Details are given in Table 1.

<table><caption>**Tab. 1.** Description of the datasets used in our experiments. Both MinION and Illumina data are available from the Genoscope's website http://www.genoscope.cns.fr/externe/nas/datasets.html.</caption><thead><tr><th rowspan="2">Dataset</th><th rowspan="2">Name</th><th colspan="2">Reference genome</th><th colspan="3">Oxford Nanopore data</th><th colspan="3">Illumina data</th></tr><tr><th>Strain</th><th>Reference sequence</th><th>Genome size</th><th># Reads</th><th>Average length</th><th>Coverage</th><th># Reads</th><th>Read length</th><th>Coverage</th></tr></thead><tbody><tr><td>A. baylyi</td><td>A. baylyi</td><td>ADPI</td><td>CR543861</td><td>3.6 Mbp</td><td>89,011</td><td>4,284</td><td>44x</td><td>900,000</td><td>250</td><td>50x</td></tr><tr><td>E. coli</td><td>E. coli</td><td>K-12 substr.</td><td>MG1655</td><td>4.6 Mbp</td><td>22,270</td><td>5,999</td><td>28x</td><td>775,500</td><td>300</td><td>50x</td></tr><tr><td>Yeast</td><td>S. cerevisae</td><td>W303</td><td>scf7180000000084-113</td><td>12.4 Mbp</td><td>205,923</td><td>5,698</td><td>31x</td><td>2,500,000</td><td>250</td><td>50x</td></tr></tbody></table>

## 5.3 Alignment-based comparison

The previously described long reads datasets were aligned with Last prior to any correction. The four different correction tools were then applied, and the obtained corrected long reads were aligned with BWA mem [21]. Results are given in Table 2 and discussed below.

We notice that, unlike the other methods, CoLoRMap output all the long reads and not only the ones it managed to correct. As the reads that could be corrected were not tagged in any way and could therefore not be extracted, it appears that CoLoRMap performed the worst correction, and did not manage to improve the accuracy of the long reads at all, except for the *E. coli* dataset. These poor results are probably due to the
<table><caption>**Tab. 2.** Runtime and statistics of the long reads, before and after correction by the different tools. NaS runtimes are omitted for the Yeast dataset because the results did not compute in 16 days. NaS corrected reads for this dataset were obtained from the Genoscope website.</caption><thead><tr><th>Dataset</th><th>Method</th><th># Reads</th><th>Average length</th><th>Cumulatize size</th><th># Aligned reads</th><th>Average identity</th><th># Error-free reads</th><th>Genome coverage</th><th>Runtime</th></tr></thead><tbody><tr><td rowspan="6">A. baylyi</td><td>Original</td><td>89,011</td><td>4,284</td><td>381,365,755</td><td>29,954</td><td>70.09%</td><td>0</td><td>100%</td><td>N/A</td></tr><tr><td>CoLoRMap</td><td>89,011</td><td>4,355</td><td>387,609,994</td><td>18,085</td><td>67.93%</td><td>2</td><td>100%</td><td>14h3min</td></tr><tr><td>Jabba</td><td>17,476</td><td>10,260</td><td>179,309,738</td><td>17,476</td><td>99.40%</td><td>16,893</td><td>99.80%</td><td>12min30</td></tr><tr><td>NaS (fast)</td><td>24,063</td><td>8,840</td><td>212,707,189</td><td>24,063</td><td>99.82%</td><td>22,984</td><td>100%</td><td>94h18min</td></tr><tr><td>NaS (sensitive)</td><td>28,492</td><td>9,530</td><td>271,526,778</td><td>28,492</td><td>99.83%</td><td>27,190</td><td>100%</td><td>128h55min</td></tr><tr><td>HG-CoLoR</td><td>23,465</td><td>11,137</td><td>261,327,970</td><td>23,461</td><td>99.44%</td><td>20,906</td><td>100%</td><td>21h08min</td></tr><tr><td rowspan="6">E. coli</td><td>Original</td><td>22,270</td><td>5,999</td><td>133,607,392</td><td>22,170</td><td>79.46%</td><td>0</td><td>100%</td><td>N/A</td></tr><tr><td>CoLoRMap</td><td>22,270</td><td>6,219</td><td>138,489,144</td><td>21,784</td><td>89.02%</td><td>152</td><td>100%</td><td>8h26min</td></tr><tr><td>Jabba</td><td>22,065</td><td>5,794</td><td>127,848,525</td><td>22,065</td><td>99.81%</td><td>21,850</td><td>99.41%</td><td>12min56</td></tr><tr><td>NaS (fast)</td><td>21,818</td><td>7,926</td><td>172,918,739</td><td>21,818</td><td>99.86%</td><td>20,383</td><td>100%</td><td>72h02min</td></tr><tr><td>NaS (sensitive)</td><td>22,144</td><td>8,307</td><td>183,958,832</td><td>22,144</td><td>99.86%</td><td>20,627</td><td>100%</td><td>81h30min</td></tr><tr><td>HG-CoLoR</td><td>22,549</td><td>5,897</td><td>132,979,813</td><td>22,549</td><td>99.59%</td><td>19,676</td><td>100%</td><td>15h15min</td></tr><tr><td rowspan="6">Yeast</td><td>Original</td><td>205,923</td><td>5,698</td><td>1,173,389,509</td><td>68,215</td><td>55.49%</td><td>0</td><td>99.90%</td><td>N/A</td></tr><tr><td>CoLoRMap</td><td>205,923</td><td>5,737</td><td>1,181,298,941</td><td>40,530</td><td>39.93%</td><td>23</td><td>99.40%</td><td>37h36min</td></tr><tr><td>Jabba</td><td>36,958</td><td>6,613</td><td>244,402,749</td><td>36,855</td><td>99.55%</td><td>34,028</td><td>93.21%</td><td>44min05</td></tr><tr><td>NaS (fast)</td><td>71,793</td><td>5,938</td><td>426,326,355</td><td>71,664</td><td>99.59%</td><td>59,788</td><td>98.70%</td><td>-</td></tr><tr><td>NaS (sensitive)</td><td>85,432</td><td>6,770</td><td>578,351,588</td><td>85,288</td><td>99.53%</td><td>69,816</td><td>99.17%</td><td>-</td></tr><tr><td>HG-CoLoR</td><td>71,518</td><td>6,604</td><td>472,306,800</td><td>71,393</td><td>99.17%</td><td>55,357</td><td>98.39%</td><td>99h16min</td></tr></tbody></table>

fact that only a few reads could be corrected, as CoLoRMap is designed to correct long reads from Pacific Biosciences, that have an error rate of about 15%, whereas the long reads used in our experiments were from Oxford Nanopore, and reached an error of at least 30% for the two other datasets.

Jabba clearly performed the best when it comes to runtime, outperforming all the other tools by several orders of magnitude. It also produced corrected long reads that aligned with a high identity, a great proportion of them aligning with no error. However, although highly accurate, these corrected long reads did not manage to completely cover any of the studied reference genomes.

When it comes to this point, only NaS and HG-CoLoR managed to cover the whole reference genomes with high identity, except for Yeast, due to the fact that even the original long reads did not cover the whole genome. Moreover, HG-CoLoR outperforming Jabba in terms of genome coverage also tends to underline the usefulness of our hybrid graph, showing that it seems to resolve the different regions of the reference genomes better than a classical de Bruijn graph, when the short reads coverage is locally insufficient.

On the three datasets, NaS yielded more corrected long reads than HG-CoLoR, both in fast and sensitive mode. The slight advantage of HG-CoLoR on the *E. coli* dataset comes from the production of fragmented corrected long reads, rather than from a greater number of processed templates. In both modes, the corrected long reads produced by NaS also aligned with a slightly higher identity than those produced by HG-CoLoR, and a greater proportion was therefore error-free. As for the average length and the cumulative size of the corrected long reads, HG-CoLoR performances were highly similar to NaS's, except on the *E. coli* dataset, where the advantage of NaS is probably due to the high quality of the original templates, and to the fact that it can recruit short reads outside of the templates, while HG-CoLoR stops once the borders are reached. However, despite its slight disadvantage on the aforementioned metrics, HG-CoLoR was at least four times faster than NaS, even in fast mode.

## 5.4 Assembly-based comparison

All the corrected long reads datasets previously described were assembled using Canu [22], without the correction and trimming steps. The following parameters were used for the assembly of all the datasets: OvlMerSize=17, MhapMerSize=17, OvlMerDistinct=0.9925, OvlMerTotal=0.9925. The correctedErrorRate parameter was tuned independently for each dataset. It was set to 0.07 for *A. baylyi*, to 0.085 for *E. coli* and to 0.125 for Yeast. Results are given in Table 3 and discussed below.

In agreement with what we observed in Table 2, the low accuracy of the long reads corrected by CoLoRMap resulted in impossible assemblies. Only the corrected long reads of the *E. coli* dataset could be assembled, due to their original high accuracy, but the generated assembly did not cover the whole genome, and displayed the worst identity among all the other assemblies.

As for Jabba, the fact that the corrected long reads did not manage to cover the whole reference genomes resulted in highly fragmented assemblies, that could not resolve large regions of the reference genomes. As a
<table><caption>**Tab. 3.** Statistics of the assemblies that were generated from the long reads, after correction by the different tools. CoLoRMap results are omitted for the *A. baylyi* and Yeast datasets, because Canu did not manage to assemble the sets of corrected long reads.</caption><thead><tr><th>Dataset</th><th>Method</th><th># Reads</th><th>Coverage</th><th># Expected contigs</th><th># Obtained contigs</th><th>Genome coverage</th><th>Identity</th></tr></thead><tbody><tr><td rowspan="5">A. baylyi</td><td>CoLoRMap</td><td>89,011</td><td>44x</td><td>1</td><td>-</td><td>89.43%</td><td>99.93%</td></tr><tr><td>Jabba</td><td>17,476</td><td>50x</td><td>1</td><td>13</td><td>100%</td><td>99.99%</td></tr><tr><td>NaS (fast)</td><td>24,063</td><td>59x</td><td>1</td><td>1</td><td>99.72%</td><td>99.98%</td></tr><tr><td>NaS (sensitive)</td><td>28,492</td><td>75x</td><td>1</td><td>2</td><td>99.97%</td><td>99.93%</td></tr><tr><td>HG-CoLoR</td><td>23,465</td><td>73x</td><td>1</td><td>1</td><td></td><td></td></tr><tr><td rowspan="5">E. coli</td><td>CoLoRMap</td><td>22,270</td><td>28x</td><td>1</td><td>29</td><td>97.74%</td><td>99.81%</td></tr><tr><td>Jabba</td><td>22,065</td><td>28x</td><td>1</td><td>41</td><td>95.76%</td><td>99.92%</td></tr><tr><td>NaS (fast)</td><td>21,818</td><td>37x</td><td>1</td><td>1</td><td>99.90%</td><td>99.99%</td></tr><tr><td>NaS (sensitive)</td><td>22,144</td><td>40x</td><td>1</td><td>2</td><td>100%</td><td>99.99%</td></tr><tr><td>HG-CoLoR</td><td>22,549</td><td>29x</td><td>1</td><td>2</td><td>99.95%</td><td>99.95%</td></tr><tr><td rowspan="5">Yeast</td><td>CoLoRMap</td><td>205,923</td><td>14x</td><td>30</td><td>-</td><td>70.52%</td><td>99.83%</td></tr><tr><td>Jabba</td><td>36,958</td><td>21x</td><td>30</td><td>134</td><td>97.44%</td><td>99.77%</td></tr><tr><td>NaS (fast)</td><td>71,793</td><td>35x</td><td>30</td><td>123</td><td>96.98%</td><td>99.80%</td></tr><tr><td>NaS (sensitive)</td><td>85,432</td><td>47x</td><td>30</td><td>123</td><td>92.19%</td><td>99.61%</td></tr><tr><td>HG-CoLoR</td><td>71,518</td><td>39x</td><td>30</td><td>108</td><td></td><td></td></tr></tbody></table>

result, long reads corrected by Jabba yielded the least covering assemblies, despite their high average length and high accuracy. This underlines the fact that, although it is extremely fast, Jabba does not seem to be adapted for correcting long reads prior to an assembly.

Surprisingly, for all the datasets, the sensitive mode of NaS produced corrected long reads that resulted in slightly less satisfying assemblies than the fast mode. However, the difference was not significant, and adapting the parameters of Canu to match the corrected long reads produced in sensitive mode addressed this issue.

Therefore, only the corrected long reads produced by NaS and HG-CoLoR could be assembled into a decent number of contigs, covering the reference genomes well, and with a high identity. However, for the Yeast dataset, none of these two tools managed to produce corrected long reads allowing to get close to the expected number of contigs, nor to the full genome coverage. This is probably due to the fact that the original long reads were of really poor quality, displaying an error rate of almost 45%, and did not cover the whole genome. They were indeed sequenced with an old chemistry, and it is more than likely that, with long reads from a more recent one as templates, both NaS and HG-CoLoR could produce corrected long reads that would greatly reduce the number of contigs and increase the genome coverage of the assembly.

## 6 Conclusion

We described HG-CoLoR, a new hybrid method for the error correction of long reads, that, like NaS, uses long reads as templates and focuses on the production of corrected long reads from assemblies of accurate short reads, rather than on the local correction of the input long reads. Our method, instead of aligning the short reads against each other in a recruiting step, like NaS, focuses on a seed-and-extend approach and introduces a brand new idea of using a hybrid structure between a de Bruijn graph and an overlap graph. This graph, which is defined from the short reads' *k*-mers, and traversed with PgSA, is used to extend and link together the seeds, which are short reads that align correctly on the input long reads, by a simple traversal, using them as anchor points. Therefore, the corrected long reads are produced by directly assembling the short reads' *k*-mers during the traversal, without using any other proper assembly tool.

We tested this new method and compared it with NaS, CoLoRMap and Jabba on Oxford Nanopore long reads from three different genomes, namely *A. baylyi*, *E. coli*, and *S. cerevisae*. On these three datasets, HG-CoLoR yielded results that compared well with NaS, while being several times faster, CoLoRMap produced corrected long reads of poor quality, and Jabba, while being the fastest tool, produced accurate corrected long reads that however did not cover the whole reference genomes. As a result, only the corrected long reads produced by NaS and HG-CoLoR could be assembled into a decent number of contigs, covering well the reference genomes, although NaS outperformed HG-CoLoR on the *S. cerevisae* dataset.

The development of this method shows that, when having anchor points, the previously introduced hybrid
graph can prove useful for hybrid error correction of long reads, and can even yield better results than a classical de Bruijn graph. For future works, it could be interesting to focus more on this graph, and directly build it instead of simulating its traversal with PgSA, in order to directly map the long reads on the graph, like Jabba, thus skipping the alignment step of the short reads on the long reads, and reducing the runtime.

## Acknowledgements

The authors would like to thank the Genoscope team for the availability of all the data used in this paper. This work has been partially supported by the Defi MASTODONS C3G project from CNRS.

## References

[1] Chen-Shan Chin, David H Alexander, Patrick Marks, Aaron A Klammer, James Drake, Cheryl Heiner, Alicia Clum, Alex Copeland, John Huddleston, Evan E Eichler, Stephen W Turner, and Jonas Korlach. Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data. *Nature Methods*, 10(6):563–569, 2013.

[2] Konstantin Berlin, Sergey Koren, Chen-Shan Chin, James P Drake, Jane M Landolin, and Adam M Phillippy. Assembling large genomes with single-molecule sequencing and locality-sensitive hashing. *Nature biotechnology*, 33(6):623–630, 2015.

[3] Kin Fai Au, Jason G. Underwood, Lawrence Lee, and Wing Hung Wong. Improving PacBio Long Read Accuracy by Short Read Alignment. *PLoS ONE*, 7(10):1–8, 2012.

[4] Thomas Hackl, Rainer Hedrich, Jörg Schultz, and Frank Förster. Proovread: Large-scale high-accuracy PacBio correction through iterative short read consensus. *Bioinformatics*, 30(21):3004–3011, 2014.

[5] Ehsan Haghshenas, Faraz Hach, S Cenk Sahinalp, and Cedric Chauve. CoLoRMap: Correcting Long Reads by Mapping short reads. *Bioinformatics*, 32(17):i545–i551, 2016.

[6] Nicolaas Govert de Bruijn. A combinatorial problem. *Koninklijke Nederlandse Akademie van Wetenschappen te Amsterdam*, 49(7):758–764, 1946.

[7] Leena Salmela and Eric Rivals. LoRDEC: Accurate and efficient long read error correction. *Bioinformatics*, 30(24):3506–3514, 2014.

[8] Giles Miclotte, Mahdi Heydari, Piet Demeester, Stephane Rombauts, Yves Van de Peer, Pieter Audenaert, and Jan Fostier. Jabba: hybrid error correction for long sequencing reads. *Algorithms Mol Biol*, 11:10, 2016.

[9] Leena Salmela, Riku Walve, Eric Rivals, and Esko Ukkonen. Accurate selfcorrection of errors in long reads using de Bruijn graphs. *Bioinformatics*, 33(6):799–806, 2017.

[10] Mohammed-Amin Madoui, Stefan Engelen, Corinne Cruaud, Caroline Belser, Laurie Bertrand, Adriana Alberti, Arnaud Lemainque, Patrick Wincker, and Jean-Marc Aury. Genome assembly using Nanopore-guided long and error-free DNA reads. *BMC Genomics*, 16:327, 2015.

[11] Tomasz Kowalski, Szymon Grabowski, and Sebastian Deorowicz. Indexing arbitrary-length k-mers in sequencing reads. *PLoS ONE*, 10(7):1–14, 2015.

[12] Andrzej Ehrenfeucht, Tero Harju, Ion Petre, David M Prescott, and Grzegorz Rozenberg. *Overlap Graphs*, pages 99–108. Springer Berlin Heidelberg, Berlin, Heidelberg, 2004.

[13] W James Kent. BLAT — The BLAST -Like Alignment Tool. *Genome research*, 12:656–664, 2002.

[14] Szymon M Kielbasa, Raymond Wan, Kengo Sato, Szymon M Kiebasa, Paul Horton, and Martin C Frith. Adaptive seeds tame genomic sequence comparison. *Genome Research*, 21(3):487–493, 2011.

[15] Nicolas Maillet, Guillaume Collet, Thomas Vannier, Dominique Lavenier, and Pierre Peterlongo. Commet: Comparing and combining multiple metagenomic datasets. In *IEEE International Conference on Bioinformatics and Biomedicine (BIBM)*, 2014.

[16] Guillaume Marçais, James A Yorke, and Aleksey Zimin. QuorUM: An Error Corrector for Illumina Reads. *PLOS ONE*, 10(6):1–13, 2015.

[17] Guillaume Marçais and Carl Kingsford. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. *Bioinformatics*, 27(6):764–770, 2011.

[18] Mark J Chaisson and Glenn Tesler. Mapping single molecule sequencing reads using basic local alignment with successive refinement (BLASR): application and theory. *BMC bioinformatics*, 13(1):238, 2012.

[19] Ole Tange. GNU Parallel - The Command-Line Power Tool. :login: *The USENIX Magazine*, 36(1):42–47, 2011.

[20] Amin Allam, Panos Kalnis, and Victor Solovyev. Karect: accurate correction of substitution, insertion and deletion errors for next-generation sequencing data. *Bioinformatics*, 31(21):3421–3428, 2015.

[21] Heng Li and Richard Durbin. Fast and accurate long-read alignment with Burrows-Wheeler transform. *Bioinformatics*, 26(5):589–595, 2010.

[22] Sergey Koren, Brian P. Walenz, Konstantin Berlin, Jason R. Miller, and Adam M. Phillippy. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. *bioRxiv*, 2016.
# PhylOligo: a package to identify contaminant or untargeted organism sequences in genome assemblies.

Ludovic Mallet¹†, Tristan Bitard-Feildel²†, Franck Cerutti¹ and Hélène Chiapello¹

¹ INRA UR875, Unité Mathématiques et Informatique Appliquées de Toulouse (MIAT),
Auzeville, 31326 Castanet-Tolosan, France

² CNRS UMR7590, Sorbonne Universités, Université Pierre et Marie Curie - Paris 6 - MNHN -
IRD - IUC, Paris, France.

Corresponding author: ludovic.mallet@inra.fr

**Reference paper:** Mallet et al. (2017) *PhylOligo: a package to identify contaminant or untargeted organism sequences in genome assemblies*. in rev.

**Abstract** Genome sequencing projects sometimes uncover more organisms than expected, especially for complex and/or non-model organisms. It is therefore useful to develop software to identify mix of organisms from genome sequence assemblies. We developed a suite of tools to tackle an often overlooked question: how to deal with sequences from untargeted or unexpected organisms (i.e. contaminants but also organelles, commensal or parasitic organisms) in genomic data?

**Availability:** The tools are written in Python3 and R under the GPLv3 Licence and can be found at https://github.com/itsmeludo/Phyloligo/.

**Keywords** multi-species, contaminant filtering, untargeted sequencing, oligonucleotide signature, metagenomics

## 1 Introduction

We present PhylOligo, a new package including tools to explore, identify and extract organism-specific sequences in a genome assembly using the analysis of their DNA compositional characteristics: the oligonucleotide signature.

Compared to existing software, PhylOligo provides several features to explore assemblies including: (i) a customisable oligonucleotide pattern (ii) an interactive cladogram-based visualisation of the contig signature similarity and cumulative size to explore the signature clusters and profile putative additional materials (iii) a species-specific sequence filtering based on a supervised learning of candidate profiles and a double threshold scan.

Our strategy includes 3 main steps:

a) **Assembly exploration** using an **interactive tree visualisation** based on oligonucleotide profiles computed from all genomic contigs. PhylOligo allows for a visual exploration of the compositional similarity distribution and structure of the contigs in an assembly based on either continuous (kmers) or spaced-pattern oligonucleotide frequencies. The oligonucleotide profile of each contig is computed and a pairwise distance matrix based is produced (Figure 1A) to generate an interactive Neighbour-Joining tree. Branch width is drawn proportional to the cumulated length of the contigs in a clade, allowing the user to track where the main part of the assembly clusters (assumed to correspond to the targeted organisms) and what significant clades branch out as hint for separate organisms (see Figure 1B). Thanks to the Ape package, sequences from a clade are interactively selected on the tree and exported to learn a prototype of their oligonucleotide profile.

b) **Oligonucleotide profile prototype learning** based on contig subsets selected by the user at nodes of the tree. ContaLocate then allows the learning of oligonucleotide profiles of the main and presumed additional organisms identified and sampled by the user at the previous step.

c) **Assembly partitioning** to locate organism-specific regions and classify contigs or segments according to the learned prototypes. The assembly is then scanned with sliding windows to locate organism-specific regions using oligonucleotide divergences computed against the targeted and the additional profiles. The distribution of the divergence against both is used to establish two thresholds best separating the different modes in the density
<figure><img src="image_22.png" /><figcaption>**Fig. 1.** Visualisation and interactive exploration of assemblies. A: Pairwise compositional divergence of contigs produced by PhylOligo. Contigs are reordered by hierarchical clustering. B: Contig tree produced by PhylOligo on the tardigrade genome. The clade in red is the current selection pointed by the user. C: Contigs clustered by HDBSCAN on oligonucleotide frequencies, Data from *Magnaporthe oryzae*. Red and blue are predicted clusters, grey are unclassified. The hyperspace is reduced to 2 dimensions with t-SNE. D: Determination of the untargeted threshold in ContaLocate based on the distribution of distances between the untargeted clade and the scanning windows over the whole assembly.</figcaption></figure>

functions (See Figure 1D). Genomic regions with a divergence simultaneously crossing respective thresholds to the targeted and to the additional profiles are labelled as part of the additional organism and exported as a GFF file.

## 2 Results

Our strategy present several advantages. 1) Unlike sequence homology based methods, PhylOligo allows the identification of putative uncharacterised and distantly related sequences in assemblies. 2) The double threshold species-specific filtration prevents the removal of HGTs and the subsequent fragmentation of the assembly. 3) Chimeric contigs are detected and split. 4) Compared to filtering unassembled reads, learning the compositional profile on contigs allows for a refined profile, thus leading to a finer filtering.

PhylOligo has been successfully applied to identify untargeted bacterial organisms in four fungi genomic datasets. The software also identified additional organisms in the scaffolds of the tardigrade assembly.

## Acknowledgements

We thank the INRA bioinformatics platforms MIGALE and GenoToul for resources and Thomas Schiex and Natalie Villa-Vialaneix for their input.
# Dynamix: Dynamic visualization by automatic selection of informative tracks from hundreds of genomic data sets

Matthias MONFORT¹, Eileen E.M. FURLONG¹ and Charles GIRARDOT¹

¹ European Molecular Biology Laboratory, Genome Biology Unit, Meyerhofstr. 1, D-69117, Heidelberg, Germany

Corresponding Author: matthias.monfort@embl.de

Visualization of genomic data is fundamental for gaining insights into genome function. Yet, co-visualization of a large number of data sets remains a challenge in all popular genome browsers and the development of new visualization methods is needed to improve the usability and user experience of genome browsers.

The Furlong Laboratory presents Dynamix [1], a genome browser plugin for JBrowse [2] that enables the parallel inspection of hundreds of genomic data sets. Dynamix takes advantage of a priori knowledge to automatically display data tracks with signal within a genomic region of interest. As the user navigates through the genome, Dynamix automatically updates data tracks and limits all manual operations otherwise needed to adjust the data visible on screen. Dynamix also introduces a new screen utilization by enabling users to independently scroll through groups of tracks.

Dynamix can be experimented navigating the 4CBrowser (http://furlonglab.embl.de/4CBrowser), a companion website released with Ghavi-Helm et al. 2014 [3], and the DynamixDemo server that demonstrates how Dynamix can be used for rich displays (http://furlonglab.embl.de/DynamixDemo). Visualization of genomics data sets remains a challenging area of research, we hope that Dynamix will encourage the development of innovative visualization methods.

Dynamix is hosted at http://furlonglab.embl.de/Dynamix under the MIT licence and has been published in *Bioinformatics* [1].

## References

[1] Matthias Monfort, Eileen E. M. Furlong, Charles Girardot; Dynamix: dynamic visualization by automatic selection of informative tracks from hundreds of genomic datasets. *Bioinformatics*, btx141, 2017

[2] Oscar Westesson, Mitchell Skinner, Ian Holmes, Visualizing next-generation sequencing data with JBrowse. *Brief Bioinform*, 14 (2): 172-177, 2013

[3] Yad Ghavi-Helm, Felix A. Klein, Tibor Pakozdi, Lucia Ciglar, Daan Noordermeer, Wolfgang Huber, Eileen E. M. Furlong, Enhancer loops appear stable during development and are associated with paused polymerase. *Nature*, 512(7512):96-100, 2014
# GeneSpy, a simple tool to explore genomic context

Pierre Simon GARCIA¹,², Frédéric JAUFFRIT¹,³, Christophe GRANGEASSE² and Céline BROCHIER-ARMANET¹

¹ LBBE Laboratoire de Biométrie et Biologie Evolutive, UCB Lyon 1 - Bât. Grégor Mendel
43 bd du 11 novembre 1918, 69622 VILLEURBANNE cedex, France

² MMSB Molecular Microbiology and Structural Biochemistry, Institut de Biologie et de
Chimie des Protéines 7 passage du Vercors, 69367 LYON7 cedex, France

³ R&D microbiology, Clinical Unit, bioMerieux, 376 Chemin de l'Orme,
69280 Marcy L'Etoile, France

Corresponding Author: pierre.garcia@etu.univ-lyon1.fr

Genomic context exploration is an important approach to provide informations in phylogenetic and genomic analysis: it brings elements to highlight evolutionary events such as duplications or horizontal transfers, functional links between genes and to assist in identification of orthologs. There are currently some available tools to visualize the genomic context of genes, such as MgcV [1], but they have some drawbacks. Indeed, most of these tools are based on non-local databases (generally made up of a set of complete genomes), do not allow dynamic navigation along genomes, provide non-publishable figures with a restricted selection of export formats (mainly pdf) and are generally slow to generate figures. In this context, we developed an application able to quickly generate visualizations of genomic context of genes, GeneSpy. This tool allows users to explore the organization of genomes and compare them through a user-friendly interface.

GeneSpy is written in python 2.7 and depends on Matplotlib and Tkinter libraries. GeneSpy uses gff files to build its local database. A download utility can be used to retrieve and prepare gff files for a set of assemblies from genbank or refseq. The list of target genes can be provided by the user via a tabular text file containing assembly identifiers of genomes and accession identifiers of proteins. The tool also supports BLAST reports from NCBI. Alternatively, users can search genes using keywords (fragment of name of strain and accession number, locus tag, name or annotation).

The context visualization in itself is heavily customizable. Many options are available such as window length, arrow size, type of coloration or verbosity of strains description. By default, the colors used for genes are procedurally generated using protein names and annotations for color consistency between contexts. Manual editing of color attribution is possible through a dedicated menu. The option to color only the genes of interest is also available.

The user can access informations relative to any displayed gene such as accession number, locus tag and predicted function. Moreover, it is possible to update the list of target genes by selecting any gene shown in the context figure. Thus, the user can simply navigate along the genome.

Figures can be exported in many formats such as png, jpg, tif,.svg, pdf and multiple pdf. An iTOL-specific export format is also available. It consists of a text file that can be imported in iTOL [2]. This format can be used to annotate a tree with relevant genomic context information, which can be very useful to study evolutionary dynamic of genomic regions. It is worth noting that database can be generated independently of GeneSpy, meaning that a pre-existing local databases can be reused with GeneSpy.

Thus, we provide a tool able to generate publishable figures of genomic context that is simple to use, fast, flexible, and adapted to any database. GeneSpy is available at https://lbbe.univ-lyon1.fr/GeneSpy/ and distributed under CeCILL licence.

This work was supported by grant from ARC1 Santé Rhône-Alpes Auvergne.

## References

[1] Overmars L, Kerkhoven R, Siezen RJ, Francke C. MGCv: the microbial genomic context viewer for comparative genome analysis. BMC Genomics; 14:209. 2013.

[2] Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res.;44(W1):W242-5. 2016.
# MCXpress: An R Package for functional interpretation of single cell RNA-Seq data using multivariate analysis.

Akira CORTAL¹, Antonio RAUSELL²

1,2Imagine InstituteClinical Bioinformatics Lab, 75015, Paris, France

Corresponding Author: akira.cortal@institutimagine.org
antonio.rausell@institutimagine.org

## Abstract

Single-cell RNA-sequencing allows unbiased transcriptome profiling of individual cells, enabling the analysis of genes expression at the cellular level. By uncovering cell heterogeneity in a given cell type, discovery of rare subpopulations of cells which could be responsible of the onset and progression of specific diseases are possible. Here we developed a comprehensible tool for the analysis of RNA-seq data using an approach based on Multiple Correspondence Analysis (MCA). MCA is a dimensionality reduction technique that allows the representation of both individuals (cells) and the variables (genes) within the same Euclidean space, thus allowing the simultaneous identification of subpopulations of cells and their gene signatures. MCA was coupled with gene set enrichment analysis; a powerful analytical method for characterizing differentially expressed genes/pathways within each groups of individuals. This singular combination allows a joint comparison of gene expression and pathway enrichment across all the groups. Combined with a shiny data visualization interface, the MCXpress R package can enhance the interpretation of both single and bulk RNA-sequencing data analysis. The tool development version is available on github at https://github.com/cbl-imagine/MCXpress (MIT License).

## References

[1] Shalek, A. K. et al. Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells, Nature 498, 236–240, 2013.

[2] Buettner F et al, Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells, Nature Biotechnology 33(2) 155, 2015.

[3] Rausell A et al, Protein interactions and ligand binding: From protein subfamilies to functional specificity, PNAS 107(5) 1995, 2010.
# Deciphering The Functional Effects of Genetic Variation With UniProt Annotations.

Benoit BELY¹, Andrew NIGHTINGALE¹ and Maria MARTIN¹

¹ Protein Function Development EMBL-EBI, Wellcome Trust Genome Campus, CB10 1SD, Hinxton, UK

Corresponding Author: benoit.bely@ebi.ac.uk

## 1 Introduction

Variants such as, substitution, frame shift and in frame deletions and insertions, in protein coding genes may cause deleterious effects manifested as debilitating inherited diseases or syndromes. Variants that cause a disease, do so by altering gene transcription or key parts of the proteins resulting in the disruption of a protein's native structure or its ability to execute functions, for example, to catalyze chemical reactions. The UniProt KnowledgeBase (UniProtKB) provides protein function information including annotation of protein altering variants, with any known functional effects and disease associations taken from scientific literature and large-scale public experimental data sets. Providing this rich knowledge at the genetic level is essential for discovering the potential deleteriousness of genetic variants.

## 2 Human proteome variation

UniProt has mapped protein annotations in the human proteome to the GRCh38 assembly of the human genome. Twenty-seven structural and functional annotations are currently provided including: enzyme active sites, modified residues, protein binding domains, protein variations etc. The complete set of human proteome sequences, including isoforms are also provided, for comparison to predicted transcripts. These mappings and related annotations have been made available as text BED files and BigBed files that are compatible with most genome browsers (http://www.uniprot.org/downloads). These files are also bundled into a public track hub that is available with the Ensembl and UCSC genome browsers through their public track hub registries or can be accessed directly from the genome track hubs registry (https://trackhubregistry.org). UniProt protein annotations can also be integrated directly into large scale genomic data programmatically via the Proteins API (http://www.ebi.ac.uk/proteins/api/doc/) [1].

## 3 Disease associated variants

To illustrate how the integration of protein function knowledge with high throughput genomic data provides unique opportunities for biomedical research. We examine some specific biological examples in disease related genes and proteins illustrating the utility of the combining protein and genome annotations for the functional interpretation of variants. Thus, showing that the UniProt genomic mappings, can help scientists to rapidly comprehend complex processes in biology.

## References

[1] Nightingale A., Antunes R., Alpi E., Bursteinas B., Gonzales L., Liu W., Luo J., Qi G., Turner E. and Martin, M. The Proteins API: accessing key integrated protein and genome information. *NAR Webservices Edition*, 2017.
Conférence invitée

John HUELSENBECK

<figure><img src="image_23.png" /></figure>

Department of Integrative Biology, UC Berkeley, USA

Bayesian inference in phylogeny for genome-scale data

Bayesian inference has permeated the field of phylogenetics. A major challenge in the field remains how to extend methods to genome-scale data. The temptation is to take short-cuts by applying fast methods that do not take full advantage of the information contained in the data. I describe several methods that may be applicable to genome-scale data. First, I describe new proposal mechanisms for better inferring large phylogenetic trees. Second, I discuss a class of models that can be used to address questions such as the identification of sites under the influence of natural selection.
# Extreme halophilic archaea derive from two distinct methanogen

## Class II ancestors

Monique Aouad¹, Najwa Taib¹, Anne Oudart¹, Michel Lecocq¹, Manolo Gouy¹ and Céline Brochier-Armanet¹

¹ Univ Lyon, Université Lyon 1, CNRS, UMR5558, Laboratoire de Biométrie et Biologie Évolutive, 43 bd du 11 novembre 1918, F-69622, Villeurbanne, France

Corresponding Author: monique.aouad@univ-lyon1.fr

**Abstract** The phylogenetic analysis of conserved core genes has disentangled most of the ancient relationships in Archaea. However, some groups remain debated, like the DPANN, a recently proposed deep-branching super-phyllum gathering various lineages of nanosized archaea with reduced genomes. Among these, the Nanohaloarchaea thrive in high-salt environments and require high-salt concentrations for growth. The discovery of Nanohaloarchaea in 2012 was significant because at the time extreme halophilic archaea were represented by a single lineage, the Halobacteria, a major class belonging to Euryarchaeota. The phylogenetic position of Nanohaloarchaea is highly debated, being alternatively proposed as the sister-lineage of Halobacteria or member of the DPANN super-phyllum. Pinpointing the phylogenetic position of extreme halophilic archaea is important to improve our knowledge of the deep evolutionary history of Archaea and decipher the underlying molecular adaptive processes and the evolutionary paths that allowed their emergence. Using comparative genomic approaches, we identified more than 250 protein markers carrying a reliable phylogenetic signal to address this issue. By combining strategies limiting the impact of biases on tree inferences, we showed that Nanohaloarchaea and Halobacteria represent two independent lineages that derived from methanogens Class II. This implies that adaption to very high salinity emerged twice independently in Archaea and that the grouping of Nanohaloarchaea within DPANN lineages is the consequence of a tree reconstruction artifact, which could challenge the existence of this group.

**Keywords** Stenosarchaea, compositional bias, long branch attraction, Slow-Fast method, rate signal

## Introduction

Recent advances in high-throughput sequencing technologies have revealed many new major uncultured environmental groups, most of them being known only through ribosomal RNA or genomic sequences [1-3]. This is for instance the case of *Nanohaloarchaea*, a group of extreme halophilic nanosized *archaea*, discovered recently in Lake Tyrell, Australia [4]. Some studies suggested that they represent the sister-lineage of *Halobacteria* [4,5], the other lineage of extreme halophilic archaea, while other analyses, based on different sets of markers, different methods and/or different taxonomic samplings, suggested instead that *Nanohaloarchaea* belong to the recently proposed DPANN super-phyllum. This deep-branching group encompasses diverse fast evolving, possibly nanosized archaea (e.g. *Diapherotrites*, *Parvarchaota*, *Micrarchaota*, *Aenigmarchaota*, *Nanoarchaota*, *Woosearchaota*, *Pacearchaota*) including *Nanohaloarchaea* [2,3,6]. Regarding *Halobacteria*, phylogenetic analyses of the RNA component of the small subunit of the ribosome and large supermatrices of conserved core genes have revealed a close relationship between *Halobacteria* and methanogens Class II, a group encompassing *Methanomicrobiales*, *Methanosarcinales* and *Methanocellales*, and indicated that *Halobacteria* could derive from a methanogenic ancestor [7]. However, the identity of the closest relative of *Halobacteria* remains debated. In fact, recently published phylogenies supported *Halobacteria* as the sister-lineage of all methanogens Class II [8-11], of *Methanomicrobiales* [12-14], or of *Methanocellales* [5,15].

Elucidating the precise position of *Halobacteria* and *Nanohaloarchaea* is particularly challenging because their proteomes harbor atypical amino acid compositions as a consequence of their extremophilic lifestyle. This can generate a compositional signal that may conflict with and dominate over the phylogenetic signal
[16], and lead to artifactual tree reconstructions where distant sequences with similar compositions are clus-
tered together [17,18]. Another source of bias could be linked to the fast evolutionary rate of nanohaloar-
chaecal and halobacterial proteomes highlighted by their very long branches in phylogenetic trees compared
to other archaeal lineages [4,5]. The phylogenetic position of fast-evolving species and long branches is par-
ticularly difficult to determine because differences in evolutionary rates among lineages can generate a rate
signal that may conflict the phylogenetic signal [16] and cause tree reconstruction artifacts such as the long
branch attraction (LBA) [19]. This well-known tree reconstruction artifact tends to group fast-evolving se-
quences/long branches and slow-evolving sequences/short branches in different parts of the trees when the
rate signal dominates over the phylogenetic signal [16]. Accordingly, we may wonder to what extent the con-
flicting positions observed for *Nanohaloarchaea* and *Halobacteria* are the consequence of tree reconstruc-
tion artifacts and if it is possible to overcome them.

To address this issue we performed an in-depth phylogenomic analysis designed to limit the impact of the non-phylogenetic signal on phylogenetic inferences. We showed that *Nanohaloarchaea* and *Halobacteria* branch robustly with *Methanocellales* and *Methanomicrobiales*, respectively, meaning that they derive from two distinct methanogen Class II ancestors. This implies also that adaptation to very high salinity occurred at least twice in archaea, and that the phenotypical similarities of *Nanohaloarchaea* and *Halobacteria* result from convergent evolutionary processes, possibly accompanied by horizontal gene transfers. Finally, our re-sults indicate also that the grouping of *Nanohaloarchaea* with other DPANN lineages is the consequence of a tree reconstruction artifact, challenging the existence of this candidate super-phyllum.

Results

The comparison of 155 proteomes from ANME-1, methanogens Class II, *Halobacteria* and their close rela-
tives: *Archaeoglobales* and *Diaforarchaea*, and from *Nanohaloarchaea* led to the delineation of 108,007
families, among which 258 presented a broad taxonomic distribution and no or very few evidences of hori-
zontal gene transfers (HGT) among these lineages and/or gene duplications. *Diaforarchaea* represent the first
diverging lineage, according to previous studies.

To test the impact of missing data on phylogenetic inference, different versions of these supermatrices were
built by gathering protein families present in more than 95% or 70% of the studied proteomes. The 68 pro-
tein markers present in the three *Nanohaloarchaea* and in more than 70% of the 155 studied proteomes (in-
cluding *Halobacteria*) were combined to build the FNANOHALO<sub>70</sub> supermatrix. The corresponding ML tree
was overall well resolved (Fig 1). In particular, the sistership between *Methanocellales* and *Methanosarci-
nales* was recovered (BV = 100%), indicating that adding *Nanohaloarchaea* and *Halobacteria* did not cause
major tree reconstruction artifacts. Regarding extreme halophiles, *Nanohaloarchaea* and *Halobacteria* did
not group together. In fact, *Halobacteria* clustered with *Methanomicrobiales* (BV = 99%, Fig 1), in agree-
ment with recent studies [12–14], whereas *Nanohaloarchaea* branched on the stem of *Diaforarchaea* (BV =
99%, Fig 1). This indicated that *Nanohaloarchaea* are not related to any methanogens Class II lineage,
ANME-I or *Archaeoglobales*, and is compatible with a deep-branching position of *Nanohaloarchaea* within
*Archaea*, as postulated by the DPANN hypothesis [2,3,6].
<figure><img src="image_24.png" /><figcaption>**Fig 1.** Maximum Likelihood phylogeny inferred with the FNANOHALO<sub>70</sub> supermatrix (18,309 positions, 155 sequences). The tree was inferred with PHYML 3.1 using the LG+I+G4+F model as proposed by IQ-TREE (BIC criteria). The scale bar corresponds to the average number of substitutions per site. Numbers at nodes correspond to bootstrap supports (100 replicates of the original dataset).</figcaption></figure>

Sequence compositional heterogeneity is an important source of bias in molecular phylogeny [16,20]. To test the impact of compositional biases on the inferred phylogenies, we removed from the supermatrices the positions most strongly responsible for the amino acid composition heterogeneity among sequences. The removal of positions with the highest compositional bias did not impact the phylogenetic position of extreme halophilic archaea. The clustering of *Halobacteria* with *Methanomicrobiales* and the branching of *Nanohaloarchaea* on the stem of *Diaforarchaea* were recovered again with high supports.

The evolutionary rate signal is one of the major causes of tree reconstruction artifacts such as the LBA [17,19]. This artifact is caused by multiple substitutions occurring at the same sites, a process which erases progressively the most ancient phylogenetic signal and results in the grouping of sequences according to their evolutionary rates (i.e. rate signal) in different parts of the inferred trees. To overcome this issue, dedicated methods have been developed.

Among them, the recoding of amino acids allows to hide substitutions among similar amino acids. To test the impact of the rate signal, we applied two different recoding schemes (dayhoff4 and dayhoff6) to all supermatrices in a Bayesian framework. The inferred BI trees confirmed the sister-ship between *Halobacteria* and *Methanomicrobiales*. Yet, surprisingly, *Nanohaloarchaea* branched with *Methanocellales* in seven out of the eight recoded supermatrices (Fig 2), suggesting that the branching of *Nanohaloarchaea* on the stem of *Diaforarchaea* could result from a tree reconstruction artefact due to the rate signal.
Tree scale: 0.1-

<figure><img src="image_25.png" /></figure>

**Fig 2.** Bayesian phylogeny inferred with the FNANOHALO<sub>70</sub> supermatrix recoded according to the Dayhoff4 scheme (18,309 positions, 155 sequences). The tree was inferred with Phylobayes using the CAT + GTR + G4 evolutionary model. The scale bar corresponds to the average number of substitutions per site. Numbers at nodes correspond to posterior probabilities.

To test this hypothesis, we used another approach to limit tree reconstruction artifacts resulting from the rate signal. This approach, called the Slow-Fast method (S-F), consists in the progressive removal of the fastest evolving sites from multiple alignments [21]. The S-F method was shown to be very efficient to reduce tree reconstruction artifacts because the fastest evolving sites are the most susceptible to be impacted by multiple substitutions [17,22]. This approach allows monitoring the support associated to a given branch of a tree throughout the removal process and thus to determine if the corresponding relationship reflects the phylogenetic or the rate signal contained in sequences [17]. To avoid biases in the estimation of evolutionary rates at each position due to missing data and/or unbalanced taxonomic sampling among lineages, the S-F method was applied to the supermatrices built with markers present in at least 95% of the studied taxa and by keeping only three to seven representative sequences for each archaeal lineage.

Removal of the fastest-evolving sites did not impact the phylogenetic position of *Halobacteria*. In fact, the grouping of *Halobacteria* with *Methanomicrobiales* was strongly supported in ML and BI trees inferred with the S-F supermatrices. This suggested that this relationship was not the consequence of the rate signal. In sharp contrast, the S-F method showed that a robust grouping of *Nanohaloarchaea* with *Methanocellales* was observed in ML and BI trees when the fastest-evolving sites were removed (Fig 3).
A-

<figure><img src="image_26.png" /></figure>

B-

<figure><img src="image_27.png" /></figure>

**Fig 3.** Effect of removal of fast-evolving positions on the phylogenetic position of *Halobacteria* and *Nanohaloarchaea* in the FNANOHALO<sub>95</sub> S-F matrices. The x axes correspond to the fraction of sites (in %) kept in the supermatrices. The removal of fastest-evolving sites proceeds from right to left. The y axes indicate the support for the grouping of *Halobacteria* with *Methanomicrobiales* (in pink), the grouping of *Nanohaloarchaea* with *Halobacteria* (in red), and the grouping of *Nanohaloarchaea* with *Methanocellales* (in blue) in ML (A) and BI trees (B) resulting from the S-F analysis. Notice that the *Halobacteria* + *Nanohaloarchaea* grouping supported by FNANOHALO<sub>95</sub> disappears when fast-evolving site are removed.
Discussion

By focusing our analysis on the part of the euryarchaeal tree that contains *Halobacteria* and methanogens Class II, we were able to assemble larger datasets of conserved markers and to use more intense taxonomic samplings compared to previous studies. More precisely, we identified 258 conserved protein families widely distributed in the 155 taxa, among which 68 were present in the three representatives of *Nanohaloarchaea* available at the time. This led us to the construction of large supermatrices containing thousands of amino acid positions. By using various methods allowing to decouple the different types of signal contained in protein sequences, we showed that the compositional signal did not significantly impact the phylogenetic positions of *Nanohaloarchaea* and *Halobacteria*, while the rate signal had a major impact on the phylogenetic position of *Nanohaloarchaea*. In fact, two independent methods allowing to reduce the impact of multiple substitutions on phylogenetic inferences, the removal of the fastest evolving sites and the recoding of amino acids, provided consistent results supporting the grouping of *Halobacteria* with *Methanomicrobiales* and of *Nanohaloarchaea* with *Methanocellales*. The robust and recurrent grouping of *Halobacteria* and methanogens Class II in many studies, suggested that they could represent a new super-class, that we propose to call *Stenosarchaea* (from the Greek stenós, meaning close/joint). Regarding *Nanohaloarchaea*, our analyses strongly suggested that they are part of the *Stenosarchaea*, and more precisely, that they represent the sister group to *Methanocellales*.

The grouping of *Nanohaloarchaea* with *Methanocellales* and that of *Halobacteria* with *Methanomicrobiales*
within *Stenosarchaea* has major implications and opens new perspectives. First, it implies that both lineages
derive from two distinct but related methanogen ancestors, which is in accordance with the fact that a few ar-
chaeal lineages, all belonging to *Euryarchaeota*, can survive at high salt concentrations [23]. These
halophilic or salt-tolerant archaea are anaerobic methanogens living in hypersaline sediments. Most are
methylotroph and belong to methanogens Class II and more precisely to *Methanosarcinaceae* (order
*Methanosarcinales*) or *Methanocalculaceae* (order *Methanomicrobiales*) [23]. It also implies that adaptation
to extreme high salt concentrations occurred at least twice independently during the evolution of *Archaea*.
Thus, the phenotypic properties shared by *Nanohaloarchaea* and *Halobacteria* should be interpreted as the
consequence of a convergent evolution that could have been facilitated by HGT. In that context, it would be
interesting to reevaluate the evolutionary history of these lineages, and the role played by HGT in the emer-
gence of *Halobacteria* and *Nanohaloarchaea* from methanogenic ancestors.

Finally, the robust grouping of *Nanohaloarchaea* with *Methanocellales* rules out alternative hypotheses for the branching of *Nanohaloarchaea*, and in particular a branching outside of *Stenosarchaea*, as expected if *Nanohaloarchaea* were part of the candidate DPANN superphylum. This challenges the existence of this group as it is currently described and questions to what extent similar artifacts could also impact the position of the other lineages composing this super-phylum, which are all fast evolving. Testing this hypothesis would require accurate and separate analyses, each focused on one lineage.

References

1. Schleper C, Jurgens G, Jonuscheit M. Genomic studies of uncultivated archaea. Nat. Rev. Microbiol. [Internet]. 2005 [cited 2016 Jun 20];3:479–88. Available from: http://www.ncbi.nlm.nih.gov/pubmed/15931166

2. Rinke C, Schwientek P, Sczyrba A, Ivanova NN, Anderson IJ, Cheng J-F, et al. Insights into the phylogeny and cod-
ing potential of microbial dark matter. Nature [Internet]. 2013 [cited 2016 Jun 20];499:431–7. Available from:
http://www.ncbi.nlm.nih.gov/pubmed/23851394

3. Castelle CJ, Wrighton KC, Thomas BC, Hug LA, Brown CT, Wilkins MJ, et al. Genomic expansion of domain ar-
chaea highlights roles for organisms from new phyla in anaerobic carbon cycling. Curr. Biol. [Internet]. 2015 [cited
2016 Jun 20];25:690–701. Available from: http://www.ncbi.nlm.nih.gov/pubmed/25702576

4. Narasingarao P, Podell S, Ugalde JA, Brochier-Armanet C, Emerson JB, Brocks JJ, et al. De novo metagenomic as-
sembly reveals abundant novel major lineage of Archaea in hypersaline microbial communities. ISME J. [Internet].
2012 [cited 2016 Jun 20];6:81–93. Available from: http://www.ncbi.nlm.nih.gov/pubmed/21716304

5. Petitjean C, Deschamps P, López-Garciá P, Moreira D. Rooting the domain archaea by phylogenomic analysis sup-
ports the foundation of the new kingdom Proteoarchaeota. Genome Biol. Evol. 2014;7:191–204.
6. Williams TA, Embley TM. Archaeal &quot;dark matter&quot; and the origin of eukaryotes. Genome Biol. Evol. [Internet]. 2014 [cited 2016 Jun 20];6:474–81. Available from: http://www.ncbi.nlm.nih.gov/pubmed/24532674

7. Forterre P, Brochier C, Philippe H. Evolution of the Archaea. Theor. Popul. Biol. [Internet]. 2002 [cited 2016 Jun 20];61:409–22. Available from: http://www.ncbi.nlm.nih.gov/pubmed/12167361

8. Gao B, Gupta RS, Woese C, Kandler O, Wheelis M, Pace N, et al. Phylogenomic analysis of proteins that are distinctive of Archaea and its main subgroups and the origin of methanogenesis. BMC Genomics [Internet]. BioMed Central; 2007 [cited 2016 Jun 20];8:86. Available from: http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-8-86

9. Yutin N, Puigbò P, Koonin E V., Wolf YI, Ramakrishnan V, Ban N, et al. Phylogenomics of Prokaryotic Ribosomal Proteins. Lespinet O, editor. PLoS One [Internet]. Public Library of Science; 2012 [cited 2016 Jun 20];7:e36972. Available from: http://dx.plos.org/10.1371/journal.pone.0036972

10. Wolf YI, Makarova KS, Yutin N, Koonin E V. Updated clusters of orthologous genes for Archaea: a complex ancestor of the Archaea and the byways of horizontal gene transfer. Biol. Direct [Internet]. 2012 [cited 2016 Jun 20];7:46. Available from: http://www.ncbi.nlm.nih.gov/pubmed/23241446

11. Nelson-Sathi S, Dagan T, Landan G, Janssen A, Steel M, McInerney JO, et al. Acquisition of 1,000 eubacterial genes physiologically transformed a methanogen at the origin of Haloarchaea. Proc. Natl. Acad. Sci. U. S. A. [Internet]. 2012;109:20537–42. Available from: http://www.pnas.org/content/109/50/20537 abstract

12. Brochier-Armanet C, Forterre P, Gribaldo S. Phylogeny and evolution of the Archaea: One hundred genomes later. Curr. Opin. Microbiol. 2011. p. 274–81.

13. Raymann K, Brochier-Armanet C, Gribaldo S. The two-domain tree of life is linked to a new root for the Archaea. Proc. Natl. Acad. Sci. U. S. A. [Internet]. 2015 [cited 2016 Jun 20];112:6670–5. Available from: http://www.ncbi.nlm.nih.gov/pubmed/25964353

14. Petitjean C, Deschamps P, López-García P, Moreira D, Brochier-Armanet C. Extending the conserved phylogenetic core of archaea disentangles the evolution of the third domain of life. Mol. Biol. Evol. 2015;32:1242–54.

15. Becker EA, Seitzer PM, Tritt A, Larsen D, Krusor M, Yao AI, et al. Phylogenetically Driven Sequencing of Extremely Halophilic Archaea Reveals Strategies for Static and Dynamic Osmo-response. Whitaker RJ, editor. PLoS Genet. [Internet]. Public Library of Science; 2014 [cited 2016 Jun 20];10:e1004784. Available from: http://dx.plos.org/10.1371/journal.pgen.1004784

16. Jeffroy O, Brinkmann H, Delsuc F, Philippe H. Phylogenomics: the beginning of incongruence? Trends Genet. 2006;22:225–31.

17. Delsuc F, Brinkmann H, Philippe H. Phylogenomics and the reconstruction of the tree of life. Nat. Rev. Genet. [Internet]. Nature Publishing Group; 2005 [cited 2016 Jun 20];6:361–75. Available from: http://www.nature.com/doifinder/10.1038/nrg1603

18. Woese CR, Achenbach L, Rouviere P, Mandelco L. Archaeal phylogeny: reexamination of the phylogenetic position of Archaeoglobus fulgidus in light of certain composition-induced artifacts. Syst. Appl. Microbiol. [Internet]. 1991 [cited 2016 Jun 20];14:364–71. Available from: http://www.ncbi.nlm.nih.gov/pubmed/11540072

19. Felsenstein J. Cases in which Parsimony or Compatibility Methods Will Be Positively Misleading. Syst. Zool. [Internet]. 1978 [cited 2016 Jun 20];27:401. Available from: http://www.jstor.org/stable/2412923?origin=crossref

20. Gribaldo S, Philippe H. Ancient Phylogenetic Relationships. Theor. Popul. Biol. [Internet]. 2002;61:391–408. Available from: http://www.scopus.com/inward/record.url?eid=2-s2.0-0036593395&partnerID=40&md5=c809464d9-ca50344aef3f1c0ecd9c871

21. Brinkmann H, Philippe H. Archaea sister group of Bacteria? Indications from tree reconstruction artifacts in ancient phylogenies. Mol. Biol. Evol. [Internet]. 1999 [cited 2016 Jun 20];16:817–25. Available from: http://www.ncbi.nlm.nih.gov/pubmed/10368959
22. Philippe H. Opinion: long branch attraction and protist phylogeny. Protist [Internet]. 2000 [cited 2016 Jun 20];151:307–16. Available from: http://www.ncbi.nlm.nih.gov/pubmed/11212891

23. Oren A. Taxonomy of halophilic Archaea: current status and future challenges. Extremophiles [Internet]. 2014 [cited 2016 Jun 20];18:825–34. Available from: http://www.ncbi.nlm.nih.gov/pubmed/25102811
# Origin and evolution of multiple haem copper oxidases in *Archaea*

Anne Oudart<fnref n="1" />, Simonetta Gribaldo<fnref n="2" /> and Céline Brochier-Armanet<fnref n="1" />

Corresponding Author: celine.brochier-armanet@univ-lyon1.fr

## Abstract

Understanding the origin and evolution of dioxygen reductases, the terminal electron acceptors of aerobic respiratory chains, can provide precious clues for the emergence of this energy process that is still debated. The haem-copper oxidases superfamily (HCO), belonging to the dioxygen reductases, contains three families of dioxygen reductases named A, B, and C according to a classification based on sequence similarity and phylogenetic analysis of the homologous catalytic subunits [1]. A phylogenomic study showed that these enzymes have very different evolutionary histories [2] with an ancient dioxygen reductase (A-HCO) present prior to the divergence of major present-day bacterial and archaeal phyla, thus before the emergence of oxygenic photosynthesis. However, this result is in contradiction with those proposed by Ducluzeau et al. 2014 [3] about the structure of the A-HCO suggesting that this dioxygen reductase would be the most recent. So the question about the origin of the haem-copper oxidases is still unresolved, and a new analysis is required. The haem-copper oxidases superfamily of the ancestor of Archaea probably used dioxygen [4]. Nevertheless, the previous conclusions are old and the available data was very limited (73 archeal complete genomes in 2009). Today, with more available data in the public database (252 archeal complete genomes) it is interesting to reassess the questions about the origin and evolution of the haem-copper superfamily for Archaea. We will present our results about subunits of the haem-copper superfamily in Archaea.

**Keywords** Aerobic respiration, LUCA, Phylogeny, Prokaryotes, Early earth

## Introduction

Dioxygen reductases (O₂-Red) represent the terminal electron-transfer enzymes of aerobic respiratory chains in the three Domains of Life (*Archaea*, *Bacteria* and *Eucarya*). These membrane-bound enzymes catalyze the reduction of dioxygen (O₂) to water by using electrons provided by either a quinol derivate or a cytochrome c. Because O₂-Red are key enzymes of aerobic respiration, understanding their origin and evolution could provide precious clues for the emergence of this key metabolism, the rise of O₂ on early Earth and its impact on the evolution of ancient microbial communities. The O₂-Red encompasses HCO divided into three subfamilies named A, B, and C according to a recently proposed classification based on sequence similarity and phylogenetic analyses of their catalytic subunits [1,5]. The O₂-Red include also the Nitric Oxide Reductases (NOR) which reduce nitric oxide (NO) to nitrous oxide (N₂O) [2].

The catalytic subunit of A-HCO is found in many bacterial lineages [2]. In archaea, it was detected in *Crenarchaeota* and *Thaumarchaeota* (two phyla belonging to the TACK super-phylum), and in *Diaforarchaea* and *Halobacteria*, two euryarchaeotal orders [2]. The catalytic subunit of B-HCO was reported in a few bacterial lineages, in *Crenarchaeota* and in *Halobacteria* [2]. The catalytic subunit of C-HCO is mainly present in *Proteobacteria* [2]. Finally, NOR catalytic subunit was mainly found in *Proteobacteria*, in some other *Bacteria*, and in a few *Archaea* [2].

The evolutionary history of HCO O₂-Red is highly debated, and three different scenarios have been proposed. First, it was proposed that the four HCO O₂-Red were present in the Last Universal Common Ances-

<footnote n="1">¹Univ Lyon, Université Lyon 1, CNRS, UMR5558, Laboratoire de Biométrie et Biologie Évolutive, 43 bd du 11 novembre 1918, F-69622, Villeurbanne, France.</footnote>

<footnote n="2">²Unité de Biologie Moléculaire du Gene chez les Extremophiles, Département de Microbiologie, Institut Pasteur, 75015 Paris, France</footnote>
tor (LUCA) [6]. More recently, it was proposed that only A-HCO was present in the ancestor of Bacteria, and possibly in LUCA, while the three other HCO O₂-Red subfamilies emerged later, in *Archaea* (B-HCO) and in *Proteobacteria* (NOR and C-HCO) [2,4]. According to these two scenarios, the LUCA and/or the ancestor of Bacteria could have been able to reduce O₂. Because these two organisms are ancestors of *Cyanobacteria*, the scenarios implied that some sources of O₂ should have been present before the emergence of oxygenic photosynthesis. In contrast, the third scenario proposed the NOR to be present in the LUCA, while HCO O₂-Red emerged during the diversification of *Bacteria* after the emergence of *Cyanobacteria* and oxygenic photosynthesis, and spread through HGT in archaeal and bacterial lineages [1,3,7]. According to this scenario, HCO O₂-Red acquired their capacity to reduce O₂ independently from distinct NOR ancestors after the emergence of *Cyanobacteria* and of oxygenic photosynthesis.

At the time, the data available were mainly limited to *Bacteria* because the number of genomes (complete or not) of *Archaea* was very restricted compared to *Bacteria*. This lack of data in *Archaea* limited considerably our capacity to understand the ancient evolutionary history of HCO O₂-Red thus of one of the most important metabolism on earth. Fortunately, the situation is changing rapidly, and now hundred genomes of *Archaea* are available in public databases, thanks to ambitious projects aiming at sequencing the cultured and the uncultured "dark matter" of the microbial diversity [8,9]. This rainfall of genomic data offers the possibility to investigate the evolutionary history of HCO O₂-Red in *Archaea*.

In *Archaea*, HCO O₂-Red have been extensively studied in *Sulfolobales*. In these *Crenarchaeota*, A-HCO and B-HCO, the catalytic subunits of O₂-Red, are associated with several proteins to form various different protein complexes (Figure 1).

<table><caption>**Fig 1.** Haem copper oxidases complexes in *Sulfolobales*.</caption><thead><tr><th rowspan="2">Subunits</th><th colspan="2">A-HCO</th><th colspan="2">B-HCO</th><th rowspan="2">Function</th></tr><tr><th>Sox</th><th>Sox</th><th>Dox</th><th>Fox</th></tr></thead><tbody><tr><td>Subunit I</td><td>SoxM</td><td>SoxB</td><td>DoxB</td><td>FoxA/FoxA</td><td>Reduce dioxygen</td></tr><tr><td>Subunit II</td><td>SoxH</td><td>SoxA</td><td>DoxC</td><td>FoxB</td><td>Transfer electron to Subunit I</td></tr><tr><td>Subunit III</td><td>Can be fused with SoxM</td><td></td><td></td><td></td><td>Assist the final step of the folding of the subunit I</td></tr><tr><td>Cytochrome b</td><td>SoxG</td><td>SoxC</td><td></td><td></td><td>Transfer electron to Rieske protein</td></tr><tr><td>Cytochrome b558/566</td><td></td><td></td><td></td><td>FoxCD</td><td>Transfer electron to multi blue copper oxidase</td></tr><tr><td>Rieske protein</td><td>SoxF</td><td>SoxL</td><td></td><td></td><td>Transfer electron to Sulfocyanin</td></tr><tr><td>Sulfocyanin</td><td>SoxE</td><td></td><td></td><td></td><td>Transfer electron to Subunit II</td></tr><tr><td>Protein of unknown function</td><td>Soxl</td><td></td><td></td><td></td><td>Unknown</td></tr><tr><td>Protein of unknown function</td><td></td><td>SoxD</td><td></td><td></td><td>Unknown</td></tr><tr><td>Protein of unknown function</td><td></td><td>SoxD'</td><td></td><td></td><td>Unknown</td></tr><tr><td>Protein of unknown function</td><td></td><td></td><td>DoxE</td><td></td><td>Unknown</td></tr></tbody></table>

The four complexes described in *Archaea* carried homologous subunits, suggesting that these respiratory protein complexes may derive from an ancestral protein complex.

## Material and methods

### Dataset assembly

The 2,610 bacterial complete proteomes and 391 archaeal complete and incomplete proteomes available at the NCBI (http://www.ncbi.nih.gov) and the JGI (http://jgi.doe.gov/) were retrieved and gathered in a local database. These corresponded to 1279 bacterial and 269 archaeal species.

Identification of the homologues of HCO subunits was performed by combining iterative sequence similarity searches and HMM profiles surveys using the BLASTP v2.2.26 [10] and the HMMER v3.1b1 programs [11] respectively.

### Classification of HCO catalytic subunit

HCO catalytic subunits were classified according to the classifier tool available at (http://www.evocell.org/hco/) [12].

### Phylogenetic analyses

For phylogenetic analyses, multiple alignments were built with MAFFT v7.123b (option L-INS-i) [13] with the BLOSUM30, BLOSUM45 and BLOSUM62 matrix. The accuracy of the resulting alignments was verified with SEAVIEW 4.5.4 [14] and compared with NorMD v1.2 [15]. The best alignments were
trimmed with BMGE v1.1 [16] with the best BLOSUM matrix identified with NorMD and used for phyloge-
netic analyses. To avoid taxonomic redundancy, only one representative strain per bacterial genus was kept
for phylogenetic reconstruction. The best fitting evolutionary models were identified by IQ-TREE (Nguyen
et al. 2015) according to the BIC criterion [17]. Maximum Likelihood (ML) trees were built using PHYML
v3.1 [18] with the NNI+SPR topology search option. The robustness of the resulting trees was estimated
with a non-parametric bootstrap procedure (100 replicates of the original datasets).

The trees were visualized with iToL [19,20].

Analysis of genomic contexts

Genomic contexts were generated with GeneSpy developed by Pierre Garcia (personal communication).

Construction of a reference phylogeny of archaeal species

A reference phylogeny of *Archaea* was inferred using ribosomal proteins. The 70 ribosomal protein fami-
lies available for 209 complete archaeal proteomes in RiboDB ([https://ribodb.univ-lyon1.fr/ribodb/ribodb-in.cgi](https://ribodb.univ-lyon1.fr/ribodb/ribodb-in.cgi)) [21] were retrieved and used as seed to query our local database. 39 on 391 available proteomes con-
tained less than 35 ribosomal protein families. To limit biases due to missing data, the corresponding strains
were not included in the phylogenetic analysis.

The 58 ribosomal proteins present in at least 70% of the 352 remaining proteomes were kept for phyloge-
netic analyses. For each ribosomal protein family a multiple alignment was built and trimmed as described
above. The resulting datasets were combined together to build a large supermatrix containing 352 sequences,
and 6,257 conserved amino acid positions. This supermatrix was used to inferred a ML tree with PHYML as
described above.

Results

At least one HCO O₂-Red catalytic subunit was detected in most bacterial and archaeal lineages. A-HCO and its accessory subunits homologues were found in many *Archaea* and *Bacteria*. In contrast, B-HCO and its functional partners were mainly present in *Archaea*, while C-HCO and NOR were almost all present in *Bacteria*. For each subunit of the four O₂-Red families, we inferred the ancestral presence of the subunit in a taxon when the gene was widely distributed in this taxon and when we could find the monophyly of the group.

NOR catalytic subunit presents a scarce taxonomic distribution, being rarely present in more than 50% of proteomes of bacterial or archaeal order. NOR sequences from different taxa appeared intermixed on the tree, indicating that the relationships among NOR sequences were inconsistent with the currently recognized prokaryotic systematics. This combined to the scarce taxonomic distribution of NOR, indicated that their evolutionary history was heavily impacted by HGT, as proposed previously [2]. In fact, based on phyloge- netic and taxonomic criteria, the ancestral presence of NOR can not be inferred in any lineage to the excep- tion of *Halobacteria*. Thus, it was very unlikely that this enzyme could have been present in LUCA as pro- posed few years ago [3,7], or even in the ancestor of Bacteria or Archaea.

C-HCO catalytic subunit was widely distributed in the *Bacteroidetes/Chlorobi* group and in *Proteobacte-
ria* that formed a well supported group, while they were scarce in other lineages. According to our data, the
possible ancestral presence of C-HCO in the ancestor of *Chlorobi* and of *Proteobacteria*, two major bacterial
lineages that didn't occupy a basal position in the phylogeny of *Bacteria* was not sufficient to infer its pres-
ence in the ancestor of Bacteria, as proposed recently [3]. The presence of these enzymes in strict anaerobic
archaea was puzzling.

Most of the B-HCO homologues are found in *Archaea*, and more precisely in *Geoarchaeota*, *Sulfolobales* and *Halobacteria*. Regarding bacterial sequences, the phylogeny of B-HCO revealed a complex evolutionary history dominated by HGT between *Bacteria* and *Archaea* and among *Bacteria*. In particular, according to phylogenetic and taxonomic criteria it was not possible to trace back B-HCO in any bacterial lineage. In contrast, its ancestral presence in three major archaeal lineages suggested that it originated in this domain of Life, in agreement with previous studies [2,3].

According to taxonomic distribution and phylogenetic criteria, the ancestral presence of A-HCO could be
inferred in *Cyanobacteria*, the *Bacteroidetes/Chlorobi* group, the *Fibrobacteres/Acidobacteria* group, *Al-
phaproteobacteria*, *Betaproteobacteria* and *Planctomycetes*. Other bacterial A-HCO appeared intermixed
suggesting that they spread through HGT in present-day bacterial lineages. Regarding Archaea, A-HCO ho-
mologues were widespread in several archaeal phyla, classes, and orders and formed monophyletic groups
on the ML phylogeny. This was for instance the case of *Thaumarchaeota, Halobacteria, Geoarchaeota,
Thermoproteales, Sulfolobales*, uncultured marine group II and Aigarchaeota, suggesting that they could be
ancestral in these groups. The grouping of sequences from *Aigarchaeota, Geoarchaeota, Thermoproteales*,
and *Sulfolobales*, close to a large *Thaumarchaeota* cluster, could indicate the ancestral presence of this en-
zyme in the ancestor of the recently proposed TACK superphylum. In contrast, the situation was less clear in
the case of *Diaforarchaea*, because they formed a paraphyletic group that also included archaeal sequences
from other taxa.

Beside catalytic subunits (subunit I), three additional subunits are shared by the A-HCO and B-HCO: the
subunit II (SoxH, SoxA and DoxC), the cytochrome *b* (SoxG and SoxC) and the Rieske protein (SoxF and
SoxL). The catalytic subunit and the subunit II co-occurred in almost all proteomes. Surprisingly, in 45% of
*Halobacteria* strains, a copy of the subunit II was not co-localized with the catalytic subunit whereas an an-
other copy did. At the exception of these Halobacteria, we can clearly distinguish the subunit II of the A-
HCO and of the B-HCO in the phylogeny. We can highlight a topology similar to the one obtained with the
catalytic subunit. Cytochrome *b* was not found in 30,7% of all proteomes of *Archaea* and complete genomes
of *Bacteria* that contain a catalytic subunit. In the B-HCO, only the Sox of *Sulfolobales* had a copy of cy-
tochrome *b* in the same genomic context as the catalytic subunit. The cytochrome *b* of A-HCO are in the
same genomic context as the catalytic subunit only in *Crenarchaeota*. Rieske proteins were not found in
42,4% of all proteomes of Archaea and complete genomes of Bacteria that contain a catalytic subunit. The
topology of Rieske protein is very similar to the topology of cytochrome *b*. It is not surprising because the
two proteins are almost always in the same genomic context. Like the cytochrome *b*, Rieske protein is not in
the same genomic context as the catalytic subunit at the exception of *Crenarchaeota*. The subunit III, subunit
of A-HCO, was not found in 14,3% of all proteomes of *Archaea* and complete genomes of *Bacteria* that con-
tain a catalytic subunit. The subunit III, in *Archaea*, is always fused with the catalytic subunit except for
some *Halobacteria*. The subunit III present a similar topology to the A-HCO catalytic subunit at the excep-
tion that *Thaumarchaeota* do not have this subunit.

Regarding other subunits, the sulfocyanin (SoxE) is associated to A-HCO. The subunit is almost found in
Archaea. More precisely, Sulfolobales sequences formed three monophyletic clusters, each of them could be
traced back in the ancestor of this archaeal order. In contrast, the presence of SoxE in other archaeal lineages
(e.g. *Thermoplasmatales, Aigarchaeota, Thermoproteales*) likely resulted from secondary HGT. We can
highlight a curiosity about this subunit because for some strains, SoxE is in the same genomic context as the
B-HCO catalytic subunit. SoxI was very likely present in the ancestors of *Geoarchaeota* and *Sulfolobales*
and maybe also present in the ancestor of *Thermoproteales* and *Aigarchaeota*. This subunit is in the same ge-
nomic context as the A-HCO catalytic subunit. Interestingly, a similar cluster existed in the phylogeny of A-
HCO, indicating that the presence of SoxI was co-opted specifically in the ancestor of this subgroup of A-
HCO and preserved along its subsequent diversification. SoxD and DoxE associated to the Sox and Dox
complexes (B-HCO), respectively, were found exclusively and in almost all *Sulfolobales*, suggesting they
emerged in the stem of this order. Finally, the SoxD' subunit (B-HCO) and the Fox complex were restricted
to a subgroup of *Sulfolobales* strains, indicating that they emerged late during the diversification of this ar-
chaeal order. Thus, they could not be traced back to the ancestor of this archaeal order and this implied that
the corresponding proteins were much younger than the other subunits.

Discussion

Studying ancient metabolisms opens windows on primitive geosphere and biosphere, and a key to deci-
pher the major transitions in the evolution of Earth [22]. Among them, the origin and rise of O₂ represent one
of the most debated questions [23], specifically regarding the divergent conclusions drawn from biological
and geological observations [23,24]. Present-day high level of O₂ in the atmosphere resulted of the oxygenic
photosynthesis performed by *Cyanobacteria* [22,25–27] (including plant chloroplast that are of cyanobacte-
rial origin), a process that started around 2.45-2.33 billion years ago [23,27,28]. Recent studies suggested
that the tectonic activity around 2.7 billion years ago could have also contributed to release significant
amounts of O₂ in the atmosphere (Lee et al. 2016; Sukumaran 2000; Jelen et al. 2016). The presence of small
amount of O₂ of abiotic origin in primitive oceans and atmosphere prior 2.7 billion years remains debated
[29–32]. In this context, studying the origin and evolution of HCO, the terminal electron acceptor of aerobic
respiratory chains, can provide precious clues for the emergence of this energetic process and the availability of $O_2$.

The survey of thousands of bacterial and archaean proteomes encompassing all major lineages allowed inferring the presence of A-HCO in the ancestor of *Cyanobacteria*, *Bacteroidetes*, *Chlorobi*, *Fibrobacteres*, *Acidobacteria*, *Alphaproteobacteria*, *Betaproteobacteria*, *Planctomycetes*, and likely in Gammaproteobacteria. This confirmed previous results suggesting that A-HCO were ancient in *Bacteria* and may have predated the emergence of *Cyanobacteria* and thus of oxygenic photosynthesis [2]. In contrast and contrarily to previous proposals [7,33,34], we didn't find any evidence suggesting that B-HCO, C-HCO and NOR could be ancient in *Bacteria*.

Regarding *Archaea*, our data revealed a very dynamic evolution. First, we confirmed that A-HCO catalytic subunit and some of the accessory subunits were ancient in this Domain of Life [2,4]. More precisely, we showed that NOR could have been present in the ancestor of *Halobacteria*, B-HCO in the ancestors of *Sulfolobales* and *Halobacteria*, and A-HCO in the ancestors of *Thaumarchaeota*, *Crenarchaeota*, *Aigarchaeota*, uncultured marine group II/III and *Halobacteria*. Worth noticing, A-HCO were not restricted to mesophilic archaea, represented here by *Halobacteria*, uncultured marine group II/III and some *Thaumarchaeota*. This contradicted the recent proposal that the acquisition of the corresponding gene from bacterial donors could have been a crucial step towards the colonization of low-temperature environments by thermophilic *Archaea* [35]. On the contrary, A-HCO and B-HCO were widespread in thermophilic/hyperthermophilic archaean lineages and we showed that B-HCO likely originated in *Archaea* and were subsequently transferred to a few *Bacteria*.

According to the rooting of the archaean tree in-between the TACK superphylum and *Euryarchaeota*, it was proposed previously that A-HCO could be ancestral in *Archaea* [2]. Alternatively, one could propose that A-HCO and B-HCO were acquired independently and secondarily by HGT after the diversification of *Archaea* and the emergence of the main archaean orders in response to the rise of $O_2$ level linked to the photosynthetic activity of *Cyanobacteria*. While, this conclusion was not contradict by our data, a new root splinting the Archaea into two clusters was recently proposed [36]. Cluster 1 corresponded to TACK, and two major euryarchaeotal lineages: the *Methanomada* (i.e. *Methanococcales*, *Methanobacteriales*, *Methanopyrales*) and the *Thermococcales*, while Cluster 2 encompassed *Diaforarchaea*, *Archaeoglobales*, *ANME-1*, *Halobacteria*, *Methanocellales*, *Methanosarcinales*, *Methanomicrobiales*. The position of this root opened the possibility that of a methanogenic and anaerobic ancestor for *Archaea* [36]. If confirmed, this would imply that the absence of any HGT of C-HCO from bacterial donor to archaea (excepted in two anaerobic methanosarcinales) was puzzling. However, independently of the root in *Archaea* our data indicate an ancient presence of A-HCO in *Archaea*.

In our study, we also carried out a thoughtful analysis of the other subunits included in the HCO complex in Archaea. Some subunits escorting the catalytic subunit, Cytochrome B, Rieske protein and the subunit II present the same evolutionary history than the catalytic subunit in Archaea and form the conserved core of A-HCO and B-HCO with the catalytic subunit. The subunit III is specific to the A-HCO complex and can be a part of the conserved core which has been lost by the ancestor of *Thaumarchaeota* and *Thermoplasmatales*. Interestingly, the other subunits have more recent origin and are a crenarchaeal innovation. Sulfocyanin (SoxE) of the A-HCO is particularly interesting because some copies are in the same genomic context as the B-HCO catalytic subunit. This result supports the hypothesis that the B-HCO could come from a duplication of the A-HCO in *Crenarchaeota* and explains that some specific subunits of A-HCO could work with B-HCO subunits. According to our results, it seems that the ancestor of TACK might have possessed two operons (A-HCO), one composed of the catalytic subunit, the subunit II and the subunit III and the other composed of the Rieske protein and the cytochrome b. During evolution, the subunit III was lost by the ancestor of *Thaumarchaeota* and *Thermoplasmatales* and the A-HCO complex has been duplicated in the ancestor of *Crenarchaeota* to give the B-HCO complex (with the lost of the subunit III) where the two operons get together. It is still unclear why Archaea and particularly Crenarchaeota innovate new subunits. Moreover, HCO was only studied in a few archaean models and there is a lot of innovations in Archaea with accumulation of new subunits, so we are wondering if we have discovered all the proteins imply in the aerobic respiration in Archaea.

An other interesting point is the presence of $O_2$Red in strains or species described as strict anaerobes, such as *Methanosarcinales*, some *Deltaproteobacteria*, some *Gammaproteobacteria*, *Bacteroidetes*, some *Chlorobi*, *Chrysiogenetes* and *Deferribacteres*. This point has already highlight for other O2Red, the Cytochrome bd that was find also in species described as strict anaerobes (some *Deltaproteobacteria* and two
*Methanosarcinales*). It is likely that they are not pseudogenes. This suggest that these enzymes could be functional and could be, maybe, involved in detoxification [37].

## Acknowledgements

We are grateful to Jose De La Torre (Department of Biology, San Francisco State University) for providing the sequences of HCO subunits of “Candidatus Nitrosocaldus yellowstonii”. We thank Najwa Taïb for discussions and helpful comments on this work. We thank also Pierre Garcia for having supplied us his tool of visualization of genomic contexts (GeneSpy).

This work has been performed using the computing facilities of the CC LBBE/PRABI.

## References

1. Pereira MM, Santana M, Teixeira M. A novel scenario for the evolution of haem-copper oxygen reductases. Biochim. Biophys. Acta [Internet]. 2001;1505:185-208. Available from: http://www.ncbi.nlm.nih.gov/pubmed/11334784

2. Brochier-Armanet C, Talla E, Gribaldo S. The multiple evolutionary histories of dioxygen reductases: Implications for the origin and evolution of aerobic respiration. Mol. Biol. Evol. 2009;26:285-97.

3. Ducluzeau A-L, Schoepp-Cothenet B, Lis R van, Baymann F, Russell MJ, Nitschke W. The evolution of respiratory O2/NO reductases: an out-of-the-phylogenetic-box perspective. J. R. Soc. Interface [Internet]. 2014 [cited 2014 Nov 7];11:20140196. Available from: http://rsif.royalsocietypublishing.org/content/11/98/20140196

4. Gribaldo S, Talla E, Brochier-Armanet C. Evolution of the haem copper oxidases superfamily: a rooting tale. Trends Biochem. Sci. [Internet]. 2009;34:375-81. Available from: http://www.ncbi.nlm.nih.gov/pubmed/19647436

5. Sousa FL, Alves RJ, Ribeiro MA, Pereira-Leal JB, Teixeira M, Pereira MM. The superfamily of heme-copper oxygen reductases: Types and evolutionary considerations. Biochim. Biophys. Acta - Bioenerg. 2012;1817:629-37.

6. Castresana J, Lübben M, Saraste M, Higgins DG. Evolution of cytochrome oxidase, an enzyme older than atmospheric oxygen. EMBO J. 1994;13:2516-25.

7. Ducluzeau A-L, van Lis R, Duval S, Schoepp-Cothenet B, Russell MJ, Nitschke W. Was nitric oxide the first deep electron sink? Trends Biochem. Sci. [Internet]. 2009;34:9-15. Available from: http://www.ncbi.nlm.nih.gov/pubmed/19008107

8. Rinke C, Schwientek P, Szczyrba A, Ivanova NN, Anderson IJ, Cheng J-F, et al. Insights into the phylogeny and coding potential of microbial dark matter. Nature [Internet]. Nature Research; 2013 [cited 2016 Dec 8];499:431-7. Available from: http://www.nature.com/doifinder/10.1038/nature12352

9. Brown CT, Hug LA, Thomas BC, Sharon I, Castelle CJ, Singh A, et al. Unusual biology across a group comprising more than 15% of domain Bacteria. Nature [Internet]. 2015;523:208-11. Available from: http://www.nature.com/doifinder/10.1038/nature14486

10. Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: A new generation of protein database search programs. Nucleic Acids Res. [Internet]. 1997 [cited 2014 Jul 9];25:3389-402. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=146917&tool=pmcentrez&rendertype=abstract

11. Eddy SR. A new generation of homology search tools based on probabilistic inference. Genome Inform. [Internet]. 2009;23:205-11. Available from: http://www.ncbi.nlm.nih.gov/pubmed/20180275

12. Sousa FL, Alves RJ, Pereira-Leal JB, Teixeira M, Pereira MM. A Bioinformatics Classifier and Database for Heme-Copper Oxygen Reductases. PLoS One [Internet]. 2011 [cited 2014 Nov 5];6:e19117. Available from: http://dx.doi.org/10.1371/journal.pone.0019117

13. Katoh K, Misawa K, Kuma K, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. [Internet]. 2002;30:3059-66. Available from: http://www.ncbi.nlm.nih.gov/pubmed/12136088
14. Gouy M, Guindon S, Gascuel O. SeaView version 4: A multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Mol. Biol. Evol. [Internet]. 2010;27:221–4. Available from: http://mbe.oxfordjournals.org/content/27/2/221.long

15. Thompson JD, Plewniak F, Ripp R, Thierry JC, Poch O. Towards a reliable objective function for multiple sequence alignments. J. Mol. Biol. [Internet]. 2001;314:937–51. Available from: http://www.sciencedirect.com/science/article/pii/S0022283601951873

16. Criscuolo A, Gribaldo S. BMGE (Block Mapping and Gathering with Entropy): a new software for selection of phylogenetic informative regions from multiple sequence alignments. BMC Evol. Biol. [Internet]. 2010 [cited 2015 Dec 7];10:1–21. Available from: http://link.springer.com/article/10.1186/1471-2148-10-210

17. Nguyen LT, Schmidt HA, Von Haeseler A, Minh BQ. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 2015;32:268–74.

18. Guindon S, Dufayard J-F, Lefort V, Anisimova M. New Alogrithms and Methods to Estimate Maximum-Likelihoods Phylogenies: Assessing the performance of PhyML 3.0. Syst. Biol. 2010;59:307–21.

19. Letunic I, Bork P. Interactive Tree Of Life (iTOL): an online tool for phylogenetic tree display and annotation. Bioinformatics [Internet]. 2007 [cited 2016 Jan 2];23:127–8. Available from: http://www.ncbi.nlm.nih.gov/pubmed/17050570

20. Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res. [Internet]. Oxford University Press; 2016 [cited 2016 Sep 19];44:W242-5. Available from: http://www.ncbi.nlm.nih.gov/pubmed/27095192

21. Jauffrit F, Penel S, Delmotte S, Rey C, de Vienne DM, Gouy M, et al. RiboDB Database: A Comprehensive Resource for Prokaryotic Systematics. Mol. Biol. Evol. [Internet]. Oxford University Press; 2016 [cited 2016 Aug 30];33:2170–2. Available from: http://mbe.oxfordjournals.org/lookup/doi/10.1093/molbev/msw088

22. Jelen BI, Giovannelli D, Falkowski PG. The Role of Microbial Electron Transfer in the Coevolution of the Biosphere and Geosphere. Annu. Rev. Microbiol. [Internet]. Annual Reviews 4139 El Camino Way, PO Box 10139, Palo Alto, California 94303-0139, USA ; 2016 [cited 2016 Aug 19];70:annurev-micro-102215-095521. Available from: http://www.annualreviews.org/doi/10.1146/annurev-micro-102215-095521

23. Farquhar J, Zerkle AL, Bekker A. Geological constraints on the origin of oxygenic photosynthesis. Photosynth. Res. [Internet]. 2011 [cited 2017 Mar 17];107:11–36. Available from: http://www.ncbi.nlm.nih.gov/pubmed/20882345

24. Shih PM. Cyanobacterial Evolution: Fresh Insight into Ancient Questions. Curr. Biol. [Internet]. 2015 [cited 2017 May 9];25:R192–3. Available from: http://www.sciencedirect.com.inee.bib.cnrs.fr/science/article/pii/S0960982214016492

25. Lee C-TA, Yeung LY, McKenzie NR, Yokoyama Y, Ozaki K, Lenardic A. Two-step rise of atmospheric oxygen linked to the growth of continents. Nat. Geosci. [Internet]. 2016 [cited 2016 Aug 19];9:417–24. Available from: http://www.nature.com/doifinder/10.1038/ngeo2707

26. Sukumaran P V. Evolution of the atmosphere and oceans: Evidence from geological records. Resonance [Internet]. Springer India; 2000 [cited 2017 Jan 19];5:6–14. Available from: http://link.springer.com/10.1007/BF02834667

27. Lyons TW, Reinhard CT 1, Planavsky NJ 1. The rise of oxygen in Earth's early ocean and atmosphere. Nat. Febr. 20, 2014 [Internet]. 2014 [cited 2014 Nov 10];506:307–15. Available from: http://gate1.inist.fr/login?url=http://ovid-sp.ovid.com/ovidweb.cgi?T=JS&CSC=Y&NEWS=N&PAGE=fulltext&D=ovfto&AN=00006056-201402200-00039

28. Luo G, Ono S, Beukes NJ, Wang DT, Xie S, Summons RE. Rapid oxygenation of Earths atmosphere 2.33 billion years ago. Sci. Adv. [Internet]. 2016 [cited 2017 Mar 13];2:e1600134–e1600134. Available from: http://advances.sciencemag.org/cgi/doi/10.1126/sciadv.1600134

29. Kendall B, Creaser RA, Reinhard CT, Lyons TW, Anbar AD. Transient episodes of mild environmental oxygenation and oxidative continental weathering during the late Archean. Sci. Adv. [Internet]. 2015;1:e1500777. Available from: http://advances.sciencemag.org/cgi/doi/10.1126/sciadv.1500777
30. Lu Z, Chang YC, Yin Q-Z, Ng CY, Jackson WM. Evidence for direct molecular oxygen production in CO2 photodissociation. Science (80-.). [Internet]. 2014 [cited 2014 Nov 5];346:61-4. Available from: http://www.ncbi.nlm.nih.gov/pubmed/25278605

31. Walker JCG. Oxygen and hydrogen in the primitive atmosphere. Pure Appl. Geophys. PAGEOPH [Internet]. Birkhäuser-Verlag; 1978 [cited 2017 Feb 28];116:222-31. Available from: http://link.springer.com/10.1007/BF01636879

32. Wang X-D, Gao X-F, Xuan C-J, Tian SX. Dissociative electron attachment to CO2 produces molecular oxygen. Nat. Chem. [Internet]. 2016 [cited 2017 Feb 28];8:258-63. Available from: http://www.ncbi.nlm.nih.gov/pubmed/26892558

33. Ducluzeau AL, Ouchane S, Nitschke W. The cbb3 oxidases are an ancient innovation of the domain bacteria. Mol. Biol. Evol. 2008;25:1158-66.

34. Ducluzeau A-L, Schoepp-Cothenet B, van Lis R, Baymann F, Russell MJ, Nitschke W. The evolution of respiratory O2/NO reductases: an out-of-the-phylogenetic-box perspective. J. R. Soc. Interface [Internet]. 2014;11:20140196. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=4233682&tool=pmcentrez&rendertype=abstract

35. López-García P, Zivanovic Y, Deschamps P, Moreira D. Bacterial gene import and mesophilic adaptation in archaea. Nat. Rev. Microbiol. [Internet]. 2015;13:447-56. Available from: http://www.ncbi.nlm.nih.gov/pubmed/26075362

36. Raymann K, Brochier-Armanet C, Gribaldo S. The two-domain tree of life is linked to a new root for the Archaea. Proc. Natl. Acad. Sci. U. S. A. [Internet]. National Academy of Sciences; 2015 [cited 2016 Aug 19];112:6670-5. Available from: http://www.ncbi.nlm.nih.gov/pubmed/25964353

37. Le Fourn C, Brasseur G, Brochier-Armanet C, Pieulle L, Brioukhanov A, Ollivier B, et al. An oxygen reduction chain in the hyperthermophilic anaerobe Thermotoga maritima highlights horizontal gene transfer between Thermococcales and Thermotogales. Environ. Microbiol. [Internet]. 2011 [cited 2017 Feb 27];13:2132-45. Available from: http://doi.wiley.com/10.1111/j.1462-2920.2011.02439.x
# Genomicus – New tools for comparative genomics and evolution in eukaryotes

Alexandra LOUIS¹, Nga Thi Thuy NGUYEN¹, Hugues ROEST CROLLIUS¹

¹ Ecole Normale Supérieure, Inserm, CNRS, IBENS, F-75005 Paris, France

Corresponding Author: alouis@biologie.ens.fr

Since 2010, the Genomicus [1,2] web server is available online at http://genomicus.biologie.ens.fr. This graphical browser provides access to comparative genomics analyses in four different phyla (Vertebrate, Plants, Fungi, and non vertebrate Metazoans). Users can manipulate evolutionary genomic information from extant species, as well as ancestral gene content and gene order for vertebrates and flowering plants.

Annotated genes are presented in a classical but highly customisable phylogenetic framework, while harnessing the added signal provided by local gene organisation. Gene presence/absence, Ka and Ks, % Id, genomic distance are available in addition to homology relationships. Genomicus is used by a broad community of users that perform >30,000 individual queries per months.

New analysis and visualisation tools have recently been implemented in Genomicus Vertebrate. Entire genomes (karyotype structures) can now be compared between multiple genomes, and synteny blocks can now be computed and visualised between any two genomes [3].

In this tutorial we will take users through the 6 different visualisation and analysis tools via a number of real test cases, aimed at answering questions on karyotype evolution, gene synteny conservation, gene duplication, gene gain, gene loss and chromosomal rearrangement.

## References

[1] Louis, A., Murat, F., Salse, J. and Roest Crollius, H. GenomicusPlants: A Web Resource to Study Genome Evolution in Flowering Plants. *Plant Cell Physiol.* 56: e4. 2015

[2] Louis, A., Nguyen, N.T., Muffato, M. and Roest Crollius, H. Genomicus update 2015: KaryoView and MatrixView provide a genome-wide perspective to multispecies comparative genomics. *Nucleic acids research* 43: D682-689. 2015

[3] Lucas, J. M., Muffato, M. & Roest Crollius, H. PhylDiag: identifying complex synteny blocks that include tandem duplications using phylogenetic gene trees. *BMC Bioinformatics* 15, 268 (2014).
Conférence invitée

Céline BROCHIER

Université Lyon 1, France

<figure><img src="image_28.png" /></figure>

The growing tree of Archaea: changing perspectives on the diversity and evolution of the third domain of life

Archaea occupy a key position in the Tree of Life, and represent a major fraction of the microbial diversity. Abundant in soils, ocean sediments and the water column, they are key players in processes mediating global carbon and nutrient fluxes, as well as important components of the animal microbiome and human body. The development of culture-independent sequencing techniques has revealed a myriad of so far inaccessible microbial lineages and filled up the archaeal tree with entirely new branches. The unprecedented access to genomic data from a large number of archaeal lineages provides the raw material for dissecting the origin of this domain, the evolutionary trajectories that have shaped its current diversity, and its relationships with Bacteria and Eukaryotes. This rainfall of data combined to cutting-edge methods allowing to disentangle the multiple signals contained in molecular sequences has shed new light on the evolutionary history of Archaea. Here I will review the major advances in the field as well as important open issues and future challenges.