Claude Université

Bernard-Lyon

M François-Noël

M Le Professeur Hamda

Ben Hadid

M Le Professeur

Germain Gillet

M Alain Helleu

Composantes Sante

J Etienne

La Mme

C Professeure

Burillon

D Bourgeois

Vinciguerra

Y Matillon

Ana Tereza

Ribeiro De Vasconcelos

Nuno Mendes

Olivier Rue

Hubert Charles

Federica Calevro

Karen Gaget

Gabrielle Duport

Susana Vinga

Nadia Pisanti

Roberto Grossi

Caio Padoan

Sá Godinho

It has been extremely rewarding to work and learn with all of them.

COMPOSANTES ET DEPARTEMENTS DE SCIENCES ET TECHNOLOGIE

CONTENTS

Introduction

This thesis mainly addresses methodological problems related to the prediction of small regulatory RNAs, specially microRNAs (miRNAs). The first topic involved the elaboration of a robust and efficient method, called Mirinho, for the prediction of pre-miRNAs in both genomic and small RNA sequencing (sRNA-seq) data of both animal and plant species. The second topic led to the development of a pipeline, called MirinhoPipe, for the treatment of small RNA sequencing data. It was specially implemented to identify the expressed miRNAs of Acyrthosiphon pisum, the pea aphid. We then moved on to solve a few problems related to the prediction and analysis of non-coding RNAs (ncRNAs) in the bacterium Mycoplasma hyopneumoniae. A method, called Alvinho, was thus developed for the prediction of targets in this bacterium, together with a pipeline for the segmentation of a numerical sequence and detection of conservation among ncRNA sequences using a k-partite graph. We finally addressed a problem related to motifs, that is to patterns, that may be composed of one or more parts, that appear conserved in a set of sequences and may thus correspond to functional elements such as DNA binding sites or miRNA families (i.e., all the isoforms of a same miRNA). We presented some clustering solutions to group the motifs that may correspond to a same such biological element, and thus to better distinguish the biologically significant ones from noise that may be present in what often are large outputs from many motif extraction algorithms.

All the methodological and biological concepts required to understand the previous topics are presented in Chapter 1. We now provide a brief introduction to each of these topics.

Given the importance and ubiquity of miRNAs in a wide range of biological processes and diseases, a plethora of methods for the prediction of miRNAs were developed. Despite all the effort put in developing them, there remained a number of issues that needed to be addressed:

1. the vast majority of the existing softwares rely on a folding algorithm of cubic time complexity to predict the characteristic hairpin structure of a pre-miRNA: this is suitable when the input is small enough, but it can become impracticable when the size of the input increases;

2. for longer pre-miRNAs (such as in plant), such folding methods moreover can produce hairpin structures different from the ones provided in miRBase [START_REF] Kozomara | mirbase: integrating microRNA annotation and deep-sequencing data[END_REF], as a consequence the miRNA may be located in a different place than a stem-arm;

3. together with folding, most methods then rely on further information that must be learned from previously validated miRNAs of closely related genomes (at a minimum within the same clade, plant or animal) for the final prediction of new miRNAs in order either to set the parameters of the model or to restrict the search to a limited space.

Mirinho was therefore developed to address all three issues. The search for pre-miRNAs is concentrated on regions with the same length as the two stem-arms separated by the length of the terminal loop. The direct application to sRNA-seq data guarantees a better quality in the prediction of the pre-miRNA structures. A quadratic time complexity algorithm improves the practical efficiency of the free energy computation. As neither of the two attributes used (length of stem-arm and terminal loop) are species-specific within the animal or the plant kingdom (they differ only between these two kingdoms), the method can easily be applied for predicting pre-miRNAs in either clade. Importantly, while the method we provide is thus much simpler, faster, and general to use, we also show for tested examples that it has a sensitivity and precision as good as other methods, in some cases even better. Moreover, we show that the secondary structures predicted by Mirinho are much closer to the ones available in miRBase than for the other compared methods. Mirinho is described in detail in Chapter 2 which is strongly based on our paper Higashi et al. (ress).

Still concerning the identification of miRNAs, however from another perspective, we treated and analysed the small RNA sequencing (sRNA-seq) data of Acyrthosiphon pisum, the pea aphid. The unique feeding habit of aphids combined with their ability to rapidly reproduce makes of them one of the most damaging pests of crops with economical importance worldwide. Considering their impact on agriculture and the role miRNAs play in gene regulation, it is imperative to better characterise and understand the function of these miRNAs. One first effort has already been made by Legeai et al. (2010a) in A. pisum, a laboratory model for the study of these pests whose genome was sequenced. It is worth noting that in Legeai's work, the miRNAs of parthenogenic females were sequenced and analysed, while we focus on the miRNAs expressed in three embryonic developmental and one larval stages. Furthermore, the potential mRNA targets of the detected miRNAs were identified by the overlapped predictions of two methods (Pita and miRanda), and correlated with the gene expression profile of the pea aphid.

To treat the data in order to guarantee a more accurate set of reads, as well as to detect the expressed miRNAs, three approaches were used: (i) MirinhoPipe, specially developed for this analysis; (ii) sRNA-PlAn, a pipeline designed for the annotation of small RNAs; and (iii) miRDeep, a classical method for the discovery of miRNAs from deep sequencing data [START_REF] Friedländer | Discovering microRNAs from deep sequencing data using mirdeep[END_REF]. The detected miRNAs were submitted to the prediction of mRNA targets. Together with such predictions, the gene expression profile of A. pisum was analysed and compared to the miRNA expression profile, leading to very interesting results. All the methodology, results and discussion are presented in Chapter 3 which is strongly based on our paper Higashi et al. (tion).

Besides miRNAs, another small regulatory molecule was investigated. This was noncoding RNAs in Mycoplasma hyopneumoniae. The bacterium M. hyopneumoniae strain 7448 is a pathogenic and obligate parasite of porcine respiratory systems. It lives adhered to the epithelium of its host respiratory tract, and together with other bacteria and viruses, it is considered one of the ethiologic agents of swine enzootic pneumonia. The disease can cause a decrease in the productivity of these animals, sometimes resulting in their death [START_REF] Byrt | Effect of enzootic pneumonia of pigs on growth performance[END_REF][START_REF] Debey | Ciliostasis and loss of cilia induced by Mycoplasma hyopneumoniae in porcine tracheal organ cultures[END_REF][START_REF] Brockmeier | Polymicrobial Diseases, chapitre 13, Porcine respiratory disease complex[END_REF]. Although some effort has already being put in understanding the infection process, the specific mechanisms relating the bacterium to the disease remain unknown [START_REF] Gardner | Detection and quantification of intergenic transcription in Mycoplasma hyopneumoniae[END_REF][START_REF] Hsu | Molecular analysis of the p97 cilium adhesin operon of Mycoplasma hyopneumoniae[END_REF][START_REF] Nicolás | Abc transporters in Mycoplasma hyopneumoniae and Mycoplasma synoviae: insights into evolution and pathogenicity[END_REF][START_REF] Siqueira | New insights on the biology of swine respiratory tract mycoplasmas from a comparative genome analysis[END_REF].

M. hyopneumoniae 7448 has only one known transcription factor (TF) and a complex gene expression pattern. The incomparability between the number of regulatory elements and the complexity of the gene expression of the bacterium, together with increasing evidences that ncRNAs are involved in this phenomenon, strongly encourage the search for ncRNAs in the genome of M. hyopneumoniae 7448. After predicting the regions with a potential to harbour ncRNA genes, additional analyses were performed in an attempt to provide more evidences to carry on with experimental validation of the ncRNAs.

The first problem that concerned us was related to the output of the pipeline for the prediction of ncRNAs: such pipeline was generating one single assembled ncRNA sequence where two or more different ncRNA candidates were in fact present. We solved this by applying a segmentation algorithm on these outputs. To then provide stronger evidence that the candidates were indeed functional, we performed the prediction of the ncRNA targets with a method, called Alvinho, that was specially developed for this purpose. Finally, to verify if conservation could play any role in the functionality of ncRNAs, the identity of intergenic regions was assessed between closely-related Mycoplasma species by means of a kpartite graph. Genomic motifs surrounding the ncRNA, such as promoters and terminators, were also verified to reinforce the functional evidence of the ncRNA candidates. All the three steps of the pipeline are available in the form of a script or a C++ implementation. All the details concerning the methods developed for the analysis of the ncRNAs are presented in Chapter 4 that is based on the paper Godinho et al. (tion).

We then looked at a problem related to structured motifs, which corresponds to a possibly complex pattern that is conserved in a sequence or a set of sequences. This is an issue that may seem unrelated to the study of miRNAs but the two may however appear combined in some studies. For instance, the motifs associated to the miRNAs that are exported from a human tissue might enable to understand what distinguishes such miRNAs from those that are not exported.

The problem of finding structured motifs was first addressed by Marsan et Sagot (2000) and implemented as a software called Smile (Structured Motifs Inference and Evaluation). Depending on the parameters given to Smile, the algorithm can generate a large output that may contain redundant information. This will happen in particular when the characteristics of the motifs are not precisely known, thus requiring that more permissive parameters are adopted in an attempt to recover them. We therefore present some clustering solutions to group together motifs that may correspond to a same biological "object", and to better identify the noise that may be present in such large outputs.

Efficiently extracting consensus sites in a set of sequences is an essential approach to identify functional elements in a genome. Examples of such elements are DNA binding sites and miRNA families (i.e., a consensus that represents all the miRNA isoforms). There are two main problems related to this identification. One is the prediction of the location of the element site, and the second is the extraction of the consensus. The algorithm Smile (Marsan et Sagot, 2000) addresses both problems: extracting and locating consensus motifs in a set of sequences. To solve this problem, Smile implements an exact algorithm for finding motifs in a set of sequences. A suffix tree is used to represent the input sequences, which together with the strategies implemented in the algorithm, result in an efficient method for the extraction of motifs. Smile requires a number of parameters, such as the number p of boxes a (structured) motif may have, the minimum number of substitutions e (one per box) between the motif and its occurrence, and the minimum number of times q (stands for quorum) the motif has to appear among the sequences.

Depending on the values of these parameters, the size of the output generated by Smile may be very large, containing redundant motifs. For example, the larger is the number e or the smaller the quorum q, the larger will be the output. In an attempt to organise such output eliminating the redundancy, we implemented an UPGMA (Unweighted Pair Group Method with Arithmetic Mean) algorithm to cluster the similar motifs according to the positions where they appear. The implementation of this algorithm was performed during the internship of Thomas Balezeau, an undergraduate student in information technology whom I co-advised together with Marie-France Sagot. Another approach that has been explored, but not yet implemented, is the use of hashing for list intersection as an estimator to find redundant motifs. All the details concerning structured motifs and the clustering approaches are presented in Chapter 5. In this chapter, we present the biological and computational backgrounds required to the comprehension of this thesis. It is certainly not possible to cover all the details about the concerned topics to provide a self-contained thesis; we therefore provide only the concepts that we find crucial for both the computational and biological sides. The chapter is divided in two sections: Section 1.1 presents the biological concepts and Section 1.2 covers the methodological concepts (computational and experimental).

The purpose of Section 1.1 is to present microRNAs (miRNA) and the involved machinery. We thus begin by defining a miRNA and by providing its historical background and current landscape to place miRNAs in a small regulatory context; from this exposition, one should be convinced of the importance of miRNAs in the different biological processes in which they are involved. We then present the miRNA biogenesis process, an important issue since any computational modelling of a miRNA is strongly based on this process. We thus address the RNA-induced silencing complex (RISC) that is responsible for the functional regulatory interaction between a miRNA and its target messenger RNA (mRNA).

The main goal of Section 1.2 is to introduce the current computational and experimental methods used to detect miRNAs and targets. We first present how the two problems, prediction of miRNAs and prediction of targets, are computationally addressed. We then introduce the experimental methods used to detect miRNAs and targets, an important aspect that may complement and validate the results obtained by the computational methods.

Biological background

MicroRNA definition, history, and landscape

A miRNA is a small non-coding regulatory molecule, present in animals, plants, and in a few viruses. It is responsible for the post-transcriptional regulation of gene expression via complementarity base-pairing with the target mRNA; frequently the result of the regulation is the silencing of the target, however, there are fewer cases in which the expression is enhanced. These transcripts of ∼22nt are derived from a precursor-miRNA (pre-miRNA) with a specific hairpin (stem-loop) structure, with small internal loops and bulges, and are located in the stem of the hairpin. A "classical" miRNA would meet all the previous features. Although in practice variations are obviously possible, the minimum requirement to classify a sequence as a miRNA is its length (∼22nt) and the presence of a hairpin loop [START_REF] Berezikov | Approaches to microRNA discovery[END_REF][START_REF] Chen | The evolution of gene regulation by transcription factors and microRNAs[END_REF][START_REF] He | MicroRNAs: small RNAs with a big role in gene regulation[END_REF].

These molecules are believed to be involved in the regulation of several basic pathways, such as in the transition of developmental stages in nematodes (lin-4 and let-7 ) [START_REF] Reinhart | The 21-nucleotide let-7 RNA regulates developmental timing in caenorhabditis elegans[END_REF][START_REF] Wightman | Posttranscriptional regulation of the heterochronic gene lin-14 by lin-4 mediates temporal pattern formation in c. elegans[END_REF][START_REF] Lee | The c. elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14[END_REF], cell proliferation and apoptosis (miRNA Bantam) [START_REF] Brennecke | bantam encodes a developmentally regulated microRNA that controls cell proliferation and regulates the proapoptotic gene hid in Drosophila[END_REF], regulation of fat metabolism (miR-14 ) [START_REF] Xu | The Drosophila microRNA mir-14 suppresses cell death and is required for normal fat metabolism[END_REF], etc. Furthermore, they also known to play a role in diseases such as autoimmune and neurodegenerative diseases, and in cancer [START_REF] Almeida | MicroRNA history: discovery, recent applications, and next frontiers[END_REF].

The first investigations concerning RNA interference (RNAi), which is a process of inhibiting gene expression (i.e., of gene silencing), started in 1990 with the efforts of two teams [START_REF] Napoli | Introduction of a chimeric chalcone synthase gene into petunia results in reversible co-suppression of homologous genes in trans[END_REF][START_REF] Van Der Krol | Flavonoid genes in petunia: addition of a limited number of gene copies may lead to a suppression of gene expression[END_REF]. The authors used a transgene in an attempt to over-express an enzyme related to the violet color of petunias. In the end, instead of obtaining darker violet petunias as it was expected, they observed white ones. They thus raised the hypothesis that the endogenous and transgenic genes were co-suppressed.

Three years later, the first miRNA was identified in the nematode Caenorhabditis elegans by [START_REF] Lee | The c. elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14[END_REF] and [START_REF] Wightman | Posttranscriptional regulation of the heterochronic gene lin-14 by lin-4 mediates temporal pattern formation in c. elegans[END_REF]. The authors cloned the gene lin-4 and discovered that it did not encode a protein but instead a small RNA of 21nt. They observed that, by partial complementarity between the miRNA and the 3 UTR of the lin-14 mRNA, the translation of the protein LIN-14 was being repressed.

In 1998, the classically established flow of information inside a cell (the so-called central dogma) became more complex as the pathway of RNAi was first described by [START_REF] Fire | Potent and specific genetic interference by double-stranded RNA in caenorhabditis elegans[END_REF], with RNAs regulating other RNAs, instead of only producing proteins. The authors discovered that, instead of a single strand RNA (ssRNA), the trigger for the gene silencing in Caenorhabditis elegans was a double strand RNA (dsRNA) [START_REF] Sen | A brief history of RNAi: the silence of the genes[END_REF]. This work introduced a new concept for the gene silencing pathway, clarifying the results of previous works [START_REF] Napoli | Introduction of a chimeric chalcone synthase gene into petunia results in reversible co-suppression of homologous genes in trans[END_REF][START_REF] Van Der Krol | Flavonoid genes in petunia: addition of a limited number of gene copies may lead to a suppression of gene expression[END_REF][START_REF] Guo | par-1, a gene required for establishing polarity in c. elegans embryos, encodes a putative ser/thr kinase that is asymmetrically distributed[END_REF], and maybe becoming one of the best known pathways for RNA silencing, since it is possible to repress the expression of a wide range of genes with just partial sequence complementarity.

In 2000, a second miRNA, namely let-7, was discovered in Caenorhabditis elegans by [START_REF] Reinhart | The 21-nucleotide let-7 RNA regulates developmental timing in caenorhabditis elegans[END_REF]. The authors noticed that the repression of let-7 caused the reappearance of larval characteristics during the adult stage, while the over-expression of let-7 caused the early expression of adult characteristics. The authors then concluded that miRNA let-7 was controlling the transition of developmental stages in Caenorhabditis elegans. A timeline with the major discoveries in gene silencing can be found in Figure 1.1.

Gradually other types of small non-coding RNAs were being discovered: Piwi-interacting RNAs (piRNA) [START_REF] Siomi | Piwi-interacting small RNAs: the vanguard of genome defence[END_REF], transcription initiation RNAs (tiRNA) [START_REF] Taft | Small RNAs derived from snoRNAs[END_REF], ./figures/figure_history1.pdf ./figures/figure_history2.pdf ./figures/figure_history3.pdf Figure 1.1: Timeline of the main discoveries in gene silencing [START_REF] Zamore | Ribo-gnome: the big world of small RNAs[END_REF][START_REF] Kunej | Cross talk between microRNA and coding cancer genes[END_REF].

nucleolar RNAs (snoRNA) [START_REF] Filipowicz | Biogenesis of small nucleolar ribonucleoproteins[END_REF][START_REF] Dieci | Eukaryotic snoRNAs: a paradigm for gene expression flexibility[END_REF], and other miRNAs [START_REF] Winter | MicroRNA biogenesis and cancer[END_REF][START_REF] Siomi | Posttranscriptional regulation of microRNA biogenesis in animals[END_REF][START_REF] De Planell-Saguer | Analytical aspects of microRNA in diagnostics: a review[END_REF] have also been identified. The full landscape of such small regulatory RNAs is presented in Figure 1.2.

./figures/figure_landscape.pdf Figure 1.2: RNA landscape and the different types of small non-coding RNAs: transcription initiation RNA (tiRNA) [START_REF] Taft | Small RNAs derived from snoRNAs[END_REF], long non-coding RNA (lncRNA) [START_REF] Geisler | RNA in unexpected places: long non-coding RNA functions in diverse cellular contexts[END_REF][START_REF] Batista | Long noncoding RNAs: cellular address codes in development and disease[END_REF], small interfering RNA (siRNA) [START_REF] Castel | RNA interference in the nucleus: roles for small RNAs in transcription, epigenetics and beyond[END_REF]Davidson et McCray, 2011), Piwi-interacting RNA (piRNA) [START_REF] Siomi | Piwi-interacting small RNAs: the vanguard of genome defence[END_REF], small nucleolar RNA (snoRNA) [START_REF] Filipowicz | Biogenesis of small nucleolar ribonucleoproteins[END_REF][START_REF] Dieci | Eukaryotic snoRNAs: a paradigm for gene expression flexibility[END_REF], small Cajal bodyspecific RNA (scaRNA) [START_REF] Darzacq | Cajal body-specific small nuclear RNAs: a novel class of 2-o-methylation and pseudouridylation guide RNAs[END_REF], X-inactivation RNA (xiRNA) [START_REF] Ogawa | Intersection of the RNA interference and x-inactivation pathways[END_REF], small RNA (sRNA) (Gottesman et Storz, 2011), small-scan RNA (scnRNA) [START_REF] Kim | Small RNAs: classification, biogenesis, and function[END_REF], promoter-associated small RNA (PASR) [START_REF] Kapranov | RNA maps reveal new RNA classes and a possible function for pervasive transcription[END_REF] (image modified from [START_REF] Ghosh | Renaissance of the regulatory RNAs[END_REF]).

MicroRNA biogenesis

Concerning the transcription of miRNA genes, these molecules can arise either from intergenic regions or from introns of spliced genes. They are either transcribed as independent units or in clusters of miRNAs by means of a polycistronic transcript. Figure 1.3 will serve as a support for all the explanation given in what follows.

The transcription of miRNA genes is mainly performed by the RNA polymerase II (Pol II). The Pol II begins the transcription in the nucleus by binding to the promoter, and its first product is a longer transcript (between 500bp and 10Kbp) called primary-miRNA (pri-miRNA), which is capped at the 5 end and polyadenylated at the 3 end. Still in the Figure 1.3: Animal miRNA biogenesis: A miRNA can be located either in an intergenic region or in an intronic region of a protein encoding gene. In the first canonical pathway (left side of the figure), the miRNA is transcribed from its own gene into a pri-miRNA, which will then be processed by Drosha into a pre-miRNA. In the second non-canonical pathway (right side of the figure), the host gene is transcribed, spliced and the miRNA comes from an intron-in this case the miRNA is called mirtron. After debranching, the sequence folds itself into a pre-miRNA, and does not require processing by Drosha. The pre-miRNA is exported by Exportin-5 to the cytoplasm, where it is further processed by the enzyme Dicer into a duplex miRNA:miRNA*; usually only one of the strands is loaded into the RISC complex and the other miRNA* is degraded.

nucleus, the pri-miRNA is processed by a microprocessor (composed of the Drosha RNase III enzyme and its cofactor DiGeorge syndrome critical region gene 8 (DGCR8)) into a ∼80nt stem-loop precursor miRNA (pre-miRNA) with a 2 nucleotides overhang at the 3 end. The pre-miRNA is then exported by the nucleocytoplasmic shuttler Exportin-5 to the cytoplasm, where another RNase III enzyme called Dicer and its cofactor transactivating response RNAbinding protein (TRBP) recognise the 2 nucleotides overhang left by Drosha and cleaves the terminal loop. The result is a short imperfect miRNA:miRNA* duplex of length ∼22nt that is unwounded, producing one functional strand (mature miRNA) and another non-functional miRNA* (miRNA star) that is usually degraded, although sometimes it can be functional too [START_REF] Petersen | 19 the biology of short RNAs[END_REF][START_REF] Yang | Widespread regulatory activity of vertebrate microRNA species[END_REF][START_REF] Okamura | The regulatory activity of microRNA species has substantial influence on microRNA and 3' utr evolution[END_REF]. The duplex configuration is known to stabilise the miRNA by protecting it from RNases degrading single-stranded miRNAs. Recently, [START_REF] Winter | Loop-mirs: active microRNAs generated from single-stranded loop regions[END_REF] provided a first evidence that single-stranded loop regions may give origin to functional regulatory miRNAs, which the authors call loop-miRNAs.

During the splicing of other genes, miRNAs can also arise from introns. After splicing, the intronic region folds into a pre-miRNA stem loop and it is then submitted to the same canonical biogenesis pathway. This kind of miRNA is called mirtron and is independent of the activity of Drosha. After all the processing, the mature miRNA is incorporated into a complex called RNA-induced silencing complex (RISC) to be further driven for target regulation, as it is detailed in the next section [START_REF] Okamura | The mirtron pathway generates microRNA-class regulatory RNAs in Drosophila[END_REF][START_REF] Meister | Argonaute proteins: functional insights and emerging roles[END_REF][START_REF] Gommans | Controlling miRNA regulation in disease[END_REF].

RNA-induced silencing complex

Once the miRNA is assembled into RISC, it anneals to its mRNA target for regulation. The RISC complex is comprised of several proteins, among them one is well known: the Argonaute (AGO) protein that is the component that links the miRNA with the complex by means of two domains, PAZ and PIWI, responsible for miRNA recognition. The RISC assembly may be divided into at least two successive steps: RISC-loading, in which miRNA duplexes are inserted into the AGO proteins; and (ii) strand dissociation, in which the two miRNA strands are separated within the AGO protein. During assembly, the AGO proteins suffer conformation changes, made by chaperones, to allow for the incorporation of the miRNA duplex. Once the duplex is incorporated, AGO releases the tension to recover its original conformation unwinding the duplex and discarding the non-functional miRNA strand called "passenger strand". The remaining functional mature miRNA (or "guide strand") is usually the one with the less stable 5 end. This mechanism of RISC assembly is mostly studied in Drosophila using AGO2-RISC as a model system. Although the exact molecular composition of RISC is unknown, a sufficient requirement for target regulation is the Argonaute protein [START_REF] Kawamata | Making risc[END_REF][START_REF] Meister | Argonaute proteins: functional insights and emerging roles[END_REF][START_REF] Scott | Diversity, overlap, and relationships in the small RNA landscape[END_REF].

Usually, the effect of RISC, miRNA and target interaction is down-regulation, either through the cleavage of the mRNA target or by repression of the translation. To cleave the target, at least two requirements are necessary: an Argonaute with catalytic activity (in humans only AGO2 has this characteristic), and a near-perfect complementarity between the guide miRNA and its target. Different from the cleavage, near-perfect complementarity is not required to repress translation. Instead, only a smaller region of 6nt, that is called seed, requires perfect complementarity. It is usually located at the 5 end of the miRNA (positions 2 to 8) and is known to be more frequent in animals [START_REF] Zheng | Animal microRNA target prediction using diverse sequencespecific determinants[END_REF]. When the pairing at the 5 end is insufficient, stronger pairing at the 3 end compensates for it [START_REF] Brennecke | Principles of microRNA-target recognition[END_REF]. Concerning the mRNA target, the interaction is usually located in the 3 untranslated region (3UTR) [START_REF] Pratt | The RNA-induced silencing complex: a versatile gene-silencing machine[END_REF]. However, studies such as [START_REF] Lytle | Target mRNAs are repressed as efficiently by microRNA-binding sites in the 5 utr as in the 3 utr[END_REF]; [START_REF] Moretti | Mechanism of translational regulation by mir-2 from sites in the 5 untranslated region or the open reading frame[END_REF]; Ørom et al. (2008); [START_REF] Qin | mir-24 regulates apoptosis by targeting the open reading frame (orf) region of faf1 in cancer cells[END_REF]; [START_REF] Fang | The impact of miRNA target sites in coding sequences and in 3 utrs[END_REF] suggested that the association can be functional in the 5 UTR also or even in the CDS.

