# Comparison of statistical methods of inference of cooccurrence networks within micro-bial ecosystems from metagenomics data

Sophie Schbath, Mahendra Mariadassou, Julie Lao

► To cite this version:

Sophie Schbath, Mahendra Mariadassou, Julie Lao. Comparison of statistical methods of inference of cooccurrence networks within micro-bial ecosystems from metagenomics data. JOBIM 2017 - Journées Ouvertes Biologie Informatique Mathématiques, Jul 2017, Lille, France., pp.234, 2017, JOBIM 2017. hal-01563604

HAL Id: hal-01563604

https://hal.science/hal-01563604

Submitted on 2 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
Journées
Ouvertes en
Biologie,
Informatique et
Mathématiques

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

Lille, 3-6 juillet 2017

Éditeurs
Cédric Lhoussaine
Hélène Touzet
# Préface

Cette année, pour sa dix-huitième édition, JOBIM (Journées Ouvertes en Biologie, Informatique et Mathématiques) fait étape à Lille. Cet événement annuel, proposé conjointement par la Société Française de BioInformatique, le GDR BioInformatique Moléculaire et l'Institut Français de Bioinformatique, est devenu au fil du temps un rendez-vous incontournable de la communauté bioinformatique et biostatistique nationale. JOBIM offre une photographie unique des travaux menés dans l'Hexagone, et parfois même un peu plus loin. C'est un lieu d'échanges et de partages, toujours convivial, où se croisent chercheurs et ingénieurs, théoriciens et praticiens, jeunes et moins jeunes.

Comme à son habitude, le programme est très riche. Il reflète la diversité et la vitalité de la discipline : analyse de données omiques, phylogénie, biologie des systèmes, bioinformatique structurale, biologie intégrative,... Le comité de programme a sélectionné 37 présentations orales, dont 22 exposés classiques et 15 démonstrations logicielles, et 100 posters. Nous avons également le plaisir d'accueillir six orateurs invités : Céline Brochier-Armanet (LBBE, Lyon), Franca Fraternali (King's college, London), John Huelsenbeck (UC Berkeley), Tobias Marschall (MPI, Saarbrücke), Julio Saez-Rodrigues (EMBL-EBI, Aachen) et Patrick Wincker (Génoscope). Vous trouverez tous les résumés dans ce volume. Enfin, se tiendront en marge de JOBIM les journées de l'association JeBiF (Jeunes Bioinformaticiens de France).

Nous vous souhaitons une belle conférence!

Pour le comité de programme,
Cédric Lhoussaine et Hélène Touzet

Pour le comité d'organisation,
Guillemette Marot et Jean-Stéphane Varré
Comité de programme

Cédric LHOUSSAINE (CRISTAL, Université de Lille)

Hélène TOUZET (CNRS, CRISTAL, Université de Lille)

Eric BAPTESTE (CNRS, AIRE, Université Pierre et Marie Curie)

Séverine BÉRARD (ISEM, Université Montpellier)

Virginie BERNARD (Institut Curie, Paris)

Samuel BLANQUART (Inria Lille)

Jérémie BOURDON (LINA, Université de Nantes)

Bastien BOUSSAU (CNRS, LBBE, Université Lyon 1)

Christine BRUN (CNRS, TAGC, Aix-Marseille Université)

Ségolène CABOCHE (TAG/CIII, Université de Lille)

Frédéric CAZALS (Inria Sophia)

Annie CHATEAU (LIRMM, Université Montpellier)

Rayan CHIKHI (CNRS, CRISTAL, Université de Lille)

Vincent CHOURAKI (Inserm U1167)

Sarah COHEN-BOULAKIA (LRI, Université Paris-Sud)

Marie-Agnès DILLIES (Institut Pasteur, Paris)

Damien EVEILLARD (LINA, Université de Nantes)

Jérôme FERET (Inria / Ecole Normale Supérieure)

Martin FIGEAC (Université de Lille)

Christine GASPIN (Inra Toulouse)

Carito GUZIOLOWSKI (LS2N, Ecole Centrale de Nantes)

Laurent JACOB (CNRS, LBBE, Université Lyon 1)

Fabien JOURDAN (INRA Toulouse)

Michaël KOPP (I2M, Aix-Marseille Université)

Vincent LACROIX (LBBE, Université Lyon 1)

Valérie LECLÈRE (Institut Charles Viollette, Université de Lille)

Claire LEMAITRE (Inria Rennes)

Marc LENSINK (CNRS, UGSF, Université de Lille)

Anne LOPES (IGM, Université Paris Sud)

Macha NIKOSLKI (CNRS, LaBRI and CBIB, Université de Bordeaux)

Loïc PAULEVÉ (CNRS, LRI, Université Paris-Sud)

Eric PELLETIER (CEA / Genoscope)

Guy PERRIÈRE (CNRS, LBBE, Université Lyon 1)

Pierre PETERLONGO (Inria Rennes)

Pierre PEYRET (MEDIS, INRA-Université Clermont Auvergne)

Céline POUX (EEP, Université de Lille)

Hugues RICHARD (LCQB, Université Pierre et Marie Curie)

Delphine ROPERS (Inria Grenoble Rhône-Alpes)

Nicolas SERVANT (Institut Curie, Paris)

Anne SIEGEL (CNRS, IRISA, Université Rennes 1)

Hédi SOULA (Université Pierre et Marie Curie)

Pascal TOUZET (EEP, Université de Lille)

Jacques VAN HELDEN (TAGC, Aix-Marseille Université)

Vincent VANDEWALLE (EA 2694, Université de Lille)

Cristian VERSARI (CRISTAL, Université de Lille)
Comité d'organisation

Marie-Bénédicte DERNONCOURT (Inria Lille)

Guillemette MAROT (EA 2694, Université de Lille)

Jean-Stéphane VARRÉ (CRISTAL, Université de Lille)

Émilie ALLART (CRISTAL, Université de Lille)

Anaïs BARRAY (bilille, Université de Lille)

Aurélien BÉLIARD (CHRU Lille, CRISTAL, Université de Lille)

Samuel BLANCK (EA 2694, Université de Lille)

Ségolène CABOCHE (TAG/CII L, Université de Lille)

Hélène CHIAPELLO (Inra Toulouse, SFBI)

Christophe DEMAY (CHRU Lille)

Gaël EVEN (Genes Diffusion)

Sophie GALLINA (CNRS, EEP, Université de Lille)

Sébastien GREC (UGSF, Université de Lille)

Isabelle GUIGON (bilille, Université de Lille)

David HOT (Institut Pasteur de Lille)

Stéphane JANOT (CRISTAL, Université de Lille)

Laetitia JOURDAN (CRISTAL, Université de Lille)

Valérie LECLÈRE (Institut Charles Viollette, Université de Lille)

Pierre MARIJON (Inria Lille)

Laurent NOÉ (CRISTAL, Université de Lille)

Pierre PERICARD (CRISTAL, Université de Lille)

Céline POUX (EEP, Université de Lille)

Maude PUPIN (CRISTAL, Université de Lille)

Tatiana ROCHER (CRISTAL, Université de Lille)

Chadi SAAD (CRISTAL, JPARC, Université de Lille)

Mikaël SALSON (CRISTAL, Université de Lille)

Olivier SAND (UMR8199, Institut de Biologie de Lille)

Pascal TOUZET (EEP, Université de Lille)
# Table des matières

## Présentations orales 1

### Lundi 3 juillet

#### Conférence invitée

Holistic metagenomics in marine plankton communities
P. WINCKER .................................................................................................................... 2

#### Article

Simka : large scale de novo comparative metagenomics
G. BENOIT, P. PETERLONGO, M. MARIADASSOU, E. DREZEN, S. SCHBATH, D. LAVENIER et C. LEMAITRE .................................................................................................................... 3

#### Article

Food-Microbiomes Transfert, a shotgun metagenomic tool and a database to analyze cheese ecosystems
T. GUIRIMAND, A.-L. ABRAHAM, S. DEROZIER, C. PAUVERT, M. MARIADASSOU, V. LOUX et P. RENAULT .................................................................................................................... 6

#### Article

Reconstruction of full-length 16S rRNA sequences for taxonomic assignment in metagenomics
P. PERICARD, Y. DUFRESNE, S. BLANQUART et H. TOUZET .................................................... 10

#### Article

FEELnc : an alignment-free tool for long non-coding RNAs annotation
V. WUCHER, F. LEGEAI, B. HÉDAN, G. RIZK, L. LAGOUTTE, E. CADIEU, A. DAVID, N. BOTHEREL, C. LE BÉGUEC, C. ANDRÉ, C. HITTE et T. DERRIEN .................................................... 18

#### Article

Analyse intégrative des ARN longs non-codant (lncRNAs) du génome canin
C. LE BÉGUEC, V. WUCHER, L. LAGOUTTE, E. CADIEU, B. HÉDAN, C. ANDRÉ, C. HITTE et T. DERRIEN .................................................................................................................... 20

#### Article

Analysis, Integration and Modeling of Cell Clustering Results in High-Dimensional Cytometry Data
G. GAUTREAU, D. PEJOSKI, R. LE GRAND, A. COSMA, A.-S. BEIGNON et N. TCHITCHEK .................................... 28
# Mardi 4 juillet

## Conférence invitée

A Guided Tour to Computational Haplotyping

<table><tr><td>T. MARSCHALL</td><td>30</td></tr></table>

## Article

GenomeOnRails : depicting microbial species diversity via a pangenome graph

<table><tr><td>G. GAUTREAU, R. PLANEL, A. PERRIN, M. TOUCHON, E. ROCHA, C. AMBROISE, C. MATIAS, S. CRUVEILLER, C. MÉDIGUE et D. VALLENET</td><td>31</td></tr></table>

## Article

Assembly of heterozygous genomes with high-order de Bruijn graphs

<table><tr><td>A. LIMASSET, C. MARCHET, P. PETERLONGO et J.-F. FLOT</td><td>32</td></tr></table>

## Article

Rapid spectra comparison with data-mining algorithms : accessing post translational modification profiles on a sample scale

<table><tr><td>M. DAVID, G. FERTIN, H. ROGNIAUX et D. TESSIER</td><td>33</td></tr></table>

## Conférence invitée

Network Models to Understand and Combat Cancer : from Clinical Genomics to Biochemical Modelling

<table><tr><td>J. SAEZ-RODRIGUEZ</td><td>35</td></tr></table>

## Article

Rare and common epilepsies converge on a shared gene regulatory network providing opportunities for novel antiepileptic drug discovery

<table><tr><td>A. DELAHAYE-DURIEZ, P. SRIVASTAVA, K. SHKURA, S. LANGLEY, L. LAANISTE, A. MORENO-MORAL, B. DANIS, M. MAZZUFERI, P. FOERCH, E. GAZINA, K. RICHARDS, S. PETROU, R. KAMINSKI, E. PETRETTO et M. JOHNSON</td><td>36</td></tr></table>

## Article

Prediction of Disease-associated Genes by advanced Random Walk with Restart on Multiplex and Heterogeneous Biological Networks

<table><tr><td>A. VALDEOLIVAS, E. REMY, L. TICHIT, G. ODELIN, C. NAVARRO, S. PERRIN, P. CAU, N. LEVY et A. BAUDOT</td><td>38</td></tr></table>

## Démonstration

CRCmapper : models of core transcriptional regulatory circuitries

<table><tr><td>V. SAINT-ANDRÉ, A. FEDERATION, C. LIN, B. ABRAHAM, J. REDDY, T. IHN LEE, J. BRADNER et R. YOUNG</td><td>44</td></tr></table>

## Article

Network Inference of Dynamic Models by the Combination of Spanning Arborescences

<table><tr><td>A. COUTANT et C. ROUVEIROL</td><td>45</td></tr></table>
## Démonstration

### Long-term Tracking of Budding Yeast Cells with CellStar

C. VERSARI, S. STOMA, K. BATMANOV, A. LLAMOSI, F. MROZ, A. KACZMAREK, M. DEYELL, C. LHOUSSAINE, P. HERSEN et G. BATT

53

## Démonstration

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

C. BECAVIN, M. KOUTERO, N. TCHITCHEK, F. CERUTTI, P. LECHAT, N. MAILLET, C. HOEDE, H. CHIAPELLO, C. GASPIN et P. COSSART

54

## Démonstration

### AskOmics, a web tool to integrate and query biological data using semantic web technologies

X. GARNIER, A. BRETAUDEAU, O. FILANGI, F. LEGEAI, A. SIEGEL et O. DAMERON

56

## Démonstration

### Regulatory and signaling network assembly through Linked Open Data

M. LEFEBVRE, J. BOURDON, C. GUZIOLOWSKI et A. GAIGNARD

57

## Article

### Scientific Workflows for Computational Reproducibility in the Life Sciences : Status, Challenges and Opportunities

S. COHEN-BOULAKIA, K. BELHAJJAME, O. COLLIN, J. CHOPARD, C. FROIDEVAUX, A. GAIGNARD, K. HINSEN, P. LARMANDE, Y. LE BRAS, F. LEMOINE, F. MAREUIL, H. MÉNAGER, C. PRADAL et C. BLANCHET

58

## Démonstration

### Sequanix : A Dynamic Graphical Interface for Snakemake Workflows

D. DESVILLECHABROL, R. LEGENDRE, C. BOUCHIER, S. KENNEDY et T. COKELAER

60

## Démonstration

### Biosphere : un portail web haut niveau pour une utilisation bioinformatique des clouds

B. BRANCOTTE, M. BEDRI, J. LORENZO, S. PERRIN, F. SÉNÉ, A. SEPOU NGAÏLO, C. BLANCHET et J.-F. GIBRAT

61

## Démonstration

### Biodjango, an open framework for bioinformatics publishing

E. GHEYOUCHE et S. TÉLETCHÉA

62

## Démonstration

### New Generation Phylogeny.fr : Refactoring Phylogeny.fr for Innovative Phylogenetic Services

D. CORREIA, V. LEFORT, O. DOPPELT-AZEROUAL, F. MAREUIL, S. COHEN-BOULAKIA et O. GASCUEL

63
## Mercredi 5 juillet

### Conférence invitée

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

F. FRATERNALI 64

### Article

Use of cross-docking simulations for identification of protein-protein interactions sites : the case of proteins with multiple binding sites.

N. LAGARDE, L. VAMPARYS, B. LAURENT, A. CARBONE et S. SACQUIN-MORA 65

### Article

In silico developments for the study of glycosylation applied to extracellular matrix proteins

C. BESANÇON, A. GUILLOT, S. BLAISE, M. DAUCHEZ, J. JONQUET, N. BELLOY et S. BAUD 66

### Article

HG-CoLoR : Hybrid Graph for the error Correction of Long Reads

P. MORISSE, T. LECROQ et A. LEFEBVRE 67

### Démonstration

PhylOligo : a package to identify contaminant or untargeted organism sequences in genome assemblies

L. MALLET, T. BITARD-FEILDEL, F. CERUTTI et H. CHIAPELLO 75

### Démonstration

Dynamix : Dynamic visualization by automatic selection of informative tracks from hundreds of genomic data sets

M. MONFORT, E. FURLONG et C. GIRARDOT 77

### Démonstration

GeneSpy, a simple tool to explore genomic context

P. SIMON GARCIA, F. JAUFFRIT, C. GRANGEASSE et C. BROCHIER-ARMANET 78

### Démonstration

MCXpress : An R Package for functional interpretation of single cell RNA-Seq data using multivariate analysis

A. CORTAL et A. RAUSELL 79

### Démonstration

Deciphering The Functional Effects of Genetic Variation With UniProt Annotations

B. BELY, A. NIGHTINGALE et M. MARTIN 80
Conférence invitée

<table><tr><td>Bayesian inference in phylogeny for genome-scale data<br>J. HUELSENBECK</td><td>81</td></tr></table>

Article

<table><tr><td>Extreme halophilic archaea derive from two distinct methanogen Class II ancestors<br>M. AOUAD, N. TAIB, A. OUDART, M. LECOCQ, M. GOUY et C. BROCHIER-ARMANET</td><td>82</td></tr></table>

Article

<table><tr><td>Origin and evolution of multiple haem copper oxidases in Archaea<br>A. OUDART, S. GRIBALDO et C. BROCHIER-ARMANET</td><td>90</td></tr></table>

Démonstration

<table><tr><td>Genomicus New tools for comparative genomics and evolution in eukaryotes<br>A. LOUIS, N. NGUYEN et H. ROEST CROLLIUS</td><td>98</td></tr></table>

Jeudi 6 juillet

Conférence invitée

<table><tr><td>The growing tree of Archaea : changing perspectives on the diversity and evolution of the third domain of life<br>C. BROCHIER</td><td>99</td></tr></table>

Article

<table><tr><td>Probing factor-dependent long-range contacts using regression with higher-order interaction terms<br>R. MOURAD, L. LI et O. CUVIER</td><td>100</td></tr></table>

Article

<table><tr><td>An integrative approach for predicting the RNA secondary structure for the HIV-1 Gag UTR using probing data<br>A. SAAIDI, Y. PONTY et B. SARGUEIL</td><td>102</td></tr></table>

Article

<table><tr><td>NCboost : a meta-classifier of pathogenic non-coding variants integrating multiple sequence conservation and unsupervised functional scores<br>B. CARON, L. SLIM et A. RAUSELL</td><td>103</td></tr></table>

Article

<table><tr><td>Context-specific prioritization of non-coding variants implicated in human diseases<br>L. MOYON, Y. CLÉMENT, C. BERTHELOT et H. ROEST CROLLIUS</td><td>104</td></tr></table>
## Posters

105

### Poster A1

Analysis workflow for smallRNA-seq data

S. NIN, S. RIALLE, E. DUBOIS et L. JOURNOT

106

### Poster A2

Analyse du réseau moléculaire impliqué dans le remodelage ventriculaire gauche postinfarctus du myocarde

M. CUVELLIEZ, C. BAUTERS, T. KELDER, M. RADONJIC, P. AMOUYEL et F. PINET

107

### Poster A3

Mitochondrial genome variability of 205 Arabian endurance horses

A. HEURTEAU, C. HOEDE, A. RICARD, D. ESQUERRÉ, M. CAROLINE, N. MACH, C. ROBERT et E. BARREY

108

### Poster A4

Epigenetic heterogeneity in multiple myeloma

J. RONDINEAU, V. GABORIT, C. GUERIN-CHARBONNEL, P. MOREAU, S. MINVIELLE et F. MAGRANGEAS

109

### Poster A5

NF-kappaB Landscape in Multiple Myeloma by high-throughput sequencing analysis

V. GABORIT, J. RONDINEAU, W. GOURAUD, J. BOURDON, S. MINVIELLE et F. MAGRANGEAS

110

### Poster A6

Efficient MinION data management

L. JOURDREN, A. BIRER, L. FERRATO-BERBERIAN, S. LEMOINE et S. LE CROM

111

### Poster A7

Mise en place d'un pipeline de contrôle de la qualité de runs MinIon

L. FERRATO-BERBERIAN, A. BIRER, A. MOHAMMAD, C. BLUGEON, F. COULPIER, S. LE CROM, L. JOURDREN et S. LEMOINE

112

### Poster A8

Impact et évolution de la correction d'erreur sur des lectures longues issues de séquençage MinION Oxford Nanopore dans un contexte transcriptomique

L. FERRATO-BERBERIAN, A. BIRER, S. LE CROM, L. JOURDREN et S. LEMOINE

113

### Poster A9

Toullig : New pipeline for nanopore data analysis

A. BIRER, L. FERRATO-BERBERIAN, S. LE CROM, S. LEMOINE et L. JOURDREN

114

### Poster A10

Extreme phenotypes define epigenetic and metabolic myeloid signatures in cardiovascular diseases

D. SEYRES, J. LAMBOURNE, P. KIRK, A. CABASSI, F. BURDEN, R. KREUZHUBER, S. FARROW, C. KEMPSTER, H. MCKINNEY, A. PARK, D. SAVAGE, J. GRIFFIN, O. STEGLE, S. RICHARDSON, K. DOWNES, W. OUWEHAND et M. FRONTINI

