# Probing factor-dependent long-range contacts using regression with higher-order interaction terms

Raphaël MOURAD¹, Lang Li² and Olivier CUVIER¹

¹ Laboratoire de Biologie Moléculaire Eucaryote (LBME), CNRS, Université Paul Sabatier (UPS), 31000 Toulouse, France

² Center for Computational Biology and Bioinformatics (CCBB), Indiana University, 46202 Indianapolis, USA

Corresponding author: raphael.mourad@ibcg.biotoul.fr

**Reference paper:** Mourad et al. (2017). Uncovering direct and indirect molecular determinants of chromatin loops using a computational integrative approach. PLOS Computational Biology, in press.

**Abstract** Chromosomal organization in 3D plays a central role in regulating cell-type specific transcriptional and DNA replication timing programs. Yet it remains unclear to what extent the resulting long-range contacts depend on specific molecular drivers. Here we propose a model that comprehensively assesses the influence on contacts of DNA-binding proteins, cis-regulatory elements and DNA consensus motifs. Using real data, we validate a large number of predictions for long-range contacts involving known architectural proteins and DNA motifs. Our model outperforms existing approaches including enrichment test, random forests and correlation, and it uncovers numerous novel long-range contacts in Drosophila and human.

**Keywords** Epigenetics; Chromatin; Hi-C; ChIP-seq; Generalized linear model

## 1 Introduction

The comprehensive analysis of 3D chromatin drivers is currently a hot topic [1]. A growing body of evidence supports the role of insulator binding proteins (IBPs) such as CTCF, and cofactors like cohesin, as mediators of long-range chromatin contacts [2]. In this work, we used a generalized linear model with interactions (GLMI) to identify the molecular determinants of loops, including protein and DNA sequence. Using this model, we uncovered numerous novel DNA loops and underlying mechanisms in *Drosophila* and human.

## 2 Results

The model is formulated as follows:

$$ 
\begin{aligned}
\log (\mathrm{E}[\mathbf{y}|\mathbf{X}]) &= \beta_0 + \beta \mathbf{X} \\
&= \beta_0 + \beta_d \mathbf{d} + \beta_B \mathbf{B} + \beta_C \mathbf{C} + \beta_g \mathbf{g}
\end{aligned}
\quad (1) $$

Variable $\mathbf{y}$ denotes the number of high-throughput chromosome conformation capture (Hi-C) contacts for any pair of bins on the same chromosome. Variable set $\mathbf{X} = \{\mathbf{d}, \mathbf{B}, \mathbf{C}, \mathbf{g}\}$ comprises several variable subsets: the log-distance variable $\mathbf{d}$ (polymer effect), the bias variables $\mathbf{B}$ (GC content, fragment length and mappability), the confounding variable set $\mathbf{C}$ and the genomic variable of interest $\mathbf{g}$. Model (1) is very general and can be developed in multiple versions depending on the variable $\mathbf{g}$ of interest. In this highlight article, we will only present the simplest version of the model, although more sophisticated versions accounting for multiple proteins mediating loops have been developed. The model is available in the R package "HiCglmi" which can be downloaded from the Comprehensive R Archive Network.

Let consider a pair of bins that we call left bin (L) and right bin (R). The attribution for left and right bins is arbitrary. Let also consider a genomic feature $F_i$ (whose binding is colored in blue in Figure 1a), that represent binding sites of the same protein. For the genomic feature $F_i$, occupancy variables $\mathbf{z}_{iL}$ and $\mathbf{z}_{iR}$ denote the occupancies of $F_i$ on left and right bins, respectively. For an occupancy variable, a value of 0/1 means absence/presence of the corresponding feature on the bin, e.g. absence/presence of the protein on the bin (a value between 0 and 1 means partial overlap of the feature). A "homologous interaction" variable $\mathbf{g} = \mathbf{n}_{ii}$ is the product of $\mathbf{z}_{iL}$ and $\mathbf{z}_{iR}$ ($\mathbf{n}_{ii} = \mathbf{z}_{iL} \times \mathbf{z}_{iR}$). The associated $\beta_{n_{ii}}$ parameter reflects the extent by which the genomic feature $F_i$ interacts with itself through chromatin contacts (Figure 1a). For instance, distant CTCF
binding sites were shown to form loops in human [2].

We compared GLMI with existing methods for their ability to identify genomic features known to be in-
volved in long-range contacts (Figure 1b). Here we used the negative binomial regression as the best spec-
ification of the GLMI in the context of Hi-C data overdispersion. We compared GLMI with (1) enrichment
test (ET) on highly confident chromatin interaction pairs as previously [3], (2) correlation (Cor) on highly con-
fident chromatin interaction pairs [4] and (3) random forests (RF) discriminating highly confident chromatin
interaction pairs from non-interacting pairs [5]. We found that GLMI outperformed the other methods to de-
tect long-range contacts between known architectural protein binding motifs. Using GLMI, we also uncovered
novel long-range contacts between architectural proteins in *Drosophila* (Figure 1c) and in human (data not
shown).

<figure><img src="image_29.png" /><figcaption>a Homologous interaction variable $n_{ii}$</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>b Known protein motifs versus other
DNA-binding protein motifs</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>C</figcaption></figure>

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

**Fig. 1.** a) Illustration of homologous interaction variable. b) Comparison between GLMI, random forest (RF), enrichment test (ET) and correlation (Cor) to detect known long-range contacts between protein motifs. c) Heatmap of long-range contacts between architectural proteins in *Drosophila*.

**3 Conclusion**

Here, we propose to use a generalized linear regression with interactions (GLMI) to study the roles of genomic features such as DNA-binding proteins, motifs or promoters to bridge long-range contacts in the genome, depending on transcriptional status or motif orientation. GLMI has multiple assets over existing approaches such as enrichment test, correlation and random forests. Compared to enrichment test or correlation that respectively assesses the protein enrichment or correlation at highly confident loops, GLMI quantitatively links the frequency of all long-range contacts to complex co-occupancies of proteins while accounting for known Hi-C biases and polymer background. Moreover, GLMI accounts for localizations among protein binding, a strong issue when analyzing protein binding sites known to largely overlap over the genome. In contrast to random forests which are efficient predictive models but sometimes poor explanatory ones, GLMI allows to identify key chromatin loop driver proteins and motifs. GLMI can also uncover numerous mechanisms behind loop formation using higher-order interaction terms and proper confounding variables (see original article). For instance, GLMI can determine if a cofactor is necessary to mediate long-range contacts between distant protein binding sites.

**References**

[1] Caelin Cubenas-Potts and Victor G. Corces. Architectural proteins, transcription, and the three-dimensional organization of the genome. *FEBS Letters*, 589(20PartA):2923-2930, 2015.

[2] Suhas S. P. Rao, Miriam H. Huntley, Neva C. Durand, Elena K. Stamenova, Ivan D. Bochkov, James T. Robinson, Adrian L. Sanborn, Ido Machol, Arina D. Omer, Eric S. Lander, and Erez L. Aiden. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. *Cell*, 159(7):1665-1680, February 2015.

[3] Ka-Chun Wong, Yue Li, and Chengbin Peng. Identification of coupling dna motif pairs on long-range chromatin interactions in human k562 cells. *Bioinformatics*, 2015.

[4] Vera Pancaldi, Enrique Carrillo-de Santa-Pau, Biola Maria Javierre, David Juan, Peter Fraser, Mikhail Spivakov, Alfonso Valencia, and Daniel Rico. Integrating epigenomic data and 3D genomic structure with a new measure of chromatin assortativity. *Genome Biology*, 17(1):1-19, 2016.

[5] Bing He, Changya Chen, Li Teng, and Kai Tan. Global view of enhancer-promoter interactome in human cells. *Proceedings of the National Academy of Sciences of the United States of America*, 111(21):201320308-E2199, May 2014.
An integrative approach for predicting the RNA secondary structure for the
HIV-1 Gag UTR using probing data

Afaf SAAIDI<fnref n="1" /> supervised by Yann PONTY<fnref n="1" /> and Bruno SARGUEIL<fnref n="2" />

Corresponding author: yann.ponty@lix.polytechnique.fr

Structure modeling is key to understand the mechanisms of RNA retroviruses such as HIV. Many *in silico*
prediction approaches suggesting structural models of moderate to good accuracies are available. However,
the prediction methods could be further improved by taking advantage of both next generation sequencing
technologies and different experimental techniques such as enzymatic and SHAPE probing data [1]. In a
published article [2], we introduce and use a structural modeling method based on the integration of many
experimental probing data to direct predictions with the aim to find the most accurate structure lying in the
intersection of disjoint sources of experiments.

**Method.** High-throughput experimental data can be derived from two sources of experimental data: SHAPE [1] and enzymatic probing. We used the stochastic sampling [3] mode of RNASubopt to sample structural models from the Boltzmann distribution in a way that favors/forces compatibility with the derived constraints. Namely, SHAPE reactivity profiles were used as soft constraints, meaning that observed reactivity values were translated into pseudo-energies. Position-specific susceptibilities to RNAses cleavage were used as hard constraints by setting arbitrary cut-offs above which specific base are forced to be paired and unpaired. Both types of constraints reduce the space of possible conformations, leading to a set of structures that are maximally compatible with the provided data.

We posited that the optimal structure(s) should be energetically stable and supported by several experimen-
tal data. Thus, for each type of probe, we generate a set of structures compatible with experimentally-derived
constraints. We merge those sets, and performed a structural distance-based clustering across experimental
conditions, to generate several sets of structural models that are well-supported by experimental data. Clusters
were scored using three criteria, namely, their stability, coherence and diversity (recurrence across structural
conditions). Within the set of clusters returned by our clustering algorithm we elect clusters on the Pareto
frontier, ie clusters that are not strictly dominated with respect to these three criteria. Finally, representative
structures (centroids), corresponding to Maximum Expected Accuracy structures are built and returned.

**Results.** The HIV-1 sequence probed in this study corresponds to the 5'UTR preceding the gag coding region from the NL-4.3 strain (Genbank: AF324493.2). A sample set of 12 000 structures, covering 6 sources of experimental data, was generated. The clustering step led us to elect two optimal clusters, whose corresponding centroids were then assessed in the light of their compatibility with specific SHAPE data. This allowed us narrow down our proposed models to a single candidate, whose base pair conservation/covariation was confirmed by comparative analysis.

**Conclusion and perspectives.** Our integrative approach allowed us to implement a consensus structure com-
patible with many different experimental probing data. As further work, we project to use our approach with
other sources of experimental probing data and to exploit the alignment to build an additional constrained
sample instead of using it for validation.

**References**

[1] K. A Wilkinson, E. J Merino, and K. M Weeks. Selective 2'-hydroxyl acylation analyzed by primer extension (SHAPE): quantitative RNA structure analysis at single nucleotide resolution. *Nature protocols*, 1(3):1610–6, 2006.

[2] J. Deforges, S. de Breyne, M. Ameur, N. Ulryck, N. Chamond, A. Saaidi, Y. Ponty, T. Ohlmann, and B. Sargueil. Two ribosome recruitment sites direct multiple translation events within HIV1 Gag open reading frame. *Nucleic Acids Research*, 2017.

[3] Y. Ding and C Lawrence. A statistical sampling algorithm for RNA secondary structure prediction. *Nucleic Acids Research*, 31(24):7280–7301, 2003.

<footnote n="1">¹ AMIBio team, Laboratoire d'informatique de l'école polytechnique, Inria Saclay, 91120, Palaiseau, France</footnote>

<footnote n="2">² Laboratoire de cristallographie et RMN Biologiques, Faculté de pharmacie Paris Descartes, 75006, Paris, France</footnote>
# NCboost: a meta-classifier of pathogenic non-coding variants integrating multiple sequence conservation and unsupervised functional scores

Barthelemy CARON¹, Lotfi SLIM¹ and Antonio RAUSELL¹

¹ Clinical BioInformatics Lab - Imagine Institute, 24 Boulevard du Montparnasse, 75015, Paris, France

Corresponding Author: barthelemy.caron@institutimagine.org

Whole Genome Sequencing (WGS) allows the identification of rare variants in non-coding regions with potential regulatory consequences. In the study of rare genetic diseases, where the number of sequenced individuals is typically low, computational approaches are needed to prioritize variants for further functional follow-up. Over the last years, a number of methods have been developed to predict pathogenicity of non-coding variants (e.g. CADD, GWAVA, FunSeq2[1]). Such methods are often based on machine-learning algorithms trained on curated disease-causing variants, and exploit features such as DNA methylation, histone modifications and transcription factor binding sites. However, their predictive performance is still low and has been shown to mainly rely on sequence conservation scores.

Here, we present NCBoost, an XGBoost-based [2] meta-classifier, trained on non-coding pathogenic variants from the Human Gene Mutation Database (HGMD). NCboost exploits a comprehensive set of purifying selection signals at three levels: the affected position, the surrounding non-coding region and the closest gene. In order to propose an integrative strategy for non-coding variants prioritization, we additionally included the functional scores of two unsupervised methods for the detection of regulatory variants: DeepSEA[3] and Eigen[4], both of them integrating epigenetic functional and conservation features.

NCboost performance was evaluated both by 10-fold cross validation and against an independent set of pathogenic variants from ClinVar database. A detailed comparative benchmark against state-of-the art methods is presented and results discussed in terms of the quality of the annotation sets and the targeted regulatory regions.

## Acknowledgements

We thank YuFei Luo for her contribution in variant annotation.

## References

[1] YaoFu et al, FunSeq2: a framework for prioritizing noncoding regulatory variants in cancer, Genome Biology, 2014.

[2] Chen, T. and Guestrin, C. (2016) XGBoost: A Scalable Tree Boosting System. ACM Press, pp. 785–794.

[3] Zhou, J. and Troyanskaya, O.G., Predicting effects of noncoding variants with deep learning-based sequence model. Nat. Methods, 12, 931–934, 2015.

[4] Iuliana Ionita-Laza et al, A spectral approach integrating functional genomic annotations for coding and non-coding variants, NatGen. 2016
# Context-specific prioritization of non-coding variants implicated in human diseases

Lambert MOYON¹, Yves CLÉMENT¹, Camille BERTHELOT¹ and Hugues ROEST CROLLIUS¹

¹ DYOGEN TEAM, IBENS, CNRS UMR8197, 46 rue d'Ulm, 75005, Paris, France

Corresponding Author: hugues.roest.crollius@ens.fr

High-quality whole genome sequences (WGS) of patients with rare or common diseases are increasingly being used to search for causative genetic variants that may explain the disease phenotype. In a large fraction of cases, no coding mutation in a known disease gene can be found, raising the possibility of at least two alternative causes: (1) a deleterious coding mutation in a new disease gene (2) a non-coding variant, for example modifying the function of an enhancer or a promoter of a disease gene. In the first case, a candidate-gene approach can be employed to investigate additional possibilities. In the second case, while methods exist to automatically annotate variants as potentially overlapping a regulatory region, investigators remain at loss for a reliable guide to efficiently prioritize the many variants that generally fall in this category, especially with respect to their impact on the patient's disease phenotype.

We present here a method aiming at prioritizing non-coding variants in a disease-oriented manner. By integrating functional, biochemical, and evolutionary information, our goal is to identify variants that are likely to be functional, with regards to the considered phenotype. This phenotype-oriented prioritization relies on the use of resources that link regulatory regions with target genes. In the context of a disease with known curated genes, this enables us to yield predictions of potentially causal variants that can be further investigated.

To evaluate our method, we are working in collaboration with the BRIDGE project from the University of Cambridge. More than 9,000 high-quality whole genome sequenceus from patients affected by rare diseases are now available, with a wide range of phenotypes, organized in 15 sub-projects. We applied our classification pipeline to cohorts of patients grouped by disease-phenotypes, and present here the specific properties of the predicted causal variants.

## References

[1] Albert, F.W., and Kruglyak, L. (2015). The role of regulatory variation in complex traits and disease. Nat Rev Genet 16, 197–212.

[2] Carss, K.J., Arno, G., Erwood, M., Stephens, J., Sanchis-Juan, A., Hull, S., Megy, K., Grozeva, D., Dewhurst, E., Malka, S., et al. (2017). Comprehensive Rare Variant Analysis via Whole-Genome Sequencing to Determine the Molecular Pathology of Inherited Retinal Disease. The American Journal of Human Genetics 100, 75–90.
Posters
# Analysis workflow for smallRNA-seq data

Sébastien NIN¹, Stéphanie RIALLE¹, Emeric DUBOIS¹ AND Laurent JOURNOT¹

¹ Plateforme MGX - Montpellier GenomiX, c/o Institut de Génomique Fonctionnelle, 141, rue de la cardonille, 34094 MONTPELLIER Cedex 05, FRANCE

Corresponding author: emeric.dubois@mgx.cnrs.fr

## Abstract

Small RNAs are small sequences, generally non-coding and playing a role in gene silencing effects, thus allowing the inhibition of the regulation or the translation of their target messenger RNA.

After the sequencing, reads obtained from small RNA libraries require a specific bioinformatics treatment in order to identify and quantify small RNAs. Here we focus on the analysis of micro RNAs.

After a quality control step, the 3' adapter sequence is trimmed from the reads (Cutadapt). Then, miRDeep2[1], a software able to quantify, identify micro RNAs (miRNAs) and discover new miRNAs through a secondary structure analysis, is used. It first maps the reads to a reference genome allowing 5 mulf-fi-mapping positions, then quantifies the miRNAs. We compared miRDeep2 to other tools and chose it for its performance in term of quantification of known miRNAs. From the quantification results, a gene differential expression analysis (edgeR[2], DESeq[3], DESeq2[4]) is done. The analysis ends with a target mRNA prediction (miRGate) and functional enrichment of Gene Ontology terms with TopGO (Bioconductor).

The workflow is used on the Montpellier GenomiX facility to provide a new analysis service.

## Citations

[1] Friedländer MR et al., miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades, Nucleic acids research 2012, Vol. 40, No 1: 37-52

[2] Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. *Bioinformatics*. 2010;26(1):139-140. doi:10.1093/bioinformatics/btp616

[3] Anders et al., Differential expression analysis for sequence count data, Genome Biology 2010, 11:R106 DOI: 10.1186/gb-2010-11-10-r106

[4] Love et al., Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, Genome Biology 2014 15:550, DOI: 10.1186/s13059-014-0550-8
# Analyse du réseau moléculaire impliqué dans le remodelage ventriculaire gauche post-infarctus du myocarde

Marie CUVELLIEZ¹, Christophe BAUTERS², Thomas KELDER³, Marijana RADONJIC³,Philippe AMOUYEL¹and Florence PINET¹

¹Inserm U1167, Institut Pasteur de Lille, 1 rue du Pr Calmette, 59019 Lille, France

²Inserm U1167, CHRU Lille, 2 avenue Oscar Lambret, 59037 Lille, France

³EdgeLeap B.V., Hooghiemstraplein 15, 3514 AX Utrecht, Pays Pas

Corresponding Author: florence.pinet@pasteur-lille.fr

Les maladies cardiovasculaires sont l'une des principales causes de mortalité dans les pays occidentaux. Suite à un infarctus du myocarde (IDM), 30% des patients développent un remodelage ventriculaire gauche (RVG) pouvant aboutir à une insuffisance cardiaque (IC).

REVE-2 est une étude multicentrique de 246 patients hospitalisés pour un premier IDM antérieur. Un suivi échocardiographique et des prélèvements sanguins sont réalisés à 4 temps après l'IDM : à l'inclusion (baseline), 1 mois, 3 mois et 1 an. A chaque temps, 25 variables moléculaires ont été dosées dans le plasma des patients (19 protéines et 6 ARNs non codants). L'expression de ces variables a été comparée entre les patients avec un RVG élevé (>20%) et un RVG faible (<20%) calculé entre l'échographie à un an et à l'inclusion.

Un réseau basé sur les interactions moléculaires connues a été construit à partir de l'ensemble des variables moléculaires dosées dans le plasma des patients de REVE-2. Pour cela, la plateforme EdgeBox (EdgeLeap), composée de 12 bases de données publiques (ENCODE, EnsemblGenes, HMDB, Microcosm, miRBase, miRRecords, miRTarBase, Reactome, STRING, TargetScan, Tfe et WikiPathways) a été utilisée. Le réseau inclut les variables REVE-2, les voisins directs de ces molécules et les molécules faisant partie des chemins les plus courts lorsque 3 interactions relient 2 variables REVE-2. Le réseau appelé REVE-2 est constitué de 1310 molécules dont 1263 protéines, 24 microARNs, 22 métabolites et un ARN long non codant. L'analyse du réseau a identifié 40 clusters qui ont été annotés pour les processus biologiques dans Gene Ontology (GO). Le réseau est visualisé sous Cytoscape (version 3.2.1). Une analyse des modules actifs a été réalisée pour chaque temps (inclusion, 1 mois, 3 mois et 1 an) à partir des variables REVE-2 significativement modulées entre les 2 groupes de patients (Pinet et al., 2017). La majorité des changements d'expression des molécules du réseau est observée en baseline et à 3 mois post-IDM, correspondant respectivement à la phase post-IDM et au développement du RVG. Une analyse de la centralité de chaque molécule a permis de déterminer leur importance dans le réseau, une centralité élevée suggérant un rôle crucial de la molécule dans le processus physiopathologique. Cette analyse a permis d'identifier de nouvelles molécules, non encore quantifiées chez les patients REVE-2, potentiellement impliquées dans le développement du RVG post-IDM. Les facteurs de transcription EP300, CTCF et ESR1 sont exprimés dans le cœur et ont été identifiés comme régulant l'activité de SOD2, protéine intervenant dans la régulation du stress oxydant dans les cardiomyocytes. Quatre ARNs non codants, les miR-26b-5p, miR-17-5p, miR-335-5p et miR-375, ont une centralité importante. Les deux premiers sont exprimés dans le cœur, contrairement aux 2 derniers. Les taux plasmatiques du miR-26b-5p sont diminués chez les patients avec une IC aigue et le miR-17-5p est impliqué dans l'apoptose des cardiomyocytes et la fibrose cardiaque.

L'analyse du réseau REVE-2 permet d'étudier les mécanismes physiopathologiques associés au RVG post-IDM au cours du temps, ainsi que d'identifier de nouveaux biomarqueurs potentiels pour détecter le RVG associé à un risque élevé d'IC.
# Mitochondrial genome variability of 205 Arabian endurance horses

Alexandre HEURTEAU¹, Claire HOEDE¹, Anne RICARD²,⁴, Diane ESQUERRE³, Caroline MORGENTHALER⁴, Nuria MACH⁴, Céline ROBERT⁴,⁵, Eric BARREY⁴,

¹ PF Bioinfo Genotoul, MIAT, Université de Toulouse, INRA, Castanet-Tolosan, France

² IFCE, Recherche et innovation, 61310 Exmes

³ GeT PlaGe, Genotoul, INRA Auzeville, Castanet Tolosan, France

⁴ GABI, BIGE, INRA, AgroParisTech, Université Paris-Saclay, Jouy-en-Josas, France

⁵ Ecole Nationale Vétérinaire d'Alfort, Université Paris-Est, Maison-Alfort, France

Corresponding Author: eric.barrey@inra.fr

## 1 Introduction

Endurance horses can run up to 160 km per race. Mitochondrial DNA variations may affect the efficiency of electron transport chain and ROS (Reactive Oxygen Species ie: peroxide) production, thus contributing to endurance performance. We studied the mitochondrial genome variability among 205 endurance horses for which racing performance were known. The objectives of this study were to propose a new strategy to call genetic variants in mitochondrial genome, and determine whether some variants are associated to endurance performance.

## 2 Material & methods

*Sampling & sequencing:* We designed 5 overlapping amplicons to specifically amplify mitochondrial genome from total DNA extracted in peripheral blood. Then we sequenced all samples with Illumina Miseq that produces 250 bp paired reads.