It has been shown that the repression is even more effective when there are multiple miRNAs binding to the same mRNA target, suggesting that the regulation is controlled by multiple miRNAs [START_REF] Fang | The impact of miRNA target sites in coding sequences and in 3 utrs[END_REF][START_REF] Krek | Combinatorial microRNA target predictions[END_REF][START_REF] Grimson | MicroRNA targeting specificity in mammals: determinants beyond seed pairing[END_REF]. One more characteristic that was found to contribute to the regulatory effect is the AU content in the 3 of the seed region [START_REF] Jing | Involvement of microRNA in au-rich elementmediated mRNA instability[END_REF][START_REF] Grimson | MicroRNA targeting specificity in mammals: determinants beyond seed pairing[END_REF]. Figure 1.4 shows a schema of the interaction between a miRNA and its target.

./figures/figure_target.pdf Figure 1.4: Mechanisms of interaction between miRNA and target. The interaction of miRNA and target can cause either mRNA degradation, when the complementarity between miRNA and target is near perfect, or repression of translation, when there is partial sequence complementarity in relation to the whole miRNA sequence. A basic feature required for the interaction is seed, a region starting at the second nucleotide of the 5' end of the miRNA with perfect base pairing. In relation to the target, the interaction frequently occurs at the 3'UTR. On the right side of this figure, the miRNAs act in a synergistic way: multiple miRNAs (the ones in red and blue) bind to the same target to cooperatively regulate it. The figure was taken from [START_REF] Sun | microRNA: a master regulator of cellular processes for bioengineering systems[END_REF].

Although regulation by miRNAs has been widely studied, a model describing in detail the mechanisms of the different modes of actions is still being debated [START_REF] Lytle | Target mRNAs are repressed as efficiently by microRNA-binding sites in the 5 utr as in the 3 utr[END_REF][START_REF] Moretti | Mechanism of translational regulation by mir-2 from sites in the 5 untranslated region or the open reading frame[END_REF][START_REF] Ørom | MicroRNA-10a binds the 5 utr of ribosomal protein mRNAs and enhances their translation[END_REF][START_REF] Qin | mir-24 regulates apoptosis by targeting the open reading frame (orf) region of faf1 in cancer cells[END_REF][START_REF] Fang | The impact of miRNA target sites in coding sequences and in 3 utrs[END_REF].

Plant microRNAs

For the sake of concision, we will highlight just the differences between animal and plant miRNAs. Starting by the transcription, it seems that the great majority of plant miRNAs are produced from their own transcription units, while animal miRNAs can also be produced from introns of spliced genes (see Section 1.1.2). Just like in animals, plant miRNAs can also appear in clusters. However, this polycistronic organisation is much more frequent in animals than in plants. In plants, instead of requiring two different enzymes to process pri-miRNA and pre-miRNA (Drosha and Dicer respectively), it seems that Dicer-Like 1 (DCL1) performs both roles in the nucleus, producing in the first step longer and more variable stemloop pre-miRNAs. Once processed, the miRNA:miRNA* duplex is exported to the cytoplasm by the transporter HASTY. It will then be loaded into the RISC complex, and by nearperfect complementarity with its target, it will induce endonucleolytic cleavage of the mRNA. Translation repression, result of a weak base paring, is yet to be explored in plants. Target sites can be located either in coding exons or in 3' UTRs [START_REF] Jones-Rhoades | MicroRNAs and their regulatory roles in plants[END_REF][START_REF] Axtell | Vive la difference: biogenesis and evolution of microRNAs in plants and animals[END_REF][START_REF] Rogers | Biogenesis, turnover, and mode of action of plant microR-NAs[END_REF][START_REF] Pasquinelli | MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship[END_REF].

Methodological background

miRBase: a reference for microRNA studies

To begin this second section, we present an important miRNA resource that serves as a reference for miRNA research. Most of the miRNA studies use miRBase as a gold standard.

To validate the miRNA predictions, either an experimental method must be performed or a gold standard must be adopted. Given that the first option is much more expensive, the great majority of the authors use miRBase as a reference for the validation. miRBase is a simple and practical database for published miRNA sequence and annotation. It is currently the main source for miRNA annotation with frequent updates. As a consequence, the number of annotated miRNAs grows exponentially (see Figure 1.5), reaching in 2014 (release 21) 28,645 annotated miRNAs. miRBase is available at http://www.mirbase.org/. A miRNA gene is identified by a prefix (3 or 4 letters) corresponding to the species, followed by the sequence type (mature or precursor miRNA), with a sequential number at the end. For instance, the identifier "hsa-miR-101" corresponds to the human (hsa) mature miRNA sequence (miR) 101. While mature sequences are identified by "miR", precursor miRNA sequences are labelled "mir". Some more informations are aggregated to the name of the mature miRNA. For instance, on the same example as before, miRNAs "hsa-miR-101" and "mmu-miR-101" are similar genes appearing in different species, while "hsa-miR-101a" and "hsa-miR-101b" are similar genes (of the same species) differing at one or two bases. Furthermore, if the same miRNA sequence arises from different pre-miRNA loci, a numbered suffix is added to the end of the miRNA name. For instance, "dme-mir-281-1" and "dmemir-281-2" are two identical miRNA sequences from Drosophila melanogaster derived from different positions of the same pre-miRNA. If two mature miRNAs are excised from both arms of a same pre-miRNA, for instance "hsa-miR-17", they are then called "hsa-miR-17-5p" (from the 5' arm) and "hsa-miR-17-3p" (from the 3' arm). The nomenclature for virus and plant miRNAs is slightly different: (i) for viruses, the genes are named according to the locus where the miRNA originates (for instance, "ebv-mir-BART1" is the miRNA from the virus Epstein Barr deriving from the BART locus); (ii) for plants, the names are in the form "ath-MIR166a", where "ath" is the plant species, "MIR166" is the name of the miRNA, and the suffix composed of one letter stands for the different loci that express the related mature miRNAs [START_REF] Griffiths-Jones | mirbase: microRNA sequences, targets and gene nomenclature[END_REF], 2008).

For each entry, several features concerning the pre-miRNA(s) and the respective mature miRNA(s) are provided. We present in what follows a few features that are worth highlighting. When it is available, the alignment of deep sequencing reads is given, showing the regions in the stem-loop that were mostly expressed. Frequently these regions correspond to the miRNA(s) loci. This is important because it proves that the miRNA was indeed transcribed and reached its mature stage. The clustered miRNAs are another feature that shows miRNAs that are close in location to the current entry, that is, < 10kb away from the current miRNA. This is relevant because it allows the user to identify miRNA genes that can be related to each other, since they are very probably being co-expressed. Finally, one functional feature is the list of predicted and validated targets. This information is pertinent since the user can go further in the analysis by knowing which genes are potentially regulated by the given miRNA.

As concerns the organisation of the data in the miRBase ftp, the information provided is separated in the following files:

• miRNA.dat: all miRNA entries in EMBL format.

• hairpin.fa: predicted pre-miRNA sequences in fasta format.

• mature.fa: mature miRNA sequences in fasta format.

• miRNA.dead: removed entries from the database.

• miRNA.diff: differences between the current and the last release.

• miFam.dat: family classification of related hairpin sequences.

In a separated directory (genomes), the gff files with the genome coordinates of the miR-NAs and pre-miRNAs are indicated.

1.2.2 Computational methods for microRNA identification [START_REF] Mendes | Current tools for the identification of miRNA genes and their targets[END_REF] classify the methods for miRNA prediction in five categories: (i) filter based approaches; (ii) machine learning approaches; (iii) target centred approaches; (iv) mixed approaches; and (v) homology search methods. In filter based approaches, the (pre-)miRNA features are verified in different filtering steps. In machine learning approaches, the methods are trained with the features of known (pre-)miRNAs to be later used for prediction. Target centred approaches are based on a more functional perspective, that is, target sequences are used to determine a potential miRNA as functional or not. Mixed approaches use high throughput experimental data and computational strategies for miRNA prediction. Finally, the homology based approach consists simply in searching for homologous miRNA sequences and/or structures, in the great majority of the cases using alignment methods. It is worth observing that many of the approaches actually use homology in one way or another (e.g. to verify sequence conservation).

In general, the methods implement characteristics originating from the biogenesis process of a miRNA to determine if a sequence is functional. The features are mainly related to the pre-miRNA hairpin, such as free energy, length of the stem-arms and terminal loop, percentage of paired nucleotides within the miRNA duplex; if small RNA sequencing data is used, the pattern of the read stacks is verified to be consistent or not with the one of an expressed miRNA. We thus provide an overview of the current methods for miRNA prediction by describing how the different methods perform this task. The described methods are summarised in Table 1.1 together with the categories to which they mainly belong.

MiRscan was one of the first methods developed for the prediction of miRNAs that is still available. It uses seven features to characterise a miRNA, such as the number of base pairs involving the miRNA candidate, conservation between related species, bulge symmetry between the two species, etc. For each of these features, the authors compute a log-odd score, and then sum them up to obtain an overall score that represents the miRNA. In that time, the authors detected 30 new genes in Caenorhabditis elegans [START_REF] Lim | The microRNAs of caenorhabditis elegans[END_REF].

Triplet-SVM implements a support vector machine (SVM) classifier trained with the features extracted from every 3 adjacent nucleotides within the hairpin structure, the authors call it triplet element features. Instead of computing the hairpin structure during the execution of the method, it requires it a priori, as an input to the software. The authors trained and tested their method on a human dataset [START_REF] Xue | Classification of real and pseudo microRNA precursors using local structure-sequence features and support vector machine[END_REF].

ProMiR is a probabilistic co-learning method based on a paired hidden HMM implemented for the prediction of miRNA with either close or distant homologs. It incorporates both sequence and structural information in a probabilistic framework and also checks for the presence of signals, such as 3' overhang, left by Drosha [START_REF] Nam | Promir ii: a web server for the probabilistic prediction of clustered, nonclustered, conserved and nonconserved microRNAs[END_REF]. As dataset, the authors used miRNAs from human chromosomes 16, 17, 18 and 19 [START_REF] Nam | Human microRNA prediction through a probabilistic co-learning model of sequence and structure[END_REF].

miRAlign identifies novel miRNAs based on sequence and structure alignment. It differentiates itself from other homology search methods because it is able to identify distant homologs, assuming little conservation of the mature miRNA. Moreover, it considers more properties of the miRNA structural conservation for the prediction of new candidates. The method was applied to Anopheles gambiae and 59 new miRNA genes were detected [START_REF] Wang | MicroRNA identification based on sequence and structure alignment[END_REF].

RNAmicro is the implementation of a SVM classifier that evaluates the information of a multiple sequence alignment. To identify the miRNAs the authors use a sliding window approach to extract segments from the genome. Then for each segment, the consensus sequence and structure are computed, and an automaton is used to evaluate the consensus secondary structure. The alignments that do not respect a few criteria are eliminated and the remaining ones are used to build the feature vector for the SVM classifier. The authors applied their method to the genomes of mammals, urochordates, and nematodes [START_REF] Hertel | Hairpins in a haystack: recognizing microRNA precursors in comparative genomics data[END_REF].

miRFinder scans whole-genomes for hairpin candidates and, by means of a SVM, evaluates the robustness of these candidates based on 18 parameters, including the minimum free energy (MFE), the frequency of the different kinds of motifs inside the hairpin structure, base pairing of the mature miRNA, etc. The search is performed in a pairwise manner, meaning that the user should provide a closely related genome in addition to the query genome. The authors applied their method to the genome pairs of chicken/human, and to Drosophila melanogaster /Drosophila pseudoobscura [START_REF] Huang | Mirfinder: an improved approach and software implementation for genome-wide fast mi-croRNA precursor scans[END_REF].

MiPred is the implementation of a random forest method that uses a hybrid feature by incorporating the local contiguous structure-sequence composition, the MFE of the secondary structure, and the P-value of a randomization test. For training and testing their method, the authors used human pre-miRNA data; real pre-miRNAs were obtained from miRBase (at that time, called the miRNA Registry database), and the pseudo pre-miRNAs were the same as those used by the authors of Triplet-SVM [START_REF] Jiang | Mipred: classification of real and pseudo microRNA precursors using random forest prediction model with combined features[END_REF].

miRank is based on a random walk ranking algorithm to characterise novel miRNAs from genomes with just a few annotated miRNAs. Differently from other machine learning approaches, this model can generalise with just a few samples. The method requires positive miRNA samples for the training step, but no negative samples are necessary. For training and validation, the authors used the genome of Anopheles gambiae, which is the vector of malaria [START_REF] Xu | MicroRNA prediction with a novel ranking algorithm based on random walks[END_REF].

SSCprofiler is a probabilistic method based on profile hidden Markov models to predict novel miRNAs. The model is trained over a set of features arising from the sequence, structure, and conservation of known miRNAs. The authors trained the model with human pre-miRNAs and applied it to cancer-associated genomic regions in search of novel miRNAs [START_REF] Oulas | Prediction of novel microRNA genes in cancer-associated genomic regions—a combined computational and experimental approach[END_REF].

HHMMIR is a de novo predictor based on a hierarchical hidden Markov model that does not require evolutionary conservation. To predict the miRNAs, the authors set a template for the structure of a typical pre-miRNA hairpin from publicly available data. They then build the HHMM model over this template that is comprised by the following regions: terminal loop, extension (area between the terminal loop and the miRNA duplex), the miRNA duplex itself, and the pri-miRNA extension. The model was trained over a human dataset and was tested on mammals, birds, fishes, worms, flies and plants [START_REF] Kadri | Hhmmir: efficient de novo prediction of mi-croRNAs using hierarchical hidden markov models[END_REF].

MIReNA finds miRNAs, given a genome and a set of known miRNAs, using a filterbased approach with no learning at a genomic scale. It uses five (physical and combinatorial) conditions to define an acceptable pre-miRNA: the miRNA cannot fold itself into a hairpin structure, there is a strong pairing between miRNA and miRNA*, the percentage of unmatched nucleotides within the hairpin, and the MFE and MFE indices are below a certain threshold. Additionally, the authors use a RepeatMasker filter, an EST data filter, and another filter that eliminates other types of RNAs. The option of using deep sequencing data is also available. To compare and validate their method, the authors used, besides the human genome, six other eukaryotic species, including Caenorhabditis elegans and Arabidopsis thaliana [START_REF] Mathelier | Mirena: finding microRNAs with high accuracy and no learning at genome scale and from deep sequencing data[END_REF].

CSHMM uses a Context-Sensitive Hidden Markov Model to represent pre-miRNA structures with estimated transition probabilities. Initially, it uses a file with the secondary structure of human pre-miRNAs to set its parameters. It is then trained with the sequences of the same positive human pre-miRNAs and with a set of negative or pseudo pre-miRNAs. Once the model is set, the authors compute the most likely sequence and its likelihood [START_REF] Agarwal | Prediction of novel precursor miRNAs using a context-sensitive hidden markov model (cshmm)[END_REF].

miRD is a webserver which runs an implementation of two independent SVM models based on two different sets of features. To combine these two models, a boosting method was used. In practice, miRD has two applications: (i) to compute the probability of a candidate pre-microRNA to be a real one; and (ii) to extract the probable pre-microRNAs from deep sequencing data. The authors predicted 92 novel pre-miRNA candidates from a small RNA sequencing dataset of the human fetal ovary [START_REF] Zhang | Prediction of novel pre-microRNAs with high accuracy through boosting and svm[END_REF].

MiRPara is an SVM implementation trained with sequences from miRBase. It makes available a script to generate the model according to the miRBase release and to the desired organism(s). The authors used a set of 77 features as input to the SVM classifier; these features were based on characteristics of the miRNA, pre-miRNA and pri-miRNA which are important to the biogenesis of a miRNA [START_REF] Wu | Mirpara: a svm-based software tool for prediction of most probable microRNA coding regions in genome scale sequences[END_REF].

miR-BAG is a set of three complementarity approaches (naive Bayes, Best First Decision tree and SVM) which employs different miRNA features such as matrices with specific miRNA guided structural profile and structural triplet density variation profiles with respect to the position of the miRNA. The prediction can be performed at both genomic scale or by using deep sequencing data. The genomes of six species human, mouse, rat, dog, nematode, and fruit fly were used by the authors. [START_REF] Jha | mir-bag: bagging based identification of microRNA precursors[END_REF].

miRDeep is a package for the discovery of miRNAs from deep sequencing data. It first eliminates reads which map to many loci in the genome, and optionally it can remove reads mapping to rRNAs, tRNAs, etc. To obtain potential pre-miRNAs, the authors use the information of the mapped reads against the genome. Pre-miRNAs with an unlikely structure are discarded and the core algorithm computes a probabilistic score related to the structure and signature of the pre-miRNA candidate. To validate their method, miRNAs from Caenorhabditis elegans and Homo sapiens were used [START_REF] Friedländer | mirdeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades[END_REF].

miRNAFold is an ab-initio method that, given a sequence as input, directly searches for pre-miRNA hairpins. The main idea of the algorithm is to find a long stem, which is then taken as an anchor to predict the hairpin structure, attempting to improve the search time. To test their method, the authors used chromosome 19 of the human genome, chromosome 2 of mouse, chromosome 4 of the zebrafish, and chromosome 7q of the sea squirt [START_REF] Tempel | A fast ab-initio method for predicting miRNA precursors in genomes[END_REF].

RNA secondary structure prediction

Since the great majority of the methods for miRNA prediction use in one way or another the secondary structure of the hairpin pre-miRNA, we focus in this section on this issue.

RNA folding consists in intra-strand base-pairing to produce a secondary structure. As concerns RNA, guanine and cytosine (GC) pair by forming a triple hydrogen bond, adenine and uracil (AU) pair by a double hydrogen bond; additionally, guanine and uracil (GU) can pair by forming a single hydrogen bond. The stability of a given secondary structure depends on: (i) the number of GC versus AU and GU base pairs (the higher the energy bonds, the more stable the structures are, eg, GC is more stabilising than AU); (ii) the number of base pairs in a stem region (longer stems result in more bonds); (iii) the number of bases in a hairpin loop region (the formation of loops with more than 10 or less than 5 bases requires more energy); and (iv) the number of unpaired bases within the structure, either interior loops or bulges (unpaired bases decrease the stability of the structure) [START_REF] Gesteland | The RNA world[END_REF]Mathews et Turner, 2006) . [START_REF] Mendes | Current tools for the identification of miRNA genes and their targets[END_REF].

Methodological background

The stability of a secondary structure is quantified as the amount of free energy released or used by forming base pairs. Positive free energy requires work to form a configuration; negative free energies release stored work. Therefore, the more negative the free energy of a structure, the more likely is the formation of that structure because more stored energy is released. Free energy changes of coupled reactions are additive, so one can determine the total free energy of a secondary structure by adding all the component free energies associated to each two consecutive base pairs (units are kilocalories per mole, kcal/mol). This is used to predict the secondary structure of a given sequence. Finding a base pair configuration with the minimum possible free energy is the aim of most secondary structure prediction algorithms [START_REF] Nelson | Dna and RNA oligomer thermodynamics: The effect of mismatched bases on double-helix stability[END_REF].

To compute the minimum free energy of a sequence, empirical energy parameters are used. These parameters summarise the free energy change (positive or negative) associated to all possible pairing configurations [START_REF] Turner | Nndb: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF]. The energy parameter depends on the place (motif) in the structure where the bases are located. An RNA secondary structure can have the following motifs: (i) a helix is the stacking of canonical base pairs (GC, AU and GU); (ii) a loop is a set of non-canonical pairs (that is, unpaired nucleotides): a terminal loop has one appended helix; an internal loop has two appended helices; a bulge loop is similar to an internal one, however, the non-canonical pairs appear just in one strand of the loop; a multibranch loop (junction) is a loop with at least three appended helices; and an exterior loop is a series of adjacent unpaired bases which are not accessible by any base pair; and (iii) a dangling end is the stacking of nucleotides at the end of helices [START_REF] Turner | Nndb: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF][START_REF] Serra | Improved parameters for the prediction of RNA hairpin stability[END_REF]. All the previously described motifs are presented in Figure 1.6a. A structure can also have a more complex motif called pseudoknot, which is formed by at least two hairpin structures, in which half of one of the stems is intercalated with the other hairpin, a hairpin being a structure formed by a stem and a terminal loop (see Figure 1 [START_REF] Akutsu | Dynamic programming algorithms for RNA secondary structure prediction with pseudoknots[END_REF]. The problem of predicting an RNA secondary structure is defined as follows. If S = s 1 s 2 ...s m is an RNA sequence of length m, then the secondary structure of S is defined as a set R of base pairs that satisfies the following criteria: (i) if s i base pairs with s j , then i < j; (ii) a base pair can only be established if the two bases are, at least, 3 nucleotides apart from one another; (iii) s i can base pair with one, and only one other base s j . Pseudoknots are usually not permitted because of the complexity that this leads to. The goal is thus to maximise the number of base pairs within R, or minimise the energy associated to the set of base pairs.

One of the first methods for RNA secondary structure prediction was described by [START_REF] Nussinov | Algorithms for loop matchings[END_REF]. The algorithm proposes the maximisation of the number of base pairs to find the best structure. For each position i in the sequence, one should verify all the possible cases: (a) i, j base pair; (b) i is unpaired; (c) j is unpaired; (d) i, j base pair with, respectively, k and k + 1. The recurrence for this algorithm is presented in Equation 1.1 [START_REF] Eddy | How do RNA folding algorithms work[END_REF]):

E(i, j) = max ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ E(i + 1, j -1) if i and j base pair E(i + 1, j) E(i, j -1) max i<k<j [E(i, k) + E(k + 1, j)] (1.1)
Clearly, filling each cell in the DP matrix takes O(n) time, and since there are O(n 2 ) cells, the complexity for the whole procedure is in O(n 3 ). However, maximising the number of base pairs is a naïve approach; a more realistic one minimises the free energy of the structure, as proposed for example in [START_REF] Mathews | Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure[END_REF]. The recurrence for the latter algorithm is presented in Equation 1.2:

E(i, j) = min ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ E(i + 1, j) E(i, j -1) min i<k<j [E(i, k) + E(k + 1, j)] P (i, j) if i and j base pair (1.2)
To minimise the free energy, one more table P is required to store the different types of motifs a structure can have, although the complexity in the worst case remains the same, namely in O(n 3 ) [START_REF] Mathews | Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure[END_REF](Mathews et al., , 2006)).

A minimum energy folding algorithm will return only one secondary structure, though there are many candidates for the natural structure. To address this problem, some algorithms (such as Zuker's mfold) are designed to provide a set of suboptimal solutions. Inferring what structure is truly representative of the natural structure requires additional information. Phylogenetic information is often used to constrain the search by identifying highly conserved motifs. Some programs allow the user to specify constraints on the secondary structure, by specifying paired, single-stranded, or non-pairable regions, or [START_REF] Gesteland | The RNA world[END_REF].

Suboptimal folding One sequence can have several different secondary structures with very similar free energies, which can also be quite close to the minimum. Instead of providing a single optimum structure, a suboptimal approach provides all the partial structures which can be later refined to complete structures. This is done during the traceback step, in which suboptimal structures are chosen. Algorithms implementing this strategy were described in [START_REF] Williams | A dynamic programming algorithm for finding alternative RNA secondary structure[END_REF] and [START_REF] Wuchty | Complete suboptimal folding of RNA and the stability of secondary structures[END_REF]. [START_REF] Mccaskill | The equilibrium partition function and base pair binding probabilities for RNA secondary structure[END_REF] aggregated more quality and robustness for the folding by using a partition function in the prediction of RNA secondary structures. A partition function considers the statistical properties of a system, in this case a secondary structure, in relation to thermodynamics. The Boltzmann factor is defined by e -ΔG • /RT , where ΔG • is the free energy of the structure, R is the constant of gas, and T is the temperature given in kelvin. Then, the probability of a given structure is defined by its Boltzmann factor divided by the partition function Z, which is defined by the sum of all the Boltzmann factors.

Partition function

Obviously, there are a number of limiting assumptions to existing folding algorithms. These include the kinetics of folding during transcription, the difficulty in predicting pseudoknots, the role of chaperone proteins in folding, and the importance of modified bases (e.g. methylated bases). Some algorithms attempt to incorporate these considerations (e.g., [START_REF] Rivas | A dynamic programming algorithm for RNA structure prediction including pseudoknots[END_REF] and [START_REF] Ruan | An iterated loop matching approach to the prediction of RNA secondary structures with pseudoknots[END_REF] for pseudoknots).

Algorithm for global sequence alignment (Needleman-Wunsch) We also present here the Needleman-Wunsch [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF] algorithm for global sequence alignment, since our method for the prediction of miRNAs makes use of it to approximate the free energy of a hairpin structure. Global alignments are applied to sequences of similar length, for which the algorithm will try to align every nucleotide in the sequences. The base case and recurrence for the algorithm are presented in Equations 1.3 and 1.4:

W (i, 0) = W (0, j) = 0, i, j ∈ 0..n (1.3) W (i, j) = max ⎧ ⎨ ⎩ W (i -1, j -1) + f (s i , s j ) W (i, j -1) + γ W (i -1, j) + γ (1.4)
where n is the length of the aligned sequences, f (s i , s j ) is the function returning the score or penalty for, respectively, a match or a mismatch, and γ is the penalty for a gap. Using this recurrence one should take, in the worst case, O(n 2 ) time to align two sequences of length n [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF]. Algorithm 1 contains this forward-filling step of the algorithm for two sequences A and B of length m and n, respectively.

Algorithme 1: Forward step of Needleman-Wunsch's algorithm. Data : Two sequences A and B Result : Dynamic matrix fulfilled

1 for i ← 0 to m do 2 F (i, 0) ← 0 3 for i ← 0 to n do 4 F (0, j) ← 0 5 for i ← 1 to m do 6 for j ← 1 to n do 7 match ← F(i-1,j-1) + f(A i , B j ) 8 delete ← F(i-1, j) + γ 9 insert ← F(i, j-1) + γ 10 F(i,j) ← max(match, insert, delete)
Once the dynamic programming matrix F is filled up, the next task consists in recovering the alignment by backtracking along the matrix. The recovery is performed by means of a recursion starting in cell F (m, n) and ending when the left or the top part of the matrix is reached, as shown in Algorithm 2.

Algorithme 2: Backtracking step of Needleman-Wunsch's algorithm. Data : Dynamic programming matrix F Result : Alignment and its score 1 AlignmentA ← "" 2 AlignmentB ← ""

3 i ← m 4 j ← n 5 while i > 0 or j > 0 do 6 if i > 0 and j > 0 and F (i, j) == F (i -1, j -1) + f (A i , B j ) then 7 AlignmentA ← A i + AlignmentA 8 AlignmentB ← B j + AlignmentB 9 i ← i -1 10 j ← j -1 11 else if i > 0 and F (i, j) == F (i -1, j) + γ then 12 AlignmentA ← A i + AlignmentA 13 AlignmentB ← "-" + AlignmentB 14 i ← i -1 15 else if j > 0 and F (i, j) == F (i, j -1) + γ then 16 AlignmentA ← "-" + AlignmentA 17 AlignmentB ← B j + AlignmentB 18 j ← j -1 19

Nearest neighbour energy model

To model the free energy change for the folding of RNAs, one can use the thermodynamic Nearest-Neighbour (NN) energy associated to each type of motif in the structure. By summing up the energy increment of each motif, it is possible to obtain a reasonable approximation of the free energy change for folding an RNA or, in other words, to obtain a measure of the stability of an RNA molecule [START_REF] Mathews | Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure[END_REF](Mathews et al., , 2006)).