115
<table>
  <tbody>
    <tr>
      <td>Poster A11<br>Analyse et co-développement d'outils de bioinformatique destinés au traitement de données NGS issues de métagénomique virale<br>E. DELPUECH, G. CROVILLE, J.-L. GUERIN, C. KLOPP et S. MAMAN</td>
      <td>116</td>
    </tr>
    <tr>
      <td>Poster A12<br>3-SMART : Bioinformatic analysis of intronic polyadenylation regulation<br>M. CADIX, I. TANAKA, P. GESTRAUD, M. SÉJOURNÉ, S. VAGNER, N. SERVANT et M. DUTERTRE</td>
      <td>117</td>
    </tr>
    <tr>
      <td>Poster A13<br>Developing and sharing reproducible bioinformatics pipelines : best practices<br>Y. LELIÈVRE, A. BIHOUÉE, E. CHARPENTIER, A. GAIGNARD, S. SOUCHET et D. VINTACHE</td>
      <td>118</td>
    </tr>
    <tr>
      <td>Poster A15<br>MICROSCOPE : an integrated platform for the Exploration and Curation of Microbial Genomes<br>D. VALLENET, A. CALTEAU, S. CRUVEILLER, M. GACHET, G. GAUTREAU, A. JOSSO, A. LAJUS, J. LANGLOIS, J. MERCIER, H. PEREIRA, R. PLANEL, J. ROLLIN ROLLIN, D. ROCHE, Z. ROUY et C. MÉDIGUE</td>
      <td>119</td>
    </tr>
    <tr>
      <td>Poster A16<br>NETSYN : NETwork SYNteny, a new tool to help functional annotation<br>B. VIART, K. BASTARD, G. REBOUL, H. PEREIRA, R. PLANEL, M. STAM, C. MÉDIGUE et D. VALLENET</td>
      <td>120</td>
    </tr>
    <tr>
      <td>Poster A17<br>Stratégie d'assemblage de génomes en cellule unique de protistes marins dans le cadre du projet Tara Oceans<br>L. D'AGATA, Y. SEELEUTHNER, J. POULAIN, P. WINCKER et J.-M. AURY</td>
      <td>121</td>
    </tr>
    <tr>
      <td>Poster A18<br>Assessing the functional impact of genomic alterations using proteogenomics<br>G. BEDRAN, Y. VANDENBROUCK, E. BONNET, J.-F. DELEUZE, D. PFLIEGER et C. BATTAIL</td>
      <td>122</td>
    </tr>
    <tr>
      <td>Poster A19<br>Deciphering early cell fate decision by Single Cell RNA-Seq and DGE-Seq<br>D. MEISTERMANN, Y. LELIÈVRE, E. CHARPENTIER, S. KILENS, T. FRÉOUR, J. BOURDON et L. DAVID</td>
      <td>123</td>
    </tr>
    <tr>
      <td>Poster A20<br>In search of W sex chromosome-specific sequences in the genome of the isopod crustacean Armadillidium vulgare<br>M. CHEBBI, T. BECKING, I. GIRAUD, B. MOUMEN, C. GILBERT, J. PECCOUD et R. CORDAUX</td>
      <td>124</td>
    </tr>
    <tr>
      <td>Poster A21<br>Exome sequencing to identify the molecular mechanism underlying chordomas pathogenesis<br>Z. TARIQ, K. DRIOUCH, V. BERNARD, I. BIECHE, V. RAYNAL, S. BAULANDE, H. MAMMAR et J. MASLIAH PLANCHON</td>
      <td>125</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>Poster A22<br>Revisiting cell lineage specification during male sex determination with single-cell RNA sequencing<br>I. STÉVANT, Y. NEIRJINCK, C. BOREL, J. ESCOFFIER, L. SMITH, S. ANTONARAKIS, E. DERMITZAKIS et S. NEF</td>
      <td>126</td>
    </tr>
    <tr>
      <td>Poster A23<br>How to analyse human or mouse Genome-scale CRISPR Knock-Out (GeCKO) datasets?<br>M. DELOGER, P. GESTRAUD, R. MARGUERON et N. SERVANT</td>
      <td>127</td>
    </tr>
    <tr>
      <td>Poster A24<br>Functional characterization of human epidermis co-expression modules<br>G. LEMOINE, M. SCOTT-BOYER, M. LECLERQ et A. DROIT</td>
      <td>128</td>
    </tr>
    <tr>
      <td>Poster A25<br>Integrated and flexible analysis of scRNA-seq data with Eoulsan scRNA-seq<br>G. BRELURUT, N. LEHMANN, C. HERNANDEZ, M. THOMAS-CHOLLIER, D. THIEFFRY, S. LE CROM et L. JOURDREN</td>
      <td>129</td>
    </tr>
    <tr>
      <td>Poster A26<br>Efficient data structure for indexing and similarity computation of nucleic sequences<br>C. MARCHET, A. MENG, L. LECOMPTE, A. LIMASSET, L. BITTNER et P. PETERLONGO</td>
      <td>130</td>
    </tr>
    <tr>
      <td>Poster A27<br>Detection of poly-adenylation sites from RNA-Seq data<br>C. FOURNIER et A. NECSULEA</td>
      <td>131</td>
    </tr>
    <tr>
      <td>Poster A28<br>Développement d'une structure pour l'indexation et la compression de multi-génomes<br>C. AGRET, M. RUIZ et A. MANCHERON</td>
      <td>132</td>
    </tr>
    <tr>
      <td>Poster A29<br>Epigenetic marks and the human transcriptome diversity<br>G. DEVAILLY, A. MANTSOKI et A. JOSHI</td>
      <td>133</td>
    </tr>
    <tr>
      <td>Poster A30<br>Microbial diversity and plant cell wall-degrading enzyme dynamics during dew-retting of flax - one of the oldest applications of biotechnology to textile transformation<br>C. DJEMIEL, S. GREC et S. HAWKINS</td>
      <td>134</td>
    </tr>
    <tr>
      <td>Poster A31<br>Splicing Lore: Speeding up the identifications of splicing factors regulating alternative exons across physiological and pathological conditions<br>H. POLVÈCHE, N. BOUCHOUICHA et D. AUBOEUF</td>
      <td>135</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>Poster A32<br>A method for DNA virus detection and quantification during pregnancy based on noninvasive prenatal testing whole genome sequencing results<br>V. CHESNAIS, A. OTT, E. CHAPLAIS, S. GABILLARD, C. VAULOUP-FELLOUS, A. BENACHI, J.-M. COSTA et E. GINOUX</td>
      <td>136</td>
    </tr>
    <tr>
      <td>Poster A33<br>RINSpector : a Cytoscape app that combines centrality analyses with DynaMine flexibility prediction<br>G. BRYSBAERT, K. LORGOUILLOUX, W. VRANKEN et M. LENSINK</td>
      <td>137</td>
    </tr>
    <tr>
      <td>Poster A34<br>Automated generation and analysis of parametric kinetic models obtained from biochemical interaction maps<br>M. BUFFARD, O. ORTEGA, C. LOPEZ et O. RADULESCU</td>
      <td>138</td>
    </tr>
    <tr>
      <td>Poster A35<br>Étude de processus de coalescence dans les paysages contemporains : bibliothèques template C++ pour le Calcul Bayésien Approché<br>A. BECHELER, C. CORON et S. DUPAS</td>
      <td>139</td>
    </tr>
    <tr>
      <td>Poster A36<br>PROqPCR : a Shiny web application for PROcessing of qRT-PCR data<br>M. SAUTREUIL et C. BÉRARD</td>
      <td>140</td>
    </tr>
    <tr>
      <td>Poster A37<br>Bivariate Negative Binomial Mixture Model for the analysis of RNA-seq data<br>M. SAUTREUIL, N. VERGNE, A. CHANNAROND, A. ROCHE, G. CHAGNY et C. BÉRARD</td>
      <td>141</td>
    </tr>
    <tr>
      <td>Poster A38<br>Quality evaluation of the mapping of small RNA-footprinting reads<br>P. FOURGOUX, E. DELANNOY, C. LURIN, G. RIGAILL et V. BRUNAUD</td>
      <td>142</td>
    </tr>
    <tr>
      <td>Poster A39<br>Medical diagnosis pipelines on the new AP-HP bioinformatics platform<br>J. BRAYET, C. BARETTE, M. BARTHELEMY, V. DESHAIES et A. LERMINE</td>
      <td>143</td>
    </tr>
    <tr>
      <td>Poster A40<br>Mathematical modeling of a genetic network controlling the regulation of Fe-S biogenesis<br>F. HAMMAMI, F. BARRAS, P. MANDIN et E. REMY</td>
      <td>144</td>
    </tr>
    <tr>
      <td>Poster A41<br>First gene-annotation enrichment analysis based on bacterial core-genome variants : Insights into mammalian and avian host adaptation of Salmonella serovars<br>M. VILA NOVA, K. DURIMEL, A. FELTEN, L. GUILLIER, M.-Y. MISTOU et N. RADOMSKI</td>
      <td>145</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>Poster A42<br>A phylogenomic network approach to decipher bacterial adaptation through horizontal gene transfer<br>D. RICHARD, V. RAVIGNÉ, A. CHABIRAND, O. PRUVOST et P. LEFEUVRE</td>
      <td>146</td>
    </tr>
    <tr>
      <td>Poster A43<br>Méthodes et outils de construction de super-arbres en phylogénie<br>M. SOULIÉ, V. LEFORT et A.-M. ARIGON CHIFOLLEAU</td>
      <td>148</td>
    </tr>
    <tr>
      <td>Poster A44<br>Implémentation d'une interface pratique pour l'évaluation de stratégies d'exploration de la diversité moléculaire par protéogénomique<br>Y. COGNE, C. ALMUNIA, O. PIBLE, D. GOUVEIA, A. CHAUMOT, O. GEFFARD et J. ARMENGAUD</td>
      <td>149</td>
    </tr>
    <tr>
      <td>Poster A45<br>Pea genetic map enrichment with Genotyping By Sequencing markers<br>A. KOUGBEADJO, G. AUBERT, M. FALQUE, D. EDWARDS, P. BAYER, J. BATLEY, M. ZANDER, S. HAYASHI, J. KREPLAK, J. BURSTIN et J. KREPLAK</td>
      <td>150</td>
    </tr>
    <tr>
      <td>Poster A46<br>Développement d'une mesure du biais d'usage des codons : application aux virus humains<br>J. BOURRET, S. ALIZON et I. BRAVO</td>
      <td>151</td>
    </tr>
    <tr>
      <td>Poster A47<br>Was the Chlamydial Adaptative Strategy to Tryptophan Starvation an Early Determinant of Plastid Endosymbiosis?<br>U. CENCI, M. DUCATEZ, D. KADOCHE, C. COLLEONI et S. BALL</td>
      <td>152</td>
    </tr>
    <tr>
      <td>Poster A48<br>OrthoInspector 3.0 : orthology en route to big data<br>Y. NEVERS, A. KRESS, R. RIPP, O. POCH et O. LECOMPTE</td>
      <td>153</td>
    </tr>
    <tr>
      <td>Poster A49<br>Déploiement automatique d'une infrastructure complexe pour le mapping des données de séquençage<br>S. PERRIN, B. BRANCOTTE, J. LORENZO, C. BLANCHET et J.-F. GIBRAT</td>
      <td>154</td>
    </tr>
    <tr>
      <td>Poster A50<br>State of the art and comparison of long reads technologies<br>A. POIRAUDEAU, M. MANNO, C. VANDECASTEELE, A. ROULET, C. ROQUES, M. VIDAL, C. ZANCHETTA, P. HEUILLARD, F. ROUX, B. MAYJONADE, J. GOUZY, Y. GUIGUEN, C. KLOPP, P. FRASSE, M. ZOUINE, C. DONNADIEU, O. BOUCHEZ, G. SALIN et C. KUCHLY</td>
      <td>155</td>
    </tr>
    <tr>
      <td>Poster A51<br>Genomic markers of species diversification in vertebrates<br>G. LOUVEL, E. LEWITUS, H. MORLON et H. ROEST CROLLIUS</td>
      <td>156</td>
    </tr>
  </tbody>
</table>
Poster A52
PREMS / ELVIS : A local plant biological resource management system
F. DUPUIS, A. LELIÈVRE, S. PELLETIER, T. THOUROUDE, J. BOURBEILLON et S. GAILLARD ... 157

Poster A53
Découverte et analyse de polymorphismes SNPs issus de RNA-seq chez le peuplier noir
O. ROGIER, S. AMANZOUGARENE, M.-C. LESAGE-DESCAUSES, S. BALZERGUE, V. BRUNAUD, J. CAIUS, A. CHATEIGNER, L. SOUBIGOU-TACONNAT, V. JORGE et V. SEGURA ... 158

Poster A54
DiNAMO : Exact method for degenerated IUPAC motifs discovery : characterization of sequence-specific errors
C. SAAD, L. NOÉ, H. RICHARD, J. LECLERC, M.-P. BUISINE, H. TOUZET et M. FIGEAC ... 159

Poster A55
Evolutionary conservation of unusual N-glycosylation sites in the human glycosyltransferase B4GALNT2
V. COGEZ, A. BARRAY, J. DE RUYCK, S. GROUX-DEGROOTE et A. HARDUIN-LEPERS ... 160

Poster A56
Comparison of statistical methods of inference of cooccurrence networks within microbial ecosystems from metagenomics data
J. LAO, M. MARIADASSOU et S. SCHBATH ... 161

Poster A57
ShRCAn : a user-friendly Shiny application for quantitative metagenomics analysis
F. THIRION, E. LE CHATELIER, N. PONS, A.-S. ALVAREZ, P. LÉONARD et D. EHRLICH ... 162

Poster A58
Genetic diversity of Anaplasma phagocytophilum among ticks and roe deer in a fragmented agricultural landscape
A. CHASTAGNER, A. PION, H. VERHEYDEN, B. LOUTRET, B. CARGNELUTTI, D. PICOT, V. POUX, J. BARD, O. PLANTARD, K. MCCOY, A. LEBLOND, G. VOURC'H et X. BAILLY ... 163

Poster A59
A transcriptional study of five fungal Mucor strains
A. LEBRETON, L. MESLET-CLADIÈRE, J.-L. JANY, G. BARBIER et E. CORRE ... 164

Poster A60
Multi-Cloud deployment for microbial genomes analysis
J. LORENZO, B. BRANCOTTE, T. LACROIX, M. BEDRI, J.-F. GIBRAT et C. BLANCHET ... 165

Poster A61
Evaluation of 15 in silico prediction tools for the classification of MED13L missense variations
T. SMOL, C. THUILLIER, S. MANOUVRIER-HANU et J. GHOUMID ... 166

Poster A62
HiFit : robust data analysis method for High-throughput qPCR
M. BAHIN, Q. VIAUTOUR, E. DIAZ, B. DUCOS et A. GENOVESIO ... 167
<table><tr><td>Poster A63<br>Is UniProtKB Missing Knowledgeable Proteins?<br>B. BELY, S. BENMOHAMMED, G. QI, N. TYAGI et M. MARTIN</td><td>168</td></tr></table>

<table><tr><td>Poster A64<br>Co-option of complex molecular system in bacterial membranes<br>R. DENISE, S. ABBY et E. ROCHA</td><td>169</td></tr></table>

<table><tr><td>Poster A65<br>High-quality, fast, and memory-efficient assembly of metagenomes and large genomes using Minia-pipeline<br>R. CHIKHI, C. DELTEL, G. RIZK, C. LEMAITRE, P. PETERLONGO, K. SAHLIN, L. ARVESTAD, P. MEDVEDEV et D. LAVENIER</td><td>170</td></tr></table>

<table><tr><td>Poster A66<br>Debugging long-read genome and metagenome assemblies using string graph analysis<br>P. MARIJON, J.-S. VARRÉ et R. CHIKHI</td><td>171</td></tr></table>

<table><tr><td>Poster A67<br>Labsquare une communauté de développeurs libres<br>S. SCHUTZ, P. MARIJON, J. ROQUET, F. PINA-MARTINS, J. SALLOU, E. TROUNE, A.-S. DENOMME et O. GUEUDELOT</td><td>172</td></tr></table>

<table><tr><td>Poster A68<br>Influence of SNP coding on the analysis of disease risk<br>H. SARTER, C. GOWER-ROUSSEAU et G. MAROT</td><td>173</td></tr></table>

<table><tr><td>Poster A69<br>Ultra High throughput, single molecule mapping of replicating DNA<br>N. MENEZES BRAGANCA, F. DE CARLI, W. BERRABAH, V. BARBE, A. GENOVESIO et O. HYRIEN</td><td>174</td></tr></table>

<table><tr><td>Poster A70<br>Data updates on Norine, the reference Non-Ribosomal Peptide knowledge base<br>Y. DUFRESNE, J. MICHALIK, A. FLISSI, V. LECLÈRE et M. PUPIN</td><td>175</td></tr></table>

<table><tr><td>Poster A71<br>NRPro : a Bioinformatics Tool for Nonribosomal Peptides Identification by Tandem Mass Spectrometry<br>E. RICART, M. CHEVALIER, M. PUPIN, V. LECLÈRE, C. FLAHAUT et F. LISACEK</td><td>176</td></tr></table>

<table><tr><td>Poster A72<br>Multi-scale characterization of symbiont diversity in the pea aphid complex through metagenomic approaches<br>C. GUYOMAR, F. LEGEAI, C. MOUGEL, C. LEMAITRE et J.-C. SIMON</td><td>177</td></tr></table>

<table><tr><td>Poster A73<br>iMOMi : a database dedicated to integration and exploration of multi-omic data<br>N. PONS, J.-M. BATTO, A. GHOZLANE, K. WEISZER, P. LÉONARD, E. LE CHATELIER, P. RENAULT et S. DUSKO EHRLICH</td><td>178</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>Poster A74<br>Sequence composition-based read binning and taxonomic profiling in infectious metagenomics diversity analyses<br>M. VANDENBOGAERT, C. BALIÈRE et V. CARO</td>
      <td>179</td>
    </tr>
    <tr>
      <td>Poster A75<br>REGOVAR, logiciel libre pour l'analyse de données de séquençage haut débit pour les maladies génétiques rares<br>A.-S. DENOMMÉ-PICHON, O. GUEUDELOT, J. ROQUET, J. SALLOU et D. GOUDENÈGE</td>
      <td>180</td>
    </tr>
    <tr>
      <td>Poster A76<br>Alignement à grande échelle pour une approche métagénomique dans le cadre du projet Tara Oceans<br>A. KOURLAIEV, C. DA SILVA, S. ENGELEN, A. BERTRAND, A. BRUNO, E. PELLETIER, P. WINCKER et J.-M. AURY</td>
      <td>181</td>
    </tr>
    <tr>
      <td>Poster A77<br>Seqenv : linking sequences to environments through text mining<br>L. SINCLAIR, U. IJAZ, L. JENSEN, M. COOLEN, C. GUBRY-RANGIN, A. CHROKOV, A. OULAS, C. PAVLOUDI, J. SCHNETZER, A. WEIMANN, A. IJAZ, A. EILER, C. QUINCE et E. PAFILIS</td>
      <td>182</td>
    </tr>
    <tr>
      <td>Poster A78<br>Predicting the ecological quality status of marine environments from eDNA metabarcoding data using supervised machine learning<br>T. CORDIER, P. ESLING, F. LEJZEROWICZ, J. VISCO, A. OUADAHI, C. MARTINS, T. CEDHAGEN et J. PAWLOWSKI</td>
      <td>183</td>
    </tr>
    <tr>
      <td>Poster A79<br>Extended HLA haplotypes in Membranous Nephropathy<br>L. CHAUVIÈRE, P. RONCO et H. DEBIEC</td>
      <td>184</td>
    </tr>
    <tr>
      <td>Poster A80<br>RSAT matrix-clustering : dynamic exploration and redundancy reduction of transcription factor binding motif collections<br>J. CASTRO-MONDRAGON, S. JAEGER, D. THIEFFRY, M. THOMAS-CHOLLIER et J. VAN HELDEN</td>
      <td>185</td>
    </tr>
    <tr>
      <td>Poster A81<br>Assemblage de novo de quasi-espèces virales basé sur un graphe de chevauchements<br>J. BAAIJENS, A. EL AABIDINE, E. RIVALS et A. SCHOENHUTH</td>
      <td>186</td>
    </tr>
    <tr>
      <td>Poster A82<br>Characterization of biological data<br>A. CHAMBON, T. BOUREAU, F. LARDEUX et F. SAUBION</td>
      <td>187</td>
    </tr>
    <tr>
      <td>Poster A83<br>Conservation of interaction energy landscapes across structural homologs through cross-docking calculations<br>H. SCHWEKE, S. SACQUIN-MORA, M.-H. MUCCHIELLI et A. LOPES</td>
      <td>188</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>Poster A84</td>
      <td>Dynamic model of Central Carbon Metabolism and electrochemical reactions of the cell</td>
      <td></td>
    </tr>
    <tr>
      <td></td>
      <td>C. MOULIN, J. DA VEIGA MOREIRA, E. BIGAN, L. SCHWARTZ, M. JOLICOEUR, L. TOURNIER et S. PERES</td>
      <td>189</td>
    </tr>
    <tr>
      <td>Poster A85</td>
      <td>Proposition d'un workflow d'analyse QIIME dans Galaxy et évaluation de trois techniques d'extraction d'ADN pour l'analyse du microbiote intestinal 16S</td>
      <td></td>
    </tr>
    <tr>
      <td></td>
      <td>S. BUFFET-BATAILLON, M. BEN ABDALLAH, P. BORDRON, E. CORRE et S. KAYAL</td>
      <td>190</td>
    </tr>
    <tr>
      <td>Poster A86</td>
      <td>BIOSPECIMENS 1.5 : a web platform to facilitate collaborative research on infectious diseases</td>
      <td></td>
    </tr>
    <tr>
      <td></td>
      <td>K. LOUIS, B. RIMBAULT, C. DELESTRE, Y. MOUSCAZ, M. ALBRIEUX, C. BOISSE, R. VILLET et G. BOISSY</td>
      <td>191</td>
    </tr>
    <tr>
      <td>Poster A87</td>
      <td>Host tropism and host-pathogen interplay of typhoidal Salmonella enterica</td>
      <td></td>
    </tr>
    <tr>
      <td></td>
      <td>L. MALLET, C. HOEDE, F. CERUTTI, A. MOISAN, C. GASPIN, I. VIRLOGEUX-PAYANT, I. SHOMER, O. GAL-MOR, T. SCHIEX et H. CHIAPELLO</td>
      <td>192</td>
    </tr>
    <tr>
      <td>Poster A88</td>
      <td>Towards a new heuristics to compute Consensus Ranking of Big Biological data</td>
      <td></td>
    </tr>
    <tr>
      <td></td>
      <td>P. ANDRIEU, L. BULTEAU, S. COHEN-BOULAKIA, A. DENISE, A. LABARRE, A. PIERROT et S. VIALETTE</td>
      <td>193</td>
    </tr>
    <tr>
      <td>Poster A89</td>
      <td>Comparaison de pipelines pour la découverte de signatures métagénomiques à partir de séquençage 16S en contexte clinique</td>
      <td></td>
    </tr>
    <tr>
      <td></td>
      <td>C. AZIZA, B. EMMA, S. IRADJ et M. DENIS</td>
      <td>194</td>
    </tr>
    <tr>
      <td>Poster A90</td>
      <td>Horizontal gene transfer from viruses in the genomes of plant-parasitic nematodes</td>
      <td></td>
    </tr>
    <tr>
      <td></td>
      <td>C. BELLIARDO, C. RANCUREL, E. DANCHIN et M. BAILLY-BECHET</td>
      <td>195</td>
    </tr>
    <tr>
      <td>Poster A91</td>
      <td>Mutation of Tyr137 of the universal Escherichia coli fimbrial adhesin FimH relaxes the tyrosine gate prior to mannose binding</td>
      <td></td>
    </tr>
    <tr>
      <td></td>
      <td>E.-M. KRAMMER, G. ROOS, M. PRÉVOST, J. BOUCKAERT et M. LENSINK</td>
      <td>196</td>
    </tr>
    <tr>
      <td>Poster A92</td>
      <td>Collapsing reads while maintaining qualities : srnaCollapser</td>
      <td></td>
    </tr>
    <tr>
      <td></td>
      <td>W. BEN SAOUD BENJERRI, C. GASPIN et M. ZYTNICKI</td>
      <td>197</td>
    </tr>
  </tbody>
</table>
## Posters partenaires 198

### Poster B1

BioMAnTM : a user-friendly interface for targeted metagenomic data visualization and analysis

P. VAISSIÉ, C. CAMUS, Y. AMOUZOU, T. CARTON, S. LE FRESNE-LANGUILLE, F. LE VACON, M. CAZAUBIEL et S. LEUILLET

199

### Poster B2

WHORMSS : un nouvel outil pour l'exploration taxonomique et fonctionnelle sans a priori des métagénomes

Y. LAURENT, J. DENONFOUX, A. BODEIN et S. FERREIRA

200

### Poster B3

CRISPR LifePipe : tools for the design of gRNAs and donor sequence required for genome editing using CRISPR/Cas9 system

V. CHESNAIS, E. CHAPLAIS, A. OTT et E. GINOUX

201

### Poster B4

Bioinfo-fr.net : présentation du blog communautaire scientifique francophone par lesEEKUS biologicus

G. BIOLOGICUS et I. STÉVANT

202

## Liste des contributeurs 203
Présentations orales
<figure><img src="image_3.png" /></figure>

Conférence invitée

Patrick WINCKER

CEA Genoscope, Evry, France

Holistic metagenomics in marine plankton communities
# Simka: large scale *de novo* comparative metagenomics

Gaëtan BENOIT¹, Pierre PETERLONGO¹, Mahendra MARIADASSOU², Erwan DREZEN¹,³, Sophie SCHBATH², Dominique LAVENIER¹ and Claire LEMAITRE¹

¹ INRIA/IRISA, Genscale team, UMR6074 IRISA CNRS/INRIA/Université de Rennes 1, Campus de Beaulieu, 35042, Rennes, France

² MaIAGE, INRA, Université Paris-Saclay, 78350, Jouy-en-Josas, France

³ CHU Pontchaillou, 35000, Rennes, France

Corresponding author: gaetan.benoit@inria.fr

**Reference paper:** Benoit et al. (2016) Multiple comparative metagenomics using multiset k-mer counting. *PeerJ Computer Science*. https://doi.org/10.7717/peerj-cs.94

Large metagenomic projects, such as Human Microbiome Project [1] or Tara Ocean Project [2], are becoming increasingly important for understanding life processes at the individual scale or at more global scales. Therefore, serious efforts are put in funding projects, collecting DNA, sequencing, and analyzing the resulting sequences. Comparative metagenomics is one of the most ubiquitous and informative of those analyzes. The purpose is mainly to estimate proximity (or distance) between two or more environmental sites at the genomic level. The comparisons are often based on identified species content. However, this approach is limited to sequences correctly assigned to known species documented in public biobanks, and this may correspond to small fractions of the datasets, in particular for environmental samples such as sea water. This limitation motivated the development of *de novo* comparison tools, such as Compareads [3] or Mash [4], based only on the non assembled read set comparisons. Compareads was for instance successfully used for the first analyses of the Tara Ocean data [5].

These reference-free methods share the use of *k*-mers as the fundamental unit used for comparing samples. Actually, *k*-mers are a natural unit for comparing communities: (1) sufficiently long *k*-mers are usually specific of a genome [6], (2) *k*-mer frequency is linearly related to genome's abundance [7], (3) *k*-mer aggregates organisms with very similar *k*-mer composition (e.g. related strains from the same bacterial species) without need for a classification of those organisms [8]. However, even if Compareads approach was designed to scale-up to large metagenomic read sets, its use on data generated by large scale projects is turning into a bottleneck in terms of time requirements. By contrast, Mash outperforms by far all other methods in terms of computational resource usage. However, this frugality comes at the expense of result quality and precision: the output distances and Jaccard indexes do not take into account relative abundance information and are not computed exactly due to *k*-mer sub-sampling. This is what motivated this work in which we propose a new *de novo* comparative metagenomic method, called Simka. Simka compares *N* metagenomic datasets based on their exact *k*-mers counts. It computes a large collection of distances classically used in ecology to compare communities, by replacing species counts by *k*-mer counts, for a large range of *k*-mer sizes, including large ones (up to 30). Simka is, to our knowledge, the first method able to rapidly compute a full range of distances enabling the comparison of any number of datasets. Simka outperforms state-of-the-art read comparison methods in terms of computational needs and result quality. For instance, Simka ran on 690 samples from the Human Microbiome Project (HMP) (totalling 32 billion reads) in less than 10 hours and using no more than 70 GB RAM.

Simka works as follows. Firstly, the *k*-mer spectrum of each dataset is computed. The *k*-mer spectrum of a dataset is the set of all its distinct *k*-mers associated with their abundance in the dataset. Secondly, *k*-mer spectra are compared in a pairwise manner to compute their distance. This comparison process basically aims at identifying which *k*-mers are shared by both spectra and which ones are not. It can be computationally very expensive because each *k*-mer spectrum can contains millions to billions of distinct *k*-mers when *k* is large (> 15). Moreover, the number of comparisons grows quadratically with the number of input datasets. To tackle this issue, we have designed a new *k*-mer counting strategy of numerous datasets, called Multiset *k*-mer Counting (MKC). MKC takes *N* datasets as input and provides an abundance vector for each distinct *k*-mer. The abundance vector of a *k*-mer consists of its *N* counts in the *N* datasets. The abundance vector generation by the MKC task is divided into two phases: (1) Sorting Count, (2) Merging Count. During the first step, the *k*-mers of each dataset are counted independently. This is performed by sorting the *k*-mers in lexicographical order. Distinct *k*-mers can thus be identified and their number of occurrences computed. This task can be very efficiently performed by popular disk-based *k*-mer counting tool such as DSK [9] or KMC2 [10]. The resulting *k*-mer spectra are written on the disk. During the second step, a Merge-Sort algorithm can be efficiently
applied on the sorted *k*-mer spectra to directly generate abundance vectors. Given those abundance vectors, the distances between each pair of datasets can be computed simultaneously. Interestingly, most of the ecological distances are additive over the distinct *k*-mers, meaning that they can be iteratively updated one abundance vector at a time. Once an abundance vector has been processed, there is thus no need to keep it on record, allowing Simka to have a very low memory footprint.