*Bioinformatics:* We chose a Arabian reference genome (GenBank ID JN398380.1) that we pseudo-circularized in silico to improve terminals' alignments. The cleaning and mapping pipeline is standard. We filtered reads with mapping quality < 60. We used the GATK package [1] to detect variations. Most parameters followed the best practices and some were fine-tuned.

*Statistics:* Association between SNPs and performance in endurance was tested using mixed model with fixed SNP Effect and random additive genetic effect with relationship matrix.

## 3 Results & conclusion

For the 205 horses, we found 590 combined variable positions (one position varies every 27 bps in average). Only 5 contain indels. 72% have already been found in previous studies. 80% of the protein coding variants are silent and the transition/transversion ratio is 22.5. Furthermore, we could observe 1.5% of non-haploid genotype (potential heteroplasmy) but some of them are probably numts (nuclear sequences of mitochondrial origin [2]). Preliminary study of statistical analysis can not yet link any variants to performance in endurance racing. Further analysis is still required.

## References

[1] Aaron McKenna, Matthew Hanna, Eric Banks, Andrey Sivachenko, Kristian Cibulskis, Andrew Kernytsky, Kiran Garimella, David Altshuler, Stacey Gabriel, Mark Daly, Mark A DePristo. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research, 20:1297-303, 2010.

[2] Solomon G. Nergadze, Manuel Lupotto, P. Pellanda, Marco Santagostino, Valerio Vitelli and Elena Giulotto. Mitochondrial DNA insertions in the nuclear horse genome. Animal Genetics, 41: 176-185, 2010.
Epigenetic heterogeneity in multiple myeloma

Jennifer RONDINEAU¹, Victor GABORIT¹,², Catherine GUERIN-CHARBONNEL¹,³, Philippe

MOREAU¹,⁴, Stéphane MINVIELLE¹,⁴, and Florence MAGRANGEAS¹,⁴

¹ CRCINA, INSERM, CNRS, Université d'Angers, Université de Nantes, Nantes, France

² LS2N, CNRS, Université de Nantes, Nantes, France

³ Institut de Cancérologie de l'Ouest, Nantes, France

⁴ CHU de Nantes, Nantes, France

Corresponding Author: jennifer.rondineau@univ-nantes.fr

Multiple myeloma (MM) is a hematological malignancy characterized by proliferation of malignant plasma cells. Despite the evolution of treatments over the last few years, MM remains an incurable disease, the vast majority of patients relapse of their disease in the years following treatment. MM represents about 20% of deaths from hematologic malignancies and 2% of cancer deaths. Although many types of genetic lesions (translocation, deletions, amplifications, mutations) have been identified, they don't fully explain the molecular mechanisms of relapse. We hypothesized that epigenetic changes contribute to the resistance of MM. Modified epigenetic genes could drive a novel mechanism of drug resistance in MM through changes in epigenetic states, subclonality and diversity. To test this hypothesis, we sought to highlight changes in the composition and diversity of epigenetic alleles during the progression of MM.

We used ERRBS (Enhanced Reduced Representation Bisulfite Sequencing) technique to analyze the DNA methylome of 17 patients at diagnosis and their relapse of MM, and 3 normal plasma cells (NPC) as control. This technique covers about 6% of CpGs of the genome and provides a single base resolution. Bisulfite reads were aligned using Bismark alignment software [1] to the bisulfite-converted hg19 genome, with non-directional model. We calculated the epigenetic changes between two stages using methclone [2]. This tool detects locus of 4 adjacent CpGs (minimum depth of 60X), called epiallele. We were particularly interested in the epiallele shift at relapse compared to diagnosis. Epigenetically shifted loci (eloci) are defined by a significant entropy shift (ΔS < -70).

We showed that epigenetic allelic diversification occurs during the initiation of the disease and also during progression and is highly variable between patients. Moreover, localization of eloci indicate important perturbations of regions with potential regulatory effects on gene expression during MM progression. These results are similar to observations made in other tumors, such as acute myeloid leukemia [3].

Acknowledgments

This work is supported by the "Fondation Française pour la Recherche contre le Myélome et les Gammapathies monoclonales".

References

[1] KRUEGER, Felix et ANDREWS, Simon R. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. bioinformatics, 2011, vol. 27, no 11, p. 1571-1572.

[2] LI, Sheng, GARRETT-BAKELMAN, Francine, PERL, Alexander E., et al. Dynamic evolution of clonal epialleles revealed by methclone. Genome biology, 2014, vol. 15, no 9, p. 472

[3] LI, Sheng, GARRETT-BAKELMAN, Francine E., CHUNG, Stephen S., et al. Distinct evolution and dynamics of epigenetic and genetic heterogeneity in acute myeloid leukemia. Nature Medicine, 2016, vol. 22, no 7, p. 792-799.

Keywords

Multiple Myeloma, epigenetic, whole genome DNA methylation, tumor heterogeneity
# NF-κB Landscape in Multiple Myeloma by high-throughput sequencing analysis

Victor GABORIT¹,², Jennifer RONDINEAU¹, Wilfried GOURAUD¹,³, Jérémie BOURDON², Stéphane

MINVIELLE¹,⁴ and Florence MAGRANGEAS¹,⁴

¹ CRCINA, INSERM, CNRS, Université d'Angers, Université de Nantes, Nantes, France

² LS2N, CNRS, Université de Nantes, Nantes, France

³ Institut de Cancérologie de l'Ouest, Nantes, France

⁴ CHU de Nantes, Nantes, France

Corresponding author: victor.gaborit@univ-nantes.fr

Nuclear factor κB (NF-κB) subunits RelA, RelB, cRel, p50 and p52 are each critical for B cell function such as regulation of cell proliferation or immune response. The NF-κB pathway is constitutively activated for more than 20% of patients with multiple myeloma (MM), a hematological malignancy characterized by an accumulation of abnormal plasma cell in bone marrow. Two different cytoplasmic mechanisms have been described [1] during NF-κB activation. First, the canonical pathway involving RelA:p50 dimer is essential for fast immune response. Second, a non-canonical pathway is activated in organogenesis, development and cell survival. It involves p52:RelB dimer. NF-κB nuclear activity and its effect in MM was studied using ChIP-seq experiments performed for each NF-κB subunit in MM cell line (MM.1S) in which the two NF-κB pathway are constitutively activated due to mutations inactivating TRAF3, a cytoplasmic inhibitor of the two NF-κB pathways of interest. Here, we present a complementary approach to characterize NF-κB nuclear activity in MM by combining different NGS datasets.

First, a ChIP-seq analysis is performed. This approach combines a quality assessment with fastqc, an alignment on hg19 human genome with bowtie2 and finally a peak calling using stringent threshold (p-value < 1e-7 and FDR < 1%) with macs2. We finally obtained 19,199 peaks for NF-κB subunit. In order to realize functional annotation of NF-κB peaks and MM.1S cell line, we used ChromHMM [2] with five histone marks to define 10 different chromatin states. We defined as *Promoters* all states showing enrichment for H3K4 trimethylation and H3K27 acetylation. *Enhancers* were defined as enrichment for H3K4 mono-methylation and H3K27 acetylation. To characterize how both pathway interacts in MM.1S, all κB regions are clustered using the k-mean method of seqMiner. It allows to identify different pattern of co-binding and correlate them to ATAC-seq enrichment. We also performed motif discovery using MEME-suite and motif enrichment determination with HOMER and ROC curves. Finally, we compared our results with DNAse-seq data of GM12878 (B cell lymphoma cell line) and ChromHMM data to compute differential enrichment of nucleosome free regions in this cell line.

Our analysis revealed that majority of NF-κB ChIP-seq peaks are located in promoters and in free accessible regions. Surprisingly, no κB motif enrichment was found for RelA and p52 subunit, but only in promoters binding RelB and in enhancers binding RelB and/or p50 subunit indicating the major role of this dimer for NF-κB specific binding. We also demonstrated that an interaction between the two NF-κB pathways exists in MM as it was shown in other malignancies [3]. Some differences on the changing of regulatory regions between those two cell lines shows the lack of enhancer location and κB motif enrichment in MM.1S cell line and implies specific mechanisms for this malignancies.

## Acknowledgements

This work is supported by the "Fondation Française pour la Recherche contre le Myélome et les Gammapathies monoclonales". Authors thank Magali Devic, Elise Douillard, Emilie Maurenton, and Nathalie Roi for excellent technical expertise.

## References

[1] Louis M Staudt. Oncogenic activation of NF-κB. Cold Spring Harbor perspectives in biology, 2(6):a000109, 2010.

[2] Jason Ernst, Pouya Kheradpour, Tarjei S Mikkelsen, Noam Shoresh, Lucas D Ward, Charles B Epstein, Xiaolan Zhang, Li Wang, Robbyn Issner, Michael Coyne, et al. Mapping and analysis of chromatin state dynamics in nine human cell types. Nature, 473(7345):43–49, 2011.

[3] Bo Zhao, Luis A Barrera, Ina Ersing, Bradford Willox, Stefanie CS Schmidt, Hannah Greenfeld, Hufeng Zhou, Sarah B Mollo, Tommy T Shi, Kaoru Takasaki, et al. The nf-κb genomic landscape in lymphoblastoid b cells. Cell reports, 8(5):1595–1606, 2014.
Efficient MinION data management

Laurent JOURDREN<fnref n="1" />, Aurélien BIRER<fnref n="1" />, Lionel FERRATO-BERBERIAN<fnref n="1" />, Sophie LEMOINE<fnref n="1" />
and Stéphane LE CROM<fnref n="1" />

Corresponding Author: jourdren@biologie.ens.fr

Keywords: High throughput sequencing, Oxford Nanopore Technologies, MinION, Data management, Demultiplexing, Quality control

Summary

In the last 12 months, the throughput of the MinION [1], a third generation sequencer provided by Oxford Nanopore Technologies (ONT), has dramatically raised. Today, with the last R9.4 chemistry more than one million of long sequence reads can be sequenced and much more are expected by the end of the year with the frequent updates (every 2-3 months) of the flowcell chemistry.

ONT choose to create one raw file (FAST5 file) for each read produced by the MinION. Moreover the
basecalling [2] of raw data is an intensive computational step as a neuronal network algorithm is used.
Hence, MinION data management and quality control is a disk space and time consuming task. Dealing
manually with the various input and output file formats from sequencer software requires many hours/days to
generate FASTQ files and QC reports.

In this poster, we present the best practices to deal with MinION data from raw data to QC reports in an
efficient manner. We now plan to create a new tool in order to automatically handle data transfer, read
demultiplexing conversion and quality control once a sequencing run has been finished as we previously
have done with Aozan [3] for Illumina sequencing.

Acknowledgements

This work was supported by the France Génomique national infrastructure, funded as part of the "Investissements d'Avenir" program managed by the Agence Nationale de la Recherche (contract ANR-10-INBS-09).

References

[1] Miten Jain *et al*. The Oxford Nanopore MinION: delivery of nanopore sequencing to the genomics community. *Genome Biol.*, 17(1):239, 2016.

[2] https://github.com/nanoporetech/nanonet

[3] Sandrine Perrin et al. Aozan: an automated post-sequencing data-processing pipeline. *Bioinformatics*, 2017.

<footnote n="1">¹ École normale supérieure, PSL Research University, CNRS, Inserm, Institut de Biologie de l'École normale supérieure (IBENS), Plateforme Génomique, 75005 Paris, France</footnote>
## Mise en place d'un pipeline de contrôle de la qualité de runs MinIon.

Lionel Ferrato-Berberian¹, Aurélien Birer¹, Ammara Mohammad¹, Corinne Blugeon¹, Fanny Coulpier¹, Stéphane Le Crom¹, Laurent Jourdren¹, Sophie Lemoine¹

¹École normale supérieure, PSL Research University, CNRS, Inserm, Institut de Biologie de l'École normale supérieure (IBENS), Plateforme Génomique, 75005 Paris, France.

Corresponding author : ferrato@biologie.ens.fr

Keywords : MinIon, nanopore, pipeline, analyse, QC

La technologie oxford nanopore MinIon est une technologie de séquençage qui permet l'obtention de lectures plus longues que la technologie illumina.

La technologie de séquençage ne repose pas sur des techniques d'imagerie comme illumina mais consiste à mesurer la différence de potentiel détecté lors du passage du brin d'ADN dans les pores.

On obtient alors des signaux électriques qui vont être convertis en lecture. Les données sont produites sous un format fast5. Les outils de suivis de runs existants ne sont plus adaptés ni aux métriques à suivre ni au format de fichier généré d'où la nécessité de développer des outils de contrôle qualité spécifiques.

Une bibliographie des outils d'analyses de la qualité des runs existants (poRe[1], poretools[2], ioniser[3], minotour[4]) a été réalisée pour prendre les meilleurs aspects de chaque outil afin de les combiner dans un nouveau pipeline d'analyse du run.

Ce poster présente le pipeline et les résultats produits par celui-ci de façon à permettre l'évaluation optimale des résultats d'un run MinIon. Notre pipeline a été développé et appliqué sur des données de runs de RNASeq MinIon comparant les transcrits WT et KO Egr2[5].

## Bibliographie

1 : https://github.com/mw55309/poRe_docs

2 : https://poretools.readthedocs.io/en/latest/

3 : http://bioconductor.org/packages/release/bioc/html/IONiseR.html

4 : http://minotour.nottingham.ac.uk/

5 : Topilko, P. et al. Krox-20 controls myelination in the peripheral nervous system. Nature 371, 796-799 (1994).
# Impact et évolution de la correction d'erreur sur des lectures longues issues de séquençage MinION Oxford Nanopore dans un contexte transcriptomique

Lionel FERRATO-BERBERIAN¹, Aurélien BIRER¹, Stéphane LE CROM¹, Laurent JOURDREN¹ and Sophie LEMOINE¹

¹ École normale supérieure, PSL Research University, CNRS, Inserm, Institut de Biologie de l'École normale supérieure (IBENS), Plateforme Génomique, 75005 Paris, France.

Corresponding author: slemoine@biologie.ens.fr

Les technologies de séquençage de 3ème génération comme celle proposée par Oxford Nanopore (ONT), permettent d'obtenir des lectures longues de plusieurs milliers de kilobases. Ces lectures longues ouvrent des champs nouveaux dans un contexte transcriptomique comme l'accès direct aux transcrits alternatifs et aux modifications post-transcriptionnelles.

Contrairement aux séquences issues des technologies de 2ème génération, les séquences longues comportent beaucoup d'erreurs, principalement des insertions-délétions. Ces erreurs peuvent être corrigées en utilisant des méthodes hybrides, les séquences courtes de type Illumina viennent corriger les lectures longues, ou des méthodes non-hybrides, basées sur une couverture importante du transcriptome. Cette étape de correction alourdit l'analyse des données menant à l'obtention d'un transcriptome de référence.

La qualité des séquences en sortie de MinION s'améliorant à chaque version de chimie, de flowcell, de logiciel d'appel de base, nous nous proposons de tester l'impact de la correction d'erreur sur différente version de séquençage MinION de façon à envisager l'opportunité ou non de la correction des lectures en RNA-Seq.
# Toullig: New pipeline for nanopore data analysis

Aurélien BIRER<fnref n="1" />, Lionel FERRATO-BERBERIAN<fnref n="1" />, Stéphane LE CROM<fnref n="1" />, Sophie LEMOINE<fnref n="1" />
and Laurent JOURDREN<fnref n="1" />.

Corresponding Author: birer@biologie.ens.fr

keywords: High throughput sequencing, Oxford Nanopore Technologies, MinION, RNA-Seq, Mapping, Quality control

## Summary

Since the release of Oxford Nanopore Technologies (ONT) MinION sequencers in 2014 [1] the number of reads produced with this new sequencing technology is still increasing. All the protocols remain in very active development (ONT provides updates of its chemistry and bioinformatics tools every 2-3 months). Hence the bioinformatical tools must be up to date throughout the development of MinION technology.

The IBENS (Institut de Biologie de l'École normale supérieure) genomic facility is currently developing a new data analysis workflow for RNA-Seq experiments using ONT sequencing output: *Toullig*. This pipeline is based on Eoulsan [2] and its bundled RNA-Seq pipeline for the Illumina reads.

The final goals of Toullig is to perform differential expression analysis from ONT long reads and produce a reference transcriptome by combining data from both Illumina and ONT technologies.

In this poster, we present our work on the Toullig pipeline with a focus on the long read mapping [3,4,5,6] and mapping quality control [7,8] steps.

The new Eoulsan modules for Toullig and the toolbox for manipulating ONT data are available on GitHub [9].

## References

[1] Miten Jain *et al*. The Oxford Nanopore MinION: delivery of nanopore sequencing to the genomics community. *Genome Biol.*, 17(1): 239, 2016.

[2] Jourdren, L., Bernard, M., Dillies, M.-A. & Le Crom, S. Eoulsan: a cloud computing-based framework facilitating high throughput sequencing analyses. *Bioinformatics* 28, 1542–1543 (2012).

[3] https://github.com/lh3/bwa.

[4] Sović, I. *et al*. Fast and sensitive mapping of nanopore sequencing reads with GraphMap. *Nat. Commun.* 7, 11307 (2016).

[5] Wu, T. D. & Watanabe, C. K. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. *Bioinforma. Oxf. Engl.* 21, 1859–1875 (2005).

[6] Kielbasa SM, Wan R, Sato K, Horton P, Frith MC. Adaptive seeds tame genomic sequence comparison. *Genome Res.* 2011 21(3):487-93.

[7] Weirather, J. L. *et al*. Comprehensive comparison of Pacific Biosciences and Oxford Nanopore Technologies and their applications to transcriptome analysis. *F1000Research* 6, 100 (2017).

[8] https://github.com/s-andrews/BamQC.

[9] https://github.com/GenomicParisCentre/toullig.

<footnote n="1">¹ École normale supérieure, PSL Research University, CNRS, Inserm, Institut de Biologie de l'École normale supérieure (IBENS), Plateforme Génomique, 75005 Paris, France.</footnote>
Extreme phenotypes define epigenetic and metabolic myeloid signatures in
cardiovascular diseases.

D. Seyres ¹,²,³, J.L. 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.H. Ouwehand ¹,²,⁸,⁹,¹⁰, M. Frontini ¹,²,¹⁰

¹ Department of Haematology, University of Cambridge, Cambridge Biomedical Campus, Cambridge,
CB2 0PT, UK

2 National Health Service Blood and Transplant (NHSBT), Cambridge Biomedical Campus,
Cambridge, CB2 0PT, UK

3 NIHR BioResource-Rare Diseases, University of Cambridge, Cambridge Biomedical Campus,
Cambridge, CB2 0QQ, UK

4 Medical Research Council Biostatistics Unit, University of Cambridge, Forvie Site,
Cambridge Biomedical Campus, Cambridge, CB2 OSR, UK

5 European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome
Campus, Hinxton, Cambridge, CB10 1SD, UK

6 Department of Clinical Biochemistry, Box 232 Addenbrooke's Hospital Hills Road Cambridge
CB2 0QQ, UK

7 University of Cambridge Metabolic Research Laboratories, Wellcome Trust-MRC Institute of
Metabolic Science, Box 289, Addenbrooke's Hospital, Cambridge CB2 0QQ, UK

8 The National Institute for Health Research (NIHR) Blood and Transplant Unit in Donor Health
and Genomics at the University of Cambridge, University of Cambridge, Strangeways Research
Laboratory, Cambridge, CB1 8RN, UK

9 Department of Human Genetics, The Wellcome Trust Sanger Institute, Wellcome Trust Genome
Campus, Hinxton, Cambridge, CB10 1SA, UK

10 BHF Centre of Excellence, Division of Cardiovascular Medicine, Addenbrooke's Hospital,
Cambridge Biomedical Campus, Cambridge, CB2 0QQ, UK

Corresponding Authors: ds777@medschl.cam.ac.uk, mf471@medschl.cam.ac.uk

Almost 6 million people are affected by cardiovascular disease (CVD) and 26% of deaths are
attributed to CVD in the UK each year. Thus, early detection of CVD onset is critical to improve quality of
life and to decrease the economic burden on healthcare providers. The identification of predictive biomarkers
and interpretable disease signatures by combining data obtained with different high-throughput omics is the
ultimate step towards reconstruction and analysis of complex multi-dimensional diseases, enabling deeper
mechanistic and medical insight. To this end, we collected data on 169 blood donors and 22 patients
representing three different diseases with high CVD risk. For each individual, monocytes and neutrophils
were isolated and underwent whole genome sequencing, ChIP-sequencing for histone modifications
representing regulatory elements (H3K4me1, H3K27ac), RNA-sequencing and DNA methylation analysis.
Additionally, plasma metabolites and lipids were quantified in all individuals.

We are currently integrating the different data types, working to select relevant covariates using penalized
likelihood approaches, followed by clustering via variational Bayes mixture models and other clustering
methods [1,2,3] according to continuous responses represented by clinical values associated with a high CVD
risk factor. Genetic contribution to a high CVD risk factor is analyzed by performing quantitative trait loci
mapping using histone modifications, methylation, transcripts expression, lipidomic and metabolic data.

Our aim is to identify metabolic and/or epigenetic signatures that could be applied for the detection of CVD
from an early onset.

References

[1] Witten & Tibshirani. A framework for feature selection in clustering. *Journal of the American Statistical Association*, 105(490), 713–726. 2010

[2] Liverani et al. PReMiuM: An R Package for Profile Regression Mixture Models Using Dirichlet Processes. Journal of Statistical Software, 64(7)). 2015

[3] Shen et al. Integrative clustering of multiple genomic data types using a joint latent variable model with application to breast and lung cancer subtype analysis. *Bioinformatics* 25 (22): 2906-2912, 2009
# Analyse et co-développement d'outils de bioinformatique destinés au traitement de données NGS issues de métagénomique virale

Emilie DELPUECH<fnref n="1" />, Guillaume CROVILLE<fnref n="2" />, Jean-Luc GUERIN<fnref n="2" />, Christophe KLOPP<fnref n="1" /> and Sarah MAMAN<fnref n="1" />

Corresponding Author: emilie.delpuech@inra.fr

Les techniques traditionnelles de détection des agents infectieux (PCR, recherche d'antigènes ou d'anticorps via la technique ELISA) par mise en culture (sur cellules ou œufs embryonnés), sont à la fois très lourdes, coûteuses et souvent pas suffisamment sensibles. Des approches dites « sans a priori » sont désormais susceptibles de les supplanter.

Les nouvelles générations de séquençage (NGS) sont devenues essentielles pour l'étude de l'identité et de la variabilité des génomes microbiens. Les techniques de type MiSeq (Illumina) et MinIon (Oxford Nanopore) sont aujourd'hui complémentaires, de par la nature et le volume de séquences générées. Une première étape d'étude bibliographique a permis de découvrir plusieurs pipelines répondant à la problématique de détection de virus. Trois d'entre eux ont été étudiés plus spécifiquement : Kraken [1], Truffle [2] et VIP [3].

Le pipeline Kraken [1] utilise une méthode d'alignement exact de K-mers pour réaliser une affiliation taxonomique des séquences d'agents pathogènes présentes dans des données NGS. Les K-mers sont recherchés dans la base de données Kraken. La correspondance exacte des 31-mers rend le pipeline stringent. Est-ce un facteur rédhibitoire pour détecter spécifiquement les séquences virales de données NGS de métagénomique virale ?

La détection de virus peut être réalisée à l'aide de l'outil Truffle [2] qui utilise des e-sondes représentantes des virus recherchés. Par conséquent, cet outil nécessite la connaissance des virus potentiellement présents dans les jeux de données. Ces e-sondes sont générées et par la suite, recherchées par alignement dans les données NGS. Un score est calculé pour chaque e-sonde. Le point fort de ce pipeline est l'obtention de scores représentatifs de la présence du virus, via une analyse statistique.

