Introduction

A Microbe or micro-organism is 'an organism that cannot be seen by the naked eye'.

Ancient and recent advances in technology have allowed scientists to study microbes at very different levels. Since their first observations through first microscopes in the 17 th century [START_REF] Payne | The Cleere Observer: A Biography of Antoni Van Leeuwenhoek[END_REF], micro-organisms became gradually the subject of many studies and experiences. In the 18 th century, the theory of spontaneous generation which states that life could arise from inanimate matter, gave rise to many debates among scientists where opponents like Francesco Redi conducted different experiments to disproved this theory [START_REF] Leikola | Francesco Redi as a pioneer of experimental biology[END_REF]. It was in the mid and late 1800s when spontaneous generation became completely disproved by Louis Pasteur through the swannecked flasks experiment. At the same time, Pasteur postulated the germ theory of disease, which states that microorganisms are the causes of infectious disease [START_REF] Schwartz | The life and works of Louis Pasteur[END_REF]. This theory drew the attention on the impact these organisms may have in human health. Since then, microbiology became an important field in science and the roles of microbes were eventually unveiled.

Nowadays, it is well accepted that microbes are essential to maintain all forms of life on earth.

With the advances of the molecular era, the importance of micro-organisms in the development of new technologies has been emphasized. They became crucial for the development of human activity as they can act as biological tools to produce nutrients [START_REF] De Vrieze | The littlest farmhands[END_REF], energy [START_REF] Tickell | From the Fryer to the Fuel Tank: The Complete Guide to Using Vegetable Oil as an Alternative Fuel[END_REF] and a great variety of molecules which can be used to fight diseases, to produce materials and many other applications.

Recent studies have focused on the capacity of micro-organisms to live in multi-species communities. Microbes create complex assemblies where they can be engaged in competition or cooperative behaviors, living in a continuous exchange of molecules [START_REF] Phelan | Microbial metabolic exchange--the chemotype-to-phenotype link[END_REF]. These communities ask to scientists to think on a broad scale, shifting their focus from 'how does an organism work?' to 'who is here and what are they doing?'. It is therefore crucial for the understanding of microbial communities, to uncover their structure and to characterize the molecular mechanisms involved in these interactions.

Micro-organisms as well as all cellular organisms have their genomic information coded in their DNA (deoxyribonucleic acid), a double-stranded molecule which is made from 4 different sub-units called nucleotides [START_REF] Alberts | Molecular Biology of the Cell[END_REF]. DNA contains specific regions called genes from which cellular products such as proteins are synthetized. Proteins are macro molecules which are present in all living organisms and they are the principal components of cells. Their functions are very diverse and they are involved in almost every cell process [START_REF] Andersen | Protein Structure, Stability, and Folding[END_REF]. Proteins along with all other gene products are the machinery preserving life and they have been shaped over million of years through evolutionary processes giving rise to the astonishing diversity of living beings we observe today.

Over the past few years, the technologies of DNA sequencing have evolved considerably. Sequencing consists in obtaining the entire sequence of nucleotides from a DNA molecule. This has allowed scientists from all over the world to determine the entire genomic sequences for a great diversity of organisms. Genomics has emerged therefore as a discipline whose aim is to decrypt the molecular information coded in genomes with the purpose of unveil the relevant mechanisms sustaining life. Genomes are typically reconstructed by assembling short sequences of DNA produced by Next Generation Sequencing (NGS) technologies, it is therefore very important to sample enough DNA fragments covering the entire genomic sequence in order to accomplish a proper reconstruction and to correct NGS errors [START_REF] Pareek | Sequencing technologies and genome sequencing[END_REF]. This technic relies therefore on the ablity of organims to be grown in a controlled environment in order to isolate enough DNA material. However, it is estimated that only less than one percent of the unicellular species can be grown in such conditions [START_REF] Pham | Cultivation of unculturable soil bacteria[END_REF]. The reason is the inability to replicate essential aspects of their environment which may include the lack of other species usually present in their environments [START_REF] Stewart | Growing unculturable bacteria[END_REF].

To overcome this limitation, a new discipline called metagenomics has emerged as a compelling tool to study microbial communities [START_REF] Bragg | Metagenomics using next-generation sequencing[END_REF]. It is based on the analysis of genomic sequences that are obtained directly from environmental samples. Driven by recent advances in NGS technologies, metagenomics allows scientists to establish a very extended catalogue of genes which are present in environmental microbial communities. In addition to metagenomics, NGS has also permitted the development of metatranscriptomics, a discipline that allows scientists to analyze changes in gene expression of the communities under specific environmental conditions.

In parallel to the recently generated metagenomics data, genomics studies have populated DNA and protein sequence databases with millions of entries. Due to this huge increment of data, the scientific community has made a great effort to classify and annotate existent genomic sequences. Experimental and algorithmic approaches have allowed scientists to create big repositories where sequences are classified by different criteria but mainly by their homology. Two sequences are said to be homologous if they share common ancestry and therefore, they are said to belong to the same family [START_REF] Dayhoff | The origin and evolution of protein superfamilies[END_REF]. This common ancestry can be inferred by sequence similarity between a pair of aligned sequences and other much more refined technics such as the characterization of a set of homologous sequences by probabilistic models [START_REF] Eddy | Profile hidden Markov models[END_REF][START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF]. It is well accepted that sequences sharing high similarity have also a high probability to share the same function [START_REF] Dayhoff | The origin and evolution of protein superfamilies[END_REF]. Therefore, in order to unveil the structure, the processes and the interactions in a microbial community it is crucial to annotate environmental DNA sequences obtained from metagenomic or metatranscriptomic studies by searching homologous sequences among the ones already referenced in well annotated databases.

DNA sequences from which proteins are synthetized are called coding sequences and they are often translated into proteins in order to achieve statistically significant alignments in homology detection. The protein alphabet [START_REF] Eddy | Accelerated Profile HMM Searches[END_REF] is larger than the DNA alphabet (4), therefore, it is easier to achieve statistically significant alignments due to the fact that it is less likely to align protein sequences by chance. DNA coding sequences can be inferred by different technics.

Ab initio technics are based on the statistical analysis of their nucleotide composition [START_REF] Goel | A comparative analysis of soft computing techniques for gene prediction[END_REF].

Once the coding sequences in an organism have been established, they can be translated into protein sequences which can be annotated by homology in order to determine the function of each protein (ideally) in the genome. In metagenomics we face a similar but much more complicated workflow as data coming from these experiments contain only small fragments from several distinct genomes. The increased complexity of the data poses computational challenges in assembling fragments from multiple organisms present in microbial communities.

Complications stem from the difficulty of annotating, and classifying the short sequence fragments typically obtained with next generation sequencing methods. So, novel computational methods are needed to address these issues and the massive amounts of sequence data that have become available.

Protein domains are functional and/or structural conserved units in a protein which are responsible for a particular function or interaction [START_REF] Janin | Structural domains in proteins and their role in the dynamics of protein function[END_REF]. Some domains may be present in distinct proteins that may have different biological functions, however, in most of the cases the detection of a domain may give insights of the function of the entire protein. During this thesis, we have developed MetaCLADE, a new methodology that improves the detection of well characterized protein domains in sequences obtained from metagenomics and metatranscriptomics studies. This method allows us to determine with more precision if a gene in a microbial community has already been characterized or whether it is a new sequence. To learn about the evolutionary origins of genes, it is very important to identify new families of proteins. It is also important in the understanding of the families of proteins to elucidate the evolutionary paths that have not been previously discovered. These paths can contribute to elucidate the functional activity of the communities.

For the development of MetaCLADE, we modified a protein domain annotation system that has been developed at the Laboratory of Computational and Quantitative Biology called CLADE (CLoser sequences for Annotations Directed by Evolution) [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF]. In general, methods for the annotation of protein domains characterize them with probabilistic models. These probabilistic models, called Sequence Consensus Models (SCMs) are built from an alignment of a set of homologous sequences belonging to different phylogenetic clades and they represent the consensus at each position of the alignment. However, when the sequences that form the homologous set are very divergent, the signals of the SCMs become too weak to be identified and therefore the annotation fails. In order to solve this problem of annotation caused by weak profiles constructed from very divergent sequences, CLADE uses an approach based on the observation that many of the functional and structural constraints in a protein are not broadly conserved among all species, but they can be found locally in phylogenetic clades. The approach consists to expand the catalogue of probabilistic models for a domain by creating new models that highlights the specific characteristics of each clade. This strategy has been applied on genomes difficult to annotate and it has been proven that the same evolutionary paths followed by domains contained in these genomes can be found in distant species.

To develop CLADE annotation system, almost 2.5 millions of probabilistic models called Clade Center Models (CCMs) have been constructed from a set of homologous sequences belonging to each domain family in the Pfam database [START_REF] Punta | The Pfam protein families database[END_REF]. These profiles were used to determine, by homology, if a new sequence can be classified into a family (group of proteins that share a common evolutionary origin). To produce a reliable detection of domains, CLADE filters the predictions of the models by consensus (SCMs) and the models characterizing the specific properties of each clade (CCMs) by means of a learning/classifying step. During the classifyng step a score is asigned to each prediction in a sequence. The score is dependent of several attributes that are calculated from the output of the set of models predicting hits in the sequence. Finally, a new algorithm: DAMA (Domain Annotation by a multi-objective approach) [START_REF] Bernardes | A multi-objective optimization approach accurately resolves protein domain architectures[END_REF] based on multiple optimization finds the most likely domain architecture for each protein. When this method is applied on the genome of Plasmodium falciparum, it predicts domains in proteins for which current methods fails to identify signficant signals (see chapter 2).

The extended library (or database) of probabilistic models generated by CLADE is used by MetaCLADE to find significant hits in environmental protein sequences. As mentioned before, the characteristics of metagenomic / metatranscriptomic sequences are quite different to sequences that have been predicted from entire genomes. For this reason, the learning/classifying step has been completely changed in MetaCLADE. Indeed, the specific properties of these environmental datasets such as sequence length, coverage and sequencing errors do not allow the use of the same attributes used by CLADE for the annotation of full length protein sequences. In MetaCLADE, these characteristics which produce a fragmentation of data, are taken into account during the learning and classifying step. In order to populate its training sets for the learning step, MetaCLADE tests exhaustively all probabilistic models in the library through the generation of a huge amount of artificial sequences. This learning step provides confidence thresholds for each domain in the library through a bi-dimensional threshold for probabilistic models hits, using the scores of the hits and the relative strength of the hit with respect to its length. This bi-dimensional threshold provides the tool with a capacity to recover significant signals from short fragments of protein domains. In addition, metagenomic datasets may contain millions of sequences. In order to manage this amount of data, MetaCLADE was developed to be used in high performance computing (HPC) machines with several cores. MetaCLADE has been tested in different metagenomics and metatranscriptomics data sets containing potentially very different species, it provides evidence of a significant increase in the quantity and the quality of the predictions compared against current methodologies for protein domain detection. In order to assess the quality of MetaCLADE annotations, a simulated metagenomics data set was constructed. Expert curated annotations for some of the simulated sequences were compared against those produced by MetaCLADE. Results from this simulation emphasize the very low rate of erroneous annotations (<1%) produced by our tool. MetaCLADE outperforms other methodologies by annotating more sequences, with more reliability in its predictions.

As mentioned before, the scope of this project is to better understand the composition of genes in microbial communities in order to depict a functional landscape of the most relevant processes that are carried out by its members. For this reason, MetaCLADE includes an automatic generation of functional profiles. This profiles characterized by histograms, are derived from the mapping of protein domain annotations to its functional annotations. These profiles enable biologists to analyze and compare functional classes which describe the metabolic preferences of the community and to observe which are the domains involved in those activities.

Through the chapters of this manuscript, we will firstly describe the biological and technical aspects that are closely related to the genomic analysis of microbial communities in order to allow the reader to completely understand the context of our work. Secondly, we will justify the redesign of CLADE approach and we will carefully describe MetaCLADE methodology. Finally, we will see the perspectives in future work that may be included in order to refine our tool. It is important to mention that this work has given rise to two publications: MetaCLADE: a multi-source annotation method for metagenomic and metatranscriptomic sequences. Ari Ugarte, Juliana Bernardes, Alessandra Carbone (2016,

Submitted)

Marine diatoms sense and respond to microscale turbulence in non-limiting nutrient conditions. Alberto Amato, Gianluca Dell'Aquila, Francesco Musacchia, Ari Ugarte, Nicolas Maille, Alessandra Carbone, Maurizio Ribera d'Alcalà, Remo Sanges, Daniele Iudicone, Maria I. Ferrante (Scientific Reports, 2016 under review)

The first publication describes in detail the MetaCLADE methodology while the second shows the usefulness of MetaCLADE when used to annotate sequences from several communities of micro-marine organisms under different micro-scale turbulent conditions.

Chapter 1

Preliminary Notions

Micro-organisms

Since the discover of micro-organisms scientists gradually became fascinated by the incredible diversity of these organisms. Adapting to environments where no life form is expected to thrive. They are able to process a great variety of nutrients, minerals and other environmental resources such as light. Moreover, they can transform them into other products that can be beneficial to other organisms in the same environment. The estimated number of these organisms is 4-6 × 10 30 making them the most abundant bio-mass on earth [START_REF] Whitman | Prokaryotes: the unseen majority[END_REF]. They have colonized every possible environment including some extreme environments where they are the only living beings that are capable to live in such conditions. One feature that help them to achieve this degree of colonization is their capability to create complex communities where the exchange of molecules is very important. The most outstanding exchanges may be considered the genetic information related to. Micro-organisms are able to import/export genomic sequences to produce many types of gene products [START_REF] Soucy | Horizontal gene transfer: building the web of life[END_REF]. This feature allows them to increase its genetic repertoire gaining new functions such as the resistance to something harmful or the ability to get energy from different sources. They are the base of the nutrients cycle in all environments, they help plants and fungus to grow through complexes exchanges in the rhizosphere [START_REF] Breidenbach | Microbial Community Structure in the Rhizosphere of Rice Plants[END_REF], they keep nutrients to make them available to other species in the oceans by being an important source of food [START_REF] Armbrust | The life of diatoms in the world's oceans[END_REF], they live in the large majority of species helping or transgressing them [START_REF] Soares | The Iron age of host-microbe interactions[END_REF]. Human activity is very related to micro-organisms since they affect our food and water sources [START_REF] Gray | Biology of Wastewater Treatment[END_REF], and they affect many of our biological processes [START_REF] Soares | The Iron age of host-microbe interactions[END_REF]. Taking all these factors into account, we can expect a great impact on the ecosystems and the living beings from the absence/presence of these organisms. It is, hence, an important task to unveil and understand as much as possible the incredible diversity in every sense that these organisms offer.

The origin of micro-organisms is still unknown. Nevertheless, their antiquity can be inferred from geological evidence. Rocks dated 3.4 billion years old, shows fossilized domelike structures, from a type of microbial sea community that became lithified (the process in which sediments compact under pressure). These structures, called stromatolites, were produced by anaerobic bacteria (primary cyanobacteria) and they are responsible for the introduction of oxygen into earth atmosphere [START_REF] Allwood | Stromatolite reef from the Early Archaean era of Australia[END_REF]. Since then, evolution has diversified their metabolism resulting in the incredibly microbial diversity we observe today.

Microbes includes mostly unicellular species from different clades but also animals have been included in this classification [START_REF] Baqai | Introduction to the Rotifera[END_REF]. Micro-organisms are very diverse and include all bacteria, archaea and most protozoa. Fungi and algae are also represented by some ecological important species such as diatoms, the most common type of phytoplankton. The most abundant group of micro-organisms are Prokaryotes [START_REF] Whitman | Prokaryotes: the unseen majority[END_REF], organisms without nucleus and they are classified in two domains: bacteria and archaea. Earlier classifications of these organisms were realized using different chemical and physical characteristics. Since the advances on the molecular era, the scientists have been able to classify them regarding its genomic sequences.

The most common method (not only for microbes but for all living organisms) is the analysis of the highly conserved sequences of ribosomal RNA (rRNA), which is a component of the ribosome, a complex found in all living cells which is essential for the protein synthesis [START_REF] Pace | A molecular view of microbial diversity and the biosphere[END_REF].

The comparison of rRNA sequences allows to unveil the evolutionary relationships among organisms. In particular, the 16S rRNA which is a component of the 30S small subunit of prokaryotic ribosomes and the 18S rRNA its eukaryotic homologous. The number of 16S rRNA sequences in current databases exceeds 4 million [START_REF] Yarza | Uniting the classification of cultured and uncultured bacteria and archaea using 16S rRNA gene sequences[END_REF] and it will continue increasing as sequences coming from uncultivated organisms will continue to be added. Figure 1 shows a phylogenetic tree constructed using 16S rRNA and 18S rRNA sequences of 64 different species. Micro-organisms are also classified by their morphology, genome size, relevance to human diseases and metabolic preferences such as their nutrient and energy sources.

Their roles in nature suggest complexes relationships as they are the primary decomposers of material. In soil they are the base of the nitrogen cycle and they play an Micro-organism do not only have a great impact on earth ecosystems, they are also very abundant inside other living beings, they stablish endosymbiotic or pathogenic relationships. It is estimated that the total amount of micro-organisms in a single human being is 7 × 10 12 , mainly at the colon (300 × 10 9 /g) [START_REF] Whitman | Prokaryotes: the unseen majority[END_REF]. The genomic material found in human gut microbiome is 100 times the human genome [START_REF] Foldout | The Inner Tube of Life[END_REF], therefore most of the genomic material found in humans come from their microbes. The estimation of microbes in all human beings is 4 × 10 23 , however, the estimation for aquatic (1 × 10 29 ) or subsurface (3.8 × 10 30 ) abundance is much more higher [START_REF] Whitman | Prokaryotes: the unseen majority[END_REF]. Micro-organisms grow and live in multi-species assemblies which can be seen as a general definition of an ecological community. It is important to unveil how these biological assemblages are structured and which are their functional activities and interactions.