One advantage of the overall Simka workflow is its high parallelism potential. During the sorting count phase of the MKC, a first parallelism level is given by the independent counts of each dataset. *N* processes can thus be run in parallel, each one dealing with a specific dataset. A second level is given by the fine grained parallelism implemented in software such as DSK or KMC2 that intensively exploit today multicore processor capabilities. As the number of distinct *k*-mers is generally huge, those tools separate the *k*-mers in *P* smaller disjoint sets that can be counted independently and thus result in *P* *k*-mer spectrum chunks per dataset. During its second step, the MKC exploit this partitionning to merge up to *P* *k*-mer spectrums chucks in parallel. Each of these merge processes generates abundance vectors from which independent contribution to the distances are computed. Since the distance computed by Simka are additive over the distinct *k*-mers, each contribution is simply accumulated and the final distance is computed. Simka implementation is based on the GATB library [11], a C++ library optimized to handle very large sets of *k*-mers. Simka is usable on standard computers and has also been entirely parallelized for grid infrastructures made of hundred of nodes, and where each node implements 8 or 16-core systems.

The quality of the distances computed by Simka were evaluated answering two questions. First, are they similar to distances between read sets computed using other de novo approaches? Second, do they recover the known biological structure of HMP samples? For the first evaluation, we show that Simka result are perfectly well correlated with Compareads results. We go further in this evaluation by showing that Simka results are highly correlated with costly but extremely accurate de novo comparison techniques relying on all-versus-all sequence alignment strategy. For the second evaluation, Simka distances were compared to taxonomic distances that are a traditional way of comparing metagenomic samples. Taxonomic distances are based on sequence assignation to taxons by mapping to reference databases. To compare Simka to such traditional reference-based methods, we used the HMP dataset. One advantage of this dataset is that it has been extensively studied, in particular the microbial communities are relatively well represented in reference databases [1,12]. We show that substituting *k*-mer counts by species counts gives admittedly different distances but that those distances are biologically relevant as they capture the same underlying biological structure and lead to the same conclusions as those based on taxonomic composition. In particular, Simka was able to retrieve two major biological results. The first one is the segregation of the HMP datasets by body sites. The second one reveals that the organisation of the gut samples is mainly driven by the relative abundances of three bacterial genera, known as enterotypes, and characterized by the relative abundances of a few genera: *Bacteroides*, *Prevotella* and genera from the *Ruminococcaceae* family. In contrast of Simka, Mash performed badly when considering HMP datasets per body site since this tool can only take into account presence/absence information and not relative abundances. As a matter of fact, differences in relative abundances are subtler signals that are often at the heart of interesting biological insights in comparative genomics studies [13, 14, 15, 16, 17].

We introduced Simka, a new method for computing a collection of ecological distances, between many large metagenomic datasets, based on their *k*-mer composition. This was made possible thanks to the Multiset *k*-mer Counting algorithm (MKC), a new strategy that counts *k*-mers of numerous datasets with state-of-the-art time, memory and disk performances.

## Acknowledgements

This work was supported by the French ANR-14-CE23-0001 Hydrogen Project.

## References

[1] Human Microbiome Project Consortium. Structure, function and diversity of the healthy human microbiome. *Nature*, 486(7402):207–214, 2012.

[2] Eric Karsenti, Silvia G. Acinas, Peer Bork, Chris Bowler, Colomban de Vargas, Jeroen Raes, Matthew Sullivan, Detlev Arendt, Francesca Benzoni, Jean Michel Claverie, Mick Follows, Gaby Gorsky, Pascal Hingamp, Daniele Iudicone, Olivier Jaillon, Stefanie Kandels-Lewis, Uros Krzic, Fabrice Not, Hiroyuki Ogata, Stéphane Pesant, Em-
manuel Georges Reynaud, Christian Sardet, Michael E. Sieracki, Sabrina Speich, Didier Velayoudon, Jean Weisenbach, and Patrick Wincker. A holistic approach to marine Eco-systems biology. *PLoS Biology*, 9, 2011.

[3] Nicolas Maillet, Claire Lemaitre, Rayan Chikhi, Dominique Lavenier, and Pierre Peterlongo. Compareads: comparing huge metagenomic experiments. *BMC bioinformatics*, 13(Suppl 19):S10, 2012.

[4] Brian D. Ondov, Todd J. Treangen, Páll Melsted, Adam B. Mallonee, Nicholas H. Bergman, Sergey Koren, and Adam M. Phillippy. Mash: fast genome and metagenome distance estimation using minhash. *Genome Biol*, 17(1):132, 2016.

[5] E. Villar, G. K. Farrant, M. Follows, L. Garczarek, S. Speich, S. Audic, L. Bittner, B. Blanke, J. R. Brum, C. Brunet, R. Casotti, A. Chase, J. R. Dolan, F. Ortenzio, J.-P. Gattuso, N. Grima, L. Guidi, C. N. Hill, O. Jahn, J.-L. Jamet, H. Le Goff, C. Lepoitivre, S. Malviya, E. Pelletier, J.-B. Romagnan, S. Roux, S. Santini, E. Scalco, S. M. Schwenck, A. Tanaka, P. Testor, T. Vannier, F. Vincent, A. Zingone, C. Dimier, M. Picheral, S. Searson, S. Kandels-Lewis, S. G. Acinas, P. Bork, E. Boss, C. de Vargas, G. Gorsky, H. Ogata, S. Pesant, M. B. Sullivan, S. Sunagawa, P. Wincker, E. Karsenti, C. Bowler, F. Not, P. Hingamp, and D. Iudicone. Environmental characteristics of agulhas rings affect interocean plankton transport. *Science*, 348(6237):1261447–1261447, may 2015.

[6] Y. Fofanov, Y. Luo, C. Katili, J. Wang, Y. Belosludtsev, T. Powdrill, C. Belapurkar, V. Fofanov, T.-B. Li, S. Chumakov, and B. M. Pettitt. How independent are the appearances of n-mers in different genomes? *Bioinformatics*, 20(15):2421–2428, apr 2004.

[7] Yu-Wei Wu and Yuzhen Ye. A novel abundance-based algorithm for binning metagenomic sequences using l-tuples. *Journal of Computational Biology*, 18(3):523–534, 2011.

[8] Hanno Teeling, Jost Waldmann, Thierry Lombardot, Margarete Bauer, and Frank O Glöckner. Tetra: a web-service and a stand-alone program for the analysis and comparison of tetranucleotide usage patterns in dna sequences. *BMC bioinformatics*, 5(1):163, 2004.

[9] Guillaume Rizk, Dominique Lavenier, and Rayan Chikhi. DSK: k-mer counting with very low memory usage. *Bioinformatics*, page btt020, 2013.

[10] Sebastian Deorowicz, Marek Kokot, Szymon Grabowski, and Agnieszka Debudaj-Grabysz. KMC 2: Fast and resource-frugal k-mer counting. *Bioinformatics*, 31(10):1569–1576, 2015.

[11] Erwan Drezen, Guillaume Rizk, Rayan Chikhi, Charles Deltel, Claire Lemaitre, Pierre Peterlongo, and Dominique Lavenier. GATB: Genome assembly & analysis tool box. *Bioinformatics*, 30(20):2959–2961, 2014.

[12] Human Microbiome Project Consortium. A framework for human microbiome research. *Nature*, 486(7402):215–221, Jun 2012.

[13] Sébastien Boutin, Simon Y. Graeber, Michael Weitnauer, Jessica Panitz, Mirjam Stahl, Diana Clausznitzer, Lars Kaderali, Gisli Einarsson, Michael M. Tunney, J. Stuart Elborn, Marcus A. Mall, and Alexander H. Dalpke. Comparison of microbiomes from different niches of upper and lower airways in children and adolescents with cystic fibrosis. *PLoS ONE*, 10(1):1–19, 01 2015.

[14] A. Shade, S. E. Jones, J. G. Caporaso, J. Handelsman, R. Knight, N. Fierer, and J. A. Gilbert. Conditionally rare taxa disproportionately contribute to temporal changes in microbial diversity. *mBio*, 5(4):e01371–14–e01371–14, jul 2014.

[15] S. Genitsaris, S. Monchy, E. Viscogliosi, T. Sime-Ngando, S. Ferreira, and U. Christaki. Seasonal variations of marine protist community structure based on taxon-specific traits using the eastern english channel as a model coastal system. *FEMS Microbiology Ecology*, 91(5):fiv034–fiv034, mar 2015.

[16] Suzanne Coveley, Mostafa S Elshahed, and Noha H Youssef. Response of the rare biosphere to environmental stressors in a highly diverse ecosystem (zodletone spring, ok, usa). *PeerJ*, 3:e1182, 2015.

[17] V. Gomez-Alvarez, S. Pfaller, J. G. Pressman, D. G. Wahman, and R. P. Revetta. Resilience of microbial communities in a simulated drinking water distribution system subjected to disturbances: role of conditionally rare taxa and potential implications for antibiotic-resistant bacteria. *Environ. Sci.: Water Res. Technol.*, 2(4):645–657, 2016.
# Food-Microbiomes Transfert, a shotgun metagenomic tool and a database to analyze cheese ecosystems

Thibaut GUIRIMAND<fnref n="1" />, Anne-Laure ABRAHAM<fnref n="2" />, Sandra DEROZIER<fnref n="2" />, Charlie PAUVERT<fnref n="3" />, Mahendra MARIADASSOU<fnref n="2" />, Valentin LOUX<fnref n="2" />, and Pierre RENAULT<fnref n="1" />

Corresponding Author: thibaut.guirimand@inra.fr

**Abstract** The manufacturing process of cheeses, as for most fermented food, involves a complex flora, composed of bacteria but also yeast and filamentous fungi. These organisms can be brought by starters added during cheese manufacturing, or by environment (milk, maturing cellar...). The exact composition of cheeses is not known. Both academic researchers and cheese manufacturers are interested to have a better insight of cheese ecosystems, and collaborate in the Food-Microbiomes Transfert project. One of the objectives was to develop a metagenomic approach with a user-friendly tool, adapted to cheese ecosystems.

**Keywords** Metagenomics, next generation sequencing, microbial genomes, Database, Web server

## 1 Cheese ecosystems characterization with a metagenomic approach

The manufacturing process of cheeses, as for most fermented food, involves a complex flora, composed of bacteria but also yeast and filamentous fungi. The wide range of final products found on the dairy market is representative of the diversity of natural starters and ripening cultures used by dairy industries or coming from the food chain, from milk to the factory. However the cheese ecosystem is not completely understood [1]. The natural starters are not constructed from pure strains and the knowledge of their exact composition remains incomplete. Classical microbiological analysis or genetic methods (qPCR, MLST ...) can be used to better understand cheese ecosystem, but these techniques are expensive and time consuming. In order to further understand cheeses ecosystems and maintain a constant quality of cheese products, there is a need for a method to characterize low abundant species and assign precisely taxonomy, in cheese samples.

Techniques based on metagenomic DNA sequencing have been developed recently to rapidly identify species in complex ecosystems. Several tools are available to manage shotgun sequencing metagenomic datasets. Some are based on marker genes (for example: MetaPhlaAn [2], MetaPhyler [3], mOTU [4]) and propose rapid approach to identify species, although the use of a small part of the genomes decreases sensitivity of these approaches. Others use different strategies to take into account all the reads, for example with k-mer approaches (CLARK [5], Kraken [6], LMAT [7], OneCodex [8]...), read mapping on reference genomes (Genometa [9], GOTTCHA [10], MEGAN [11], MicrobeGPS [12], Sigma [13]...), assembly and functional annotation (EBI metagenomic web server [14], MG-RAST [15]...). Very few are available for biologists.

We are working in partnership with dairy manufacturers to develop a metagenomic approach based on shotgun sequencing of the cheese samples adapted to cheese ecosystems. Cheese ecosystem contain a reduced number of species (less than one hundred in most of the cases) and lots of reference genomes have been sequenced. However, there is a need to assign taxonomy of the present organisms, up to the strain level if possible, and to identify low abundance species. As different strains can have different property on the cheese manufacturing, it is important to have a tool able to identify genes present in the ecosystems.

We have developed a method based on the mapping of metagenomic reads (first 35 nucleotides) on a set of reference genomes with Bowtie [16], with 3 mismatches allowed. These parameters have been chosen to allow

<footnote n="1">¹ Micalis Institute, INRA, AgroParisTech, Université Paris-Saclay, 78350 Jouy-en-Josas, France</footnote>

<footnote n="2">² MaIAGE, INRA, Université Paris-Saclay, 78350 Jouy-en-Josas, France</footnote>

<footnote n="3">³ BIOGECO, INRA, Université de Bordeaux, 69 Route d'Arcachon, 33612 Cestas, France</footnote>
detecting microorganism displaying up to 10% divergence with reference genomes. A first filter based on read distribution on CDS allows discarding false positives. We then propose a method to identify among reference genomes of the same species those that are the closest to ecosystem micro-organisms. It is based on the number of CDS detected on each genome, a comparison between observed and expected CDS coverage [17], and informations about mismatches in the alignments. We have implemented this analysis method under a python pipeline named GeDI.

## 2 Food-Microbiomes Transfert, a specific database and an interface to analyze cheese ecosystems

Food-Microbiomes Transfert aims to provide a user-friendly tool to analyze cheeses ecosystems. To create the most accessible tool, the project offers a web interface to GeDI and a cheese specific genomes database. This interface will allow users to analyze their own metagenomes (or public metagenomes).

The genomes database has been created using PostgreSQL and currently contains 99 cheese specific microorganisms (we expect soon 300 genomes, with at least one genome for all the genus that have been described in dairy products). These genomes have been extracted from public databases by dairy products ecosystems experts and will be enriched with ecological metadata using text-mining tools and the Ontobiotope ontology [18]. In addition to these public genomes, the user is able to add his own private genomes to perform analysis.

A metagenome database allows to store metagenomics raw data and GeDI results. For this purpose, a metadata model representing cheese manufacturing, sampling method and cheese classification has been created, in partnership with cheese manufacturers and researchers to identify the most accurate and accessible model.

These two parts of the database have been conceived using the Minimum Information about a Genome/Metagenome Sequence (MIGS and MIMS [19]). It allows to have a standard and reusable set of data.

The client side is developed using JavaScript/HTML5/CSS3/RDFa and interacts with the server using AJAX queries. The aim is to create a dynamic interface with minimal user interaction needs and an easy way to perform/manage analysis and data.

The user can upload data, share them with other users, manage genomes lists to reuse for a later analysis, and perform analysis with a minimal steps amount: i) select the metagenome ii) select genomes from public or private lists iii) start GeDI with default parameters.

A results page allows the user to visualize the mapping summary and to download the results (mapping tables, charts as shown in Fig 1).

<figure><img src="image_4.png" /><figcaption>**Fig 1.** Here is an example of graphics generated by the tool by mapping a 10 million reads metagenome on *Streptococcus thermophilus* JIM 8232. We can observe some uncovered CDS marking differences between the reference genome and the strain present in the ecosystem.</figcaption></figure>
The server, hosted by the Migale platform, is based on two specific technologies. The web server uses the Python Django framework to manage web client requests, databases and users. The GeDI computation is done on a cluster using a Galaxy [20] instance called by the Python Bioblend library [21]. The use of Galaxy facilitate the reproducibility of research because of the possibility of exporting the histories and tools. It also allows to easily link the web server to the analysis pipeline because of the use of the same language: Python. The actual computation time is about 10 hours for a 1 million reads metagenome and a hundred of reference genomes and about 5 days for a 10 million reads metagenome.

## 3 Prospects

We are working on the improvement of GeDI tool: validation on several datasets, computation time... The genome database will be enriched with new genomes and expert annotations especially with text-mining tools. The metadata of the metagenomic database will be added. We are also working on the interface improvement in order to make analysis even more intuitive, and to provide tools to perform cross-comparisons between analyses. GeDI tool will be soon available with command line and Galaxy wrappers. The database and interface will be available via the Migale platform.

## Acknowledgements