The motifs forming an RNA structure are determined by the base-pairs AU, GC and GU. The arrangement of these base pairs can shape into the different types of motifs, such as helices, bulge loops, and internal loops. The stabilising motifs are: the Watson-Crick helix represented by the stacking of at least two base pairs; and a dangling end which is a single base at the end of a helix. The destabilising motifs are of three types: the hairpin loop which is composed of non-canonical base pairs closed by one canonical base pair; the bulge loop which is an arrangement of unpaired nucleotides in one of the strands of a helix; and finally, the internal loop which includes unpaired nucleotides in both strands of a helix. There exist three more types of motifs which are the multi-branch loop, the exterior loop, and pseudoknots. However, they are not present in a pre-miRNA stem-loop structure, and will therefore not be explored in detail here (Mathews et al., 2006;[START_REF] Turner | Nndb: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF].

As mentioned before, to compute the free energy of an RNA structure, it is necessary to sum the increments according to the type of the motif. The equations presented hereafter describe how to compute the free energy associated to each kind of motif.

The energy of a dangling end depends only on the base-pair before the dangling nucleotide

Chapter 1. Background and on this latter. For all the other types of motifs, the equations are given below. The energy of an internal loop is computed by means of Equation 1.5:

ΔG Internal = ΔG i (n) + (ΔG a * |n 1 -n 2 |) + ΔG m1 +ΔG m2 + (ΔG ru * λ) (1.5)
where ΔG i (n) is the initiation energy to form an internal loop of n ≤ 30 unpaired nucleotides; ΔG a = 0.6 is the asymmetry penalty multiplied by the absolute value of the difference between the number of unpaired nucleotides in each strand; ΔG m1 and ΔG m2 are the energy of the first and the last mismatches in the internal loop; and ΔG ru = 0.7 is the penalty for an RU closure, where R = {A, G} and λ is the lambda function which returns 0 or 1, corresponding, respectively, to the presence or absence of, in this case, the RU closure.

For the bulge loops, one should use Equation 1.6:

ΔG Bulge(n=1) = ΔG i (1) + ΔG C + ΔG s -RT ln(t) + (ΔG ru * λ) ΔG Bulge(n>1) = ΔG i (n) (1.6)
where ΔG i (n) is the energy required to form a bulge with n ≤ 30 unpaired nucleotides; if the bulge is comprised of the nucleotide C only, and there is at least one more C not in the bulge (meaning, it is paired with a G), one should add the C bulge penalty ΔG C = -0.9 kcal/mol; ΔG s is the base pair stacking around the bulge; t is the number of possible loop conformations with identical sequence; R = 8.3144621 J/mol K is the gas constant and T = 310.15 K is the temperature in kelvin. Notice that for bulges and helices, ΔG ru = 0.45 and is referred to as the penalty for a RU end (and not closure as for internal loops). For bulges and internal loops larger than 30 nucleotides (n > 30), Equation 1.7 should be applied instead:

ΔG n>30 = ΔG i (30) + 1.75 × RT × ln(n/30) (1.7)
Finally, for a helix, one should apply Equation 1.8:

ΔG helix = ΔG stck + ΔG sym + (ΔG ru * λ) (1.8)
where ΔG stck is the stacking energy of each two consecutive base pairs; ΔG sym is the symmetry correction for self complementarity duplexes; and ΔG ru = 0.45 is, as mentioned before, the RU end penalty.

All the thermodynamic NN energies used in this work, as well as the equations described above, were obtained in the Nearest Neighbor Database (NNDB) [START_REF] Turner | Nndb: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF][START_REF] Zuker | Algorithms and thermodynamics for RNA secondary structure prediction: A practical guide[END_REF].

Experimental methods for microRNA detection and quantification

As mentioned in Section 1.1.1, a sequence must fulfill three criteria to characterise a miRNA: (i) the mature miRNA must be expressed as a transcript of ∼22nt (the expressed transcript is detected by means of experimental techniques such as northern blot, small RNA sequencing, etc.); (ii) the mature miRNA must derive from a precursor miRNA with a typical hairpin structure containing small bulges and internal loops; and (iii) the pre-miRNA should be processed by Dicer, as an increased accumulation of the precursor is noticed when Dicer is absent. The experimental methods for discovering miRNAs are based on these definitions with variations among them [START_REF] Berezikov | Approaches to microRNA discovery[END_REF]. A brief description of these methods, such as PCR based methods, microarray, northern blot, and RNA sequencing is given in what follows. It is important to observe that the first three methods-PCR, microarray, and northern blot-require the miRNA sequence, that is, these methods validate specific sequences known a priori while for RNA sequencing it is not necessary to know the miRNA sequence [START_REF] Chaudhuri | MicroRNA detection and target prediction: integration of computational and experimental approaches[END_REF].

PCR based methods

The reverse transcription polymerase chain reaction (RT-PCR) protocol is based on the reverse transcription of the small RNA to cDNA. The reverse transcription is generated either by the addition of a poly A queue or by means of a stem-loop primer. The reverse transcript cDNA is then submitted to a PCR for amplification and quantification. The accumulation of the reaction product can be monitored in real time (it is then called real time RT-PCR).

The method is widely used because of its ease of incorporation; consequently, it is a very well established method with a good sensitivity and specificity. The disadvantages are the medium-throughput concerning the number of samples processed per day, and the inability to detect novel miRNAs [START_REF] Pritchard | MicroRNA profiling: approaches and considerations[END_REF][START_REF] Aldridge | Introduction to miRNA profiling technologies and cross-platform comparison[END_REF]. The available assay/platforms for RT-PCR are: TaqMan by ABI, miRCURY LNA qPCR by Exiqon, Biomark HD system by Fluidigm, SmartChip human microRNA by Wafergen, and miScript miRNA PCR array by SABiosciences/ Qiagen.

MicroRNA microarray

A microarray is a chip composed of several microscopic spots. Each spot is filled with DNA/RNA molecules for the measurement of their expression level. These molecules are specific oligonucleotide sequences (e.g., miRNAs, genes, etc.) known as probes. These probes will then be submitted to hybridisation with specific cDNA/cRNA targets under specific conditions. The occurrence of hybridisation between probe and target will be detected by a label that is linked to the target and will emit a fluorescent light. The different light spectra will quantify the level of expression [START_REF] Yin | Profiling microRNA expression with microarrays[END_REF][START_REF] Pritchard | MicroRNA profiling: approaches and considerations[END_REF].

The advantage of microarray experiments is that it is not an expensive method, while allowing for the parallel profiling of a large number of molecules. On the other hand, microarray technology has low specificity if the miRNAs are similar, and the absolute quantification of miRNA expression is not easily performed, while it is better suited for detecting the relative abundance of specific miRNAs in 2 different states [START_REF] Pritchard | MicroRNA profiling: approaches and considerations[END_REF][START_REF] Liu | MicroRNA expression profiling using microarrays[END_REF].

The available platforms for microarrays are: Geniom Biochip miRNA by CBC, GeneChip miRNA array by Affymetrix, GenoExplorer by Genosensor, MicroRNA microarray by Agilent, miRCURY LNA microRNA array by Exiqon, NCode miRNA array by Invitrogen, nCounter (not a microarray but hybridization-based) by Nanostring, OneArray by Phalanx Biotech, and Sentrix array matrix and BeadChips by Illumina.

Northern blot

Northern blot is a technique to identify a specific RNA from a bunch of RNAs. The total RNA is denatured and after the addition of an agent, the RNA remains unfolded in its linear conformation. The collection of RNAs is then sorted by size by means of an electrophoresis gel and moved to a nitrocellulose filter in which the RNAs are attached. A labelled probe is added to the filter, and it is finally submitted to autoradiography that will quantify the expression level of the given RNA [START_REF] Lodish | Molecular cell biology[END_REF].

The majority of the methods are specialised in the detection of smaller mature miRNAs. The advantage of northern blot is that it can detect a wide range of sizes from primary miRNAs to mature ones. However, the approach is low throughput and has low sensitivity, besides being time consuming and requiring a large amount of total RNA [START_REF] Pritchard | MicroRNA profiling: approaches and considerations[END_REF][START_REF] Liu | MicroRNA expression profiling using microarrays[END_REF].

RNA sequencing

Since we produced and analysed small RNA sequencing data during this thesis, we give more details about this experimental technique, focusing on the Ilumina platform. Details about the sRNAseq data produced and analysed are provided in Chapter 3.

RNA sequencing makes use of NGS technology to verify the presence of and to quantify the RNAs from a genome in a specific condition. For miRNA sequencing, the input library is enriched for this kind of molecule. The advantages of using NGS for miRNA profiling are identification of known and novel miRNAs and precision in identifying very similar miRNAs, such as isomiRs of different length and miRNAs differing by a single nucleotide. On the other hand, small RNA sequencing can produce several putative small RNAs of novel sequence, and they are not necessarily bona fide miRNAs [START_REF] Pritchard | MicroRNA profiling: approaches and considerations[END_REF].

The first step of the experiment consists in the isolation of the total RNA from a sample by means of a specific reagent, which depends on the kit used. The total RNA is then filtered by size using a polyacrylamide gel which is submitted to a process of electrophoresis. In the case of miRNAs, the selected size range is from 17nt to 25nt. Adapters are thus ligated on both 5' and 3' ends of the RNA sequence to act as the binding sites for the primers used during the next step (RT-PCR). Since sequencing technology is designed to sequence only DNA, it is necessary to convert the RNA into cDNA by means of a reverse transcriptase. The total amount of cDNA is amplified with a PCR, and it is finally submitted for sequencing. Figure 1.7 presents a summary of the above-mentioned process.

For sequencing, different kinds of platforms exist, we mention the most used ones: sequencingby-synthesis on the Illumina1 platform, pyrosequencing on the 454 Life Sciences2 platform, and ABI Solid Sequencing3 platform. Since for this thesis, the miRNAs of the Acyrthosiphon pisum (pea aphid) were sequenced on an Illumina platform, we provide more details about it.

The cDNA library is given as an input to the Illumina sequencer. The first step occurs in a device called Cluster Station over a 8-channel flow cell where amplification of the reads occurs. Oligos, that are complementary to the adaptors ligated to the cDNAs, are attached to this flow cell. The cDNA fragments will thus bind to these oligos and the DNA polymerase will produce approximately one million copies of the original fragments, which is a sufficient amount to generate the required signal intensity of the incorporated bases. After that, the single cDNA fragments are replaced by clusters of fragments. The four nucleotides, enriched with a unique fluorescent label, are then added to the channels of the flow cell, together with the DNA polymerase, to be incorporated into the clustered fragments. Each base incorporation is followed by an imaging step that scans the emitted light associated to each base. Each base incorporation corresponds to a cycle; consequently, the number of cycles is equivalent to the length of the fragments.

./figures/figure_ilumina_sequencing.pdf Figure 1.7: Preparation of a miRNA Illumina sequencing library. The steps are described in order as follows: (i) isolation of total RNA from the sample; (ii) size fractionating of total RNA using denaturing PAGE, and selection of small RNA by size (17-25 nt); (iii) 3' and 5' adapter ligation; (iv)

reverse transcription of RNA sequences, and PCR amplification; (v) flow cell attachment, bridge amplification, annealing of sequencing primers and base extension, base calling till the number of cycles is finished. Figure taken from Wikipedia.

As mentioned before, a plethora of putative novel miRNAs is produced, and they are not necessarily bona fide miRNAs. In order to retain the real miRNAs, further criteria must be applied for the annotation of small RNA sequences as a miRNA: length of ∼22nt, hairpin structure of the corresponding precursor miRNA, sequenced reads aligning to both arms, -3p and -5p of the precursor, and, when a close species is available, conservation across species. All the details concerning the sRNAseq data analysis of the pea aphid is provided in Chapter 3.

Computational methods for target prediction

As briefly mentioned before in Section 1.1.3, the exact mechanism used by miRNAs to regulate target gene expression is still uncertain. Cases have been reported of target mRNA cleavage, translation repression, and also of activation of gene expression. There are even evidences, in both plants and animals, that miRNAs can reduce protein (and not mRNA) levels. More specific mechanisms are not clear, the decrease of gene expression can be associated to the prevention of translation initiation or elongation, and also to the proteolysis of peptides [START_REF] Liu | Identifying miRNAs, targets and functions[END_REF][START_REF] Pasquinelli | MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship[END_REF].

Although the mechanism of miRNA target regulation is not clear, the problem is modelled focusing either on characteristics that are specific of some cases (such as cleavage of mRNA), or on features that are in principle common to all the cases.

Most of the computational methods incorporate features related to the base pairing interaction between miRNA and target. These include the presence of perfect complementarity of the seed region located at the 5' end of the miRNA (nucleotides 2-7), and 3' UTR for the target mRNA. Accessibility of the 3' UTR in its secondary structure is also verified and it is associated to AU content in the flanking regions. Target conservation is used to eliminate false positives. Even if these general rules have been successful in many predictions, a substantial part of the methods diverge in their results, with levels of false predictions that are not easy to evaluate. One of the causes of the previous mentioned problem is the lack of experimentally validated miRNA-mRNA interactions [START_REF] Witkos | Practical aspects of microRNA target prediction[END_REF][START_REF] Pasquinelli | MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship[END_REF].

As mentioned before, methods for target prediction basically employ two features: (i) base pairing between miRNA and target, specially considering the seed region; and (ii) target conservation across related species. The currently available methods are TargetScanS, miRanda, DIANA-microT, PicTar, and RNAHybrid.

TargetScanS requires a seed region of length 6nt from positions 2 to 7 in relation to the miRNA. It also demanded target site conservation across all the five genomes the authors studied: human, mouse, rat, dog, and chicken. The authors observed the presence of conserved adenosines flanking the seed region in the target mRNA, suggesting that these nucleotides can play a decisive role in the recognition of miRNA targets [START_REF] Lewis | Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets[END_REF].

miRanda looks for sequence match between miRNA and target, permitting GU wobble pairs and moderate insertions and deletions, while giving a stronger weight to complementarity at the 5' end of the miRNA. The free energies of the duplexes are then computed with the Vienna package [START_REF] Lorenz | ViennaRNA package 2.0[END_REF]. Conservation is thus verified according to three factors: (i) a miRNA should match orthologous UTRs in the three species the authors studied (i.e., Drosophila melanogaster, Drosophila pseudoobscura and Anopheles gambiae); (ii) the target sequences in all three species should respect a threshold identity with each other; and (iii) the positions of both target sites are equivalent according to a cross-species UTR alignment [START_REF] Enright | MicroRNA targets in drosophila[END_REF].

DIANA-microT requires the interactions to meet two criteria. The first is high-affinity measured on the basis of free energy. The second criterion considers the proteins associated to the interaction between miRNA and target; the authors verify it by analysing the position and sizes of the loops (bulges and internal loops) within the miRNA:target duplex that are imposed by the associated proteins [START_REF] Kiriakidou | A combined computational-experimental approach predicts human microRNA targets[END_REF].

PicTar defines a seed as a perfect base pair of length 7nt starting at position 1 or 2 at the miRNA 5' end. Insertions or mutations are allowed only if the free energy of the duplex does not increase and does not form GU wobbles. A combined score is computed for the mRNA target; it is comprised of the maximum likelihood of the given target to be regulated by a set of miRNAs, plus a few other features observed in the results of experimental interactions. If a miRNA seed aligns to overlapping positions of the UTR sequences across the different species, conservation is considered to be verified [START_REF] Krek | Combinatorial microRNA target predictions[END_REF].

RNAHybrid is an adaptation of the classical RNA secondary structure prediction described by [START_REF] Zuker | Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information[END_REF]. Instead of folding a single sequence in the energetically most favourable conformation, it determines the most favourable hybridisation site between two sequences. The presence of a seed is also verified by this method, however the properties (such as length, position) are set by the user [START_REF] Krüger | RNAhybrid: microRNA target prediction easy, fast and flexible[END_REF].

PITA predicts target sites by verifying the presence of seed regions (allowing for GU wobbles and mismatches). It next uses target accessibility, the core of their algorithm, a concept that is strictly related to the secondary structure of the target transcript. The hypothesis is based on the fact that the mRNA structure is an important factor in the recognition of the target, by thermodynamically favouring or not the interaction. The free energy gained from the formation of the miRNA-target duplex, and the energetic cost of unfolding the target to make it accessible to the miRNA are computed, and the tradeoff between these two measures is assessed to classify an interaction as functional or not [START_REF] Kertesz | The role of site accessibility in microRNA target recognition[END_REF].

Experimental methods for microRNA target identification

This section is included for the sake of completeness, since we did not performed any wet experiments involving target identification. It is then a concise section providing only an overview of the methods available.

Interactions between miRNA and target are often validated by fusing the target site to a reporter gene and verifying, in the presence or absence of the miRNA, if regulation occurs. In this case, the original cellular context is lost. Nevertheless, a recent technology called crosslinking immunoprecipitation sequencing (CLIP-seq) allows the identification of endogenous target sites by means of the sequencing of those targets that co-immunoprecipitate with RISC components [START_REF] Pasquinelli | MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship[END_REF].

The use of microarray experiments, after miRNA overexpression or knockdown, consists in another method for identifying genes regulated by miRNAs. Since miRNAs reduce the levels of gene transcripts, measuring the expression of a given mRNA after an abnormal miRNA expression provides an effective manner to verify functional interactions [START_REF] Thomas | Desperately seeking microRNA targets[END_REF].

Using a stable isotope labeling with amino acids (SILAC) in cell culture followed by mass spectrometry based proteomics, one can evaluate the effect of down or overexpression of miRNAs on global protein expression [START_REF] Thomas | Desperately seeking microRNA targets[END_REF].

There are mainly two types of regulation performed by miRNAs: translation repression and mRNA cleavage. Parallel analysis of RNA ends (PARE), also known as degradome sequencing, detects interactions originating from the second case. It detects the products of mRNA cleavage. An RNA adaptor is ligated on the mRNA 3' fragments resulted from the Argonaute-mediated cleavage. These fragments are submitted to RT-PCR for enrichment, followed by deep sequencing [START_REF] German | Construction of parallel analysis of RNA ends (pare) libraries for the study of cleaved miRNA targets and the RNA degradome[END_REF].

Other methods exist, which are mainly improvements of the mentioned ones; we do not detail them here as this would be beyond the scope of this thesis. For a complete survey of the current methods, see [START_REF] Thomson | Experimental strategies for microRNA target identification[END_REF] This chapter is strongly based on the paper Higashi et al. (ress). Here, we address the problem of prediction of miRNAs, focusing on three main issues: (i) efficiency of the algorithm for free energy due to the use of a quadratic algorithm (instead of a cubic one as used so far by other methods); (ii) high quality hairpin secondary structure when sRNA-seq data is available; and (iii) dependence on as little information as it is possible to compute the free energies. These items were defined in details and implemented in a software called Mirinho, which is available at http://mirinho.gforge.inria.fr. Besides the better time complexity of the algorithm itself, a speed-up was implemented during the Master internship of a bioinformatics student, Cyril Fournier, whom I co-advised together with Marie-France Sagot.

Introduction

Given the ubiquity of miRNAs and their functional importance, it became crucial to develop methods for the prediction and analysis of miRNAs. As a consequence a plethora of such methods have been developed (as shown in Chapter 1). Despite all the effort put in developing them, there remain a number of issues that need to be addressed: (i) to predict the characteristic hairpin structure of a pre-miRNA, the vast majority of the existing softwares rely on a folding algorithm of cubic time complexity which is suitable when the input is small enough, but it can become impracticable when the size of the input increases; (ii) for longer pre-miRNAs (such as in plant), such folding methods moreover can produce hairpin structures different from the ones provided in miRBase [START_REF] Kozomara | mirbase: integrating microRNA annotation and deep-sequencing data[END_REF], which uses sRNA-seq data to do so; (iii) together with folding, most methods then rely on further information that must be learned from previously validated miRNAs of closely related genomes (at a minimum within the same clade, plant or animal) for the final prediction of new miRNAs in order either to set the parameters of the model or to restrict the search to a limited space.

We therefore developed Mirinho to address all three issues. The search for pre-miRNAs is concentrated on regions with the same length as the two stem-arms separated by the length of the terminal loop. The direct application to sRNA-seq data guarantees a better quality in the prediction of the pre-miRNA structures. A quadratic time complexity algorithm improves the practical efficiency of the free energy computation. As neither of the two attributes used (length of stem-arm and terminal loop) are species-specific within the animal or the plant kingdom (they differ only between these two kingdoms), the method can easily be applied for predicting pre-miRNAs in either clade.

Importantly, while the method we provide is thus much simpler, faster, and general to use, we also show for tested examples that it has a sensitivity and precision as good as other methods, in some cases even better. Moreover, we show that the secondary structures predicted by Mirinho are much closer to the ones available in miRBase than for the other compared methods.

Material and methods

Algorithm

Screening the genome to identify potential pre-miRNAs A pre-treatment of the data is performed in order to identify all the inputs for our algorithm. This was done without loss of information at this step, meaning that all sequences that are potential candidates for being a pre-miRNA were selected.

We set a sliding window w of length = 25, that is the mean length of a miRNA sequence. For each stem-arm represented by st 1 = [w i , ..., w i+ -1 ], we looked for its putative stem-arm pair st 2 = [w i+ +n-1 , ..., w i+2 * +n-1 ], n nucleotides away from the first one, with n between 5 and 20. We thus have that w represents the potential stem-arms and n the length of the terminal loop, as shown in Figure 2.1. Each such pair (w, n) will be an input for the alignment algorithm described next. Assessing the potential pre-miRNAs Each pair of putative stem-arms screened in the previous step was given to an alignment algorithm in order to evaluate whether it is a stable stem-loop structure. For that, we implemented the Needleman & Wunsch global alignment algorithm [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF]) (Section 1.2.2) with a scoring strategy based on the Nearest Neighbour energy model (Section 19). Instead of using the sum of the integer penalties for gaps, matches and mismatches, the alignment is assessed according to the free energy related to each two consecutive nucleotides in the alignment. We define the alphabet Σ = {M xy , S xy , I xy , D xy }, where the symbols correspond, respectively, to Match, Mismatch, Insertion and Deletion, and x, y ∈ {A, U, C, G, -}. The definition of an alignment of two putative stem-arms, st 1 and st 2 , is a vector comprised by the symbols in Σ, such that align(st

W[i] W[i+l-1] W[i+l-1+n] W[i+2*l-2+n]
1 , st 2 ) = v and v = [v i , v i+1 , ..., v n ], where v i ∈ Σ.
To determine the stability of a pre-miRNA stem-loop, we go through vector v and sum up the free energy of each pair (v i , v i+1 ) according to the type t of the motif it is inserted in. For that, we use Equation 2.1 below to compute the energy of each motif in the structure:

(t) = k(t, m) + m-1 i e(v i , v i+1 ) (2.1)
where t is the motif type that can be an internal loop, a bulge loop or a helix. The value k(t, m) accounts for the penalties associated to the motif t, which appears m times in the structure. For example, for a motif of type t = helix, one should consider the symmetry correction ΔG sym for self-complementary duplexes (see Equation 1.8). Finally, the function e returns the energy associated to the pair

(v i , v i+1 ).
We then sum all the energies related to the different types of motifs to obtain the final free energy E of the structure using Equation 2.2:

E = (t) (2.2)
where t is again the different types of motifs a given structure can have.

Alignment speed-up Considering that a stable hairpin structure should not contain very large bulges neither internal loops, an ideal alignment should be concentrated around the main diagonal of the dynamic programming (DP) matrix. Instead of using the whole matrix, the user can therefore constrain the alignment to this diagonal and prune parts of the bottom-left and top-right corners of the matrix, thus saving time in the computation of the free energies with a small loss, as shown in Figure 2.2. A parameter dw (diagonal width) is established that depends on the length of the aligned sequences and on a compromise between sensitivity and precision in relation to the version that uses the full matrix (see the Section 2.3.2 to determine how to set an appropriate value for this parameter).

Dataset

To set an appropriate energy threshold for Mirinho, we chose chromosomes from six different metazoan genomes with the respective miRBase miRNA annotations.

• Chromosome 25 from Bos taurus (27 miRNAs) To compare our method to other pre-miRNA predictors, we applied it to:

• the prediction of plant pre-miRNAs: we used the sequence of chromosome 4 of Arabidopsis thaliana (version 2.0) as well as 340,114 reads of high-throughput small RNA sequencing data from Arabidopsis thaliana (GEO accession number GPL3968).

• three animal chromosomes for which the miRNAs are well characterised:

-Chromosome III of Caenorhabditis elegans (44 miRNAs)

-Chromosome 2R of Drosophila melanogaster (92 miRNAs)

-Chromosome 19 of Homo sapiens (234 miRNAs)

In the latter case, as two of the softwares to which Mirinho was compared are too slow, the predictions were performed on smaller sets of sequences obtained in the following way: for each of the three chromosomes (III in Caenorhabditis elegans, 2R in Drosophila melanogaster, and 19 in Homo sapiens), 10 miRNAs were randomly chosen together with 100nt both up and downstream. Each fragment (miRNA+extension) of length n was flanked by sequences of the same length, which were generated based on the nucleotide distribution of the given chromosome. In the end, we obtained three different sequences of ∼4265nt that were given as input to CSHMM, MIReNA, Mirinho, and miRPara (mentioned in the Section 2.2.3).

For computing sensitivity and precision in a genomic scale, we used the genomes of three insects that are of special interest for our group:

• Acyrthosiphon pisum genome assembly version 2 (123 miRNAs)

• Culex quinquefasciatus genome assembly version 1 (120 miRNAs)

• Heliconius melpomene genome assembly version 1.1 (101 miRNAs) All the chromosomes, genomes, and sRNA-seq data were obtained from the NCBI. The annotations concerning the known (pre-)miRNAs were obtained from mirBase (release 20) [START_REF] Kozomara | mirbase: integrating microRNA annotation and deep-sequencing data[END_REF].