Communities are never static, therefore it is also important to understand how these structures and functional interactions changes in space and time. [START_REF] Konopka | What is microbial community ecology[quest][END_REF].

The number of different prokaryote species is unknown but some estimations conclude that there are between 10 7 and 10 9 [START_REF] Curtis | Estimating prokaryotic diversity and its limits[END_REF][START_REF] Dykhuizen | Santa Rosalia revisited: why are there so many species of bacteria?[END_REF]. The number of different microbial species in a community depends on the environment and different estimations for different environments has been proposed. To make a concrete example, the estimations of the number of different species per gram of soil vary from 2 000 up to 830 000. [START_REF] Gans | Computational improvements reveal great bacterial diversity and high metal toxicity in soil[END_REF][START_REF] Schloss | Toward a census of bacteria in soil[END_REF]. Even with the lower estimations, taking into account that the average size of a soil bacterial genome is 6.3 MB, [START_REF] Trevors | One gram of soil: a microbial biochemical gene library[END_REF] we can expect therefore to have at least 12.6 GB of DNA information per gram of soil which is almost 4 times the size of the human genome (3.3 GB). It is clear that they have a great genomic potential that is reflected in their ability to adapt to different conditions. Two decades have passed since the first bacterial genome from bacterium Haemophilus influenzae Rd was published, [START_REF] Fleischmann | Whole-genome random sequencing and assembly of Haemophilus influenzae Rd[END_REF] since then, many projects for complete bacteria sequencing have been realized and today it is estimated that there are more than 30 000 complete sequenced prokaryotic genomes available in sequence databases [START_REF] Land | Insights from 20 years of bacterial genome sequencing[END_REF]. This has led to great discoveries from genomic analysis in the latest years and a great number of databases containing genomic sequences has been established. This has allowed scientists to better understand the microorganisms functional characteristics as individual but there is still a lack in the comprehension of micro-organism working as a community. In relation with the total estimation of 10 7 and 10 9 [START_REF] Curtis | Estimating prokaryotic diversity and its limits[END_REF][START_REF] Dykhuizen | Santa Rosalia revisited: why are there so many species of bacteria?[END_REF] different prokaryotic species, we are still far away from the completion of all existent genomic sequences. One of the main limitations to obtain its sequences is the unavailability to grow many of these organisms under laboratory conditions. Recent advances in next generation sequencing has allowed researchers to explore microbe diversity directly from nature, genomics specialists are now retrieving information from meta-omics data sets and computational methods are used to infer the taxonomic composition of a community and to describe its functional activity. This information is also important for ecologists as it can be related to different environmental conditions, we are therefore in a context where genomics, ecology and micro-biology intersects.

Metabolic Pathways

Metabolism refers to all chemical reactions that occur in a cell to maintain it alive.

Metabolism can be divided in two types of reactions:

1. Catabolic reactions where complexes macro-molecules are broken down into simple molecules creating energy 2. Anabolic reactions where synthesis of complex molecules is made, requiring energy Catabolic reactions produce Adenosine Triphosphate ATP, the molecule that stores energy in cells. ATP is composed of an adenosine nucleoside and three phosphates. The energy in ATP is stored in the bonds between the phosphate groups, when these bonds are broken usable energy is released leaving a molecule of adenosine diphosphate (ADP or APP). ATP is continually being formed from lower-energy molecules of APP and AMP. The biosynthesis of ATP is achieved throughout processes such as substrate-level phosphorylation, oxidative phosphorylation, and photophosphorylation, all of which facilitating the addition of a phosphate group to an APP molecule [START_REF] Alberts | Molecular Biology of the Cell[END_REF].

Metabolic pathways refer to the different pathways or series of reactions by which essentially molecules are formed or broken down. Anabolic and metabolic reactions in metabolic pathways are the opposites of each other. These pathways are driven by a great diversity of enzymes which catalyze chemical reactions. Enzymes are glob-like proteins or complexes which sizes ranges between ~ 60 and ~2500 residues [START_REF] Chen | 4-Oxalocrotonate tautomerase, an enzyme composed of 62 amino acid residues per monomer[END_REF][START_REF] Smith | The animal fatty acid synthase: one gene, one polypeptide, seven enzymes[END_REF]. Each step in a metabolic pathway requires a different enzyme. Metabolic processes produce metabolites which are intermediate and products of enzymatic reactions. The Glycolisis metabolic pathway is an example of a catabolic metabolic pathway. This pathway which is present in all living organisms converts glucose into pyruvate. This process uses 2 ATPs during the first steps of the path but it generates 4 ATPs, so the net gain is of 2 ATPs. This pathway can also start from other input points such as F6P (fructose). It also produces other useful metabolites such as dihydroxyacetone phosphate (DHAP) which can be used to form fat. Microbes are mainly lithotrophic organisms, they capture their energy from inorganic chemicals in the environment and they are the most abundant form of life. Lithotrophic organisms are classified by its cellular respiration. Aerobic respiration requires oxygen to generate ATP while anaerobic respiration uses other compounds [START_REF] Alberts | Molecular Biology of the Cell[END_REF].

Metabolic pathway components can be detected on microbial communities through the identification of enzymes. Metagenomics and metatrascriptomics studies can provide evidence of the most active pathways in a community. Studies of the alteration in metabolic pathways of the human microbiota under specific conditions have given evidence that some pathways can became more or less actives affecting human health [START_REF] Zheng | Metagenomic sequencing reveals altered metabolic pathways in the oral microbiota of sailors during a long sea voyage[END_REF]. Resources providing a well curated mapping between current domain databases and metabolic pathway databases will improve considerably the characterization of the complexes activities occurring in communities.

Next Generation Sequencing (NGS)

We have shown the importance of environmental studies by revealing its genomic potential and the dynamics it can produce. Technologies producing DNA sequencing are the ones who had driven these studies to produce significant results. Manufacturers provide data with distinct characteristics which are suitable to specific analyses (e.g. genome assembly, gene expression, structural variations) [START_REF] Liu | Comparison of Next-Generation Sequencing Systems[END_REF]. Through this section of the manuscript, we will describe current NGS technologies in order to understand advantages and limitations of each one of them.

Sequencing technologies have undergone great changes since Frederick Sanger developed the Sanger sequencer [START_REF] Sanger | DNA sequencing with chain-terminating inhibitors[END_REF], a technology based on the chain-termination method.

Sanger sequencer is known as the first commercial sequencer as the method had low radioactivity and provided high efficiency. Since then, different approaches have emerged and automatic sequencing instruments have been produced. Software embedded in these machines became the standard among constructors as it improves the quality of delivered data. At the beginning, the great majority of sequencer were based in Sanger technologies. It is only in recent years that completely different approaches to Sanger starts to emerge, these technologies represent the second sequencing generation [START_REF] Liu | Comparison of Next-Generation Sequencing Systems[END_REF]. As technologies evolve, sequencing became cheaper and the data they can produce became massive this is shown in the cost per nucleotide since second generation sequencing appeared (Figure 3). Manufacturers have slightly modified its base approaches in order to provide data with different characteristics, depending on the type of the study, scientists may prefer read length or to priorize throughput. In general, larger reads imply lower throughput.

DNA sequencing preparation

Once DNA has been extracted it is broken into smaller pieces. Fragmentation of DNA is an early step in next generation sequencing workflows, as well in the creation of DNA inserts for expression libraries. Methods of DNA fragmentation include: Acoustic shearing, nebulization and sonication [START_REF] Liu | Comparison of Next-Generation Sequencing Systems[END_REF].

After DNA fragmentation, artificial DNA oligonucleotides (oligos) are added to the denatured fragments in vitro. One kind of these sequences are the PCR primers, which are used to prime DNA replication reactions. In addition to PCR primers used in the amplification step, different kind of adapters can be added to DNA molecules in order to tag with a specific sequence some of the molecules. In paired end reads. Once adapters have been added, DNA fragments are amplified. There are two different technics for DNA amplification: Bead Based Emulsion PCR and Solid phase bridge amplification [START_REF] Liu | Comparison of Next-Generation Sequencing Systems[END_REF].

In bead based emulsion PCR (polymerase chain reaction) [START_REF] Williams | Amplification of complex gene libraries by emulsion PCR[END_REF], DNA adapters in fragments are complementary to fixed oligonucleotides in the capture beads. Both along with enzyme reagent in a water mixture are put into small cylindrical containers containing a synthetic oil. Then, vigorous shaking makes water to create droplets around the beads.

Typically, droplets should contain only one DNA fragment. Complexes in the droplets amplify isolated DNA fragments into millions of copies (up to 10M). After DNA amplification, beads are cleaned and those beads that do not hold DNA are eliminated. Beads containing more than one DNA fragment will be filtered during the sequencing step (See figure 4A).

In Solid phase bridge amplification [START_REF] Kawashima | Patent: Method of nucleic acid amplification[END_REF], the process where each fragment is isothermically amplified is called clustering. In the clustering step. The flow cell is a glass slide with lanes. Each lane is a channel coated with a lawn composed of two types of oligos.

Hybridization is enabled by the first of the two types of oligos on the surface. This oligo is complementary to the adapter region on one of the fragments strands. A polymerase creates a compliment at the hybridized fragment, the double stranded molecule is denatured and the original template is washed away. The strands are clonally amplified to bridge amplification.

In this process the strands pull over and the adapter region hybridizes to the second type of oligo on the flow cell. Polymerases generate the complementary strand, forming a double stranded bridge. This bridge is denatured, resulting in 2 single-stranded copies of the molecule that are tethered to the flow cell. The process is then repeated over and over and occurs simultaneously to millions of clusters resulting in clonal amplification up all the fragments (See figure 4B).

Systems

MinION (Oxford Nanopore) [START_REF]Oxford Nanopore Presents Details on New High-throughput Sequencer, Improvements to MinIon[END_REF]. The MinIon device works with a proprietary nanopore (a nano scale hole). The nanopore is inserted into an electric resistant membrane created from synthetic polymers, a potential (ion current) is applied across the membrane resulting in a current flowing through the aperture (hole) of the nanopore. Single molecules that enter the nanopore cause characteristic disruptions in the current. This event is the nanopore signal. By measuring that signal, the molecule can be identified (for example the four standar DNA bases). This approche, is suitable for single molecule analysis such as DNA, RNA and proteins. Inside the minion device there is a flow cell containing an array of microscaffolds (a sensor array which is a collection of electrodes and micro supports). Each microscaffold is connected to a channel on the Application-Specific Integrated Circuits providing real-time data.

Once a DNA strand has been completely sequenced, a new one will enter the nanopore. This means that the experiment can go on as long as the user required it. The number of nanopores PicoTiterPlate is loaded with capture beads containing the amplified DNA fragments, it is placed into the 454 sequencing instrument. During the cyclic flowing process, each of the deoxynucleotides (dNTPs or bases) are flown sequentially and in the same order. Up to 400 cycles are accomplished per run. When the flowed nucleotide is complementary to the template strand, the polymerase extends the primer and pauses. Pyrosequencing uses luciferase (generic oxidative enzymes that produce bioluminescence) to generate a chemi-luminiscent signal after the base has been incorporated. This signal is detected by a CCD camera included in the instrument. Then current dNTPs are degradated before the inclusion of the next dNTPs. The intensity of the light generated is proportionally dependent to the number of consecutive incorporated nucleotides. The order and the intensity of the light signal are recorded as flowgrams and interpreted as the resulting DNA sequence by the system's software producing millions of sequenced bases per hour.

NGS Performance

Next Generation Sequencing technologies are usually benchmarked by their limitations in term of the type and the size of the sequences, the amount or throughput of the sequences, the accuracy and the error profile. Other variables to take into account are the run time and naturally the cost of the equipment. Translation [START_REF] Alberts | Molecular Biology of the Cell[END_REF] is the synthesis of proteins from mRNA. The ribosome, a complex composed by several proteins and RNAs links amino acids together following the sequence in the mRNA. The mRNA is read by sequences of three nucleotides called codons. Codons determine which amino acid will be incorporated by the ribosome into the protein.

As mentioned before, amino acids form long chains linked by peptide-bonds. The amino acid sequence of a protein is called its primary structure. Secondary structure, refers to hydrogen bonding dependent regular local sub-structures which are formed within a protein.

The most common secondary structures are alpha helices and beta sheets; different groups of amino acids may tend to be structured in one or the other conformation. The tertiary structure is the final three-dimensional shape of a folded protein. The quaternary structure of proteins is the 3d structured formed by the assembly of many proteins in larger complexes. When proteins bind with other molecules, their conformation may change [START_REF] Pauling | The structure of proteins; two hydrogenbonded helical configurations of the polypeptide chain[END_REF].

Proteins do most of the cellular task in a cell, some types include enzymes which catalyze the biochemical reactions that occur in cells, structural proteins which maintain cell shape, transport/storage proteins which bind and carry atoms and small molecules within cells and messenger/receiver which produces/reacts to different metabolic signals. The function of a protein/complex depends on its tertiary or quaternary structure which is in turn, depends on its primary structure or amino acid sequence [START_REF] Pauling | The structure of proteins; two hydrogenbonded helical configurations of the polypeptide chain[END_REF]. As mentioned before, proteins sharing significant regions of its amino acid sequence are thought to share also the same function. This similarity indicates a common evolution origin. Proteins are also composed of functional/structural sub-units called domains. Domains are conserved regions in proteins which can evolve and exist independently of the rest of the proteins. These structural sub-units which mean size is ~100 amino acids [START_REF] Richardson | The anatomy and taxonomy of protein structure[END_REF], can have functional roles such as the assisting in protein-protein interactions that may be useful in proteins having different biological functions [START_REF] Janin | Structural domains in proteins and their role in the dynamics of protein function[END_REF]. Domains are therefore combined in different ways in order to exploit its individual functional characteristics by generating an architecture which defines the entire protein function [START_REF] Bernardes | A multi-objective optimization approach accurately resolves protein domain architectures[END_REF]. Domains may give therefore, insights of the biological processes in which a protein is involved. Moreover, some domains are present only in certain types of proteins which turns them into fingerprints which can be used to infer the role of a protein. Highly conserved domain regions may characterize functional sites which are the region of domains (i.e. proteins) that are involved in a specific function such as ligandbinding or protein-protein interaction. In the case of enzymes, these regions are the active sites, the regions where molecules undergo chemical reaction.

exchanges data with 2 other databases: The European Nucleotide Archive [START_REF] Leinonen | The European Nucleotide Archive[END_REF] (ENA) part of the EMBL bank [START_REF] Stoesser | The EMBL Nucleotide Sequence Database[END_REF] and the DNA Databank of Japan (DDBJ) [START_REF] Kodama | The DNA Data Bank of Japan launches a new resource, the DDBJ Omics Archive of functional genomics experiments[END_REF]. The daily exchanges between the databases ensures uniform and comprehensive sequence information among them.

The three databases are part of the International Nucleotide Sequence Database Collaboration [START_REF] Cochrane | The International Nucleotide Sequence Database Collaboration[END_REF] (INSDC). More than 260 000 species are represented in the GenBank sequences.

Sequences submitted to the GenBank are functional divided regarding to the sequencing strategies used to obtained them such as: High-throughput genomic (HTG) and high-throughput cDNA (HTC) sequences, expressed sequence tags (ESTs), Whole-Genome Shotgun Sequences (WGS), Environmental sample sequences (ENV) and transcriptome Shotgun Assembly (TSA) sequences. GeneBank also classify its sequences in 12 taxonomic divisions. GenPept [START_REF] Benson | GenBank[END_REF] is the database that contains the translated coding sequences from the GenBank database.

RefSeq [START_REF] O'leary | Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation[END_REF] Protein Data Bank (PDB) [START_REF] Berman | The Protein Data Bank[END_REF] is a database of structural information about macromolecules such as proteins and nucleic acids. The data is obtained by several experimental approaches such as X-ray crystallography, NMR spectroscopy, or, increasingly, cryo-electron microscopy. PDB have more than 115 000 structures in its database.

UniProt [START_REF]UniProt: a hub for protein information[END_REF] is a resource providing high-quality protein sequences and functional annotations. It was created by the the Protein Information Resource (PIR) [START_REF] Wu | The Protein Information Resource: an integrated public resource of functional annotation of proteins[END_REF], the EBI [START_REF] Brooksbank | The European Bioinformatics Institute's data resources[END_REF] (European Bioinformatics Institute) and the SIB [START_REF] Stockinger | Fifteen years SIB Swiss Institute of Bioinformatics: life science databases, tools and support[END_REF] (Swiss Institute of Bioinformatics) through a grant from the National Institute of Health [START_REF] Collins | Exceptional opportunities in medical science: a view from the National Institutes of Health[END_REF] (NIH). UniProt derives its protein sequences from different external sources such as GenBank, RefSeq, PDB and other resources.

All sequences are included in the UniParc sequence archive, which is a database containing most of the publicy available amino acid sequences. UniParc is a non-redundant database, providing a unique identifier (UPI). Sequences composing the UniParc database are included in the UniProt Knowledgebase (UniProtKB) which is a collection of information on proteins.

The UniProtKB is divided in two sections. The unreviewed part: TrEMBL (translated EMBL-Bank), which contains computationally generated annotations and the reviewed part (Swiss-Prot), which contains manually annotated and reviewed data. Swiss-Prot is considered to be the golden standard of protein annotations, and almost every software attempting to annotate proteins should use this data base as reference to evaluate its performance. UniProtKB annotations includes references to other sequence, family and domain databases. Based on the sequences included in the knowledgebase, UniProt provides a set of taxon specific proteomes and a set of sequence clusters at different thresholds called UniRef. A great amount of supporting data is also available through its portal.