We are grateful to the INRA Migale bioinformatics platform (http://migale.jouy.inra.fr) for providing help and support. We also would like to thanks Claire Nédellec, Robert Bossy and Estelle Chaix from the INRA Bibliome team for their work and their support to use Ontobiotope.

## References

[1] Almeida M, Hébert A, Abraham AL, Rasmussen S, Monnet C, Pons N, Delbès C, Loux V, Batto JM, Leonard P, Kennedy S. Construction of a dairy microbial genome catalog opens new perspectives for the metagenomic analysis of dairy fermented products. BMC genomics, 13:15(1):1101, 2014 Dec.

[2] Segata N, Waldron L, Ballarini A, Narasimhan V, Jousson O, Huttenhower C. Metagenomic microbial community profiling using unique clade-specific marker genes. Nature methods, 1:9(8):811-4, 2012 Aug.

[3] Liu B, Gibbons T, Ghodsi M, Pop M. MetaPhyler: Taxonomic profiling for metagenomic sequences. InBioinformatics and Biomedicine (BIBM), pages 95-100, 2010 IEEE International Conference, 2010 Dec 18.

[4] Sunagawa S, Mende DR, Zeller G, Izquierdo-Carrasco F, Berger SA, Kultima JR, Coelho LP, Arumugam M, Tap J, Nielsen HB, Rasmussen S. Metagenomic species profiling using universal phylogenetic marker genes. Nature methods, Dec 1;10(12):1196-9, 2013.

[5] Ounit R, Wanamaker S, Close TJ, Lonardi S. CLARK: fast and accurate classification of metagenomic and genomic sequences using discriminative k-mers. BMC genomics, 25:16(1):236, 2015 Mar.

[6] Wood DE, Salzberg SL. Kraken: ultrafast metagenomic sequence classification using exact alignments. Genome biology, Mar 3;15(3):R46, 2014.

[7] Ames SK, Hysom DA, Gardner SN, Lloyd GS, Gokhale MB, Allen JE. Scalable metagenomic taxonomy classification using a reference genome database. Bioinformatics, 15:29(18):2253-60, 2013 Sep.

[8] Minot SS, Krumm N, Greenfield NB. One Codex: A Sensitive and Accurate Data Platform for Genomic Microbial Identification. bioRxiv, 1:027607, 2015 Jan.

[9] Davenport CF, Neugebauer J, Beckmann N, Friedrich B, Kameri B, Kokott S, Paetow M, Siekmann B, Wieding-Drewes M, Wienhöfer M, Wolf S. Genometa-a fast and accurate classifier for short metagenomic shotgun reads. PloS one, 21;7(8):e41224, 2012 Aug.

[10] Freitas TA, Li PE, Scholz MB, Chain PS. Accurate read-based metagenome characterization using a hierarchical suite of unique signatures. Nucleic acids research, 12:gkv180, 2015 Mar.

[11] Huson DH, Auch AF, Qi J, Schuster SC. MEGAN analysis of metagenomic data. Genome research, 1;17(3):377-86, 2007 Mar.

[12] Lindner MS, Renard BY. Metagenomic profiling of known and unknown microbes with MicrobeGPS. PloS one, 2;10(2):e0117711, 2015 Feb.

[13] Ahn TH, Chai J, Pan C. Sigma: Strain-level inference of genomes from metagenomic analysis for biosurveillance. Bioinformatics, 29:btu641, 2014 Sep.

[14] Hunter S, Corbett M, Denise H, Fraser M, Gonzalez-Beltran A, Hunter C, Jones P, Leinonen R, McAnulla C, Maguire E, Maslen J. EBI metagenomics—a new resource for the analysis and archiving of metagenomic data. Nucleic acids research, 1;42(D1):D600-6, 2014 Jan.

[15] Meyer F, Paarmann D, D'Souza M, Olson R, Glass EM, Kubal M, Paczian T, Rodriguez A, Stevens R, Wilke A, Wilkening J. The metagenomics RAST server–a public resource for the automatic phylogenetic and functional analysis of metagenomes. BMC bioinformatics, 19;9(1):386, 2008 Sep.
[16] Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology, 2009 Mar.

[17] Lander ES, Waterman MS. Genomic mapping by fingerprinting random clones: A mathematical analysis. Genomics, 2, 3, 231-239, 1988.

[18] Bossy R, Golik W, Ratkovic Z, Bessières P, Nédellec C. Bionlp shared task 2013–an overview of the bacteria biotope task. In: Proceedings of the BioNLP Shared Task pages 161-169, 2013 Workshop, 2013 Aug.

[19] Yilmaz P, Kottmann R, Field D, Knight R, Cole JR, Amaral-Zettler L, Gilbert JA, Karsch-Mizrachi I, Johnston A, Cochrane G, Vaughan R. Minimum information about a marker gene sequence (MIMARKS) and minimum information about any (x) sequence (M1xS) specifications. Nature biotechnology, 1:29(5):415-20, 2011 May.

[20] Goecks J, Nekrutenko A, Taylor J. Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome biology, 25:11(8):R86, 2010 Aug.

[21] Sloggett C, Goonasekera N, Afgan E. BioBlend: automating pipeline analyses within Galaxy and CloudMan. Bioinformatics, 1:29(13):1685-6, 2013 Jul.
# Reconstruction of full-length 16S rRNA sequences for taxonomic assignment in metagenomics

Pierre PERICARD¹,², Yoann DUFRESNE¹,², Samuel BLANQUART²,¹ and Hélène TOUZET¹,²

¹ CRISTAL (UMR CNRS 9189, Université Lille 1), 59655 Villeneuve d'Ascq Cedex, France

² Inria Lille Nord Europe, 59650, Villeneuve d'Ascq, France

Corresponding author: pierre.pericard@gmail.com, helene.touzet@univ-lille1.fr

**Abstract** Advances in the sequencing of uncultured environmental samples, raise a growing need for accurate taxonomic assignment. Accurate identification of organisms present within a community is essential to understanding even the most elementary ecosystems. However, current high-throughput sequencing technologies generate short reads which partially cover full-length marker genes and this poses difficult bioinformatic challenges for taxonomy identification at high resolution. We designed MATAM, a software dedicated to the fast and accurate targeted assembly of short reads sequenced from a genomic marker of interest. The method implements a stepwise process based on construction and analysis of a read overlap graph. It is applied to the assembly of 16S rRNA markers and is validated on simulated, synthetic and genuine metagenomes. We show that MATAM outperforms other available methods in terms of low error rates and recovered genome fractions and is suitable to provide improved assemblies for precise taxonomic assignments.

**Keywords** Metagenomics, 16S rRNA, assembly, taxonomic assignment, algorithm.

## 1 Introduction

Shotgun metagenomic sequencing provides an unprecedented opportunity to study uncultured microbial samples, with multiple applications ranging from the human microbiome to soil or marine samples, for which the vast majority of microorganisms diversity remains unknown [1].

A major goal of metagenomic studies is to characterize the microbial diversity and ecological structure. This is often achieved by focusing on one of several phylogenetic marker genes [2,3], that are ubiquitous in the taxonomic range of interest and exhibit variable discriminative regions. For bacterial communities, the gold standard marker is the 16S ribosomal RNA (rRNA, ~1500bp avg. length), for which millions of sequences are available in curated reference databases, such as Silva [4], RDP [5] or GreenGenes [6]. Traditional approaches such as amplicon sequencing are limited to the analysis of small portions of the marker sequences. This leads to strong technological limitations for organisms identification at sufficiently precise taxonomic levels, typically beyond genus [7]. To assign marker sequences to species, or even strains, we need to be able to recover full length rRNA with less than a few errors per kilobase. Metagenomic assemblers are not suitable for this task, because they are optimized to deal with whole genomes, and struggle to differentiate between very similar sequences [8]. To this respect, marker-oriented methods such as EMIRGE [9] and REAGO [10] were recently developed in order to assemble metagenomic read subsets into full length 16S rRNA contigs, thus aiming to improve the taxonomic assignment accuracy of environmental samples. EMIRGE uses a Bayesian approach to iteratively reconstruct 16S rRNA full length sequences. REAGO identifies rRNA reads using Infernal [11], and then constructs an overlap graph by searching for exact overlaps between reads using a suffix/prefix array. However, such tools still show some limitations in terms of recovery error rates as well as dealing with low abundance species.

In this work, we present MATAM, a new approach based on the construction and exploitation of an overlap graph, carefully designed to minimize the error rate and the risk of chimera formation. MATAM was validated on both simulated and actual sequencing data. It is able to reconstruct nearly full length 16S rRNAs and is robust to variations in the sequencing depth as well as community complexity.

## 2 Methods

### 2.1 Overview of MATAM

The MATAM (Mapping-Assisted Targeted-Assembly for Metagenomics) pipeline takes as input a set of shotgun metagenomics short reads and a reference database containing the largest possible set of sequences
from a given target marker gene. MATAM identifies reads originating from that marker, and assembles nearly full length sequences of it. It is composed of four major steps. Although this method should work for any conserved and widely surveyed gene, we will focus on the 16S rRNA for the remainder of the article.

## 2.2 Reference database construction

The availability of a reference database for the marker gene is an essential feature of the method, because it allows us to model the target sequences. For applications to 16S rRNA assembly, MATAM utilizes Silva 128 SSU Ref NR database [4]. From this reference database that we denote as *complete*, we also build a *clustered* reference database, that provides a coarse-grained representation of the taxonomic space. For that task, we use *Sumaclust* [12,13] (http://metabarcoding.org/sumaclust) using a 95% identity threshold.

## 2.3 rRNA reads identification and mapping

In the first step, reads are mapped against the clustered reference database using SortMeRNA [14,15]. This step allows to quickly sort out 16S rRNA reads from the whole set of reads, providing high quality alignments. For each read, we keep up to ten best alignments against the reference database. Moreover, this mapping step yields a broad classification of the 16S rRNA reads. Indeed, reads coming from distantly related species are aligned against their respective closest known references, which nest in distant lineages of the taxonomy, while reads from closely related species are aligned against closely related references.

## 2.4 Construction of the overlap graph

The identified 16S rRNA reads are then organized into an *overlap graph* defined as follows: graph nodes are reads, and an undirected edge connects two nodes if the two reads overlap with a sufficient length and with a sufficient identity to assert that they originated from a common sampled taxon. The standard approach to build such an overlap graph requires comparison of each read with each other, which is time-consuming. Here, we use alignment information to sort through candidate read pairs in a very efficient manner. For each pairing, we consider only reads that share alignments with at least one common reference sequence and for which the alignments are overlapping on more than 50 nucleotides with 100% identity. This strict criterion allows us to reduce the risk of connecting reads from unrelated taxa, which would in turn produce chimeras. By doing so, we discard reads containing sequencing errors in their overlap, which is bearable considering the nowadays very low sequencing error rates of short reads.

## 2.5 Extracting contigs from the overlap graph

Although the *overlap graph* appears very bushy, it also reveals some general trends. While it exhibits highly connected subgraphs, it also displays disjoint paths. We simplify the graph by performing a breadth first traversal starting from a random node to annotate the nodes with their depth. All nodes with equal depth that are connected in a single connected component are collapsed into a single *compressed node* and outgoing edges are merged into a *compressed edge*. Low support compressed nodes containing a single read, and compressed edges representing a single overlap are removed. The resulting graph, called the *compressed graph*, is several order of magnitude smaller than the initial *overlap graph*. We partition this graph in three categories of subgraphs: *hubs*, that are nodes with an degree strictly greater than two, *specific paths* that are sequences of nodes of degree two or one, and *singletons* that are non-connected nodes. Intuitively, hubs correspond to the highly connected subgraphs in the *overlap graph*, and are likely to contain mainly reads coming from conserved regions shared in many species, thus overlapping without error even for distantly related taxa. Specific paths tend to contain reads originating from variable regions of the 16S gene, that are specific to one or few closely related species. For each subgraph in the *compressed graph* (hubs, specific paths, singletons), we extract the underlying sets of reads and build an individual assembly using the genomic assembler SGA [16]. Note that any other state-of-the art genomic assembler could be used here. As a result, we obtain one or more contigs for each subgraph.

## 2.6 Contigs scaffolding

We use a greedy algorithm to scaffold the contigs obtained in the previous step. For that task, contigs are first mapped against the complete reference database, and all alignments within the 1% range of suboptimal scores are kept. We then select contigs by increasing number of matches and decreasing lengths. By doing so, a long contig with a unique alignment will be selected for scaffolding before a short contig exhibiting a large number of alignments. Such long contig can be assigned non-ambiguously to a single species, while the short contig with multiple matches rather corresponds to a conserved region of the marker and is used to fill in the blanks between the specific contigs. Contigs matching against the same reference sequence are then
<table><caption>**Tab. 1.** Results for the simulated dataset with varying sequencing depth. We provide averaged metrics for the five sequencing depths. ACL is the average contig length.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Chimera (%)</th><th colspan="2">TAL/TL (%)</th><th colspan="2">ER (%)</th><th colspan="2">Ns (%)</th><th colspan="2">ACL</th></tr><tr><th>mean</th><th>SD</th><th>mean</th><th>SD</th><th>mean</th><th>SD</th><th>mean</th><th>SD</th><th>mean</th><th>SD</th></tr></thead><tbody><tr><td>MATAM</td><td>1.28</td><td>0.55</td><td>99.3</td><td>0.2</td><td>0.03</td><td>0.02</td><td>0.00</td><td>0.00</td><td>1252</td><td>116.9</td></tr><tr><td>EMIRGE</td><td>36.89</td><td>9.42</td><td>79.9</td><td>11.6</td><td>0.62</td><td>0.16</td><td>0.55</td><td>0.36</td><td>1436</td><td>15.4</td></tr><tr><td>REAGO</td><td>42.11</td><td>10.36</td><td>91.5</td><td>0.8</td><td>0.31</td><td>0.13</td><td>0.00</td><td>0.00</td><td>1333</td><td>298.9</td></tr><tr><td>SPAdes</td><td>21.23</td><td>9.05</td><td>73.5</td><td>15.9</td><td>0.60</td><td>0.49</td><td>0.02</td><td>0.04</td><td>966</td><td>47.4</td></tr><tr><td>MEGAHIT</td><td>23.81</td><td>2.85</td><td>80.3</td><td>4.9</td><td>0.36</td><td>0.18</td><td>0.00</td><td>0.00</td><td>962</td><td>87.6</td></tr></tbody></table>

merged into a single consensus scaffold. Redundant scaffolds included in larger ones are removed. Finally, only scaffolds larger than 500bp are retained. This yields the final MATAM output which could be used for the purpose of taxonomic assignment.

## 3 Implementation

MATAM was implemented in Python 3, except for the overlap graph building and compression steps that were written in C++11 using the SeqAn library [17], and is available *via* Docker and Conda. MATAM is distributed under the GNU Affero GPL v3.0 licence and the source code is freely available at the following URL: https://github.com/bonsai-team/matam. All MATAM runs presented in this article were performed using MATAM v0.9.9.

## 4 Results

MATAM performance was compared with those of two general-purpose metagenomic assemblers, SPAdes [18 and MEGAHIT [19], as well as with two methods specialized in 16S rRNA assembly, EMIRGE [9] and REAGO [10]. The five tools were run on three different datasets, chosen for their complementarity and the possibility to validate the reconstructed candidate 16S rRNA sequences: a simulated dataset [20], a synthetic microbial community [21], and two environmental samples from human gut and mouth providing amplicon based taxonomic assignments [22]. SortMeRNA was used to extract 16S rRNA reads from these datasets before assembling them with SPAdes and MEGAHIT. Complete command-lines and parameters are available in the Supplementary Results.

In order to compare the five methods on a common ground, the same validation procedure was applied for all experiments. Only reconstructed sequences with lengths exceeding 500bp were considered, and chimeric sequences were filtered out by the UCHIME algorithm [23] implemented in VSEARCH [24] and querying the Silva 128 SSU Ref Nr99 database. For each experiment, we indicate the proportion of chimeric contigs (% chimeras, which is the total size of all chimeric contigs divided by the assembly total size). All the following measures were then computed on the remaining assemblies. When the sequences present in the sample are actually known (see Sections 4.1 and 4.2), the assembly quality assessment was performed with MetaQuast [25] by aligning the contigs against the original sample sequences, and considering the following metrics: the number of contigs (#contigs), which is the total number of contigs of lengths greater than 500bp; the total length (TL), which is the total number of bases in the contigs; the total aligned length (TAL), which is the total number of aligned nucleotides in the contigs; the genome fraction (GF), which stands for the total number of nucleotides from the original sample sequences covered with contigs divided by the total size of the sample sequences; the error rate (ER), which consists in the percentage of observed mismatches and indels with respect to the closest matched sequence in the original sample. Finally, taxonomic assignments were carried out with the RDP Classifier [26].

### 4.1 Simulated metagenomic datasets with varying sequencing depth

In the first experiment, we evaluated the ability of methods to correctly reconstruct the 16S rRNA sequences in the context of low sequencing depth. For that, we used a selection of 122 genomes providing a realistic taxomical diversity [20,27], that contains 287 distinct 16S rRNA copies. We generated five datasets with varying sequencing depths: 50x, 20x, 10x, 5x and 2x per genome. Illumina reads were simulated with the ART simulator [28], using the HiSeq2500 built-in error profile, 101bp read length, and 250bp fragment length with a 30bp standard-deviation (SD). In this simulation, all species are equally distributed, which corresponds to the high complexity community introduced in [20].
<figure><img src="image_5.png" /><figcaption>Fig. 1. Effect of sequencing depth on the assemblies genome fractions.</figcaption></figure>

<table><caption>Tab. 2. Results for the synthetic community.</caption><thead><tr><th></th><th>Chimera (%)</th><th>#contigs</th><th>TL</th><th>TAL</th><th>GF (%)</th><th>ER (%)</th><th>Ns (%)</th></tr></thead><tbody><tr><td>MATAM</td><td>3.2</td><td>101</td><td>139220</td><td>130654</td><td>83.1</td><td>0.05</td><td>0</td></tr><tr><td>EMIRGE</td><td>17.4</td><td>82</td><td>117138</td><td>102856</td><td>50.7</td><td>0.17</td><td>1.12</td></tr><tr><td>REAGO</td><td>15.5</td><td>59</td><td>90269</td><td>81297</td><td>42.8</td><td>0.06</td><td>0</td></tr><tr><td>SPAdes</td><td>5.5</td><td>59</td><td>70229</td><td>59988</td><td>39.9</td><td>0.11</td><td>0.05</td></tr><tr><td>MegaHit</td><td>3.0</td><td>61</td><td>77251</td><td>68904</td><td>44.3</td><td>0.18</td><td>0</td></tr></tbody></table>

Table 1 shows the results averaged over the five datasets (*mean* metrics and their respective standard deviation, SD). More than 99% of the MATAM sequences were aligned by MetaQuast to one of the 287 16S rRNA sequences from the initial sample (mean TAL/TL), while among other methods, this proportion reached at best 91%, with REAGO. Congruently, MATAM sequences obtained the lowest average error rate (ER=0.03%), which represents more than a ten-fold accuracy gain compared to the other assemblers, and a twenty-fold improvement over EMIRGE. Furthermore, EMIRGE sequences contained 0.5% of unknown nucleotides (Ns), bringing its effective ER above 1%. Additionally, MATAM recovered about thirty times less chimeras than REAGO and EMIRGE did.

For each of the five tools, we reported the recovered genome fraction (GF) with respect to increasing sequencing depth (Figure 1). MATAM recovered from 76% to 85% of the reference sequences for sequencing depths greater than 10x, while EMIRGE recovered less than 55% of the reference sequence, and the GF for other methods is lower than 22%. MATAM also achieved the best performance facing a low sequencing depth of 2x, reaching a GF of 33%, while GFs ranged between 5% and 10% with all other assemblers.

## 4.2 Synthetic archaeal and bacterial community

Inching toward more realistic applications, a second dataset provides Illumina reads extracted from a synthetic microbial community composed of 16 archaeal species from 12 genera, as well as 48 bacterial species from 36 genera (accession SRR606249; [21]). As emphasized by the authors, the selected organisms cover a wide range of environmental conditions and adaptation strategies. In contrast to the previous simulated dataset (Section 4.1), the proportion of each species in the sample is not uniform, which results in individual genome average sequencing depth varying from 9x to 318x. The number of 16S rRNA paralogs per genome appears also highly diverse, ranging from 1 to 10 copies per genome. Altogether, this dataset represents a total amount of 106 distinct 16S rRNA sequences with pairwise sequence identities ranging from 59.64% to 99.93%.

The organisms were sequenced on Illumina HighSeq2000, providing 109 million 101bp paired-end reads with an average fragment size of 250bp. We quality cleaned the reads using Prinseq Lite [29], removed adapter sequences using Cutadapt [30], filtered out short reads (<50bp), and obtained a total number of 67.6 million reads, which were analyzed with MATAM and EMIRGE. The uncleaned raw dataset was provided to REAGO, considering that the method could not handle reads with varying lengths. Finally, for SPAdes and MEGAHIT, the 16S rRNA reads were extracted from the cleaned dataset using SortMeRNA, which provided 108,560 16S rRNA reads to assemble.
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Alignment of the reference sequences with the assembled contigs shows MATAM ability to differentiate between very close sequences. MATAM, EMIRGE and REAGO contigs are shown respectively in blue, red and green. In an ideal setting, each software should produce contigs that cluster closely to each reference (black) sequence. Contigs followed by a star, and drawn in a darker color, were considered as chimeric by VSEARCH.</figcaption></figure>

Results are shown in Table 2. Confirming the trends observed on the simulated dataset, MATAM is able to recover the highest number of sequences together with the highest GF (83%). Most importantly, with lower ER than achieved by the other tested methods, the MATAM assembly appears highly accurate. While EMIRGE is the second best approach in terms of recovered GF, it also yields the greatest ER and Ns over all the compared tools. Moreover, a RDP classification of MATAM and EMIRGE sequences indicates that while MATAM missed one expected genus only, EMIRGE missed 4 genera out of 48.

Inspection of the MetaQuast alignments of the assemblies against the original 16S rRNAs revealed that all methods accurately assembled the genes sharing less than 90% sequence identity with their closest relatives within the sample. However, performances significantly dropped when attempting to assemble the closely related genes in the dataset. This especially concerned the paralogous 16S rRNA copies sharing around 99% sequence identity.

We selected sequences from a representative subset of four related species possessing one to three such paralogous copies. Those 16S rRNAs and their corresponding assembled candidate sequences were selected for a phylogenetic tree reconstruction. The obtained tree (Figure 2) demonstrates that MATAM correctly assembled all the different paralogs with nearly no error, while EMIRGE and REAGO only managed to recover one candidate sequence per species. Thus, EMIRGE and REAGO merged into a single candidate sequence the reads issued from distinct paralogs, resulting in erroneous assemblies with high ER and underestimated GF. Indeed, each of the sequences assembled with REAGO, as well as one EMIRGE sequence over four, appear to cluster at a slight distance from their respective targeted paralogs. Those distances simply account for the methods reconstruction errors. Consistently, in two cases, the candidates assembled by EMIRGE and REAGO were identified as chimeras by VSEARCH.

## 4.3 Human Microbiome Project

Finally, we used two metagenomic samples from the Human Microbiome Project (gut: SRS011405, and mouth: SRS016002, [22]) in order to validate MATAM on real metagenomic datasets sequenced from genuine environments. The reads were already quality cleaned and trimmed, and no additional filtering was performed. Hence, reads having different lengths, we were not able to run REAGO on these datasets. Results obtained
<table><caption>**Tab. 3.** Results for the gut and mouth HMP datasets. The column #classes indicates the total number of taxonomic classes found with RDP from the assemblies, with the number of these classes validated with the QIIME OTUs (in parentheses). The column #genera gives the same information at the genus level.</caption><thead><tr><th></th><th></th><th>Chimera (%)</th><th>#contigs</th><th>TL #classes</th><th>#genera</th></tr></thead><tbody><tr><td rowspan="2">SRS011405</td><td>MATAM</td><td>3.37%</td><td>218 187710</td><td>5 (4)</td><td>21 (17)</td></tr><tr><td>EMIRGE</td><td>43.04%</td><td>273 393152</td><td>2 (2)</td><td>12 (8)</td></tr><tr><td rowspan="2">SRS016002</td><td>MATAM</td><td>4.92%</td><td>353 320748</td><td>13 (13)</td><td>31 (28)</td></tr><tr><td>EMIRGE</td><td>46.01%</td><td>282 394087</td><td>12 (12)</td><td>25 (23)</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** % identity distribution of best matches against Silva 128 SSU Ref NR.</figcaption></figure>

with SPAdes and MEGAHIT using the following protocol appeared highly inaccurate and therefore, they are not further commented in this work. Thus, we only present the results obtained with EMIRGE and MATAM.

For these two datasets, the exact ground truth is unknown. Thus we could not perform the same validation procedure as in the two previous examples and we had to resort to alternative strategies. First, we took advantage of the availability of OTU sequences inferred through a QIIME analysis of the V1-V3 hypervariable regions for the same biological samples (available from the SRS accession numbers). We compared the assignments obtained from assemblies, calculated with RDP, with these of amplicon OTUs (Table 3). For both samples, MATAM identified more classes and genera than EMIRGE did, and most of these taxa were validated by the amplicon OTUs. Interestingly, we observed that in the two samples, three genera were recovered both by MATAM and EMIRGE, but not by the amplicon approach: *Odoribacter*, *Peptococcus*, and *Bergeyella*. Since some species from these genera are known to be adapted to the human gut and mouth environments, it is plausible that they were missed by the amplicon approach while being accurately recovered by MATAM and EMIRGE from the metagenomic samples.

Moreover, we evaluated assembly quality by aligning MATAM and EMIRGE sequences against the complete Silva 128 SSU Ref NR database, using BLAST. The rationale for this experiment is that most of the species in these human gut and mouth samples are possibly already known, and therefore should be found in Silva. We observed that nearly all MATAM sequences matched with a known 16S rRNA in Silva with more than 99% identity, among which a majority matched with 100% identity (Figures 3a and 3b), which suggests that MATAM sequences could possibly be assigned at the species or even the strain level. On the other hand, EMIRGE sequences provided a discordant picture. In the case of the human mouth sample, most of the EMIRGE sequences obtained a match above 97% identity, but only a slight proportion of them matched with 100% identity against a known 16S rRNA (Figure 3b). The observation is even more pronounced with the human gut sample, where only 43% of the EMIRGE sequences obtained a match above 97% identity against a Silva 16S rRNA sequence (Figure 3a). Thus, conversely to MATAM, EMIRGE sequences would suggest that only a slight proportion of the human gut and mouth diversity has a known isolate registered in Silva. However, considering our previous conclusions on controlled datasets, we assume that part of this diversity inferred with EMIRGE might in fact corresponds to reconstruction artifacts.
## 5 Discussion

Taxonomic assignments of environmental samples is a strikingly difficult task which suffers from inherent limitations of high-throughput sequencing technologies. In this respect, we designed MATAM as an alternative to existing software helping to better understand the taxonomic structures of shotgun metagenomic samples. Our experimental results show that MATAM outperforms other available tools providing phylogenetic marker assemblies. Reconstructing full length 16S rRNAs allows to reach a higher precision of taxonomic assignments than individual read analysis or amplicon sequencing do, because the reconstructed sequences effectively contain stronger phylogenetic signal. Moreover, metagenomic shotgun sequencing is naturally immune against the primer and amplification biases attached to the amplicon sequencing technology, and therefore is more adequate to sequence unknown species.

Our approach opens up several new perspectives. Although we have focused this work on the assembly of 16S rRNA genes, MATAM was designed to deal with any marker of taxonomic interest. Indeed, there is currently an emerging trend to consider a combination of universal (single-copy) marker families, such as provided in the recently published database proGenomes [31]. Sequences from this database, or from any other customized one, could be used with MATAM to target a variety of markers, and thus provide improving taxonomic assignments. MATAM could also be used in combination with other types of sequencing data. Long read sequencing is able to produce fragments that cover large regions of the DNA molecules, up to several thousands of bases. When long reads are available, they could serve as a guide in the scaffolding step of MATAM and concomitantly, MATAM low-error contigs could be used to correct them. Finally, targeted gene capture, that allows to sequence at high depth captured DNA regions of interest from an environmental sample [32], could also prove to be an exciting application field for MATAM.

## References

[1] Kenneth J. Locey and Jay T. Lennon. Scaling laws predict global microbial diversity. *Proceedings of the National Academy of Sciences*, 113(21):5970–5975, May 2016.

[2] Bo Liu, Theodore Gibbons, Mohammad Ghodsi, Todd Treangen, and Mihai Pop. Accurate and fast estimation of taxonomic profiles from metagenomic shotgun sequences. *BMC Genomics*, 12(Suppl 2):S4, July 2011.

[3] Nicola Segata, Levi Waldron, Annalisa Ballarini, Vagheesh Narasimhan, Olivier Jousson, and Curtis Huttenhower. Metagenomic microbial community profiling using unique clade-specific marker genes. *Nature Methods*, 9(8):811–814, 2012.

[4] Christian Quast, Elmar Pruesse, Pelin Yilmaz, Jan Gerken, Timmy Schweer, Pablo Yarza, Jörg Peplies, and Frank Oliver Glöckner. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. *Nucleic Acids Research*, 41(D1):D590–D596, January 2013.

[5] James R. Cole, Qiong Wang, Jordan A. Fish, Benli Chai, Donna M. McGarrell, Yanni Sun, C. Titus Brown, Andrea Porras-Alfaro, Cheryl R. Kuske, and James M. Tiedje. Ribosomal Database Project: data and tools for high throughput rRNA analysis. *Nucleic Acids Research*, 42(D1):D633–D642, January 2014.

[6] T. Z. DeSantis, P. Hugenholtz, N. Larsen, M. Rojas, E. L. Brodie, K. Keller, T. Huber, D. Dalevi, P. Hu, and G. L. Andersen. Greengenes, a Chimera-Checked 16s rRNA Gene Database and Workbench Compatible with ARB. *Applied and Environmental Microbiology*, 72(7):5069–5072, January 2006.

[7] Rachel Poretsky, Luis M. Rodriguez-R, Chengwei Luo, Despina Tsementzi, and Konstantinos T. Konstantinidis. Strengths and Limitations of 16s rRNA Gene Amplicon Sequencing in Revealing Temporal Microbial Community Dynamics. *PLoS ONE*, 9(4):e93827, April 2014.

[8] Alexander Sczyrba, Peter Hofmann, Peter Belmann, David Koslicki, Stefan Janssen, Johannes Droege, Ivan Gregor, Stephan Majda, Jessika Fiedler, Eik Dahms, Andreas Bremges, Adrian Fritz, Ruben Garrido-Oter, Tue Sparholt Jorgensen, Nicole Shapiro, Philip D. Blood, Alexey Gurevich, Yang Bai, Dmitrij Turaev, Matthew Z. DeMaere, Rayan Chikhi, Niranjan Nagarajan, Christopher Quince, Lars Hestbjerg Hansen, Soren J. Sorensen, Burton K. H. Chia, Bertrand Denis, Jeff L. Froula, Wang, and et al. Critical Assessment of Metagenome Interpretation - a benchmark of computational metagenomics software. *bioRxiv*, page 099127, January 2017.

[9] Christopher S. Miller, Brett J. Baker, Brian C. Thomas, Steven W. Singer, and Jillian F. Banfield. Emirge: reconstruction of full-length ribosomal genes from microbial community short read sequencing data. *Genome Biology*, 12(5):R44, 2011.

[10] Cheng Yuan, Jikai Lei, James Cole, and Yanni Sun. Reconstructing 16s rRNA genes in metagenomic data. *Bioinformatics*, 31(12):i35–i43, June 2015.

[11] Eric P. Nawrocki, Diana L. Kolbe, and Sean R. Eddy. Infernal 1.0: inference of RNA alignments. *Bioinformatics*, 25(10):1335–1337, May 2009.
[12] Celine Mercier, Frederic Boyer, Aurelie Bonin, and Coissac Eric. Sumatra and sumaclust: fast and exact comparison and clustering of sequences, 2013.

[13] Evguenia Kopylova, Jose A. Navas-Molina, Céline Mercier, Zhenjiang Zech Xu, Frédéric Mahé, Yan He, Hong-Wei Zhou, Torbjørn Rognes, J. Gregory Caporaso, and Rob Knight. Open-source sequence clustering methods improve the state of the art. *mSystems*, 1(1), 2016.

[14] Evguenia Kopylova, Laurent Noé, and Hélène Touzet. Sortmerna: fast and accurate filtering of ribosomal rnas in metatranscriptomic data. *Bioinformatics*, 28(24):3211-3217, 2012.

[15] Evguenia Kopylova, Laurent Noé, Pierre Pericard, Mikael Salson, and Hélène Touzet. Sortmerna 2: ribosomal rna classification for taxonomic assignation. In *Workshop on Recent Computational Advances in Metagenomics*, ECCB 2014, 2014.

[16] Jared T. Simpson and Richard Durbin. Efficient de novo assembly of large genomes using compressed data structures. *Genome Research*, 22(3):549-556, 2012.

[17] Andreas Döring, David Weese, Tobias Rausch, and Knut Reinert. Seqan an efficient, generic c++ library for sequence analysis. *BMC Bioinformatics*, 9(1):1-9, 2008.

[18] Anton Bankevich, Sergey Nurk, Dmitry Antipov, Alexey A. Gurevich, Mikhail Dvorkin, Alexander S. Kulikov, Valery M. Lesin, Sergey I. Nikolenko, Son Pham, Andrey D. Prjibelski, Alexey V. Pyshkin, Alexander V. Sirotkin, Nikolay Vyahhi, Glenn Tesler, Max A. Alekseyev, and Pavel A. Pevzner. SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing. *Journal of Computational Biology*, 19(5):455-477, May 2012.

[19] Dinghua Li, Chi-Man Liu, Ruibang Luo, Kunihiko Sadakane, and Tak-Wah Lam. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. *Bioinformatics*, 31(10):1674-1676, May 2015.

[20] Konstantinos Mavromatis, Natalia Ivanova, Kerrie Barry, Harris Shapiro, Eugene Goltsman, Alice C. McHardy, Isidore Rigoutsos, Asaf Salamov, Frank Korzeniewski, Miriam Land, Alla Lapidus, Igor Grigoriev, Paul Richardson, Philip Hugenholtz, and Nikos C. Kyrpides. Use of simulated data sets to evaluate the fidelity of metagenomic processing methods. *Nature Methods*, 4(6):495-500, June 2007.

[21] Migun Shakya, Christopher Quince, James H. Campbell, Zamin K. Yang, Christopher W. Schadt, and Mircea Podar. Comparative metagenomic and rRNA microbial diversity characterization using Archaeal and Bacterial synthetic communities. *Environmental microbiology*, 15(6):1882-1899, June 2013.

[22] The Human Microbiome Project Consortium. Structure, function and diversity of the healthy human microbiome. *Nature*, 486(7402):207-214, June 2012.

[23] Robert C. Edgar, Brian J. Haas, Jose C. Clemente, Christopher Quince, and Rob Knight. UCHIME improves sensitivity and speed of chimera detection. *Bioinformatics*, 27(16):2194-2200, August 2011.

[24] Torbjørn Rognes, Tomáš Flouri, Ben Nichols, Christopher Quince, and Frédéric Mahé. VSEARCH: a versatile open source tool for metagenomics. *PeerJ*, 4:e2584, October 2016.

[25] Alla Mikheenko, Vladislav Saveliev, and Alexey Gurevich. MetaQUAST: evaluation of metagenome assemblies. *Bioinformatics*, 32(7):1088-1090, April 2016.

[26] Qiong Wang, George M. Garrity, James M. Tiedje, and James R. Cole. Naïve Bayesian Classifier for Rapid Assignment of rRNA Sequences into the New Bacterial Taxonomy. *Applied and Environmental Microbiology*, 73(16):5261-5267, August 2007.

[27] Miguel Pignatelli and Andrés Moya. Evaluating the Fidelity of De Novo Short Read Metagenomic Assembly Using Simulated Data. *PLOS ONE*, 6(5):e19984, May 2011.

[28] Weichun Huang, Leping Li, Jason R. Myers, and Gabor T. Marth. ART: a next-generation sequencing read simulator. *Bioinformatics*, 28(4):593-594, February 2012.

[29] Robert Schmieder and Robert Edwards. Quality control and preprocessing of metagenomic datasets. *Bioinformatics*, 27(6):863-864, March 2011.

[30] Marcel Martin. Cutadapt removes adapter sequences from high-throughput sequencing reads. *EMBnet.journal*, 17(1):pp. 10-12, May 2011.

[31] Daniel R. Mende, Ivica Letunic, Jaime Huerta-Cepas, Simone S. Li, Kristoffer Forslund, Shinichi Sunagawa, and Peer Bork. proGenomes: a resource for consistent functional and taxonomic annotations of prokaryotic genomes. *Nucleic Acids Research*, 45(D1):D529-D534, January 2017.

[32] Cyrielle Gasc, Eric Peyretaillade, and Pierre Peyret. Sequence capture by hybridization to explore modern and ancient genomic diversity in model and nonmodel organisms. *Nucleic Acids Research*, 44(10):4504-4518, June 2016.
# FEELnc: an alignment-free tool for long non-coding RNAs annotation

Valentin WUCHER, Fabrice LEGEAI<sup>2,3</sup>, Benoît HÉDAN<sup>1</sup>, Guillaume RIZK<sup>3</sup>, Lætitia LAGOUTTE<sup>1</sup>, Édouard CADIEU<sup>1</sup>, Audrey DAVID<sup>2</sup>, Nadine BOTHEREL<sup>1</sup>, Céline LE BÉGUEC<sup>1</sup>, Catherine ANDRÉ<sup>1</sup>, Christophe HITTE<sup>1</sup>, Thomas DERRIEN<sup>1</sup>

<sup>1</sup>IGDR, CNRS, UMR6290, University Rennes1, Rennes, Cedex 35043, France

<sup>2</sup>IGEPP, BIPAA, INRA, Campus Beaulieu, Le Rheu 35653, France

<sup>3</sup>INRIA, Genscale, Campus Beaulieu, Rennes 35042, France

Corresponding Author: tderrien@univ-rennes1.fr

**Paper Reference:** Wucher et al. (2017) FEELnc: a tool for long non-coding RNA annotation and its application to the dog transcriptome. *Nucleic Acids Research*. https://doi.org/10.1093/nar/gkw1306

**Abstract** Whole transcriptome sequencing (RNA-seq) has become a standard for cataloguing and monitoring RNA populations. One of the main bottlenecks consists in correctly identifying the different classes of RNAs among the plethora of reconstructed transcripts, particularly those that will be translated (mRNAs) from the class of long non-coding RNAs (lncRNAs). Here, we present FEELnc (FIExible Extraction of LncRNAs), an alignment-free program that accurately annotates lncRNAs based on a Random Forest model trained with general features such as multi k-mer frequencies and relaxed open reading frames. Benchmarking versus five state-of-the-art tools shows that FEELnc achieves similar or better classification performance on gold standard annotation datasets. Importantly, the program provides specific modules that enable the user to fine-tune classification accuracy, to formalize the annotation of lncRNAs and to identify lncRNAs even in the absence of an lncRNA training set. FEELnc moves beyond conventional coding potential classifiers by providing a standardized and complete solution for annotating lncRNAs. FEELnc is freely available at https://github.com/tderrien/FEELnc.

**Keywords** Long non-coding RNAs, Machine learning, Non-model species, Classification.

## 1 FEELnc description

FEELnc is an all-in-one solution from the filtering of non-lncRNA-like transcript models, to the computation of a new coding potential score and the automation of the definition of lncRNA classes. Based on a relaxed definition of ORFs and a very fast analysis of k-mer frequencies, the program implements an alignment-free strategy using Random Forests [1] to classify lncRNAs and mRNAs. Particularly, we developed FEELnc to be used on organisms for which no set of lncRNAs is available by deriving species-specific lncRNA models from mRNA sequences and automatically computing the coding potential score cut-off that maximizes classification performances.

<figure><img src="image_8.png" /><figcaption>**Fig 1.** FEELnc workflow and its 3 modules (from left to right): Filter, Coding\_potential, Classifier.</figcaption></figure>

The FEELnc workflow (Fig. 1) starts with a module, the filter, which aims at identifying non-lncRNA transcripts from the reconstructed transcript models given by transcriptome reconstruction tools [2]. The
second FEELnc module aims at computing a coding potential score given the assembled sequences and thus discriminates lncRNAs/mRNAs. Finally, the third FEELnc module, the classifier, employs a sliding window strategy around each lncRNAs to report all the reference biotypes/transcripts located within the window.

## 2 FEELnc benchmarks

In order to compare the performance of FEELnc with state-of-the art methods, we benchmarked tools based on human and mouse gold-standard GENCODE annotations [3]. In summary, FEELnc shows similar or better performance compared to five alignment-based and alignment-free tools (Table 1).

<table><caption>**Table 1.** Tools performance on the human data sets (similar results were observed in mouse). Bold values are the highest for each metrics. Programs are ranked by MCC (Matthews Correlation Coefficient).</caption><thead><tr><td>Programs</td><td>Alignment</td><td>Sensitivity</td><td>Specificity</td><td>Precision</td><td>Accuracy</td><td>MCC</td></tr></thead><tbody><tr><td>FEELnc</td><td>no</td><td><strong>0.923</strong></td><td>0.915</td><td>0.916</td><td><strong>0.919</strong></td><td><strong>0.838</strong></td></tr><tr><td>CPAT</td><td>no</td><td>0.899</td><td>0.924</td><td>0.922</td><td>0.912</td><td>0.823</td></tr><tr><td>CNCI</td><td>no</td><td>0.829</td><td>0.979</td><td>0.975</td><td>0.904</td><td>0.817</td></tr><tr><td>PLEK</td><td>no</td><td>0.732</td><td><strong>0.985</strong></td><td><strong>0.981</strong></td><td>0.858</td><td>0.741</td></tr><tr><td>PhyloCSF</td><td>yes</td><td>0.906</td><td>0.802</td><td>0.820</td><td>0.854</td><td>0.712</td></tr><tr><td>CPC</td><td>yes</td><td>0.699</td><td>0.739</td><td>0.728</td><td>0.719</td><td>0.438</td></tr></tbody></table>

## 3 Training FEELnc without known lncRNAs

One issue when using machine-learning approach is the requirement of both a positive and a negative set (here mRNA and lncRNA) to train the model. While the former is often available for most organisms, the latter is usually not, especially for many organisms [4]. To simulate lncRNAs, we also assessed two strategies called "shuffle" and "cross-species". The shuffle strategy is based on the paradigm that lncRNAs are derived from 'debris' of protein-coding genes [5]. To this end, we shuffled mRNA sequences from a reference annotation using the Ushuffle program [6], while preserving a given k-mer frequency of the input sequences. The cross-species strategy makes use of lncRNA sets annotated in other species to extract non-coding predictors and train the Random Forest model.

As expected, we showed that FEELnc performance for the cross-species strategy decreases with the time of speciation between the targeted species and the species providing lncRNAs ($Sp_{Pho} = -0.85$, $p_{val} < 10^{-4}$) and the shuffle strategy is a valuable alternative when no lncRNAs are annotated in closely related species.

## 4 Conclusion

We present FEELnc, a new program to annotate lncRNAs based on RNA-Seq assembled transcripts. In addition to providing good performance metrics, FEELnc allows to be self-trained by own users' datasets and also can be used for non-model organisms for which no set of lncRNAs is annotated.

## References

[1] Breiman L. Random Forests. Machine Learning. Kluwer Academic Publishers; 2001;45: 5-32. doi:10.1023/A:1010933404324

[2] Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology. 2010;28: 511-515. doi:10.1038/nbt.1621

[3] Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, et al. GENCODE: The reference human genome annotation for The ENCODE Project. Genome Res. 2012;22: 1760-1774. doi:10.1101/gr.135350.111

[4] Tagu D, Colbourne JK, Nègre NN. Genomic data integration for ecological and evolutionary traits in non-model organisms. BMC Genomics. 2014;15: 1-16. doi:10.1186/1471-2164-15-490

[5] Duret L, Chureau C, Samain S, Weissenbach J, Avner P. The Xist RNA gene evolved in eutherians by pseudogenization of a protein-coding gene. Science. 2006;312: 1653-1655. doi:10.1126/science.1126316

[6] Jiang M, Anderson J, Gillespie J, Mayne M. uShuffle: A useful tool for shuffling biological sequences while preserving the k-let counts. BMC Bioinformatics; 2008;9: 192-11. doi:10.1186/1471-2105-9-192
# Analyse intégrative des ARN longs non-codants (lncRNAs) du génome canin

Céline LE BÉGUEC¹, Valentin WUCHER¹,², Lætitia LAGOUTTE¹, Edouard CADIEU¹, Benoît HÉDAN¹, Catherine ANDRÉ¹, Christophe HITTE¹, Thomas DERRIEN¹

¹ Institut de Génétique et Développement de Rennes, CNRS UMR6290, Université Rennes1, 35000 Rennes, France. ² Adresse actuelle : Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Dr. Aiguader 88, Barcelona 08003, Spain.

Auteur pour correspondance : Thomas.derrien@univ-rennes1.fr

**Résumé** : Les ARN longs non codants (lncRNAs) constituent une famille d'ARN hétérogènes qui jouent un rôle majeur dans de nombreux processus biologiques. Nous avons développé récemment un outil d'annotation des lncRNAs appelé FEELnc, basé sur une approche sans alignement, et utilisant une stratégie 'Random Forest' entraînée sur les fréquences de multiples k-mer. Dans le cadre du consortium européen LUPA, FEELnc a permis d'identifier plusieurs milliers de nouveaux lncRNAs du génome du chien mais leur annotation fonctionnelle demeure dans son ensemble mal caractérisée. Dans ce travail, nous avons produit et analysé les profils d'expression des lncRNAs canins (n > 10000) sur 26 RNA-seq représentant une grande diversité de tissus. Nous présentons une caractérisation fonctionnelle des lncRNAs canins portant sur l'identification de leur signature d'expression spécifique par tissu, l'analyse de leurs roles potentiels comme régulateurs transcriptionnels, leurs niveaux de conservation et leurs contenus en éléments transposables. Nous avons identifié 4600 lncRNAs présentant un patron d'expression tissu-spécifique, avec près de 63 lncRNAs exprimés spécifiquement par tissu, suggérant un rôle essentiel dans la genèse et le maintien de l'intégrité des tissus. Nous avons construit un réseau de co-expression de l'ensemble des paires lncRNA:mRNA pour analyser des relations de cis-régulation. L'analyse statistique des paires a permis de déterminer des corrélations significatives (p-adjust BH <0,01) qui suggèrent un rôle régulateur pour plus de 900 lncRNAs. L'analyse des mRNAs cibles identifiés par les corrélations, basée sur les termes GO, a permis d'identifier des annotations fonctionnelles significativement enrichies (p <1e-4, FDR <0.05) correspondant aux processus de développement tels que le développement d'organes sensoriels et la croissance cellulaire. L'analyse de la co-occurrence des éléments transposables et des lncRNAs canins a montré que 84% des lncRNAs contiennent des SINEs spécifiques aux canidés. Ces co-occurences SINE:lncRNA suggèrent l'importance des éléments transposables dans l'expression tissu-spécifique des lncRNAs et dans la biogénèse des lncRNAs identifiés chez le chien.

**Mots clés** : lncRNAs, transcriptome, co-expression, éléments transposables, chien.

## 1. Introduction

Avec l'avancée des nouvelles technologies de séquençage haut débit, les analyses du transcriptome (RNA-seq) permettent d'identifier l'ensemble des classes d'ARN dont la classe des ARN longs non codants (lncRNAs) [1,2]. Le transcriptome correspond ainsi à l'ensemble des molécules d'ARN transcrites, avec ou sans capacité de coder des protéines, pour un temps et une condition ou tissu donnés [3]. Arbitrairement définis selon un critère de taille (généralement plus de 200 nucléotides), les lncRNAs possèdent des caractéristiques similaires à celles des ARN codant pour des protéines (mRNAs), c'est-à-dire qu'ils peuvent être épissés, et posséder (ou non) une queue de polyadénylation, mais ils se différencient par une absence de cadre de lecture ouvert fonctionnel. A la suite du séquençage de l'ensemble des transcrits ARN, l'annotation et la classification des différents ARN consiste à reconstruire les modèles de transcrits, à partir desquels il est crucial de définir les classes fonctionnelles des nouveaux transcrits identifiés. Les lncRNAs annotés peuvent
être localisés soit dans des régions intergéniques (lincRNAs), soit chevauchant et anti-sens des mRNAs (AS-lncRNAs) et représentent de bons candidats pour assurer un rôle de régulateur notamment des mRNAs proximaux.

Le modèle canin a émergé récemment comme une ressource nouvelle pour étudier la base génétique des traits complexes, incluant la morphologie, la physiologie et le comportement [4]. De multiples ressources génétiques et une séquence génomique de haute qualité positionnent désormais le chien comme une espèce modèle importante pour comprendre l'évolution du génome, la génétique de la population canine et les gènes sous-jacents aux traits phénotypiques complexes. Cependant, si les ressources génomiques nouvellement développées ont élargi notre compréhension du génome canin, l'annotation exhaustive des éléments fonctionnels tels que les ARN régulateurs de l'expression des gènes demeure nécessaire pour faciliter l'identification des relations génotype-phénotype [5]. Une des particularités du génome du chien réside dans la présence d'une famille d'éléments transposables spécifiques (SINEC_Cf) [6]. Ces rétrotransposons dégénérés (éléments auto-réplicants utilisant un ARN intermédiaire) sont dérivés d'un ARnt-Lys et ont été impliqués dans de nombreuses maladies génétiques [7,8] ou différences phénotypiques entre races de chiens [9,10]. Chez l'homme et la souris, les éléments transposables sont fréquemment retrouvés dans les lncRNAs [11,12] et de nombreuses études soulignent l'importance des éléments transposables dans la régulation de l'expression des gènes (voir pour revue récente [13]).

Nous avons récemment développé un outil appelé FEELnc qui analyse l'ensemble des transcrits et permet d'extraire, d'annoter et de classifier les ARN longs non codants [14]. Nous avons utilisé FEELnc dans le cadre du programme européen LUPA (http://eurolupa.org/), un consortium dédié au modèle canin, pour étendre l'annotation du répertoire des lncRNAs à partir de 26 RNA-seq de tissus distincts canins. Nous présentons dans cette étude une vaste analyse d'expression des >10000 lncRNAs chez le chien, leur corrélation d'expression avec les ARN messagers, et l'analyse de leur conservation et de leur contenu en éléments transposables. Nos résultats précisent la classification et les caractéristiques d'expression des lncRNAs au regard des mRNAs, la nature spécifique de leur expression et identifient un sous-ensemble de lncRNAs conservés entre l'homme et le chien.

## 2. Analyse du patron d'expression des lncRNAs

Nous avons produit et utilisé les données de l'annotation canine 'canFam3.1-plus' comportant 10444 gènes lncRNAs et 21810 gènes mRNAs [14]. Selon le protocole bioinformatique décrit dans Djebali et al. [15], nous avons aligné les lectures de 26 échantillons RNA-seq (>30 millions/RNA-seq), séquencés par HiSeq2000 à partir de banque cDNAs orientés, sur le génome de référence canin et à l'aide de l'annotation canFam3.1-plus avec le programme STAR (STAR_2.5.0a) [16] et avons déterminé l'expression des gènes et isoformes lncRNAs et mRNAs pour les 26 RNA-seq avec le programme RSEM (RSEM-1.2.25) [17]. L'analyse des 26 échantillons canins, représentant une grande diversité de tissus, a mis en évidence que 75% (7764/10444) des gènes lncRNAs étaient exprimés dans au moins un tissu (Transcripts Per Million i.e TPM >1). En comparant les distributions des niveaux d'expression entre gènes lncRNAs et les mRNAs, nous observons, comme chez l'homme [2] ou la souris [18] que les lncRNAs ont une plus faible expression (au moins 20 fois inférieure) pour tous les tissus analysés, à l'exception notable du tissu testiculaire ou la différence de moyennes d'expressions entre lncRNAs et mRNAs est moins importante (6 fois) (Fig. 1A).
<figure><img src="image_9.png" /><figcaption>**Fig 1.** (A) Analyse comparée des niveaux de transcription en log (TPM+1) entre gènes mRNAs (rouge foncé) et lncRNAs (bleu clair) dans 26 tissus. (B) Proportions de gènes lncRNAs (bleu clair) et mRNAs (rouge foncé) tissu-spécifique (le tissu testiculaire est représenté dans l'encadré étant donné la forte proportion de gènes exprimés spécifiquement dans ce tissu) (C) Clustering hiérarchique des 26 tissus basé sur les corrélations de Spearman mesurées à partir des données d'expression lncRNAs.</figcaption></figure>

Nous avons déterminé le niveau de tissu-spécificité par le calcul d'un score de spécificité tissulaire, "tau" [19], pour chaque gène lncRNAs et mRNAs (dont le niveau d'expression TPM >1) et avons ensuite mené une étude comparative de la spécificité tissulaire pour les 26 tissus. Ce score de spécificité, compris entre 0 et 1 (tau = 1 pour les gènes très spécifiques) est défini selon l'équation suivante :

$$ \tau = \frac{\sum_{i=1}^{n} (1 - \hat{x}_i)}{n-1}, \quad \hat{x}_i = \frac{x_i}{\max_{1 \le i \le n} (x_i)} $$

avec n correspondant au nombre de tissus analysés et $x_i$, l'expression du gène dans le tissu i. Dans une récente étude comparative, ce score a été classé parmi les plus robustes pour annoter les gènes tissu-spécifiques (TS) [20]. En imposant un filtre strict tau ≥ 0,95 (qui correspond en moyenne à un ratio $max(x_i) / max_{2nd}(x_i)$ supérieur à 4), un total de 4600 lncRNAs tissu-spécifiques et 3600 mRNAs tissu-spécifiques ont été identifiés. Nous observons en moyenne 63 lncRNAs exprimés spécifiquement par tissu, en excluant le tissu testiculaire. Ces sous-ensembles de profils d'expression, majoritairement lncRNAs, révèlent de véritables signatures associées à chaque tissu qui suggèrent des fonctions spécifiques dans les
processus physiologiques des tissus (Fig. 1B). Comme chez l'homme [1,2], le tissu testiculaire est particulièrement enrichi en gènes tissus-spécifiques (~3000 lncRNAs), ce qui met en évidence la singularité de ce tissu, probablement dû à l'état plus relâché de sa chromatine et à la présence de nombreux types cellulaires [21]. Par extension, l'utilisation des lncRNAs tissu-spécifique dans le cadre d'une analyse de la variabilité d'expression ou du différentiel d'expression pourra permettre d'investiguer leur implication fonctionnelle liée à l'état pathologique ou à l'environnement.

À partir des corrélations de Spearman mesurées sur les données d'expression lncRNAs, nous avons construit une matrice de distance et réalisé un clustering hiérarchique (méthode d'agrégation de ward) (Fig. 1C), qui met à nouveau en évidence la spécificité du tissu testiculaire par rapport aux 25 autres tissus. De plus, le clustering permet aussi de visualiser le regroupement cohérent des tissus ayant des origines histologiques communes et/ou des proximités fonctionnelles comme les tissus tégumentaires, musculaires ou nerveux (Fig. 1C). Par exemple, nous identifions le lncRNA RLOC 00033166 exprimé dans les 5 tissus nerveux (TPM >2) et non transcrit dans tous les autres tissus. Ce lncRNA est par ailleurs localisé en anti-sens du gène NRG3 (Neuregulin 3) qui est impliqué dans la différentiation cellulaire des neuroblastes et représente donc un candidat potentiel à sa régulation.

### 3. Construction et analyse du réseau de co-expression

La majorité des lncRNAs qui ont pu être caractérisés chez d'autres espèces possèdent une fonction de régulation de l'expression d'un gène en coordonnant des processus épigénétiques, transcriptionnels ou post-transcriptionnels [22]. Des analyses statistiques de corrélations positives ou négatives des valeurs d'expression entre mRNAs et lncRNAs peuvent permettre de mettre en valeur des lncRNAs régulant l'expression de gènes codant pour des protéines sachant que la plupart des interactions lncRNA:mRNA validées fonctionnellement concernent des gènes qui se localisent dans une certaine proximité génomique (généralement <1 Mb) [23]. Par conséquent, la caractérisation positionnelle génomique [24] des lncRNAs vis-à-vis des mRNAs est une étape initiale et essentielle pour analyser des relations de cis-régulation potentielles. En utilisant le module de classification du programme FEELnc pour identifier l'ensemble des paires lncRNA:mRNA, nous avons pu identifier 9615 interactions classées selon 2 types (i) AS-lncRNA:mRNA défini par les lncRNAs chevauchant des mRNAs transcrit en anti-sens (n = 4531) et (ii) lincRNAs:mRNA défini par les lncRNAs intergéniques localisés à moins d'1 Mb d'un mRNA (n = 5083). En utilisant les données d'expression issues des 26 RNA-seq, nous avons ensuite appliqué le concept de guilt-by-association (coupable par association), principe qui repose sur l'observation que les transcrits co-exprimés sont plus susceptibles d'être co-régulés, de partager des fonctions similaires ou de participer à des processus biologiques similaires [25]. Les analyses statistiques de corrélations des 9615 données de co-expression, ont permis de déterminer des corrélations significatives (p Adjustment BH <0,01) pour 492 paires de type lincRNA:mRNA et pour 411 paires AS-lncRNA:mRNA.

Pour attribuer des fonctions potentielles aux lncRNAs, nous avons mené une étude d'enrichissement des termes GO (Biological Process) des mRNAs avec lesquels ils sont co-exprimés. Un ensemble de 12 termes GO ont été retrouvés significativement enrichis (p <1e-4; FDR <0,05) et correspondent aux processus de développement tels que 'sensory organ development' (GO:0007423, 28 gènes) et 'cell growth' (GO:0016049, 25 gènes). Ces résultats permettent de proposer une assignation fonctionnelle aux lncRNAs, bien que les prédictions fondées sur ces enrichissements sont dépendants des annotations sous-jacentes.

### 4. Analyse de la conservation des lncRNAs canins avec l'homme

Pour identifier des lncRNAs orthologues entre l'homme et le chien qui vont suggérer des fonctions potentiellement conservées entre espèces, nous avons utilisé la base de données Compara d'Ensembl [26] qui recense les régions synténiques orthologues via un alignement multiple complet de génomes de plusieurs espèces. Cet alignement, qui ne prend pas en compte le sens de transcription, implique de restreindre l'analyse aux lincRNAs (n = 5651) et ne pas considérer les AS-lncRNAs. L'analyse a donc consisté à cartographier les positions des lincRNAs canins sur le génome humain via Compara, et a permis d'identifier
727 (12,86%) lincRNAs possédant un orthologue humain. Par comparaison, l'analyse pour les mRNAs a déterminé que 66,02% des mRNAs ont un orthologue humain selon la même méthodologie.

L'analyse des relations d'orthologie nous a permis d'étudier la distribution des scores de spécificité tissulaire (tau) en fonction de leur conservation au cours de l'évolution. Nous avons ainsi observé que les gènes sans relation d'orthologie sont significativement plus tissu-spécifique que les gènes conservés (test wilcoxon ; p <2,2<sup>-16</sup>). Ces résultats montrent que les gènes possédant une expression localisée dans un type cellulaire ou un tissu sont également des gènes moins conservés au cours de l'évolution [27]. Ces observations suggèrent que ces gènes, en grande majorité les lincRNAs évoluent rapidement en terme de séquence, et ainsi se différencient pour devenir des acteurs spécifiques du développement, de la signalisation et de la régulation des principaux types et fonctions cellulaires.

## 5. LincRNAs et éléments transposables canins

Chez l'homme, 80% des lincRNAs chevauchent au moins un élément transposable (TE) et près de 40% des séquences lincRNAs sont dérivées de TEs [11,12] suggérant un rôle essentiel des TEs dans la genèse et l'évolution des lincRNAs [28]. Une des particularités du génome du chien réside dans la présence d'une famille spécifique de TEs, les SINEC_Cf, qui ont été montrés comme liés à la diversité phénotypique observée entre les races de chiens [6,9,10] et à l'origine de maladies [7,8].

<figure><img src="image_10.png" /><figcaption>**Fig 2.** Analyse des éléments transposables (TEs). (A) Couverture du génome canin, des lincRNAs et mRNAs pour différentes classes et familles de TEs. (B) Distribution du score de spécificité tissulaire en fonction de la présence ou non de TEs dans les lincRNAs (bleu clair) et mRNAs (rouge foncé).</figcaption></figure>

Nous avons croisé l'annotation RepeatMasker [29] de 4 classes (et leurs familles principales) d'éléments transposables canins (DNA transposons, LTRs, LINEs et SINEs) et montrons que ~84% (8793/10444) des lincRNAs canins contiennent au moins un TE et que ~20% des séquences exoniques des lincRNAs sont composées exclusivement de TEs (Fig. 2A). Cette proportion en TEs est inférieure à celle observée dans la totalité du génome du chien (37,7%) mais 2,5 fois plus importante que celle des séquences mRNAs (7,9%) mettant en évidence la forte prévalence des TEs dans les séquences lincRNAs. Comparé au génome, la plupart des familles de TEs sont sous-représentées dans les lincRNAs à l'exception des rétrovirus ERVL-MaLR (Fig. 2A) qui sont aussi significativement enrichis dans les lincRNAs humains [12]. L'analyse des données d'expression sur les 26 tissus de cette étude montrent que les gènes contenant des TEs sont transcrits de manière plus tissu-spécifique (Fig. 2B) que ceux ne contenant aucun TEs. Ces observations qui montrent la corrélation entre l'insertion de TEs et une expression tissu-spécifique, suggèrent que les
modifications rapides au niveau de la séquence des gènes par insertion de TEs vont concourir à une spécialisation fonctionnelle notamment des lncRNAs.

## 6. Conclusions

Nous reportons dans cette étude l'analyse d'expression de plus de 10000 lncRNAs canins dans 26 tissus distincts. Les résultats montrent que les lncRNAs canins, comme pour les autres espèces, sont exprimés plus faiblement que les mRNAs mais de manière spécifique à un tissu. Ainsi, nous avons identifié ici les signatures transcriptionnelles spécifiques aux 26 tissus analysés. La forte spécificité tissulaire des lncRNAs peut constituer un facteur limitant pour annoter le répertoire des lncRNAs d'une espèce puisque leur annotation va ainsi dépendre de la disponibilité de nombreux types cellulaires, de l'analyse de nombreux tissus et de multiples conditions de temporalité, pour définir de manière exhaustive le catalogue complet des lncRNAs et de leur isoformes.

L'étude des réseaux de co-expression révèle plus de 900 corrélations d'expression avec les gènes codant pour des protéines proches et met en évidence de possibles processus de cis-régulation géniques [30]. Par une approche de génomique comparative, nous identifions plus de 700 lincRNAs avec un orthologue humain et l'analyse du contenu en éléments transposables montre que ~84% des lncRNAs canins contiennent au moins un TE. Nous montrons que ces caractéristiques de conservation et de contenu en TE corrèlent avec une expression tissu-spécifique. Ces observations suggèrent que les lncRNAs canins évoluent rapidement en terme de séquence (insertion de SINEC_Cf) et ainsi se différencient pour devenir des éléments fonctionnels spécifiques du développement et de la régulation de fonctions cellulaires spécialisées.

Cette étude indique que les lncRNAs appartiennent à de multiples classes fonctionnelles et suggère des fonctions potentielles aux lncRNAs, dans des processus fondamentaux tels que la spermatogenèse et le développement, mais aussi dans des mécanismes plus spécifiques tels que le développement d'organes sensoriels et la croissance cellulaire.

## Remerciements

Nous remercions l'ensemble des collègues de l'équipe Génétique du Chien de l'IGDR et Sarah Djebali de l'INRA GenPhyse de Toulouse. Nous remercions Cani-DNA (ANR-11-INBS-0003) pour les échantillons canins, le consortium LUPA (http://eurolupa.org) et le BROAD Institute pour la partie séquençage.

## Références

[1] M. N. Cabili, C. Trapnell, L. Goff, M. Koziol, B. Tazon-Vega, A. Regev, and J. L. Rinn, "Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses.," *Genes & Development*, vol. 25, no. 18, pp. 1915-1927, Sep. 2011.

[2] T. Derrien, R. Johnson, G. Bussotti, A. Tanzer, S. Djebali, H. Tilgner, G. Guernec, D. Martin, A. Merkel, D. G. Knowles, J. Lagarde, L. Veeravalli, X. Ruan, Y. Ruan, T. Lassmann, P. Carninci, J. B. Brown, L. Lipovich, J. M. Gonzalez, M. Thomas, C. A. Davis, R. Shiekhattar, T. R. Gingeras, T. J. Hubbard, C. Notredame, J. Harrow, and R. Guigo, "The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression.," *Genome Res*, vol. 22, no. 9, pp. 1775-1789, Sep. 2012.

[3] Z. Wang, M. Gerstein, and M. Snyder, "RNA-Seq: a revolutionary tool for transcriptomics.," *Nature Publishing Group*, vol. 10, no. 1, pp. 57-63, Jan. 2009.

[4] E. K. Karlsson and K. Lindblad-Toh, "Leader of the pack: gene mapping in dogs and other model organisms," *Nature Publishing Group*, vol. 9, no. 9, pp. 713-725, Sep. 2008.

[5] L. Andersson, A. L. Archibald, C. D. Bottema, R. Brauning, S. C. Burgess, D. W. Burt, E. Casas, H. H. Cheng, L. Clarke, C. Couldrey, B. P. Dalrymple, C. G. Elsik, S. Foissac, E. Giuffra, M. A. Groenen, B. J. Hayes, L. S. Huang, H. Khatib, J. W. Kijas, H. Kim, J. K. Lunney, F. M. McCarthy, J. C. McEwan, S. Moore, B. Nanduri, C. Notredame, Y. Palti, G. S. Plastow, J. M. Reecy, G. A. Rohrer, E. Sarropoulou, C. J. Schmidt, J. Silverstein, R. L. Tellam, M. Tixier-Boichard, G. Tosser-Klopp, C. K. Tuggle, J. Vilkki, S. N. White, S. Zhao, H. Zhou, FAANG Consortium, "Coordinated international action to accelerate genome-to-phenome
with FAANG, the Functional Annotation of Animal Genomes project.," *Genome Biol*, vol. 16, no. 1, p. 57, Mar. 2015.

[6] W. Wang, "Short interspersed elements (SINEs) are a major source of canine genomic diversity," *Genome Res*, vol. 15, no. 12, pp. 1798-1808, Dec. 2005.

[7] L. Lin, J. Faraco, R. Li, H. Kadotani, W. Rogers, X. Lin, X. Qiu, P. J. de Jong, S. Nishino, and E. Mignot, "The sleep disorder canine narcolepsy is caused by a mutation in the hypocretin (orexin) receptor 2 gene," vol. 98, no. 3, pp. 365-376, Aug. 1999.

[8] M. Pelé, L. Tiret, J.-L. Kessler, S. Blot, and J.-J. Panthier, "SINE exonic insertion in the PTPLA gene leads to multiple splicing defects and segregates with the autosomal recessive centronuclear myopathy in dogs.," *Human Molecular Genetics*, vol. 14, no. 11, pp. 1417-1427, Jun. 2005.

[9] B. Hédan, S. Corre, C. Hitte, S. Dréano, T. Vilboux, T. Derrien, B. Denis, F. Galibert, M.-D. Galibert, and C. André, "Coat colour in dogs: identification of the merle locus in the Australian shepherd breed.," *BMC Vet. Res.*, vol. 2, p. 9, 2006.

[10] H. G. Parker, B. M. Vonholdt, P. Quignon, E. H. Margulies, S. Shao, D. S. Mosher, T. C. Spady, A. Elkahloun, M. Cargill, P. G. Jones, C. L. Maslen, G. M. Acland, N. B. Sutter, K. Kuroki, C. D. Bustamante, R. K. Wayne, and E. A. Ostrander, "An expressed fgf4 retrogene is associated with breed-defining chondrodysplasia in domestic dogs.," *Science*, vol. 325, no. 5943, pp. 995-998, Aug. 2009.

[11] A. Kapusta, Z. Kronenberg, V. J. Lynch, X. Zhuo, L. Ramsay, G. Bourque, M. Yandell, and C. Feschotte, "Transposable Elements Are Major Contributors to the Origin, Diversification, and Regulation of Vertebrate Long Noncoding RNAs," *PLoS Genet*, vol. 9, no. 4, p.e1003470, Apr. 2013.

[12] D. Kelley and J. Rinn, "Transposable elements reveal a stem cell-specific class of long noncoding RNAs.," *Genome Biol*, vol. 13, no. 11, p. R107, Nov. 2012.

[13] E. B. Chuong, N. C. Elde, and C. Feschotte, "Regulatory activities of transposable elements: from conflicts to benefits.," *Nature Publishing Group*, vol. 18, no. 2, pp. 71-86, Feb. 2017.

[14] V. Wucher, F. Legeai, B. Hédan, G. Rizk, L. Lagoutte, T. Leeb, V. Jagannathan, E. Cadieu, A. David, H. Lohi, S. Cirera, M. Fredholm, N. Botherel, P. A. J. Leegwater, C. Le Beguec, H. Fieten, J. Johnson, J. Alföldi, C. André, K. Lindblad-Toh, C. Hitte, and T. Derrien, "FEELnc: a tool for long non-coding RNA annotation and its application to the dog transcriptome.," *Nucleic Acids Res*, p. gkw1306, Jan. 2017.

[15] S. Djebali, V. Wucher, S. Foissac, C. Hitte, E. Corre, and T. Derrien, "Bioinformatics Pipeline for Transcriptome Sequencing Analysis.," *Methods Mol. Biol.*, vol. 1468, pp. 201-219, 2017.

[16] A. Dobin, C. A. Davis, F. Schlesinger, J. Drenkow, C. Zaleski, S. Jha, P. Batut, M. Chaisson, and T. R. Gingeras, "STAR: ultrafast universal RNA-seq aligner.," *Bioinformatics*, vol. 29, no. 1, pp. 15-21, Jan. 2013.

[17] B. Li and C. N. Dewey, "RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome.," *BMC Bioinformatics*, vol. 12, no. 1, p. 323, Aug. 2011.

[18] D. D. Pervouchine, S. Djebali, A. Breschi, C. A. Davis, P. P. Barja, A. Dobin, A. Tanzer, J. Lagarde, C. Zaleski, L.-H. See, M. Fastuca, J. Drenkow, H. Wang, B. Pei, S. Balasubramanian, J. Monlong, A. Harmanci, M. Gerstein, M. A. Beer, C. Notredame, R. G. Oacute, and T. R. Gingeras, "Enhanced transcriptome maps from multiple mouse tissues reveal evolutionary constraint in gene expression," *Nat Commun*, vol. 6, pp. 1-11, Jan. 2015.

[19] I. Yanai, H. Benjamin, M. Shmoish, V. Chalifa-Caspi, M. Shklar, R. Ophir, A. Bar-Even, S. Horn-Saban, M. Safran, E. Domany, D. Lancet, and O. Shmueli, "Genome-wide midrange transcription profiles reveal expression level relationships in human tissue specification.," *Bioinformatics*, vol. 21, no. 5, pp. 650-659, Mar. 2005.

[20] N. Kryuchkova-Mostacci and M. Robinson-Rechavi, "A benchmark of gene expression tissue-specificity metrics.," *Briefings in Bioinformatics*, Feb. 2016.

[21] F. Chalmel and A. D. Rolland, "Linking transcriptomics and proteomics in spermatogenesis.," *Reproduction*, vol. 150, no. 5, pp. R149-57, Nov. 2015.

[22] A. C. Mallory and A. Shkumatava, "LncRNAs in vertebrates: Advances and challenges.," *Biochimie*, Mar. 2015.

[23] A. R. Bassett, A. Akhtar, D. P. Barlow, A. P. Bird, N. Brockdorff, D. Duboule, A. Ephrussi, A. C. Ferguson-Smith, T. R. Gingeras, W. Haerty, D. R. Higgs, E. A. Miska, and C. P. Ponting, "Considerations when investigating lncRNA function in vivo.," *Elife*, vol. 3, p. e03058, Aug. 2014.

[24] I. Ulitsky, "Evolution to the rescue: using comparative genomics to understand long non-coding RNAs.," *Nature Publishing Group*, vol. 17, no. 10, pp. 601-614, Oct. 2016.

[25] M. Guttman, I. Amit, M. Garber, C. French, M. F. Lin, D. Feldser, M. Huarte, O. Zuk, B. W. Carey, J. P. Cassady, M. N. Cabili, R. Jaenisch, T. S. Mikkelsen, T. Jacks, N. Hacohen, B. E. Bernstein, M. Kellis, A. Regev, J. L. Rinn, and E. S. Lander, "Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals.," *Nature*, vol. 458, no. 7235, pp. 223-227, Mar. 2009.

[26] J. Herrero, M. Muffato, K. Beal, S. Fitzgerald, L. Gordon, M. Pignatelli, A. J. Vilela, S. M. J. Searle, R. Amode, S. Brent, W. Spooner, E. Kulesha, A. Yates, and P. Flicek, "Ensembl comparative genomics
resources.," *Database* (Oxford), vol. 2016, p. bav096, 2016.

[27] C. Kutter, S. Watt, K. Stefflova, M. D. Wilson, A. Goncalves, C. P. Ponting, and D. T. Odom, "Rapid Turnover of Long Noncoding RNAs and the Evolution of Gene Expression," *PLoS Genet*, vol. 8, no. 7, p. e1002841, Jul. 2012.

[28] C. P. Ponting, P. L. Oliver, and W. Reik, "Evolution and functions of long noncoding RNAs.," vol. 136, no. 4, pp. 629-641, Feb. 2009.

[29] A. Smit, R. Hubley, and P. Green, "RepeatMasker Open-3.0.," 1996.

[30] U. A. Ørom and R. Shiekhattar, "Long noncoding RNAs usher in a new era in the biology of enhancers.," vol. 154, no. 6, pp. 1190-1193, Sep. 2013.
# Analysis, Integration and Modeling of Cell Clustering Results in High-Dimensional Cytometry Data

Guillaume GAUTREAU, David PEJOSKI, Roger LE GRAND, Antonio COSMA, Anne-Sophie BEIGNON and Nicolas TCHITCHEK

CEA - Université Paris Sud 11 - INSERM U1184, Immunology of viral infections and autoimmune diseases, Fontenay-aux-Roses, France.

Corresponding author: nicolas.tchitchek@gmail.com

**Reference paper:** Gautreau et al. (2016) SPADEVizR: an R package for Visualization, Analysis and Integration of SPADE results. *Bioinformatics*. http://doi.org/10.1093/bioinformatics/btw708

**Abstract** Cytometry is an experimental technique measuring cell marker expressions at the single cell level. The recent increase in the number of markers simultaneously measurable has led to the development of new automatic gating algorithms. Especially, the SPADE algorithm has been proposed as a novel way to identify clusters of cells having similar phenotypes. While SPADE or other cell clustering algorithms are powerful approaches, complementary analysis features are needed to characterize better the identified cell clusters.

We have developed SPADEVizR, an R package designed for the visualization, analysis, and integration of cell clustering results. The available statistical methods allow highlighting cell clusters with relevant biological behaviors or integrating them with additional biological variables. Moreover, several visualization methods are available to better characterize the cell clusters. SPADEVizR can also generate mathematical models to predict biological variables, based on the cell cluster abundances.

These analysis features are essential to interpret properly the behaviors and phenotypes of the identified cell clusters.

**Keywords** Mass Cytometry, Automatic Gating, Statistics, Visualization, Predictive models

## 1 Introduction

Cytometry is an experimental technique used to characterize cell properties at the single cell level. Thanks to mass cytometry, the number of simultaneously measurable cell markers has increased up to 50 [1]. This increase of measurable cell markers has led to the development of new automatic gating algorithms to identify group of cells, also named cell clusters, having similar expressions for selected markers.

The SPADE algorithm [2] was developed to identify cell clusters in the context of mass cytometry data. SPADE is a hierarchical clustering-based algorithm combined to a density-based down-sampling procedure.

While SPADE is a powerful approach, the interpretation of the behaviors or phenotypes of the identified cell clusters can be challenging, in particular in the scope of a whole dataset. For instance, SPADE has no methods allowing to highlight cell clusters with a cell abundance statistically different between two biological conditions or associated with an additional biological variable. Moreover, SPADE lacks of visualization methods to deeply characterize the phenotypes of the cell clusters in the whole dataset.

We have developed SPADEVizR, an R package to visualize, analyze and integrate results provided by SPADE. This package extends the original SPADE outputs with techniques such as volcano plots, streamgraphs, parallel coordinates, heatmaps, or distograms. Moreover, several statistical methods allow the identification of clusters with important biological behaviors. SPADEVizR also has features allowing the quantification and the visualization of the quality of clustering results and can be used with results generated by algorithms different from SPADE.

We illustrated the capabilities of SPADEVizR in the context of a vaccine study to identify new B cell populations altered by MVA immunizations in macaques [3].

## 2 Statistical methods

SPADEVizR allows the identification of three types of relevant cell clusters obtained from automatic gating algorithms. Abundant Clusters (AC) correspond to clusters having a cell abundance statistically greater than
a specific threshold for a given set of samples, identified using one-sample t-tests. Differentially Abundant Clusters (DAC) correspond to clusters having a cell abundance statistically different between two biological conditions, identified using two-samples t-tests. Correlated Clusters (CC) correspond to clusters having a cell abundance statistically correlated with an additional biological variable, identified using the Pearson or Spearman coefficients of correlation. These clusters can be visualized using scatter plot or volcano plot (Fig. 1A) representations.

## 3 Visualization methods

Boxplot (Fig. 1B) and kinetic representations available in SPADEVizR allows efficient visualizations and comparisons of cluster abundances between different samples and conditions. Moreover, streamgraph representations can display simultaneously absolute and relative cell abundances for a set of clusters (Fig. 1C).

Phenotypical characterization of the cell clusters can be performed using categorical heatmaps or parallel coordinates (Fig. 1D). While heatmaps provide global overviews, parallel coordinates provide more details by highlighting the homogeneity of marker expressions between the samples. SPADEVizR can generate multidimensional scaling representations to visualize the similarities between samples or clusters.

<figure><img src="image_11.png" /><figcaption>**Fig. 1.** Selected visualization representations available in SPADEVizR. (A) Volcano plot showing Differentially Abundant Clusters (DAC). (B) Boxplot showing the cell abundances for a given cluster in each sample and each condition. (C) Streamgraph showing absolute and relative abundances for a set of clusters across all the samples. (D) Parallel coordinates showing the phenotype of a given cluster.</figcaption></figure>

## 4 Conclusion

SPADEVizR constitutes a powerful approach for interpreting clustering results from several automatic gating algorithms. The available methods are very valuable to analyze properly high-dimensional cytometry data.

## Acknowledgements

French government program: "Investissement d'avenir: Equipements d'Excellence" (EQUIPEX) - 2010 FlowCyTech, Grant number: ANR-10-EQPX-02-01. Grant sponsor: Infrastructures Nationales en Biologie et Santé (INBS) - 2011 Infectious Disease Models and Innovative Therapies (IDMIT); Grant number: ANR-11-INBS-0008.

## References

[1] Sean C Bendall, Erin F Simonds, Peng Qiu, El-ad D Amir, Peter O Krutzik, Rachel Finck, Robert V Bruggner, Rachel Melamed, Angelica Trejo, Olga I Ornatsky, Robert S Balderas, Sylvia K Plevritis, Karen Sachs, Dana Pe'er, Scott D Tanner, and Garry P Nolan. Single-cell mass cytometry of differential immune and drug responses across a human hematopoietic continuum. *Science* (New York, N.Y.), 332(6030):687-96, May 2011.

[2] Peng Qiu, Erin F Simonds, Sean C Bendall, Kenneth D Gibbs, Robert V Bruggner, Michael D Linderman, Karen Sachs, Garry P Nolan, and Sylvia K Plevritis. Extracting a cellular hierarchy from high-dimensional cytometry data with SPADE. *Nat. Biotechnol.*, 29(10):886-91, October 2011.

[3] David Pejoski, Nicolas Tchitchek, André Rodriguez Pozo, Jamila Elhmouzi-younes, Rahima Yousfi-bogniaho, Christine Rogez-kreuz, Pascal Clayette, Yves Lévy, Antonio Cosma, Roger Le Grand, and Anne Sophie Beignon. Identification of Vaccine-Altered Circulating B Cell Phenotypes Using Mass Cytometry and a Two-Step Clustering Analysis. *J. Immunol*, 196:4814-4831, 2016.
Conférence invitée

Tobias MARSCHALL

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

Algorithms for Computational Genomics, Max-Planck-Institut für Informatik, Saabrücken, Germany

A Guided Tour to Computational Haplotyping

Humans and many other species are diploid. Every individual inherits two versions of each autosomal chromosome, called haplotypes, one from its mother and one from its father. Moving from (sequences of) genotypes to haplotypes is known as phasing or haplotyping. The knowledge of haplotypes is critical for addressing a variety of important questions in fundamental and clinical research. In this talk, I will highlight both algorithmic and experimental aspects of reconstructing haplotypes, with a special emphasis on recent technological advancements and their impact on the computational problems to be solved. I will briefly touch on population-based and pedigree-based phasing method, but will mostly focus on direct experimental methods that allow to reconstruct haplotypes for single individuals. Haplotype reconstruction from sequencing reads is most commonly formalized as the Minimum Error Correction (MEC) problem. Recent advances on fixed-parameter tractable (FPT) algorithm allow us to (quickly) solve practically relevant instances of this NP-hard problem optimally. I will present experimental results from five different platforms (PacBio, Oxford Nanopore, Hi-C, StrandSeq, and 10X Genomics) and highlight how combinations of these technologies allow to accurately reconstruct dense chromosome-length human haplotypes at manageable costs.
GenomeOnRails: depicting microbial species diversity via a pangenome graph

Guillaume GAUTREAU<sup>1,2,3,4</sup>, Rémi PLANEL<sup>1,2,3,4</sup>, Amandine PERRIN<sup>5,6</sup>, Marie TOUCHON<sup>5,6</sup>, Eduardo ROCHA<sup>5,6</sup>, Christophe AMBROISE<sup>2,7,8</sup>, Catherine MATIAS<sup>9</sup>, Stéphane CRUVEILLER<sup>1,2,3,4</sup>, Claudine MEDIGUE<sup>1,2,3,4</sup> and David VALLENET<sup>1,2,3,4</sup>

<sup>1</sup> CEA, Genoscope, 91000 Evry, France

<sup>2</sup> Université d'Evry, 91000 Evry, France

<sup>3</sup> CNRS-UMR-8030, 91000 Evry, France

<sup>4</sup> Université Paris-Saclay, 91000 Evry, France

<sup>5</sup> Microbial Evolutionary Genomics Unit, Institut Pasteur, 75724 Paris, France

<sup>6</sup> CNRS-UMR-3525, 75015 Paris, France

<sup>7</sup> Laboratoire de Mathématiques et Modélisation d'Evry (LaMME), 91000 Evry, France

<sup>8</sup> UMR-CNRS-8071, 91000 Evry, France

<sup>9</sup> Laboratoire de Probabilités et Modèles Aléatoires (LPMA), 75252 Paris, France

Corresponding Author: ggautrea@genoscope.cns.fr

Comparative genomics approaches in microbiology now use thousands of genomes to analyze a given species in different environmental or medical contexts. By collecting and comparing these genomic sequences, many studies are focused on the overall gene content of a species (i.e. the pangenome) to understand its evolution in terms of core and accessory parts. The core genome is defined as the set of genes shared by all the organisms of a taxonomic unit (generally a species). Accessory part (variable regions) is crucial to understand the adaptive potential of bacteria and contains genomic regions that are exchanged between strains by horizontal gene transfer (i.e. the mobilome). As recently suggested [1], a consensus representation of multiple genomes would provide a better analysis framework than using individual reference genomes. Here, we introduce an extension of this concept, giving it a formal mathematical representation using a graph model built up from genes clustered into families.

Pangenomes are generally stored in a binary matrix denoting the presence or absence of each gene family across organisms. However, this structure does not handle the genomic organization of gene families in each organism. In our approach, we propose a graph model where nodes represent families and edges chromosomal neighborhood information. Indeed, it is known that core gene families share conserved organizations whereas variable regions are rather randomly distributed along genomes.

Based on this data structure, our method classifies gene families through an Expectation/Maximization algorithm based on Bernoulli mixture model. Moreover, in order to take in account the genomic context of gene families, we smooth the classification with neighborhood information using Markov random field model [2]. This approach splits pangenomes in three groups: (1) *persistent genome*, equivalent to a relaxed core genome (genes conserved in all but a few genomes); (2) *shell genome*, genes having intermediate frequencies corresponding to moderately conserved genes potentially associated to environmental adaptation capabilities; (3) *cloud genome*, genes found at very low frequency.

Pangenomics is a relevant paradigm for very large scale comparative genomics. Further development of this tool should provide solid bases for efficient study of the dynamics of pangenome species and almost exhaustive references for metagenomic studies.

References

[1] Chan A. P., Sutton G., DePew J., Krishnakumar R., Choi Y., Huang X-Z., Beck E., Harkins D. M., Kim M., Lesho E. P., Nikolich M. P. and Fouts D. E. A novel method of consensus pan-chromosome assembly and large-scale comparative analysis reveal the highly flexible pan-genome of Acinetobacter baumannii. *Genome Biology*, 2015.

[2] Ambroise C., Dang M. and Govaert G. Clustering of spatial data by the EM algorithm. *geoENV-I-Geostatistics for environmental applications*, pages 493-504, 1997.
# Assembly of heterozygous genomes with high order De Bruijn graph

Antoine LIMASSET<fnref n="1" />, Camille MARCHET<fnref n="1" />, Pierre PETERLONGO<fnref n="1" /> and Jean-François FLOT<fnref n="2" />

Corresponding author: antoine.limasset@gmail.com

## 1 Introduction

Heterozygous genome assembly constitutes a complex task for which no satisfying solution exists at the present time. Phasing chromosomes in diploid or polyploid species or in metagenomes is still an open problem, despite the emergence of new long-read technologies and other dedicated approaches. In particular, intra-species or inter-species variations are usually discarded and/or result in highly fragmented assemblies. Several approaches tried to produce phased contigs using regular or paired-end short reads. On such genomes, String graph approaches like MIRA[1] produce large contigs but are intractable on large genomes. De Bruijn graph approaches (Spades[2], platanus[3], discovarDenovo[4]) are limited to very high or very low heterozygosity rates since they produce unsatisfactory fragmented contigs in the intermediary rates.

## 2 BWISE, a high order De Bruijn graph assembler

In this work we propose a novel de Bruijn graph-based assembler called BWISE that performs the construction of a very high order De Bruijn graph. This assembly fully takes advantage of read lengths and paired-end relationships between reads. Thereby correct paths are allowed to be determined resolving haplotypes and genomic repeats.

Results on simulated datasets show that BWISE is able to produce order of magnitude longer contigs than state of the art methods on highly heterozygous genomes. We also show that BWISE is comparable to state of the art assemblers for "regular" haploid genomes, and has the potential to scale up to very large genomes as the human one.

Beyond paired-end reads, the proposed framework allows in principle the integration of long-range information (mate pairs, 3Cseq/Hi-C, 10X, Single Molecule Real-Time (SMRT) as well as Nanopore reads) to determine accurate long paths in the assembly graph, with the ultimate goal of generating a single high-quality contigs per chromosome.

## References

[1] Bastien Chevreux, Thomas Wetter, Sándor Suhai, et al. Genome sequence assembly using trace signals and additional sequence information. In *German conference on bioinformatics*, volume 99, pages 45–56. Heidelberg, 1999.

[2] Anton Bankevich, Sergey Nurk, Dmitry Antipov, Alexey A Gurevich, Mikhail Dvorkin, Alexander S Kulikov, Valery M Lesin, Sergey I Nikolenko, Son Pham, Andrey D Prjibelski, et al. Spades: a new genome assembly algorithm and its applications to single-cell sequencing. *Journal of computational biology*, 19(5):455–477, 2012.

[3] Rei Kajitani, Kouta Toshimoto, Hideki Noguchi, Atsushi Toyoda, Yoshitoshi Ogura, Miki Okuno, Mitsuru Yabana, Masayuki Harada, Eiji Nagayasu, Haruhiko Maruyama, et al. Efficient de novo assembly of highly heterozygous genomes from whole-genome shotgun short reads. *Genome research*, 24(8):1384–1395, 2014.

[4] Neil I Weisenfeld, Shuangye Yin, Ted Sharpe, Bayo Lau, Ryan Hegarty, Laurie Holmes, Brian Sogoloff, Diana Tabbaa, Louise Williams, Carsten Russ, et al. Comprehensive variation discovery in single human genomes. *Nature genetics*, 46(12):1350–1355, 2014.

<footnote n="1">¹ IRISA, 263 Avenue Général Leclerc, 35000 Rennes, France</footnote>

<footnote n="2">² ULB, Avenue Franklin Roosevelt 50, 1050 Bruxelles, Belgium</footnote>
# Rapid spectra comparison with data-mining algorithms: accessing post translational modification profiles on a sample scale

Matthieu DAVID¹², Guillaume FERTIN¹, Hélène Rogniaux² and Dominique TESSIER²

¹ LS2N UMR CNRS 6004, Université de Nantes, F-44300, Nantes, France

² INRA UR1268 Biopolymères Interactions Assemblages, F-44316 Nantes, France

Corresponding author: matthieu.david@univ-nantes.fr

**Reference paper:** David et al. (2016) *International Workshop on Algorithms in Bioinformatics*. https://doi.org/10.1007/978-3-319-43681-4_6

**Abstract** Peptide identification from mass spectrometry experiments remains a challenging task, often leading to numerous missing or erroneous assignments. Because of their prohibitive computational cost, almost all peptide identification algorithms reduce their search space and in consequence, limit their ability to identify peptides displaying post-translational modifications. To resolve this difficulty, we developed a novel algorithm able to handle the computational complexity of all-to-all comparisons of spectra in the context of large volumes. We illustrated the interest of our approach with results obtained from the analysis of a public dataset and showed its capability to resolve the drawbacks usually accompanying spectra identifications without any a priori filter (e.g. computation time and decreased sensibility).

**Keywords** Proteomics, Mass Spectrometry, Spectra Comparisons, Algorithms, Data Mining

## 1 Introduction

Tandem mass spectrometry is used to identify the proteins present in a mixture. To this end, peptides obtained from the digestion of the proteins by an enzyme are fragmented into a mass spectrometer. All the masses of the fragments are measured, and this series of masses generates an *experimental mass spectrum*.
The protein identification process relies, in a first step, on the interpretation of experimental spectra in terms of peptides. In a second step, identified peptides are used to infer the list of proteins present in the mixture. In order to interpret spectra as peptides, the most widely used methods compare experimental spectra with *theoretical spectra* inferred from in silico digestion of a protein database. A scoring function is used to evaluate the similarity between pairs of spectra, whose definition varies depending on the method; however, every such scoring function is based on the count of shared masses between spectra, as a high number of shared masses is largely accepted to lead to more reliable identifications. The number of shared masses between two given spectra is the scoring function that we used in this study.

Currently, sensibility and reliability of the protein identification process remain limited. In a classical tandem mass spectrometry experiment, up to 75% of the spectra are indeed still unassigned [1]. Achieving pairwise comparisons of tens of thousands experimental spectra against hundreds of thousands of theoretical spectra is impractical. In order to avoid an explosion of the computation time, each experimental spectrum is only compared with theoretical spectra exhibiting almost the same mass, thus reducing drastically the search space. However, since proteins display a high rate of post-translational modifications (PTMs), and because these PTMs alter the mass of the peptides, this search space reduction hinders the peptide identification process – the subset of theoretical spectra considered for identifying the peptide may not contain the correct spectrum. PTMs are therefore incriminated as the most frequent cause for missing assignments [2]. Even though traditional peptide identification engines may add some PTMs in their searches, considering too many PTMs leads once again to excessive computation time. Recent Open Modification Search (OMS) algorithms [3] resolve the loss of identifications due to modified peptide masses with a loosened mass filter. The processing time however remains a significant bottleneck, associated to a degradation of results quality with more false positives and false negatives [4].

This work proposed a novel OMS approach inspired from data-mining techniques, providing access to the PTMs profiles of the sample while avoiding the computation time explosion pitfall.
## 2 Methodology

We designed a new data-structure, *SpecTrees*, to identify spectra rapidly without a priori filtering. *SpecTrees* contains by construction the number of shared masses between any two spectra. Inherent property of the data layout, this information is however distributed inside the structure and requires retrieval through the specific procedure *SpecXtract*. Let $s_e$ (resp. $s_t$) be an experimental (resp. theoretical) spectrum. For each ($s_e, s_t$) displaying a sufficiently high number of shared masses, a supplementary module analyzes this pair of spectra in order to identify peptides, possibly containing PTMs. We used a traditional statistical validation to ensure the results correctness, by means of evaluation of the False Discovery Rate (FDR).
Since the reference paper cited above was published, we extended our tests to larger datasets. We downloaded a tandem mass spectrometry dataset from the PRIDE database (PXD001428) and obtained an experimental dataset containing 37,703 spectra. The protein database is originated from the Ensemble genome assembly and contains 500,685 theoretical spectra after replication of the digestion enzyme behavior.

## 3 Results

Altogether, *SpecTrees* and *SpecXtract* compute the number of shared masses between any pair of experimental and theoretical spectra ($s_e, s_t$). We achieved on the larger datasets a running time under 10 minutes for different sets of parameters, clearly outperforming traditional identification tools or other OMS approaches. In-memory occupation of the program did not exceed 3 GigaBytes, and compatibility with the current throughput of proteomic analyzes is therefore ensured.

Our experiment demonstrated the ability of our software to identify 11,404 spectra, approximately 28% of the total, with a false detection rate of about 1%. Among those identifications, some classical cases of missed-cleavage (absence of digestion from the enzyme at a given location between two peptides), carbamylation, deamidation, dioxidation, formylation, oxidation and isotopic peptides were found, therefore validating the consistency of our approach. Rare modifications were also highlighted with a sufficient confidence, some of which not yet described in PTMs databases. Our approach additionally confers more reliability to the FDR for modified peptides using search space stratification, where traditional algorithms fail to provide an homogeneous FDR [5].

## 4 Conclusion

The *SpecTrees* data-structure is designed to compare important collections of spectra. Without any preliminary mass filter in spectra comparisons, *SpecTrees* can highlight proteins with original PTMs useful to direct further biological experiments. Moreover, flexible design enables the use of our method to detect peak patterns within any set of masses. Finally, the execution time of the method is compatible with a routine use in mass spectrometry laboratories.

## Acknowledgements

This project is partly funded by the Région Pays de la Loire (France) GRIOTE program (2013-2018).

## References

[1] J. Griss, Y. Perez-Riverol, S. Lewis, D. L. Tabb, J. A. Dianes, N. Del-Toro, M. Rurik, M. W. Walzer, O. Kohlbacher, H. Hermjakob, R. Wang, and J. A. Vizcaino. Recognizing millions of consistently unidentified spectra across hundreds of shotgun proteomics datasets. *Nat. Methods*, 13(8):651–656, Aug 2016.

[2] B. Bogdanow, H. Zauber, and M. Selbach. Systematic Errors in Peptide and Protein Identification and Quantification by Modified Peptides. *Mol. Cell Proteomics*, 15(8):2791–2801, Aug 2016.

[3] E. Ahrne, M. Muller, and F. Lisacek. Unrestricted identification of modified proteins using MS/MS. *Proteomics*, 10(4):671–686, Feb 2010.

[4] S. Na, N. Bandeira, and E. Paek. Fast multi-blind modification search through tandem mass spectrometry. *Mol. Cell Proteomics*, 11(4):M111.010199, Apr 2012.

[5] G. Hart-Smith, D. Yagoub, A. P. Tay, R. Pickford, and M. R. Wilkins. Large Scale Mass Spectrometry-based Identifications of Enzyme-mediated Protein Methylation Are Subject to High False Discovery Rates. *Mol. Cell Proteomics*, 15(3):989–1006, Mar 2016.
Conférence invitée

Julio SAEZ-RODRIGUEZ

RWTH, Aachen University, Germany and EMBL-EBI

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

Network Models to Understand and Combat Cancer: from Clinical
Genomics to Biochemical Modelling

Large-scale genomic studies are providing unprecedented insights into the molecular basis of cancer, but it remains challenging to leverage this information for the development and application of therapies. We have performed an integrated analysis of the molecular profiles of 11,215 primary tumours and 1,001 cancer cell lines, along with the response of the cell lines to 265 anti-cancer compounds. This analysis finds alterations in tumours that can confer drug sensitivity or resistance, and sheds light on which data types are most informative to prioritize treatment. Integration of this data with various sources of prior knowledge, in particular signaling pathways and transcription factors, points at molecular processes involved in resistance mechanisms, and offer hypotheses for novel combination therapies. Our own analysis as well as the results of a crowdsourcing effort (DREAM challenge) reveals that prediction of drug efficacy is far from accurate, implying important limitations for personalised medicine. I will argue than an important missing aspect is the dynamics of signaling networks, and show how applying logic models, trained with phosphoproteomic measurements upon perturbations, can further improve our understanding of the molecular basis of drug resistance, thereby providing new treatment opportunities not noticeable by static molecular characterization.
# Rare and common epilepsies converge on a shared gene regulatory network providing opportunities for novel antiepileptic drug discovery

Andrée DELAHAYE-DURIEZ<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Prashant SRIVASTAVA<fnref n="1" />,*, Kirill SHKURA<fnref n="1" />,*, Sarah R. LANGLEY<fnref n="1" />,<fnref n="5" />,*, Liisi LAANISTE<fnref n="1" />, Aida MORENO-MORAL<fnref n="2" />,<fnref n="5" />, Bénédicte DANIS<fnref n="6" />, Manuela MAZZUFERI<fnref n="6" />, Patrik FOERCH<fnref n="6" />, Elena V. GAZINA<fnref n="7" />, Kay RICHARDS<fnref n="7" />, Steven PETROU<fnref n="7" />,<fnref n="8" />,<fnref n="9" />, Rafal M. KAMINSKI<fnref n="6" />, Enrico PETRETTO<fnref n="2" />,<fnref n="5" /> and Michael R. JOHNSON<fnref n="1" />

Corresponding Authors: andree.delahaye@inserm.fr
enrico.petretto@duke-nus.edu.sg
m.johnson@imperial.ac.uk

**Paper Reference:** Delahaye-Duriez et al. (2016) Rare and common epilepsies converge on a shared gene regulatory network providing opportunities for novel antiepileptic drug discovery. Genome Biology.
http://dx.doi.org/10.1186/s13059-016-1097-7

## Abstract

### Background

The relationship between monogenic and polygenic forms of epilepsy is poorly understood and the extent to which the genetic and acquired epilepsies share common pathways is unclear. Here, we use an integrated systems-level analysis of brain gene expression data to identify molecular networks disrupted in epilepsy.

### Results

We identified a co-expression network of 320 genes (M30), which is significantly enriched for non-synonymous *de novo* mutations ascertained from patients with monogenic epilepsy and for common variants associated with polygenic epilepsy. The genes in the M30 network are expressed widely in the human brain under tight developmental control and encode physically interacting proteins involved in synaptic processes. The most highly connected proteins within the M30 network were preferentially disrupted by deleterious *de novo* mutations for monogenic epilepsy, in line with the centrality-lethality hypothesis. Analysis of M30 expression revealed consistent downregulation in the epileptic brain in heterogeneous forms of epilepsy including human temporal lobe epilepsy, a mouse model of acquired temporal lobe epilepsy, and a mouse model of monogenic Dravet (SCN1A) disease. These results suggest functional disruption of M30 via gene mutation or altered expression as a convergent mechanism regulating susceptibility to epilepsy broadly. Using the large collection of drug-induced gene expression data from Connectivity Map, several drugs were predicted to preferentially restore the

<footnote n="1">¹ Division of Brain Sciences, Imperial College Faculty of Medicine, London, UK</footnote>

<footnote n="2">² MRC Clinical Sciences Centre, Imperial College London, London, UK</footnote>

<footnote n="3">³ Université Paris 13, Sorbonne Paris Cité, UFR de Santé, Médecine et Biologie Humaine, France</footnote>

<footnote n="4">⁴ PROTECT, INSERM, Université Paris Diderot, Sorbonne Paris Cité, Paris, France</footnote>

<footnote n="5">⁵ Duke-NUS Medical School, 8 College Road 169857 Singapore, Republic of Singapore</footnote>

<footnote n="6">⁶ Neuroscience TA, UCB Pharma, S.A. Allée de la Recherche, 60 1070 Brussels, Belgium</footnote>

<footnote n="7">⁷ The Florey Institute of Neuroscience and Mental Health, The University of Melbourne, Parkville, Victoria 3052, Australia</footnote>

<footnote n="8">⁸ The Centre for Neural Engineering, The Department of Electrical Engineering, The University of Melbourne, Parkville, Victoria 3052, Australia</footnote>

<footnote n="9">⁹ The Australian Research Council Centre of Excellence for Integrative Brain Function, Parkville, Victoria 3052, Australia</footnote>

<footnote>* These authors participated equally to this work as joint second authors</footnote>
downregulation of M30 in epilepsy toward health, most notably valproic acid, whose effect on M30 expression was replicated in neurons.

The code of R functions newly created for this study is provided as a convenient R-package at
https://github.com/adelahay/BrainCell (DOI:10.5281/zenodo.164147).

## Conclusions

Taken together, our results suggest targeting the expression of M30 as a potential new therapeutic strategy in epilepsy.

**Keywords** Epilepsy, Co-expression, Regulatory network, Protein-protein interactions, Valproic acid.
# Prediction of Disease-associated Genes by advanced Random Walk with Restart on Multiplex and Heterogeneous Biological Networks

Alberto VALDEOLIVAS¹,², Elisabeth REMY¹, Laurent TICHIT¹, Gaëlle ODELIN², Claire NAVARRO², Sophie PERRIN², Pierre CAU³, Nicolas LEVY³ and Anaïs BAUDOT¹

¹ Aix-Marseille Université, CNRS, Centrale Marseille, I2M UMR 7373, Marseille, France.

² ProGeLife, 8 Rue Sainte Barbe 13001, Marseille, France.

³ Aix-Marseille Université, INSERM, UMR_S910, Faculté de Médecine, France.

Corresponding author: alberto.valdeolivas@etu.univ-amu.fr, anais.baudot@univ-amu.fr

**Abstract** Rare monogenic diseases globally affect millions of persons, but many causative genes remain to be discovered. Several computational approaches have been developed to predict disease-associated genes. Guilt-by-association strategies on protein interaction networks, in particular, postulate that proteins lying in a close network vicinity are functionally-related and implicated in similar phenotypes.

However, current network approaches are limited as they do not exploit the richness of biological networks, which are both multiplex (i.e., containing different layers of physical and functional interactions between genes and proteins), and heterogeneous (i.e., containing both interactions between genes/proteins, and interactions between diseases). In the present study, we extended the Random Walk with Restart algorithm to leverage these complex biological networks.

We compared our algorithm to classical random walks thanks to a leave-one-out strategy. The Random Walk with Restart on multiplex and heterogeneous networks takes advantage of data pluralism and shows increased performances to predict known disease-associated genes. We finally applied it to predict candidate genes for the Wiedemann-Rautenstrauch Syndrome.

**Keywords** Random Walks with Restart, Biological Networks, Multiplex and Heterogeneous Networks, Disease-Gene Prioritization.

## 1 Introduction

Rare monogenic diseases are often opposed to common diseases, but they jointly affect millions of persons. Overall, the causative gene(s) are often unknown, many patients remain undiagnosed, and no treatment exists for most of them. The disease phenotypes are resulting not from perturbations of isolated genes or proteins, but of complex networks of molecular interactions [1,2,3]. Proteins, for instance, do not act in isolation, but rather interact with each other to perform their functions in signalling pathways or metabolic reactions. Thanks to the scaling of the experimental techniques allowing interaction discovery, recent years have witnessed the accumulation of interaction datasets. For instance, protein-protein interactions (PPI) are nowadays screened at the proteome scale revealing thousands of physical interactions between proteins. Interaction data are commonly represented as networks, in which the nodes correspond to the genes or proteins, and the edges to their interactions.

In this context, network-mining approaches are applied to study human diseases, and in particular rare monogenic diseases. The rationale underlying these network approaches for human diseases is the clustering of proteins participating to the same cellular functions or biological processes in close network vicinity. Consequently, mutated genes coding for network-related proteins will lead to the same or similar phenotypes [4]. Following this idea, the mapping of the human protein-protein interactome network, about 10 years ago, was used to reveal new disease-associated genes, or allowed prioritizing candidate disease genes [5,6,7]. It also unveiled unsuspected interactions between disease-causing proteins, such as between proteins coded by genes mutated in ataxias [8]. More globally, interaction networks can also help deciphering the etiology and physiopathology of diseases [3], and their comorbidity relationships through their distances in molecular networks [9].

Among the various network-based study of human genetic diseases, Random Walk with Restart (RWR) appears as one of the state-of-the-art guilt-by-association approaches to predict new candidate disease genes. It was initially applied to explore the surroundings of disease-associated protein *seed(s)* in a PPI network. Every
protein in the global network is ranked according to its affinity to the seed(s), thereby allowing the prioritization of new candidate disease proteins [10]. The RWR algorithm was then extended, in particular to leverage phenotypic information [11,12,13,14]. For instance, Li and Patra RWR algorithm [11] considers jointly a PPI network and a network of phenotypic similarities between diseases. The two networks are connected by bipartite protein-disease associations, and form a heterogeneous network, i.e. a network containing nodes and edges of different nature connected by bipartite associations.

However, a common feature and limitation of these approaches is that they ignore the rich variety of information on physical and functional relationships between genes and proteins. Indeed, not only PPI are nowadays described on a large-scale: affinity purification followed by mass-spectrometry experiments inform on the *in vivo* molecular complexes, pathways interaction data are cured and stored in dedicated databases. In addition, functional interactions can be derived, for instance by constructing a co-expression network from transcriptionics expression data. Overall, the exploitation of this diversity of interaction data is lagging behind.

Sets of networks sharing the same nodes, but in which edges belong to different categories or represent interactions of different nature, are known as multiplex (aka multi-layer or multi-slice) networks [15]. In a biological multiplex network, each layer contains a different category of physical and functional interactions between genes or proteins. The combination of the different interaction sources, each having its own features and bias, provides a complementary view on genes and protein cellular functioning.

We present here the extension of the RWR algorithm to multiplex and heterogeneous biological networks (RWR-MH). We demonstrate the increased performance of this algorithm when compared with classical and current RWR approaches. Finally, as a real-case biological example, we applied the RWR-MH algorithm to predict candidate genes for the Wiedemann-Rautenstrauch syndrome (WRS), whose responsible gene(s) remain unknown.

## 2 Methods

We constructed a multiplex network as described in [16], but updated from downloads on November 2016 (Tab 1). The multiplex network is composed of 3 layers of physical and functional interactions between genes and proteins: protein-protein interactions (PPI), Pathway interactions extracted from pathway databases and Co-Expression interactions derived from RNA-seq expression data. The network nodes correspond equally to genes or proteins.

Additionally, we built a disease-disease similarity network, in which the edges between 2 diseases correspond to significant phenotype similarities. Briefly, we retrieved diseases and their associated phenotypes from the Human Phenotype Ontology (HPO) [17]. We then computed the phenotype similarity between diseases, by measuring the relative information content of the common phenotypes of every disease pair. We thereby assumed that rare phenotypes in the HPO database are more informative than frequent ones, as proposed by [18]. Finally, the disease-disease similarity network is constructed by linking every disease to its 5 most similar diseases according to their phenotype similarity scores, as proposed in [11].

The multiplex and heterogeneous network is constructed by linking every layer of the multiplex with the disease-disease similarity network using known gene-disease bipartite associations extracted from OMIM [19].

<table><caption>Tab. 1. Size of the networks used in this work.</caption><thead><tr><td>Network</td><td>Number of nodes</td><td>Number of edges</td></tr></thead><tbody><tr><td>PPI</td><td>12 621</td><td>66 971</td></tr><tr><td>Pathways</td><td>10 534</td><td>254 766</td></tr><tr><td>Co-expression</td><td>10 458</td><td>1 337 347</td></tr><tr><td>Disease-disease similarity</td><td>6 947</td><td>28 246</td></tr></tbody></table>

We extended the RWR algorithm to consider multiplex and heterogeneous networks (RWR-MH). In a nutshell, starting from initial seed node(s), the RWR progresses following the graph topology, with a non-zero probability to jump back to the initial seed node(s) at each step. After a sufficiently large number of iterations, RWR reaches a stationary state. In this stationary state, each nodes is associated to a score reflecting its proximity or pertinence with respect to the initial seed(s). Our extended RWR-MH algorithm has the capacity to explore one network layer, but also to jump between the different layers of the multiplex, because the same
nodes are present in the different layers. In addition, it can jump to the disease-disease similarity network thanks to the gene-disease bipartite associations. At every time, the walk can return to the initial protein and/or disease seed node(s), with a defined probability that we set to 0.7, as in previous studies [10,11,12].

We applied a leave-one-out cross-validation (LOOCV) strategy to compare RWR-MH with other RWR approaches. For every disease in OMIM [19] associated to 2 genes or more, each disease-associated gene is removed one-by-one (we will later refer to this removed gene as the left-out gene). The remaining disease-associated genes and the disease itself are used as seed nodes, and the RWR algorithms are applied. The Cumulative Distribution Functions (CDFs) are used to evaluate and compare the performances of the different approaches. They display the percentage of left-out genes that are ranked within the top k genes.

# 3 Results

## 3.1 RWR-MH outperforms current RWR approaches

We first compared different algorithms: i) the classical RWR on a monoplex PPI network, ii) the RWR on a heterogeneous network (built with the PPI and the disease-disease similarity network, RWR-H, as proposed by [11]), iii) the RWR on a multiplex network (composed of 3 layers, namely PPI, pathways and co-expression, RWR-M), and finally iv) the RWR on the multiplex-heterogeneous network (RWR-MH).