Compared methods

To compare the accuracy of our method with other predictors, we first made an extensive search of the available ones (see Table 1.1 in Section 1.2.2). We put aside the predictors that required other kinds of input files than just the fasta sequence and/or sRNA-seq data, as well as those incompatible with the Unix system. Web-servers were also disregarded because there always is a restriction to the length of the sequence that may be input. The methods that remained were CSHMM, MIReNA, and miRPara. Notice that as one of our main contributions is the efficiency in the prediction of pre-miRNAs in relation to other methods that use cubic complexity algorithms, it was natural to compare Mirinho to methods that adopt this kind of algorithm. However, we also included in the comparison a method such as CSHMM which does not use the same cubic algorithm for the prediction of miRNAs.

Since the set of input parameters differs for each method, it is not a trivial task to set them accordingly to the data, and at the same time be fair in the comparison. We then applied the methods with default parameters. However, we adapted one aspect that was common to all the methods: the set of known (pre-)miRNAs. All the methods were trained, when required, with animal (pre-)miRNA sequences. The description of each method, and how they were trained and performed is given below.

To set the initial parameters for CSHMM, we used the secondary structures of the kingdom metazoan that are available in miRBase release 20. To generate the likelihood score, the same metazoan hairpin sequences were given as the positive training set, and as the negative instances the sequences used by the authors were employed.

MIReNA provides different starting points for the prediction based on the different kinds of input files. We then chose the one that allows genomic inputs (-M option), and the set of known mature miRNAs required was from the same metazoan kingdom, taken from miRBase release 20.

miRPara makes available a script to generate the model according to the miRBase release and to the desired organism(s)/clade. In our case, we chose the model trained with metazoan pre-miRNAs of the latest release 20.

To analyse the quality of the predicted structures, we used RNAfold [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF] and miRNAFold [START_REF] Tempel | A fast ab-initio method for predicting miRNA precursors in genomes[END_REF]. The first is a classical method for predicting an RNA secondary structure through energy minimisation. If one has access to GPU facilities, we may cite two papers that implemented algorithms for such a kind of technology: [START_REF] Rizk | Gpu accelerated rna folding algorithm[END_REF] and [START_REF] Steffen | Gpu parallelization of algebraic dynamic programming[END_REF]. The second is a method for predicting a hairpin structure that takes into account specific criteria (such as length of the stem, percentage of nucleotides, size of terminal loops) related to known hairpins from miRBase, and verifies if these are present in the query structure. miRNAFold is moreover, as far as we know, the Chapter 2. Mirinho: Efficient precursor miRNA predictor only other method that has quadratic complexity for predicting a miRNA structure. For more details on how each of the methods were used, see Section 2.3.3.

Measuring sensitivity and precision

To evaluate the performance of each method in the prediction of pre-miRNAs, we used as measures sensitivity and precision (besides the stem-loop structures available in miRBase). The first is the proportion of true pre-miRNAs that are correctly predicted while the second is the fraction of predicted pre-miRNA candidates that are real pre-miRNAs:

Sensitivity = T P T P + F N (2.
3)

P recision = T P T P + F P (2.4)
where TP stands for True Positive, FP for False positive, and FN for False Negative.

To compute the number of true pre-miRNAs predicted by each method, we do the following. For each of the six species, there is a control set C = {c j , c j+1 , ...} of the miRNAs that are considered to be true miRNAs following according to miRBase, where j ∈ 1..n and n is the number of true miRNAs for a given species. Ideally, to compute the number of TPs, one should compare a predicted miRNA pm with the control miRNAs cm j that has at least one position in common with it. However, not all the softwares provide the exact coordinate of the predicted miRNA. Instead, all of them give the coordinates of the respective predicted pre-miRNA ppm. In order to compute the number of true miRNAs for a given species, we therefore verified, for each ppm, whether it fully covered a control cm j . If that was the case, we accounted for one TP. If the same ppm covered more than one control miRNA, we considered just the one with the best prediction score according to each method.

Results and discussion

Regression analysis of the free energies

To verify how close we get to the algorithms based on a secondary structure prediction, we present a regression analysis between the energies of the pre-miRNAs corresponding to the true positive pre-miRNAs predicted by Mirinho and their energies when predicted by RNAfold [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF].

Figures 2.3a-2.3f shows the relationship of the energies for the true positive pre-miRNAs of chromosome 25 of Bos taurus, chromosome I of Caenorhabditis briggsae, chromosome 2R of Drosophila simulans, chromosome 25 of Gallus gallus, chromosome 22 of Gorilla gorilla, and chromosome 19 of Mus musculus. We consider as the dependent variable the energies of Mirinho and as the independent variable the energies of RNAfold. As we can see, the energies are quite close to each other with, in general, bigger energies predicted by Mirinho. It is expected that RNAfold produces energies that are more negative than Mirinho since it minimises the free energy while the algorithm used by Mirinho maximises the number of base pairs. This provides reasonable evidence that our method approximates well the free energy of hairpins.

Time efficiency

As mentioned (in Section 2.2.1), we further improved the alignment algorithm by pruning the DP matrix and focusing on its diagonal only.

To establish the size of the diagonal portion of the DP matrix we should compute, we assessed different values for the parameter dw (diagonal width). The values for dw were evaluated empirically; they varied from 4 to 6 (see Table 2.1). A very small value for dw means to constrain the alignment to a very limited space around the diagonal part of the DP matrix, that is to permit a few or almost no bulges nor internal loops. This situation would not represent the real structure of a stem-loop and that is why we chose as minimum value dw = 4. On the other hand, a very large value for dw would not achieve the goal of the pruning strategy, that is time efficiency. In our experiments, the best results were obtained when using dw = 5 or dw = 6, which corresponds to the maximum number of unpaired nucleotides in the stem formed by both strands. The default value for the dw parameter was then set to 6. The user of Mirinho is given the freedom to compute the whole matrix instead of only its diagonal for a given value of dw. In this case, dw should be set equal to the length of the stem-arm (option -a).

Using this pruning strategy, the region exploited by the alignments is much smaller and the method performs, in general, 30% faster than the original version. Sensitivity and precision remain similar between the original and the optimised versions, in the great majority of the cases it remained the same.

Time efficiency is even more evident when comparing our method to other predictors, such as CSHMM, MIReNA, and miRPara. Table 2.2 presents the computation times for the prediction of putative pre-miRNAs in a sequence of length 4,951nt, running under a Mac OS X 10.6.8, 2.7 GHz Intel. As one can see, our method is indeed much faster than the others, making the prediction of pre-miRNAs much more feasible. To show that Mirinho is much more applicable, we compared the time of prediction of Mirinho, CSHMM, and miRpara. To facilitate the comparison of the predicted pre-miRNA candidates, we used the human chromosome 19, as the authors of CSHMM did. All three softwares were then submitted in a cluster queue of 29 hours (maximum job time without special bureaucratic request). Mirinho finished its job after 5 hours, while the other two exceeded the 29 hours without finishing their prediction, with no reported result. Clearly, one can fragment the input in smaller pieces to finish the prediction with CSHMM and MIReNA, however the message here is to show that no fragmentation is required for long sequences since Mirinho can finish its prediction in a smaller amount of time.

Method

miRNA hairpin structure prediction in sRNA-seq of plant

To obtain a high quality structure, Mirinho needs the information on the length of the stemarms and terminal loop. It is clear that, when the search is made at a genomic scale, the precise information about length is unknown. However with sRNA-seq data, the length of the stem-arms and terminal loop may be naturally inferred from the alignment of the reads against the genome. This characteristic of our method thus allows its direct application to sRNA-seq data, enriching the prediction and quality of the hairpin structures.

To demonstrate this, we started by mapping the 340,114 reads of high-throughput small RNA sequencing data from Arabidopsis thaliana (GEO accession number GPL3968) to chromosome 4 of Arabidopsis thaliana using Bowtie2 [START_REF] Langmead | Ultrafast and memory-efficient alignment of short dna sequences to the human genome[END_REF]. We considered only the mapped regions that verified the expression profile of a pre-miRNA: high coverage on (at least one of) the stem arms and lower coverage in the terminal loop. It is easy to see that the length l of a stem-arm and the length t of the terminal loop can be naturally inferred from these alignments. We then gave l, t, and the respective pre-miRNA sequence as input to Mirinho, RNAfold, and miRNAFold.

For Mirinho, we set the stem-arm length to l (option -a), and the minimum and maximum length of the terminal loop to t (options -n and -x respectively). Given that RNAfold is a method for predicting the secondary structure of an RNA in general, we used the option -C to force the structure to be a hairpin. We then required that the stem-arm regions, each of length l, were paired, and that the terminal loop region of length t was unpaired. For miRNAFold, we gave as the sliding window parameter the length of the whole pre-miRNA, that is, l+t+l.

As miRBase is the basis for miRNA studies, we took its hairpin structures as a gold standard. In order to compare the structures predicted by the three methods, we then considered three criteria: (i) number of internal loops and bulges within the stem; (ii) length of the predicted stem-arm; (iii) length of the predicted terminal loop. For each predicted structure, we verified which method produced the best result. This corresponded to the predicted structure that produced values that are closest to those of the structure in miRBase. For example, if the miRBase structure s has 3 bulges, and RNAfold predicted a structure with 2 bulges while Mirinho predicted one with 1 bulge, the first method would be considered the best one.

From the set of 50 pre-miRNAs of chromosome 4 of Arabidopsis thaliana, we randomly chose 10 structures from miRBase, for which such structures were predicted with the three methods from the sequences. In the end, Mirinho obtained the closest structure in 80% of the cases, RNAfold was the second with 50%, and miRNAFold the third with 40%.

Figures 2.4, 2.5, and 2.6 show, respectively, cases in which the closest structure was found by RNAfold, Mirinho, and miRNAFold. As we can see, even in the cases where Mirinho was not the best, it was very close to the best.

Sensitivity and precision

To determine an appropriate energy threshold for the prediction of pre-miRNAs, we used randomly generated genomes. The reasoning behind this strategy is that, if the energy model is robust enough, there should exist a certain energy that is able to differentiate the stable hairpin structures from the randomly generated ones in which the base pairs would be established by chance.

To choose the different genomes for setting the threshold, we mainly considered the GC content as it plays an important role in determining a hairpin structure. We thus chose chromosomes with different percentages of GC varying from 37% to 54%, as shown in Table 2.3.

For each of the genomes, the nucleotide frequency distribution was used to generate the respective random genomes. After that, the prediction of the pre-miRNAs was performed in both versions (original and random) of each genome. We then chose as threshold the biggest energy for which the number of true miRNAs remains zero in the random genome, as can be seen in Figures 2.7a-2.7f. To define a true positive miRNA in the random genome, we simply verified if a given true miRNA in the original genome was present in the respective random region. Using this approach, the selected genomes had the thresholds presented in Table 2. [START_REF]Number of target interactions found for the 40 miRNAs known in Acyrthosiphon pisum, using the methods for target prediction Pita and miRanda. Chapter 4 Prediction of non-coding RNAs and targets in Mycoplasma hyopneumoniae Contents 4.1 Introduction[END_REF].

We provide to the user of Mirinho an "automatic" way to set the threshold. In addition to the query genome G q , the user should give as input a similar genome G s , and an annotation file with the coordinates of the respective (true) miRNAs. Mirinho will then generate a random genome G r based on the nucleotide distribution of G s , predict its pre-miRNAs, and compute the energy threshold. If the user chooses not to provide these additional files, the default energy is set to -20.6 kcal/mol, that is the mean of the previously mentioned energies

Species

Energy threshold Chromosome GC% generated with the same approach. One should also remember that the only characteristics used by Mirinho in the prediction of pre-miRNAs are the length of the terminal loop and stem-arms and the width of the diagonal. The other methods apply additional criteria that are based on other attributes, such as AU content, sequence homology, number of unpaired nucleotides, etc. Despite this, Mirinho performs as well as the other compared methods and is at least 100 times faster than the quickest one (miRPara).

To analyse the sensitivity and precision at a genomic scale, we used the genomes of three insects, one of which, Acyrthosiphon pisum, is of particular interest to us. The results are shown in Table 2.4. Notice that the prediction is often far from being perfect for all methods; in particular, there is as usual a delicate choice to be made between sensitivity and precision, in as much as we are currently capable of accurately measuring the latter. The low precision for all the methods may be due to two reasons. One is that the model used for predicting (pre-)miRNAs needs refinement. The other is that the precise definition of a FP miRNA is completely dependent on the known miRNAs, which could represent just a small fraction of those that really exist.

Organism

Conclusion

With Mirinho, we propose a faster and flexible method for the prediction of pre-miRNAs, using minimal information about known pre-miRNAs. Concerning the prediction results, we obtain very reasonable sensitivity and precision similar to the other tested methods, and in some cases even better. As concerns the quality of the predicted structures, the hairpins predicted by Mirinho are much closer to the ones available in miRBase than the ones predicted by RNAfold and miRNAFold.

Our method is faster because we employ a quadratic time complexity algorithm to predict the free energy of the hairpin, instead of the so used cubic algorithm for prediction of RNA secondary structure. We are flexible in two aspects. First, as concerns the input type we accept both whole genome sequence and sRNA-seq data. Second, Mirinho may be used for the prediction of either plant or animal pre-miRNAs, with a minimal adjustment (of the length of the stem-arm and terminal loop only). Finally, the only a priori knowledge we use is the length of the stem-arm, the length of the terminal loop, and the width of the diagonal. The persons involved in the wet experiments were Gabrielle Duport, Karen Gaget, Federica Calevro, and Hubert Charles from the SymTrophique team at BF2I (Biologie Fonctionnelle, Insectes et Interactions, UMR0203). To treat the data in order to guarantee a more accurate set of reads, as well as to detect the expressed miRNAs, three approaches were used: (i) Chapter 3. MicroRNA expression profile during embryonic development in A.

Chapter 3

MicroRNA expression profile during embryonic development in

pisum: combining deep sequencing data and Mirinho to identify miRNAs MirinhoPipe, specially developed for this analysis; (ii) sRNA-PlAn, a pipeline designed for the annotation of small RNAs; and (iii) miRDeep, a classical method for the discovery of miRNAs from deep sequencing data [START_REF] Friedländer | Discovering microRNAs from deep sequencing data using mirdeep[END_REF]. The detected miRNAs were submitted to the prediction of mRNA targets. Together with such predictions, the gene expression profile of Acyrthosiphon pisum was analysed and compared to the miRNA expression profile, leading to very interesting results.

Introduction

The unique feeding habit of aphids combined with their ability to rapidly reproduce makes of them one of the most damaging pests of crops with economical importance worldwide.

Considering their impact on agriculture and the role miRNAs play in gene regulation, it is imperative to better characterise and understand the function of these miRNAs.

One first effort has already been made by Legeai et al. (2010a) in Acyrthosiphon pisum (the pea aphid), a laboratory model for the study of these pests whose genome was sequenced. The authors combined small RNA sequencing data from parthenogenetic females and bioinformatics approaches to identify 103 Acyrthosiphon pisum miRNAs. It is worth noting that in Legeai's work, the miRNAs of parthenogenic females were sequenced and analysed, while we focus on the miRNAs expressed in three embryonic developmental and one larval stages; all the details concerning their methodology is presented in Section 3.3.2. Furthermore, the potential mRNA targets of the detected miRNAs were identified by the overlapped predictions of two methods (Pita and miRanda).

Another effort published by [START_REF] Hansen | Widespread expression of conserved small RNAs in small symbiont genomes[END_REF], that is not directly related to small RNAs in Acyrthosiphon pisum but instead to small RNAs in its symbiont Buchnera aphidicola, provides evidence of protein regulation and of a reasonable number of conserved small RNA and UTR sequences among different Buchnera strains. The authors predicted small RNAs involved in the post-transcriptional mechanisms of the bacterium, as well as other types of mechanisms, for instance involving proteases at the same post-transcriptional level.

Material and methods

Aphid rearing and embryo isolation

A long-established parthenogenetic clone (LL01) of Acyrthosiphon pisum was maintained at 21 • C, with a 16 hour photoperiod, on Vicia faba (L. cv. Aquadulce). In order to have a supply of synchronised aphids and embryos, around one hundred mass-reared winged adults were maintained on young plants and removed after 24 h. The resulting apterous insects were maintained on Vicia faba plants for a nine-day period, until they reached the adult stage. Embryos were dissected from synchronised parthenogenetic viviparous adult aphids, removing the ovariole sheath in a buffer kept on ice. We used an RNase-free buffer composed of 35 mM Tris-HCl (pH 7.5), 25 mM KCl, 10 mM MgCl 2 , 250 mM sucrose, in 0.1% diethyl pyrocarbonate water. Following a stereoscopical analysis (Olympus IX-81, Olympus, France), embryos were classified according to their length and morphological characteristics into 3 groups (see Table 3.1 and Figure 3.1): early embryos (EE) (≤ 0.4 mm), intermediate embryos (IE) (0.4 to 0.8 mm), and late embryos (LE) (> 0.8 mm) corresponding, respectively, to the developmental stages ≤ 15, 16-18 and 19-20 as described by [START_REF] Miura | A comparison of parthenogenetic and sexual embryogenesis of the pea aphid acyrthosiphon pisum (hemiptera: Aphidoidea)[END_REF]. For L1 aged from 0 to 24h, viviparous adults were maintained on young plants for 24 hours and the resulting L1s were collected. The first column presents the four developmental groups of Acyrthosiphon pisum, followed by their abbreviation, and the developmental stage itself as described by [START_REF] Miura | A comparison of parthenogenetic and sexual embryogenesis of the pea aphid acyrthosiphon pisum (hemiptera: Aphidoidea)[END_REF]. The size of the organism is subsequently presented together with the morphological features associated to the developmental stage; the external morphological features can be observed in Figure 3.1. 3.1 for details).

RNA extraction

Total RNA was prepared using the mirVana TM miRNA Isolation Kit (Ambion, Austin, TX, USA). Three independent extractions were prepared for each group starting with 60 embryos for the EE group, 30 embryos for both the IE and LE groups, and 30 larvae for the L1 group (0-24h). The extraction was followed by a step of DNase treatment using DNA-free TM DNase Treatment and Removal Reagents (Ambion, Austin, TX, USA). Total RNA concentration and quality were initially checked using the NanoDrop R ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and samples had to meet the following quality parameters: A260/A280 ≥ 1.8 and A260/A230 ≥ 1.8, in order to be used in the subsequent analysis.

The RNA samples were then run using the Agilent RNA 6000 Nano Kit on the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) to check their integrity. Degraded samples appeared as significantly lower intensity traces, with the main peak area shifted to the lower molecular weights, and they typically exhibited much more noise on the trace. Only good quality samples were sent for sequencing.

Next-generation Illumina Sequencing

Total RNA was shipped to ProfilExpert Genomic Platform (Université de Lyon, France). RNA concentration was verified using the RiboGreen R Assay Thermo Scientific, Wilmington, DE, USA) for precise quantification before sequencing. Barcoded small RNA libraries were created from 1 ug of total RNA for each sample according to Illumina TruSeq small RNA Sample Preparation Guide (Illumina, San Diego, CA, USA): adaptor ligation was followed by RT-PCR amplification. The small RNA libraries were gel purified and they were validated on the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) to check for quality and size. The 12 samples were sequenced with a Single Read 50 cycles run on one lane of the flow cell v3 (150 million raw reads) of the Illumina HiSeq-2500 (Illumina, San Diego, CA, USA).

Treatment of the small RNA sequencing data

To identify the expressed miRNAs in the four developmental stages of the pea aphid, three methods were used: MirinhoPipe, sRNA-PlAn, and miRDeep. The first was developed specially for the purpose of analysing Acyrthosiphon pisum sRNAseq data by the author of this thesis at the BAMBOO-BAOBAB team (head Dr Marie-France Sagot) of the LBBE-UMR5558. sRNA-PlAn was developed as a pipeline for small RNA annotation by one of our collaborators, Oliver Rue, at UBIA & PF GenoToul Bioinfo (head Dr Christine Gaspin). miRDeep2 is a classical method for the discovery of miRNAs from deep sequencing data; it uses the read stacks, that are consistent to the ones of an expressed miRNA, to select the best candidates to further verify other characteristics, such as the free energy [START_REF] Friedländer | mirdeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades[END_REF]. To document and describe the technical details of the pipelines, a wiki was created and is available at http://mirinho.gforge.inria.fr/mirinhopipe.html.

To present the details on how the methods were developed and/or performed, we use Figure 3.2 as a reference guide. As shown in the figure, the first three steps were common to MirinhoPipe, sRNA-PlAn and miRDeep2. From raw data, we used Cutadapt (version 1.4.1) to trim the adapters from the 3' end (option -a) of the reads, and to filter out reads with less than 16nt (option -m). The redundancy was removed by collapsing redundant reads within each of the four samples, followed by the copy number computation of the unique reads. Only reads with a copy number greater than ten (10X) were kept. As mentioned before, for each sample three biological replicates were made. If a miRNA transcript is indeed expressed, it should thus appear in all the three replicates. Based on that, only the reads that appeared in all the three replicates were considered for subsequent analysis.

MirinhoPipe

Figure 3.2: Flowchart describing the steps of the methods for the treatment of small RNA sequencing (sRNA-seq) data. The information flows from top to bottom: (i) trimming the adapters from the 3' end and filtering out reads smaller than 16nt using Cutadapt (version 1.4.1); (ii) collapsing redundant reads and computing their copy number; (iii) only reads appearing in the three replicates remained for subsequent analysis; (iv) on the left: mapping the reads to the genome with Bowtie2; (v) computing coverage of mapped regions with genomeCov and excising potential pre-miRNAs sequences; (vi) computing the free energies of pre-miRNA hairpins with MirinhoPipe. The set of unique reads appearing in the three replicates were also given as input to miRDeep (on the right). The details concerning this last method are provided in the text.

After the cleaning process previously mentioned, the more accurate set of reads were mapped to the genome of Acyrthosiphon pisum (assembly version 2) using Bowtie2 (version 2.1.0). It was required that the reported reads mapped to at most 5 different loci in the genome (option -k), and only alignments with at most 1 mismatch were permitted (field "XM:i:<N >" from the Bowtie2 output represents the number N of mismatches). The very same set of accurate reads was also given as input to miRDeep2. The details concerning the subsequent steps of each method are provided in what follows.

MirinhoPipe

It is worth noting that all the preceding steps are included in MirinhoPipe; however, to maintain a structured presentation we split the description of the steps. From the filtered pisum: combining deep sequencing data and Mirinho to identify miRNAs Bowtie mappings, we first computed the coverage of each position in the genome using a tool called genomeCov from the toolset Bedtools (version 2.17). If a region has a coverage of a minimum height of one and a minimum length of 20nt (length of a miRNA), it is considered as a region with potential to harbour a (pre-)miRNA. To guarantee that the whole pre-miRNA is identified, flanking portions are also taken into account: if the region is smaller than 70nt, a flanking portion of 60nt down and upstream is considered and the final pre-miRNA locus is extracted.

These potential pre-miRNA loci are then given as input to Mirinho for the computation of their secondary structures and free energies. The energy threshold used is -20.6 kcal/mol, which is set as described in Section 2.3.4.

sRNA-PlAn

As for MirinhoPipe, the reads mapping to a same locus in the reference genome are assembled into a longer region resulting in a potential miRNA locus. Each locus is submitted to the annotation process and prediction of miRNA(s). To annotate a locus, non-coding RNA (ncRNA) databases are used to assign one or more putative function(s). All the loci are then submitted to a prediction step that will determine if it is a potential miRNA or not. For each potential miRNA, the up and downstream flanking portions are accounted to extract the pre-miRNA sequence. A glocal alignment is then computed between the most represented read (the putative miRNA) and its 5' or 3' neighbour region, in order to mimic the hybridisation between both strands of the potential hairpin. Each alignment is scored according to a few criteria related to the expression profile. The top ranked miRNAs are thus classified in three classes: (i) miRNA-annotated/predicted; (ii) other-function-annotated/predicted; (iii) annotation-orphan/predicted. These putative candidates can be sorted according to their score to be further submitted to experimental validation. The details about this pipeline were omitted in this thesis manuscript because it is not yet publish.

miRDeep2

The miRDeep package is composed of two modules, Mapper and miRDeep2. The first module maps the reads to the genome with Bowtie (version 1), keeping only the alignments with 0 mismatches (option -n) in the seed region. The seed region, set to 18nt (option -l), is defined as the n first nucleotides of a read. A maximum of 2 mismatches (option -e 80) occurring after the seed region were allowed (option -n ). Only reads that do not map to more than five different loci in the genome were kept (option -m). Option -best-strata was used to order the mappings (from best to worse) according to the strata definition of Bowtie [START_REF] Langmead | Ultrafast and memory-efficient alignment of short dna sequences to the human genome[END_REF][START_REF] Friedländer | mirdeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades[END_REF].

In the second module, potential miRNA precursors are excised from the genome using the read mappings as guidelines. Then the two genome strands of each genome sequence are scanned separately, from 5' to 3' . Excision is initiated when a stack of reads (height one or more) is encountered. If there is a higher read stack within 70nt downstream of the current read stack, then this is chosen instead. In this way, the highest local read stack is identified. Then the sequence covered by the highest local read stack is excised twice, once including a 70nt upstream and a 20nt downstream flanking sequence, and once including a 20nt upstream and a 70nt downstream flanking sequence. The second step of the module is to prepare the signature file. The Bowtie-build tool is used with default options to build a Burrows-Wheeler transform index of the excised potential precursors. Then the set of sequencing reads is mapped to the index, using Bowtie (version 0.12.7). The set of known mature miRNAs for the reference species is also mapped to the index. The RNA secondary structures of the potential precursors are then predicted with RNAfold with default options. Finally, the potential precursors are individually scored or discarded by the core algorithm of miRDeep2.

MicroRNA expression profile

Before verifying the expression profile of the pre-miRNAs, a normalisation of the read counts is necessary. To that purpose, the RPM (reads per million) number was computed according to Equation 3.1 (modified from RPKM-reads per kilobase per million-equation described in [START_REF] Ammar | Rpkm: The rough possibilistic k-modes[END_REF]):

RP M = r/(R s * 10 -6 ) (3.1)
where r is the number of reads mapped to the given transcript (in this case the pre-miRNA), and R s is the total number of reads from sample s that mapped to any locus of the genome. The normalisation originally includes the length of the transcript; however, we did not consider it here because miRNAs are roughly of the same length (i.e., ∼22nt).