Le pipeline VIP [3] utilise quant à lui, une méthode d'assemblage de novo. Après un contrôle qualité et la suppression des séquences hôtes éventuelles, les lectures sont ensuite assemblées en contigs puis assignées via un alignment nucléotidique contre des bases de données virales spécifiques (ViPR/IRD), ou contre les bases de données virales du NCBI. A partir de cette assignation un assemblage de novo permet de reconstruire le génome des virus détectés. En complément, une analyse phylogénétique est proposée. Ce pipeline développé pour la détection de virus humains dans le cadre clinique, peut être appliqué à la détection de virus d'intérêt vétérinaire.

Cette étude a permis une prise en main de ces pipelines pour rendre chacun d'eux compatible, avec des jeux de données provenant de séquençages de tissus animaux. Les pipelines ont été testés et comparés afin d'optimiser leurs paramétrages. D'autre part nous réalisons une étude de faisabilité de l'intégration du pipeline VIP à la plateforme web d'analyse de données, Galaxy [4]. La mise à disposition d'un pipeline d'identification et de découverte de virus, est importante pour le diagnostic microbiologique, la surveillance de virus d'intérêt en santé publique et plus largement, la découverte de nouveaux virus.

## References

[1] D.E. Wood, S.L. Salzberg, Kraken: ultrafast metagenomic sequence classification using exact alignments, Genome Biol. 15 (2014) R46. doi:10.1186/gb-2014-15-3-r46.

[2] M. Visser, J.T. Burger, H.J. Maree, Targeted virus detection in next-generation sequencing data using an automated e-probe based approach, Virology. 495 (2016) 122–128. doi:10.1016/j.virol.2016.05.008.

[3] Y. Li, H. Wang, K. Nie, C. Zhang, Y. Zhang, J. Wang, P. Niu, X. Ma, VIP: an integrated pipeline for metagenomics of virus identification and discovery, Sci. Rep. 6 (2016) 23774. doi:10.1038/srep23774.

[4] B. Giardine, C. Riemer, R.C. Hardison, R. Burhans, L. Elnitski, P. Shah, Y. Zhang, D. Blankenberg, I. Albert, J. Taylor, W. Miller, W.J. Kent, A. Nekrutenko, Galaxy: A platform for interactive large-scale genome analysis, Genome Res. 15 (2005) 1451–1455. doi:10.1101/gr.4086505.

<footnote n="1">¹ Plateforme Bioinformatique Toulouse Midi-Pyrénées, UBIA, INRA, Chemin de Borde Rouge, 31326 Castanet-Tolosan, France</footnote>

<footnote n="2">² IHAP, Université de Toulouse, INRA, ENVT, 23 Chemin des Capelles, 31300 Toulouse, France</footnote>
3-SMART: Bioinformatic analysis of intronic polyadenylation regulation

Mandy CADIX¹,², Iris TANAKA², Pierre GESTRAUD¹, Marine SÉJOURNÉ¹,², Stéphan VAGNER², Nicolas SERVANT¹ and Martin DUTERTRE²

¹ Bioinformatics and Computational Systems Biology of Cancer (U900) - Institut Curie - PSL Research University - MINES ParisTech - Inserm, 26 rue d'Ulm, 75248 Paris cedex 05, France

² CNRS UMR3348 - Institut Curie - PSL Research University, Centre Universitaire Bâtiments 110, 91405 Orsay, France

Corresponding Author: mandy.cadix@curie.fr

The 3' end maturation of precursor mRNAs includes 3' end cleavage and addition of the polyadenylated tail. In addition, alternative polyadenylation (APA) occurs in about half of the genes in mammals. There are two types of alternative polyadenylation. The first one is located upstream of the last exon of the gene and is called intronic polyadenylation (IPA). The second is located within the last exon of the gene and is called polyadenylation in tandem. APA is widely regulated in oncogenic transformation and cell response to genotoxic (DNA-damaging) agents. 3'-seq (high-throughput sequencing of 3'-ends of polyadenylated transcripts) is a particular type of RNA sequencing technique. This methods allows sequencing the 3' end of the transcript upstream of the polyadenylated tail, in order to evaluate APA regulation.

The aim of this work was to develop a bioinformatics pipeline for differential 3'-seq analysis, called 3-SMART (3'-seq Mapping Annotation and Regulation Tool). Our pipeline includes 6 steps : 1/ Quality control and trimming of reads ; 2/ Reads mapping to genome ; 3/ Identification of peaks ; 4/ Annotation of peaks (IPA or last exon of genes) ; 5/ Filtering out of artefactual peaks corresponding to internal priming of oligo(dT) primers during reverse-transcription ; and 6/ Differential analysis of each IPA between two conditions. We applied this approach to a set of experimental data obtained in our lab to analyze IPA regulation in lung cancer cell response to cisplatin, a genotoxic anticancer drug. We found that cisplatin up-regulates IPA in many genes. The 3-SMART pipeline was implemented in bash and required a linux system. The pipeline is available at GitHub (https://github.com/InstitutCurie/3-SMART).
# Developing and sharing reproducible bioinformatics pipelines: best practices

Yohann LELIEVRE<fnref n="1" />, Audrey BIHOUÉE<fnref n="2" />, Eric CHARPENTIER<fnref n="2" />, Alban GAIGNARD<fnref n="2" />,<fnref n="4" />, Simon SOUCHET<fnref n="3" /> and
Damien VINTACHE<fnref n="1" />

Corresponding author: Yohann.Lelievre@univ-nantes.fr

## 1 Introduction

Life-sciences are nowadays conducted in multi-disciplinary and multi-centric studies. In this context, the same software components must be deployed in multiple environments for reproducibility and scalability issues. In addition, data analysis pipelines are usually composed of multiple components, continuously evolving, which leads to maintenance and long-term support challenges. To promote FAIR<fnref n="1" /> principles, providing controlled software environments becomes mandatory. We propose a set of best practices taking advantage of proven or promising tools: Git, Conda, SnakeMake[1], Jenkins and Docker.

## 2 Motivations

Bio-informaticians and software developers need to build data analysis pipelines in controlled environments to ensure long-term re-execution and better reproducibility. From an end-user point of view, typically a biologist, data analysis pipelines should be automatically installable in a local or dedicated computing infrastructure, including any software or data dependency. Pipelines should be launched in three steps: i) environment setup/activation, ii) parameters tuning, and iii) pipeline execution.

## 3 Approach and Results

The BiRD pipeline registry results from applying these guidelines in the context of Exome sequencing and RNAseq (variant calling, differential gene expression, gene fusion detection, single-cell). These pipelines are described in a GitLab web portal. GitLab allows i) to document the pipeline and its usage, and ii) to host and version the associated source code. To ease installation and dependency management, we packaged and deployed the executable software components through the Conda package manager in a dedicated repository<fnref n="2" />. To assess their long-term re-execution, workflows and associated software environments are nightly assembled into minimal Docker images through a Jenkins continuous integration system.

## 4 Conclusion and perspectives

The best practices hereby proposed aim at promoting findable and accessible data analysis pipelines through web-based resources. This process allows to re-package and re-execute pipelines in the long run, and to adapt to continuously evolving environments. Our future works include two main directions: i) handling data resources as part of the pipeline distribution process (e.g. BioMaj), and ii) studying how to promote interoperability between multiple systems and infrastructures. To enhance trust for end-users and to encourage reuse, provenance metadata and controlled vocabularies (e.g. EDAM) offer interesting perspectives to associate produced/analyzed with large-scale bio-resource registries such as BioTools.

## Acknowledgements

This work was supported by the BiRD core facility, the SyMeTRIC project and the GRIOTE project.

## References

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

<footnote n="1">¹ LS2N, UMR CNRS 6004, IMT Atlantique, ECN, Université de Nantes, Nantes, France.</footnote>

<footnote n="2">² l'institut du thorax, INSERM, CNRS, Université de Nantes, Nantes, France.</footnote>

<footnote n="3">³ Angers Academic Hospital, CHU d'Angers, France</footnote>

<footnote n="4">⁴ Nantes Academic Hospital, CHU de Nantes, France</footnote>

<footnote n="1">1. Findable - Accessible - Interoperable - Reusable</footnote>

<footnote n="2">2. https://anaconda.org/BiRD</footnote>
# MICROSCOPE: an integrated platform for the Exploration and Curation of Microbial Genomes

David VALLENET¹, Alexandra CALTEAU¹, Stéphane CRUVEILLER¹, Mathieu GACHET¹, Guillaume GAUTREAU¹, Adrien JOSSO¹, Aurélie LAJUS¹, Jordan LANGLOIS¹, Jonathan MERCIER¹, Hugo PEREIRA¹, Rémi PLANEL¹, Johan ROLLIN¹, David ROCHE¹, Zoé ROUY¹, Claudine MÉDIGUE¹
CEA/Genoscope/LABGeM, Université d'Evry Val-d'Essonne, CNRS-UMR 8030, Université Paris-Saclay, 2 rue Gaston Crémieux, 91057, Evry, France

Corresponding author: vallenet@genoscope.cns.fr

## 1 Introduction

The analysis of genomes from NGS platforms needs to be automated and fully integrated. However, maintaining consistency and accuracy in annotation is a challenging task because millions of proteins in databases are not assigned reliable functions.