We applied a LOOCV strategy to evaluate the performances of the different RWR algorithms (Methods) for the prediction of known disease-associated genes (Fig 1). The multiplex RWR-M shows a superior performance than the classical RWR on a monoplex network. It is however comparable to the heterogeneous RWR-H. The multiplex and heterogeneous RWR-MH is remarkably better, since more than 45% of the left-out genes are ranked within the top 20.

<figure><img src="image_14.png" /><figcaption>**Fig. 1.** Cumulative Distribution Functions of the rank of each left-out gene retrieved by the LOOCV strategy. The algorithms are the classical RWR applied to the PPI monoplex network, the RWR-M applied to the 3-layer multiplex network, the RWR-H applied to the heterogeneous network built from the PPI and the disease-disease similarity networks, and the RWR-MH applied to the multiplex-heterogeneous network.</figcaption></figure>

## 3.2 RWR-MH prediction of candidate genes for the Wiedemann-Rautenstrauch Syndrome

The Wiedemann-Rautenstrauch Syndrome (WRS; MIM code: 264090), also called neonatal progeroid syndrome, is a disorder characterized by intrauterine growth retardation with subsequent failure to thrive and short stature [20]. In addition, patients display a progeroid appearance at birth and during the infancy, decreased subcutaneous fat, hypotrichosis and macrocephaly [21]. Only a few case reports have been documented, and no gene(s) has been described as causative of the syndrome yet.
To illustrate our approach, we applied our RWR-MH algorithm using as seed the WRS disease node. For visualization purposes, we displayed only the top 25 ranked diseases and top 25 ranked genes scored by the RWR-MH algorithm (Fig 2).