The NR database [START_REF]Database resources of the National Center for Biotechnology Information[END_REF] is a protein database maintained by the NCBI. NR stands for non-redundant so every pair of sequences in the database are non-identical. The NR database integrates sequences from a great variety of sources such as UniProt, PDB and RefSeq. It provides a great coverage of the entire set of protein sequences publicy available. It can be compared in coverage to UniParc database (see UniProt above).

Homology detection software

The BLASTp (Basic Local Alignment Search Tool) [START_REF] Altschul | Basic local alignment search tool[END_REF] algorithm search for similar sequences to a sequence query in a protein database. The search results consist of the hits satisfying a threshold, along with a scoring scheme to describe significance of the match. The BLASTp algorithm can be described in three phases (Figure 10)

1. For protein searches, BLASTp compiles a preliminary list of pairwise alignments called word pairs.

2. The algorithm scans a database for word pairs that meet some threshold score T.

When this occurs, such hits are extended using ungapped and gapped alignments.

BLAST extends the word pairs to find those that surpass a cutoff score S, at which point those hits will be reported to the user. Scores are calculated from scoring matrices (such as BLOSUM62 [START_REF] Henikoff | Amino acid substitution matrices from protein blocks[END_REF]) along with gap penalties.

3. A trace-back procedure is performed in which the locations of insertions, deletions and mismatches are assigned.

In the first phase, the BLASTp algorithm parses the query sequence and obtain a list of words of a fixed length w. For each of these words, BLASTp creates a new set of words matching it. A threshold value T is established for the score of aligned words. Those words having a score S >= T are kept and used to find matches in the database matches. For protein searches the word size typically has a default value of 3 (8000 possible words).

In the second phase, the BLASTp algorithm scans a database for matches against the pre-computed word list. Then it searches for two separate word pairs within a certain distance D from each other. It then generates an ungapped extension of these hits. BLASTp then extend hits to produce high-scoring segment pairs (HSPs). If HSPs scores are above a threshold T', a gapped extension is done. The extension process is terminated when a score falls below a cutoff.

In the third phase, a trace-back is performed in which the locations of insertions, deletions, and matches are assigned. Composition-based statistics are applied. As a heuristic algorithm, BLASTp offer both speed and sensitivity.

Figure 10 The BLASTp algorithm In phase 1 a list of words derived from the input with a score above a threshold T is generated. In phase 2 the input database is scanned in order to generate high scoring pairs by performing gap and ungapped extensions of the initial hits. In phase 3 scores are calculated and the final gapped alignment is generated. http://www.bioinfbook.org/wiley/3e/chapter4/JWST573_Pevsner-chapter-

04-figures.ppt

When two sequences are aligned, they produce a score. We are therefore, interested in determine if this score is likely to have occurred by chance. The problem is to find the probability of a best score above a threshold s using BLASTP. For the comparison of a query sequence to a database of random sequences of uniform length, the highest scores can be plotted and shown to have the shape of an extreme value distribution (Figure 11). Therefore, for two random sequences of length and , the cumulative distribution function of an extreme value distribution of scores is defined:

< = exp ( -( ) ) (1)
For ungapped alignments, the parameter is dependent on the lengths of the sequences being compared and is defined:

= (2)
Where K is a constant to scale the search space and is a constant to scale the scoring. it is also the unique solution to the equation:

, = 1 (3)
Where is the frequency of amino acid in the query and is the frequency of the amino acid in the subject and is the score of aligning an , pair.

When combining equations 1 and 2 we have:

> = 1 -exp ( - ) (4)
The term for the expected number of ungapped alignments with a score of at least by chance is:

= (5)
For gapped alignments however, λ, cannot be calculated analytically, but they can be estimated from pre-computed values.

Raw scores are calculated from the scoring system and gap penalty parameters. The bit score S′ is calculated from the raw score by normalizing with the statistical variables that define a given scoring system. Raw scores do not have units and have little meaning alone. Bit scores account for the scoring system that was used and allow scores to be compared between different database searches, even if different scoring matrices are employed.

=

The e-value corresponding to is:

= ×2 (7)
PSI-BLAST (Position-Specific Iterative Basic Local Alignment Search Tool) [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF] is a version of the BLASTp algorithm that constructs PSSMs (Position Specific Score Matrix) in each iteration. A PSSM is a type of scoring matrix constructed from a multiple sequence alignment and can be seen as a model of evolution as it describes the amino acids substitutions Each domain family is composed by a set of sequences called the full set. A sub-set of representative sequences from the full set called the seed set is used to generate a multiple sequence alignment and to construct a sequence global consensus model using a pHMM. Pfam uses HMMer3 to find hits against it pHMM library. For each family in the database they set a bit score gathering (GA) threshold by hand, such that all sequences scoring at or above this threshold appear in the full alignment.

Gene3d [START_REF] Lees | Gene3D: Multi-domain annotations for protein sequence and comparative genome analysis[END_REF] is a database of protein domains that are directly mapped from the structures of the CATH [START_REF] Sillitoe | CATH: comprehensive structural and functional annotations for genome sequences[END_REF] (Class, Architecture, Topology Homology) database or predicted using a library of pHMMs representing each of the 2 737 CATH superfamilies. he largest domain superfamilies contain sequences that have greatly diverged in molecular function.

Recently work has been carried out to help improve the functional purity of domain assignments by dividing the domain superfamilies into smaller functionally coherent groups termed Functional Families or FunFams (3). These FunFams greatly improve the ability to interpret the functions of an experimentally uncharacterized protein based on its domain assignments.

FunFams are constructed using an in-house protocol [START_REF] Janin | Structural domains in proteins and their role in the dynamics of protein function[END_REF], which involves profile-profile-based clustering of domain sequences in each superfamily to identify functional families. We have recently improved the speed and assignment quality of the method so that it is possible to identify FunFams in all Gene3D superfamilies. After family identification, an HMM is built for each FunFam using HMMER, in conjunction with a model-specific bitscore threshold based on the score attained by the most remote member sequence. Gene3d contains 53,479,436 nonstructural protein domain entries.

TIGRFam [START_REF] Haft | TIGRFAMs and Genome Properties in 2013[END_REF] is a database of curated protein families including sequence alignments, pHMMs, Gene Ontology assignments (see Gene Ontology Section) and pointers to Pfam and InterPro databases. TIGRFam makes special restriction in thresholds and sequences used to construct pHMM models. They exploit information from sequence similarity, evidence of function from literature, phylogenetics, species specific metabolic contexts in order to select sequences to be included in the final alignment to create probabilistic models. It introduces the concept of equivalogs, which are proteins that have conserved the same function since their last common ancestor. In TIGRFam a super family is a set of proteins having global homology over their whole length and it members may vary greatly in function. Subfamilies are clades within super families, as they share specific functions or properties. TIGRFAM contains models of full-length proteins and shorter regions at the levels of superfamilies, subfamilies and equivalogs. The majority of the 4284 models are of type equivalogs, giving TIGRFam models great accuracy at a functional level. TIGRFam is focused in a broad coverage of prokaryotic proteins but also include some families of interest from eukaryotes and plant. TIGRFam as other pHMM databases is powered by HMMer3.

PANTHER [START_REF] Mi | PANTHER in 2013: modeling the evolution of gene function, and other gene attributes, in the context of phylogenetic trees[END_REF] (protein annotation through evolutionary relationship) is a large database covering ~90 % of mammalian protein-coding genes. Panther database is composed of protein sequences from more than 100 complete genomes taken from the UniProt Reference Proteomes dataset (a data set generated manually and algorithmically in order to span as much as possible the entire UniProt dataset). Sequences are clustered into families, and for each family a phylogenetic tree is build in order to detect sub-families by close similarity.

Evolutionary relationships (such as speciation or duplication) are inferred by PANTHER for all sequences in the tree. The phylogenetic data generated for families in PANTHER is used to generate functional models matching experimental data and characterising the gain and loss of functions in specific parts of the tree. These functional models can be used to identify functions of uncharacterised sequences. Each family and sub-family is represented by a pHMM from the alignment of member sequences. PANTHER contains more than 12 000 families.

SMART [START_REF] Schultz | SMART, a simple modular architecture research tool: Identification of signaling domains[END_REF]

Functional Annotation

The Gene Ontology (GO) [START_REF]Gene Ontology Consortium: going forward[END_REF] project consists of three structured ontologies to describe gene products. The three ontologies describe in a species independent manner, three different domains of gene products: the biological processes, the cellular components and the molecular functions. Biological processes are series of steps that are done by one or more organized assemblies of molecular functions. Molecular function is related to the activities or jobs that a gene product does such as binding, transporting, etc. Cellular component describes locations, at the levels of subcellular structures and macromolecular complexes but it does not Chapter 2

CLADE system and performance on simulated and real metagenomics datasets

CLADE

CLADE [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF] (CLoser sequences for Annotations Directed by Evolution) is a domain annotation tool for entire proteins which combines predictions from almost 2.5 M of probabilistic models. The outcomes obtained from models are processed and transformed into features used to train a meta-classifier (Support Vector Machine [START_REF] Scholkopf | Advances in kernel methods: support vector learning[END_REF]), that assigns a confidence score to each domain prediction. Based on this score and domain co-occurrence, CLADE finds the most probable architecture for each protein sequence using the DAMA [START_REF] Bernardes | A multi-objective optimization approach accurately resolves protein domain architectures[END_REF] software (Domain Annotation by a Multi-Objective Approach).

CLADE Extended Library.

As seen in chapter 1, a domain can be characterized by a probabilistic model. These probabilistic models are build from the multiple sequence alignment of a set of homologous sequences. Current annotation systems are based mainly in two approaches:

1. Given a set of homologous sequences S i associated to a particular domain family D i . A subset S i* of representative sequences is used to construct a Sequence Consensus Model (SCM) 2. Given a set of homologous sequences S i associated to a particular domain D i .

Closely related subsets {S i0* … S iN* } of S i are used to construct several models characterizing a domain. We shall speak of a multi-source strategy.

SCMs are used in the Pfam database, the biggest collection of known protein domains with more than 15 000 different domains. Methodologies using the second approach are rather restricted to a small set of domain families and they also produce a small number of probabilistic models for each domain, for instance, SUPERFAMILY [START_REF] Wilson | The SUPERFAMILY database in 2007: families and functions[END_REF] However, current methodologies cover a small fraction of the entire Pfam database and moreover, the number of probabilistic models characterizing a single domain is rather small, having on average less than 10 models per domain [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF].