The LABGeM team of the Genoscope sequencing center focuses its research activities on the development and application of new methods for genome analysis. These tools are then made available through MicroScope (http://www.genoscope.cns.fr/agc/microscope) [1], an integrated platform dedicated to microbial genome annotation and comparative analysis, which is being developed in our group since 2004.

## 2 Methods

The resource provides data from complete and ongoing genome projects together with post-genomic experiments (i.e. transcriptomics, re-sequencing of evolved strains) allowing users to improve the understanding of gene functions. We will present an overview of the MicroScope analysis pipelines and illustrate the use of several new functionalities in the context of data discovery and expert annotation, which concern:

- comparative genomics with synteny computations and pan-genome analyses,

- the prediction of virulence and antimicrobial resistance genes,

- the detection and annotation of genomic regions of interest, like, secretion systems, integrons and secondary metabolite biosynthesis gene clusters,

- and metabolic network reconstruction assisted by the GROOLS expert system (https://github.com/grools) [2].

## 3 Conclusions

To date, MicroScope contains data for about 7,000 microbial genomes, part of which are manually curated and maintained by microbiologists ( > 3,200 personal accounts in March 2017). The platform enables collaborative work in a rich comparative genomic context and improves community-based curation efforts.

## References

[1] David Vallenet et al. Microscope in 2017: an expanding and evolving integrated resource for community expertise of microbial genomes. *Nucleic Acids Res.*, (45(D1)):D517–D528, 2017.

[2] Jonathan Mercier et al. Grools: reactive graph reasoning for genome annotation through biological processes. *bioRxiv*, 2017.
NETSYN : NETwork SYNteny, a new tool to help functional annotation

Benjamin VIART<fnref n="1" />, Karine BASTARD<fnref n="1" />, Guillaume REBOUL<fnref n="2" />, z, Hugo PEREIRA<fnref n="1" />, Remi PLANEL<fnref n="1" />, Mark STAM<fnref n="1" />, Claudine MEDIGUE<fnref n="1" /> and David VALLENET<fnref n="1" />

Corresponding author: bviart@genoscope.cns.fr

# 1 Abstract

About 25% of the protein families collected in the Pfam database are annotated with unknown functions [1]. Besides, in databanks in which proteins are mainly annotated based on sequence similarity only (e.g. TrEMBL, NR), the level of misannotation reaches about 60% for superfamilies [2]. To detect possible errors or to suggest hypothetical functions for these proteins, we developed a bioinformatics approach, called NetSyn, based on genomic neighborhood.

NetSyn aims at exploring conserved genomic contexts (i.e. syntenies) to classify proteins within a family. It consists of (i) a pairwise comparison of gene chromosomal organization to find conserved syntenies, (ii) the representation of the results using a non-oriented weighted graph where nodes represent proteins and edges the average score of the number of genes involved in the synteny, and (iii) nodes are clustered using the weighted Markov Cluster Algorithm [3]. Each cluster is supposed to gather iso-functional proteins. The software can be used either to study multi-functional enzyme families or to find interactions between several families.

As a proof of concept, NetSyn was tested on two families. In the β-keto acid cleavage enzymes family, 7 profiles of active sites have been described, each one has been linked experimentally with specific enzymatic activities [4]. NetSyn was able to retrieve the proteins associated with the 7 known activities and to disclose 4 other clusters of proteins that might have new functions. In the amidinotransferases family [5], proteins annotated as arginine deiminase were split in two different clusters. One of them probably contains proteins that might not be correctly annotated as these proteins do not share similar genomic context with characterized arginine deiminase. Our method points out another cluster in which half the proteins are annotated as "Unknown function". Because they are highly connected with enzymes experimentally characterized as dimethylarginine dimethylaminohydrolase, these proteins might share the same enzymatic activities. The program is available to the scientific community on GitHub.

# References

[1] R. Mudgal, S. Sandhya, N. Chandra, and N. Srinivasan. De-DUFing the DUFs: Deciphering distant evolutionary relationships of Domains of Unknown Function using sensitive homology detection methods. *Biol. Direct*, 10:38, Jul 2015.

[2] A. M. Schnoes, S. D. Brown, I. Dodevski, and P. C. Babbitt. Annotation error in public databases: misannotation of molecular function in enzyme superfamilies. *PLoS Comput. Biol.*, 5(12):e1000605, Dec 2009.

[3] S. van Dongen. A cluster algorithm for graphs. *Technical Report INS R0010, National Research Institute for Mathematics and Computer Science in the Netherlands*, May 2000.

[4] K. Bastard, A. A. Smith, C. Vergne-Vaxelaire, A. Perret, A. Zaparucha, R. De Melo-Minardi, A. Mariage, M. Boutard, A. Debard, C. Lechaplais, C. Pelle, V. Pellouin, N. Perchat, J. L. Petit, A. Kreimeyer, C. Medigue, J. Weissenbach, F. Artiguenave, V. De Berardinis, D. Vallenet, and M. Salanoubat. Revealing the hidden functional diversity of an enzyme family. *Nat. Chem. Biol.*, 10(1):42–49, Jan 2014.

[5] H. Shirai, T. L. Blundell, and K. Mizuguchi. A novel superfamily of enzymes that catalyze the modification of guanidino groups. *Trends Biochem. Sci.*, 26(8):465–468, Aug 2001.

<footnote n="1">¹ LABGeM, CEA, Genoscope, Université d'Evry, CNRS-UMR8030, Université Paris-Saclay, 2 rue Gaston Cremieux, 91057 Evry, France</footnote>

<footnote n="2">² Unité d'écologie systématique et évolution, CNRS UMR 8079, Université Paris Sud</footnote>
# Stratégie d'assemblage de génomes en cellule unique de protistes marins dans le cadre du projet Tara Oceans

Léo D'AGATA<fnref n="1" />, Yoann SEELEUTHNER<fnref n="1" />, Julie POULAIN<fnref n="1" />, Patrick WINCKER<fnref n="2" />, and Jean-Marc AURY<fnref n="1" />

Corresponding Author: ldagata@genoscope.cns.fr

Le plancton des océans est responsable de plus de la moitié de l'oxygène produit sur la planète. Il n'en demeure pas moins très peu étudié et l'on connait encore mal l'impact que pourrait avoir le réchauffement climatique et la pollution sur ce dernier. Son rôle dans l'absorption du CO2 et la régulation du climat pourrait en effet être directement affecté. C'est dans le but de répondre à ces questions que l'expédition Tara, à laquelle le Genoscope participe, fut lancée entre 2009 et 2013 sous le nom de *Tara Oceans* [1].

Un des objectifs est l'étude précise des génomes de protistes marins. Ces derniers n'étant pour la plupart pas cultivables en laboratoire, la solution mise en place fut de recourir au séquençage en cellule unique. Le protocole de ce dernier consistait à extraire l'ADN d'une seule cellule, à l'amplifier par MDA [2] et enfin à le séquencer. Les étapes d'extraction et d'amplification génèrent toutefois de nombreux problèmes rendant l'assemblage de leur génome difficile. Certaines régions peuvent en effet être non capturées lors de l'extraction de l'ADN tandis que d'autres ne sont pas amplifiées. L'amplification appliquée à ces protistes génère de plus une couverture verticale très irrégulière.

Des outils adaptés sont donc nécessaires pour assembler ces génomes. Un workflow spécifique a ainsi été mis en place au Genoscope pour traiter ces données. Une idée proposée par ce workflow consiste notamment à se servir de la synergie de plusieurs cellules pour résoudre au mieux les problèmes de régions non couvertes. Cette technique dite de « co-assemblage » nécessite néanmoins de savoir à l'avance si ces cellules appartiennent au même organisme. L'ADN ribosomique 18S (gène marqueur chez les eucaryotes) peut apporter un début de réponse. Cette séquence n'est cependant pas toujours très résolutive et ne reflète pas totalement la diversité des génomes.

Pour répondre à cette problématique nous avons utilisé différentes méthodes et outils associés. Deux méthodes principales ont été retenues. La première consiste à se servir de la technique des graphs de De Bruijn colorés sur un premier co-assemblage [3]. Chaque couleur correspondant à une cellule, il nous est possible d'observer par quelles cellules sont recouverts les contigs. Les régions génomiques partagées entre les cellules nous informent ainsi sur leur proximité ou leur éloignement. La deuxième méthode consiste quant à elle à découper les lectures en mots chevauchants de longueur k (k-mers) et à mesurer la proportion de ces mots qui sont partagés entre deux cellules séquencées de façon indépendante [4]. Plus le taux de k-mers partagés par deux cellules est élevé, plus on peut supposer que celles-ci appartiennent à la même espèce. L'utilisation conjointe de ces méthodes nous aide donc à choisir les cellules pouvant être ou non co-assemblées.

## References

<footnote n="1">¹ Commissariat à l'Energie Atomique (CEA), Institut de Génomique (IG), Genoscope, F-92057 Evry, France.</footnote>

<footnote n="2">² Genoscope, CEA, CNRS UMR 8030, Université d'Evry, France.</footnote>

<footnote>[1] Eric Karsenti, Silvia G Acinas, et al. A holistic approach to marine eco-systems biology. *PLoS Biol* 9(10): e1001177, 2011.</footnote>

<footnote>[2] Frank B Dean, Seiyu Hosono, et al. Comprehensive human genome amplification using multiple displacement amplification. *Proc Natl Acad Sci USA* 99(8):5261-5266, 2002.</footnote>

<footnote>[3] Narjes S Movahedi, Mallory Embree, et al. Efficient Synergistic Single-Cell Genome Assembly. *Front Bioeng Biotechnol* 4: 42, 2016.</footnote>

<footnote>[4] Gaëtan Benoit, Pierre Peterlongo, et al. Multiple comparative metagenomics using multiset k-mer counting. *PeerJ Computer Science* 2:e94, 2016.</footnote>
Assessing the functional impact of genomic alterations using proteogenomics

Georges Bedran¹,²,³, Yves Vandenbrouck¹, Eric Bonnet², Jean-François Deleuze²,
Delphine Pflieger¹,⁴* and Christophe Battail¹,²*

¹ CEA/DRF/BIG/BGE/EDYP, 17 rue des Martyrs, 38054, Grenoble, France

² CEA/DRF/JACOB/CNRGH, 2 rue Gaston Crémieux, 91057, Evry, France

³ UFR Sciences et Techniques - Université de Rouen, Place Emile Blondel, 76821, Mont-Saint-Aignan, France

⁴ CNRS, 25 rue des Martyrs, Grenoble, 38042, France

* Co-senior authors

Corresponding Author: georges-bedran@etu.univ-rouen.fr

Proteomic data is obtained using a combination of liquid chromatography and tandem mass spectrometry (MS/MS) where peptides are most commonly identified by matching MS/MS spectra against theoretical spectra of all candidate peptides represented in a generalist protein sequence reference database[1]. The limitation of this approach is that variant peptides are missing from this reference database, thus they cannot be detected.

Proteogenomics[2] is an alternative approach where the reference database is replaced by a customized protein sequence database generated using genomic information extracted from RNA-seq data[3]. This strategy allows the detection of novel SAAVs, INDELs and splice junctions.

Using exploratory proteogenomics on a colorectal cancer cell line, we identified over 134 Single Amino Acid Variants (SAAV), from which 88 were also found by CPTAC colorectal cancer study[4]. However, in such analysis a significant number of alterations found by RNA-seq are not detected at the protein level (3605 non-synonymous mutations); therefore, we are aiming to address them using a targeted approach [5].

Our perspective is to develop a bioinformatic methodology capable of assessing the impact of the genetic alterations on the proteome in a cancer context.

Keywords

Proteogenomics, RNA-sequencing, colorectal cancer, SAAVs (single amino acid variants), splice junctions, INDELs (insertions / deletions), exploratory proteomics, targeted proteomics.

References

1. Nesvizhskii AI. A survey of computational methods and error rate estimation procedures for peptide and protein identification in shotgun proteomics. J. Proteomics. 2010;73:2092-2123.

2. Nesvizhskii AI. Proteogenomics: concepts, applications and computational strategies. Nat Meth. 2014;11:1114-25.

3. Wang X, Zhang B. customProDB: an R package to generate customized protein databases from RNA-Seq data for proteomics search. Bioinformatics [Internet]. 2013. Available from:
http://bioinformatics.oxfordjournals.org/content/early/2013/09/20/bioinformatics.btt543 abstract

4. Zhang B, Wang J, Wang X, Zhu J, Liu Q, Shi Z, et al. Proteogenomic characterization of human colon and rectal cancer. Nature. 2014;513:382-7.

5. Yeom J, Kabir MH, Lim B, Ahn H-S, Kim S-Y, Lee C. A proteogenomic approach for protein-level evidence of genomic variants in cancer cells. Sci. Rep. 2016;6:35305.
# Deciphering early cell fate decision by Single Cell RNA-Seq and DGE-Seq

Dimitri MEISTERMANN¹,², Yohann LELIEVRE², Eric CHARPENTIER³, Stéphanie KILENS¹,
Thomas FREOUR¹,⁴, Jérémie BOURDON² and Laurent DAVID¹,⁵

¹ CRTI, INSERM UMR1064, Université de Nantes, F-44000, Nantes, France

² LS2N, UMR 6004, Université de Nantes, F-44000, Nantes, France

³ l'institut du thorax, INSERM UMR1087, CNRS UMR6291, Université de Nantes, F-44000, Nantes, France

⁴ CHU de Nantes, Service de Biologie de la Reproduction, F-44000, Nantes, France

⁵ iPSC core facility, INSERM UMS 016, SFR Francois Bonamy, F-44000, Nantes, France;
Corresponding author: dimitri.meistermann@univ-nantes.fr

The major goal of our laboratory is to understand and quantify human preimplantation development, from fertilization to implantation in the uterus, and to predict its outcome. Specifically, we aim at deciphering the molecular mechanism driving cell fate during this first step of our existence. Understanding human preimplantation is therefore critical to improve assisted reproductive technologies (ART) and broaden the use of human pluripotent stem cells in regenerative medicine. It is during this timeframe that embryonic cells make their first choice of cellular fate, moving from one totipotent cell in the zygote to an embryo stratified by three cell types in the mature blastocyst. Moreover, the main objective of in vitro fertilization is to support the development of the zygote into a blastocyst, before transfer in infertile patients. To discover how early cell fate specification is regulated, our team develops several strategies to model embryos. To reach this goal, we employed two transcriptomics approaches: single-cell RNA-Seq and DGE RNA-Seq.

Firstly, we used *Single-Cell RNA-Seq* to identify unknown subpopulation of preimplantation embryo cells and associated genes network. For this purpose, we designed a single-cell analysis pipeline written with *snakemake*, a new and handy workflow tool. Biological data in single-cell are highly heterogenous and require unsupervised approaches. Thus, the pipeline answer to these constraints and includes complete analysis from FASTQ to identification of cell subpopulations with dimension reduction methods. We also analyzed cell trajectories with pseudo-time methods (*Slingshot*, *Monocle2* [2]) to pinpoint critical moment of cell fate decisions. Then, we determined specific gene networks associated with subpopulations, embryo lineage or embryo stage with Weighted Gene Correlation Analysis (*WGCNA* [3]). We observed that *WGCNA* is a particularly suitable tool for single-cell analysis, giving coherent results from noisy and heterogenous data. Here we will present the pipeline and show how it has solved the problem of heterogeneity and lack of biological material inherent to human embryo. Secondly, a novel, cost-efficient RNA-Seq method, digital expression RNA-Seq (*DGE-Seq*) was employed to sequence human induced naive pluripotent stem cells (*hiNPSC*), a recently characterized state of pluripotency in human [4]. As a main result of this work, we show that hiNPSC have a more similar metabolism to preimplantation embryo than usual hiPS cells. This makes hiNPSC an appropriate model of preimplantation development.

Here we present methods of computational biology that guided our research strategy, and promoted the use of transcriptomics in Nantes. Both approaches contributed significantly to our understanding of preimplantation development, opening new avenues of research in the fields of ART and regenerative medicine fields.

## Acknowledgements

This work was supported by the BiRD bioinformatics facility.

## References

[1] J. Rossant et P. P. L. Tam, « New Insights into Early Human Development: Lessons for Stem Cell Derivation and Differentiation », *Cell Stem Cell*, vol. 20, n° 1, p. 18-28, janv. 2017.

[2] C. Trapnell et al., « The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells », *Nat. Biotechnol.*, vol. 32, n° 4, p. 381-386, mars 2014.

[3] P. Langfelder et S. Horvath, « WGCNA: an R package for weighted correlation network analysis », *BMC Bioinformatics*, vol. 9, p. 559, 2008.

[4] G. Guo et al., « Naive Pluripotent Stem Cells Derived Directly from Isolated Cells of the Human Inner Cell Mass », *Stem Cell Rep.*, vol. 6, n° 4, p. 437-446, avr. 2016.
# In search of W sex chromosome-specific sequences in the genome of the isopod
crustacean *Armadillidium vulgare*

Mohamed Amine CHEBBI, Thomas BECKING, Isabelle GIRAUD, Bouziane MOUMEN,

Clément GILBERT, Jean PECCOUD and Richard CORDAUX

Laboratoire Ecologie et Biologie des interactions UMR CNRS 7267, Equipe Ecologie Evolution
Symbiose, Université de Poitiers, Bât. B8, 5 rue Albert Turpin, TSA 51106, 86073 Poitiers Cedex
9 France

Corresponding Author: mohamed.amine.chebbi@univ-poitiers.fr

## Abstract

In the isopod crustacean *Armadillidium vulgare*, genetic sex determination follows female heterogamety (ZZ males and ZW females). Z and W sex chromosomes show no apparent heteromorphy. WW female individuals are viable and genetic sex determination in the closely related species *Armadillidium nasatum* follows male heterogamety (XY males and XX females). These observations suggest that the evolution of *A. vulgare* sex chromosomes is at an incipient stage of the specialization of a pair of ancestral autosomes carrying sex determinants. To test this hypothesis and identify W-specific sequences, we sequenced and assembled the genome of *A. vulgare* by combining Illumina and PacBio sequencing technologies. The 1.7 Gb genome of *A. vulgare* is highly repeated, as simple repeats and transposable elements represent ~70% of the genome. To identify sex linked sequences, we obtained male and female Illumina sequencing reads independently and used two different bioinformatics approaches: Chromosome Quotient (CQ) and Y Genome screening (YGS) methods. These analyses confirmed the homomorphy of *A. vulgare* sex chromosomes and allowed us to identify candidate W-specific scaffolds.

## Keywords:

female heterogamety, homomorphy, W-specific sequences, hybrid de novo genome assembly, highly repeated, Chromosome Quotient (CQ), Y Genome screening (YGS)
# Exome sequencing to identify the molecular mechanism underlying chordomas pathogenesis

Zakia Tariq<fnref n="1" />, Keltouma Driouch<fnref n="1" />, Virginie Bernard<fnref n="2" />, Ivan Bieche<fnref n="1" />, Virginie Raynal<fnref n="2" />, Sylvain Baulande<fnref n="2" />, Hamid Mammar<fnref n="3" />, Julien Masliah Planchon¹⁴

Corresponding Author: zakia.tariq@curie.fr

Chordomas are rare bone tumors, often observed in cranial, spinal and sacral sites. With a slow evolution, these tumors are clinically diagnosed at a late stage, decreasing patients’ life expectancy [1]. However, these tumors remain poorly described. This project aims to identify and characterize somatic alterations involved in this pathology. To this end, a series of tumors were analyzed by means of high-throughput exome sequencing.

On Illumina platform, paired-end 100x100 exome-Seq was performed for eight primary tumors and matching germline DNAs obtained from blood samples. After quality control and mapping files cleaning, somatic variants search was performed. Therefore, three variant callers were used in parallel to increase the sensitivity and the confidence of predicted variations (SNV). After Annovar annotation, variants were considered as somatic when absent in the germline samples and if their reported frequency in 1000Genome database was below 0.1%. Finally, variants were validated by IGV visualization discarding sequencing and alignment errors. In another hand, copy number variations (CNV) were predicted using Facets and Sequenza tools, to determine putative copy number alterations.

The number of identified somatic mutations ranged from 12 to 32 per sample, except one tumor exhibiting 166 somatic variants. Among the mutated genes, known drivers, such as *KIT* and *PIK3CA*, were detected. Mutated genes highlighted recurrent pathways involved in biologic processes such as chromatin organization, and epigenetic modifications. With CNV analysis, we identified a chromothripsis phenomenon as formerly described [2]. In our study, such events were observed in 2 out of 8 tumors. Furthermore, we characterized a high frequency of *CDKN2A/B* homozygous deletion in 80% of samples; a much higher rate than previously reported [3]. Strikingly, the most extensively mutated sample, presented a homozygous deletion encompassing *MLH1* gene, suggesting a DNA mismatch repair deficiency in this tumor.

The integrated analysis of SNV and CNV profiles emphasized several molecular pathways involved in chordomas. In particular, the *MEFC2* gene, regulating bone differentiation, showed a biallelic inactivation i.e. a stop gain variation coupled with a loss of heterozygosity suggesting a complete loss of the protein function.

In conclusion, we observed different variations already described in the literature, but we also characterized novel putative genes involved in chordoma, not linked to the pathology yet. To validate our hypotheses, exome-Seq will be performed on a new sample series and RNA-Seq data will be generated to complete the study.

## References

[1] Walcott BP, et al. Chordoma: Current concepts, management, and future directions. *The Lancet Oncology*, 13(2):e69-76, 2012.

[2] Stephens PJ, et al. Massive Genomic Rearrangement Acquired in a Single Catastrophic Event during Cancer Development. *Cell*, 144(1):27-40, 2011.

[3] Hallor KH, et al. Frequent deletion of the CDKN2A locus in chordoma: analysis of chromosomal imbalances using array comparative genomic hybridisation. *Br J Cancer*, 98: 434-42, 2008.

<footnote n="1">¹ Service de génétique, Institut Curie, 26 rue d’Ulm 75248, Paris, France.</footnote>

<footnote n="2">² Next generation sequencing platform, ICGex, Institut Curie, Paris, France.</footnote>

<footnote n="3">³ Service radiothérapie, Institut Curie, Paris, France.</footnote>

<footnote n="4">⁴ INSERMU830, Département biologie des tumeurs, Institut Curie, Paris, France.</footnote>
Revisiting cell lineage specification during male sex determination with single-cell RNA sequencing

Isabelle STÉVANT<sup>1,2,3</sup>, Yasmine NEIRJINCK<sup>1</sup>, Christelle BOREL<sup>1</sup>, Jessica ESCOFFIER<sup>1</sup>, Lee B. SMITH<sup>4,5</sup>, Stylianos E. ANTONARAKIS<sup>1,2</sup>, Emmanouil T. DERMITZAKIS<sup>1,2,3</sup>, Serge NEF<sup>1,2</sup>

<sup>1</sup> Department of Genetic Medicine and Development, University of Geneva, 1211 Geneva, Switzerland;

<sup>2</sup> iGE3, Institute of Genetics and Genomics of Geneva, University of Geneva, 1211 Geneva, Switzerland;

<sup>3</sup> SIB, Swiss Institute of Bioinformatics, University of Geneva, 1211 Geneva, Switzerland;

<sup>4</sup> MRC Centre for Reproductive Health, University of Edinburgh, Edinburgh EH16 4TJ, UK;

<sup>5</sup> School of Environmental and Life Sciences, University of Newcastle, Callaghan, NSW 2308, Australia.

Corresponding Author: isabelle.stevant@unige.ch

Among the genetic characteristics that influence our physical identity all along our life, the most important is surely our sex. The ultimate goal of male and female sex differentiation and development is to provide organisms the necessary attributes for sexual reproduction. In mammals, the genetic sex of individuals is cast at fertilisation with the bring of an X or a Y chromosome from a spermatozoa to a bearing X oocyte. Phenotypic sex is only revealed during foetal development, when the gonads start to differentiate as ovaries or testes (around 6 embryonic weeks in human, and at embryonic day 11.5 in mouse). Consequently to the ovary or testis differentiation, the whole embryo will adopt the secondary sexual characteristics such as male and female reproductive tracts and appendixes.

Current knowledge of sex determination was built on gene by gene knock-in/knock-out in mice and transcriptomic analysis on pool of purified cells. The first method suffers from a very low throughput and produces a partial knowledge of the function of one or a limited set of genes. The second method caveat resides in the lack of specific reporters for the different cell types within the gonad and makes difficult the appreciation of the cell type heterogeneity, especially before sex determination. Moreover, a pool of purified cells results in an average message of non synchronous differentiating cells and thus blur the chronology of gene expressions.

In this study, we used transgenic mice carrying the Nr5a1-GFP reporter to isolate by FACS the gonadal somatic cells at five key stages of sex determination and gonad development in male. We proceeded to the RNA-sequencing of 391 single cells and identified the different cell types present in the developing testis and reconstructed their cell lineages. At E10.5, we detected one homogeneous progenitor cell population of NR5A1⁺ cells expressing epithelial and stem cell marker genes, consistently with their bi-potential state before sex determination. From E11.5, one fraction of these cells activated a strong genetic program and initiated their differentiation as Sertoli cells and the other fraction of progenitor cells evolved transcriptionally and progressively start to express steroidogenic lineage markers such as *Arx* and *Pdgfra* and are restricted to the interstitial compartment of the testis. From E12.5, we detected few fetal Leydig cells differentiating from the interstitial progenitors.

This study represents the most granular transcriptomic database of gonadal somatic cells during early testis development. With these data, we contribute in the deeper understanding of the cell differentiation during sex determination.
# How to analyse human or mouse Genome-scale CRISPR Knock-Out (GeCKO) datasets ?

Marc DELOGER<fnref n="1" />,*, Pierre GESTRAUD<fnref n="1" />,*, Raphaël MARGUERON<fnref n="2" /> and Nicolas SERVANT<fnref n="1" />

Corresponding Author: marc.deloger@curie.fr

## 1 Introduction

Cas9 (Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR) associated protein 9) target specificity is determined by a short 20bp sequence of a single guide RNA (sgRNA). Consequently, large scale oligo synthesis of guide sequences suggests a new way to interrogate gene function at a genome-wide scale. It has been shown that infecting cell with Cas9:sgRNAs libraries using lentiviruses can facilitate both positive and negative loss-of-function screening in mammalian cells [1, 2]. When Cas9-induced double-strand breaks are introduced into coding sequences, error-prone repair machinery will often introduce « indel » potentially leading to loss-of-function allele [3]. In order to check the Genome-scale CRISPR Knock-Out (GeCKO) screening modification efficiency, off-target modification rate, consistency between unique sgRNAs targeting the same genes and validation rate of screen hits, we decided to propose a full R package that goes from FASTQ files to an HTML report giving all quality control and statistical analysis results of your GeCKO experiments taking into account process data difficulties such as results ranking.

## 2 Data pre-processing, Quality Control and Statistical analysis

Original human GeCKOv2 library is composed of 123'411 sgRNAs (~ 6 per gene + 2000 negative controls). First, we propose to clean the library by merging redundant sequences and removing the ones that have multi-hits on hg38 genome, obtaining 113'761 guides. These sequences were detected (with or without mismatches) and counted on each sample sequenced reads (FASTQ files) in order to generate a count table as output for further analysis. Nine criteria are used to assess the experiment quality. An R markdown notebook is automatically generated with all the corresponding values/figures.

From the count table, data are normalised using classical RNAseq methods from R package edgeR. We propose to compute the normalisation factor either on the whole dataset or on the non-targeting guides only. Differential abundance of guides is estimated with the limma/voom framework. From the moderated t-statistics estimated by limma, we derived a one-sided p-value to test explicitly the depletion or the enrichment of guides. A gene-level score is then computed by the RRA method [4] and p-values are derived from a null distribution estimated from random genes sampled into non-targeting guides scores.

## Acknowledgements

This project was done thanks to SiRIC-Curie funding and Raphaël MARGUERON's datasets

## References

[1] Shalem O, Sanjana NE, Hartenian E, Shi X, Scott DA, Mikkelsen TS, Heckl D, Ebert BL, Root DE, Doench JG, Zhang F. Genome-scale CRISPR-Cas9 knockout screening in human cells. Science. 2014 Jan 3;343(6166):84-7.

[2] Sanjana NE, Shalem O, Zhang F. Improved vectors and genome-wide libraries for CRISPR screening. Nat Methods. 2014 Aug;11(8):783-4.

[3] Rodgers K, McVey M. Error-Prone Repair of DNA Double-Strand Breaks. J Cell Physiol. 2016 Jan;231(1):15-24.

[4] Kolde R, Laur S, Adler P, Vilo J: Robust rank aggregation for gene list integration and meta-analysis. Bioinformatics. 2012, 28: 573-580. 10.1093/bioinformatics/btr709

<footnote n="1">¹ Institut Curie, PSL Research University, INSERM, U 900, MINES, ParisTech, F-75005, Paris, France</footnote>

<footnote n="2">² Institut Curie, PSL Research University, CNRS, UMR 3215, INSERM, U 934, F-75005, Paris, France</footnote>

<footnote>* Equally contributed</footnote>
Functional characterization of human epidermis co-expression modules

Gwenaelle LEMOINE¹, Marie Pier SCOTT-BOYER¹, Mickael LECLERQ¹, Arnaud DROIT¹

¹ Computational Biology Laboratory - CHU de Québec - Université Laval, 2705 Boulevard Laurier, G1V 4G2, Québec, Canada

Corresponding Author: lemoine.gwenaelle@gmail.com

Underlying molecular mechanisms in human tissues tend to be more documented everyday for a certain number of them [1]. Usually studied in pathological context, deepening the understanding of molecular mechanisms lead to a better prevention and treatment strategies of related diseases. Still, many (patho)physiological reactions remain unpredictable [2,3] if solely observed at the molecular level.

However, the recent surge of high-throughput omics technologies (e.g. transcriptomics, proteomics, metabolomics) allowed the emergence of new discoveries by system-wide analyses. A popular approach to analyze these large data is the modeling of gene similarity based networks. Also called co-expression networks, they group highly dependent genes by modules. Those modules can mostly be related to and used to further understanding of phenotypic traits information. Besides summarizing gene expression, this method allow to functionally annotate unknown genes.

Here, we emphasis on epidermal data which remains one of the easiest tissues to explore because of its ease of access for sampling (unlike tissues from internal organs). Helped by the R package WGCNA, first step of the pipeline involves the network creation. It is based on previous defined similarity and a threshold consisting of numerical set value (called hard-thresholding), or a function with a defined parameter value (called soft-thresholding). Module detection inside the network is then performed in order to group highly co-expressed genes. This process is then repeated on several pertinent data sets in order to identify epidermis-specific modules and compare them against each other.

Because those modules are more likely to mediate phenotypic traits of interest [4], a further investigation aims at comparing the activity of epidermis-specific modules under different conditions. On one hand, skin-linked factors like UV exposure, skin diseases or scaring will be studied. On the other hand, we will look into more general element such as sex or age impact. Most data comes from available public databases, essentially GTEx [5] or GEO [6], and internal sequencing data of *in-vivo* skin samples will complete them. Finally, the main frame is to lead to a refining of module functional characterization, and therefore allow some supervised learning for ulterior prediction models.

Since transcriptomic impacts other omics stratum in many ways, the incoming challenge is the association of those outlined epidermis-specific gene modules with other biological networks such as protein-protein interactions networks or metabolic pathways in the same tissue [7].

References

[1] Glass et al. *Gene expression changes with age in skin, adipose tissue, blood and brain*. Genome Biology, 2013

[2] Sudharsana Sundarrajan, Mohanapriya Arumugam. *Weighted gene co-expression based biomarker discovery for psoriasis detection*. GeneVolume 593, Issue 1, 15 November 2016, Pages 225–234

[3] Kuehne et al. *An integrative metabolomics and transcriptomics study to identify metabolic alterations in aged skin of humans in vivo Andreas*. BMC Genomics, 2017

[4] Sipko van Dam, Urmo Võsa, Adriaan van der Graaf, Lude Franke, João Pedro de Magalhães. *Gene co-expression analysis for functional classification and gene-disease predictions*. Brief Bioinform, 2017

[5] The GTEx Consortium. *The Genotype-Tissue Expression (GTEx) project*. Nature Genetics 45, 580–585, 2013

[6] Edgar R, Domrachev M, Lash AE. *Gene Expression Omnibus: NCBI gene expression and hybridization array data repository*. Nucleic Acids Res., 2007, 34, 2002

[7] Serin Elise A. R., Nijveen Harm, Hilhorst Henk W. M., Ligterink Wilco. *Learning from Co-expression Networks: Possibilities and Challenges*. Frontiers in Plant Science, 2015
# Integrated and flexible analysis of scRNA-seq data with Eoulsan scRNA-seq

Geoffray BRELURUT<fnref n="1" />*, Nathalie LEHMANN<fnref n="1" />*, Céline HERNANDEZ<fnref n="1" />, Morgane THOMAS-CHOLLIER<fnref n="1" />, Denis THIEFFRY<fnref n="1" />, Stéphane LE CROM<fnref n="2" /> and Laurent JOURDREN<fnref n="2" />

Corresponding Author: jourdren@biologie.ens.fr

Keywords: Single-Cell RNA sequencing, distributed computing, high-throughput sequencing, pipeline.

## Summary

Single cell RNA sequencing (scRNA-seq) provides new opportunities to characterize cell transcriptomic heterogeneity, e.g. the discovery of new cell populations and the reconstruction of cell lineages. Sequence reads are produced from individual cells using two main approaches: full transcript amplification or partial amplification with single molecule tagging using Unique Molecular Identifiers (UMI) [1]. These different protocols combined with various computational analyses result into numerous possible computational workflows.

Here we present a new workflow dedicated to scRNA-seq data analysis. Still under development, this workflow is based on the *Eoulsan* framework [2], starts with raw FASTQ files, and encompasses various quality control, normalisation, read mapping, and more advanced steps.

The first steps of the analysis (from read filtering to expression estimation) can be distributed on computer clusters to greatly reduce processing time. Most of the popular frameworks or job schedulers like *Hadoop*, *TORQUE* or *HTCordor* are currently supported by *Eoulsan*. Noteworthy, we take full advantage of the *Docker-Galaxy* framework introduced in *Eoulsan 2* [3] to propose a series of modules allowing the analysis of both read-based and UMI-based data. Downstream steps include a wide range of tools for gene differential expression (*SCDE* [4], *SCDD* [5]), cell clustering and lineage reconstruction (*Monocle 2* [6], *Destiny 2* [7]).

The experimental design ofr an *Eoulsan scRNA-seq* workflow is stored in simple text file, while parameters are stored in a second xml file, ensuring flexibility and traceability. Keeping track of each successful step, this approach allows to swiftly resume large analyses upon trouble-shooting, and further ensure reproducibility.

In conclusion, *Eoulsan scRNA-seq* provides an integrated workflow for scRNA-seq data analysis on standalone workstations or on computer clusters. With its modular structure and distributed data processing, it can handle large amounts of data in a reproducible, yet flexible manner.

## References

[1] Christoph Ziegenhain et al. Comparative Analysis of Single-Cell RNA Sequencing Methods. *Molecular Cell*, 65: 631-43, 2017.

[2] Laurent Jourdren et al. Eoulsan: A Cloud Computing-Based Framework Facilitating High Throughput Sequencing Analyses. *Bioinformatic*, 28:1542-3, 2012.

[3] http://www.outils.genomique.biologie.ens.fr/eoulsan2/

[4] Peter V. Kharchenko, Lev Silberstein, and David T. Scadden. Bayesian approach to single-cell differential expression analysis. *Nature Methods*, 11:740-2, 2014.

[5] Keegan D. Korthauer et al. A statistical approach for identifying differential distributions in single-cell RNA-seq experiments. *Genome Biology*, 17:222, 2016.

[6] Cole Trapnell et al. Pseudo-temporal ordering of individual cells reveals dynamics and regulators of cell fate decisions. *Nature Biotechnology*, 32:381-6, 2014.

[7] Philipp Angerer et al. Destiny: diffusion maps for large-scale single-cell data in R. *Bioinformatics*, 32: 1241-3, 2015.

<footnote n="1">¹ École normale supérieure, PSL Research University, CNRS, Inserm, Institut de Biologie de l'École normale supérieure (IBENS), Equipe Biologie Computationnelle des Systèmes, 75005 Paris, France.</footnote>

<footnote n="2">² École normale supérieure, PSL Research University, CNRS, Inserm, Institut de Biologie de l'École normale supérieure (IBENS), Plateforme Génomique, 75005 Paris, France.</footnote>

<footnote>* These authors contributed equally to this work</footnote>
# Efficient data structure for indexing and similarity computation of nucleic sequences

Camille MARCHET<fnref n="1" />, Arnaud MENG<fnref n="2" />, Lolita LECOMPTE<fnref n="1" />, AntoineLIMASSET<fnref n="1" />, LucieBITTNER<fnref n="2" /> and Pierre PETERLONGO<fnref n="1" />

Corresponding author: camille.marchet@irisa.fr

## 1 Introduction

Sequencing of (meta-)genomes and (meta-)transcriptomes generates huge sets of sequences, that are often chopped in voluminous sets of k-mers for their further analysis. Such instances constitute challenges for high performance computation. To extract relevant pieces of information from the large data sets generated by current sequencing techniques, one must rely on extremely scalable methods and solutions. In this work we present a straightforward indexing structure called quasi-dictionary that scales to billions of elements and we propose direct applications based on k-mer diversity to explore (meta-)genomics/transcriptomics data sets.

## 2 Quasi-dictionary applications to sequencing data

The quasi-dictionary structure is tailored to handle huge quantities of k-mers to index. It relies on a MPHF library [1] combined with a fingerprint system. We present two applications of such a data structure for short reads data: SRC_COUNTER and SRC_LINKER. Both can be used to compare two samples or sample versus itself in particular when traditional methods do not scale. SRC_COUNTER links any read to its estimated abundance in a collection of samples. SRC_LINKER connects any read from a given sample to similar reads in the data set [2]. Furthermore we present an application of SRC_LINKER to long reads. Contrary to the previous applications, the input are long, erroneous reads from last generation of sequencing (PacBio, Nanopore).

## 3 Results

We present the quasi-dictionary data structure performances and its applications results and practical use cases, notably at work on marine data sets. We state the memory and time performances of the quasi-dictionary as well as the practical impact of its false positives. We also compare SRC_LINKER with state of the art tools for short and long reads for retrieving similarities between reads. As for short reads, we demonstrate our gain in scaling, that applies in particular in meta-genomics field. With long reads, we show the benefits of having a lightweight data structure to deal with high error rates with robustness, combined to scalability.

We highlight the fact that SRC_LINKER is one of the few tools that can be applied successfully to explore short as well as long reads sequencing data. Finally the data structure and its applications remain simple, which make them easily adaptable for further challenges.

## 4 Citations

### References

[1] Antoine Limasset, Guillaume Rizk, Rayan Chikhi, and Pierre Peterlongo. Fast and scalable minimal perfect hashing for massive key sets. *arXiv preprint arXiv:1702.03154*, 2017.

[2] Veronika B Dubinkina, Dmitry S Ischenko, Vladimir I Ulyantsev, Alexander V Tyakht, and Dmitry G Alexeev. Assessment of k-mer spectrum applicability for metagenomic dissimilarity analysis. *BMC Bioinformatics*, 17(1):38, dec 2016.

<footnote n="1">¹ IRISA Inria Rennes Bretagne Atlantique, GenScale team</footnote>

<footnote n="2">² Sorbonne Universités, UPMC Univ Paris 06, CNRS, Evolution Paris Seine - Institut de Biologie Paris Seine (EPS - IBPS), 75005 Paris, France</footnote>
# Detection of poly-adenylation sites from RNA-Seq data

Cyril FOURNIER and Anamaria NECSULEA

Univ Lyon, Laboratoire de Biométrie et Biologie Evolutive, Université Claude Bernard Lyon1, 69622 Villeurbanne cedex, France

Corresponding author: cyril.fournier@univ-lyon1.fr

## Submission

Transcriptome sequencing technologies (RNA-seq) are increasingly used to accurately and sensitively measure gene expression levels, across a wide variety of biological samples [1]. The utility of RNA-seq goes beyond the simple assessment of gene expression levels, as it allows detection and quantification of alternative transcript isoforms, as well as annotation of new genes and isoforms. Here, we propose a method to annotate cleavage and poly-adenylation sites of messenger RNAs and long non-coding RNAs [2,3] from RNA-Seq data.

Our computational pipeline comprises three steps:

* First, we efficiently search for polyA tails in RNA-Seq data, allowing for sequencing and biological errors. This step results in a complete set of reads with polyA tracts, which includes genuine 3' ends of transcripts but also encompasses polyA repeats present in the genome.

* Second, to filter out these repeats from the dataset and to enrich for genuine polyA tails, the reads are artificially cleaved to remove the polyA tail and the remainder is aligned on the reference genome. We then select the reads for which the polyA tract does not map on the genome, as expected given that polyA tails are added post-transcriptionally to messenger RNAs.

* Third, the predicted poly-adenylation and cleavage sites are assessed by checking the flanking nucleotides and the presence of a poly-adenylation signal, and finally compared to known cleavage sites.

The running time of this method with a sample of 100 millions reads of 100 bp is approximately of 20 minutes (test run on Ubuntu 16.04 LTS with 16Go RAM and i7-6700HQ CPU @ 2.60GHz). Thus, this method can be easily applied to the vast quantities of RNA-Seq data that are accumulating today. Our method is thus considerably faster than other programs that detect polyA sites from RNA-Seq data, such as *KLEAT* [4] and *ContextMap 2* [5].

We applied our method to a collection of RNA-Seq data, derived from four major organs (brain, kidney, liver and testis) and five developmental stages (from early organogenesis to adult and aged individuals), for mouse and rat. For a standard RNA-Seq sample, the fraction of reads with polyA tail detected is approximately 3/1000 reads. We annotated over 100,000 putative polyA sites in each species and we confirmed about 10,000 sites annotated in the Ensembl database and 3,000 sites using polyA-sequencing data, from *Derti et al* [6].

Ongoing analyses include comparative analyses of the predicted polyA sites across species, organs and developmental stages. We also plan to use our method to refine the 3' ends of long non-coding RNAs, whose annotation is currently challenging.

## References

[1] Zhong Wang, Mark Gerstein, and Michael Snyder. RNA-Seq: a revolutionary tool for transcriptomics. *Nature reviews Genetics*, 10, jan 2009.

[2] Diana F. Colgan and James L. Manley. Mechanism and regulation of mRNA polyadenylation. *Genes & Development*, 11:2755-2766, 1997.

[3] Bin Tian and Joel H. Graber. Signals for pre-mRNA cleavage and polyadenylation. *Wiley Interdisciplinary Reviews: RNA*, 3(3):385-396, 2012.

[4] Birol I, Raymond A, Chiu R, Nip KM, Jackman SD, and Kreitzman M et al. Kleat: cleavage site analysis of transcriptomes. *Pac Symp Biocomput*, 20:347-358, 2015.

[5] Thomas Bonfert and Caroline C. Friedel. Prediction of poly(a) sites by poly(a) read mapping. *PLOS ONE*, 12(1):1-32, 01 2017.

[6] Adnan Derti, Philip Garrett-Engele, Kenzie D. MacIsaac, Richard C. Stevens, Shreedharan Sriram, Ronghua Chen, Carol A. Rohl, Jason M. Johnson, and Tomas Babak. A quantitative atlas of polyadenylation in five mammals. *Genome Research*, 22(6):1173-1183, 2012.
# Développement d'une structure pour l'indexation et la compression de multi-génomes / New structure to index and compress multi-genomes

Clément AGRET<fnref n="1" />,<fnref n="2" />, Manuel RUIZ<fnref n="1" /> and Alban MANCHERON<fnref n="2" />

Corresponding Author: clement.agret@gmail.com

Alors qu'il a fallu treize ans et plusieurs millions de dollars pour séquencer et assembler les trois milliards de nucléotides composant le génome humain [1], un tel séquençage nécessite aujourd'hui seulement quelques jours et à peine quelques milliers de dollars. L'International Rice Research Institute (IRRI), dans le cadre du consortium GRISP (Global Rice Science Partnership), a initié un programme de séquençage de l'ensemble des variétés de riz, et aujourd'hui plus de 3000 génomes sont déjà disponibles [2].

Pour analyser et stocker efficacement cette importante masse de données, il est nécessaire de représenter les génomes dans une forme permettant leur consultation rapide, tout en économisant le plus possible l'espace nécessaire à leur stockage.

La structure de données qui fait actuellement le succès des méthodes d'analyse de génomes est le FM-index [3]. Il s'agit d'une structure compressée exploitant les propriétés de réorganisation des données de la Transformée de Burrows-Wheeler (BWT) [4] appliquées sur le génome à indexer. Cependant le FM-index n'est pas optimisé pour compressor une collection de génomes similaires. Ce qui veut dire que si l'on souhaite analyser les 3000 génomes de riz, on va devoir créer 3000 index pour pouvoir ensuite interroger chaque index l'un après l'autre.

Nous avons exploré et comparé les méthodes existantes (PanTools, TwoPaco, CHICO, etc.) permettant de construire efficacement un index commun aux 3000 génomes. Nous développons aussi une méthode basée sur un découpage des génomes par k-mers. Pour cela, nous avons dans un premier temps étudié l'évolution du nombre de k-mers communs entre différents génomes complets de riz. La représentation des ensembles de k-mers communs entre les différents génomes sous forme de diagrammes de Venn nous a confirmé la pertinence de notre approche. En raison du volume de données qui seront indexées, nous avons défini une représentation de l'ensemble des k-mers et de leurs présence/absence basée sur des structures succinctes. Les problématiques qui émergent sont d'une part la mise à jour dynamique de la structure, par exemple lorsqu'un génome est ajouté ou retiré de l'ensemble des génomes déjà indexés; et d'autre part la formulation et l'optimisation des requêtes que notre structure d'index doit permettre. Enfin, les problématiques liées à la représentation et l'exploration visuelle de l'index demeurent un sujet d'étude qu'il conviendra de traiter par la suite.

## References

[1] International Human Genome Sequencing Consortium: Initial sequencing and analysis of the human genome. Nature, 409(6822):860–921, 2001.

[2] Li JY, Wang J, Zeigler RS: The 3,000 rice genomes project: new opportunities and challenges for future rice research. Gigascience 3, 8, 2014.

[3] Paolo FERRAGINA et Giovanni MANZINI: Opportunistic Data Structures with Applications. In Proceedings of the 41st Annual Symposium on Foundations of Computer Science (FOCS), pages 390–398.

[4] Michael BURROWS et David WHEELER: A Block-Sorting Lossless Data Compression Algorithm. Rapport technique, 1994.

<footnote n="1">¹ UMR AGAP, CIRAD, Avenue Agropolis, 34398, Montpellier, FRANCE</footnote>

<footnote n="2">² LIRMM, Campus St Priest, 161 Rue Ada, 34090 Montpellier, FRANCE</footnote>
# Epigenetic marks and the human transcriptome diversity

Guillaume DEVAILLY¹, Anna MANTSONKI¹ and Anagha JOSHI¹

Developmental biology, The Roslin Institute, The University of Edinburgh, Easter Bush,
Midlothian, EH25 9RG, Scotland, UK

Corresponding author: guillaume.devailly@roslin.ed.ac.uk

Links between epigenetic marks and transcription regulation have been increasingly explored in recent years. Major technical advancements allowed the generation of massive, under-exploited datasets by consortia and individual laboratories. We present a systematic integrative re-analysis of Roadmap Epigenomics RNA-seq, WGBS, DNAseI and histone marks ChIP-seq data in 33 human cell types. We combined two complementary approaches to link epigenetic marks at location of interests (transcription start sites, transcription end sites, exons) to transcription features (transcription level, exon inclusion ratio, cryptic transcription start sites): 1- for all genes/exons within a cell, and 2- for every single gene/exon across all 33 cell types.

Results generated include the confirmation of many known relationships (i.e. promoter DNA methylation is negatively correlated to gene expression, while H3K4me3 at promoter is positively regulated to gene expression) at a large scale, but also the confirmation or rejection of less accepted correlations. For example, while we did find a slight decrease of DNA methylation density at non-included region within a cell, we could not find any trend linking changes of exon methylation to changes of inclusion ratio when carrying the analysis across cells. A web portal will be developed in the near future to allow easy exploration of the results.

## Acknowledgements

AJ is a Chancellor's fellow and AJ labs is supported by institute strategic funding from Biotechnology and Biological Sciences Research Council (BBSRC, BB/J004235/1). GD is funded by the People Programme (Marie Curie Actions FP7/2007-2013) under REA grant agreement No PCOFUND-GA-2012-600181.

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

**Fig. 1.** Promoter DNA methylation and gene expression level in pancreas. A. Gene expression level in all 44114 genes annotated by Gencode. First side bar indicates the gene type (yellow, protein coding genes), the second side bar indicates the 5 bins used in panel F-J (purple: highly expressed genes, green: lowly expressed genes). Protein coding genes are globally more expressed than other types of genes. B-E. Stacked profiles of CpG density (B), the mCpG ratio (mCpG/CpG) (C), the mCpG density (D), and the WGBS coverage (E) at promoter, sorted according to their expression level. Promoter CpG density is positively correlated with expression level, mCpG ratio and mCpG density are negatively correlated to expression level. Coverage of CpG rich region is lower than at CpG poor regions. F. Boxplot of gene expression level in each of the 5 bins defined in A, and used in G-J. G-J. Average profile of CpG density (G), mCpG ratio (H), mCpG density (I) and WGBS coverage (J) +/- SEM for each bin of promoters.
# Microbial diversity and plant cell wall-degrading enzyme dynamics during dew-retting of flax - one of the oldest applications of biotechnology to textile transformation

Christophe DJEMIEL¹, Sébastien GREC¹ and Simon HAWKINS¹

¹ UGSF - Unité de Glycobiologie Structurale et Fonctionnelle, Université de Lille, 59000, Lille, France

Corresponding Author: christophe.djemiel@gmail.com

The existence of spun and colored flax fibers dating from the Upper Paleolithic suggests that man has long exploited the biotechnological natural process known as dew-retting to extract these fibers from the flax plant for textile production [1,2]. This process is achieved directly on the soil surface of the field [3,4]. Despite many studies of this process to evaluate the degree of retting [5], relatively little is known about (i) the composition and the evolution of the microflora population during retting, (ii) the kinetics of the microbial communities colonizing the plant material and (iii) the composition and the evolution of Carbohydrate Active enZymes degrading plant cell walls [6].

To improve our understanding of dew-retting, we first used a metabarcoding approach to identify the membership and structure of the microbial communities (focusing on bacteria and fungi). This approach also allowed us to identify (i) some potential bacterial major enzymatic functions related to carbohydrate degradation based on functional prediction using PICRUST (http://picrust.github.io/picrust/) and (ii) a strong pattern of fungal trophic modes [7,8]. In a second step we developed a metatranscriptomic approach to access of the evolution of the exogenous (soil) and endogenous (plant) enzymatic arsenal potentially across the Tree of Life, involved in the degradation of carbohydrate and aromatic substances in decomposing plant matter.

The methodology used to explore microbial and enzymatic diversity using High Throughput Sequencing (Illumina system) will be described. We will then correlate colonization complexity dynamics to the progress in plant cell wall degradation. Finally, the microbial ecology of the retting process will be compared to other natural plant material (e.g. forest litter) degradation processes.

## Acknowledgements

This work is funded within the framework of the collaborative French 'Future project' SINFONI.

## References

[1] Roland, J.-C., Mosiniak, M. & Roland, D., 1995. Dynamique du positionnement de la cellulose dans les parois des fibres textiles du lin (Linum usitatissimum). Acta Botanica Gallica, 142(5), pp.463-484.

[2] Kvavadze, E. et al., 2009. 30,000-Year-Old Wild Flax Fibers. Science, 325(5946), pp.1359-1359.

[3] Md. Tahir, P. et al., 2011. Retting process of some bast plant fibres and its effect on fibre quality: A review. BioResources, 6(4), pp.5260-5281.

[4] Akin, D.E., 2013. Linen most useful: Perspectives on structure, chemistry, and enzymes for retting flax. ISRN Biotechnology, 2013, p.23.

[5] Martin, N. et al., 2013. Influence of the degree of retting of flax fibers on the tensile properties of single fibers and short fiber/polypropylene composites. Industrial Crops and Products, 49, pp.755-767.

[6] Lombard, V. et al., 2014. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic acids research, 42(Database issue), pp.D490-5.

[7] Langille, M.G. et al., 2013. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat Biotechnol, 31(9), pp.814-821.

[8] Nguyen, N.H. et al., 2016. FUNGuild: An open annotation tool for parsing fungal community datasets by ecological guild. Fungal Ecology, 20, pp.241-248.
# Splicing Lore: Speeding up the identifications of splicing factors regulating alternative exons across physiological and pathological conditions.

HELENE POLVÈCHE¹, NAHED BOUCHOUICHA¹ and DIDIER AUBOEUF¹

¹ LBMC, 46 Allée d'Italie, 69364 Lyon Cedex 07, France

Corresponding Author: helene.polveche@ens-lyon.fr

Alternative splicing (AS) leads to the production of different transcripts from each gene and relies on splicing factors (SFs) site recognition. Massive splicing variations are observed in many diseases. However, identifying which mechanisms that are responsible for these variations is challenging. To address this concern, we developed "Splicing Lore" which aims at predicting the SF(s) regulating the inclusion rate of exons from a selected list of alternative exons.

Splicing Lore contains the list of exons regulated by 64 SFs, which corresponds to 94 publicly datasets. These datasets were analyzed with an homemade pipeline (FaRLine) allowing to quantify the effect of SF downregulation on the inclusion rate of human exons. The results were stored in a `MySql` database, Splicing Lore DB that will be available for download.

An user-friendly *interface* allows to enter a list of exons and interrogate Splicing Lore DB. This allows users to know which SF(s) can regulate at least part of the exons of interest. Splicing Lore provides support toward the identification of the molecular mechanisms driving splicing variations across physiological or pathological situations.

Keyword : *alternative splicing*, *splicing factor*, *RNAseq*, *database*, *interface*
# A method for DNA virus detection and quantification during pregnancy based on noninvasive prenatal testing whole genome sequencing results

Virginia CHESNAIS<fnref n="1" />, Alban OTT<fnref n="1" />, Emmanuel CHAPLAIS<fnref n="1" />, Samuel GABILLARD<fnref n="1" />, Christelle VAULOUP-FELLOUS<fnref n="2" />, Alexandra BENACHI<fnref n="3" />, Jean-Marc COSTA<fnref n="4" />, and Eric GINOUX<fnref n="1" />

Corresponding Author: vchesnais@lifeandsoft.com

Human cytomegalovirus (HCMV) is a DNA β-herpevirus of critical importance to human health during pregnancy. HCMV primary infection of pregnant women could lead to congenital infection of fetus and could have severe clinical complication in child [1]. With 10–100 billion fragments per milliliter of plasma, circulating cell-free DNA is an information-rich window into human physiology, with rapidly expanding applications in genetic prenatal diagnosis. The whole genome sequencing (WGS) of cell-free plasma DNA is classically used to diagnose fetal aneuploidy during pregnancy. Because WGS has also become a standard tool for pathogen discovery in biological samples [2], the purpose of this study is to propose a new method to detect and quantify circulating viral DNA during pregnancy using the same sequencing results as noninvasive prenatal testing whole genome sequencing data.

Our approach used the human unmappable reads to search viral specific reads and quantify the species viral load on low depth sequencing data, by combining multiple steps of alignment on reference genome to filter host reads and mapped exogenous reads on in-house validated reference genome of different viral species. Different quality steps checked that sequencing results had a sufficient quality to be analyzed and that the alignment on targeted viral genome is homogenous. We used sequenced ranged samples to calibrate our method and determine a mean depth threshold to classify infected and non-infected samples with a sensitivity of 100% [88%-100] and a specificity of 100% [94%-100%]. Because the mean depth was linearly correlated to the theoretical viral concentration of each sample, we constructed a model able to determine the viral load of each sample.

Our method was then applied to a cohort of 538 pregnant women to validate our approach with real clinical samples. In this validation cohort, we found two positive samples for HCMV with very low viral load. The serologic status of these two samples revealed that the patients were immunized against HCMV at the beginning of their pregnancy, suggesting a possible viral reactivation or secondary infection. In the same cohort, we also found three samples positive for HBV. For two of them the serologic status was known: these women were positive for HBV infection, confirming that our pipeline can detect several targeted viral species.

In summary, this study demonstrates the whole potential benefit of WGS based monitoring for pregnant women to perform complete prenatal diagnosis based on a single test.

## References

[1] W. Britt, Manifestations of human cytomegalovirus infection: proposed mechanisms of acute and chronic disease, Curr. Top. Microbiol. Immunol. 325 (2008) 417–470.

[2] Q. Wang, P. Jia, Z. Zhao, VirusFinder: software for efficient and accurate detection of viruses and their integration sites in host genomes through next generation sequencing data, PloS One. 8 (2013) e64465. doi:10.1371/journal.pone.0064465.

<footnote n="1">¹ Life&Soft, 8b avenue Descartes, 92350, Plessis-Robinson, France</footnote>

<footnote n="2">² Hôpital Paul Brousse, 12 avenue Paul Vaillant Couturier, 94800, Villejuif, France</footnote>

<footnote n="3">³ Hôpital Antoine-Béclère, 157 rue de la Porte de Trivaux, 92140, Clamart, France</footnote>

<footnote n="4">⁴ Laboratoire CERBA, 11 rue de l'Équerre, 95310, Saint-Ouen-l'Aumône, France</footnote>
# RINSpector: a Cytoscape app that combines centrality analyses with DynaMine flexibility prediction

Guillaume Brysbaert<fnref n="1" />, Kevin Lorgouilloux<fnref n="1" />, Wim Vranken<fnref n="2" /> and Marc F Lensink<fnref n="1" />

Corresponding Author: guillaume.brysbaert@univ-lille1.fr

Deciphering the function of a protein is an intricate task that needs integrated approaches. In silico methods based on sequence and structure give good insights to solve this problem but few permit to integrate their results. An approach based on residue interaction networks analyses and in particular centrality analyses showed to give clues to the involvement of key residues in function and folding of a protein - e.g. [1, 2]. These networks are built from a protein structure where nodes are residues and edges are detected interactions between residues. Next to these methods, flexibility of a backbone and changes upon mutation provide additional information that help to select such key residues in the design of mutagenesis experiments and to unravel the function of a protein [3, 4].

We developed an app called RINSpector for the Cytoscape network analysis and visualization software [5] to analyze residue interaction networks and visualize flexibility predictions associated to a protein chain. This app performs centrality analyses based on shortest path lengths, with associated Z-scores, and queries the DynaMine server [6, 7] to retrieve a prediction of the dynamic of a protein chain. Results can be visualized on an interactive flexibility graph which permits selection of residues that can be mutated to compare the new graph with the wild type. A connection is made between the flexibility graph, the residue interaction network and the protein structure (if structureViz app/Chimera installed [8, 9]) to select a residue simultaneously in the three representations thus to directly visualize this residue in its context. Our tool can help with the rapid identification of key residues for protein function and stability.

RINSpector is available in the Cytoscape app store.

## References

[1] del Sol, A., et al., Residues crucial for maintaining short paths in network communication mediate signaling in proteins. Mol Syst Biol, 2006. 2: p. 2006 0019.

[2] Amitai, G., et al., Network analysis of protein structures identifies functional residues. J Mol Biol, 2004. 344(4): p. 1135-46.

[3] Teague, S.J., Implications of protein flexibility for drug discovery. Nat Rev Drug Discov, 2003. 2(7): p. 527-41.

[4] Golovanov, A.P., et al., Structural consequences of site-directed mutagenesis in flexible protein domains: NMR characterization of the L(55,56)S mutant of RhoGDI. Eur J Biochem, 2001. 268(8): p. 2253-60.

[5] Shannon, P., et al., Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res, 2003. 13(11): p. 2498-504.

[6] Cilia, E., et al., From protein sequence to dynamics and disorder with DynaMine. Nat Commun, 2013. 4: p. 2741.

[7] Cilia, E., et al., The DynaMine webserver: predicting protein dynamics from sequence. Nucleic Acids Res, 2014. 42(Web Server issue): p. W264-70.

[8] Morris, J.H., et al., structureViz: linking Cytoscape and UCSF Chimera. Bioinformatics, 2007. 23(17): p. 2345-7.

[9] Pettersen, E.F., et al., UCSF Chimera--a visualization system for exploratory research and analysis. J Comput Chem, 2004. 25(13): p. 1605-12.

<footnote n="1">¹ University of Lille, CNRS UMR8576 UGSF, F-59000 Lille, France</footnote>

<footnote n="2">² Interuniversity Institute of Bioinformatics in Brussels (IB)², ULB-VUB and Structural Biology Brussels, VUB, Brussels, Belgium</footnote>
# Automated generation and analysis of parametric kinetic models obtained from biochemical interaction maps

Marion BUFFARD<fnref n="1" />,<fnref n="2" />, Oscar O. ORTEGA<fnref n="3" />, Carlos F. LOPEZ<fnref n="3" /> and Ovidiu RADULESCU<fnref n="1" />

Corresponding Author: marion.buffard@etu.umontpellier.fr

The biochemical interaction map database, KEGG [1], contains approximately 500 Pathway maps. Creating parametric kinetic models manually for all of the interactions contained within all of the available maps would be tedious and would require the knowledge of tens of thousands of unknown kinetic parameters. A few approaches have been proposed to perform this task automatically such as for instance Path2Models [2]. However, these approaches do not cope with parametric problem. We propose a new method to automate this process including the generation of semi-quantitative kinetic laws. The analysis of such large dynamic models could allow for increased comprehension of the different phenotype or behavior of the cell and could elucidate new key proteins.

Using PySB [3], a framework for building mathematical models of biochemical pathways, we develop a tool that generates biochemical kinetic models from static interaction maps. The models are built as Python programs using PySB libraries for mechanistic interactions. PySB allows one to divide models into modules and to call libraries of reusable elements that encode standard biochemical actions. Species in the PySB model can be described with multiple sites and states, which reduce considerably the numbers of entities to declare.

To each interaction we associate a semi-quantitative parameter (integer) representing the order of magnitude of the interaction timescale. The values of these integers rank the interactions according to their timescales. Each interaction of the static model provides a system of rules. These rules are parametrized consistently with the interaction timescale orders. We consider that all internal processes, with exception to the limiting step, are more rapid than the timescale of the interaction. The choice of the limiting step can correspond either to quasi-equilibrium or to quasi-stationary conditions. The parametrized kinetic models are further analyzed. We are interested in classifying possible states and transitions between these states. We use tropical equilibration branches, a concept recently introduced in [4], as proxies for metastable states of the biochemical system.

As case studies we consider signaling models MAPK and TRAIL to compare our model to known models created using the classical approach.

## Acknowledgements

This participation in the 18th JOBIM has been funded with thanks to Labex Numev (http://www.lirmm.fr/numev/).

## References

[1] Hiroyuki Ogata, Susumu Goto, Kazushige Sato, et al. KEGG: Kyoto encyclopedia of genes and genomes. *Nucleic acids research*, vol. 27, no 1, p. 29-34, 1999.

[2] Finja Büchel, Nicolas Rodriguez, Nei Swainstonl, et al. Path2Models: large-scale generation of computational models from biochemical pathway maps. *BMC systems biology*, vol. 7, no 1, p. 116, 2013.

[3] Carlos F. Lopez, Jeremy L. Muhlich, John A. Bachman, et al. Programming biological models in Python using PySB. *Molecular systems biology*, vol. 9, no 1, p. 646, 2013.

[4] Satya Swarup Samal, Aurélien Naldi, Dima Grigoriev, et al. Geometric analysis of pathways dynamics: application to versatility of TGF-β receptors. *Biosystems*, vol. 149, p. 3-14, 2016.

<footnote n="1">¹ DIMNP, 163 rue Auguste Broussonnet, 34090, Montpellier, France</footnote>

<footnote n="2">² Master of Sciences and Digital Technologies for Healthcare BioInformatics, Montpellier University, 34090, Montpellier, France</footnote>

<footnote n="3">³ Center for Quantitative Sciences, Vanderbilt University, 2220 Pierce Ave, Tennessee 37232, Nashville, United States</footnote>
# Étude de processus de coalescence dans les paysages contemporains : bibliothèques template C++ pour le Calcul Bayésien Approché.

Arnaud BECHELER¹, Camille CORON² et Stéphane DUPAS¹

¹ EGCE, Bâtiment 13, Campus CNRS, 1 Avenue de la Terrasse, 91198 Gif-sur-Yvette, France

² LMO, Faculté des Sciences d'Orsay Université Paris-Sud, F-91405 Orsay Cedex, France

Correspondant : Arnaud.Becheler@egce.cnrs-gif.fr

## 1 Introduction

Les modalités d'évolution récente des populations contemporaines deviennent une question majeure à l'heure de quantifier les impacts des changements globaux. Par exemple, dans le cadre d'une espèce invasive comme le frelon à pattes jaunes (Vespa velutina) en France, l'analyse d'un jeu de données génétique peut permettre la compréhension des réactions de l'espèce à son nouvel environnement.

En couplant des modèles démographiques spatialement explicites à des modèles de coalescence, il est possible d'inférer les paramètres des lois de croissance et de dispersion de l'espèce par Calcul Bayésien Approché (Approximate Bayesian Computation, ABC) [1]. Des modèles environnementaux permettent de présenter les lois de croissances ou de dispersion comme une fonction des patrons paysagers locaux [2].

Toutefois, les méthodes de simulation antérieures reconstruisent le coalescent jusqu'à l'ancêtre commun le plus récent (MRCA) des gènes échantillonnés [3]. Or, le MRCA étant situé dans une fenêtre spatiotemporelle lointaine, cela force à renseigner des processus historiques anciens peu informés et/ou peu informatifs. La nouvelle méthode présentée ici permet d'éviter les coûts d'une prise en compte de l'histoire ancienne (coût en hypothèses, en calcul et en données) en recentrant l'analyse sur les processus de coalescence très récents, c'est à dire ceux qui permettront d'informer l'histoire invasive du frelon asiatique.

## 2 Méthode

La nécessité d'envisager de nombreux modèles de différents niveaux de complexité, conjointement au besoin de performance imposé par l'ABC, nous a mené à développer des bibliothèques template C++ pour la simulation de processus de coalescence. Ces bibliothèques génériques, modulaires et extensibles, sont destinées à être rendues libres et offrent dans la définition du modèle simulatoire une liberté et une efficacité à notre connaissance sans équivalents actuels.

Lors d'une simulation, la coalescence peut être interrompue brutalement à la date où les données deviennent trop insuffisantes pour renseigner le processus, menant ainsi à une forêt de généalogies de gènes partitionnant le jeu de données génétique. Forêt de coalescents simulée et jeu de données observé sont alors convertis en partitions floues sans perte d'information, suite à quoi la procédure ABC permet d'accepter les valeurs de paramètres ayant servi à la simulation si la distance de transfert floue [4] calculée entre la partition simulée et la partition observée est inférieure à un certain seuil.

## Références

[1] Estoup, Arnaud, et al. Combining genetic, historical and geographical data to reconstruct the dynamics of bioinvasions: application to the cane toad Bufo marinus. *Molecular ecology resources* (10.5):886-901, 2010.

[2] He, Qixin, Danielle L. Edwards, and L. Lacey Knowles. Integrative testing of how environments from the past to the present shape genetic structure across landscapes. *Evolution* (67.12): 3386-3402, 2013.

[3] Currat, Mathias, Nicolas Ray, and Laurent Excoffier. SPLATCHE: a program to simulate genetic diversity taking into account environmental heterogeneity. *Molecular Ecology Notes* (4.1): 139-142, 2004.

[4] Campello, Ricardo JGB. Generalized external indexes for comparing data partitions with overlapping categories. *Pattern Recognition Letters* (31.9): 966-975, 2010.
# PROqPCR : a Shiny web application for PROcessing of qRT-PCR data

Mathilde SAUTREUIL<fnref n="1" /> and Caroline BÉRARD<fnref n="2" />

Corresponding author: mathilde.sautreuil@etu.univ-rouen.fr

**Summary:** On this poster, we present PROqPCR, a user friendly tool enabling the PROcessing of quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) data. The qRT-PCR is a common-used technique, which allows to quantify the transcriptome for pre-defined specific genes [1]. Despite the emergence of next generation sequencing, the qRT-PCR remains widely used in the laboratories. Indeed, one use of this experiment is, for instance, the mandatory validation of the results when high-throughput experiment is performed.

When a qRT-PCR experiment is performed, many files are created by the quantitative PCR instrument, and there is a long and repetitive treatment to do to carry out the analysis (averaging over replicates, normalization). That's why we automated it in a web application. This application allows biologists to perform easily the processing of qRT-PCR data by providing only the experimental design and the files created by the software of the quantitative PCR instrument.

At the end of the analysis, the biologist can visualize and explore the results from several graphs. Three different graphs are proposed and can be downloaded in png format:

1. Barplot of all conditions for a given gene;

2. Barplot of all genes for a given condition;

3. Comparison of barplots for selected genes in all conditions.

Moreover, two graphs are also provided in order to facilitate the comparison between the results of qRT-PCR with that of a RNA-seq experiment.

PROqPCR is written in R using the Shiny library, which provides access to powerful R-based functions and libraries through a simple user interface.

PROqPCR is free, open source and is available at https://qpcrapp.shinyapps.io/proqpcr/.

## References

[1] C. A. Heid, J. Stevens, K. J. Livak, and P. M. Williams. Real time quantitative PCR. *Genome Research*, 6(10):986–994, January 1996.

<footnote n="1">¹ LMRS UMR 6085, Normandie Université, Avenue de l'Université, BP.12, F76801 Saint-Étienne-du-Rouvray, France</footnote>

<footnote n="2">² LITIS EA 4108, Normandie Université, Avenue de l'Université, F76801 Saint-Étienne-du-Rouvray, France</footnote>
Bivariate Negative Binomial Mixture Model for the analysis of RNA-seq data

Mathilde SAUTREUIL¹, Nicolas VERGNE¹, Antoine CHANNAROND¹, Angelina ROCHE³, Gaëlle CHAGNY¹
and Caroline BÉRARD²

¹ LMRS UMR 6085, Normandie Université, Avenue de l'Université, BP.12, F76801
Saint-Étienne-du-Rouvray, France

² LITIS EA 4108, Normandie Université, Avenue de l'Université, BP.12, F76801
Saint-Étienne-du-Rouvray, France

³ CEREMADE UMR CNRS 7534, Université Paris Dauphine, Place du Maréchal de Lattre de
Tassigny, 75016 Paris, France

Corresponding author: mathilde.sautreuil@etu.univ-rouen.fr

**Summary:** RNA sequencing constitutes a method of choice to quantify gene expression. This experiment provides discrete read counts assigned to target genome regions measuring the expression level. Our goal is to perform differential analysis from these data, that is to say compare the counts of a given region between two conditions. Many methods of differential analysis were developed to compare two conditions using mainly statistical tests [1,2]. Here we propose a mixture model to search the differentially expressed genes. The idea is to classify the genes into three groups : one group where expression is the same between the two conditions and two groups where expression is higher in a condition than in the other one.

We developed a Bivariate Negative Binomial Mixture Model. Let $X_t = (X_{1t}, X_{2t})$ be the counts for each gene $t$ in the conditions one and two, and $Z_t$ the group this gene belongs. As $X_t$ is bivariate, we have proposed to work with a mixture model in two dimensions. Moreover, the RNA-seq data are modeled using a Negative Binomial distribution because these data are discrete and overdispersed.

As written in Shi and Valdez [3], $\bar{X} = (X_1, X_2)$ follows a Bivariate Negative Binomial distribution if there exists the independent variables $Y_1, Y_2, Y_3$ such as $Y_i \sim NB(\theta_i, \alpha_i)$ with $i = 1, 2, 3$ and

$$ \begin{cases} X_1 = Y_1 + Y_3 \\ X_2 = Y_2 + Y_3. \end{cases} $$

Thus, $X_i \sim NB(\theta_i + \theta_3, \alpha_i + \alpha_3)$, with $\theta_i = \frac{\theta_3\alpha_i}{\alpha_3}$.

The Bivariate Negative Binomial Mixture Model is written as :

$$ \mathbb{P}(x|\psi) = \sum_{k=1}^{K} \pi_k MNB(x|(\theta_k, \alpha_k)), \text{ with} $$

- K the number of groups (here $K = 3$ as written before)

- MNB is the Multivariate Negative Binomial distribution

- $\psi = (\pi_1, ..., \pi_K, \theta_1, ..., \theta_K, \alpha_1, ..., \alpha_K)$, where $\theta_k = (\theta_{k1}, \theta_{k2}, \theta_{k3})$ et $\alpha_k = (\alpha_{k1}, \alpha_{k2}, \alpha_{k3})$

- $0 < \pi_k \le 1, k = 1, ..., K \text{ et } \sum_k \pi_k = 1$.

The parameters are estimated with the EM algorithm [4]: $\pi$ and $\theta$ have explicit estimators but $\alpha$ is estimated numerically. Then the genes are classified in one of the three groups according to the Maximum A Posteriori rule based on the posterior probabilities.

Through a simulation study, we will compare the obtained results with a Bivariate Gaussian Mixture Model, a Bivariate Poisson Mixture Model and the model we have developed. Moreover, we will compare our model with methods commonly used like DESeq2 [1]. Finally, we will apply the model on real data from an RNA-seq experiment on the *Phaeodactylum Tricornutum* diatom. In order to handle with replicates, we propose to average the replicates which is rounded at the closer integer.

References

[1] Michael I. Love, Wolfgang Huber, and Simon Anders. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biology*, 15:550, 2014.

[2] Paul L. Auer and Rebecca W Doerge. A Two-Stage Poisson Model for Testing RNA-Seq Data. *Statistical Applications in Genetics and Molecular Biology*, 10(1), 2011.

[3] Peng Shi and Emiliano A. Valdez. Multivariate Negative Binomial Models for Insurance Claim Counts. SSRN Scholarly Paper ID 2175226, Social Science Research Network, Rochester, NY, November 2012.

[4] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. *Journal of the Royal Statistical Society, Series B*, 39(1):1–38, 1977.
Quality evaluation of the mapping of small RNA-footprinting reads

Pauline Fourgoux¹, Etienne Delannoy¹, Claire Lurin¹, Guillem Rigaill¹, Véronique Brunaud¹

¹ Institute of Plant Science Paris-Saclay (IPS2), Bât. 630 Plateau du Moulon rue Noetzlin, 91192, Gif sur Yvette cedex, FRANCE

Corresponding Author: pauline.fourgoux@u-psud.fr

Pentatricopeptide repeat (PPR) proteins are modular proteins which bind RNA [1]. They are found in every eukaryotes with a high number in Plants. There are about 500 PPRs in Arabidopsis. PPRs are key regulators in mRNA processing of organelles. For example they are involved in splicing, editing and stability[1]. The goal of our project is to detect RNA binding sites of PPRs using RNA-footprinting experiments. Arabidopsis thaliana will be the plant model for this study. PPR binding sites are small sequences of 8-30bp. Thus reads obtained from RNA-footprinting experiments are small. Correctly mapping those reads on chloroplast and mitochondrial genomes is difficult [2] yet crucial for our project.

Indeed it is important to correctly map those small reads to identify real PPR binding sites. Mapping small RNAs can be complicated, because they can match at different positions on the genome. Default parameters of mapping tools do not allow multiple hits and they can not be used in this case. In the case of PPRs, there are different causes to explain multiple hits. Firstly, some regions are repeated regions. Such regions lead to artefactual multiple hits. Secondly, some PPR proteins can bind to several identical RNA [4] at different position on the genome. In that case, we get “real” multiple hits.

The goal of our study is to identify the best tool to map our small RNA reads. We used pIRS [3] to simulate data and assess the mapping quality. pIRS allows to include SNP and InDel error rate parameters. Those parameters are important, because the PPR recognition code is inaccurate [1]. The analysis of RNA-Seq simulation showed two results. Firstly, they enabled to identify complex regions prone to multiple hits. We excluded those regions from the rest of our analysis. Secondly, they allowed us to calibrate quality threshold to use in case of multiple hits.

To conclude this precise quality assessment of mapping tools for small reads will facilitate the identification of real PPR binding sites using peak detection approach mostly developed for ChIP-Seq. Finding those sites is the main goal of our project. In particular, some of these sites can be related to cytoplasmic male sterility. Finding these sites will provide easy ways of producing F1 hybrids which are of high agronomical value.

References

[1]Manna, S. (2015). An overview of pentatricopeptide repeat proteins and their applications. *Biochimie*, **113**, 93-99.

[2]Johnson, N. R., Yeoh, J. M., Coruh, C., & Axtell, M. J. (2016). Improved placement of multi-mapping small RNAs. *G3: Genes* | *Genomes* | *Genetics*, **6**(7), 2103-2111.

[3]Hu, X., Yuan, J., Shi, Y., Lu, J., Liu, B., Li, Z., ... & Yue, Z. (2012). pIRS: Profile-based Illumina pair-end reads simulator. *Bioinformatics*, **28**(11), 1533-1535.

[4]Barkan, A., Rojas, M., Fujii, S., Yap, A., Chong, Y. S., Bond, C. S., & Small, I. (2012). A combinatorial amino acid code for RNA recognition by pentatricopeptide repeat proteins. *PLoS Genet*, **8**(8), e1002910.
Medical diagnosis pipelines on the new AP-HP bioinformatics platform

Jocelyn BRAYET¹, Camille BARETTE¹, Mathieu BARTHELEMY¹, Vivien DESHAIES¹ and Alban LERMINE¹

¹ Bioinformatic platform of AP-HP, 5 rue Santerre (Hôpital Rothschild), 75012, Paris, France

Corresponding Author: jocelyn.brayet@aphp.fr

The Assistance Publique – Hôpitaux de Paris (AP-HP) is a teaching hospital groupment with a European dimension globally recognized. The AP-HP is organized into twelve hospital groups, for a total of 39 hospitals localized in Paris and its region. Currently, those hospitals attend each year 8 million patients.

One year ago, a new bioinformatics platform was created for multiple missions: the progressive centralization of the storage of genomic data produced by hospitals, their analyses in controlled and standardized workflows and the provisioning of tools for results exploitation.

In this abstract, we present our solution to develop and standardize medical diagnosis pipelines. The goal of this project is to propose a development environment for tools and parameters testing and an analysis environment with fixed pipelines and for results reproducibility. We decided to use *Docker* [1], *Galaxy* [2], and *Snakemake* [3] technologies. Docker allows to eliminate tool dependencies problems and sets a version tool in an image. We created a Docker image per tool per version. Galaxy is a web-based platform, designed to provide easy access to a versatile toolbox for biological users. Using Galaxy interface, scientists can test tools, tool versions and parameters but also design new workflows. When the pipeline is validated, Snakemake rules are written, versioned and tagged. Snakemake is a workflow management system with implicit rule implementation (input and output logic). The advantage over Make is the capability to allow Python to be interspersed through the pipeline and thereby reducing a lot of ambiguity [4]. As other pipeline frameworks, error recovery, automatic parallelization and workflow integrity features are included in Snakemake. With this project, on one hand, those Docker images are executed by a Galaxy instance. On the other hand, those same Docker images are run by Snakemake pipelines. To ease medical diagnosis routine, AP-HP's scientists are able to execute tagged workflows with their data through a web interface.

Furthermore, every Docker tools and Snakemake pipelines are tested with GitLab-CI [5]. GitLab is a web-based repository manager. GitLab-CI provides continuous integration features, so that when a new tool version is added in a Docker image, GitLab will build, run the container and automatically test the tool.

Key words: Galaxy, Docker, Snakemake, diagnosis pipelines and continuous integration

References

[1] Web site: https://www.docker.com/

[2] Enis A, Dannon B, Marius VDB, Daniel B, Dave B, Martin C, John C, Dave C, Nate C, Carl E, Björn G, Aysam G, Jennifer HJ, Greg VK, Eric R, Nicola S, Nitesh T, James T, Anton N, and Jeremy G. The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2016 update. *Nucleic Acids Research* (2016) 44(W1): W3-W10 doi:10.1093/nar/gkw343

[3] Johannes Köster, Sven Rahmann; Snakemake—a scalable bioinformatics workflow engine. *Bioinformatics* 2012; 28 (19): 2520-2522. doi: 10.1093/bioinformatics/bts480

[4] Jeremy Leipzig; A review of bioinformatic pipeline frameworks. *Brief Bioinform* 2017; 18 (3): 530-536. doi: 10.1093/bib/bbw020

[5] Web site: https://about.gitlab.com/gitlab-ci/
# Mathematical modeling of a genetic network controlling the regulation of Fe-S biogenesis

Firas HAMMAMI<fnref n="1" />,<fnref n="2" />, Frédéric BARRAS<fnref n="1" />, Pierre MANDIN<fnref n="1" /> and Elisabeth REMY<fnref n="2" />

Corresponding Authors: elisabeth.remy@univ-amu.fr or pmandin@imm.cnrs.fr

Iron-sulfur (Fe-S) clusters are cofactors conserved in all domains of life. A-Type Carriers (ATC) proteins play a crucial role in Fe-S biogenesis by delivering clusters to their targets [1]. In *E. coli*, the ErpA ATC protein is essential in aerobic growth, whereas NfuA plays an important role in stress conditions, such as oxidative stress or iron deficiency.

Recent work in our laboratory showed that *erpA* expression is regulated by IscR, the main Fe-S cluster homeostasis regulator, and by the non-coding RNA RyhB, expressed in iron deficient conditions [2]. While both regulators repress *erpA* and *nfuA* expression in opposite conditions in regards to iron concentration, these genes present different expression profiles.

In order to understand the mechanisms underlying *erpA* and *nfuA* regulation, we used a mathematical modeling approach relying on logical formalism. This model takes into account environmental conditions that perturb Fe-S biogenesis such as iron and oxygen levels. The model was validated with experimental data, allowing us to make predictions in different mutants and to infer qualitative network properties. For instance, the model suggests different inhibition thresholds for RyhB.

(Work in progress)

## References

[1] Roche B, Aussel L, Ezraty B, Mandin P, Py B, Barras F. Iron/sulfur proteins biogenesis in prokaryotes: formation, regulation and diversity. *Biochim Biophys Acta*, Mar;1827(3):455-69, 2013

[2] Mandin P, Chareyre S and Barras F. A Regulatory Circuit Composed of a Transcription Factor, IscR, and a Regulatory RNA, RyhB, Controls Fe-S Cluster Delivery. *Mbio*, Sep 20;7(5), 2016

<footnote n="1">¹ Laboratoire de Chimie Bactérienne, Aix Marseille Université-CNRS, UMR 7283, Institut de Microbiologie de la Méditerranée, CNRS, 31 Chemin Joseph Aiguier, 13009 Marseille, France</footnote>

<footnote n="2">² Institut de Mathématiques de Marseille, Aix-Marseille Université, CNRS, Centrale Marseille, I2M UMR 7373, Marseille, France.</footnote>
# First gene-annotation enrichment analysis based on bacterial core-genome variants: Insights into mammalian and avian host adaptation of *Salmonella* serovars

Meryl VILA NOVA¹*, Kevin DURIMEL¹*, Arnaud FELTEN¹, Laurent GUILLER¹, Michel-Yves MISTOU¹ and Nicolas RADOMSKI¹

¹ Laboratory for food safety (Anses), 14 rue Pierre et Marie Curie, 94700, Maisons-Alfort, France
*Poster presenters

Corresponding Author: nicolas.radomski@anses.fr

## 1 Introduction

Most of bacterial genomic studies exploring the host adaptation focus on the accessory genome describing how gain and loss of genes explain evolution processes leading colonization of new habitats. In this context, we propose a robust phylogenetic inference based on single nucleotide polymorphisms (SNPs) and recombination events, identification of fixed SNPs and small insertions/deletions (InDels) distinguishing homoplastic and non-homoplastic core-genome variants, and gene-annotation enrichment analyses in order to describe the main metabolism pathways impacted by these fixed variants during adaptation of *Salmonella* enterica subsp. Enterica to multi-host (S. Enteritidis), mammalian (S. Dublin), and avian (S. Pullorum and S. Gallinarum) hosts [1].

## 2 Results

The developed workflow ‘VARCall-GO’ produced a robust phylogenetic inference based on SNPs. The monophyletic clade S. Dublin diverged to the polyphyletic clade S. Enteritidis which includes the divergent clades S. Pullorum and S. Gallinarum. Firstly, this workflow gave the opportunity to detect intragenic and non-homoplastic fixed variants supporting the phylogenetic reconstruction [2]. Secondly, it identified representative metabolic pathways related to the host adaptation using the first gene-annotation enrichment analysis [3] based on bacterial core-genome variants. The host adaptation of *Salmonella* serovars were driven by fixed variants impacting mainly biosynthetic and metabolic processes of carbon sources alternative to glycose, amino acids, and ion transport, especially potassium.

## 3 Conclusion

We propose a new core-genome approach coupling identification of fixed SNPs and InDels with regards to the inferred phylogenetic clades, and gene-annotation enrichment analysis in order to describe the adaptation of *Salmonella* serovars Dublin, Enteritidis, Pullorum, and Gallinarum. All these generic bioinformatic tools can be applied on other bacterial genera without additional developments.

## Acknowledgements

This work was supported by the COMPARE European project aiming to provide genomic tools in the context of outbreak detection among humans and animals worldwide. *Horizon 2020 research and innovation programme* under grant agreement No. 643476

## References

[1] Langridge GC, Fookes M, Connor TR, Feltwell T, Fasey N, Parsons BN, et al. Patterns of genome evolution that have accompanied host adaptation in *Salmonella*. Proc. Natl. Acad. Sci. 2015;112:863–8.

[2] Aberer AJ, Pattengale ND, Stamatakis A. Parallelized phylogenetic post-analysis on multi-core architectures. J. Comput. Sci. 2010;1:107–14.

[3] Rachael P Huntley, Tony Sawford, Maria J Martin, Claire O’Donovan; Understanding how and why the Gene Ontology and its annotations evolve: the GO within UniProt. *Gigascience* 2014; 3 (1): 1-9. doi: 10.1186/2047-217X-3-4
# A phylogenomic network to decipher bacterial adaptation through horizontal gene transfer

Damien RICHARD¹,²,³, Virginie RAVIGNE¹, Aude CHABIRAND², Olivier PRUVOST¹ and Pierre LEFEUVRE¹

¹ UMR PVBMT, Cirad, F-97410, St Pierre, Réunion, France

² Plant Health Laboratory, ANSES, F-97410, St Pierre, Réunion, France

³ Université de la Réunion, F-97410, St Pierre, Réunion, France

Corresponding Author: damien.richard@cirad.fr

Horizontal gene transfer (HGT) is often reported as being the motor of bacterial adaptation. Among the various processes involved in HGT, plasmid transfers were shown to be of prime importance for many bacterial families [1,2]. Plasmid-encoded accessory genes enrich the species' gene pool and provide new adaptive traits in response to environmental modifications, such as the use of new antimicrobials or the colonization of a new ecological niche [3,4]. Therefore, in order to further our understanding of bacterial adaptation and ecology, it is pivotal to improve our knowledge on (1) plasmid genomic structure, (2) mechanisms involved in plasmid transfers and (3) their limits. Phylogenomic networks are useful tools to picture HGT among prokaryotes [1]. Indeed, each edge of the network represent an HGT event between two genomes (i.e. nodes of the network). Using this approach, it was demonstrated that HGT is limited by the taxonomic relatedness of the organisms, the characteristics of the host genomes (such as GC content or codon usage) and the ecological niche they occupy.

Copper-based antimicrobial compounds are widely used to control plant bacterial pathogens. *Xanthomonas citri* pv. *citri* (Xcc), a major citrus pathogen, has adapted in response to this selective pressure [5,6] and it was demonstrated that its copper-resistance was plasmid-borne. In order to understand the genetic basis of this adaptation, we fully sequenced 13 copper-resistant strains from the Xanthomonadaceae family: six Xcc strains [7], six *Xanthomonas* strains pathogenic to solanaceous species [8] and one commensal *Stenotrophomonas* strain [9]. Genome comparison between the six Xcc strains revealed that the copper-resistance genes were encoded on an adaptive transposon located on a conjugative plasmid. Most strains of other *Xanthomonas* species encoded a highly conserved copy of the entire plasmid, but one that displayed a similar adaptive transposon inserted in a distinct plasmid. The commensal *Stenotrophomonas* strain presented with a chromosomally encoded copy of the transposon. These findings highlighted the existence of two overlapping levels of mobility (transposon and plasmid) [10].

To further investigate the spread of copper-resistance genes among the Xanthomonadaceae family, we searched for homologues of every plasmid's genes in the public NCBI databases NR and WGS. Based on sequence identity, we constructed and analysed phylogenomic networks. It revealed that homologues of the plasmid present in Xcc were only identified from Xcc and few *Xanthomonas* species pathogenic to solanous species. In contrast, genes homologous to the adaptive transposon were detected from 14 species included in five genera (*Xanthomonas*, *Stenotrophomonas*, *Pseudoxanthomonas*, *Pelomonas* and *Pseudomonas*) and grouping in three families. Globally, homologues of genes encoded on the transposon were found further apart in the taxonomy than plasmid backbone homologues.

A similar approach using all available plasmids from the Xanthomonadaceae family was applied to more finely define the reach of gene exchange for bacteria in this family. The analyses of the resulting networks would certainly help uncover the global context of gene exchange, their limits and specificities. It may also reveal specific trends for genes involved in adaptation within the particular context of agricultural settings.

## References

[1] O. Popa and T. Dagan. Trends and barriers to lateral gene transfer in prokaryotes. *Curr Opin Microbiol*, (14):615-623, 2011.

[2] S. Halary, J. W. Leigh, B. Cheaib, P. Lopez and E. Bapteste. Network analyses structure genetic diversity in independent genetic worlds. *Proceedings of the National Academy of Sciences, USA*, (107):127-132, 2010.

[3] J. L. Hobman and L. C. Crossman. Bacterial antimicrobial metal ion resistance. *J Med Microbiol*, (64):471-497, 2015.

[4] H. Ochman, J. G. Lawrence and E. A. Groisman. Lateral gene transfer and the nature of bacterial innovation. *Nature*, (405):299-304, 2000.
[5] D. Richard, C. Boyer, S. Javegny, K. Boyer, P. Grygiel, O. Pruvost, A. L. Rioualec, V. Rakotobe, J. Iotti, R. Picard, C. Vernière, C. Audusseau, C. François, V. Olivier, A. Moreau and A. Chabirand. First Report of Xanthomonas citri pv. citri Pathotype A Causing Asiatic Citrus Canker in Martinique, France. *Plant Disease*, (100):1946-1946, 2016.

[6] D. Richard, N. Tribot, C. Boyer, M. Terville, K. Boyer, S. Javegny, M. Roux-Cuvelier, O. Pruvost, A. Moreau, A. Chabirand and C. Vernière. First Report of Copper-resistant Xanthomonas citri pv. citri Pathotype A Causing Asiatic Citrus Canker in Réunion, France. *Plant Disease*, (101):503, 2016.

[7] D. Richard, C. Boyer, C. Vernière, B. I. Canteros, P. Lefeuvre and O. Pruvost. Complete Genome Sequences of Six Copper-Resistant Xanthomonas citri pv. citri Strains Causing Asiatic Citrus Canker, Obtained Using Long-Read Technology. *Genome Announcements*, (5):2017.

[8] D. Richard, C. Boyer, P. Lefeuvre, B. I. Canteros, S. Beni-Madhu, P. Portier and O. Pruvost. Complete Genome Sequences of Six Copper-Resistant Xanthomonas Strains Causing Bacterial Spot of Solaneous Plants, Belonging to X. gardneri, X. euvesicatoria, and X. vesicatoria, Using Long-Read Technology. *Genome Announcements*, (5):2017.

[9] D. Richard, C. Boyer, P. Lefeuvre and O. Pruvost. Complete Genome Sequence of a Copper-Resistant Bacterium from the Citrus Phyllosphere, Stenotrophomonas sp. Strain LM091, Obtained Using Long-Read Technology. *Genome Announcements*, (4):2016.

[10] D. Richard, V. Ravigné, A. Rieux, B. Facon, C. Boyer, K. Boyer, P. Grygiel, S. Javegny, M. Terville, B. I. Canteros, I. Robène, C. Vernière, A. Chabirand, O. Pruvost and P. Lefeuvre. Adaptation of genetically monomorphic bacteria: evolution of copper resistance through multiple horizontal gene transfers of complex and versatile mobile genetic elements. *Molecular Ecology*, Special Issue: Microbial Local Adaptation, 2017.
# Méthodes et outils de construction de super-arbres en phylogénie

Morgan SOULIÉ², Vincent LEFORT¹ et Anne-Muriel Arigon CHIFOLLEAU¹

¹ LIRMM, UMR 5506, CNRS et Université Montpellier, 161 rue Ada, 34095, Montpellier, France

² Parcours Bioinformatique, Connaissances, Données (BCD) du Master Sciences et Numérique pour la Santé (SNS)

Auteur référent: morgan.soulie@etu.umontpellier.fr

## Résumé

Un super-arbre est un arbre phylogénétique construit à partir d'une collection d'arbres, dits arbres sources, partageant complètement ou partiellement un même ensemble d'espèces. Les informations portées par chaque arbre source peuvent être complémentaires et contradictoires entre elles. Un super-arbre représente une synthèse d'hypothèses relationnelles qui couvre l'ensemble des organismes présents dans chacun des arbres sources. Il permet de révéler les liens de parenté non visibles dans les arbres sources pris séparément. Afin de résoudre les problèmes de compatibilité entre les parties communes des arbres sources, de nombreuses méthodes de construction de super-arbre ont été développées. La plateforme ATGC (www.atgc-montpellier.fr) de l'IFB (Institut Français de Bioinformatique), adossée à l'équipe MAB (Méthodes et Algorithmes pour la Bioinformatique) du LIRMM, est dédiée à la bioinformatique pour la génomique évolutive comparative et fonctionnelle. Nous souhaitons développer un portail web convivial dédié aux analyses phylogénomiques intégrant les outils de référence du domaine, et notamment les outils de construction de super-arbres. Ce portail permettra de guider les biologistes dans leurs analyses afin de leur fournir des services adaptés à leurs problématiques. Pour ce faire, nous avons réalisé un état des lieux des méthodes et outils de construction de super-arbre existants et identifié les plus utilisés. Nous avons sélectionné ceux qui pourraient être directement intégrés à la plateforme.

Plusieurs approches sont utilisées dans la construction de super-arbres : l'approche libérale ou de "vote" qui résout les conflits en "faisant voter" les arbres sources et en optant pour l'alternative topologique qui maximise un critère d'optimisation, celui-ci variant d'une méthode à l'autre (MRP - représentation matricielle avec parcimonie, similarité, bipartitions, quadruplets, SPR - subtree pruning and regrafting, ...); l'approche consensus qui utilise des méthodes issues de la théorie des graphes, l'approche "veto" dont la topologie des super-arbres inférés respectent la topologie de chacun des arbres sources; ou de nouvelles approches basées, par exemple, sur le maximum de vraisemblance ou sur l'analyse des relations entre paires de taxons. Les outils actuels se basent sur ces approches et implémentent certaines de ces méthodes. Notre étude nous a permis de préselectionner un ensemble d'outils implémentant les principales méthodes : PhySIC IST (Scornavacca et al., 2008 ; approche libérale et approche "veto"), L.U.St (Akanni et al., 2014 ; approche par maximum de vraisemblance), COSPEDTree (Bhattacharyya et al., 2016 ; approche par paire de taxons), SPRsupertrees (Whidden et al., 2014 ; approche libérale) et Clann (Creevey et al., 2005 ; approche libérale et approche consensus). Parmi les méthodes de super-arbre, la méthode MRP est la méthode la plus populaire (Olaf R.P. Bininda-Emonds, 2004), et les articles la décrivant (Baum, 1992 et Ragan, 1992) ont respectivement été cités 459 et 529 fois. Par la suite, nous porterons donc une attention particulière à l'implémentation et l'utilisation de cette méthode.

## Remerciements

Cette participation à la 18ème dition des Journes Ouvertes en Biologie, Informatique et Mathmatiques (JOBIM) a t' finance par le Labex Numev (http://www.lirmm.fr/numev/).
# Implémentation d'une interface pratique pour l'évaluation de stratégies d'exploration de la diversité moléculaire par protéogénomique

YANNICK COGNE¹, CHRISTINE ALMUNIA¹, OLIVIER PIBLE¹, DUARTE GOUVEIA², ARNAUD CHAUMOT², OLIVIER GEFFARD² and JEAN ARMENGAUD¹

¹ CEA-Marcoule LI2D, PRAE Marcel Boiteux, F-30200, Bagnols-sur-Cèze, France

² IRSTEA - Centre de Lyon, 5 rue de la Doua, 69616, Villeurbanne, France

Corresponding Author: jean.armengaud@cea.fr

Le projet ANR Proteogam propose de définir des marqueurs moléculaires d'un organisme sentinelle afin de tester rapidement son état de santé pour sonder la qualité biologique de l'eau d'une rivière. L'organisme sentinelle choisi est le gammare, un petit amphipode d'eau douce, présent dans tous les cours d'eau d'Europe. Les marqueurs moléculaires sont des protéines dont l'abondance peut varier en fonction de la présence de polluants dans les eaux de rivière. Typiquement, les concentrations de protéines présentes dans les gonades peuvent signer la présence de perturbateurs endocriniens.

Le Gammare est une espèce dite « non-modèle » car il n'existe, pour cette espèce, pas de séquence génomique de référence permettant l'attribution des spectres MS/MS enregistrés en protéomique. Les génomes de référence les plus proches sont ceux de *Hyalella azteca* et *Daphnia pulex* qui restent cependant trop éloignés pour une attribution suffisante. Pour pallier ce manque d'information, une base de données de référence peut être constituée par assemblage de novo de lectures de transcrits matures. Pour sonder la diversité des populations de *Gammarus fossarum*, la stratégie protéogénomique mise en œuvre nécessite la combinaison des outils bioinformatiques d'assemblage de données RNAseq et d'attribution protéomique [1]. Cette stratégie requiert un environnement de travail souple et reproductible, permettant la gestion simple et efficace de données volumineuses et de formats multiple.

L'environnement de travail choisi est Galaxy avec une solution portable aisée à l'aide de Docker, un programme de gestion de containers permettant son exécution sous Windows en mimant UNIX. L'interface proposée utilise donc Galaxy docker développé par Björn Grüning [2]. Une utilisation de docker dans le Galaxy docker a été choisie afin d'utiliser un container par outil. Cette séparation par outils permet via l'utilisation d'un gestionnaire d'installation (conda) d'apporter une gestion des versions et une mise à jour simple de chaque outil indépendamment ainsi qu'une forte reproductibilité des variables d'environnement au lancement de chaque outil.

La présentation détaillera l'implémentation pour l'évaluation de stratégies d'exploration de la diversité moléculaire par protéogénomique, notamment la réalisation de la gestion des données issues de nouvelles technologies de séquençage ARN et leurs assemblages.

## References

[1] Armengaud, J., Trapp, J., Pible, O., Geffard, O., Chaumot, A., & Hartmann, E. M. (2014). Non-model organisms, a species endangered by proteogenomics. *Journal of proteomics*, *105*, 5-18

[2] Björn Grüning: https://github.com/bgruening/docker-galaxy-stable

[3] Afgan, E., Baker, D., Van den Beek, M., Blankenberg, D., Bouvier, D., Čech, M., ... & Grüning, B. (2016). The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2016 update. *Nucleic acids research*, *gkw343*.
# Pea genetic map enrichment with Genotyping By Sequencing markers

Ayité Kougbeadjo¹, Grégoire Aubert¹, Mathieu Falque², David Edwards³, Philipp Bayer³,
Jacqueline Batley³, Manuel Zander³, Satomi Hayashi³, Jonathan Kreplak¹ and Judith Burstin¹

¹ INRA UMR Agroecologie, 17 Rue Sully, 21000, Dijon, France

² INRA Ferme du Moulon, 91190, Gif sur Yvette, France

³ School of biological Sciences, University of Western Australia

Corresponding Author: jonathan.kreplak@inra.fr

Pea (Pisum sativum L.) is a model plant and an important protein crop used in human and animal feed.

The laboratory developed a 13.2K SNP Infinium genotyping array and built individual and consensus genetic maps for [1] which gave information on the genome structure. We developed a pipeline on GBS (Genotyping by sequencing) data to enrich the consensus map with more markers to get a more accurate map.

New sequencing technologies as GBS (Genotyping by Sequencing) allow to access genetic information like markers cost-effectively. GBS is a simple and reproducible genome reduction approach using enzymes of restriction sensitive to methylation before sequencing. It allows reducing the number of sequence reads and can allow to avoid repetitive regions[2]. A Genotype by Sequencing (GBS) experiment was conducted by Warkentin et al. on a recombinant inbred line *Cameor x Melrose* produced by INRA.

Our pipeline first mapped the reads on the temporary reference assembly using bwa mem[3] and kept unique mapping with MAPQ (PHRED MAPping Quality) higher or equal to 30. Optical duplicates were removed with the picard tools[4]. We then called the variants using SAMtools mpileup and BCFtools[5]. The results showed that abnormal numbers of heterozygous were detected for all sites where the depth of sequencing was less than 3 reads.

To support low depth sites, we first called SNPs from the 7x resequencing data of the parents. We filtered them to only keep SNPs that are non-variant between *Cameor* resequencing reads and the *Cameor* genome assembly, and variant between *Melrose* resequencing reads and the *Cameor* genome assembly. This defines the set of possible SNPs in the progeny. Only GBS polymorphic sites that are common with parents polymorphic sites are taken into account. All GBS sites with a read depth below 3 was then called again using an in-house Python script.

We finally obtained a 473 587 SNPs of which 468 448 SNPs were placed on the genetic map. This map combined with a bionano map and a whole genome profiling map will allow us to rearrange scaffolds in our assembly and to do further analysis. We are working to encapsulate the pipeline in a reusable package.

## References

1. Tayeh N, Aluome C, Falque M, Jacquin F, Klein A, Chauveau A, et al. Development of two major resources for pea genomics: the GenoPea 13.2K SNP Array and a high-density, high-resolution consensus genetic map. Plant J.,(84): 1257-1273, 2015

2. Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, et al. A Robust, Simple Genotyping-by-Sequencing (GBS) Approach for High Diversity Species. Public Library of Science, 2011

3. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 2009

4. Picard [Internet]. Available: http://picard.sourceforge.net/

5. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and
# Développement d'une mesure du biais d'usage des codons :
application aux virus humains

Jérôme BOURRET¹,², Samuel ALIZON¹ et Ignacio G BRAVO¹

¹ MIVEGEC (UM CNRS 5290, IRD 224, UM), 911 avenue Agropolis, 34394, Montpellier, France
2 Parcours Bioinformatique, Connaissances, Données (BCD) du Master Sciences et Numérique pour la Santé (SNS), place Eugène Bataillon, 34090, Montpellier, France

Adresse mail : bourret.jerome@gmail.fr

Au cours de la synthèse d'une protéine à partir d'un ARN messager, la machinerie ribosomale traduit 61 codons en 20 acides aminés. Le code génétique est de ce fait « dégénéré » : seuls 2 acides aminés sont codés par des codons uniques alors que les 18 autres sont codés par des combinaisons de 2, 3, 4 ou 6 codons. Même si ces derniers sont généralement appelés « codons synonymes », leur utilisation ne se fait pas de manière aléatoire. Les préférences d'usage de certains codons varient entre *organismes*, entre *gènes* à l'intérieur d'un génome et même entre *les positions d'un même gène*. On parle de *biais mutationnel* : lors de la synthèse ou de la réparation de l'ADN, tous les nucléotides ne sont pas utilisés de manière identique, ce qui peut influencer la composition nucléotidique et donc biaiser l'usage du code [2]. Le second moteur du biais d'usage des codons est la * sélection traductionnelle* : elle repose sur l'hypothèse que certains codons confèrent une traduction plus efficace lors de la synthèse protéique [2].

Pour estimer le degré d'usage non-aléatoire des codons synonymes, plusieurs mesures ont été développées depuis les années 1980. Ces mesures appartiennent principalement à deux classes : l'une où les fréquences des codons synonymes d'une séquence requête sont comparées à celles d'un jeu de données de référence [3] et l'autre où l'hypothèse nulle correspond à un usage aléatoire des différents codons synonymes [4].

Notre objectif est de développer une nouvelle mesure qui combine ces deux approches : elle calcule le biais d'usage des codons d'une séquence requête face aux préférences d'usage des codons dans un jeu de données de référence tout en normalisant selon une hypothèse d'usage aléatoire des codons. Nous avons effectué le développement de cette mesure et l'avons implémentée au sein d'un outil codé en Python, R et C++. Nous avons ensuite réalisé une étude comparative entre notre mesure et les autres outils et mesures de biais d'usage des codons.

Dans un second temps, notre outil sera utilisé pour mesurer le biais d'usage des codons d'un large panel de virus humains en utilisant comme référence celui de leur hôte humain. D'éventuelles différences de biais d'usage des codons détectées par notre outil permettront d'aborder des phénomènes tels que le contrôle spatio-temporel différentiel de l'expression génique à l'intérieur d'un même génome ou l'exposition différentielle au système immunitaire des protéines virales en fonction des interactions hôte-parasite [5].

## Remerciements

Cette participation à la 18ème édition des **JOBIM** a été financée grâce au laboratoire d'excellence Numev (http://www.lirmm.fr/numev/).

## Références

[1] R. Grantham, C. Gautier, M. Gouy, R. Mercier and A. Pavé. *Codon catalog usage and the genome hypothesis*. Nucleic Acids Research 8(1):r49-r62, 1980.

[2] A. Roth, M. Anisimova and G. M. Cannarozzi. *Measuring codon usage bias*. In: Cannarozzi G.M, Schneider A, editors. Codon Evolution Mechanism and Models:University of Bern.University of Utrecht:189-198, 2012.

[3] P. M. Sharp and W. H. Li. *The codon Adaptation Index-a measure of directional synonymous codon usage bias, and its potential applications*. Nucleic Acids Research 15(3):1281-1295, 1987.

[4] F. Wright. *The 'effective number of codons' used in a gene*. Gene 87(1):23-29, 1990.

[5] M. Félez-Sánchez, J. Trösemeier, S. Bedhomme, M.I. González-Bravo, C. Kamp, and I.G. Bravo. *Cancer, Warts, or Asymptomatic Infections: Clinical Presentation Matches Codon Usage Preferences in Human Papillomaviruses*. Genome Biology and Evolution, vol. 7, pp. 2117-2135, 2015.
# Was the Chlamydial Adaptative Strategy to Tryptophan Starvation an Early Determinant of Plastid Endosymbiosis?

Ugo CENCI<fnref n="1" />, Mathieu DUCATEZ<fnref n="1" />, Derifa KADOUCHE<fnref n="1" />, Christophe COLLEONI<fnref n="1" /> and Steven G. BALL<fnref n="1" />

Corresponding Author: steven.ball@univ-lille1.fr

## Abstract

Chlamydiales were recently proposed to have sheltered the future cyanobacterial ancestor of plastids in a common inclusion [1,2]. The intracellular pathogens are thought to have donated those critical transporters that triggered the efflux of photosynthetic carbon and the consequent onset of symbiosis [3,4]. Chlamydiales are also suspected to have encoded glycogen metabolism TTS (Type Three Secretion) effectors responsible for photosynthetic carbon assimilation in the eukaryotic cytosol [1]. We now turn our attention to the reasons underlying other chlamydial lateral gene transfers evidenced in the descendants of plastid endosymbiosis. In particular, we show that half of the genes encoding enzymes of tryptophan synthesis in Archaeplastida are of chlamydial origin. Tryptophan concentration is an essential cue triggering two alternative modes of replication in Chlamydiales. In addition, sophisticated tryptophan starvation mechanisms are known to act as antibacterial defenses in animal hosts [5,6]. We propose that Chlamydiales have donated their tryptophan operon to the emerging plastid to ensure increased synthesis of tryptophan by the plastid ancestor. This would have allowed massive expression of the tryptophan rich chlamydial transporters responsible for symbiosis. It would also have allowed possible export of this valuable amino-acid in the inclusion of the tryptophan hungry pathogens. Free-living single cell cyanobacteria are devoid of proteins able to transport this amino-acid. We therefore investigated the phylogeny of the Tyr/Trp transporters homologous to *E. coli* TyrP/Mre and found yet another LGT from Chlamydiales to Archaeplastida thereby considerably strengthening our proposal.

## References

[1] Ball, S.G., Subtil, A., Bhattacharya, D., Moustafa, A., Weber, A.P.M., Gehre, L., Colleoni, C., Arias, M.-C., Cenci, U., and Dauvillée, D. Metabolic Effectors Secreted by Bacterial Pathogens: Essential Facilitators of Plastid Endosymbiosis? *Plant Cell* 25, 7–21, 2013.

[2] Facchinelli, F., Colleoni, C., Ball, S.G., and Weber, A.P.M. Chlamydia, cyanobiont, or host: who was on top in the ménage à trois? *Trends Plant Sci.* 18, 673–679, 2013.

[3] Facchinelli, F., Pribil, M., Oster, U., Ebert, N.J., Bhattacharya, D., Leister, D., and Weber, A.P.M. Proteomic analysis of the Cyanophora paradoxa muroplast provides clues on early events in plastid endosymbiosis. *Planta* 237, 637–651, 2013.

[4] Karkar, S., Facchinelli, F., Price, D.C., Weber, A.P.M., and Bhattacharya, D. Metabolic connectivity as a driver of host and endosymbiont integration. *Proc. Natl. Acad. Sci. 112*, 10208–10215, 2015.

[5] Bonner, C.A., Byrne, G.I., and Jensen, R.A. Chlamydia exploit the mammalian tryptophan-depletion defense strategy as a counter-defensive cue to trigger a survival state of persistence. *Front. Cell. Infect. Microbiol.* 4, 2014.

[6] Lo, C.-C., Xie, G., Bonner, C.A., and Jensen, R.A. The Alternative Translational Profile That Underlies the Immune-Evasive State of Persistence in Chlamydiaceae Exploits Differential Tryptophan Contents of the Protein Repertoire. *Microbiol. Mol. Biol. Rev.* 76, 405–443, 2012.

<footnote n="1">¹ Unité de Glycobiologie Structurale et Fonctionnelle, UMR8576 Centre National de la Recherche Scientifique, Université des Sciences et Technologies de Lille, Villeneuve d'Ascq, France</footnote>

<footnote n="2">² Laboratory, Address, zip code, Town, Country</footnote>
# OrthoInspector 3.0: orthology en route to big data

Yannis NEVERS¹, Arnaud KRESS¹, Raymond RIPP¹, Olivier POCH¹ and Odile LECOMPTE¹

¹ CSTB team, UMR7357 - ICube, 4 rue Kirschleger 67085 Strasbourg, France

Corresponding Author: yannis.nevers@etu.unistra.fr

Homologs, genes that derive from a common ancestor, can be separated in two classes: paralogs -deriving from an ancestral gene by a duplication event, and orthologs - deriving from an ancestor by a speciation event. Orthologs generally perform a similar function in different species with possible species-specific adaptations. Prediction of correct orthology relationship is critical in many functional and evolutionary applications: phylogenetic tree inferences, gene function predictions, genome annotation, selection of a relevant model organism for experimental studies... OrthoInspector [1] is one of the leading algorithm [2] for pairwise orthologous relationship predictions. In addition to the independent software package, the OrthoInspector website (www.lbgi.fr/orthoinspector) currently offers precomputed databases, with a sampling of 259 eukaryotes, 1568 bacteria and 120 archaea.

Here, we present our current developments to build a new release of OrthoInspector databases. This major update aims at presenting a comprehensive set of species. As a starting point, we used a selection of non-redundant proteome representing a wide variety of taxa: the Uniprot Reference Proteomes (Nov 2012) [3]. Statistical analyses of proteome contents were used to identify and filter low-quality proteomes. The 4752 selected proteomes (representing 87% of the original set) were used as our dataset for the BLAST all-versus-all [4] searches, first step of OrthoInspector calculation. With more than 20 million BLASTP searches, this step is computationally expensive. Thus, we segmented the BLAST database to parallelize this step and to perform it on the EGI computing grid [5]. Computing of orthology relationships were in turn parallelized on our local infrastructure. The final computed databases of orthology relationships provide a major breakthrough in terms of covered species with 711 eukaryotes, 3862 bacteria and 179 archaea. It makes OthoInspector database the most comprehensive orthology database available to date.

In parallel, we are developing new features to the OrthoInspector suite to improve the possibilities offered by OrthoInspector package and website. This includes support for SQLite databases, automatic update procedures to keep our resources up-to-date with Uniprot Reference Proteomes and more importantly, a definition of ortholog families. Families will be available in complement to the pairwise relationships currently supported. We plan to provide evolutionary characterization of these gene family through multiple sequence alignment, domain conservations and integration of cross-references toward external biological resources.

## References

[1] B. Linard et al., "OrthoInspector 2.0: Software and database updates," *Bioinforma. Oxf. Engl.*, vol. 31, no. 3, pp. 447-448, Feb. 2015.

[2] A. M. Altenhoff et al., "Standardized benchmarking in the quest for orthologs," *Nat. Methods*, vol. 13, no. 5, pp. 425-430, May 2016.

[3] UniProt Consortium, "UniProt: a hub for protein information," *Nucleic Acids Res.*, vol. 43, no. Database issue, pp. D204-212, Jan. 2015.

[4] S. F. Altschul, W. Gish, W. Miller, E. W. Myers, and D. J. Lipman, "Basic local alignment search tool," *J. Mol. Biol.*, vol. 215, no. 3, pp. 403-410, Oct. 1990.

[5] European Grid Infrastructure. http://www.egi.eu
# Déploiement automatique d'une infrastructure complexe pour le mapping des données de séquençage

Sandrine PERRIN¹, Bryan BRANCOTTE¹, Jonathan LORENZO¹, Christophe BLANCHET¹
and Jean-Francois GIBRAT¹

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

Corresponding Author: sandrine.perrin@france-bioinformatique.fr

L'accroissement des données de NGS implique que les équipes disposent de solutions performantes, souples et reproductibles pour le traitement de données volumineuses. De telles infrastructures sont très lourdes pour être proposées au niveau local. L'Institut Français de Bioinformatique (IFB)[1] est une infrastructure nationale de service en bio-informatique pour les sciences de la vie. Elle met à disposition un cloud académique et anime un réseau national de plateformes bioinformatiques. Cette e-infrastructure de calcul est l'environnement adéquat pour répondre aux besoins de traitement des données biologiques. Elle met à disposition des machines virtuelles pré-configurées, des clusters sur-mesure et prochainement des gestionnaires de données volumineuses dont les données de référence.

Le projet européen CYCLONE H2020 [2] a pour mission d'aider aux développements de solutions pour déployer facilement des infrastructures complexes, il apporte notamment des utilitaires pour la construction de VPN et une authentication sécurisée.

La solution présentée ici répond à un cas d'utilisation proposé par CYCLONE sur une solution prenant en charge l'étape de mapping, qui est consommatrice en ressources mémoire (chargement des index), en ressources de calcul et surtout en temps (taille des données d'entrée). Cette étape est préalable à tout traitement sur les données de séquençage. Elle est facilement automatisable. La solution développée par l'IFB propose un déploiement en un clic d'un cluster configurable en fonction des besoins de l'analyse à réaliser et prêt à l'emploi, elle comprend :

* un cluster swarm docker [3] extensible et configurable ;

* un accès à un gestionnaire de données (tel que iRODS), performant avec des volumes importants, à la charge de l'utilisateur de télécharger ces données dans son espace sur le cloud ;

* une interface graphique de gestionnaire de conteneurs Docker, l'analyse sera faite par des outils dockerisés référencés, par exemple dans BioShadock [4], BioContainer ou Docker-hub. Dans notre cas, l'outil TopHat2 [5] est utilisé pour comparer les temps de traitement avec les résultats obtenus classiquement. D'autres étapes peuvent être intégrées : par exemple un pipeline de contrôle qualité ;

* un accès à une banque de données de référence gérée avec BioMaj [6] ;

* un réseau virtuel privé (VPN) pour isoler le cluster et l'accès aux données, la solution CNSMO a été développée dans le cadre du projet CYCLONE ;

* un accès sécurisé au cluster grâce à la fédération d'identités eduGAIN pour les accès SSH et Web ;

* au terme de l'étape de mapping, l'analyse peut être poursuivie avec d'autres outils mis à disposition dans le cloud au sein du cluster ou dans d'autres machines virtuelles, disponibles dans le catalogue RAINBio [7].

L'IFB-core travaille parallèlement à la mise en place d'une fédération de clouds avec les plateformes bioinformatiques françaises partenaires. Les utilisateurs pourront alors rapprocher au mieux les environnements de calcul de leurs données.

[1] Institut Français de Bioinformatique : http://www.france-bioinformatique.fr/

[2] Projet CYCLONE (EU H2020 - 644925), http://www.cyclone-project.eu/

[3] Docker & Swarm : https://docs.docker.com/engine/swarm/ ; https://www.docker.io/

[4] Moreews F, Sallou O, Ménager H *et al*. BioShaDock: a community driven bioinformatics shared Docker-based tools registry. *F1000Research* 4:1443 2015.

[5] Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. *Genome Biology* (14:R36) 2013.

[6] BioMAJ, http://biomaj.genouest.org

[7] RAINBio. https://biosphere.france-bioinformatique.fr/catalogue/
State of the art and comparison of long reads technologies

Anais POIRAUDEAU¹, Maxime MANNO¹, Céline VANDECASTEELE¹, Alain ROULET¹, Celine ROQUES¹, Marie VIDAL¹, Catherine ZANCHETTA¹, Pauline HEUILLARD¹, Fabrice ROUX², Baptiste MAYJONADE², Jérôme GOUZY², Yann GUIGUEN³, Christophe KLOPP⁴, Pierre FRASSE⁵, Mohamed ZOUINE⁵, Cécile DONNADIEU¹, Olivier BOUCHEZ¹, Gérald SALIN¹ and Claire KUCHLY¹

1 INRA, US 1426, GeT-PlaGe, Genotoul, Castanet-Tolosan, France
2 LIPM, UMR 441 INRA, 31326, Castanet-Tolosan, France
3 INRA, UR 1037 LPGP Laboratoire de Physiologie et Génomique des Poissons, 35000, Rennes, France.
4 MIAT, UR875 INRA, Plateforme bioinformatique, 31326, Castanet-Tolosan, France
5 ENSAT, UMR 990 INRA/INPT-ENSAT, Laboratoire Génomique et Biotechnologie des Fruits, 31326, Castanet-Tolosan, France

Corresponding Author: maxime.manno@inra.fr

Genomic issues such as complex genome assembly, structural variant discovery or phasing can be addressed by Long Read technologies.

Thanks to its experience on short reads sequencing using the Illumina technologies, the GeT-PlaGe core facility began to evaluate and use long read technologies since the beginning of 2015 : Pacific BioSciences RSII, Oxford Nanopore Technology MinION and 10XGenomics Chromium. Results from different genomes sequenced by those platforms and sequences obtained on a PacBio Sequel from another France Genomic core facility will be presented and compared on this poster.

As DNA quality is the most important requirement to obtain an efficient sequencing, sample requirements for each technology, and quality controls performed on GeT-PlaGe will be detailed in a first point. Secondly, we will compare sequencing data in terms of file format, sequence metrics, barcoding solutions and needed IT resources. Current projects will also be presented concerning de novo assembly results obtained using Long Read technologies, for several genomes (bacteria, plant and fish).

Acknowledgements

We thank the GenoToul bioinformatics facility for its support in computing resources and data storage, and for helping us processing the data and making available all the needed software and infrastructure.
Genomic markers of species diversification in vertebrates

Guillaume LOUVEL¹, Eric LEWITUS¹, Hélène MORLON¹ and Hugues ROEST CROLLIUS¹

Institut de biologie de l'École Normale Supérieure, 46 rue d'ulm, 75005, Paris, France

Corresponding author: guillaume.louvel@ens.fr

Abstract

Several mechanisms have been hypothesized to explain species divergence, from genomic incompatibilities to divergent selection pressures [1]. Given the current availability of full genomes for many non-model organisms sampling various branches of the vertebrate phylogeny, we can now combine genomic data with patterns of speciation from more complete phylogenies [2]. As Ohno [3] initially postulated, duplicated genes are good candidates for generating functional novelty and adaptation. We first dated duplications using dS (synonymous substitution rate) calculations, in order to obtain a fine estimation of the rate of gene duplication through time and lineages. This method however is sensitive to multiple bias (fast evolving branches, quality of gene alignments, etc). We are currently working on improving these dating method, and aim towards more sophisticated models of gene evolution that could estimate duplication ages (either adapting existing models or developping one). Our broader aim is to compare duplication rates with the diversification rates of taxons, and to assess the role of gene duplication in evolution.

References

[1] J A Coyne and H A Orr. The evolutionary genetics of speciation. *Philosophical transactions of the Royal Society of London. Series B, Biological sciences*, 353(1366):287-305, feb 1998.

[2] Eric Lewitus and Hélène Morlon. Natural Constraints to Species Diversification. *PLOS Biology*, 14(8):e1002532, aug 2016.

[3] S Ohno. *Evolution by Gene Duplication*. Springer-Verlag, 1970.
# PREMS / ELVIS : A local plant biological resource management system

Fabrice DUPUIS¹, Aurélie LELIÈVRE¹, Sandra PELLETIER¹, Tatiana THOUROUDE¹, Julie BOURBEILLON¹ and
Sylvain GAILLARD¹

¹IRHS, Agrocampus-Ouest, INRA, Université d'Angers, SFR 4207 QuaSaV, 49071, Beaucouzé, France

Corresponding Author: sylvain.gaillard@inra.fr

The "Institut de Recherche en Horticulture et Semences" (IRHS) hosts some collections of biological resources for bacteria, fungi and plants. These collections are part of certified Biological Resource Centers (BRC): the "Collection Française de Bactéries associées aux Plantes" (CFBP) for microbes and more recently the "Collection FRuits À PÉpins et Rosiers" (FrAPeR) and the "Collection Apiaceae" for the plants. The management of these last two has leaded to the development of a tool for tracking the material and logging phenotyping and characterization of these resources.

We have developed a webservice oriented database (ELVIS [1]): Expérimentations en Laboratoire Végétal Information System) coupled with a dynamic web interface (PREMS [2]: Plant REsources Management System). These tools take care of specificities linked to the material maintenance (seeds for annual / biannual plants or vegetative reproduction for perennial plants) as well as confidentiality for breeding program.

## Featured functionalities:

- Orchard / field location (apple tree, rose...) or geolocation (wild material characterization) through multi-referential tagging

- Material characterization (plant phenotyping in field, greenhouse, laboratory)

- Material input / output tracking

- Material sample and exploitation tracking (sample, location, transformation)

- Document association (MTA, protocols)

As a webservice oriented application, ELVIS can implement API to communicate with other systems like SIReGal [3], the national portal of INRA for plant genetic resources and through this be a part of the European network for genetic resource management. Moreover, the local installation of PHIS, the PHENOME [4] information system, in the laboratory allows us to build interoperability bridges.

ELVIS is built with Python and PostgreSQL and is freely available under open source CeCILL license.
The webservice API uses the JSON-RPC protocol.

PREMS is powered by Qooxdoo [5] and also available under open source CeCILL license.

## References

[1] https://sourcesup.renater.fr/projects/elvis/

[2] https://sourcesup.renater.fr/projects/prems/

[3] http://prodinra.inra.fr/record/216501

[4] https://www.phenome-fppn.fr/

[5] http://www.qooxdoo.org/
# Découverte et analyse de polymorphismes SNPs issus de RNA-seq chez le peuplier noir

Odile ROGIER¹, Souhila AMANZOUGARENE¹, Marie-Claude LESAGE-DESCAUSES¹, Sandrine BALZERGUE²,³,
Véronique BRUNAUD²,³, José CAIUS²,³, Aurélien CHATEIGNER¹, Ludivine SOUBIGOU-TACONNAT²,³,
Véronique JORGE¹ and Vincent SEGURA¹

¹ INRA, UR0588, UAGPF Amélioration Génétique et Physiologie Forestières, F-45075 Orléans,
France

2 Institute of Plant Sciences Paris Saclay IPS2, CNRS, INRA, Université Paris-Sud,
Université Evry, Université Paris-Saclay, Bâtiment 630, 91405 Orsay, France

3 Institute of Plant Sciences Paris-Saclay IPS2, Paris Diderot, Sorbonne Paris-Cité,
Bâtiment 630, 91405 Orsay, France

Corresponding Author: odile.rogier@inra.fr

Le peuplier noir (*Populus nigra*) est une espèce majeure de la forêt alluviale en Europe occidentale et un support de la biodiversité écosystémique. Pour explorer les composantes de la variabilité génétique au sein de cette espèce parentale de peupliers hybrides commerciaux, nous avons initié une approche intégrative combinant des données génomiques, transcriptomiques et phénotypiques dans une vaste collection d'individus échantillonnés à partir de populations naturelles et évalués dans un même site expérimental. De par leur grande abondance et facilité d'accès par séquençage, les SNPs (Single Nucleotide Polymorphisms) sont des marqueurs moléculaires de choix pour les études de génétique quantitative (incluant les études d'association), de génétique des populations ou de sélection génétique. Chez le peuplier noir, si les premières études de polymorphismes de séquence se sont focalisées sur le reséquençage de quelques gènes candidats [1,2,3] ou sur l'identification de variants rares [4], des travaux plus récents se sont intéressés à une analyse plus globale avec le développement d'une puce de génotypage à partir de SNPs détectés par séquençage de génomes complets [5]. Nous rapportons ici l'analyse de polymorphismes SNPs dans des séquences transcrites produites par RNA-seq. Plus précisément, nous présentons un pipeline de détection et de typepage de SNPs à partir de données RNA-seq ainsi que sa validation grâce à des données déjà existantes.

L'ARN a été extrait à partir de jeune xylème et de cambium recueillis sur 24 arbres correspondant à 2 répétitions de 12 génotypes originaires de 6 populations naturelles représentatives de l'aire de distribution de l'espèce en Europe occidentale. Après quantification de l'ARN par tissu, les ARNs de xylème et de cambium du même arbre ont été mélangés et soumis à un séquençage à haut débit sur un HiSeq2000 d'Illumina ( séquences pairées de 100 pb). Nous avons mis en place un pipeline de détection de SNPs à partir de données RNA-seq. Les séquences ont d'abord été nettoyées puis alignées sur le génome de référence (*Populus trichocarpa*). Ensuite, nous avons effectué un post-traitement des données d'alignement et nous avons lancé la découverte de variants à l'aide 3 outils en plusieurs modes, combinés à plusieurs paramètres de filtration. Les résultats obtenus par ces différentes techniques ont été comparés. Enfin, les SNPs détectés par les différents outils ont été annotés. Pour valider la technique, nous avons comparé les génotypes trouvés à des résultats de génotypage issus d'une puce Illumina Infinium 12k fournissant 7 903 SNPs polymorphes et qui avait été utilisée pour génotyper une collection de peupliers noirs incluant nos 12 génotypes étudiés [5]. Le taux moyen de similarité de génotypage aux positions communes entre les SNPs de la puce et ceux trouvés par RNA-seq (d'un même individu) est supérieur à 98%. Ces résultats démontrent la faisabilité et l'intérêt du RNA-seq pour la découverte et le typepage de SNPs dans des populations naturelles.

## References

[1] Fabio Marroni *et al*. Large-Scale Detection of Rare Variants via Pooled Multiplexed next-Generation Sequencing: Towards next-Generation Ecotilling: Detection of Rare Variants via next-Generation Sequencing. *The Plant Journal* 67 (4): 736–45, 2011.

[2] Fabio Marroni *et al*. Nucleotide Diversity and Linkage Disequilibrium in *Populus nigra* Cinnamyl Alcohol Dehydrogenase (CAD4) Gene. *Tree Genetics & Genomes* 7 (5): 1011–23, 2011.

[3] Fernando P. Guerra *et al*. Association Genetics of Chemical Wood Properties in Black Poplar (*Populus nigra*). *New Phytologist* 197 (1): 162–76, 2013.

[4] Bartel Vanholme *et al*. Breeding with Rare Defective Alleles (BRDA): A Natural *Populus nigra* HCT Mutant with Modified Lignin as a Case Study. *New Phytologist* 198 (3): 765–76, 2013.

[5] Patricia Faivre-Rampant *et al*. New resources for genetic studies in *Populus nigra*: genome-wide SNP discovery and development of a 12k Infinium array. *Molecular Ecology Resources* 16 (4): 1023–1036, 2016.
# DiNAMO: Exact method for degenerate IUPAC motifs discovery, characterization of sequence-specific errors

Chadi SAAD<fnref n="1" />,<fnref n="2" />, Laurent NOË<fnref n="2" />, Hugues RICHARD<fnref n="3" />, Julie LECLERC<fnref n="1" />, Marie-Pierre BUISINE<fnref n="1" />, Hélène TOUZET<fnref n="2" /> and Martin FIGEAC<fnref n="4" />

Corresponding author: Chadi.Saad@univ-lille1.fr

Next generation sequencing technologies are still associated with relatively high error rates, about 1%, which correspond to thousands of errors in the scale of a complete genome. Each region needs therefore to be sequenced several times and variants are usually filtered based on depth criteria. The significant number of artifacts, in spite of those filters, shows the limit of conventional approaches and indicates that some sequencing artifacts are recurrent. This recurrence underlines that sequencing errors can depend on the upstream nucleotide sequence context. Our goal is to search for overrepresented motifs that tend to induce sequencing errors.

Previous studies showed that some motifs, such as GGT [1,2], induce sequencing errors in the Illumina tech-
nologies. However, these studies were dedicated to exact motifs, and did not take into account approximate
motifs, limiting the statistical power of such approaches. On the other hand, some tools, such as FIRE [3],
DREME [4] and Discrover [5], were developed to search for degenerate motifs over the 15-letter IUPAC al-
phabet in the context of chip-seq studies. However, these tools use greedy algorithms, implying a lack of
sensitivity. So we developed an exact algorithm to search for degenerate motifs by enumerating all possible
IUPAC motifs. This algorithm is based on mutual information and uses hashtables with graphs data structure
to store the motifs. It is independent from the sequencing technology.

Experimental results on real data show that there are many overrepresented motifs upstream of sequencing artifacts. These latter are identified through the strand bias between forward and reverse reads. The homopoly- mer of length 3 CCC seems to be sufficient to induce errors on IonTorrent. On Illumina, motifs are mainly composed of GGC followed by GGT (like: TGGCNGGT) or homopolymers. We have also noticed a base quality fall after the detected motifs. Our exact algorithm requires less than one minute (Intel® Core™ i5-4570 CPU, 3.20GHz), and less than 2GB of RAM to search for full degenerate motifs of length 6 on a dataset of approximately 24000 sequences, extracted from 11 exomes sequenced on IonTorrent Proton.

Availability: https://github.com/bonsai-team/DiNAMO

Acknowledgements

This work is supported by Lille University Hospital and Hauts-de-France region

References

[1] Manuel Allhoff, Alexander Schönhuth, Marcel Martin, Ivan G Costa, Sven Rahmann, and Tobias Marschall. Discovering motifs that induce sequencing errors. In *BMC bioinformatics*, volume 14, page S1. BioMed Central Ltd, 2013.

[2] Frazer Meacham, Dario Boffelli, Joseph Dahbi, David IK Martin, Meromit Singer, and Lior Pachter. Identification and correction of systematic error in high-throughput sequence data. *BMC bioinformatics*, 12(1):451, 2011.

[3] Olivier Elemento, Noam Slonim, and Saeed Tavazoie. A universal framework for regulatory element discovery across all genomes and data types. *Molecular cell*, 28(2):337–350, 2007.

[4] Timothy L Bailey. Dreme: motif discovery in transcription factor chip-seq data. *Bioinformatics*, 27(12):1653–1659, 2011.

[5] Jonas Maaskola and Nikolaus Rajewsky. Binding site discovery from nucleic acid sequences by discriminative learning of hidden markov models. *Nucleic acids research*, 42(21):12995–13011, 2014.

<footnote n="1">¹ JPARC (UMR1172 Inserm, Lille 2 University and Lille University Hospital ), 59000, Lille, France</footnote>

<footnote n="2">² CRISTAL (UMR CNRS 9189 Lille 1 University and Inria Lille), Team BONSAI, 59000, Lille, France</footnote>

<footnote n="3">³ LCQB (UMR 7238 CNRS Pierre Marie Curie University ), 75006, Paris, France</footnote>

<footnote n="4">⁴ Functional and Structural Genomic Platform, Lille 2 University, 59000, Lille, France</footnote>