<figure><img src="image_15.png" /><figcaption>**Fig. 2.** Multiplex and Heterogeneous Network linking the top 25 genes and top 25 diseases obtained when the RWR-MH algorithm is applied using WRS as seed (yellow node). Grey elliptical nodes represent diseases, while turquoise rectangles represent genes. Black edges account for the bipartite gene-disease associations; Grey edges are the similarity links in the disease-disease similarity network; Blue edges correspond to PPI interactions; Red edges represent co-expression relationships; Orange edges represent pathway interactions.</figcaption></figure>

We predicted the top scored genes as top candidates for being involved in WRS. Many of them, such as *FIG4*, *RNF113A* or *LMNA*, are implicated in diseases directly connected to WRS from phenotypic similarities. For instance, mutations in *LMNA* are responsible for Hutchinson-Gilford progeria syndrome (HGPS; MIM code: 176670) and other premature aging syndromes such as Mandibuloacral Dysplasia with type A Lipodysrophy (MAD-A; MIM code: 248370). However, no mutations were found targeting *LMNA* sequence by gene sequencing analyses of WRS patients [21,22]. Additionally, the RWR-MH algorithm evidenced *ZMPSTE24*, which is known to be responsible of severe progeroid syndromes such as restrictive dermopathy (RD, MIM code: 275210) [23]. This peptidase acts during the post-translation modifications of the prelamin A, coded by *LMNA*, to undergo the complete maturation to lamin A. The direct interaction between the products of *LMNA* and *ZMPSTE24* is missing in the databases we used to construct our multiplex network. The *ZMPSTE24* nodes is however retrieved through different trajectories in the random walk. Once again, no mutations were found in the *ZMPSTE24* gene among the 5 WRS patients [22].