To visualise the pre-miRNA gene expression in each sample, we used MeV (MultiExperiment Viewer), a software that was originally developed for microarray data analysis. MeV incorporates algorithms for clustering, visualisation, classification, statistical analysis, and biological theme discovery from single or multiple experiments [START_REF] Howe | Mev: multiexperiment viewer[END_REF]. The HCL option in MeV allows for the visualisation of the datasets by means of a heatmap, a graphical representation of the expression data in which the values in the matrix are represented by colors, in an organised manner, via a dendrogram, to look for emergent trends.

To build the dendogram, a hierarchical clustering is implemented in MeV. It is based on the average-linkage method developed by [START_REF] Sokal | A statistical method for evaluating systematic relationships[END_REF] for clustering correlation matrices. The algorithm assembles all elements, in this case the pre-miRNAs, into a single tree. For a set of n pre-miRNAs, an upper-diagonal similarity matrix is computed by using the Pearson correlation as a metric to score all pairs of pre-miRNAs. The pair of pre-miRNAs with the most similar expression profile is determined by finding the largest value in the matrix (i.e., the best correlation). A new cell is created by joining the two pre-miRNAs, and a new expression profile (normalised read count) is computed for the cell by averaging the expression of the joined elements. The similarity matrix is updated with this new cell replacing the two joined elements, and the process is repeated n -1 times until only a single element remains. We used a similar strategy to cluster motifs as shown in Section 5.2.2.

Target prediction

The 3'UTR target sequences were obtained at AphidBase (Legeai et al., 2010b), requiring a minimum length of 50nt for a target sequence. As a consequence, from the 40,336 sequences, 32,127 remained. Together with the targets, the set of detected miRNAs were given as input to two algorithms for the prediction of targets, Pita and miRanda. We also used RNAhybrid, however, as it produced too many interactions and none of them were common to the ones of the two other, we decided to put is aside. The main difference between miRanda and Pita, is that the latter uses accessibility to predict functional interactions; more details about each method are provide below.

miRanda uses classical features such as sequence matching, highly scoring matches at the miRNA 5' end (seed location), free energy, and target site conservation in three insect Chapter 3. MicroRNA expression profile during embryonic development in A.

pisum: combining deep sequencing data and Mirinho to identify miRNAs species Drosophila melanogaster, Drosophila obscura, and Anopheles gambiae; for validation the authors used targets from Caenorhabditis elegans and Caenorhabditis briggsae [START_REF] Enright | MicroRNA targets in drosophila[END_REF].

Pita models the target site accessibility by defining a score, ΔΔG, which is computed as the difference between the energy gained to form the duplex ΔG duplex and the energetic cost to unpair the target secondary structure ΔG open [START_REF] Kertesz | The role of site accessibility in microRNA target recognition[END_REF]. To validate their method, the authors used a quantitative luciferase assay in Drosophila melanogaster tissue. More details about the methods are presented in Section 1.2.4.

Results and discussion

Statistical summary of the sequenced reads

As mentioned in the previous sections, the sRNA-seq data were submitted to a series of processing steps before the miRNA detection itself. Table 3.2 shows the evolution of the read counts after each processing step. Filtering reads out by length (column "Cutadapt") eliminated an average of ∼36.8% of the total reads-31.7% for EE, 30.3% for IE, 35.5% for LE, and 49.7% for L1. Collapsing the redundant reads (column "Unique reads") resulted in an average removal of ∼97.2% of the reads (in relation to the previous step). Disregarding reads with a copy number smaller than 10X eliminated ∼82.2% of the reads. Finally, considering only the reads appearing in all the three replicates and only the reads that mapped to the genome, removed respectively ∼30% and ∼40,5% of the total reads. Initially the dataset was comprised of 187,357,260 reads; after all the processing and pre-treatment of this dataset, a more accurate set of 352,061 reads remained. present the read length distribution across the four samples EE, IE, LE, and L1. The length distribution is analysed from four different perspectives: (i) considering all the reads (Figure 3.3a), i.e., the reads obtained after step "Trim adapters" from Figure 3.2; (ii) considering collapsed unique reads (Figure 3.3b), i.e., the reads acquired after the step "Remove redundancy" from the same figure; (iii) considering all the reads that appear in all the three replicates (Figure 3.3c) -from the reads used in step i, only the ones appearing in all the three replicates remained; and (iv) considering collapsed unique reads that appear in all the three replicates. These four perspectives were chosen to first verify how the "noise" of redundant reads could bias the distribution, and how a more accurate set of reads (i.e. the ones appearing in all three replicates) could affect the same distribution. In the distributions with no redundancy, it is easier to see the points (which correspond to the different miRNA lengths) where the reads are concentrated and to see differences between each sample. The same occurs when we consider only the reads in the three replicates, mainly for the case "all" reads (left side of Figure 3.3). We can notice that the majority of the reads are of length 22nt, which is indeed the mean length of a miRNA. This peak at 22nt occurs for samples EE, IE and LE. For sample L1, differently from what is observed for the other samples, the peak is at 20nt, with a number of reads similar to the ones at 21nt and 22nt. Since the larval stage is developmentally farther from the three embryo stages, it is natural to expect a different behaviour for L1 in relation to the three others, and a similar behaviour within the three embryo stages. This heterogeneity of the miRNAs in L1 may be explained by the heterogeneity of the organisms from which the RNA was extract. The organisms in L1 live in a different environment than the embryos: L1 larvae are exposed to an in vivo medium, feeding on the plants, while embryos are inside their progenitor in a more stable environment. This difference in the environment results in a variation in the organism, and as a consequence (c) Considering all the reads that appear in all the 3 replicates.

MicroRNAs Expressed in Acyrthosiphon pisum

In 2010, Legeai and co-workers identified 103 mature miRNAs in Acyrthosiphon pisum using three different approaches. In the first approach, the authors blasted insect miRNAs from miRBase (release 14) against the genome of Acyrthosiphon pisum (assembly version 1.0). The second approach consisted in sequencing small RNAs from a mixed generation sample of Acyrthosiphon pisum parthenogenetic females, and mapping the 850,000 unique reads against the same genome used in the first approach. The mappings were thus given as input to miRDeep. In the third approach, the authors implemented a machine learning classifier trained with 30 pea aphid miRNAs. The ensemble of these three methods initially produced 149 mature miRNAs that were deposited in miRBase. From miRBase release 14 to the current release 21, 46 pea aphid miRNAs were removed from the database, thus remaining 103 miRNAs of Acyrthosiphon pisum.

While Legeai et al. (2010a) sequenced parthenogenic females, we sequenced the small RNAs extracted at four different developmental stages of Acyrthosiphon pisum: early embryo (EE), intermediate embryo (IE), late embryo (LE), and larvae (L1) stage. Moreover, to guarantee quality and consistency of the sRNA-seq data, for each sample three biological replicates were made (for the details about the experimental procedure, see Section 3.2). After treating the sequenced reads, the expressed miRNAs discovered with our methodology were classified in three categories: (i) miRNAs known in Acyrthosiphon pisum; (ii) miRNAs known in other species but not present in Acyrthosiphon pisum; and (iii) potential novel miRNAs. Table 3.3 summarises the number of miRNAs in each of these categories identified using the three analysis methods.

Known in Acyrthosiphon pisum

MicroRNAs known in Acyrthosiphon pisum

In miRBase (release 21), there are currently 103 pea aphid mature miRNAs deriving from 123 precursors. To identify the known miRNAs in our data, we used Blastn (version 2.2.28+) to align the 103 mature miRNA sequences against the pre-miRNAs identified by the three methods. To define a miRNA as a known miRNA in the pea aphid, we used the following criteria: (i) glocal (global+local) alignment required, local for the pre-miRNA and global for the miRNA, that is, the miRNA must be fully covered by the pre-miRNA; and (ii) maximum of one mismatch in the alignment.

From the 103 known miRNAs, MirinhoPipe retrieved 70, miRDeep 65, and sRNA-PlAn 56 miRNAs. As we can see in Figure 3.4, the number of strict consensus miRNAs between each two methods is close to the number of miRNAs predicted when considering each of the two methods separately; for example miRDeep and sRNA-PlAn find 51 strict consensus miRNAs, while miRDeep and sRNA-PlAn alone predicts, respectively, 65 and 56 miRNAs. It means that the methods are converging in their predictions and are consistent in their results. The combination MirinhoPipe+miRDeep detected 53 and MirinhoPipe+sRNA-PlAn 43 strict consensus miRNAs. The combination of all the three methods resulted in a high confidence set of 40 miRNAs known in Acyrthosiphon pisum.

From the 103 Acyrthosiphon pisum miRNAs in miRBase, 100 miRNAs have an annotation with an experimental evidence ("Evidence: experimental; Illumina"), while three miRNAs (api-miR-1923, api-miR-281, api-miR-iab-4) are annotated with an evidence obtained by similarity ("Evidence: by similarity"). It is worth noting that from these three miRNAs, miRNA api-miR-281 appears in our high confidence list (i.e., identified by the three methods), while miRNA api-miR-iab-4 was retrieved by both MirinhoPipe and miRDeep; miRNA api-miR-1923 was not recovered by any method. One difference between these three miRNAs that is important to highlight is that the high confidence miRNA api-miR-281 is expressed in 41 different species and miRNA api-miR-iab-4 in 24, while the miRNA api-miR-1923 that was not found in our predictions appears in only one species Bombyx mori. Based on these results and considering that conservation is a strong argument, the annotation field "Evidence" in miR-Base, for the miRNAs api-miR-iab-4 and api-miR-281 should be updated to experimental, as we have obtained it from sRNA-seq data.

MicroRNAs known in other species and not (yet) identified in Acyrthosiphon pisum

There are currently 32,488 miRNAs in miRBase (release 21) that are known in other species and were not (yet) identified in Acyrthosiphon pisum. From these miRNAs, 26 were identified in our data by MirinhoPipe, 21 by miRDeep, and 21 miRNAs by sRNA-PlAn, as shown in Figure 3.5. It is worth noting that to compute the preceding numbers, the unique mature miRNA sequence was considered instead of miRNA families. Although miRNAs miR-2a and miR-2b are very similar (differing in one or two bases), they are two unique miRNA sequences. The miRNAs in this category are of special interest, since they were not known to be expressed in Acyrthosiphon pisum before. To proceed with further analyses, we focused on the miRNAs identified by all the three methods (see Table 3.5), and with no family member in the pea aphid. A "family member" is a miRNA very similar in sequence; for example, miRNA dme-miR-184 was identified as "known in other species", however, the pea aphid expresses miRNA miR-184b that differs in two bases in relation to dme-miR-184. Using this definition, the only miRNA with no family member in the pea aphid is miR-79 (see Table 3.5).

We now focus on the expression profile of miRNA miR-79: first the 17 precursors giving rise to this miRNA were recovered and then their read coverage was computed. A expression profile is consistent with Dicer and Drosha processing if a few criteria are met. According to miRBase, a sequence must meet the criteria below to be annotated with high confidence:

1. At least 10 reads must map with no mismatches to each of the two possible mature miR-87-3p 5 miR-87 5 miR-87 api-miR-87a, api-miR-87b miR-9-1 1 miR-9-2 1 miR-9-3p 2 miR-9-5p 24 miR-9 17 miR-9 api-miR-9a, api-miR-9b Table 3.5: List of the miRNA genes known in other species and found in Acyrthosiphon pisum by all the three methods. The number of species expressing the miRNA is given (No. of species), followed by the corresponding miRNA families and their members present in the pea aphid. The line in bold represents the only miRNA, miR-79-5p, with no family member in Acyrthosiphon pisum.

Results and discussion

57 microRNAs derived from the hairpin precursor.

2. The most abundant reads from each arm of the precursor must pair in the mature microRNA duplex with 0-4nt overhang at their 3' ends.

3. At least 50% of the reads mapping to each arm of the hairpin precursor must have the same 5' end.

4. The predicted hairpin structure must have a folding free energy of <-0.2 kcal/mol/nt.

5. At least 60% of the bases in the mature sequences must be paired in the predicted hairpin structure.

To verify if the criteria applied to the 17 pre-miRNAs, their indexes were first built with Bowtie (2.2.0), and the reads from the four samples EE, IE, LE, and L1 were mapped to the precursors. It is worth noting that the packages of reads used were the ones obtained after trimming out adapters and filtering out reads < 16nt; unique reads were not used because the real expression profile would be "hidden" by the removal of the copies. Our set was then comprised of 34,170,694 reads for EE, 33,195,700 reads for IE, 30,084,568 reads for LE, and 21,732,072 reads for L1.

The EE reads mapped to 13 precursors, the IE reads to 13 too, the LE reads to 11, and the L1 reads to 9 precursors. Filtering out these precursors according to the criteria mentioned above, only two remained: one belonging to contig GL350203 (471,709..471,793) and the other to contig GL349650 (1,158,472..1,158,559). To make it simple, we call mir-79-GL350203 the first precursor that appeared in all the samples, and mir-79-GL349650 the second one that appeared only in sample LE.

The three first criteria are related to the pattern of the mapped reads, while the two last are related to the precursor sequence. Table 3.6 presents the different values for these criteria considering the four samples. The three first criteria apply to sample EE only, while criterion 3 did not apply to the IE, LE, and L1 samples. Furthermore, the EE sample has the largest number of reads aligning to the precursor mir-79-GL350203, as shown in Figure 3.6. To check criteria 4 and 5, which refer to characteristics of the secondary structure, we use Figure 3.7) as a guide. The first precursor mir-79-GL350203 has a secondary structure with a free energy of -27.6 kcal/mol, thus -27.6/84nt = -0.328 kcal/mol/nt (criterion 4), and 68% of the nucleotides of the miRNA duplex were paired (criterion 5). The second precursor mir-79-GL349650 has a secondary structure with a free energy of -24.52 kcal/mol, so -24.52/84nt = -0.291 kcal/mol/nt, and 86% of paired nucleotides in the duplex. As one can notice, only precursor mir-79-GL350203 together with the reads of sample EE fulfilled all the criteria.

These facts provide a strong evidence that the miRNA api-miR-79 derives from the precursor mir-79-GL350203 since it verifies all the biological criteria. Moreover, when considering the number of reads mapping to this precursor, the highest stacks are obtained with the reads from sample EE. This means that the expression profile of the miRNA api-miR-79 is more prominent during the early embryo stage of Acyrthosiphon pisum. This make us believe that api-miR-79 has an important function in the developmental process of early embryos. Although only the precursor mir-79-GL350203 (more strongly expressed in the EE stage) fulfilled all the criteria, we do not discard the hypothesis that the same precursor is being expressed in other stages, since the great majority of the criteria also applied to the IE, LE and L1 stages. Based on that, we present in Figures 5. , the expression profile of precursor mir-79-GL350203, and in Figures 5.7-5.9 expression profile of precursor mir-79-GL349650 during the three remaining stages. pisum: combining deep sequencing data and Mirinho to identify miRNAs

Novel precursor microRNAs in Acyrthosiphon pisum

Only the potential novel pre-miRNAs retrieved by all the three methods were considered for downstream analysis. To verify if a given pre-miRNA was a consensus among two or three methods, we applied a global sequence alignment (instead of a local as used for known miRNAs), since we are comparing two sequences with similar lengths (pre-miRNA sequences).

For that, we used the tool ggsearch36 implemented in the FASTA package. A pre-miRNA sequence was considered as a consensus between n methods if the alignment between the sequences outputted by the n different methods had no mismatches. Note that the sequences are not necessarily identical since there can exist gaps. Using this strategy, the number of consensus pre-miRNAs was computed for each two methods and for all the three, as shown in Figure 3.8. As we can see, MirinhoPipe obtains a larger number of potential novel pre-miRNAs because we consider all the mapped regions for prediction, while miRDeep and sRNA-PlAn eliminate unlikely regions before the prediction. This means that MirinhoPipe will also consider low expressed miRNAs while the two other methods will preferentially detect the highly expressed ones. miRDeep uses the pattern of the read stack to constrain the prediction to a smaller region while sRNA-PlAn uses a database of other kinds of RNAs (e.g. ribosomal RNAs) to eliminate "non-miRNAs".

Considering that there is no a priori knowledge about novel pre-miRNAs, an additional criterion was used: only pre-miRNAs holding a pattern of read coverage consistent with the Dicer and Drosha processing were deemed as strong candidates. As pointed out in Kozomara et Griffiths-Jones (2011), a typical pre-miRNA would meet the five criteria mentioned in the previous section: at least 10 mapped reads, 3' end overhangs, 50% of the reads with a same 5' end, a minimum free energy, and 60% of paired bases. We thus selected only the pre-miRNAs respecting these criteria. For that, we used the same approach as for the miRNA miR-79. First, the indexes of the pre-miRNAs were built with Bowtie2 (version 2.2.0) and the reads from samples EE, IE, LE, and L1 were mapped to the pre-miRNAs. From the 23 potential novel pre-miRNAs, 14 had EE read mappings, 16 had IE read mappings, 13 had LE read mappings, and 6 had L1 read mappings, all of them with more than 10 reads aligning to each mature miRNA (criterion 1). After applying the other criteria (2-5), 10 precursors remained for EE, 15 for IE, 12 for LE, and 5 for L1, giving a total of 14 unique precursor sequences, as shown in Tables 3.7 and 3.8. The first table presents the developmental stages during which the pre-miRNA is being expressed, while the second table presents the sequence of the precursor miRNA.

Precursor EE IE LE L1 1 X 2 X X X 3 X X 4 X 5 X X 6 X 7 X X 8 X X X X 9 X 10 X X X 11 X 12 X X 13 X 14 X
Table 3.7: The 14 novel precursor-miRNA sequences organised by the developmental stages in which they are expressed. To recover the precursor sequence, see Table 3.8. Table 3.8: The 14 collapsed sequences corresponding to the new precursor-miRNAs identified by all the three methods and applying to all the five criteria for a high confidence annotation. The sequences are in a 5' to 3' orientation, and the mature miRNAs are highlighted in red. Although sequences 6-10 appear to be the same, they differ from each other by one or two nucleotides. That is for instance the case of sequences 7 and 8 for which in position 41 (underlined) there is a an "A" for 7 and a "T" for the other sequence.

Chapter 3. MicroRNA expression profile during embryonic development in A. pisum: combining deep sequencing data and Mirinho to identify miRNAs

MicroRNA gene expression profile

All the previous analyses considered the ensemble of discovered miRNAs disregarding the samples from which they originated. In this section, we focus on the differential miRNA expression in each sample (EE, IE, LE, and L1). It is important to identify the miRNA genes specific to a certain stage to understand which miRNAs are biological determinants in the development of the pea aphid. To identify these genes, the expression profile of the miRNAs found with our methodology was computed with Mev (MultiExperiment Viewer), which uses an unsupervised hierarchical clustering-generated by an average linkage method with euclidean distance and no leaf order optimisation [START_REF] Howe | Mev: multiexperiment viewer[END_REF].

To measure the expression, we first normalised the read counts using Equation 3.1, which is simply the ratio between the number of reads, specific to one sample, that aligned to the given pre-miRNA and the total number of reads that aligned to the genome. These normalised counts are then submitted to Mev, and using the option "HCL", the result is a heatmap of the expression of the miRNAs; the expression profiles are arranged in clusters within a dendogram, as shown in Figure 3.9.

The expression profiles are organised in the heatmap by miRNAs vs. samples: the lines are the different miRNAs found with our methodology, and the columns are the samples (including the replicates) in which these miRNAs were expressed. In Figure 3.9, the clustering shows that our samples can be classified based on the miRNA expression levels: there are four clusters of expression profiles that precisely agree with the four different samples (see the top of the figure). We can also observe that the profiles of the stages IE and LE are clustered together, meaning that they are biologically closer to each other. Moreover, the height of the branches represents the differences between the clusters, i.e., the more distant are the samples the longer are the branches. As one can notice in the dendogram, the branches of cluster L1 are longer, and this reflects the biological conditions in which the samples were obtained. As mentioned before, embryos live inside their progenitor in a stable environment, while L1 larvae are exposed to an in vivo environment feeding on plants. The samples from the L1 stage are thus more heterogeneous (longer branches) than the ones extracted from embryos, due to this environmental conditions and/or to their age (varying from 0 to 24 hours after their birth).

When we consider the clusters from the perspective of the miRNAs, using a distance cutoff of 0.396 results in eight clusters, from which four have one single miRNA gene, while the others have, respectively, 3, 10, 14 and 39 miRNAs, as shown in the same Figure 3.9 (left part of the dendogram in blue). The miRNAs within each of the different clusters may be regulating a specific gene (or a specific set of genes), since synergism is known to play an effective role in the regulation [START_REF] Xu | MiRNA-miRNA synergistic network: construction via co-regulating functional modules and disease miRNA topological features[END_REF][START_REF] Lutter | Intronic microRNAs support their host genes by mediating synergistic and antagonistic regulatory effects[END_REF]. We will confirm this with the expression and functional analysis of the predicted targets.

MicroRNA target prediction

The three sets of miRNAs, comprised of 40 known miRNAs, 1 known miRNA in other species but not present in the pea aphid, and 14 novel miRNAs, were submitted to the prediction of targets by two methods, Pita and miRanda (see Section 1.2.4 for an introduction of these softwares). We set an energy threshold of -10 kcal/mol for miRanda, while for Pita, a negative ΔΔG was required. As mentioned in Section 3.2.6, the accessibility is measured by the difference between the free energy gained from the formation of the microRNA-target duplex and the energetic cost of unpairing the folded target to make it accessible to the microRNA (ΔΔG = ΔG duplex -ΔG open ). The ideal situation would be a small cost to open Figure 3.9: An unsupervised hierarchical clustering, generated by an average linkage method with euclidean distance and no leaf order optimisation, of the number of reads mapping to each one of the 81 identified miRNAs. The colour chart indicates expression intensities using a base 2 logarithmic scale: blue and red represent, respectively, lower (2.0) and upper ( 16) expression intensities. This expression profile was computed with MeV (MultiExperiment Viewer) [START_REF] Howe | Mev: multiexperiment viewer[END_REF].

Chapter 3. MicroRNA expression profile during embryonic development in A.

pisum: combining deep sequencing data and Mirinho to identify miRNAs

Conclusion

The first result of this work is a pipeline for the analysis of small RNA sequencing data. We preferred to develop a more flexible method, in which all the steps may be adjusted according to the data, rather than completely automating the process, and leave it like a "black box". The documentation of MirinhoPipe is available at http://mirinho.gforge. inria.fr/mirinhopipe.html.

The combination of MirinhoPipe, together with sRNA-PlAn and miRDeep, allowed us to analyse, for the first time, the miRNAs on the pea aphid parthenogenesis, that revealed several novel miRNAs with potential to play key roles in the transcription regulation during the development of this insect.

From the miRNAs discovered in this work, forty were known in Acyrthosiphon pisum. Two among these, api-miR-iab-4 and api-miR-281, were annotated with an evidence of "by similarity" in miRBase. We thus suggest that their annotation should be changed to "experimental". We found a miRNA, api-miR-79, that was not known to be expressed in the pea aphid, mainly during the early embryo developmental stage, suggesting that it may play an important role at such stage. We do not discard the possibility that api-miR-79 is also expressed during the other stages, IE, LE, and L1, since it was detected during these stages; however, the expression pattern did not fulfill all the five criteria for a (pre-)miRNA high confidence annotation [START_REF] Kozomara | mirbase: integrating microRNA annotation and deep-sequencing data[END_REF]. Twenty-three further potentially novel (pre-)miRNAs were found in our data, out of which 14 were annotated with high confidence (based on the criteria mentioned above). A few were specific to certain stage(s), while others were common to the four stages.

A clustering of the normalised expression profiles of the detected miRNAs allowed to verify the quality of the samples. Those belonging to a same stage remained clustered in sub-groups, while the samples obtained during closer stages (i.e., IE and LE) were also found within a same sub-group. The sample clusters are a reflection of development and this result indicates indirectly an important role of miRNAs in the pea aphid development.

Target prediction using two methods (miRanda and Pita), resulted in 68,787 interactions between the 40 pea aphid miRNAs and the 3'UTR sequences; 980 interactions between the miRNA api-miR-79 and its putative targets; and finally 204,163 interactions for the 23 potentially novel miRNAs.