CLADE tackles this problem through a significant extension of the Pfam database for which each domain is now characterized in the CLADE library by: the global consensus pHMM model generated by Pfam and up to 350 probabilistic models named Clade Center Models (CCMs). CCMs models are not build from internal partitions of a set of homologs (FULL) associated to each domain in the Pfam database. Instead, for each domain, CLADE defines a subset of sequences of the Pfam FULL set. Then, for each of those sequences it finds closely related sequences in the NR database (Figure 16) and constructs a Clade Center Model (CCM) model using PSI-BLAST with a threshold of 1e-3 which is reported to avoid corrupt PSSMs [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF]. The criteria to select the subset of sequences used as queries to construct CCM models is based on their taxonomy. CLADE selects sequences from species that are uniformly spread in the phylogenetic tree of life. The selection guarantees that species belong to different phylogenetic clades and that the phylogenetic tree is well represented. Figure 17 shows the distribution of the sequences grouped by phylogenetic source that were used to construct the CCMs. CLADE library is composed of 2 389 235 CCMs and 14 831 SCMs (one for each

The CLADE pipeline:

The CLADE pipeline is divided in three different steps:

1.-Finding hits. Each model in the CLADE library is used as query to find hits in a database of protein sequences. HMMer is used to find matches for the SCMs and PSI-BLAST is used to find matches for the CCMs. From these 5 attributes, a meta classifier (SVM) has been trained to discriminate between the true and false predictions (Figure 18). More precisely, a one-vs-rest SVM i . The positive training set for each domain i is defined as the S set excluding the sequences used for the construction of models. The negative training set, is composed from random shuffle of the 2-mer sequences in each set S i , where i = 1…14 831 (the number of Pfam27 domains). The method Platt [] was used to map the output of the SVM to a probability a posteriori by using a logistic function. For each domain D i , a specific SVM probability cut-off T Di was defined. data is much more difficult to accomplish due to the loss of domain signal.

One method to overcome fragmentation in short reads is to assemble reads into bigger contigs, however the succeed of such task is closely related to the abundance and taxon distribution in the sample and also to the phylogenetic distance between them [START_REF] Charuvaka | Evaluation of short read metagenomic assembly[END_REF]. We should remember that metagenomics samples can contain thousands of different species (See chapter 1). Therefore, we can expect to have reads belonging to poorly abundant species that cannot be assembled as a result of low coverage regions. Non-abundant species can contain also important functions for the community and identification of their genes is also important for the understanding of the community. Fragmentation caused by insertion/deletion errors can be avoided by using error models during the detection of CDS (coding sequences) in reads in order to correct them. As we saw in chapter 1, NGS technologies producing the longest reads have up to 80 percent of insertion/deletion errors.

Metagenomics studies are considered big data analyses, the amount of samples and replicates per study are increasing as a result of lower prices per run in NGS (see Chapter 1).

Coding sequences identification is very important in order to reduce the size of the data to be analyzed. As mentioned before, CLADE uses an extended probabilistic model library resulting in a higher run time during the search step. Therefore, it was important for us to evaluate available prediction tools which are designed to identify coding sequences in short reads. This, in order to reduce the number of sequences taken as input. We were also interested to analyze CLADE performance in metagenomics datasets, this was a crucial step in MetaCLADE development. It was important to assess not only whether the extended library of probabilistic models provide an increase in the signal of fragmented data, but also to assess the limitations of the learning strategy in CLADE when applied to metagenomics data.

Simulation

In order to provide evidence of the performance of the CLADE library, we generated a simulated metagenomic dataset. For this simulation, we produced short reads from different organisms for which its genomic sequences are available. This allows us, to retrieve information about the regions coding for proteins that are present in simulated reads. Hence, we can determine whether a read contains domain fragments or not. In order to evaluate if CLADE probabilistic model library provide a significant increase of domain hits in simulated reads, its results were compared against those produced by HMMer. The input for both methods was generated as follows:

1.-The dataset is composed of simulated reads from 56 sea microorganisms for which genomes were downloaded from the ebi genomes site [START_REF] Brooksbank | The European Bioinformatics Institute's data resources[END_REF] (see Appendix A). Only bacteria and archaea organisms were included in the simulation. 500 000 fragments from the entire genomic sequences were sampled from random positions using MetaSim, [START_REF] Richter | MetaSim: a sequencing simulator for genomics and metagenomics[END_REF] a sequencing simulation tool for genomics and metagenomics which provide a generator of random sequences from the input. Using this generator, we produced short DNA sequences sampled from the genomes of sea prokaryote organisms. These sequences referenced as clones were produced with a mean size of 1000 nt and a standard deviation of 100 nt (Figure 20). In order to sample most of the coding sequences present in the simulated dataset, species were set to be equally abundant. This assures that all regions of the genomes are sampled through small sequences.

2.

-FlowSim [START_REF] Balzer | Characteristics of 454 pyrosequencing data-enabling realistic simulation with flowsim[END_REF] a tool to emulate the sequencing of roche's 454 pyrosequencer was used to produce simulated reads of 454 titanium GS FLX. The 500 000 clones generated in the previous step were used as input to produce reads with size (mean size of 523 base pairs) and error rates as the 454 titanium technology (Figure 21). FlowSim can produce reads with insertion/deletion errors according to the 3 models of roche's 454 pyrosequencers that are available. In order to validate that errors were produce with similar characteristics to the 454 titanium technology, we aligned reads and clones using blast. In Figures 22,23 we can observe that the data produced has the desired characteristics having a mean error rate below 1% and sequencing errors which tend to occur at the end of the read. Figure 22 Error rate in simulated reads. Figure 23 Error position distribution in simulated reads.

3.-The third step consist in determine the coding sequence composition of each read.

To perform this task, start and end positions of coding sequences (CDS) from microbial genomes were compared against start and end positions from which clones and reads were generated. The overlapping (> 20 nucleotides) of these positions defines the content of CDS fragments in each read (yellow regions in simulated reads figure 24). By the end of the simulation, we obtained a set of 500 000 reads for which 681 702 CDS fragments have been mapped. On average, there are 1,37 CDS Fragments per read.

of sequences with a huge amount of spurious CDS below 150 nt. MetaGeneMark seems to produce a similar distribution to the one generated by the ORF in the range 150-200 nucleotides but it discriminates a lot of short ORFs. However, it still produces a lot of spurious predictions bellow the 100 nt. The CDS fragments dataset and the three pCDS datasets obtained by the gene detector tools were translated to amino acid sequences and taken as input to search domains using the Pfam database. HMMScan was set to report the best non-overlapping domain hits having an evalue less than 1e-3, the distribution of the e-values of identified domains in CDS fragments is

A first evaluation of CLADE and HMMScan in metagenomics simulated data

The 583 849 pCDS produced by FragGeneScan were annotated using CLADE annotation system and HMMScan using an e-value threshold of 1. Best hits in non-overlapping positions were reported from HMMScan predictions. Results at overlapping positions were then compared and agreements and disagreements were quantified. In total, HMMScan identifies 335 745 domains while CLADE identifies 434 522.

We count as agreement when two predictions overlap and both have the same domain or same pfam clan (see Chapter 1). From the 335 745 identified domains by HMMScan, CLADE agrees in 323 879 domains that were predicted with the same accession number or same clan at overlapping positions. We can observe that in general, the E-values associated to CLADE predictions are much more significant than the ones obtained with HMMer Figure 33.

This is due to predictions coming from CCM models which seems to amplify domains signals.

However, some domains were annotated with more significant E-values by HMMScan. This is due to the domain co-occurrance algorithm and the SVM assigning high scores to non significant predictions (see disagreement)

Disagreed predictions. We identified 11 866 domains for which CLADE disagrees with HMMscan annotation. CLADE misses 6,917 (Plotted over the x-axis for visualization in Figure 34) domains that were annotated by HMMScan. CLADE uses an algorithm (DAMA) that was designed to find the most likely domain architecture for entire CDS based in co-occurrence. When two non-significant predictions are co-occurrents they are preferred over single predictions having much more significant E-values. Due to this, DAMA may predict 2 or more non-significant co-occurrent domains instead of one with a much more significant evalue. The spurious short signals may result in false predictions. We identified 4,949 overlapping predictions for which CLADE found domains belonging to a different clan than the ones predicted by HMMScan. For most of them, CLADE found more significant hits, nevertheless this can't assure that CLADE predictions are correct. There are properties in hits other than the scores that should be taken into account, for instance the identity against the probabilistic models (see Chapter 3).

Finally, we observe an increase of 110 643 domains which are identified exclusively by CLADE. The increase in the quantity of identified domains is the result of the domain library extension. In addition, some of the E-values associated to these hits are below 1e-50 (a very significant value) are completely missed by the SCMs (pHMMs) models of the Pfam database.

This means that some signals completely disappear in SCMs and can be detected only by the CCMs models. However, almost 40 % of the 110 643 new predictions have an E-value greater than 1e-8. In Table 3, we can observe the percent new identified domains with CLADE compared to HMMer when annotating:

1. All the proteins of the genomes present in simulation.

The exact (without errors) CDS Fragments present in reads

3. Predicted CDS in reads.

We can observe that CLADE identifies 15% more domains in entire proteins (full length CDS in figure 9), and it identifies 17% more domains when only fragments of those sequences are analyzed (yellow fragments in Figure 29). The 33% increase which correspond to the new 110 643 identified domains when pCDS (green sequences in Figure 29) seems to be disproportioned with respect to the increase of fragmented CDS. One can conclude, that there is an important amount of spurious detections included in the 110 643 detections, most probably belonging to the 40% that are above the 1e-8 threshold. In the next chapter, we will see how MetaCLADE proposes a new strategy to discriminate predictions with low confidence scores.

very short reads will not be detected by the tool as the expected coverage will be poor. The same problem may arise in large reads/contigs containing many insertion/deletion errors, the domain signal can only be detected when the correct reading frame is being analyzed.

Metagenomic samples may contain a lot of unclassified species, for which evolutionary pathways may vary from the ones present on current databases. It is then possible, that only a few probabilistic models on the CCM model library matches these 'rare' sequences. Due to this, the fourth meta-feature, which is the percent of the models supporting a prediction, may ignore sequences that are related with only a few probabilistic models in the library. The last meta-feature is calculated as the percent of models within the same phylogenetic clade that supports the match. For this feature, taxonomic assignment must be performed, however, taxonomic assignment of millions of sequences even at the phylogenetic clade level is a time expensive task and therefore it can rise the computational time of the pipeline. In addition, taxonomic assignment can contain errors because some sequences are simply too short to be correctly classified or too similars even belonging to different clades. Missing annotations can result in a significant bias in the quantification of expressed genes. This is more evident in metatranscriptomics studies. This domain has highly significant hits and is found in many transcripts with a coverage exceeding 50%. However, its score of SVM is low because there are only a few models supporting this prediction.

Conclusions

We provided evidence that the learning step introduced by CLADE is not suitable for metagenomics/metatranscriptomics datasets. The fragmentation of data and the presence of sequences which may be identified only by a few very divergent models can result in low SVM scores. Moreover, in simulated data, we can observe an increase of spurious detections above the 1e-8 E-value threshold and therefore, a much more refined way to discriminate false positives coming from non-significant hits should be introduced. During the next chapter we

Calssification of environmental sequences

Ecosystem changes are often correlated with the presence of new communities entering them and disturbing their stability by importing new metabolic activities [START_REF] Jackson | Changes in community properties during microbial succession[END_REF][START_REF] Tyson | Community structure and metabolism through reconstruction of microbial genomes from the environment[END_REF][START_REF] Freilich | Competitive and cooperative metabolic interactions in bacterial communities[END_REF][START_REF] Johnson | Metabolic specialization and the assembly of microbial communities[END_REF]. Very often, such communities, their functional behavior and their mutual interactions are hard to identify and analyse [START_REF] Thompson | The geographic mosaic of coevolution[END_REF][START_REF] Whitham | A framework for community and ecosystem genetics: from genes to ecosystems[END_REF][START_REF] Chisholm | Theory predicts a rapid transition from niche-structured to neutral biodiversity patterns across a speciation-rate gradient[END_REF][START_REF] Dini-Andreote | Disentangling mechanisms that mediate the balance between stochastic and deterministic processes in microbial succession[END_REF][START_REF] Hand | Landscape community genomics: understanding eco-evolutionary processes in complex environments[END_REF]. Unraveling the structure of the community, what is functionally done by the community and its subcommunities is crucial for understanding their metabolic dynamics and activities. Any computational study contributing to improve the detection of the functional preferences of environmental communities and pin-point their interactions is important to our understating of the ecosystem changes [START_REF] Handelsman | Metagenomics: application of genomics to uncultured microorganisms[END_REF][START_REF] Allen | Community genomics in microbial ecology and evolution[END_REF][START_REF] Tyson | Cultivating the uncultivated: a community genomics perspective[END_REF][START_REF] Delong | Community genomics among stratified microbial assemblages in the ocean's interior[END_REF][START_REF] Eisen | Environmental shotgun sequencing: its potential and challenges for studying the hidden world of microbes[END_REF][START_REF] Van Straalen | An introduction to ecological genomics[END_REF].

Large scale environmental effects, induced by nutrients or temperature, guide the interest of developing tools to zoom in the metabolic activities of communities and to compare environments in detail. The aim is to relate genetic information with environmental factors and to understand how these factors affect the genetic material and the dynamics of the expression from one environment to another. Here, the word "environment" refers to individuals in a population (e.g. patients), to geographical sites (e.g. lakes), to environmental niches (e.g.

corals).

One of the most challenging aspects in the classification of environmental sequences is the loss of genetic signal due to fragmentation. This means that coding sequences (CDS) might be only partially sampled and genetic signals might become hard to detect. To overcome this difficulty, one possibility is to assemble reads into contigs/scaffolds. However, the effectiveness of the assembly relies on the community structure and its species diversity.

Communities with many closely related strains will lead to a poor assembly performance [START_REF] Charuvaka | Evaluation of short read metagenomic assembly[END_REF].

MetaCLADE approach for classifying environmental sequences

The approach we use is to annotate reads directly through their domain identification.

This allows us to deal with functional units much shorter than CDSs, and yet sufficiently precise to inform us on the potential functional activity of the communities. We should notice that, even though domain sizes vary from a few tens up to several hundreds of amino acids, 90% of the known domains are smaller than 200 residues and the mean size of a domain is 100 residues [START_REF] Richardson | The anatomy and taxonomy of protein structure[END_REF][START_REF] Janin | Structural domains in proteins and their role in the dynamics of protein function[END_REF][START_REF] Xu | Favorable domain size in proteins[END_REF]. Therefore, it is reasonable to look for functional signals in short reads of at least 150 nucleotides, with the hope of finding significant evidence for a domain.

Hence, the domains identified on reads can be quantified and classified into functional classes. Under the hypothesis that the most populated classes define the community preferences, we shall infer the community functional importance. Differences in domain counts for specific classes can be used to characterize and compare environments. Here, we wish to improve domain annotation, provide a better resolution of the functional activity of the community, and clarify, as much as possible, preferences and missing functional features of a community. We shall show several examples of comparative analysis of metagenomic data in different environments that illustrate the advantage of applying a precise domain annotation method to highlight specific characteristics of the environments.

We introduce MetaCLADE, the next generation metagenomics and metatranscriptomics annotation method based on the multi-source strategy. Multi-source annotation was introduced in [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF], where it was applied to fully sequenced genomes. Some seminal ideas on the usage of multiple models for domain representation were introduced in [START_REF] Rehmsmeier | Phylogenetic information improves homology detection[END_REF][START_REF] Gough | Assignment of homology to genome sequences using a library of hidden Markov models that represent all proteins of known structure[END_REF][START_REF] Callebaut | Prediction of the general transcription factors associated with RNA polymerase II in Plasmodium falciparum: conserved features and differences relative to other eukaryotes[END_REF][START_REF] Lees | Gene3D: a domain-based resource for comparative genomics, functional annotation and protein network analysis[END_REF] and are found in systems like SUPERFAMILY [START_REF] Gough | Assignment of homology to genome sequences using a library of hidden Markov models that represent all proteins of known structure[END_REF] and Gene3D [START_REF] Lees | Gene3D: a domain-based resource for comparative genomics, functional annotation and protein network analysis[END_REF][START_REF] Lees | Gene3D: Multi-domain annotations for protein sequence and comparative genome analysis[END_REF]. Multi-source annotation is based on the idea that protein evolution pathways are bound to be few due to the numerous structural and functional constraints that a protein undergoes. This means that the evolutionary constraints driving a protein evolution in a specific species and the corresponding (conserved, structural, physico-chemical) signals identifiable in a sequence, might be more easily detectable by looking closely at the evolutionary solutions found by some other species. In fact, we hypothesized that some species share their evolutionary solutions with the species to be annotated, even if they are phylogenetically distant from it. Hence, we used a few hundred clade-centered models (CCMs) associated to each single protein domain. These models span regions of the protein sequence space that are not well represented in a model based on sequence consensus (SCM). They might highlight motifs, structural characteristics or physico-chemical properties shared by a reference sequence and a pool of sequences similar to it. Hence, if the original set of reference sequences for a domain is made of divergent homologs, CCMs are expected to describe properties that could be missed by the SCM representing global consensus.

The great improvement in annotation obtained with the multi-source strategy compared to the mono-source strategy, employed by the two most commonly used annotation tools HMMer [START_REF] Eddy | Accelerated Profile HMM Searches[END_REF] and HHblits [START_REF] Soeding | Protein homology detection by HMM-HMM comparison[END_REF][START_REF] Remmert | HHblits: lightning-fast iterative protein sequence searching by HMM-HMM alignment[END_REF], was proven in [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF] for genomes. The usage of CCMs is likely to be even more important in metagenomics, where domain divergence might be very large and the number of identifiable evolutionary solutions might be greater. With MetaCLADE, we go beyond species identification constituting the communities, which provide an important landmark signature for different environments, and characterise communities with respect to their functional activities.

From CLADE to MetaCLADE

MetaCLADE has been designed with the purpose of annotating metagenomics and metatranscriptomics reads. It is based on CLADE [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF], a multi-source domain annotation tool adapted to full genomes. The characteristics of metagenomic and metatranscriptomic reads are their short lengths and the fact that they contain multiple sequencing errors compared to full length CDSs. They demand the design of a special computational protocol taking into account the particular nature of the data

The multi-source annotation strategy and the CLADE library

The main idea shared by MetaCLADE and CLADE is the multi-source annotation strategy: several probabilistic models, instead of a single model generated from the consensus of a set of homologous sequences, are used to represent a protein domain. The mono-source strategy typically performs well when the sequences are highly conserved. For this case, the consensus model captures the most conserved features in domain sequences and it can be successfully used to find new domains in databases of sequences, sharing the same features as the original sequence. However, when sequences have highly diverged, consensus signals become too weak to generate a useful probabilistic representation and models constructed by global consensus do not characterize domain features properly. To overcome this fundamental bottleneck, CLADE introduced clade-centered models (CCMs). For each Pfam domain D i , it considers the FULL set of homologous sequences S i in Pfam [START_REF] Finn | Pfam: the protein families database[END_REF] associated to D i , and for some representative sequences s j in S i (see below), it constructs a model by retrieving with PSI-Blast a set of sequences similar to s j from the NR database. The probabilistic model generated in this way displays features that are characteristic of the sequence s j , and that might be very different for other sequences s k in S i . The more divergent the homologous domain sequences s j and s k are, the more models constructed from these sequences are expected to display different features. It is therefore important for a domain D i to be represented by several models that can characterise its different pathways of evolution within different clades. The multi-source annotation strategy has proven more efficient than the mono-source annotation strategy when applied to full genomes [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF].

The CLADE pipeline

For each domain, the CLADE library includes the Pfam sequence consensus models (SCM) and at most 350 clade-centered models (CCM), with an average of 161 models per domain. The representative sequences associated to these models are selected in order to span most of the tree of life, the underlying idea being that evolutionary patterns can be found in species that are very far apart in the tree. This amounts to almost 2.5 million probabilistic models (constructed in about 4 months of computer time on a 250 nodes cluster). It should be noticed that the clade-centered domain library has been built on Pfam database (version 27), but it can be extended to other domain databases. The CLADE pipeline combines the output of its rich database of probabilistic models with a meta learning strategy in order to determine a set of best predictions for each domain sequence. Then DAMA is used to find the best domain architecture, by using information on domain co-occurrence and by exploiting multi-objective optimisation criteria [START_REF] Bernardes | A multi-objective optimization approach accurately resolves protein domain architectures[END_REF].

The MetaCLADE's pipeline

MetaCLADE's pipeline is illustrated in Figure 37. It is based on two main steps and a precomputed learning step.

MetaCLADE domain hits identification and the precomputed learning step

The first step takes as input a set of metaG/metaT reads where CDS/ORFs were previously identified, and the CLADE model library. For each sequence, the CDS region is scanned with the model library and all domain hits are identified. For each hit, two scores are considered, the bit-score for the entire hit and the mean-bit-score (i.e. the bit-score of the hit divided by its length).

The MetaCLADE precomputed learning step estimates domain specific gathering thresholds for combinations of bit-scores and mean-bit-scores, that best separate positive hits from negative hits. It is realized with a learning procedure that divides the two-dimensional sequence space (bit-score and mean-bit-score) associated to a domain into several regions, each with an associated probability. below). Hence, each CDS sequence analysed in MetaCLADE first step will be represented by a set of domain hits, where each domain hit is described by its bit-score, mean-bit-score, predicted state (positive or negative) and probability of the predicted state. The bit-score and the mean-bit-score are provided by the matching of a model in the CLADE library. The predicted state of the hit and its probability are provided by the learning step.

MetaCLADE domain selection

MetaCLADE second step filters the set of hits as follows:

1. All pairs of overlapping hits associated to the same domain with an overlap region that covers at least 85% of the length of each hit, are filtered from the list of CDSs.

For such pair of hits, we eliminate the hit that has the lower bit-score. The filtering is realized independently for CCMs and for SCMs. This filtering step eliminates multiple overlapping hits coming from the same model (either a SCM or a specific CCM).

Based on the parameter estimation obtained with the Naive Bayes classifier applied

to each Pfam domain (see "A Naive Bayes classifier sets two-dimensional thresholds for fragmented domains" below), MetaCLADE filter accepts only the domain hits whose probability p for being a positive hit is p > 0.9 and whose bitscore is greater than a domain sensitive lower bound identified by the classifier (i.e.the smallest bit-score of the negative sequences considered by the classifier).

3. Hits are sorted and filtered by a ranking score based on the bit-score and the percentage of identity against the model: for each hit a value between 0 and 1 representing the statistical significance of the bit-score (see Table 4) is multiplied by the percentage identity of the match (represented as a value between 0 and 1).

Namely, domain hits are ordered by higher ranking scores, then iteratively eliminated if they share at least 10 residues with some domain with higher score.

Table 4 Ranking scores for selecting hits on the third step of domain selection

The output of this filtering step is the CDS annotated with non-overlapping domain hits.

Due to the short length of the reads, one expects an annotation of one or two domains per read, possibly flanked by domain fragments on the right and/or the left. Consequently, and in contrast to CLADE, there is no reconstruction of the best architecture with DAMA. Due to the amount of data to be annotated, MetaCLADE provides exactly one annotation for each read. Also, note that the first filter is used to reduce the size of the set of domain hits to select, possibly huge at the beginning due to redundant predictions. The second filter is used to identify hits with a high probability to be positive hits. The third filter is used to identify the best solution. As a consequence of the construction of the probability space of sequences for a domain, note that the second filter asks for domain hits to have a bit-score greater than the smallest bit-score of the negative sequences in the space. This is because negative sequences considered by the classifier are a selected sampling of the space of negatives (see "Generation of negative sequences" below). Namely, among all negatives generated by the algorithmic procedure, we selected those that lie further away from the origin, and that, in consequence, have the highest statistical significance. These selected negative sequences tend to group together further from the origin of the space, and to lie at the borderline of regions characterized by positive sequences. Hence, one should properly evaluate the acceptance threshold against this specificity.

The precomputed step in MetaCLADE defines two-dimensional gathering thresholds

Metagenomic and metatranscriptomic samples demand to annotate domain fragments, possibly of small length. To explicitly distinguish small hits from long ones, MetaCLADE explicitly estimates the likelihood for a small hit to be a positive sequence by considering the bit-score of the hit and also its mean-bit-score. Namely, it defines a two dimensional gathering threshold (GA) for each domain by combining bit-score and mean-bit-score and by identifying multiple regions in the two-dimensional sequence space that, with a high probability, provide reliable annotations for short sequences. Probabilities are estimated with a Naive Bayes classifier and the statistical procedure is explained below. 3. by constructing a Markov model of order 3 for each domain and by using it to generate random sequences with positional probabilities.

Construction of positive and negative training sets