An interesting result is the small subnetwork composed of the genes *IGF2*, *INS*, *INSR* and *RPS6KA3*, which all participate to the insulin pathway. We retrieved these genes as top candidates due to their associations with two different diseases linked to WRS. This pathway is suspected to play a role in WRS [22]. Similarly, a cluster of proteins related to the cell cycle and DNA repair is connected to WRS through the Wolf-Hirschhorn syndrome (MIM code: 194190). DNA repair defects are also suspected to be involved in WRS [22]. The next step will be to validate these predictions, for instance using exome-sequencing data. Overall, our extended guilt-by-association RWR-MH algorithm could be integrated in analysis pipelines to help predicting candidate genes for rare diseases.
## 4 Discussion

Random Walk with Restart (RWR) is one of the state-of-the-art guilt-by-association approaches to prioritize candidate disease genes. We here extended the classical RWR algorithm in order to navigate multiplex and heterogeneous networks. We also demonstrated the increased performance of the RWR multiplex and heterogeneous strategy by leave-one-out cross validations. This improvement is due to the ability of our algorithm to extract and integrate the information from many interaction sources. Other types of networks could be integrated in the future, for instance to include interactions with non-coding RNAs. In addition, it would be interesting to explore the impact on the results of different disease-disease network topologies, i.e. taking a different criteria to build the disease-disease phenotype similarity network.