The study of miRNAs showing differential expression in different stages, and a more detailed analysis of the predicted targets (including a comparison with the mRNA microarray based profiles [START_REF] Rabatel | Tyrosine pathway regulation is host-mediated in the pea aphid symbiosis during late embryonic and early larval development[END_REF], will allow to characterise the underlying regulatory network.

This chapter is strongly based on the paper Godinho et al. (tion). In the context of a collaboration, made possible through LIRIO, an International Associated Laboratory (LIA -Laboratoire International Associé) between the LBBE-UMR5558, and notably the BAMBOO-BAOBAB team (head Dr Marie-France Sagot), and the Laboratório de Bioinformática of the LNCC/MCT, Brazil (head Dr Ana Tereza Vasconcelos), the Master student Caio Padoan de Sá Godinho came to Lyon to work on problems related to the prediction of non-coding RNA (ncRNA) in Mycoplasma hyopneumoniae (that was the main topic of his project). The problems included the segmentation of numerical sequences that represented a predicted ncRNA, the prediction of ncRNA targets, and the analysis of conservation between intergenic regions in Mycoplasmas. Although the problems were not directly related to the regulation in eukaryotes, as miRNAs are, it is important to understand how the regulation in bacteria works since one of the perspectives for future works is to understand the regulatory interactions between the partners in a symbiotic relationship, for example the interaction model between the bacterium Buchnera aphidicola and its eukaryotic host Acyrthosiphon pisum (the pea aphid), and between the bacterium Mycoplasma hyopneumoniae and its host Suc scrofa (the swine). It is this latter case that interested Caio in his Master, and that will therefore concern us also in this chapter.

Introduction

The bacterium Mycoplasma hyopneumoniae strain 7448 is a pathogenic and obligate parasite of porcine respiratory systems. It lives adhered to the epithelium of its host respiratory tract, and together with other bacteria and viruses, it is considered one of the ethiologic agents of swine enzootic pneumonia. The disease can cause a decrease in the productivity of these animals, sometimes resulting in their death [START_REF] Byrt | Effect of enzootic pneumonia of pigs on growth performance[END_REF][START_REF] Debey | Ciliostasis and loss of cilia induced by Mycoplasma hyopneumoniae in porcine tracheal organ cultures[END_REF][START_REF] Brockmeier | Polymicrobial Diseases, chapitre 13, Porcine respiratory disease complex[END_REF].

Although some effort has already being put on understanding the infection process, the specific mechanisms relating the bacterium and the disease remain unknown. Between the different sequenced strains of the same species, only Mycoplasma hyopneumoniae J (ATCC 25934) was deemed non-pathogenic [START_REF] Gardner | Detection and quantification of intergenic transcription in Mycoplasma hyopneumoniae[END_REF][START_REF] Hsu | Molecular analysis of the p97 cilium adhesin operon of Mycoplasma hyopneumoniae[END_REF][START_REF] Nicolás | Abc transporters in Mycoplasma hyopneumoniae and Mycoplasma synoviae: insights into evolution and pathogenicity[END_REF][START_REF] Siqueira | New insights on the biology of swine respiratory tract mycoplasmas from a comparative genome analysis[END_REF].

Mycoplasma hyopneumoniae 7448 has only one known transcription factor (TF) and a complex gene expression pattern. The incomparability between the number of regulatory elements and the complexity of the gene expression of the bacterium, together with increasing evidences that ncRNAs are involved in this phenomenon, strongly encourage the search for ncRNAs in the genome of Mycoplasma hyopneumoniae 7448.

After predicting the regions with a potential to harbour ncRNA genes, additional analyses were performed in an attempt to provide more evidences to carry on with experimental validation of the ncRNAs. The first problem was related to the output of the pipeline for the prediction of ncRNAs: the pipeline was generating one single assembled ncRNA sequence where two or more different ncRNA candidates were in fact present. We solved this by applying a segmentation algorithm on these outputs. To then provide stronger evidence that the candidates were indeed functional, we performed the prediction of the ncRNA targets with a method, called Alvinho, that was specially developed for this purpose. Finally, to verify if conservation could play any role in the functionality of ncRNAs, the identity of intergenic regions was assessed between closely-related Mycoplasma species by means of a k-partite graph. Genomic motifs surrounding the ncRNA, such as promoters and terminators, were also verified to reinforce the functional evidence of the ncRNA candidates. All the three steps of the pipeline are available in the form of a script or a C++ implementation.

Material and methods

Prediction of non-coding RNAs

The main component of the pipeline for the prediction of ncRNAs is a method called Single Genome ncRNA Search (SIGRS), developed by [START_REF] Larsson | De novo search for non-coding RNA genes in the at-rich genome of dictyostelium discoideum: performance of markov-dependent genome feature scoring[END_REF]. The method uses a set of known ncRNAs, provided by the user, to guide the search for new ncRNAs with a similar nucleotide composition profile. If an annotation file is also provided, the coding regions of the genome are masked according to the annotation of known genes, and the search is concentrated in smaller regions with lesser noise.

In this work, we used a set of 816 ncRNAs from species of the class Mollicutes (the class to which the bacterium Mycoplasma hyopneumoniae 7448 belongs to) and the gene annotation of the organism. To focus the search in a more specific space, all the regions containing an annotated gene were masked, that is, the nucleotides in these positions were replaced by X's. The ncRNA sequences were obtained at the Bacterial Small RNA Database (BSRD) [START_REF] Li | BSRD: a repository for bacterial small regulatory RNA[END_REF] and the annotation file at the NCBI. The set of known ncRNAs and the masked genome are provided to SIGRS, which creates a scoring system, based on the nucleotide composition of the known ncRNAs, to transform the genome in a numerical sequence. The segments with a high cumulative sum are thus considered as an ncRNA candidate.

Scoring system

In (SIGRS) [START_REF] Larsson | De novo search for non-coding RNA genes in the at-rich genome of dictyostelium discoideum: performance of markov-dependent genome feature scoring[END_REF], a scoring system is first built and then evaluated. It considers the dependency between two consecutive nucleotides. The nucleotides of a given sequence are said to be independent if the corresponding dinucleotide frequency does not differ significantly from the one generated by chance; the test used to compute it is called G-test [START_REF] Zar | Biostatistical analysis[END_REF]. The scoring system is computed according to the result of the G-test: either the nucleotides are independent, and in this case a model M0 is used to compute the scores, or the nucleotides are dependent on the preceding adjacent neighbour, and as a consequence a model M1 is used.

To build the scoring system, the frequencies of the (di)nucleotides of both the ncRNA sequences and the masked genome sequences must be computed. We denote by f α the frequency of a 1-mer word α, with α ∈ N = {A, T, C, G}. The same notation stands for a 2-mer word αβ ∈ N × N = {AA, AC, AG, . . . , T T }, for which the frequency is represented by f αβ .

These frequencies are then used in SIGRS for the construction of a stochastic model that enables to compute the score. The random variable X t represents an element of N at time t, with probability P (X t = α), with α ∈ N . At every instant t, the element X t is concatenated with its preceding elements X 0 X 1 X 2 . . . X t-1 , thus forming a chain of nucleotides. To compute the probability P (X t = α), two stochastic models may be used, M0 and M1. Model M0 assumes that the probabilities P (X t = α) are constant and do not depend on earlier events. The vector p(α) = [p A , p C , p T , p G ] is then defined and the chain of nucleotides can be built in an iterative manner. For the model M0, it is clear that the transition probability p associated with state α is simply equal to f α . The scores s are thus assigned to each state α as shown in Equation 4.1:

s α = 10 log 2 p nc α p gf α = 10 log 2 f nc α f gf α . (4.1)
where f nc α is the frequency of the nucleotide α in the ncRNA sequence and f gf α is the frequency of the same nucleotide in the genome sequence. As for model M1, the probabilities P (X t = α) are conditioned to the previous event X t-1 , i.e., P (X t = α|X t-1 = β). This is exactly the Markov property, which is described in Equation 4.2: Chapter 4. Prediction of non-coding RNAs and targets in Mycoplasma hyopneumoniae

P (X t = α|X t-1 = β, X t-2 = γ, . . . , X 0 = ω) = P (X t = α|X t-1 = β). (4.2)
These additional criteria also apply for the definition of a Markov chain:

(i) The initial probabilities of the states P (X 0 = α) = p α are given by the vector p(α), to all α ∈ N ;

(ii ) The conditional probabilities of all other states P (X t = α|X t-1 = β) = T βα , t > 0 are determined by the transition matrix T(β, α), for all (β, α) ∈ N × N ;

(iii ) The sum of all probabilities of the same conditional status should result in α∈N T βα = 1.

For the model M1, the transition probabilities T αβ are also related to f αβ , however, they are normalised with respect to f α . The scores s αβ are therefore computed using Equation 4.3.

s αβ = 10 log 2 T nc αβ T gf αβ = 10 log 2 f nc αβ /f nc α f gf αβ /f gf α . (4.3)
It is worth noting that the scoring system is built in such a way that: (i) at least one value of s is positive and the transition probability p is not null; (ii) positive scores are assigned to profiles similar to known ncRNAs, and negative scores are assigned otherwise; and (iii) the average of the scores is negative. A formalism of the previous can be found [START_REF] Karlin | Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes[END_REF]; [START_REF] Karlin | Limit distributions of maximal segmental score among markov-dependent partial sums[END_REF]. An example of scoring S 0 and S 1 , generated respectively by the models M0 and M1, is presented in

Computing the scores of ncRNA candidates

Once the query genome is converted to a numerical sequence Λ of length n, SIGRS can identify the subsequences with high cumulative score, which is obtained by the partial sum H j i :

4.2 Material and methods 71 (4.4) where Λ x is the score in position x of the sequence; variable H j i then provides the cumulative score of the region from positions i to j.

H j i = j x=i Λ x , 0 ≤ i ≤ j ≤ n,
To determine a high representative value for a partial sum, a probability density function (PDF) is used. A PD function gives the probability to find a sequence with a larger or equal score to the one generated randomly, given the score S and the associated transition probabilities p. This probability is represented by a classical measure called e-value.

Segmentation of the outputs of SIGRS One problem with this scoring system is that, instead of outputting two regions (for instance 0-2800 and 3400-5000 in Figure 4.3) each associated with a distinct ncRNA candidate, it considered the whole region (in this case 0-5000) as a unique ncRNA candidate, which would provide the wrong answer. One solution to this problem requires the segmentation of the numerical sequence that represents the ncRNA candidate in order to identify the largest local slopes in a given sequence. In order to do this, the algorithm by Kadane [START_REF] Bentley | Programming pearls: algorithm design techniques[END_REF] designed to identify the largest cumulative sums, was adapted to find these slopes, thereby allowing for the fragmentation of the output into the correct number of candidates. The adapted algorithm is defined in what follows. Let Λ be a numerical sequence, with Λ i ∈ R, which represents the scoring system shown above. It is easy to see that the product k × n i=0 Λ i may represent the most negative score a sequence can reach, where n is the length of the numerical sequence and k a parameter optimised as follows. To measure how the different values of k modified the segmentation, we used two other parameters related to the proportion of nucleotides in accordance with known ncRNAs. The first parameter α is associated to the incorrect rejection of a true nucleotide (i.e., a nucleotide that should be in the ncRNA sequence and was discarded); in statistical hyopneumoniae hypothesis testing it is called "type I error". The second parameter β is related to the failure to reject a false nucleotide (i.e., a nucleotide that should not be in the ncRNA sequence, but it was not discarded); this measure is called "type II error". It is clear that the ideal situation is to commit no error, i.e., α = β = 0. To thus choose the best k, we varied its value from (0, 1] with a spacer of 0.01, and at each interaction the euclidean distance to the point (0, 0), that represents "no error", was computed.

Alvinho: An algorithm for the prediction of non-coding RNA targets

Alvinho was initially developed for the prediction of miRNA targets, which in turn was inspired in Mirinho (see Section 2) (Higashi et al., ress). Since base pair interaction is a common characteristic of the regulatory system of both eukaryotes and prokaryotes, the software was adapted to detect base pair interactions in bacterial systems. Although it may seem a much too simplified approach, this is the only characteristic that is precisely known and well defined. Moreover, in bacteria the regulation mediated by small RNAs is in general performed by different mechanisms involving proteins that are specific to certain types of bacteria, such as the RNA chaperone Hfq that is present only in gram-negative bacteria (which is not the case of Mycoplasma hyopneumoniae 7448) [START_REF] Storz | Regulation by small RNAs in bacteria: expanding frontiers[END_REF]. CRISPR (clustered regularly interspaced short palindromic repeats) is another mechanism of regulation in bacteria; however, this mechanism has not yet been described in Mycoplasma hyopneumoniae 7448 [START_REF] Hale | RNA-guided RNA cleavage by a crispr RNA-cas protein complex[END_REF]. Based on these facts, we decided, at least in a first step, to use base pair interaction as the only feature.

As mentioned before, Alvinho is based on Mirinho, and therefore, for the sake of concision, we will focus only on the differences between the two methods. The main one is that Mirinho is designed to compute the free energy of one single sequence that folds with itself (i.e., the free energy of hybridising the stem-arms of a hairpin) while in the case of target identification, we are dealing with the interaction between two different sequences (i.e. the ncRNA sequence, and the mRNA target sequence).

As a consequence, the alignment algorithm used in Mirinho has to be modified leading to a new one in Alvinho. In the case of Mirinho, the aligned sequences (the two stem-arms) have the same length, and a global alignment must be used. When the aligned sequences are of different lengths, which is the case for targets, a local alignment is applied. Concerning the algorithmic aspect, there are two main differences between these two approaches: (i) the base conditions; and (ii) the starting and ending point of the backtracking step. The base condition of a local alignment is presented in Equations 4.5 and 4.6 below, and the base case for a global one is presented in Equation 1.3.

W (i, 0) = i k=0 γ (4.5) W (0, j) = j k=0 γ (4.6)
where γ is the penalty for gaps. The recurrence of the local alignment is the same as for a global one (see Equation 1.4). Given a DP matrix of size m × n, to recover a local alignment, the starting point in the backtracking step is the largest value in the row i = m -1 and the ending point is the first cell with a zero value W (i, j) = 0 or W (1, 1). The starting point of a global alignment is the cell W (m -1, n -1) and the ending point is necessarily the cell W (1, 1). compsub is the set to be extended by a new vertex or shrunk by one vertex on travelling along a path of the graph. The points that are eligible to extend compsub, i.e., that are connected to all points in compsub, are collected recursively in the remaining two sets; (ii) the set candidates contains all vertices that will in due time serve as an extension to the current configuration of compsub; and (iii) the set not is the set of all vertices that have at an earlier stage already served as an extension of the current configuration of compsub and are now explicitly excluded. The algorithm generates all the extensions of a given configuration of compsub using the elements in the set candidates that are not contained in the set not. The algorithm can be summarised in five steps:

1. Select a vertex candidate.

2. Add the candidate to compsub.

3. Create new sets candidates and not from the old sets by removing all the vertices not connected to the selected candidate, keeping the old sets intact.

4. Call the extension operator to extends the formed sets.

5. Upon returning, remove the selected candidate from compsub and add it to the old set not.

A necessary condition to have a clique is that the set candidates be empty, otherwise compsub could still be extended. This condition, however, is not sufficient, because if not is non-empty, the current configuration of compsub is contained in another and is therefore not maximal. compsub is thus a clique as soon as both not and candidates are empty.

Results and discussion

Identified ncRNA candidates

After applying SIGRS to the genome of Mycoplasma hyopneumoniae 7448, and segmenting the output when it was necessary, 48 regions susceptible of harbouring ncRNA genes were identified. From these 48 regions, 36 resulted from the segmentation process of 25 potential new ncRNAs; the remaining 12 regions were known ncRNAs. Table 4.1 presents the characteristics of the 48 ncRNA candidates, including the ncRNA identifier, start, end, and length of the ncRNA sequence, the GC content, the free energy of the secondary structure of the ncRNA, and the strand from which it originated. The free energy of a folded sequence is the sum of the energies associated to each base base; frequently, the methods implement an approach to minimise this energy, since the most negative this energy, the more stable a molecule is. The free energy of the sequences were computed with RNAfold and were normalised by the length of each sequence.

Predicted non-coding RNA targets

The interaction between a ncRNA and its target may occur mainly in the UTR region; however, a fewer cases have also been observed in the coding regions. Based on this, the whole CDS was considered for target prediction. For each annotated gene, a flanking portion of 150nt downstream of the start codon and 50nt upstream of the stop codon were taken into account. The 48 putative ncRNAs together with the 698 annotated genes were then given as input to Alvinho. From the outputted interactions, only the best ones were considered, that is, for each ncRNA the interaction with the most negative free energy was taken.

From these interactions, 41,7% are associated to proteins annotated as hypothetical; Tables 4.3-4.5 present the characteristics of the identified interactions. This large percentage of hypothetical proteins agrees with the number of annotated genes (294) with the same classification, that is 42,12% (294/698) of the annotated genes. The proteins are related to the following biological functions: hybridisation, RNA translation, ABC transporters, carbohydrate metabolism, adhesins, and lipoproteins. All of these biological functions are of extreme importance to the survival of Mycoplasma hyopneumoniae, and some of them can be directly related to its pathogenicity, such as adhesins. Adhesins are cell-surface components of a bacterium that facilitate adhesion to other cells. The regulation of adhesins is thus very susceptible to be related to the process of infection of the bacterium in the swine [START_REF] Madsen | Transcriptome changes in mycoplasma hyopneumoniae during infection[END_REF]. Lipoproteins are known to be related to the immune evasion system in the swine [START_REF] Kelesidis | The cross-talk between spirochetal lipoproteins and immunity[END_REF]. The regulation of lipoproteins may then be relevant to pathogenicity [START_REF] Razin | The genus mycoplasma and related genera (class mollicutes)[END_REF]. These results sustain the hypothesis of the existence of ncRNAs as regulatory elements in the studied bacterium with fundamental roles in its survival and pathogenesis. hyopneumoniae

ID

Conserved ncRNAs

Using the approach described in Section4.2.3, only four IGRs were observed to be conserved, representing 3% of the intergenic content of Mycoplasma hyopneumoniae 7448. From the four conserved IGRs, only one was found in the predicted ncRNAs. This result is coherent with other studies showing that for other gram-positive bacteria, the level of conservation is close to null, even for closely related species [START_REF] Bibliography Acebo | Identification of 88 regulatory small RNAs in the tigr4 strain of the human pathogen streptococcus pneumoniae[END_REF][START_REF] Richter | Accessibility and conservation: General features of bacterial small RNA-mRNA interactions[END_REF].

Conclusion

The work described in this chapter is still ongoing, the results obtained in silico having now to be validated experimentally, however, from the results obtained so far, we may already conclude a few points. Using the approach implemented in SIGRS, that considers the nucleotide composition to detect potential ncRNA genes, 48 putative ncRNA were discovered in Mycoplasma hyopneumoniae 7448. The segmentation approach allowed the detection of fragmented ncRNAs that were "hidden" within longer sequences. Genes related to the life and pathogenicity of the bacterium were found to be interacting with the putative ncRNAs, an important additional evidence that reinforces the idea that the ncRNAs are indeed playing a regulation role in the bacterium. Finally, very few strong conservation was found between the IGRs of closely-related Mycoplasma species, something that was is in agreement with previous studies.

Chapter 5

Cluster analysis of structured motifs In this chapter, we present a problem related to structured motifs, which is basically a pattern, that may be composed of one or more parts separated by a certain distance, that one may look for in a sequence or a set of sequences (a more formal definition is provided in the next sections). This is an issue that may seem unrelated to the study of miRNAs but the two may however appear combined in some studies. For instance, during the internship of an undergraduate student in the team, Evgueni Jacob, the motifs associated to the miRNAs that were exported from a human tissue were analysed and classified according to their statistical significance.

The problem of finding structured motifs was first addressed by Marsan et Sagot (2000) and implemented as a software called Smile (Structured Motifs Inference and Evaluation). Depending on the parameters given to Smile, the algorithm can generate a large output that may contain redundant information. For instance, if the characteristics of the motifs are not precisely known, one should choose more permissive parameters in an attempt to recover such motifs. Here we present some clustering solutions to group motifs that may correspond to the same biological "object", and to better identify the noise that may be present in such large outputs.

Introduction

Efficiently identifying biological sites or features in a set of sequences is an essential approach to identify functional elements in a genome. Example of such elements are DNA binding sites and miRNA families (i.e., all the isoforms of a same miRNA). There are two main problems related to this identification. One is the inference of a consensus sequence for such elements, the other is the prediction of the location of the sites or features that represent true positive representatives of the corresponding elements in the set of sequences. The algorithms for the prediction of location often use the results produced by the consensus extraction methods to establish all true positive positions along a genome, although the two can also be extracted together. Indeed, this is the case of the algorithm Smile (Marsan et Sagot, 2000): it simultaneously infers consensus motifs and and locates the corresponding elements in a set of sequences. The software is available at https://team.inria.fr/bamboo/en/softwares/ smile/.

Smile implements an exact algorithm for finding motifs in a set of sequences. A suffix tree is used to represent the input sequences, which together with the strategies implemented in the algorithm, result in an efficient method for the extraction of motifs. Smile requires a number of parameters, such as the number p of parts, called boxes, that a (structured) motif may have, the minimum number of substitutions e (one per box) between the motif and its occurrence, and the minimum number of times q (which stands for quorum) that the motif has to appear among the sequences. Depending on the values of these parameters, the size of the output generated by Smile may be very large, and may contain redundant motifs, or motifs that overlap and may be considered as one single functional element. For example, the larger is the value of e, or inversely the smaller the value of q, the larger will be the output. In an attempt to organise such output by eliminating the redundancy or by grouping together motifs that correspond to a same functional element, we implemented an UPGMA (Unweighted Pair Group Method with Arithmetic Mean) algorithm to cluster the motifs that are similar according to the positions in the sequences where these motifs appear. The implementation of this algorithm was performed during the internship of Thomas Balezeau, an undergraduate student in computer science, whom I co-advised together with Marie-France Sagot. Another approach that has been explored, but not yet implemented, is the use of hashing for list intersection as a quicker estimator to find redundant motifs, or motifs that represent a single biological entity.

Materials and methods

A brief reminder on Smile

Basic definitions

A motif is a pattern that "appears" in a set of sequences. Each such "appearance" is called an occurrence. An occurrence is thus a word in a sequence, while a motif may be seen as a "representation" of a set of occurrences. Motifs thus serve to both locate and to describe certain words, their occurrences, in a set of sequences.

More formally, a structured motif (or simply motif ) is defined as an ordered set of p ≥ 1 "box(es)", with p maximum error rates (one for each box), and p -1 intervals of distance (one for each pair of successive boxes). Let Σ be the alphabet of nucleotides A, C, G, T . An element m ∈ Σ+ is said to be a motif, if there is at least one occurrence u in s such that: (i) s = xuy for x, y ∈ Σ * , and (ii) the Hamming distance (i.e., minimum number of substitutions) between u and m is no more than e, a non-negative integer. Given N sequences s 1 , ..., s N ∈ Σ * and an integer 1 ≤ q ≤ N , an element m ∈ Σ+ is said to be a valid motif if it has at least one occurrence in a quorum q of distinct sequences. From now on, we will call simply motif any valid one, given a quorum q. Notice that if e is strictly greater than 0, a motif may never appear exactly in any of the sequences of the set.

Parameters and output

The user of Smile has the option to generate a generic parameter file by providing the number of required boxes, or to manually specify the parameters in the command line. These parameters include: the name of the input and output files, the alphabet for the motifs that may be the same as for the sequences from which they are inferred or an extended IUPAC one, the quorum q (minimum number of sequences where the motifs have to appear), minimum and maximum length of the motif, number of substitutions e, and number of boxes. If motifs with more than one box are sought, additional informations must be provided, such as the minimum and maximum length of the spacer between the boxes (that is the minimum and maximum size of the interval separating the two boxes), together with the value of delta if this is strictly greater than zero (otherwise, delta does not need to be specified). Once the motifs are found and their occurrences extracted, a statistical measure is used to check whether the motifs may be considered potentially significant or not. Notice that statistical significance does not necessarily imply biological significance, but may be seen as a first filter for the latter. To that purpose, the authors compute a X 2 test (with one degree of freedom) on two contingency tables, one corresponding to what is observed, the other to what was expected under the null hypothesis. To determine what would be expected under the null hypothesis, the idea is to shuffle the original sequence(s) from which the motifs were extracted, and to count how many times the motifs found in the original dataset are present, considering a Hamming distance with the same value of e, in the shuffled dataset. The user of Smile is required to provide the number of shufflings to be performed and the size of the k-mer to be conserved when shuffling the sequences.

The output of Smile is composed of the parameters summarised in the header, the sequence of the motif, followed by a numerical encoding of the motif sequence and the number of sequences in which it appeared. The source sequence and the positions of the occurrences are listed below the motif, and finally the total number of occurrences is presented. One example of output is presented in Figure 5.1.

As mentioned before, depending on how permissive are the input parameters, this output can be very large, possibly producing motifs that are redundant in the sense that many correspond to a same functional element. This procedure may also allow to reveal motifs that are clearly noise. We made a first attempt to address this problem by implementing an UPGMA algorithm, described as follows, to organise in clusters the motifs that are very similar.

Unweighted Pair Group Method with Arithmetic Mean

UPGMA is an agglomerative hierarchical clustering method that was initially proposed by [START_REF] Sneath | Numerical taxonomy. The principles and practice of numerical classification[END_REF] and improved by [START_REF] Murtagh | Complexities of hierarchic clustering algorithms: State of the art[END_REF][START_REF] Murtagh | A survey of recent advances in hierarchical clustering algorithms[END_REF]) into an (O(n 2 ) time and O(n 2 ) space) algorithm. As the name indicates, it is a method that is unweighted (all pairwise distances contribute equally), pair group (groups are combined in pairs, dichotomies only), and arithmetic mean (the pairwise distance to a group is the mean of all the distances to each member of that group).

From a distance matrix that provides the distances (e.g., euclidean) between the pairwise points, the algorithm first finds the smallest distance. The two corresponding points are inserted as leaves in a rooted tree (that will represent the structure of the pairwise matrix). The same two points are agglomerated in a single cell (a cluster) in the matrix and the distance of this new cell to all the other points is computed as the mean of the distances to all the members of the given cluster. For example, if the distance between the points a and b is the The header is simply a summary of the input parameters provided. For each motif, its sequence and numerical encoding are presented, together with the number of sequences where the motif had occurrences. The source sequence and start positions of the occurrences are presented together, followed by the total number of occurrences.

smallest in the matrix, these two points will form a single clustered cell ab. If the distance from a to c is 0.1 and the distance from b to c is 0.2, the distance between the new cluster ab and c will be 0.15 (i.e., (0.1+0.2)/2). The algorithm can be simply summarised in the following steps:

1. Determine all interpoints dissimilarities.

2. Form a cluster from the two closest points or clusters.

3. Redefine dissimilarities between the new cluster and the other points or clusters (all the other interpoint dissimilarities remaining unchanged).

Return to

Step 2 until all points are in one single cluster.

Metric definition and matrix construction

To compute the dissimilarities between the points, in this case the motifs, we used a metric that we called "motif co-occurrence metric" and that takes into account the overlapping positions in the original sequence of two motifs. The idea is that similar motifs would more probably co-occur in the same positions. The clustering method would thus enable to group occurrences of a given motif. We used Equation 5.2 to compute the metric:

d = 1 - ∩(m 1 , m 2 ) ∪(m 1 , m 2 ) (5.2)
which is simply the ratio between the number of overlapping positions between motifs m 1 and m 2 and the sum of the lengths of both motifs. Figure 5.2 shows an example of the motif co-occurrence metric.

Figure 5.2: An example of the co-occurrence metric for motifs m 1 and m 2 . The distance between the two motifs is equal to d = 1 -∩(m 1 ,m 2 ) ∪(m 1 ,m 2 ) = 10 12 = 0.83, that is, motifs m 1 and m 2 are 83% similar to each other.

To compute the dissimilarity matrix between the motifs of an output of Smile, one needs only to parse the output file recovering all the motifs per sequence, and for each pair of motifs, to compute the dissimilarity measure shown in Equation 5.2. Once the matrix is built, it can be given as input to the UPGMA algorithm.

Initial results and discussion