Note that in 3, the space of 160 000 (204) 4-tuples is evaluated by assigning a probability to appear in a domain sequence to each 4-tuple. This is done with a pseudo-count, by considering each 4-tuple to appear at least once, and by counting the number of occurrences n of the 4-tuple in the SEED sequences of the domain. The probability of a 4-tuple is set to [START_REF] Eddy | Profile hidden Markov models[END_REF] where n is the number of occurrences for this particular 4-tuple and N is the total number of 4-tuple occurrences in the SEED set associated to this domain. The Markov model of order 3 is defined on these probability estimations Globally, we ensure the full training set is comprised of roughly 50% of positive sequences and 50% of negative ones (Figure 39B). This proportion varies from domain to domain and depends on the difficulty to generate correctly annotated random sequences. In Figure 39, we report the proportions of negative sequences for CCMs and SCMs generated by the first two methods (Figure 39A) and compare them to the distributions of sequences generated all three methods (Figure 39B). Clearly, the third method contributes the largest number of negative sequences for each domain and establishes the expected numerical balance between the two training sets, all domains confounded. Definition of positive sequences. The training set of positive sequences was constructed as follows. For each domain D i and for each sequence in the Pfam SEED set of homologous sequences for D i , we created a set of prefixes and suffixes of the sequence to simulate small domain portions coming from the beginning or the end of the domain sequence that may be found in metagenomics reads. The maximum size of prefixes and suffixes was set to 30% of the entire domain sequence length and to a maximum of 100 aa. Fragments were determined by sliding a window of a given length with a fixed step size. For domains comprising between 15 and 75 aa the sliding window step was set at 5 aa, for sizes > 75 aa it was set to 10 aa, and for sizes < 15 aa it was set to 1 aa. For large domains, > 270 aa (this corresponds to one standard deviation away from the mean in the distribution of domain model sizes as reported in Appendix C), we expect that reads may lie somewhere in the middle of the domain and therefore we extracted random sequences from the original sequence that were not already covered by small fragmentations of the extremes. Fragment positions were set by randomly choosing their first position along the middle part of the sequence and fragment lengths were randomly picked from a normal distribution with mean 50 and standard deviation 25.

A Naive Bayes classifier sets two-dimensional thresholds for fragmented domains