RWR has been mainly employed in biology to predict disease-associated genes [10, 11, 12, 14, 13]. But it has also been applied to address other biological problems, such as the prediction of drug-target potential interactions [24], the identification of clusters from PPI networks [25], or the prediction of adverse drug reactions [26]. Multiplex approaches are likely to boost the results of all these different applications, and could also be adapted to study cellular functioning as a whole.

The next step will be to resolve the degree bias of the algorithm. Indeed, RWR algorithms and other network propagation methods are biased towards networks hubs [27]. Therefore, poorly-connected and not well-characterized proteins, which can also be related to diseases, are hard to detect. Biased random walks can tackle this issue by relating the probability of transition to the degree of the nodes [28].

## Acknowledgements

This work was supported by JOBIM 2017.

## References

[1] Hannah Carter, Matan Hofree, and Trey Ideker. Genotype to phenotype via network analysis. *Current opinion in genetics & development*, 23(6):611–21, dec 2013.

[2] Marc Vidal, Michael E Cusick, and Albert-László Barabási. Interactome networks and human disease. *Cell*, 144(6):986–98, mar 2011.

[3] Trey Ideker and Roded Sharan. Protein networks in disease. *Genome research*, 18(4):644–52, apr 2008.

[4] M Oti and H G Brunner. The modular nature of genetic diseases. *Clinical Genetics*, 71(1):1–11, jan 2007.

[5] Miguel Angel Pujana, Jing-Dong J Han, Lea M Starita, Kristen N Stevens, Muneesh Tewari, Jin Sook Ahn, Gad Rennert, Víctor Moreno, Tomas Kirchhoff, Bert Gold, Volker Assmann, Wael M Elshamy, Jean-François Rual, Douglas Levine, Laura S Rozek, Rebecca S Gelman, Kristin C Gunsalus, Roger a Greenberg, Bijan Sobhian, Nicolas Bertin, Kavitha Venkatesan, Nono Ayivi-Guedehoussou, Xavier Solé, Pilar Hernández, Conxi Lázaro, Katherine L Nathanson, Barbara L Weber, Michael E Cusick, David E Hill, Kenneth Offit, David M Livingston, Stephen B Gruber, Jeffrey D Parvin, and Marc Vidal. Network modeling links breast cancer susceptibility and centrosome dysfunction. *Nature genetics*, 39(11):1338–49, nov 2007.

[6] Insuk Lee, U Martin Blom, Peggy I Wang, Jung Eun Shim, and Edward M Marcotte. Prioritizing candidate disease genes by network-based boosting of genome-wide association data. *Genome research*, 21(7):1109–21, jul 2011.

[7] Kasper Lage, E Olof Karlberg, Zenia M Størling, Påll I Olason, Anders G Pedersen, Olga Rigina, Anders M Hinsby, Zeynep Tümer, Flemming Pociot, Niels Tommerup, Yves Moreau, and Søren Brunak. A human phenome-interactome network of protein complexes implicated in genetic disorders. *Nature Biotechnology*, 25(3):309–316, mar 2007.

[8] Janghoo Lim, Tong Hao, Chad Shaw, Akash J Patel, Gábor Szabó, Jean-François Rual, C Joseph Fisk, Ning Li, Alex Smolyar, David E Hill, Albert-László Barabási, Marc Vidal, and Huda Y Zoghbi. A protein-protein interaction network for human inherited ataxias and disorders of Purkinje cell degeneration. *Cell*, 125(4):801–14, may 2006.

[9] Jörg Menche, Amitabh Sharma, Maksim Kitsak, Susan Dina Ghiassian, Marc Vidal, Joseph Loscalzo, and Albert-László Barabási. Uncovering disease-disease relationships through the incomplete interactome. *Science*, 347(6224), feb 2015.

[10] Sebastian Köhler, Sebastian Bauer, Denise Horn, and Peter N Robinson. Walking the Interactome for Prioritization of Candidate Disease Genes. *AJHG*, 82(April):949–958, 2008.

[11] Yongjin Li and Jagdish C. Patra. Genome-wide inferring gene-phenotype relationship by walking on the heterogeneous network. *Bioinformatics*, 26(9):1219–1224, 2010.

[12] Yongjin Li and Jinyan Li. Disease gene identification by random walk on multigraphs merging heterogeneous genomic and phenotype data. *BMC genomics*, 13 Suppl 7(Suppl 7):S27, 2012.
[13] MaoQiang Xie, YingJie Xu, YaoGong Zhang, TaeHyun Hwang, and Rui Kuang. Network-based phenome-genome association prediction by bi-random walk. *PLoS ONE*, 10(5):1-18, 2015.

[14] Zhi Qin Zhao, Guo Sheng Han, Zu Guo Yu, and Jinyan Li. Laplacian normalization and random walk on heterogeneous networks for disease-gene prioritization. *Computational Biology and Chemistry*, 57:21-28, 2015.

[15] Federico Battiston, Vincenzo Nicosia, and Vito Latora. Structural measures for multiplex networks. *Physical Review E - Statistical, Nonlinear, and Soft Matter Physics*, 89(3):1-16, 2014.

[16] Gilles Didier, Christine Brun, and Anaïs Baudot. Identifying Communities from Multiplex Biological Networks. *PeerJ*, pages 1-9, 2015.

[17] Sebastian Köhler, Sandra C. Doelken, Christopher J. Mungall, Sebastian Bauer, Helen V. Firth, Isabelle Bailleul-Forestier, Graeme C M Black, Danielle L. Brown, Michael Brudno, Jennifer Campbell, David R. Fitzpatrick, Janan T. Eppig, Andrew P. Jackson, Kathleen Freson, Marta Girdea, Ingo Helbig, Jane A. Hurst, Johanna Jähn, Laird G. Jackson, Anne M. Kelly, David H. Ledbetter, Sahar Mansour, Christa L. Martin, Celia Moss, Andrew Mumford, Willem H. Ouwehand, Soo Mi Park, Erin Rooney Riggs, Richard H. Scott, Sanjay Sisodiya, Steven Van Vooren, Ronald J. Wapner, Andrew O M Wilkie, Caroline F. Wright, Anneke T. Vulto-Van Silfhout, Nicole De Leeuw, Bert B A De Vries, Nicole L. Washingthon, Cynthia L. Smith, Monte Westerfield, Paul Schofield, Barbara J. Ruef, Georgios V. Gkoutos, Melissa Haendel, Damian Smedley, Suzanna E. Lewis, and Peter N. Robinson. The Human Phenotype Ontology project: Linking molecular biology and disease through phenotype data. *Nucleic Acids Research*, 42(D1):966-974, 2014.

[18] Sarah K Westbury, Ernest Turro, Daniel Greene, Claire Lentaigne, Anne M Kelly, Tadbir K Bariana, Ilenia Simeoni, Xavier Pillois, Antony Attwood, Steve Austin, Sjoert Bg Jansen, Tamam Bakchoul, Abi Crisp-Hihn, Wendy N Erber, Rémi Favier, Nicola Foad, Michael Gattens, Jennifer D Jolley, Ri Liesner, Stuart Meacham, Carolyn M Millar, Alan T Nurden, Kathelijne Peerlinck, David J Perry, Pawan Poudel, Sol Schulman, Harald Schulze, Jonathan C Stephens, Bruce Furie, Chris Van Geet, Augusto Rendon, Keith Gomez, Michael A Laffan, Michele P Lambert, Paquita Nurden, Willem H Ouwehand, Sylvia Richardson, and Andrew D Mumford. Human phenotype ontology annotation and cluster analysis to unravel genetic defects in 707 cases with unexplained bleeding and platelet disorders. *Genome Medicine*, 7:36, 2015.

[19] Ada Hamosh, Alan F. Scott, Joanna S. Amberger, Carol A. Bocchini, and Victor A. McKusick. Online Mendelian Inheritance in Man (OMIM), a knowledgebase of human genes and genetic disorders. *Nucleic Acids Research*, 33(DATABASE ISS.):514-517, 2005.

[20] H. V. Toriello. Syndrome of the month: Wiedemann-Rautenstrauch syndrome. *J. Med. Genet.*, pages 256-257, 1990.

[21] Aslihan Kiraz, Samim Ozen, Filiz Tubas, Yusuf Usta, Ozgur Aldemir, and Yasemin Alanay. Wiedemann-Rautenstrauch syndrome: Report of a variant case. *American Journal of Medical Genetics, Part A*, 158 A(6):1434-1436, 2012.

[22] Jia Woei Hou. Natural Course of Neonatal Progeroid Syndrome. *Pediatrics and Neonatology*, 50(3):102-109, 2009.

[23] Claire L. Navarro, Pierre Cau, and Nicolas Lévy. Molecular bases of progeroid syndromes. *Human Molecular Genetics*, 15(SUPPL. 2):151-161, 2006.

[24] Hui Liu, Mengmeng Guo, Ting Xue, Jihong Guan, and Libo Luo. Screening lifespan-extending drugs in Caenorhabditis elegans via label propagation on drug-protein networks. *BMC Systems Biology*, 10(Suppl 4), 2016.

[25] Kathy Macropol, Tolga Can, and Ambuj Singh. RRW: repeated random walks on genome-scale protein networks for local cluster discovery. *BMC Bioinformatics*, 10(1):283, 2009.

[26] Xiaowen Chen, Hongbo Shi, Feng Yang, Lei Yang, Yingli Lv, Shuyuan Wang, Enyu Dai, Dianjun Sun, Wei Jiang, K. M. Giacomini, M. Roy, R. Dumaine, A. M. Brown, E. Lounkine, L. Yang, J. Chen, L. He, L. Yang, J. B. Pan, M. Kuhn, M. Rarey, B. Kramer, T. Lengauer, G. Klebe, M. Campillos, M. Kuhn, A. C. Gavin, L. J. Jensen, P. Bork, L. Brouwers, M. Iskar, G. Zeller, V. van Noort, P. Bork, F. Napolitano, E. Bresso, M. Kuhn, M. Campillos, I. Letunic, L. J. Jensen, P. Bork, Z. L. Ji, J. X. Zhang, Z. Gao, X. Chen, Z. L. Ji, Y. Z. Chen, D. Szklarczyk, Y. Li, J. C. Patra, X. Chen, M. X. Liu, G. Y. Yan, S. Kohler, S. Bauer, D. Horn, P. N. Robinson, Q. Jiang, M. Duran-Frigola, P. Aloy, W. Jiang, M. Zhou, Y. Lv, J. Pinero, C. J. Zheng, G. Bindea, J. Turkson, R. Jove, F. A. Zouein, M. N. Richard, J. F. Deniset, A. L. Kneesh, D. Blackwood, G. N. Pierce, H. Kobori, M. Nangaku, L. G. Navar, A. Nishiyama, A. E. Cain, R. A. Khalil, and R. A. Khalil. Large-scale identification of adverse drug reaction-related proteins through a random walk model. *Scientific Reports*, 6(August):36325, 2016.

[27] Sinan Erten, Gurkan Bebek, Rob M Ewing, and Mehmet Koyutürk. DADA: Degree-Aware Algorithms for Network-Based Disease Gene Prioritization. *BioData mining*, 4(1):19, 2011.

[28] Federico Battiston, Vincenzo Nicosia, and Vito Latora. Efficient exploration of multiplex networks. *New Journal of Physics*, 18(4):043035, 2016.
CRCmapper: models of core transcriptional regulatory circuitries

Violaine Saint-André¹,², Alexander J Federation³, Charles Y Lin³, Brian J Abraham¹, Jessica Reddy¹,⁴, Tong
Ihn Lee¹, James E Bradner²,⁵and Richard A Young¹,⁴

¹ Whitehead Institute for Biomedical Research, Cambridge, Massachusetts 02142, USA

² Institut Curie, CNRS UMR3244/UPMC, 26 rue d'Ulm, 75005 Paris, France

³ Department of Medical Oncology,
Dana-Farber Cancer Institute, Boston, Massachusetts 02215, USA

⁴ Department of Biology, Massachusetts Institute
of Technology, Cambridge, Massachusetts 02139, USA

⁵ Department of Medicine, Harvard Medical School, Boston,
Massachusetts 02115, USA

Corresponding Author: violaine.saint-andre@curie.fr

There is considerable evidence that the control of gene expression programs is dominated by a small number of transcription factors (TFs). In embryonic stem cells and a few other cell types, these core TFs collectively regulate their own gene expression, forming an interconnected auto-regulatory loop that is considered the core transcriptional regulatory circuitry (CRC). There is limited knowledge of core TFs for most cell types. We recently discovered that genes encoding known core TFs forming CRCs are driven by super-enhancers, which provides an opportunity to systematically predict CRCs in poorly studied cell types through super-enhancer mapping. We have developed a tool, CRCmapper [1], which enables users to generate transcriptional regulatory circuitry maps for any sample for which they can provide suitable ChIP-seq data. Applying CRCmapper to a large set of human samples we have generated core circuitry models for 75 human cell and tissue types. These core circuitry models recapitulate and expand on what on previous CRCs and should prove valuable for further investigating cell-type-specific transcriptional regulation in healthy and diseased cells.

References

[1] Saint-André V, Federation AJ, Lin CY, Abraham BJ, Reddy J, Lee TI, Bradner JE, Young RA. Models of human core transcriptional regulatory circuitries. *Genome Res.* 2016 Mar;26(3):385-96.
# Network Inference of Dynamic Models by the Combination of Spanning Arborescences

Anthony COUTANT and Céline ROVEIROL

Laboratoire d'Informatique de Paris Nord (LIPN), UMR CNRS 7030, Université Paris 13
99 avenue Jean-Baptiste Clément, 93430 Villetaneuse, France

Corresponding author: firstname.lastname@lipn.univ-paris13.fr

**Abstract** In this paper, we tackle the problem of generative learning of dynamic models from "fai" time series data (high #variables/#individuals ratio), leading to a high sensitivity of learned models to the dataset noise. To overcome this problem, we propose a method computing a mixture of many highly biased but optimal spanning arborescences obtained from many perturbed versions of the original dataset, introducing variance to counterbalance the strong arborescence bias. The method is theoretically at the boundary between structure oriented Bayesian model averaging and recent work on density estimation using mixtures of poly-trees through a perturb and combine framework, transposed to a dynamic setting. In practice, preliminary results on the recent DREAM D8C1 challenge are promising.

**Keywords** Network inference, Ensemble Learning, Model Averaging, Spanning Arborescences

## 1 Introduction

Accurate network inference, or generative structure learning, is a major problem in bioinformatics for the comprehension of systems of different kinds, such as regulatory networks [1] or cell signaling pathways [2]. Although a difficult problem in any dataset, generative structure learning on biological contexts is additionally challenged by the scarcity of the available datasets, obtaining such data often being very expensive and time consuming.

A major issue occuring in structure learning algorithms from data is the one of *data fragmentation*, where the computation of frequencies from data is not reliable enough to robustly estimate statistics and avoid over-fitting. This problem occurs more particularly when the dataset is scarce and suffers from a very high variables/samples ratio, such as in many biological contexts where the number of available sets of measurements does not exceed several tens in a good situation, for a living organism made of hundreds or thousands of genes and involving many more molecules.

Biological systems are also often characterized by their dynamic properties, and many biological phenomenon under study are actually time evolving processes, e.g. the diauxic shift of the *Saccharomyces cerevisiae* yeast [3], for which available datasets can take the form of time series. Together, the data scarcity and the objective of finding dynamic models reinforce the data fragmentation issue, since considering several time points in space for each variable in the model both increases the number of actually considered variables in the network to learn, and reduces the samples size of the original dataset, due to the underlying used sliding window.

A possible strategy to prevent overfitting in such situation is to reduce the learning algorithm variance by reducing the number of possible models. This can be achieved for example by constraining the search space so that the resulting subspace have good properties, or by constraining the search algorithm so that a subset of possible models is reachable. Used alone, this strategy can find a good local, even global, optimum, but relatively to a potentially inadequate space where a good solution for the overall learning problem is unexistent.

Another possibility is to find an asymptotic structure which is the result of a consensus between different models learned from the dataset [4,5]. This way, the lack of sufficient statistics on data is partly offset by an attempt to compute sufficient statistics on potentially very noisy intermediate models.

In this article, we use both solutions. The general behaviour of the algorithm we propose is to compute an expected *composite model* using Bayesian Model Averaging [4] theory and a set of expected features, the expected edge existence, computed from the learning of many *component models*. These components are themselves highly biased models, more precisely spanning arborescences [6], with the interesting property for