As this issue was addressed at the end of this thesis, we present here only initial results. First, Smile was applied over a set of miRNA sequences derived from a human tissue during the internship mentioned at the beginning of this chapter. The mature miRNA sequences were extended 500nt up and downstream, and were given as input to Smile. Four different configurations of the algorithm were used, all of them requiring motifs with one box (p = 1), a maximum of one substitution (e = 1), over the alphabet "ATCG". The minimum and maximum lengths were respectively: 6nt and 11nt, 10nt and 11nt, 15nt (for both minimum and maximum), 14nt (for both minimum and maximum), one for each of the four configurations. The UPGMA algorithm was then fed with the four outputs of Smile, generated as described above. The number of motif occurrences per sequence and the respective number of clusters are presented in Figure 5.3 for the first configuration, that we call "configuration 611". As one may notice, when the number of occurrences increases, the number of clusters decreases, meaning that the variability is more apparent when the number of occurrences is smaller. This may be expected since the distance between the clusters is computed by the mean of the distances of their components.

Conclusion

Although this is the beginning of a study on clustering motifs, we may present a few conclusions. The number of clusters grouping the different occurrences seems to be coherent to what was expected. To verify the consistency of these clusterings, and more importantly, to determine if the grouped occurrences are biologically functionally related, we will explore different datasets for which the biological motifs are precisely described. As concerns the method and its performance, a substantial improvement, either in the implementation or in the method itself, must still be performed since it currently is time consuming. For instance, for an input of 1.6Mb, it took ∼16 minutes to compute the clusters, running under a Mac OS X 10.6.8, 2.7 GHz Intel.

Conclusion and perspectives

The most important contribution of this thesis was the development of a reliable, flexible, and much faster method for the prediction of pre-miRNAs. Mirinho predicts pre-miRNAs as well as the other tested methods, however it is orders of magnitude faster. Our method was used as the basis for other issues addressed during this thesis. It is at the heart of the pipeline MirinhoPipe for the treatment of sRNAseq data and was adapted inside the method Alvinho for the prediction of ncRNA targets. Moreover, Mirinho is currently been used in other projects of the team, for example that involve the prediction of pre-miRNAs in swines.

The efficiency and reliability of our method creates new perspectives related to the "miRNA world". The incorporation of a larger number of features for the detection of miRNAs is now possible due to the speed of our method. Such features have already been defined in [START_REF] Kozomara | mirbase: annotating high confidence mi-croRNAs using deep sequencing data[END_REF] and appear to be very precise in determining a positive (pre-)miRNA. To our knowledge, they have not yet been incorporated in any software. Another characteristic we have been explored but need to develop further is the use of targets to eliminate false pre-miRNAs. Besides possibly providing a more accurate set of pre-miRNAs, this approach would enrich the results by providing a functional overview of such molecules.

The direct application of Mirinho to sRNAseq data allows the processing of millions of reads in a more feasible time. As NGS is a constantly evolving technology, the quantity of such type of data can only increase. The efficient extraction of knowledge from such data is an essential task to provide a richer comprehension of how regulation is influencing species evolution. One point that deserves special attention as concerns sRNAseq is the large number of identified pre-miRNAs when low expression must be considered. One possible solution to this problem would be the incorporation of features in the efficient prediction of pre-miRNAs, such as the ones associated to the structure of the hairpin (e.g., minimum free energy per base) and to the location of the miRNA within the structure (e.g., in the stem with an overhang of ∼2nt at each 3' end of the miRNA duplex). As noticed during this thesis (Chapter 3), such features can be powerful in discriminating true from false pre-miRNAs.

Moving now to the sRNAseq data that we analysed: the miRNAs identified in the pea aphid, together with their putative targets, open perspectives that need to be addressed. One crucial task is the identification of the miRNAs that are being differentially expressed between the different developmental stages. If we are able to address this problem, we will be able to precisely determine which miRNA is playing a key role in each stage. Another question that needs to be treated is the huge number of interactions that were found. To this purpose, a functional analysis of the targets will be performed together with an analysis of the correlation between the expression of the miRNAs and the respective predicted targets. This will provide more accurate evidences for the potential functional interactions, and shed some light on the consequences on the development of the pea aphid.

The last issue addressed in this thesis was related to the clustering of motifs. To verify the consistency of the clusterings, and more importantly, to determine if the grouped occurrences are biologically functionally related, different datasets for which the biological motifs are better described need to be explored. A starting point would be the datasets described in [START_REF] Vanet | Promoter sequences and algorithmical methods for identifying them[END_REF] and [START_REF] Vanet | Inferring regulatory elements from a whole genome. an analysis of Helicobacter pylori σ 80 family of promoter signals[END_REF]. Once the consistency and relevance of the clusters are verified, an interesting application of the clustering approach would be in providing additional evidence for the predicted targets. In this case, the motif is the predicted miRNA, for which we know a few characteristics such as its length and the number e of substitutions between its isoforms. The sequences where the motif would be searched are mRNAs, the occurrences being potential targets for the miRNAs. This could lead to an extra verification of the interaction between miRNA and target reinforcing the target prediction results. As concerns the method and its performance, a substantial improvement, either in the implementation or in the method itself, needs to be performed since it is currently time consuming. Alternatively, other approaches, using for instance a hashing intersection list, will be investigated in future and implemented to verify which approach gives the best performance.

• miRAlign http://bioinfo.au.tsinghua.edu.cn/miralign/

• miRNAFold http://evryrna.ibisc.univ-evry.fr/miRNAFold/ As mentioned in Section 3.3.2, we present the expression profile of the two precursors that give rise to miRNA api-miR-79. Figures 5. 4-5.6 present the expression profiles of the reads mapping to the precursor mir-79-GL350203, during the three remaining developmental stages IE, LE, and L1. Figures 5.7-5.9 show the expression profiles for the precursor mir-79-GL349650 during the same three stages. For the latter precursor, the read counts in the figure are smaller than 10 (criterion 1), this is because these are the mappings of unique collapsed reads, when the reads are expanded the counts are larger than 10, hence fulfilling criterion 1. Résumé: La principale contribution de cette thèse est le développement d'une méthode fiable, robuste, et rapide pour la prédiction des pré-miARNs. Deux objectifs avaient été assignés : efficacité et flexibilité. L'efficacité a été rendue possible au moyen d'un algorithme quadratique. La majorité des prédicteurs publiés utilisaient un algorithme de complexité polynomiale de degré 3 pour évaluer la structure en épingle à tige-boucle des pré-miARNs, conduisant à des temps de calculs excessifs pour des données volumineuses. La flexibilité repose sur deux aspects, la nature des données expérimentales et la position taxonomique de l'organisme (en particulier plantes ou animaux). Mirinho accepte en entrée des séquences de génomes complets mais aussi les très nombreuses séquences résultant d'un séquençage massif de type NGS de "RNAseq". "L'universalité" taxonomique est obtenu par la possibilité de modifier les contraintes sur les tailles de la tige (double hélice) et de la boule terminale. Dans le cas de la prédiction des miARN de plantes la plus grande longueur de leur pré-miARN conduit à des méthodes d'extraction de la structure secondaire en tige-boule moins précises. Mirinho prend en compte ce problème lui permettant de fournir des structures secondaires de pré-miARN plus semblables à celles de miRBase que les autres méthodes disponibles. Mirinho a été utilisé dans le cadre de deux questions biologiques précises l'une concernant des RNAseq l'autre de l'ADN génomique. La première question a conduit à le traitement et l'analyse des données RNAseq de Acyrthosiphon pisum, le puceron du pois. L'objectif était d'identifier les miARN qui sont différentiellement exprimés au cours des quatre stades de développement de cette espèce et sont donc des candidats à la régulation des gènes au cours du développement. Pour cette analyse, nous avons développé un pipeline, appelé MirinhoPipe. La deuxième question a ermis d'aborder les problèmes liées à la prévision et l'analyse des ARN non-codants (ARNnc) dans la bactérie Mycoplasma hyopneumoniae. Alvinho a été développé pour la prédiction de cibles des miRNA autour d'une segmentation d'une séquence numérique et de la détection de la conservation des séquences entre ncRNA utilisant un graphe k-partite. Nous avons finalement abordé un problème lié à la recherche de motifs conservés dans un ensemble de séquences et pouvant ainsi correspondre à des éléments fonctionnels. L'originalité de la méthode réside dans la complexité des motifs recherchés qui peuvent être constitué de sous motifs séparés. Cela avait déjà été abordée dans une méthode robuste appelé Smile mais conduisant à des sorties très volumineuses et difficilement interprétables. Nous avons développé des solutions utilisant des méthodes de classification pour regrouper les motifs pouvant correspondre à un même élément biologique. Cette approche permet de mieux distinguer les motifs biologiquement pertinents de séquences apparaissant de manière aléatoire.

Mots-Clefs : pre-microARN; programmation dynamique; modèle de énergie du plus proche voisin; prédiction; sequecançage des petit ARNs; puceron du pois; cibles de ARN non-codants; motifs Title: MiRNA and co: Methodologically exploring the world of small RNAs Abstract: The main contribution of this thesis is the development of a reliable, robust, and much faster method for the prediction of pre-miRNAs. With this method, we aimed mainly at two goals: efficiency and flexibility. Efficiency was made possible by means of a quadratic algorithm. Since the majority of the predictors use a cubic algorithm to verify the pre-miRNA hairpin structure, they may take too long when the input is large. Flexibility relies on two aspects, the input type and the organism clade. Mirinho can receive as input both a genome sequence and small RNA sequencing (sRNA-seq) data of both animal and plant species. To change from one clade to another, it suffices to change the lengths of the stem-arms and of the terminal loop. Concerning the prediction of plant miRNAs, because their pre-miRNAs are longer, the methods for extracting the hairpin secondary structure are not as accurate as for shorter sequences. With Mirinho, we also addressed this problem, which enabled to provide pre-miRNA secondary structures more similar to the ones in miRBase than the other available methods. Mirinho served as the basis to two other issues we addressed. The first issue led to the treatment and analysis of sRNA-seq data of Acyrthosiphon pisum, the pea aphid. The goal was to identify the miRNAs that are expressed during the four developmental stages of this species, allowing further biological conclusions concerning the regulatory system of such an organism. For this analysis, we developed a whole pipeline, called MirinhoPipe, at the end of which Mirinho was aggregated. We then moved on to the second issue, that involved problems related to the prediction and analysis of non-coding RNAs (ncRNAs) in the bacterium Mycoplasma hyopneumoniae. A method, called Alvinho, was thus developed for the prediction of targets in this bacterium, together with a pipeline for the segmentation of a numerical sequence and detection of conservation among ncRNA sequences using a kpartite graph. We finally addressed a problem related to motifs, that is to patterns, that may be composed of one or more parts, that appear conserved in a set of sequences and may correspond to functional elements. This had already been addressed in a robust method called Smile. However, depending on the input parameters, the output may be too large to be tractable, as was realized in other works of the team. We then presented some clustering solutions to group the motifs that may correspond to a same biological element, and thus to better distinguish the biologically significant ones from noise that may be present in what often are large outputs from many motif extraction algorithms. Keywords: pre-microRNA; dynamic programming; nearest neighbor energy model; prediction; small RNA sequencing; pea aphid; non-coding RNA target; motifs

Faculté

  de l'Univers de Lyon Directeur : M. B. GUIDERDONI Polytech Lyon Directeur : M. P. FOURNIER Ecole Supérieure de Chimie Physique Electronique Directeur : M. G. PIGNAULT Institut Universitaire de Technologie de Lyon 1 Directeur : M. C. VITON Ecole Supérieure du Professorat et de l'Education Directeur : M. A. MOUGNIOTTE Institut de Science Financière et d'Assurances Directeur : M. N. LEBOISNE

Figure 1 . 6 :

 16 Figure 1.6: Types of RNA secondary structures.

Figure 2

 2 Figure 2.1: Stem-loop coordinates and representation. The black lines represent the stem-arms, and the stripped line represents the terminal loop.

Figure 2

 2 Figure 2.2: Pruned dynamic programming matrix according to the parameter dw (diagonal width). The alignment is concentrated only in the diagonal portion of the matrix. Alignments touching the border of the diagonal portion are disregarded.

  Chromosome 25 of Bos taurus (ρ = 0.8742886). Chromosome I of Caenorhabditis briggsae (ρ = 0.7408282). Chromosome 25 of Gallus gallus (ρ = 0.9178415).

  Chromosome 22 of Gorilla gorilla (ρ = 0.6336104).

  Chromosome 19 of Mus musculus (ρ = 0.8320502).

Figure 2 . 3 :

 23 Figure 2.3: Regression analysis of the energies predicted by Mirinho and RNAfold for six different species.

Figure 2 . 4 :

 24 Figure 2.4: From top to bottom: standard secondary structure of the pre-miRNA in miR-Base (with miRNA coloured in red), and structures respectively predicted by Mirinho, miRNAfold, and RNAfold. RNAfold obtained the best prediction for the pre-miRNA MI0019239, with the closest values of stem length, terminal loop length, and number of bulges and internal loops as in miRBase.

Figure 2 . 5 :

 25 Figure 2.5: From top to bottom: standard secondary structure of the pre-miRNA in miR-Base (with miRNA coloured in red), and structures respectively predicted by Mirinho, miRNAfold, and RNAfold. Mirinho obtained the best prediction for the pre-miRNA MI0002409, with the closest values of stem length, and number of bulges and internal loops as in miRBase.

Figure 2 . 6 :

 26 Figure 2.6: From top to bottom: standard secondary structure of the pre-miRNA in miR-Base (with miRNA coloured in red), and structures respectively predicted by Mirinho, miR-NAfold, and RNAfold. miRNAFold obtained the best prediction for the pre-miRNA MI0005382, with the closest values of terminal loop length, and number of bulges and internal loops as in miRBase.

  The energy threshold for Mus musculus is -19 kcal/mol.

Figure 2 . 7 :

 27 Figure 2.7: Number of TP miRNAs predicted when using the original and the random genomes for the different species. The vertical line represent the energy threshold that better distinct true from false pre-miRNAs.

Figure 3 . 1 :

 31 Figure 3.1: Micro-photographs of the four stages, the scale bar represents 200 μm in all photographs to allow for size comparison. The microphotographs show just one embryo stage among those belonging to the corresponding groups (see Table3.1 for details). Figure taken from[START_REF] Rabatel | Tyrosine pathway regulation is host-mediated in the pea aphid symbiosis during late embryonic and early larval development[END_REF].

  Figure 3.1: Micro-photographs of the four stages, the scale bar represents 200 μm in all photographs to allow for size comparison. The microphotographs show just one embryo stage among those belonging to the corresponding groups (see Table3.1 for details). Figure taken from[START_REF] Rabatel | Tyrosine pathway regulation is host-mediated in the pea aphid symbiosis during late embryonic and early larval development[END_REF].

  collapsed unique reads that appear in all the 3 replicates.

Figure 3 . 3 :

 33 Figure 3.3: Read length distribution across the four samples EE, IE, LE, and L1 from four different perspectives.

Chapter 3 .Figure 3 . 4 :

 334 Figure 3.4: Venn diagram of the miRNAs (known in Acyrthosiphon pisum) recovered by Mir-inhoPipe (70), miRDeep (65), and sRNA-PlAn (56). The preceding numbers in parenthesis represent the total number of miRNAs in each set.

Figure 3 . 5 :

 35 Figure 3.5: Venn diagram of the miRNAs (known in other species but not (yet) identified in Acyrthosiphon pisum) recovered by MirinhoPipe (26), miRDeep (21), and sRNA-PlAn (21). The numbers in parenthesis represent the total number of miRNAs in each set.

  .................. ...................... ...................... ...................... ...................... ...................... .................... ....................... ....................... ........................ ....................... ....................... ..................... ......................... ...................... ......................... ..................... ...................... ........................ ...................... .................. ..................... ...................... ...................... ..................... ...................... ...................... ...................... ................ ...................... .......................

Figure 3 . 6 :

 36 Figure 3.6: EE reads mapped against the precursor mir-79-GL350203 (miRNA miR-79). For the sake of the presentation, the reads used in this figure were the unique reads of the EE sample because otherwise it would not fit in the page.

Chapter 3 .Figure 3 . 8 :

 338 Figure 3.8: Venn diagram of the novel pre-miRNAs identified by MirinhoPipe (4908), miRDeep (454), and sRNA-PlAn (826). The numbers in parenthesis represent the total number of pre-miRNAs identified by each method.

Figure 4 .

 4 1 shows an example for model M0.

Figure 4 . 1 :

 41 Figure 4.1: An example of the stochastic model M0 for the construction of a nucleotide sequence.

Figure 4

 4

Figure 4 . 2 :

 42 Figure 4.2: An example of the transformation of a nucleotide sequence into a numerical sequence. The system at the top of the figure assumes that the nucleotides are independent (model M0) from adjacent neighbours. Each nucleotide is replaced by a score according to the scoring system S 0 . The system at the bottom of the figure assumes that the nucleotides are dependent (model M1) and the scoring scheme S 1 is used instead.

Figure 4 . 3 :

 43 Figure 4.3: Scoring system of SIGRS: the genome coordinates are presented in the x-axis, and the cumulative sum of the scoring is shown in the y-axis. The striped red lines represent the largest local slope.

Figure 4 . 5 :

 45 Figure 4.5: Illustrated graph concepts.

Figure 5 . 1 :

 51 Figure 5.1: An example of output of Smile. The header is simply a summary of the input parameters provided. For each motif, its sequence and numerical encoding are presented, together with the number of sequences where the motif had occurrences. The source sequence and start positions of the occurrences are presented together, followed by the total number of occurrences.

Figure 5 . 3 :

 53 Figure 5.3: Regression curve between the number of motif occurrences (per sequence) and the respective number of clusters grouping these occurrences. The regression was performed on the clustering of the Smile's output obtained with the following configuration: p = 1 box, at most e = 1 substitution, and minimum and maximum motif lengths of 6nt and 11nt.

  >GL350203:471709:471793:premirna_12522:85:+ CAATGTTGATCTCTTTGGTACTTTAGCTGTAGGTATATTTTAAAGAGACGCCCTAAAGCTTCTGTACCAATGTTATTGGCAATTG ..................... ...................... ...................... ...................... .................... ....................... ....................... ........................ ....................... ....................... ..................... ..................... ......................... ......................... ........................ ...................... ..................... ...................... ..................... ...................... ..................... ...................... ...................... ...................... ................ .......................... ...................... ................... ....................... ...................... ..........................

Figure 5 . 4 :

 54 Figure 5.4: IE reads mapping to the precursor of miRNA mir-79-GL350203. For the sake of the presentation, the reads used in this figure were the unique reads of the IE sample because otherwise it would no fit in the page.

Figure 5 . 5 :

 55 Figure 5.5: LE reads mapping to the precursor of miRNA mir-79-GL350203. For the sake of the presentation, the reads used in this figure were the unique reads of the LE sample because otherwise it would no fit in the page.

Figure 5 . 6 :

 56 Figure5.6: L1 reads mapping to the precursor of miRNA mir-79-GL350203. For the sake of the presentation, the reads used in this figure were the unique reads of the L1 sample because otherwise it would no fit in the page.

  ........................ ............................. ............................ ............................ ............................ ........................... ......................... ........................... ............................ .......................... ........................... ............................ .......................... ...........................

Figure 5 . 7 :

 57 Figure 5.7: IE reads mapping to the precursor of miRNA mir-79-GL349650. For the sake of the presentation, the reads used in this figure were the unique reads of the IE sample because otherwise it would no fit in the page.

Figure 5 . 8 :

 58 Figure 5.8: LE reads mapping to the precursor of miRNA mir-79-GL349650. For the sake of the presentation, the reads used in this figure were the unique reads of the LE sample because otherwise it would no fit in the page.

Figure 5 . 9 :

 59 Figure5.9: L1 reads mapping to the precursor of miRNA mir-79-GL349650. For the sake of the presentation, the reads used in this figure were the unique reads of the L1 sample because otherwise it would no fit in the page.

  Biological background . . . . . . . . . . . . . . . . . . . . . . . . . .

	Chapter 1
	Background
	Contents
	1.1

2 1.1.1 MicroRNA definition, history, and landscape . . . . . . . . . . . . . 2 1.1.2 MicroRNA biogenesis . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.3 RNA-induced silencing complex . . . . . . . . . . . . . . . . . . . . . 6 1.1.4 Plant microRNAs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Methodological background . . . . . . . . . . . . . . . . . . . . . . 8 1.2.1 miRBase: a reference for microRNA studies . . . . . . . . . . . . . . 8 1.2.2 Computational methods for microRNA identification . . . . . . . . . 10 1.2.3 Experimental methods for microRNA detection and quantification . 18 1.2.4 Computational methods for target prediction . . . . . . . . . . . . . 22 1.2.5 Experimental methods for microRNA target identification . . . . . . 23

  The full links to the corresponding websites are presented in the appendix of this thesis.

	Method	Type	Category Website*	Author/Paper
	MiRscan	Webserver / Standalone on demand i	MiRscan website	Lim et al. (2003)
	miRAlign	Webserver	i, v	miRAlign website	Wang et al. (2005)
	ProMiR	Standalone	ii	ProMiR website	Nam et al. (2005)
	Triplet-SVM	Standalone	ii	Triplet-SVM website	Xue et al. (2005)
	RNAmicro	Standalone for Linux	i, ii	RNAmicro website	Hertel et Stadler (2006)
	miRFinder	Standalone for Windows	ii	miRFinder website	Huang et al. (2007)
	MiPred	Webserver	ii	MiPred website	Jiang et al. (2007)
	miRank	Standalone for Windows / Linux	ii	miRank website	Xu et al. (2008)
	HHMMIR	Standalone for Linux	ii	HHMMIR website	Kadri et al. (2009)
	SSCprofiler	Webserver	ii	SSCprofiler website	Oulas et al. (2009)
	CSHMM	Webserver / Standalone on demand ii	CSHMM website	Agarwal et al. (2010)
	MIReNA	Standalone for Linux	i, iv, v	MIReNA website	Mathelier et Carbone (2010)
	miRD	Webserver	ii, iv	miRD website	Zhang et al. (2011)
	MiRPara	Standalone for Linux	ii	MiRPara website	Wu et al. (2011)
	miR-BAG	Weserver / Standalone	ii	miR-BAG website	Jha et al. (2012)
	miRDeep	Standalone for Linux	iv	miRDeep website	Friedländer et al. (2012)
	miRNAFold	Webserver	i	miRNAFold website	Tempel et Tahi (2012)

*

Table 1

 1

	.1: Main information on the current methods for the prediction of miRNAs. It includes
	the type of the prediction software (standalone or webserver), the respective website and ref-
	erence, and the category in which the method was classified: (i) filter based approaches; (ii)
	machine learning approaches; (iii) target centred approaches; (iv) mixed approaches; and (v)
	homology search methods. Such classification scheme was defined by

  .6b).

		Chapter 1. Background
	./figures/figure_motifs.pdf
		./figures/figure_pseudoknot.pdf
	(a) The different motifs a
	RNA secondary structure can	(b) A pseudoknot is formed by
	have: helix, internal loop,	at least two hairpin structures,
	bulge loop, hairpin loop, ex-	in which half of one of the
	terior loop, multi-branch loop.	stems is intercalated with the
	Figure taken from Gesteland et	other hairpin. Figure taken
	Atkins (1993).	from

  . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2 Material and methods . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.1 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.2 Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2.3 Compared methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2.4 Measuring sensitivity and precision . . . . . . . . . . . . . . . . . . . 30 2.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.1 Regression analysis of the free energies . . . . . . . . . . . . . . . . . 30 2.3.2 Time efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.3.3 miRNA hairpin structure prediction in sRNA-seq of plant . . . . . . 33 2.3.4 Sensitivity and precision . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

	Chapter 2
	Mirinho: Efficient precursor miRNA
	predictor
	Contents
	2.1

Table 2

 2

		Time (in sec)
	Mirinho	0.998
	miRPara v6.2	68.008
	MIReNA	989.958
	CSHMM	1824.474

.2: Running time comparison. Running time (in seconds) for the prediction of putative pre-miRNAs in a sequence of

length 4,951nt, on a Mac OS X 10.6.8, 2.7 GHz Intel Core i7.

Table 2 .

 2

	Caenorhabditis briggsae	-16	I	37, 76
	Mus musculus	-19	19	42, 73
	Gorilla gorilla	-19	22	47, 74
	Drosophila simulans	-21	2R	43, 93
	Gallus gallus	-24	25	54, 96
	Bos taurus	-25	25	46, 96
	Caenorhabditis elegans	-	III	35, 75
	Drosophila melanogaster	-	2R	41, 84
	Homo sapiens	-	19	50, 06

3: Energy threshold obtained with the methodology mentioned in this section, and the GC% of the different chromosomes, including the ones for test (three last lines).

Table 2 .

 2 4 presents a comparison between the different methods and Mirinho with the mean energy threshold. As we can see, in humans Mirinho has the best sensitivity (70%) and precision (50%) together with CSHMM. As concerns Drosophila melanogaster, Mirinho also has the best sensitivity (80%), while MIReNA gets the best precision (75%). For Caenorhabditis elegans, CSHMM obtains the best sensitivity (70%), and MIReNA the best precision (44.44%).

	CSHMM MIReNA miRPara Mirinho	32202.854s 918.588s 110.261s 1.667s	Precision Sensitivity Precision Sensitivity Precision Sensitivity Precision Sensitivity	23.08 60.00 50.00 10.00 13.00 60.00 50.00 70.00	26.92 70.00 75.00 30.00 08.00 60.00 61.54 80.00	29.17 70.00 44.44 40.00 04.00 20.00 35.71 50.00
				Homo sapiens	Drosophila melanogaster	Caenorhabditis elegans

Table 2 . 4 :

 24 Comparison of the sensitivity, precision, and computing time of CSHMM, MIReNA, Mirinho, The values for sensitivity and precision are given in percentage. Values in italic represent the best result for the given measure. The low precision for all the methods may be due to two reasons. One is that the model used for predicting (pre-)miRNAs needs refinement. The other is that the precise definition of a FP miRNA is completely dependent on the known miRNAs, which could represent just a small fraction of those that really exist.

	using as input
	and miRPara

Table 2

 2

		Method	Sensitivity Precision
		Mirinho	69.92	0.52
	Acyrthosiphon pisum	CSHMM miRPara	23.58 36.59	0.05 0.14
		MIReNA	24.39	3.42
		Mirinho	69.17	0.25
	Culex quinquefasciatus	CSHMM miRPara	48.51 28.33	0.10 0.07
		MIReNA	18.33	2.00
		Mirinho	78.22	0.94
	Heliconius melpomene	CSHMM miRPara	48.51 58.42	0.10 0.23
		MIReNA	31.68	7.88