Positive and negative sequences were put together and analyzed to obtain best separation parameters for CCMs and SCMs. Namely, we use a discrete version of the Naive Bayes classifier [START_REF] Hall | The WEKA Data Mining Software: An Update[END_REF] (downloadable from www.cs.waikato.ac.nz/ml/weka/citing.html) to construct learning models for each Pfam domain. The discrete version of the Naive Bayes classifier provides a finite partition of the sequence space and an estimation of the probability for a sequence to be a positive or a negative hit. Notice that we realize two different analyses, one on CCMs (generated by PSI-Blast) and the other on SCMs pHMMs generated by HMMer), because we cannot immediately compare their bit-scores. By so doing, we determine two distinct separation spaces and appropriated parameters for the two model predictions. In particular, only one probability space is estimated for all CCMs of a domain. All positive and negative sequences, generated for all CCMs, are considered in the same sequence space and the associated probability space is computed. Figure 40A illustrates an example of separation of the spaces of positive and negative hits for the CCMs and SCM of a Pfam domain, analyzed with the Naive Bayes classifier. Note that short fragments have small bit-scores with a possibly large mean-bit-score, and that negative sequences are characterized by small bit-scores and small mean-bit-scores. Also, the identification of fragmented coding regions (especially important for the annotation of real negative training sets"). In the CCM space, negative sequences are more clearly separated by both bit-scores and mean-bit-scores from positive ones than in the SCM space. In fact, since

CCMs are "closer" to sequences than SCMs, one expects their scores to be higher for positive sequences in CCMs than in SCMs. Also, the usage of a two dimensional sequence space, determined by bit-scores and mean-bit-scores, improves the separation of positive and negative sequences in MetaCLADE compared to HMMer (hmmscan). In the plot describing the SCM space (Figure 40A, right), the GA threshold excludes from hmmscan predictions a large majority of Bacteriorhodopsin-like protein sequences detected by CCMs. More generally, for all Pfam domains, we computed the difference between the GA threshold associated to the SCM and the mean of the bit-scores for the 5 negative sequences with higher bit-score identified by the SCM. The distribution of the differences, displayed in Appendix K, shows a small standard deviation suggesting that MetaCLADE and Pfam/HMMer estimation of the (one dimensional) cut-off is similar. This control shows that Naive Bayes classification produces reasonable thresholds when projected in one dimension. We generated a set of fragmented sequences from a set of fully sequenced genomes made of 11 archaea and 44 bacteria. The list of genomes is reported in Appendix A. The genomes have been fragmented with MetaSim [START_REF] Richter | MetaSim: a sequencing simulator for genomics and metagenomics[END_REF] and the outcoming clones have been parsed

with FlowSim [START_REF] Balzer | Characteristics of 454 pyrosequencing data-enabling realistic simulation with flowsim[END_REF], to simulate realistic insertion and deletion errors expected during DNA sequencing. This resulted in about 500 000 reads for which FragGeneScan [START_REF] Rho | FragGeneScan: predicting genes in short and error-prone reads[END_REF] predicted 602 968 coding sequences (pCDS) that were given as input to MetaCLADE.

MetaSim was run to generate clones, based on a normal distribution with a clone mean size of 800nt and a standard deviation of 100nt. We assumed all species be equally abundant.

Simulated clones were "sequenced" with FlowSim by using the 454 Titanium generation model. FragGeneScan was set to use the 454 Titanium model for indel error correction. For this simulated metagenomic dataset, we selected positive domain hits with a probability threshold of 0.85 estimated by the Naive Bayes classifier.

The use of clans and Interpro families

To evaluate MetaCLADE on the dataset of simulated sequences, we use Pfam domain clans [START_REF] Finn | Pfam: clans, web tools and services[END_REF] and Interpro families [START_REF] Mitchell | The InterPro protein families database: the classification resource after 15 years[END_REF]. Clans are groups of proteins for which common ancestry can be inferred by similarity of sequence, structure or profile-HMM [START_REF] Finn | Pfam: clans, web tools and services[END_REF]. The list of Pfam clans was retrieved at ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/. Interpro families represent groups of evolutionarily related proteins that share common functions. Such entries tend to be near full length and typically do not undergo recombination, in contrast to domains [START_REF] Mitchell | The InterPro protein families database: the classification resource after 15 years[END_REF]. The list of all Interpro families and its relationships was retrieved from the ebi portal at https://www.ebi.ac.uk/interpro/download.html

Real metagenomic and metatranscriptomic datasets

To validate MetaCLADE on real data, we analyzed eleven metagenomics and metatranscriptomic samples. The characteristics of these eleven datasets, such as number of reads, average read size, sequencing technique used to generate the dataset, whether it is a metagenomic (metaG) or a metatranscriptomic (metaT) dataset, are provided in Appendix D.

Available websites for download are given in Appendix E. Five metatranscriptomic samples come from different geographic locations in the oceans, Antarctic (ANT), North Pacific (NPAC), Equatorial Pacific (EPAC), Arctic (ARC) and North Atlantic (NATL) [START_REF] Toseland | The impact of temperature on marine phytoplankton resource allocation and metabolism[END_REF]. We have identified domains on the reads by analyzing the translation of the six reading frames and annotated them with HMMer and MetaCLADE. Four published metagenomics datasets come from very different environments: soil, ocean, ancient bones and guts. For the gut environment, we also considered a metatranscriptomic sample. These five sets of reads were previously analyzed for CDS identification by EBI with FragGeneScan [START_REF] Rho | FragGeneScan: predicting genes in short and error-prone reads[END_REF]. CDS sequences have been annotated by EBI based on five different domain databases found in InterPro [START_REF] Mitchell | The InterPro protein families database: the classification resource after 15 years[END_REF]: Pfam [START_REF] Finn | Pfam: the protein families database[END_REF],

TIGRFAM [START_REF] Haft | TIGRFAMs and genome properties in 2013[END_REF], Gene3D [START_REF] Lees | Gene3D: Multi-domain annotations for protein sequence and comparative genome analysis[END_REF], PRINTS & ProSite [START_REF] Attwood | The PRINTS database: a fine-grained protein sequence annotation and analysis resource -its status in 2012[END_REF][START_REF] Sigrist | New and continuing developments at PROSITE[END_REF]. The search was realised with InterproScan [START_REF] Mitchell | EBI metagenomics in 2016 -an expanding and evolving resource for the analysis and archiving of metagenomic data[END_REF] as the final step of the EBI Metagenomics pipeline (www.ebi.ac.uk/metagenomics/about). These systems are accessible from www.ebi.ac.uk/services/proteins. The O'Connor lake metagenomic dataset was downloaded from the EBI Metagenomics portal (www.ebi.ac.uk/metagenomics/projects/ERP009498). It was realised with Illumina HiSeq 2000 technology and contains 1 315 435 very short reads, with an average length of 123nt. CDSs were identified by EBI. On all these datasets, MetaCLADE selected positive domain hits with a probability threshold of 0.9

Domain abundance

The functional analysis of a metagenomic/metatranscriptomic sample is realized by characterizing domain abundance within a functional class with a normalized value between 0 and 1. This normalization is done by dividing the number of domains detected in a functional class by the total number of domains belonging to the most represented class in the environmental sample. We speak about "normalised abundance". A second kind of normalisation is realised with respect to multiple environments and it is used for comparing domain abundance within the same functional class across these environments. A normalised domain abundance , where S is the sample and I is the domain, is computed as the product of the actual domain abundance per megabase by the average size of all samples. By multiplying by the average size of all samples, we provide an indication of the expected number of domains if all environments had the same size and can compare environments with respect to such estimations.

Functional analysis of annotated real datasets

To validate MetaCLADE on all real metatranscriptomic and metagenomics datasets, we associated a function based on GO classification to both the domains identified with MetaCLADE and the domains identified with HMMer (hmmscan). We used Pfam2GO [START_REF] Camon | An evaluation of GO annotation retrieval for BioCreAtIvE and GOA[END_REF] and annotated biological process terms of GO-slim [START_REF] Mitchell | The InterPro protein families database: the classification resource after 15 years[END_REF][START_REF] Hunter | EBI metagenomics -a new resource for the analysis and archiving of metagenomic data[END_REF]. Pfam2GO was retrieved from geneontology.org/external2go/pfam2go and GO-slim classification for metagenomics was retrieved from geneontology.org/page/download-ontology. To highlight the differences between MetaCLADE and HMMer, we compared domain abundance in all GO-term classes.

For this, we normalized domain count in each metatranscriptomic dataset with respect to the size of the sample as described above.

Software used to analyse the data and to compare performances

To compare MetaCLADE, we used several systems and domain model libraries:

Pfam/HMMer [START_REF] Finn | Pfam: the protein families database[END_REF], TIGRFAM [START_REF] Haft | TIGRFAMs and genome properties in 2013[END_REF], Gene3D [START_REF] Lees | Gene3D: Multi-domain annotations for protein sequence and comparative genome analysis[END_REF], PRINTS & ProSite [START_REF] Attwood | The PRINTS database: a fine-grained protein sequence annotation and analysis resource -its status in 2012[END_REF][START_REF] Sigrist | New and continuing developments at PROSITE[END_REF].

Annotations produced by these systems were downloaded from available annotation files provided by the EBI metagenomics pipeline https://www.ebi.ac.uk/metagenomics/pipelines/ (versions 1.0 or 2.0 depending on the dataset).

The EBI metagenomics pipeline uses InterProScan (v5.0 for pipeline 1.0 and 5.9-50 for pipeline 2.0) as tool to annotate predicted CDS (using FragGeneScan 1.15) using the domain libraries above. HMMer v3.0 was downloaded from hmmer.org/ and HMMScan was run with default parameters and curated inclusion thresholds. The option --cut ga, for model-specific thresholding (using profile's GA gathering cutoffs to set all thresholding), was used.

MetaCLADE software

The pipeline is implemented in Python 2.7 and is available at http://www.lcqb.upmc.fr/metaclade. This includes the annotation system (MetaCLADE two main steps in Figure 37) and the program pre-computing domain specific gathering thresholds (MetaCLADE precomputed step in Figure 37). The CLADE model library used in MetaCLADE was constructed based on Pfam database v27 and was released with CLADE [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF].

It can be found at http://www.lcqb.upmc.fr/CLADE.

Results

Metagenomics and metatranscriptomics datasets can be explored to learn about the functional activity of the community. We show that it is possible to reach a very fine degree of accuracy in annotation though a state-of-the-art multi-source domain annotation method.

MetaCLADE is especially designed to handle metagenomic and metatranscriptomic data. The multi-source annotation strategy encoded in MetaCLADE and the predictive power of its multiple models have been already proven for annotation of full genomes [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF]. Here, we demonstrate that the multi-source annotation strategy and a careful handling of short reads can provide a very high predictive power in metagenomic datasets. On a simulated dataset, we test the performance of the learning procedure in MetaCLADE and, on eleven environmental samples satisfying very different characteristics, we demonstrate that MetaCLADE annotation highly improves current methods

A simulated metagenomic dataset: reads annotation based on clans

We defined a set of reads starting from 56 fully sequenced genomes from archaea and bacteria (Appendix A). We applied MetaSim to generate 500 000 clones and, with FlowSim, we broke up the clones into reads by adding suitable errors to the sequences simulating sequencing behavior induced by 454 Titanium last generation sequencing. Then, with

FragGeneScan, we predicted CDS in our simulated reads and annotated them with MetaCLADE.

In order to evaluate MetaCLADE performance, we identify domains at the clan or Interpro family level. This is done because sequence similarity within domains in the same Pfam clan is usually high and genome annotation is often misleading by domains belonging to the same clan. This is even more true in metagenomics datasets, where one often needs to annotate fragments of a domain displaying a weaker signal due to the reduced length.

We compare the annotation obtained by MetaCLADE against the annotation given by MetaCLADE shows that the larger amount of domains it detects falls coherently in functional classes of interest for specific environments. A functional annotation of domains in the five oceans metatranscriptomic data reported in [START_REF] Toseland | The impact of temperature on marine phytoplankton resource allocation and metabolism[END_REF], sharply shows this point. The mapping of domains identified with MetaCLADE and HMMer (hmmscan) allowed us to show that MetaCLADE reaches a much better resolution of significant terms among all GO-Slim functional classes. In Figure 42A, we report the relative abundance of reads found in the 5 metatranscriptomic samples. Certain functional classes, such as "translation", are over represented, as expected, for both MetaCLADE and HMMer. Others are characteristic of certain environmental conditions, and they are only detected by MetaCLADE. A striking example is the "ion transport" GO-slim functional class for the EPAC and ANT samples. In this specific case, HMMer annotation completely missed the large presence of bacteriorhodopsin-like domains in EPAC, as illustrated in Figure 42B. In other environments, such as the ANT sample, there is a much weaker bias towards these specific domains but their presence is nevertheless captured: the bacteriorhodopsin-like domains are classified by GO as "ion transport" and they are represented by a red color of the associated box in Figure 42D compared to the yellow color in Figure 42C. Hence, even though there is no particular bias towards a specific domain, MetaCLADE annotation is much finer than HMMer annotation. In general, this is illustrated by the hierarchical tree-graphs of the GO-slim functional class "ion transport" reported in Figure 42C for HMMer and in Figure 42D for MetaCLADE. The MetaCLADE tree-graph is much more detailed and precise in the annotation of domains. It contains seven nodes more than the corresponding HMMer tree-graph. Note that "metal ion transport" for instance, represented by just one node of 44 identified domains in the graph for the HMMer analysis (Figure 42C), is detailed by a more complex MetaCLADE sub-tree-graph of 165 identified domains, associated to iron ion, nickel cation, cobalt ion and ferrous iron transport (Figure 42D). This association to specific functional roles of the identified domains can help the biologist to better characterize the metabolic regimes of the sample. Overall, MetaCLADE uniformly annotates more domains and with a more specific functional association than HMMer.

In Figure 42B, some functional classes appear as the most represented in exactly one environmental sample. This is the case for the pyrophosphates in NATL, the transferrin and the ammonium transporter in NPAC. Others, are shared by several samples. They might be present in the remaining samples as well, but relatively less represented. This comparative information is crucial for zooming in the functional activity of an environment.

Finally, one should notice the distribution of species used by MetaCLADE to annotate domains in the five oceanic samples (Figure 42E). These eukaryotic read sequences were mostly annotated with CCMs generated by Metazoan and Alveolata domain sequences. A large contribution from other organisms, as Bacteria, is also present as expected.

MetaCLADE identification of divergent domains by conserved small motifs

MetaCLADE multi-source annotation strategy is used with the purpose to identify very divergent domain sequences lying in reads. CCMs are probabilistic models that describe closely specific sequences and they can capture conserved patterns that are specific of homologs niches and that are missed by SCMs. As a consequence, several CCMs for a domain have the possibility to describe domain sequences in greater detail and span a greater space of homologous sequences, including a larger number of very divergent ones. For instance, in Figure 43A, we consider the conservation profile of the sequence alignment associated to a CCM used in the annotation of the Rhodopsin-like domain in metagenomic fragments, missed by HMMer as discussed above. With this and other CCM models, MetaCLADE could annotate 360 sequences in EPAC, that could not be detected by HMMer in [START_REF] Toseland | The impact of temperature on marine phytoplankton resource allocation and metabolism[END_REF] due to the strong sequence divergence (Figure 43D). The conservation profile of the alignment of the 360 environmental sequences is reported in Figure 43B. It is very conserved and corresponds to a portion of the rhodopsin-like domain, whose expression is expected in the Equatorial Pacific.

This conserved pattern makes the third of the length of the entire domain. The rest of the sequence is divergent and remains with no annotation. One can visually appreciate the stronger similarity of the CCM profile (Figure 43A) to the metagenomic sequences (Figure 43B) compared to the Pfam SCM profile (Figure 43D) of the bacterial-like rhodopsin. Note that the motif identified by MetaCLADE in the eukaryotic metagenomic sample was recently identified in the dinoflagellate Prorocentrum donghaiense [START_REF] Sanger | DNA sequencing with chain-terminating inhibitors[END_REF] (Figure 43C) with an alignment comprised by homologs from Oxyrrhis marina and bacteria. The conserved positions characteristic of the dinoflagellate sequence [START_REF] Shi | Rhodopsin gene expression regulated by the light dark cycle, light spectrum and light intensity in the dinoflagellate Prorocentrum[END_REF] are recovered in the alignment of our metagenomic sequences, confirming MetaCLADE functional annotation.

Time complexity

The most expensive phase is the domain hit identification, where all probabilistic models (CCMs and SCMs) are used as query to search against a database of protein sequences.

This step can be run in several cores resulting in a total time of a couple of days for a dataset of ~200 MB in a cluster with 64 cores if a gene detection tool such as FragGeneScan is used to reduced the search space of coding sequences. As the cost of high performance computing is decreasing, we shall think that our pipeline is still suitable for a typical metagenomics study.

The first step of filtering, which needs to sort identified hits in order to remove redundancy, is also time consuming, but approximatively 10x less than the domain hit identification step. The pre-computed step allows MetaCLADE to determine the Naïve Bayes probability associated to a prediction (second filter) in log2(n) where n is the number of probability regions associated to a domain (training set), typically ~20 regions per training set. The third step of the filtering runs in a single core but it only takes a few minutes. Table 8 shows the total time per core to compute each of the steps of the pipeline for several datasets. investigation of the geographical, genetic, temporal and physiological characteristics of an environment. To capture common and rare entities in a given environment, functional annotation methods need to be as precise as possible in identifying remote homology.

We have shown that MetaCLADE outperforms existing domain annotation tools. In particular, we report that MetaCLADE identifies key domains, such as the rhodopsin-like domains for the Equatorial Pacific Ocean, that are missed by most used annotation tools.

Rhodopsin-like domains are not particularly highly expressed but they are present in a specific environment as among the most represented domains compared to others. This differential information turns out to be crucial for the understanding of environmental features and differences in metabolic activities of different environments.

A meta learning step and a domain selection designed for reads. MetaCLADE has been especially designed to annotate metagenomic and metatranscriptomic data. In this respect, it presents several differences compared to CLADE. First, the training step defining a twodimensional probability space for each domain was designed in such way that the tool is able to determine the significance of a match taking into account its size. This is very important because domain fragmentation is always present in metagenomics data. Even in sequencing technologies producing long reads, the high rate of insertion/deletion errors causes frame shifts resulting in non-continuous signals.

Based on the probability space constructed for each domain, we modified the way to combine together hits obtained from different CLADE models (CCMs and SCM). Namely, MetaCLADE filters domain predictions first by selecting overlapping hits of the same domain with best bit-score, then by selecting all domains that are positive sequences with a high probability, and finally, by filtering overlapping domain hits keeping those with best bit-score and percentage of identity. Note that, at times, domains may reach high bit-scores even by correctly matching parts of the models where conservation is weak, therefore favouring conserved parts is a crucial step to give an accurate final prediction. Note that CLADE considers other features for hit selection, like the overlapping of multiple hits coming from different models of the same domain, the E-value of domain hits, the phylogenetic origin of the models identifying hits. Domain coexistence also helps to screen predictions and filter out false positives. Such features can be used for genome annotation but are hardly fitting reads annotation.

Some perspective to improve MetaCLADE. MetaCLADE was especially designed to consider the partial information contained in domain fragments, localised in reads. For this, we defined a powerful two dimensional domain dependent gathering threshold and we use multiple models to represent each domain, possibly characterising small conserved motifs for the domain. In future development we foresee to improve the system in several ways (see also [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF]):

More domains and new models for an improved MetaCLADE annotation. An obvious improvement in MetaCLADE will be the extension of the reference domain library, based on Pfam domains, with the set of domains included in Gene3D and TIGRFAM. The motifs represented in PRINTS and ProSite could be also considered and the associated profiles handled in MetaCLADE. Constructing a library of conserved small motifs. At the moment, MetaCLADE does not have any specific handling of motifs but this extension will be considered in the future. Search for sequence motifs in an environmental sample might be realised with a computationally costly all against all read comparison. Alternatively, starting from the most conserved patterns comprised in CCMs, we shall generate a repertoire of significant motifs specific of each domain in order to improve hit selection criteria. A systematic classification of these motifs might lead to environmental datasets of motifs that could be used as environmental signatures of metabolic activities.

A suitable encoding of these "environmental patterns" can be used to find new domains in environmental samples with MetaCLADE. The advantage in this search, compared to a "all against all" strategy, beside the computational time, is the fact that patterns were constructed starting from an original domain search, possibly functionally annotated, and that this annotation could be used to associate a potential functional role to new domains discovered through the pattern.

Annotation of longer sequences. Another important extension of MetaCLADE concerns the possibility to work on long reads, assembled contigs. In this case one would expect to annotate large stretches of genes and possibly all coding sequences. As for CLADE, one could run a new step for the construction of the best architecture, by using DAMA [START_REF] Bernardes | A multi-objective optimization approach accurately resolves protein domain architectures[END_REF]. DAMA could be added as the last step of the pipeline.

Reduction of domain number in MetaCLADE. Some of the probabilistic models in MetaCLADE library are expected to be redundant and a suitable handling of these models, after their clustering, should help to increase the speed of the method and to preserve the same predictive power. Future development of MetaCLADE will ask for a reduction of the domain number.

New criteria to filter overlapping hits in MetaCLADE. Different domain hits could be selected by exploiting further the characteristics of the two dimensional space of sequences precomputed for the domains. For instance, one could privilege the domain hits with larger bitscore/mean-bit-score distance from the closest negative in the space. Filtering conditions of this kind could improve the annotation and need to be tested at large scale.

Extension of CCMs spreading through a larger number of species. New CCMs could be added to the library with the hope to reach novel and unrepresented evolutionary solutions for a domain. Note that, MetaCLADE provides the program to pre-compute gathering thresholds for all domain models. This allows the user to run the same thresholds for the domain on the new CCMs.

Chapter 4

MetaCLADE applications in environmental studies Our aim in this chapter is to describe how MetaCLADE may help scientists to analyze the huge amount of functional data that is produced by the analysis of metagenomics and metatranscriptomics data sets. This part of MetaCLADE is completely integrated with the rest of the pipeline and can offer a zooming on the functional landscape of microbial communities.

By means of several reports, MetaCLADE provides a simple but proper visualization of the functional classes associated to protein domains, allowing users to compare the functional preferences of different environments.

As mentioned in the introduction, MetaCLADE was used in a study of differential gene expression in sea diatoms under micro scale turbulences. As readers may want to retrieve the whole article from the peer-review publication, here we will briefly resume the main points of the article and we will explain the contribution of MetaCLADE to the study.

MetaCLADE functional annotation

As said before, the scope of MetaCLADE is to provide a characterization of the functional preferences of micro-organism communities through an enhanced domain annotation. The Gene Ontology (GO) database (see Chapter 1) provide an extended classification of well referenced gene products and locations. This database, being a public resource, has been used for scientists to develop its functional analysis. In metagenomics, the use of GO term enrichment analysis may highlight the functional differences among different environments. Moreover, as seen in the previous chapter the InterPro consortium maintains an InterPro2GO mapping between integrated domain signatures from databases such as Pfam and Gene Ontology terms. This allows us to directly map identified domains to GO terms and to construct functional profiles.

The directed acyclic graph (DAG) generated by the 3 different ontologies: Biological Processes, Cellular Component and Molecular Function provide more than 40 000 terms [START_REF]Gene Ontology Consortium: going forward[END_REF].

This DAG which can also be seen as a tree like structure, contains a hierarchical structure in which terms (nodes of the DAG) are much more specific as the distance increase from the root node. Since the amount of terms involved in a study can be overwhelming it is hence crucial to regroup terms at a certain level of the DAG in order to provide a much broad view of the content, this in order to summarize the most important functional classes without enter in the full detail of much more specific terms. For this reason, several groups have defined high-level subsets of gene ontology terms called the GO-Slim subsets. These subsets allow users to describe in a broad way the results of GO annotation of a genome, microarray, or cDNA collection. The metagenomics GO-Slim [START_REF] Hunter | EBI metagenomics -a new resource for the analysis and archiving of metagenomic data[END_REF] set, which is maintained by the EBI, was generated taking into account more than 22 billion of GO terms annotations associated to publicly available datasets submitted to the ebi Metagenomics portal. This slim set provide a high-level description of the main protein functions that may be present in a sample. Terms on it, were selected in order to cover most of the gene-ontology annotations that may be present in a sample, using the data from the 22 billion of GO terms they have accomplished to define a sub-set which covers more than 90% of these annotations. This set of ontologies can be seen as the most important functional classes that may be present in environmental data sets and they provide a high-level functional profile of microbial communities. MetaCLADE functional annotation has been therefore accomplish using information from both, the high level metagenomic GO Slim terms, and all its much more specific terms.

MetaCLADE not only provides the list of go terms associated to identified domains in a sample, but also a visual representation of data that can be useful to better understand functional activity in microbial communities. It generates different types of visualizations, depending at which level we want to analyze the data.

Marine diatoms sense and respond to microscale turbulence in nonlimiting nutrient conditions

MetaCLADE was used to refine domain annotation in an experiment of the morphological and gene expression responses of the diatom Chaetoceros Decipiens to microscale turbulences. In this part of the manuscript we will describe the main points of the publication and our contribution using MetaCLADE to refine domain annotation in metatranscriptomic sequences.

Diatoms are a fundamental microalgal phylum that thrives in turbulent environments. Despite several experimental and numerical studies, if and how diatoms may profit from turbulence is still an open question. These turbulences and in general the movements of water affect key factors for the survival of phytoplankton, such that the available light, nutrient concentrations, the abundance of prey and the temperature of the water [START_REF] Mckiver | Resonant plankton patchiness induced by large-scale turbulent flow[END_REF]. The movement of water introduces kinetic energy in the system and turbulence is the way by which kinetic energy is transmitted, through a cascade of dissipation, over several eddy-like structures down to the smallest scale. Bellow this scale, the energy is dissipated as heat via the friction of the viscosity, water motion cannot prevail over molecular diffusion but can control it by changing local gradients [START_REF] Stocker | Marine microbes see a sea of gradients[END_REF][START_REF] Taylor | Trade-offs of chemotactic foraging in turbulent water[END_REF]. This is particularly important for unicellular phytoplankton which are surrounded by a fluid boundary layer where molecular diffusion is the dominant process and only solute (nutrient) chemical gradients assure cell provisioning. A distortion of the boundary layer would change these gradients, hence nutrients would diffuse more rapidly, enhancing the uptake rate [START_REF] Pasciak | Transport limitation of nutrient uptake in phytoplankton[END_REF]. For non-motile cells, like diatoms, the distortion of the boundary layer can be produced only by sinking or by the shear generated by the decay of turbulent kinetic energy. Microturbulence may also favour unicellular autotrophs, bringing them into the upper illuminated layer of the ocean, the euphotic zone [START_REF] Margalef | Life-forms of phytoplankton as survival alternatives in an unstable environment[END_REF]. The arguments above, together with cell size and the ability to produce chains [START_REF] Arin | Phytoplankton size distribution and growth rates in the Alboran Sea (SW Mediterranean): short term variability related to mesoscale hydrodynamics[END_REF][START_REF] Machado | Influence of microscale turbulence on the phytoplankton of a temperate coastal embayment, Western Australia[END_REF][START_REF] Peters | Effects of small-scale turbulence on the growth of two diatoms of different size in a phosphorus-limited medium[END_REF][START_REF] Romero | Dynamic forcing of coastal plankton by nutrient imbalances and match-mismatch between nutrients and turbulence[END_REF][START_REF] Pahlow | Impact of cell shape and chain formation on nutrient acquisition by marine diatoms[END_REF], are frequently invoked to explain why diatoms should be favoured in turbulent environments [START_REF] Margalef | Turbulence and marine life[END_REF][START_REF] Wyatt | Margalef's mandala and phytoplankton bloom strategies[END_REF].

Mechanistic studies [START_REF] Pasciak | Transport limitation of nutrient uptake in phytoplankton[END_REF][START_REF] Karp-Boss | Motion of diatom chains in steady shear flow[END_REF] predict that diatoms could profit from turbulent pulses, even without any physiological adjustment. We note that following these studies some of the necessary conditions for an impact of microturbulence on diatoms, i.e., intense turbulence, nutrient depletion and grazing pressure, are infrequently met in the oceans [START_REF] Doney | From genes to ecosystems: the ocean's new frontier[END_REF][START_REF] Giraud | Importance of coastal nutrient supply for global ocean biogeochemistry[END_REF]. Thus, on the basis of current theories, diatoms would not specifically be adapted to microscale turbulence [START_REF] Barton | The impact of fine-scale turbulence on phytoplankton community structure[END_REF]. Nonetheless, that diatoms can sense mechanical stimuli was demonstrated by shaking a suspension of a Phaeodactylum tricornutum aequorin transformant with a needle and observing cytosolic calcium increases after 1-2 seconds after application of the stimulus and declining soon after [START_REF] Falciatore | Perception of environmental signals by a marine diatom[END_REF]. Cytosolic calcium waves trigger the activation of signal transduction [START_REF] Tuteja | Calcium signaling network in plants: an overview[END_REF] i.e.

a response to a perceived stimulus. This raises the question if turbulence produces changes on some proximal environmental variables which would initiate diatom response similarly to what happens for light [START_REF] Bailleul | An atypical member of the light-harvesting complex stress-related protein family modulates diatom responses to light[END_REF][START_REF] Coesel | Diatom PtCPF1 is a new cryptochrome/photolyase family member with DNA repair and transcription regulation activity[END_REF][START_REF] Depauw | Exploring the molecular basis of response to light in marine diatoms[END_REF] or if turbulence acts as a signal, carrying information about the environmental context that diatoms exploit to rearrange their physiology.

To address this question, controlled levels of turbulence have been applied on the diatom Chaetoceros decipiens using a generator of turbulence (TURBOGEN [START_REF] Amato | TURBOGEN: Computer-controlled vertically oscillating grid system for small-scale turbulence studies on plankton[END_REF]). The conditions are therefore similar to those we could find in the ocean. A Microscopic inspection and a transcriptome analysis have been carried out to determine if there is a response, morphological or in gene expression, to micro-scale turbulences in diatoms.

Results

RNA-Seq data were produced with the technology Illumina and used to generate a reference transcriptome for Chaetoceros decipiens, composed by 28 000 contigs and 21 000 unique transcripts, a summary of the sequences and annotations is shown in the Table 9. In order to assess whether C. decipiens cells exposed to turbulence experienced changes in their gene expression profile, a low density culture was aliquoted in the six TURBOGEN cylinders, turbulence was applied to three cylinders and samples were collected from still and turbulent cylinders at 48 and 72 hours from the beginning of the experiment (time points T2 and T3, respectively). About 12% (T2) and ∼2% (T3) of the total C. decipiens transcripts were differentially expressed, of which the vast majority were unique to T2 (Figure 51a). RNA-seq results were validated by qPCR using independent RNA samples. In KO [START_REF] Kanehisa | KEGG for linking genomes to life and the environment[END_REF] (KEGG ORTHOLOGY) pathway (level 1, Figure 51b) and GO-term (Figure 51c) enrichment analyses the frequency of a given term in the entire transcriptome was compared with frequency of the same term in the differential expressed (DE) transcript dataset, showing that a number of functions were significantly enriched in the latter. This is an indication that these functions might be involved in turbulence perception and/or response.

In addition to the DE analysis, cultures have been collected for a microscopic inspection and for a count of the cells every 24 hours from the beginning of the experience. TURBOGEN has also been used to measure the physiological response to micro turbulences in other two (ochre); biological processes (blue); molecular function (pink). On the x-axis, the percentage of transcripts associated to a given KO (b) or to a specific GO term (c) calculated over the total number of transcripts in the given dataset is reported.

# reads 2,3 x 10 Bioinformatics analyses made by the main authors. The fastq files from RNAseq experiments were inspected using FastQC [START_REF] Andrews | FastQC: a quality control tool for high throughput sequence data[END_REF] tool and further cleaned and trimmed using Trimmomatic [START_REF] Bolger | Trimmomatic: a flexible trimmer for Illumina sequence data[END_REF]. Trinity software [START_REF] Haas | De novo transcript sequence reconstruction from RNA-Seq: reference generation and analysis with Trinity[END_REF] (ver. trinity_201407) was used for the assembling of the reads. Quantification of transcript expression levels was done by mapping reads against the assembled transcriptome using Bowtie [START_REF] Langmead | Aligning short sequencing reads with Bowtie[END_REF] (ver. 1.1). To count the reads mapped Samtools [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF] (ver. 0.1.19-44428cd) was used. The resulting reference transcriptome was annotated using Annocript software [START_REF] Musacchia | Annocript: a flexible pipeline for the annotation of transcriptomes able to identify putative long noncoding RNAs[END_REF] that aligns transcripts against known proteins, domains and non-coding RNAs. Differential expression analysis. EdgeR [START_REF] Robinson | edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[END_REF] was used to select transcripts differentially (DE) between still and turbulent conditions. Transcripts were considered as DE if the false discovery rate (FDR) was smaller or equal to 0.05 and the fold change greater than 2. Enrichment analyses for GO and Pathways terms were performed exploiting the Fisher exact test [START_REF] Fisher | On the interpretation of χ2 from contingency tables, and the calculation of P[END_REF] and the Benjamini and Hochberg correction of the pvalues [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. GO terms and Pathways were considered enriched when associated to at least 10 DE transcripts with an adjusted p-value smaller than 0.1.

Our contribution. In order to refine domain annotation, we used META-CLADE, a software designed on the basis of CLADE with the purpose to annotate metagenomics and metatranscriptomics reads. We mapped Pfam domains into contig sequences, annotated the contigs and made the analysis (counting) on each sample. The importance of a functional class is highlighted by the abundance of the domains within the class. To compare the estimations obtained on each sample, we normalised the abundance with respect to the size of each sample (annotation per megabase). META-CLADE has provided annotations to 892 unannotated (by the other tools) contigs that are differential expressed under turbulence conditions.

Discussion

Diatoms are organisms that thrive in turbulent waters, incapable of active movements, they are universally believed to do better in turbulent waters only because turbulence mobilises nutrients and distorts the boundary layer around the cells, enhancing nutrient gradients towards the cell wall and therefore uptake. This implies that in nutrient repletion (i.e. when nutrients are never limiting for growth) diatoms should not be affected by turbulence. However, the results of this study suggest that there is a response both physiological and the expression of genes that was not expected. The results show evidence of the existence of mechanisms which allows diatoms to respond to mechanical stimuli.

Conclusion

Microbial communities are a vast source of genomic material. Next generation sequencing technologies provide an unmatched opportunity to explore the genomic potential of the 'unseen majority'. However, a paradigm shift in the analysis of environmental communities is mandatory as we are getting into the examination of potentially completely unknown mechanisms. In this context, it is suitable to take advantage of all available information rather than trying to reach an understanding of communities with a general abstraction of genomic information. For instance, when sequence consensus models attempt to characterize domain families that are very divergent, signals become too weak as they are trying to abstract common information that is not available. This does not imply that computational approaches providing abstraction/reduction of information are not suitable for the exploration of environmental communities, but rather that the information should be divided/regrouped and then extracted to provide a more meaningful characterization of the sequences. One should not underestimate the diversity that may be present in a family of closely related sequences. Until now, sequences have been classified in a 'broad manner', providing a certain number of known families and scientist trying to fit new sequences on those classifications. With the increasing number of 'new' sequences, we are understanding that this classification is not enough as families are being divided into subfamilies with slightly but important functional differences. The more sequences we explore the more this classification become much more specific.

We shall also, analyze information with the idea that the great diversity which is present in microbial communities has arisen not only from the 'individual' adaptation of the member organisms but also from the 'combined' adaptation of the communities as a whole. Molecular mechanisms may evolve in an orchestrated manner in which interactions between member organisms play a crucial role.

From a computational view, scientists should also change their view with respect to the resources that are needed to analyze environmental data. One cannot expect to try to unveil all the genomic diversity without the assistance of high performance computing (HPC). In parallel to the reduction of the cost per nucleotide of next generation sequencing technologies, HPC has also decreased its cost and many solutions are available at a reasonable cost. This does not replace the good practices in software development and the eternal search for algorithms or heuristics with low complexity however, one cannot dismiss approaches such as the multisource strategy, which increases the complexity but can provide a significant improvement in the quality of annotations.

During this work, we have been able to see the importance of the annotation of protein domains in the characterization of microbial communities. Some metabolic functions are coded in sequences that cannot be detected by the probabilistic models constructed from a set of very divergent sequences. Moreover, when predictions of models are combined with a trusted threshold taking into account the fragmentation of domains, we are able to discriminate and recover true signals coming from closely related sequences for which abundance may for instance, in metatranscriptomics studies, unveil an important preference in a community.

We have provided evidence that the approach multi-source introduced by CLADE amplify the signal of the local characteristics that are present in families composed of very divergent homolog sequences. Hence, the extension of the Pfam domain database by the addition of new probabilistic models (CCMs) makes the approach very sensitive. A proof of this, was shown during the construction of the training sets for the learning step in MetaCLADE, where these models are able to detect much more artificial sequences than the pHMMs. The artificial sequences used to set sensitive thresholds, have a highly biased 4-mer composition generated from the set of positive sequences but they are still completely random as this bias defines only the probability emissions and not the order of the 4-mer. CCM models detects these sequences and allows us to discriminate between the scores which can be achieved just by chance in a highly biased amino acid composition and the true homolog sequences. We shall speak of two factors when significant scores are achieved by random sequences, the first one is the composition of amino acids or n-mers that a random sequence must have in order to achieve a high score, and the second one, the required length to achieve this score. MetaCLADE exploits the second feature based on the hypothesis that true positives are more likely to achieve significant scores more directly (i.e. in less combinations or positions). This is true for sequences which are closer in size and composition to the models. The amount of CCM models introduced in CLADE provides a good coverage of the evolutionary pathways which may be found in a domain family. Hence it provides a good characterization of the different sizes and compositions of the homolog sequences in a domain. By a simple calculation of the mean bitscore, we can distinguish between short hits that are more likely to be positives and negative hits which achieve a high score only by the combination of its amino acid composition and a sufficient large region to produce such a score. This technic, has allowed MetaCLADE to detect highly conserved motifs in very short sequences that are usually dismissed due to trusted threshold which do not take into account the length of the hit.

Giving a critical reading of our tool, it is clear that the extension of the probabilistic model library comes with an increase in the complexity. MetaCLADE was conceived to run in HPC (High Performance Computing) machines as the first step of the pipeline in which it searches for domain hits in the input sequences is computationally expensive. Nevertheless, options like a reduction of redundant models are under work. It would be also important to develop software which could extract specific signatures from models in order to accelerate the search of hits by indexing high scoring words extracted from the models. This step which is done by PSI-BLAST every time a model is used as a query, could be pre-computed. Moreover, we could use this high scoring words of the entire library to detect HSP as we read the input sequences in a process which could be easily parallelized. A list of high scoring words will also unveil significant motifs which may be present in the CLADE library allowing us to improve the quality of the hits by helping to decide the most likely prediction if one of these fingerprints are present.

Due to the increased number of predictions its is important to give significant thresholds for each domain. This problem was tackled in MetaCLADE by the generation of balanced training sets during the learning step. It has provided evidence that it has similar threshold scores as the Pfam database in which the mean GA score was recreated for the pHMM models and therefore we could expect that for most of the CCMs the significant thresholds were correctly fixed. However, the generation of negative sequences is somehow problematic. As seen in chapter 3, some domains need highly biased random sequences in order to produce significant hits for its probabilistic models. With our method, MetaCLADE have assured an exhaustive test of our models through the generation of billion of highly biased artificial sequences. A strategy of meta learning has been used to establish thresholds that separate the better the positive sequences of the negative.

All these approaches have been used in a set of simulated data which has allowed to assess MetaCLADE performance in the prediction of domains. We have given evidence, that the method increases the number of domains identified by other methodologies. Through an enhance domain detection, our tool is capable of amplifying the panorama of functional activities in microbial communities. Its usefulness has been tested in a study aimed at the transcriptional response of marine diatoms to turbulence, annotating more of transcripts than the other methods. The results showed that MetaCLADE has been able to annotate most of the differential expressed sequences under turbulence conditions.

This first version of the tool has many opportunities for improvement by incrementing its predictive power, its run time and its precision. Nevertheless, it is already able to refine protein domain annotations and to provide biologists with an extended landscape of the genetic diversity that is present in the microbial communities. MetaCLADE may provide a powerful tool to classify environmental sequences in current databases and also to retrieve new sequences for which it would be of interest to unveil its evolutionary pathways. To unveil the complexes metabolic processes in communities, it will be essential also, to provide a mapping between our models and databases with well annotated genes that can be directly mapped to metabolic pathways.

MetaCLADE has emerged therefore, as the first multi-source tool wich provides high coverage of the evolutionary pathways of referenced domains in metagenomics and metatranscriptomics samples.

[Combining Machine Learning And Evolution For the Annotation of Metagenomics Data]

Abstract :

Metagenomics is used to study microbial communities by the analyze of DNA extracted directly from environmental samples. It allows to establish a catalog very extended of genes present in the microbial communities. This catalog must be compared against the genes already referenced in the databases in order to find similar sequences and thus determine their function. In the course of this thesis, we have developed MetaCLADE, a new methodology that improves the detection of protein domains already referenced for metagenomic and metatranscriptomic sequences. For the development of MetaCLADE, we modified an annotation system of protein domains that has been developed within the Laboratory of Computational and Quantitative Biology clade called (closer sequences for Annotations Directed by Evolution) [START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF]. In general, the methods for the annotation of protein domains characterize protein domains with probabilistic models. These probabilistic models, called sequence consensus models (SCMs) are built from the alignment of homolog sequences belonging to different phylogenetic clades and they represent the consensus at each position of the alignment. However, when the sequences that form the homolog set are very divergent, the signals of the SCMs become too weak to be identified and therefore the annotation fails. In order to solve this problem of annotation of very divergent domains, we used an approach based on the observation that many of the functional and structural constraints in a protein are not broadly conserved among all species, but they can be found locally in the clades. The approach is therefore to expand the catalog of probabilistic models by creating new models that focus on the specific characteristics of each clade. MetaCLADE, a tool designed with the objective of annotate with precision sequences coming from metagenomics and metatranscriptomics studies uses this library in order to find matches between the models and a database of metagenomic or metatranscriptomic sequences. Then, it uses a pre-computed step for the filtering of the sequences which determine the probability that a prediction is a true hit. This pre-calculated step is a learning process that takes into account the fragmentation of metagenomic sequences to classify them. We have shown that the approach multi source in combination with a strategy of meta-learning taking into account the fragmentation outperforms current methods.

Keywords : [metagenomic, metatranscriptomic, domain annotation, machine learning, protein annotation, genomic, probabilistic model, random sequences, microbial community]

Figure 2

 2 Figure 2 The Glycolysis metabolic pathway. During steps 1-3 two ATPs are used as input. During steps 7-10 2 ATPs are produced, (4 in total as it happens twice per glucose molecule).

Figure 3 DNA

 3 Figure 3 DNA Sequencing Costs: Data from the NHGRI Genome Sequencing Program (GSP) Available at: www.genome.gov/sequencingcostsdata. Accessed [2016].

Figure 4

 4 Figure 4 DNA sequencing library preparation4A.Amplification by emulsion PCR. 4B. Solid phase amplification. 4C. Inmobilization of a single-molecule by a primer 4D Inmobilization of a single-molecule by a template. 4E Inmobilization of a polymerase. http://www.nature.com/nrg/journal/v11/n1/fig_tab/nrg2626_F1.html

Figure 8

 8 Figure 8 The protein structure phases.

Figure 11

 11 Figure[START_REF] Dayhoff | The origin and evolution of protein superfamilies[END_REF] Gumbel Distribution[START_REF] Gumbel | Statistics of extremes[END_REF] Figure12Gumbel Cumulative Distribution[START_REF] Gumbel | Statistics of extremes[END_REF]

= ( 10 )

 10 scores for each position in the sequence. A PSSM assumes independence between the positions of the sequence, because it calculates a score at each position independently of the amino acids at other positions. The scores are represented by positive and negative numbers. They are basically calculated as the log 2 of the observed substitution frequency at a given position divided by the expected substitutions frequency at that position. Positive scores indicate favored substitutions. During the first iteration, the BLASTp algorithm search for similar hits of a query sequence above a certain score or below an e-value score and other parameters such as gap open/extension penalties and substitution matrix. Hits may cover the entire query or just locally Where are the observed frequencies, are the background frequencies and are the target frequencies from the substitution matrix. Next the target frequencies are calculated as follows: Where α and β are relative weights assigned to the observed frequencies and the pseudocount residue frequencies PSI-BLAST raw scores are normalized with equation 6 in order to transform them into bit-scores.HMMER3[START_REF] Eddy | Accelerated Profile HMM Searches[END_REF] is a software suite to detect homology between a sequence and a Profile Hidden Markov Model (pHMM). A pHMM is a probabilistic model that is built from a multiple sequence alignment. A pHMM shows an alignment by means of a sequence of nodes, generally a node per position in alignment and a start and end state. Each node representing a position is composed by three states: correspondence (M), insertion (I) and deletion (D) represented by the mutation rate. The states M model the conserved regions in alignment while the states I and D represent insertions or deletions. The HMM profiles have probabilities on two events: A transition from one state to another and the probability that a state emits a specific amino acid. Only states M and I generate characters. D states are silent. There is therefore, a transmission probability distribution for all states M and I (Figure14). For proteins in specific there are 20 entries, one for each amino acid. The emission probabilities are obtained by amino acid frequency count for each position in the alignment. In the case of a hidden Markov model (HMM), we cannot observe the states directly. However, we do have observations from which we can infer the hidden states. In the case of molecular sequences, the observed states are the positions of amino acids (or nucleotides) in a multiple sequence alignment. The hidden states are the match states, insert states, and delete states. Together, such states define a model for the sequence of that protein or nucleotide family. Different architectures for pHMMs have been proposed. HMMer uses the Plan7[START_REF] Eddy | Accelerated Profile HMM Searches[END_REF] architecture which allows transitions between the end and the start states, inserting a intermediate state which can emit unimportant amino acids.Using this architecture HMMer modelises gapped aligments independently of the score system. that have arisen from a single evolutionary origin into super families called clans. Pfam version 30.0 is composed of 16 306 families and it uses a sequence database called Pfamseq which is based on UniProt release 2016_02 and it was produced at the European Bioinformatics Institute.

2 .

 2 -Combining Model Outputs. CLADE uses a combination of criteria, ultimately converted into a score provided by a Support Vector Machine (SVM), to specifically deal with false positives and eliminate them. This issue is fundamental in domain prediction. The SVM discriminates potential domains by evaluating which prediction is more probable among those displaying: a small E-value, a sufficiently long domain hit, the phylogenetic proximity between the taxon of the sequence to be annotated and the reference species generating the CCMs leading to annotation, and a large agreement among models leading to the prediction. More formally described, given C = {C , … , C , C } the whole of the models characterizing an arbitrary domain D . Given s a sequence query that we want to compare against all models present in C and given C * a subset of C where each C ∈ C * is a model with the best score for a match with a exclusive (non-overlapping) segment of s . For s clade generates x meta-examples (each meta-example contains 5 attributes) where x is the size of C * . To put the emphasis on the individual results of the models, CLADE excerpt three attributes of the output of C : The E-value (1), the length of the hit expressed as the proportion of the model that is covered by the hit (2) and an attribute with a binary value that indicates if the e-value is smaller than a threshold T' (3). To give an indication of the performance of all models, CLADE defines two attributes concerning the percentage of models in C which supports the prediction of C . A model C support the prediction of a model C if: their hits occur in overlapping positions, the percentage of the overlapping region is higher than 50% and if the correspondence with C has a e-value smaller than a threshold T'' (4). The last attribute represents the percentage of models that supports C and which belong to the same phylogenetic clade of s (5).

Figure 18

 18 Figure 18 The SVM Classifier in CLADE. Attributes are calculated from the models outcomes. Then the SVM classifier assigns an score to each prediction.

  3.-Filtered hits are used by DAMA to produce the most likely architecture. Finally, the DAMA software [19] uses the pre-computed list of domain pairs presenting strong cooccurrence in known domain architectures and the list of domain architectures extracted from UniProtKB in order to provide a final prediction.

Figure 20

 20 Figure 20 MetaSim clone length distribution.Figure21FlowSim read length distribution.

Figure 25 .

 25 Figure 25. Length Distribution ORF finder. 5 444 191 open reading frames were identified by the open reading frame (ORF) finder.

Figure 26 Length

 26 Figure 26 Length Distribution MetaGeneMark (MGM).709 627 pCDS were identified by MetaGeneMark

Figure 30 CDS

 30 Figure 30 CDS Fragments vs Open Reading Frame (ORF) Finder. -log10 Evalues of overlapping predictions

  Figures 35 and 36 show the SVM scores obtained from CLADE when applied to a real metatranscriptomics dataset from marine diatoms that will be fully introduced in next chapter. The Y-axis represents the SVM score and the X-axis a -log 10 transformation of the associated e-values. Yellow dots (triangle and cross) that generate a curve from the bottom back to the top, are hits of the domain S-adenosyl-L-homocysteine hydrolase.

Figure 35 SVMFigure 36 A

 3536 Figure 35 SVM Scores of the hits produced by CLADE when applied to a metatranscriptomics dataset of marine diatoms. Points represents identified hits by the CLADE library. Colors represents the average length of the domain model, red [0-100 aa], orange [100-200aa], blue[200-300] and green [>300aa]. The Figure represents the computed coverage, + have a coverage of at least 70% of the sequence and Δ between 50% and 70%

3 MetaCLADE 3 . 1

 331 will introduce MetaCLADE, a tool which combines the enhanced detection of the new probabilistic models introduced by CLADE (CCMs) with a new learning strategy which takes into account the fragmentation of protein domains.Chapter MetaCLADE scopeBiochemical and regulatory pathways have until recently been thought and modelled within one cell type, one organism, one species. This vision is being dramatically changed by the advent of whole microbiome sequencing studies, revealing the role of symbiotic microbial populations in fundamental biochemical functions. The new landscape we face requires the reconstruction of biochemical and regulatory pathways at the community level, as the result of function integration in complex symbiotic interactions of distinct cell types. To reach this complex level of description, a fine domain annotation is of paramount importance, to estimate the abundance of domains involved in the most important functional classes within the community,to recognise the absence of specific domains in an environment, and to highlight enzymatic biochemical functions associated to specific environments. There are no adapted methodologies for this problem nowadays while data are accumulating rapidly. The multisource annotation strategy recently proved to improve the annotation of fully sequenced genomes. Here, we propose MetaCLADE, a novel domain annotation pipeline based on the multi-source domain annotation strategy and designed for metagenomics and metatranscriptomics data. Meta-CLADE exploits multiple probabilistic models representing protein domains to reach high accuracy in the annotation of divergent sequences, and introduces a two dimensional domain specific gathering threshold to handle short domain fragments in reads.

For

  each domain, MetaCLADE estimates bit-score and mean-bit-score domain sensitive thresholds. To do it, it constructs a sequence space for each domain, by defining a set of positive sequences and by generating a set of negative sequences (i.e. sequences wrongly annotated with the domain). Ideally, for each domain, one would like to have a training set comprised of a comparable number of positive sequences and negative sequences. Generation of negative sequences. To define a set of negative sequences for each model (CCM or SCM) associated to a domain, we generate a large amount of decoy sequences and select as negatives those where the original domain is identified by the model (with an Evalue < 1 for CCMs and a positive bit-score for SCMs). The algorithm generates first sequences with two different methods: 1. a random shuffling of the 2-mers of each SEED sequence, 2. the reversal of SEED sequences and checks whether they are negatives or not. If the number of negatives reaches at least the 50% of the positive sequences, then the algorithm stops the search. Otherwise, new sequences are generated with a third method:

Figure 39

 39 Figure 39 Distribution of negative sequences generated or not by Markov models. The impact of the construction of Markov models is illustrated by the distribution of negative sequences in the training sets of CCMs (bottom) and SCMs (pHMMs; top) constructed by either reshuffling of 2-mers or inversion only (A) compared to the generation of negative sequences constructed with the three methodological approaches, i.e. reshuffling, inversion and Markov models (B). Note that about 50% of the domains are characterised by training sets containing at least the 50% of negative sequences.

Figure 3 . 5 A

 35 Figure 40B illustrates the general behavior of the probability spaces for all CCMs and SCMs, all domains confounded. It shows the coherence of the spaces across models of the same domain and highlights bit-scores and mean bit-scores intervals defining rejecting and accepting regions. One observes large regions associated to high probability values accepting true positives

  SwissProt. Namely, predicted CDS overlapping regions containing only SwissProt annotated CDS were taken into account. A total of 65 816 SwissProt annotated Pfam domains were contained in the set of overlapping CDS. MetaCLADE recovers 86.2% of these domains: 56 763 domains fall in the same clan or the same Interpro family, 54 833 of which have the same accession number as SwissProt. The 13.8% of domains that are not detected by MetaCLADE are very small fragments with an average length of 38aa. The distribution of E-values for domains annotated by MetaCLADE is plotted in Figure 41A. (If TrEMBL is used instead of SwissProt, the distribution of E-values is reported in Figure 41B.)

  Figure 42B. See also Figure 42C.) Most importantly, MetaCLADE demonstrated that its

Figure 1 12 Figure 2 15 Figure 3 16 Figure 4 18 Figure 5 19 Figure 6 20 Figure 7 23 Figure 8 25 Figure 9 26 Figure 10 29 Figure 11 30 Figure 12 30 Figure 13 32 Figure 14 34 Figure 15 38 Figure 16 CLADE 41 Figure 17 A 41 Figure 18 43 Figure 19 CLADE 44 Figure 20 47 Figure 21 47 Figure 22 47 Figure 23 47 Figure 24 48 Figure 25 . 49 Figure 26 49 Figure 27 49 Figure 28 49 Figure 29 50 Figure 30 CDS 51 Figure 31 CDS 51 Figure 32 CDS 51 Figure 33 53 Figure 34 53 Figure 35 SVM 55 Figure 36 A 56 Figure 37 62 Figure 38 67 Figure 39 68 Figure 40 70 Figure 41 . 76 Figure 42 78 Figure 43 80 Figure 44 82 Figure 45 91 Figure 46 92 Figure 47 93 Figure 48 94 Figure 49 94 Figure 50 95 Figure 51

 1122153164185196207238259261029113012301332143415381641174118431944204721472247234724482549264927492849295030513151325133533453355536563762386739684070417642784380448245914692479348944994509551 Figure 1 Phylogenetic tree ........................................................................................................12 Figure 2 The Glycolysis metabolic pathway ............................................................................15 Figure 3 DNA Sequencing Costs ..............................................................................................16 Figure 4 DNA sequencing library preparation ..........................................................................18 Figure 5 Single strand nanopore sequencing of a DNA enzyme complex ...............................19 Figure 6 High-powered magnification of an Illumina run ........................................................20 Figure 7 The 21 amino acids grouped by their physiochemical properties ..............................23 Figure 8 The protein structure phases. ......................................................................................25 Figure 9 Multiple Sequence Alignment ....................................................................................26 Figure 10 The BLASTp algorithm ............................................................................................29 Figure 11 Gumbel Distribution .................................................................................................30 Figure 12 Gumbel Cumulative Distribution .............................................................................30 Figure 13 The PSI-BLAST workflow .......................................................................................32 Figure 14 pHMM architecture ..................................................................................................34 Figure 15 The Gene Ontology DAG .........................................................................................38 Figure 16 CLADE Library construction ...................................................................................41 Figure 17 A. Distribution of species in CLADE models ..........................................................41 Figure 18 The SVM Classifier in CLADE ...............................................................................43 Figure 19 CLADE Results for PlasmoDB ................................................................................44 Figure 20 MetaSim clone length distribution ...........................................................................47 Figure 21 FlowSim read length distribution .............................................................................47 Figure 22 Error rate in simulated reads .....................................................................................47 Figure 23 Error position distribution in simulated reads ..........................................................47 Figure 24 Metagenome Simulation ...........................................................................................48 Figure 25. Length Distribution ORF finder ..............................................................................49 Figure 26 Length Distribution MetaGeneMark (MGM) ..........................................................49 Figure 27 Length Distribution FragGeneScan (FGS) ...............................................................49 Figure 28 Length Distribution mapped coding sequences (CDS) ............................................49 Figure 29 Produced datasets from simulation ...........................................................................50 Figure 30 CDS Fragments vs Open Reading Frame (ORF) Finder ..........................................51 Figure 31 CDS Fragments vs MetaGeneMark (MGM) ............................................................51 Figure 32 CDS Fragments vs FragGeneScan (FGS) ................................................................51 Figure 33 CLADE vs HMMer -log 10 E-values of overlapped agreed predictions ...................53

Table 1

 1

	shows a comparative among current NGS

  is another database which is maintained by the NCBI. It is not part of the INSDC but almost all its sequences come from the INSDC database members. RefSeq provides only non-redundant curated data that includes in addition to the sequences, descriptive information, publications and other type of annotations that may come from very different sources. Protein sequence records are translations of the annotated open reading frame.

(Simple Modular Architecture Research Tool) is

  PROSITE[START_REF] Sigrist | PROSITE, a protein domain database for functional characterization and annotation[END_REF] is a database of patterns and profiles characterising protein domains, families and functional sites. Each entry in the PROSITE database is called a signature and it is linked to a document containing different information such as taxonomic occurrence, function, 3D structures and other useful additional information. PROSITE is mainly used to detect precise functional annotation for protein sequences and it is complemented by

	ProRules[86], a set of manually created rules that improves the annotation of PROSITE
	signatures. The latest version of PROSITE contains 1308 patterns, 1039 profiles and 1041
	ProRules.
	PRINTS [87] is a database of conserved motif among protein families and
	structural/functional domains. Protein fingerprints (containing one or more motifs) can be used
	CDD (Conserved Domain Database) [85] is a collection of PSSM models from ancient
	domains and entire protein sequences. Families are divided into sub families only when
	phylogenetic construction of member sequences suggests an origin by gene duplication
	occurring at least 5 billion years in the past. For the alignment of the sub families, it uses

a web-resource for domain identification and the exploration of protein architectures. SMART is composed of manually curated HMMs covering around 1200 different protein domains. SMART database is based on the complete UniProt database combined with predicted proteins from all stable ENSEMBL genomes. After a clusterisation per-species step, SMART database is composed of 1.3 million multiprotein clusters. structural information to improve the alignment. CDD data base also includes models coming from other databases such a Pfam, SMART, TIGRFAM and NCBI databases but it includes only non redundant models. to distinguish sequences at different levels for example, it can be used to distinguish the true domain for sequences having significant hits for more than one domain belonging to the same clan or family. The latest version of PRINTS includes 2156 fingerprint encoding 12 444 individual motifs. PRINTS as PROSITE contains an extensive documentation of each of the fingerprints included in the database. InterPro [88] is one of the most used resources for protein sequence annotation. It is composed by signatures (predictive models characterizing domains, repeats and sites) which are classified into family signatures. Interpro integrates signatures from different databases such as all previous described databases in this sections. It provides also mapping to external functional resources such as Gene Ontology. Intepro scope is to provide a central resource of non redundant signatures, it is powered by InterProScan, a program that allows to combine different signatures to scan input sequences. The latest version of InterPro contains 29608 entries coming from 14 different databases.

2 Performance in simulated metagenomics datasets

  Metagenomics data generated by NGS technologies contains specific properties which must be taken into account. The most relevant of these features is the fragmentation of data.Depending on the sequencing technology, fragmentation can be caused either by the length of the reads or by insertion/deletion errors. As mentioned in chapter 1, the size of reads produced by NGS technologies producing high throughput may vary from ~50 up to 700 nucleotides (nt) while the mean size of a domain is ~300 nt (see Chapter 1). One can expect therefore, to sample only fragments of domains. Insertion/Deletions errors produced by NGS with low throughput but longer reads, result in a constant shift among the reading frames, resulting in the fragmentation of the coding sequences by mixing parts of it with non-coding regions from the other 2 reading frames. As a result, annotation through domain identification in metagenomics

	distinguished a posteriori analyses of the set of species generating the CCMs that helped
	CLADE predictions (Figure 19C). First, we observed that the 54% of the contribution is
	provided by homologs belonging to the Alveolata clade and that the 46% of homologs belongs
	to other clades and, among them, Metazoa, Fungi, Viridiplantae appear to be the most

Domain annotation of all P. falciparum proteins. Over the 5 542 proteins of PlasmoDB

[START_REF] Bahl | PlasmoDB: the Plasmodium genome resource. A database integrating experimental and computational data[END_REF]

, HMMScan identifies 6 037 domains but leaves 2 068 proteins with no identified putative domains. CLADE drastically reduced this number to 1 544, providing 25% improvement and a global annotation of 7 841 domains. In many domain predictions, CLADE exploits CCMs in an exclusive manner: about 87% of CLADE domain predictions obtained with E-value ≤ 1e-60 are contributed by CCMs, and a total of 5 630 domains are identified by CCMs against 2 211 identified by SCMs at e-value ≤ 1e-3. Also, more than a half of the domains predicted by CLADE are co-occurring domains and this is true for all predictions, independently on the E-values (Figure

19A

). To understand how the large space of available sequences is exploited by CLADE to attaint P. falciparum predictions, we performed two represented Eukaryotic clades. A non negligible contribution is also recorded from Viruses, Bacteria and Archaea homologs. (Compare Bacteria with Viridiplantae in the inset of Figure

19C

.) This finding agrees with the intuition that best sequence similarity most likely appears 2.

Table 7 Conserved

 7

motifs used to identify new domains in the O'Connor lake metagenomic dataset

Table 9

 9 8 Transcriptome statistics and differential analyses outcomes

	Mean read length	50 bp
	Unique transcripts	21 224
	Trinity contigs	27 923
	Average length	1426 bp
	Minimum length	201 bp
	Maximum length	11631 bp
	MetaCLADE domains	20 523
	HMMScan domains	15 101
	T2 Up regulated contigs	1 421
	T2 Down regulated contigs	1 107
	T3 Up regulated contigs	192
	T3 Down regulated contigs	272

Table 1

 1 Next Generation Sequencing Summary .......................................................................21

Table 2

 2 Number of annotated domains by CLADE in CDS Fragments ...................................50

Table 3

 3 Annotation summary for CLADE and HMMScan. .....................................................53

Table 4

 4 Ranking scores for selecting hits .................................................................................64

Table 5

 5 Number of CDS (# of CDS) found in reads with FranGeneScan and annotated by different tools ............................................................................................................................81Table 6 Most abundant new domains identified in the O'Connor lake metagenomic dataset .83 Table 7 Conserved motifs ........................................................................................................84 Table 8 Computational time for the analysis of the datasets (using one core) .........................85

Table 9

 9 Transcriptome statistics and differential analyses outcomes .......................................99Ugarte Ari -La combinaison de l'apprentissage statistique et de l'évolution pour l'annotation des données métagénomiques -2016 Résumé : La métagénomique sert à étudier les communautés microbiennes en analysant de l'ADN extrait directement d'échantillons pris dans la nature, elle permet également d'établir un catalogue très étendu des gènes présents dans les communautés microbiennes. Ce catalogue doit être comparé contre les gènes déjà référencés dans les bases des données afin de retrouver des séquences similaires et ainsi déterminer la fonction des séquences qui le composent. Au cours de cette thèse, nous avons développé MetaCLADE, une nouvelle méthodologie qui améliore la détection des domaines protéiques déjà référencés pour des séquences issues des données métagénomiques et métatranscriptomiques. Pour le développement de MetaCLADE, nous avons modifié un système d'annotations de domaines protéiques qui a été développé au sein du Laboratoire de Biologie Computationnelle et Quantitative appelé CLADE (CLoser sequences for Annotations Directed by Evolution)[START_REF] Bernardes | Improvement in Protein Domain Identification Is Reached by Breaking Consensus, with the Agreement of Many Profiles and Domain Cooccurrence[END_REF]. En général les méthodes pour l'annotation de domaines protéiques caractérisent les domaines connus avec des modèles probabilistes. Ces modèles probabilistes, appelés Sequence Consensus Models (SCMs) sont construits à partir d'un alignement des séquences homologues appartenant à différents clades phylogénétiques et ils représentent le consensus à chaque position de l'alignement. Cependant, quand les séquences qui forment l'ensemble des homologues sont très divergentes, les signaux des SCMs deviennent trop faibles pour être identifiés et donc l'annotation échoue. Afin de résoudre ce problème d'annotation de domaines très divergents, nous avons utilisé une approche fondée sur l'observation que beaucoup de contraintes fonctionnelles et structurelles d'une protéine ne sont pas globalement conservées parmi toutes les espèces, mais elles peuvent être conservées localement dans des clades. L'approche consiste donc à élargir le catalogue de modèles probabilistes en créant de nouveaux modèles qui mettent l'accent sur les caractéristiques propres à chaque clade. MetaCLADE, un outil conçu dans l'objectif d'annoter avec précision des séquences issues des expériences métagénomiques et métatranscriptomiques utilise cette libraire afin de trouver des correspondances entre les modèles et une base de données de séquences métagénomiques ou métatranscriptomiques. En suite, il se sert d'une étape pré-calculée pour le filtrage des séquences qui permet de déterminer la probabilité qu'une prédiction soit considérée vraie. Cette étape pré-calculée est un processus d'apprentissage qui prend en compte la fragmentation de séquences métagénomiques pour les classer.Nous avons montré que l'approche multi source en combinaison avec une stratégie de méta apprentissage prenant en compte la fragmentation atteint une très haute performance.Mots clés : [metagenomique, metatranscriptomique, annotation de domaine, apprentissage statistique, annotation de protéine, genomique, modèle probabiliste, sequences aleatoires, communauté microbienne]

Acknowledgment

I have crossed during my internships. I acknowledge the support from the LabEx CALSIMLAB constituting a part of the "Investissements d'Avenir" program for providing financial support for my thesis and for providing access to the HPC resources of the Institute for Scientific Computing and Simulation, Excellence Program "Investissement d'Avenir". I thank Thomas Mock (University of East Anglia, UK) and Angela Falciatore for providing us the samples from [124], and Alberto Amato (Stazione Zoologica Anton Dohrn, Italy) and his collaborators for proposing

Only generated sequences whose original domain has been correctly identified by PSI-Blast (for CCMs) with an E-value < 1 or by HMMer (for SCMs) with a positive bit-score are considered as negative sequences for the MetaCLADE models (CCMs or SCM, respectively) and are included in their training sets. The usage of different threshold for the two systems, PSI-Blast and HMMer, is due to the observation that it is easier to produce negatives with PSI-BLAST than with HMMer; Therefore, an E-value threshold < 1 is much more selective than a positive bit-score. The statistical significance and the impact of these thresholds on the space of positive and negative sequences is discussed below (see "A Naive Bayes classifier sets twodimensional thresholds for fragmented domains").

The algorithm estimates the number of decoy sequences that should be generated to obtain the 50% of negative sequences and stops when this estimated number of sequences is generated. For example, suppose to have 100 positive sequences for a domain, then we seek to generate at least 50 negative sequences. If random reshuffling and reversal generate only 10 negative sequences, a Markov model is expected to generate 40 sequences. Since most decoys generated by the Markov model will not be selected as negatives, one estimates the number of decoys that should be generated by the Markov model to obtain 40 negatives and stops the algorithm after such number. The estimation has been realized based on the observation that false positives are found after a very large number of decoy generations: roughly one expects to obtain 1-10 false positives out of 10 000 decoys for SCMs and out of 1 000 decoys for CCMs.

CCMs lie very close to actual sequences and for this reason we expect them to be much more effective in recognizing a domain in a random sequence generated by a Markov model of that domain than a SCM. If a domain contains only a few sequences in its SEED set, n is too small produce a significant bias in the 4-mer probability emission. Therefore, n is multiplied by a factor W with an initial value of 10 and incremented by one until the domain reaches the same rate of 1-10 negative sequences out of 10 000 decoys for SCMs (We use only SCMs to estimate the weight). Therefore, for domains having a small set of SEED sequences we have emission probabilities × × [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF] where W′ is the final weight obtained. Note that each decoy is tested against both SCM and CCMs associated to the domain, and that it is considered as negative if at least one model identifies the domain in it. If too many negative sequences were produced, then x of those that are the most distant from the origin of the sequence space (this is the Euclidean distance from the two coordinates, bit-score and mean bit-score of a sequence, to the origin of the two-dimensional space; note that sequences that are most distant from the origin are those with higher statistical significance) are retained, where n is the number of Appendix Appendix A Archaeal and bacterial genomes used for generating simulated metagenomic data

List of Figures

List of Tables