.5: Sensitivity and precision of three insect genomes. The energy threshold used in Mirinho was e = -20.6. Values are given in percentage, and the ones in italic represent the best value for the given measure.

  A. pisum: combining deep sequencing data and Mirinho to identify miRNAs Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . MicroRNA gene expression profile . . . . . . . . . . . . . . . . . . . 62 3.3.4 MicroRNA target prediction . . . . . . . . . . . . . . . . . . . . . . . 62 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 This chapter is strongly based on the paper Higashi et al. (tion) (in preparation). It presents an analysis of the small RNA sequencing (sRNA-seq) data of Acyrthosiphon pisum, the pea aphid, which were obtained during this thesis at ProfilExpert Genomic Platform (Université de Lyon, France) under the supervision of Dr Stefano Colella and with financial support of Dr Marie-France Sagot-European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement no. [247073]10.

44 3.2 Material and methods . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.2.1 Aphid rearing and embryo isolation . . . . . . . . . . . . . . . . . . 44 3.2.2 RNA extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.2.3 Next-generation Illumina Sequencing . . . . . . . . . . . . . . . . . . 46 3.2.4 Treatment of the small RNA sequencing data . . . . . . . . . . . . . 46 3.2.5 MicroRNA expression profile . . . . . . . . . . . . . . . . . . . . . . 49 3.2.6 Target prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3.1 Statistical summary of the sequenced reads . . . . . . . . . . . . . . 50 3.3.2 MicroRNAs Expressed in Acyrthosiphon pisum . . . . . . . . . . . . 52 3.3.3

Table 3 .

 3

1: Description of embryonic and larval stages used for the extraction of the total RNA (subsequently submitted to Ilumina sequencing).

Table 3 .

 3 The number of unique reads and the respective copy number is computed. Reads with copy number smaller than 10X are discarded. Only the reads appearing in all the three replicates remain. Finally, the number of mapped reads, with at most 1 mismatch, is presented. The number in parenthesis is the percentage of reads that remained in relation to the column "Cutadapt".Chapter 3. MicroRNA expression profile during embryonic development in A.pisum: combining deep sequencing data and Mirinho to identify miRNAs more variation in the miRNA transcripts.

		Raw Cutadapt Unique reads	10X 3 replicates Mapped reads
	EE-1 17,009,637 12,510,677
	EE-2 15,740,057 10,338,438
	EE-3 17,214,015 11,321,579
	EE	49,963,709 34,170,694	1,065,896 192,172	186,354 105,271 (0.031%)
	IE-1 14,857,272 10,743,632
	IE-2 14,639,000	9,920,251
	IE-3 18,103,230 12,531,817
	IE	47,599,502 33,195,700	1,045,952 182,050	176,256 105,582 (0.032%)
	LE-1 17,260,430 10,881,355
	LE-2 16,128,653 10,579,881
	LE-3 13,227,401	8,623,332
	LE	46,616,484 30,084,568	903,452 161,709	157,979	92,622 (0.031%)
	L1-1 15,509,610	8,235,329
	L1-2 13,611,983	7,130,465
	L1-3 14,055,972	6,366,278
	L1	43,177,565 21,732,072	448,818	80,125	77,318	48,586 (0.011%)

2: Read counts at the different steps of the treatment workflow for the four samples EE, IE, LE, and L1-the number after the dash represents the replicate. From raw reads, the first step is to trim the adapters from the 3' end and to remove the reads smaller than 16nt with Cutadapt.

Table 3 . 3

 33 Known in other species *Potential Novel

	MirinhoPipe	70	26	4908
	miRDeep	65	21	454
	sRNA-PlAn	56	21	826
	Predicted by all the 3	40	16	23

: Summary of the miRNAs predicted by MirinhoPipe, miRDeep, and sRNA-PlAn organised in three categories: (i) miRNAs known in Acyrthosiphon pisum; (ii) miRNAs known in other species but not present in Acyrthosiphon pisum; and (iii) potentially novel miRNAs. In this table, we refer to all the discovered miRNAs disregarding the sample(s) from which they originated. The last line contains the number of strict consensus miRNAs (i.e. predicted by all the three methods). *These are the predicted miRNAs that did not fit into any of the two categories i and ii, for the final list of potential novel miRNAs more criteria were verified (see

Section 3.3.2)

.

Table 3 . 6 :

 36 Summary of the criteria of a high confident precursor, for precursor mir-79, together with the reads from the four samples. The second column (from left to right), for example, represents the EE reads mapped to the precursor mir-79 originated from the contig GL350203 of Acyrthosiphon pisum. The Criterion 1 stands for the number of reads aligning to, respectively, the 5p-arm and the 3p-arm. Criterion 2 refers to the number of overhanging nucleotides in the 3' ends. Criterion 3 is related to the number of reads in each of the arms that have the same end. On the left side of the slash are the values relative to the 5p-arm and on the right side the values relative to the 3p-arm of the precursor. Criterion 4 stands for free energy associated to each nucleotide. Criterion 5 refers to the percentage of mature miRNA,

	>GL350203:471709:471793:premirna_12522:85:+:-27.62:forward
	CAATGTTGATCTCTTTGGTACTTTAGCTGTAGGTATATTTTAAAGAGACGCCCTAAAGCTTCTGTACCAATGTTATTGGCAATT
	C A	CT T	TTT	G	TATATTTT
	A TGTTGAT C TTGGTAC	AGCT TAGG	A
	| ||||||| | |||||||	|||| ||||
	T ACGGTTA G AACCATG	TCGA ATCC	A
	T A	TT T	TCT	A	CGCAGAGA
	>GL349650:1158472:1158559:premirna_3079:88:+:-24.52:forward
	GCGTTGTGTTATCTGGCTGTTGACTTTTTCCGAAACATTCAGCCTGGTTTTTCGGAAAATCAACGGGCTCGGTGCTGTGAAAAA
	-GCG	T	T -		CT	CATTCA
	TTGTG TATC GG CTGTTGA TTTTCCGAAA	G
	||||| |||| || ||||||| ||||||||||
	AGTGT GTGG TC GGCAACT AAAAGGCTTT	C
	AAAA	C	C G		--	TTGGTC
	Figure 3.7: In the upper box is shown the secondary structure of the first precursor mir-79-
	GL350203 with a free energy of -27.6 kcal/mol, thus -27.6/84nt = -0.328 kcal/mol/nt, and
	68% of paired nucleotides in the miRNA duplex; this precursor appeared in all the samples.

In the bottom part is shown the secondary structure of the other precursor mir-79-GL349650 with a free energy of -24.52 kcal/mol, so -24.52/84nt = -0.291 kcal/mol/nt, and 86% of paired nucleotides in the duplex; this precursor appeared only in the LE sample. Secondary structures and respective free energies were computed with Mirinho. within the hairpin stem, that is paired. The criteria are precisely described in the beginning of this section. Values in gray are the ones that did not reach the minimum threshold for the given criterion.

Table 4 . 1 :

 41 Characteristics of the 48 ncRNA candidates including the normalised free energy.

		Start	End	%GC Length (nt) ΔG (kcal/mol) Strand
	SIGRS_34 138331 138427 47,42	97	-0,229	-
	SIGRS_82 139192 139291	33	100	-0,191	-
	SIGRS_80 139364 139448 38,82	85	-0,136	-
	SIGRS_36 139841 140042 41,09	202	-0,222	+
	SIGRS_35 220605 220672 42,65	68	-0,224	-
	SIGRS_15 303056 303334 39,43	279	-0,248	-
	SIGRS_26 353281 353424 45,83	144	-0,162	-
	SIGRS_3	353927 354072 44,52	146	-0,114	-
	SIGRS_59 354179 354261 38,55	83	-0,222	-
	SIGRS_30 371715 371810 39,58	96	-0,182	-
	SIGRS_20 388732 388832 43,56	101	-0,144	+
	SIGRS_69 389273 389361	38,2	89	-0,175	-
	SIGRS_66 389384 389513 33,85	130	-0,139	+
	SIGRS_6	389727 389789 55,55	63	-0,287	+
	SIGRS_12 407854 407952 38,38	99	-0,192	+
	SIGRS_38 407972 408050 39,24	79	-0,125	+
	SIGRS_72 421281 421362 41,46	82	-0,233	+
	SIGRS_16 427966 428039 39,19	74	-0,143	-
	SIGRS_8	428042 428151 42,73	110	-0,141	+
	SIGRS_18 434330 434389 53,33	60	-0,29	-
	SIGRS_23 488344 488437 43,62	94	-0,185	+
	SIGRS_40 488512 488582 36,62	71	-0,075	+
	SIGRS_17 489408 489488 45,68	81	-0,202	+
	SIGRS_31 515596 515739 38,89	144	-0,177	+
	SIGRS_14 516015 516096 53,66	82	-0,257	-
	SIGRS_75 516107 516212 34,91	106	-0,177	-
	SIGRS_33 516286 516705	0,35	420	-0,183	-
	SIGRS_11 517984 518076 45,61	93	-0,266	+
	SIGRS_100 523316 523410 31,58	95	-0,205	+
	SIGRS_118 523517 523617 32,67	101	-0,126	-
	SIGRS_27 569512 569661	38	150	-0,173	-
	SIGRS_7	570424 570484 49,18	61	-0,251	+
	SIGRS_5	571166 571236 46,48	71	-0,168	-
	SIGRS_43 574204 574274 39,44	71	-0,155	+
	SIGRS_64 574341 574513 32,37	173	-0,148	+
	SIGRS_9	583898 584023 38,89	126	-0,275	+
	SIGRS_1	585178 585304 38,58	127	-0,184	+
	SIGRS_19 585272 585351	40	80	-0,239	+
	SIGRS_22 585384 585474 40,66	91	-0,218	-
	SIGRS_29 585779 586036 38,76	258	-0,227	+
	SIGRS_13 585954 586072 42,86	119	-0,252	+
	SIGRS_52 586154 586255 32,35	102	-0,172	-
	SIGRS_25 605409 605471 52,38	63	-0,162	+
	SIGRS_4	624366 624583 42,66	218	-0,234	+
	SIGRS_24 624684 624829 41,78	146	-0,142	+
	SIGRS_32 624876 625140 39,62	265	-0,296	-
	SIGRS_10 625172 625595 35,85	424	-0,223	-
	SIGRS_2	625621 625731 40,54	111	-0,167	+

Table 4 .

 4 2: Identified interactions between the predicted ncRNAs and target genes in the forward strand. hyopneumoniae

	ΔG	(kcal/mol)	-26,2	-25,9	-25,8	-25,7	-24,95	-24,85	-24,8	-24,55	-24,15	-24,1	-23,8	-23,6	-23,2	-22,4	-22,15	-21,6	-21,5	-21,1	-20,25	-19,95
	Target end	894522c(3529)	145548c(988)	903662(190)	821667c(2373)	640064c(112)	237777(452)	134781(785)	543278(595)	119568c(656)	248057(278)	459872(1652)	581494c(479)	783353c(37)	111085(317)	248057(1053)	603236c(438)	60071(451)	182604(162)	313704(279)	1648(498)
	Target start	890363c(3516)	142634c(960)	903202(162)	817365c(2350)	637216c(84)	236186(436)	131792(755)	541342(574)	116682c(642)	246749(262)	457488(1630)	574659c(460)	781617c(1)	110613(303)	246749(1035)	602635c(424)	59315(432)	181851(143)	313085(249)	57(483)
	Target gene	AAZ54025.1|Adhesin like-protein	AAZ53482.1|Protein P102-copy 2	AAZ54028.2|hypothetical	AAZ53979.1|rpoC|DNA-directed RNA	polymerase subunit beta	AAZ53855.1|hypothetical	AAZ53581.2|lysS|Lysine-tRNA ligase	AAZ53477.2|hypothetical	AAZ53796.1|tkt|Transketolase	AAZ53468.2|uvrA|Excinuclease ABC	subunit A	AAZ53589.2|oppF-1|Oligopeptide ABC	transporter ATP-binding protein	AAZ53737.2|lipoprotein	AAZ53812.2|hypothetical	AAZ53952.1|thrS|threonyl-tRNA ligase	ABP01100.1|hypothetical	AAZ53589.2|oppF-1|Oligopeptide ABC	transporter ATP-binding protein	AAZ53827.1|rplK|50S ribosomal protein	L11	AAZ53424.2|atpB|ATP synthase subunit	a	AAZ53514.2|hemK| Protoporphirogen	oxidase	AAZ53633.2|tmk|Thymidylate kinase	AAZ53378.2|dnaA|Chromosomal	replication initiator protein dnaA
	ncRNA end	14	32	48	34	49	17	30	22	15	1 6	37	20	293	15	22	1 4	20	2 0	36	29
	ncRNA start	1	4	20	12	21	1	1	1	1	1	16	1	259	1	4	1	1	1	4	14
	ncRNA ID	SIGRS_19	SIGRS_4	SIGRS_33	SIGRS_6	SIGRS_12	SIGRS_22	SIGRS_64	SIGRS_35	SIGRS_34	SIGRS_31	SIGRS_82	SIGRS_2	SIGRS_10	SIGRS_13	SIGRS_9	SIGRS_80	SIGRS_24	SIGRS_59	SIGRS_66	SIGRS_43

Table 4 . 3 :

 43 Identified interactions between the predicted ncRNAs and target genes in the forward strand (continuation).

	ΔG	(kcal/mol)	-130	-120,8	-88,1	-83,6	-64	-59,3	-46,85	-37,7037	-32,15	-32,1	-31,25	-31,0037	-30,9		-29,9		-29,85	-29,35	-29,35	-28,85	-28,8	-28,8	-28,15	-27,4	-27	-26,8	-26,6	-26,5	-26,05
	Target end	220750c(446)	585823(416)	627274c(46)	585823(45)	304284c(30)	372231c(31)	707072(226)	366043(121)	125805(147)	765922(346)	386274(1929)	298126c(124)	602637c(541)		764123c(1252)		707072(297)	382327c(807)	767297(843)	752539c(252)	357600c(889)	821667c(2373)	116398(120)	262955c(933)	68391(213)	582533(184)	674299(1023)	122183c(1551)	426748c(221)
	Target start	220227c(379)	585408(350)	625686c(1)	585408(1)	303305c(1)	371780c(1)	704311(182)	364482(87)	125114(123)	765105(329)	383903(1906)	294444c(92)	601890c(521)		762317c(1229)		704311(272)	381480c(789)	766420(824)	751968c(224)	356707c(863)	817365c(2347)	115482(89)	260688c(916)	67405(194)	582095(142)	672369(1002)	119575c(1531)	425927c(192)
	Target gene	AAZ53569.1|rpsJ|30S ribosomal protein	S10	AAZ53814.2|hypothetical	AAZ53845.1|atpA|ATP synthase alpha	chain	AAZ53814.2|hypothetical	AAZ53623.1|smf|DNA processing protein	SMF	ABP01106.1|hypothetical	AAZ53894.2|gyrA|DNA gyrase subunit A	AAZ53675.1|permease	AAZ53473.2|tpx|thiol	AAZ53939.1|rpsD|30S ribosomal protein	S4	AAZ53687.1|ABC transporter	ATP-binding protein	AAZ53620.2|conserved hypothetical	protein	AAZ53826.1|rplA|50S ribosomal protein	L1	AAZ53937.2|PTS system,	N-acetylglucosamine-specific II ABC	component	AAZ53894.2|gyrA|DNA gyrase subunit A	AAZ53684.2|conserved hypothetical	protein	AAZ53941.1|fpg| Foramidopyrimidine	DNA gycosylase	AAZ53929.1|hypothetical	AAZ53669.1|rpsF|30S ribosomal protein	S6	AAZ53979.1|rpoC|DNA-directed RNA	polymerase subunit beta	AAZ53467.2|conserved hypothetical	protein	AAZ53595.2|nrdE|	Ribonucleoside-diphosphate reductase	AAZ53433.1|tsf|Elongation factor Ts	ABP01126.1|hypothetical	AAZ53873.1|pdhD|dihydrolipoamide	dehydrogenase	AAZ53469.2|conserved hypothetical	protein	AAZ53714.1|hypothetical
	ncRNA end	1	1	66	214	250	66	48	5	1	1	8	10	1		5		20	14	4	4	19	8	1	1	1	1	1	4	1
	ncRNA start	68	67	111	258	279	96	92	40	25	18	31	44	21		28		45	32	24	32	45	33	32	18	19	40	22	25	29
	ncRNA ID	SIGRS_35	SIGRS_22	SIGRS_2	SIGRS_29	SIGRS_15	SIGRS_30	SIGRS_12	SIGRS_32	SIGRS_5	SIGRS_4	SIGRS_72	SIGRS_14	SIGRS_16		SIGRS_17		SIGRS_38	SIGRS_1	SIGRS_23	SIGRS_11	SIGRS_25	SIGRS_31	SIGRS_33	SIGRS_18	SIGRS_9	SIGRS_24	SIGRS_19	SIGRS_59	SIGRS_100

Table 4 . 4 :

 44 Identified interactions between the predicted ncRNAs and target genes in the reverse strand. hyopneumoniae

	ΔG	(kcal/mol)	-26,05	-25,15	-25,1	-24,8	-24,7	-24,65	-24,1	-24	-23,9	-23,2	-23,1	-22,6	-22,55	-22,3	-22	-21,35	-21,15	-21,1	-19,95	-19	-16,5
	Target end	497377(179)	309525c(798)	854536c(211)	645269(191)	9075(548)	290365(2509)	154424c(716)	713967c(1777)	556944c(1572)	426748c(593)	119568c(443)	833657c(1595)	784486c(158)	412861c(140)	145548c(1397)	739677c(22)	37220(485)	237777(769)	172346(22)	148856c(242)	116398(16)
	Target start	496624(154)	308509c(781)	853368c(193)	644284(168)	7934(523)	287568(2488)	153256c(704)	709474c(1763)	554613c(1554)	425927c(578)	116682c(422)	831808c(1579)	783349c(136)	412029c(121)	142634c(1384)	739396c(1)	36517(472)	236186(756)	172033(1)	145527c(231)	115482(4)
	Target gene	AAZ53764.2|DNA methylase	AAZ53628.2|lplA|Lipoate-protein ligase A	AAZ53997.1|gcp|tRNA N6-adenosine	threonylcarbamoyltransferase	AAZ53860.2|mannose-6-phosphate	AAZ53385.1|ftsY|Cell recognition particle	receptor FtsY	AAZ53615.1|secD|Protein-export	membrane protein	AAZ53486.1|pfkA|6-phosphofructokinase	AAZ53898.2|polC|DNA polymerase III	polC-type	AAZ53807.2|conserved hypothetical	protein	AAZ53714.1|hypothetical	AAZ53468.2|uvrA|Excinuclease ABC	subunit A	AAZ53985.1|dam|DNA adenine methylase	AAZ53953.2|trpS|Tryptophanyl-tRNA	ligase	AAZ53705.1|lipoprotein	AAZ53482.1|Protein P102-copy 2	AAZ53919.1|PTS system	galactitol-specific enzyme IIB component	AAZ53407.2|hypothetical	AAZ53581.2|lysS|Lysine-tRNA ligase	AAZ53507.1|rpmA|50S ribosomal protein	L27	AAZ53483.1|Protein P97-copy 2	AAZ53467.2|conserved hypothetical	protein
	ncRNA end	1	5	1	1	1	1	24	1	1	1	1	1	1	1	1	147	9	4	48	4	82
	ncRNA start	25	23	19	24	27	22	36	16	20	16	22	18	25	19	15	168	22	17	69	15	94
	ncRNA ID	SIGRS_20	SIGRS_6	SIGRS_13	SIGRS_69	SIGRS_43	SIGRS_75	SIGRS_27	SIGRS_8	SIGRS_34	SIGRS_118	SIGRS_64	SIGRS_26	SIGRS_52	SIGRS_40	SIGRS_7	SIGRS_36	SIGRS_3	SIGRS_10	SIGRS_80	SIGRS_82	SIGRS_66

Table 4 . 5 :

 45 Identified interactions between the predicted ncRNAs and target genes in the reverse strand (continuation).

	4.4 Conclusion

http://www.illumina.com/

http://www.454.com/

http://solid.appliedbiosystems.com/

Chapter 2. Mirinho: Efficient precursor miRNA predictor

(a) An example of an undirected graph. (b) An example of a k-partite graph, in this example k = 3. (c) Example of a clique of size n = 4 (in orange).

Acknowledgements

the target structure, combined with a strong base pairing between the duplex (more negative energy). Based on that, we set the threshold to ΔΔG to be at least negative. Table 3.9 presents the total number of interactions found for each Acyrthosiphon pisum miRNA by the two methods, and the number of common predictions between both.

To guarantee a more reliable set of interactions, we used a similar strategy as the one used for finding miRNAs: only those interactions predicted by the two methods were kept. While miRanda and Pita separately predicted respectively 80,191 and 1,408,997; 68,787 interactions were found by both methods for the 40 pea aphid miRNAs. Considering the only miRNA known in other species, api-miR-79, miRanda predicted 1,336 and Pita 8,967, while the overlap consisted in 980 interactions for this miRNA. For the 23 potential novel miRNAs, miRanda found 358,360 interactions and Pita 2,217,052, with 204,163 in common.

Once the local alignment applied, the same nearest neighbour energy model can be used over each two consecutive base pairs to recover the corresponding free energy; all the details concerning the energy model is presented in Section 1.2.2. As concerns the energy model, there is one small difference between the one used here and the one used for a global alignment. Since the hybridisation is given between two different sequences, here the energy model does not account for the symmetry correction for self complementary duplexes ΔG sym in Equation 1.8. The implementation of Alvinho produces an output like the one shown in Figure 4.4.

It is available at https://sourceforge.net/p/alvinho/code/ci/master/tree/.

Conservation analysis

To verify if conserved ncRNAs were more susceptible to be functional in Mycoplasma hyopneumoniae 7448, a conservation analysis was performed.

Before describing how this was done, a few definitions are necessary. A graph G = (V, E) composed of a set V of vertices and a set E of edges is said to be undirected if the edges have no direction, that is, the relations between pairs of adjacent vertices are symmetric. Two vertices are said to be adjacent if there is an edge connecting them. A k-partite graph G = (V, E) is a graph whose vertices can be decomposed into k disjoint sets so that a pair of vertices is adjacent if and only if the two vertices belong to two different sets. A clique in an undirected graph G = (V, E) is a subset C of V such that the subgraph G of G induced by C is complete, that is, for every two vertices in C, there exists an edge connecting them in G (and thus in G). A k-partite graph may have cliques of size at most k (i.e., having k vertices). Figures 4.5a,4.5b,and 4.5c show examples of respectively a graph, a k-partite graph, and a clique in a graph.

To identify the conserved ncRNAs, four species were considered: Mycoplasma hyopneumoniae 7448, Mycoplasma hyorhinis HUB1, Mycoplasma synoviae 53, andMycoplasma agalactiae PG2. The set of intergenic regions (IGRs) where the ncRNAs may be found in each species was composed of, respectively, 567, 518, 511, and 630 IGRs. Each of the four IGR sets represents one subset (one partition) of the set of vertices of the k-partite graph; in this case k = 4. Two vertices u and v belonging to two different subsets (partitions) are adjacent, if and only if an identity of I(u, v) > 70% between the two IGR sequences labelling the vertices was verified. To compute the identity between the sequences, Blast was used [START_REF] Altschul | Basic local alignment search tool[END_REF]. A sequence was considered as conserved if and only if a clique of size four was associated to it. The algorithm of [START_REF] Bron | Algorithm 457: finding all cliques of an undirected graph[END_REF] was used to list all the cliques of size k = 4.

The authors first define three sets that are essential for the core algorithm: (i) the set

Algorithm

The algorithm implemented in Smile solves the problem of identifying motifs and is described as follows. Given a set of N sequences s 1 , ..., s N , a non-negative integer e, and a positive integer q ≤ N , the goal is to find all the motifs ((m 1 , .

Notice that we have a single motif when p = 1 and d min = d max = 0 (in this case, by default, δ 1 = 0), otherwise we have a structured motif composed of p boxes.

As mentioned, a suffix tree T is used to represent the set of sequences s 1 , ..., s N . Suffix trees were introduced by [START_REF] Mccreight | A space-economical suffix tree construction algorithm[END_REF], and modified by [START_REF] Gusfield | Algorithms on strings, trees and sequences: computer science and computational biology[END_REF] and [START_REF] Bieganski | Generalized suffix trees for biological sequence data: Applications and implementation[END_REF] to consider N ≥ 1 sequences. To extract all the valid single motifs m ∈ Σ k≥1 with a number e of substitutions allowed and appearing in at least q (quorum) sequences, Marsan et Sagot (2000) implemented an algorithm that traverses simultaneously and recursively the lexicographic trie M of all possible motifs of length k and the suffix tree T of the sequences. The algorithm is based on a recurrence that is stated by the following lemma:

and only if, one of the following two conditions is verified:

(match) A pair (parent(v), e v ) is a node-occurrence of m and the label of the arc from parent(v) to v is α;

(subst.) A pair (parent(v), e v -1) is a node-occurrence of m and the label of the arc from parent(v) to v is β = α.

As for structured motifs, the lemma above together with extensions described in Marsan et Sagot (2000) are used.

Algorithm 1

Here we describe the procedure used in Marsan et Sagot (2000) to find structured motifs of the type ((m 1 , m 2 ), (d min , d max )), i.e. with p = 2 and δ 1 = 0. In other words, we want to find a structured motif with two boxes separated by a fixed interval (that can be a fixed length if

Using the suffix tree T , the first motif of length k can be found together with its set V 1 of T -node-occurrences (which are nodes located at level k in T ). Once an occurrence of motif m 1 is found to finish at node v of the tree T , a "jump" from level(v) to level(w), with d min ≤ level(w)level(v) ≤ d max , is performed. The node w corresponds to the potential starts of node-occurrences of w of motif m 2 , with w ∈ V 2 , such that:

From a node-occurrence v in V 1 of motif m 1 , a jump is thus made in T to all potential start node-occurrences w of m 2 . If the nodes v in V 1 and the nodes w in V 2 satisfy the recurrence formula given in lemma 1, the structured motif ((m 1 , m 2 ), (d min , d max )) is verified.

To find structured motifs with p > 2 and δ 1 > 0, the authors extended the algorithm accordingly. For a detailed description, see Marsan et Sagot (2000).

Appendix A

List of the websites of the corresponding methods for miRNA prediction.