Flujo de Trabajo para el Análisis Exhaustivo de Metagenomas
Abstract
El desarrollo de nuevas tecnologías de adquisición de datos ha propiciado una enorme disponibilidad de información en casi todos los campos existentes de la investigación científica, permitiendo a la vez una especialización que resulta en desarrollos software particulares. Con motivo de facilitar al usuario final la obtención de resultados a partir de sus datos, un nuevo paradigma de computación ha surgido con fuerza: los flujos de trabajo automáticos para procesar la información, que han conseguido imponerse gracias al soporte que proporcionan para ensamblar un sistema de procesamiento completo y robusto. La bioinformática es un claro ejemplo donde muchas instituciones ofrecen servicios específicos de procesamiento que, en general, necesitan combinarse para obtener un resultado global. Los ‘gestores de flujos de trabajo’ como Galaxy [1], Swift [2] o Taverna [3] se utilizan para el análisis de datos (entre otros) obtenidos por las nuevas tecnologías de secuenciación del ADN, como Next Generation Sequencing [4], las cuales producen ingentes cantidades de datos en el campos de la genómica, y en particular, metagenómica. La metagenómica estudia las especies presentes en una muestra no cultivada, directamente recolectada del entorno, y los estudios de interés tratan de observar variaciones en la composición de las muestras con objeto de identificar diferencias significativas que correlacionen con características (fenotipo)de los individuos a los que pertenecen las muestras; lo que incluye el análisis funcional de las especies presentes en un metagenoma para comprender las consecuencias derivadas de éstas. Analizar genomas completos ya resulta una tarea importante computacionalmente, por lo que analizar metagenomas en los que no solo está presente el genoma de una especie sino de las varias que conviven en la muestra, resulta una tarea hercúlea. Por ello, el análisis metagenómico requiere algoritmos eficientes capaces de procesar estos datos de forma efectiva y eficiente, en tiempo razonable. Algunas de las dificultades que deben salvarse son (1) el proceso de comparación de muestras contra bases de datos patrón, (2) la asignación (m apping ) de lecturas (r eads ) a genomas mediante estimadores de parecido, (3) los datos procesados suelen ser pesados y necesitan formas de acceso funcionales, (4) la particularidad de cada muestra requiere programas específicos y nuevos para su análisis; (5) la representación visual de resultados ndimensionales para la comprensión y (6) los procesos de verificación de calidad y certidumbre de cada etapa. Para ello presentamos un flujo de trabajo completo pero adaptable, dividido en módulos acoplables y reutilizables mediante estructuras de datos definidas, lo que además permite fácil extensión y customización para satisfacer la demanda de nuevos experimentos.
Full text
2
ESCUELATÉCNICASUPERIORDEINGENIERÍAINFORMÁTICA GRADUADOENINGENIERÍAINFORMÁTICA FlujodeTrabajoparaelAnálisisExhaustivodeMetagenomas ComputationalWorkflowfortheFineGrainedAnalysisof MetagenomicSamples Realizadopor ESTEBANPÉREZWOHLFEIL Tutorizadopor OSWALDOTRELLESSALAZAR Departamento DEPARTAMENTODEARQUITECTURADECOMPUTADORES UNIVERSIDADDEMÁLAGA MÁLAGA,Juniode2016 Fechadefensa: ElSecretariodelTribunal 3
4
RESUMEN El desarrollo de nuevas tecnologías de adquisición de datos ha propiciado una enorme disponibilidad de información en casi todos los campos existentes de la investigación científica, permitiendo a la vez una especialización que resulta en desarrollos software particulares. Con motivo de facilitar al usuario final la obtención de resultados a partir de sus datos, un nuevo paradigma de computación ha surgido con fuerza: los flujos de trabajo automáticos para procesar la información, que han conseguido imponerse gracias al soporte que proporcionan para ensamblar un sistema de procesamiento completo y robusto. La bioinformática es un claro ejemplo donde muchas instituciones ofrecen servicios específicos de procesamiento que, en general, necesitan combinarse para obtener un resultado global. Los ‘gestores de flujos de trabajo’ como Galaxy [1], Swift [2] o Taverna [3] se utilizan para el análisis de datos (entre otros) obtenidos por las nuevas tecnologías de secuenciación del ADN, como Next Generation Sequencing [4], las cuales producen ingentes cantidades de datos en el campos de la genómica, y en particular, metagenómica. La metagenómica estudia las especies presentes en una muestra no cultivada, directamente recolectada del entorno, y los estudios de interés tratan de observar variaciones en la composición de las muestras con objeto de identificar diferencias significativas que correlacionen con características (fenotipo) de los individuos a los que pertenecen las muestras; lo que incluye el análisis funcional de las especies presentes en un metagenoma para comprender las consecuencias derivadas de éstas. Analizar genomas completos ya resulta una tarea importante computacionalmente, por lo que analizar metagenomas en los que no solo está presente el genoma de una especie sino de las varias que conviven en la muestra, resulta una tarea hercúlea. Por ello, el análisis metagenómico requiere algoritmos eficientes capaces de procesar estos datos de forma efectiva y eficiente, en tiempo razonable. Algunas de las dificultades que deben salvarse son (1) el proceso de comparación de muestras contra bases de datos patrón, (2) la asignación (mapping ) de lecturas (reads ) a genomas mediante estimadores de parecido, (3) los datos procesados suelen ser pesados y necesitan formas de acceso funcionales, (4) la particularidad de cada muestra requiere programas específicos y nuevos para su análisis; (5) la representación visual de resultados ndimensionales para la comprensión y (6) los procesos de verificación de calidad y certidumbre de cada etapa. Para ello presentamos un flujo de trabajo completo pero adaptable, dividido en módulos acoplables y reutilizables mediante estructuras de datos definidas, lo que además permite fácil extensión y customización para satisfacer la demandadenuevosexperimentos. Palabrasclave:flujodetrabajo;análisisdemetagenomas;cienciadebigdata;asignaciónde lecturas;genómicacomparativa 5
ABSTRACT The development of new data acquisition systems has promoted an enormous quantity of available data in nearly all scientific fields, thus enabling a more specificoriented research that is resulting in concrete software developments. In the aim of making it simpler for the final user to obtain processed results from raw data, a new computational paradigm is arising: the use of automated workflows for data processing. These have managed to take the lead because of their potential to build complete and reliable processing systems. In particular, the field of bioinformatics is a clear example where lots of worldwide institutions are offering specific services that need to be combined to produce valuable results. Automated workflow managers such as Galaxy, Swift or Taverna are constantly being used for data analysis, especially since the development of new sequencing technologies such as Next Generation Sequencingthatarecapableofproducingvastsamountsofgenomicandmetagenomicdata. Metagenomics aims to study the collection of species present at their original environment. Within the metagenomics field, different types of studies can be carried out; among them, comparative studies of the changes in the composition of samples in the aim of identifying significant differences that correlate phenotype and the individuals, which includes functional analysis of the present species to determine the derived consequences of their presence. While genomics still offers computational difficulties, the science of metagenomics proposes a multiple analysis of unknown species, and therefore requires efficient algorithms capable of processing the enormous amounts of data in reasonable time. Some of the difficulties to be faced are (1) the sequence comparison process with a reference database, (2) the process of mapping reads to species (3) the large size of processed results requires functional storage to allow posterior use in an efficient way, (4) the specificity of each sample demands new processing methods, (5) representation and visualization of ndimensional aspects and (6) qualitycontrolandverificationmethodstoensurethecorrectnessoftheprocess. Therefore we present a complete workflow divided into software modules connected by open and reusable datafiles, which, in conjunction, enables future extension and customization to satisfythedemandofnewexperiments. Keywords:workflow;metagenomeanalysis;bigdatascience;readsmapping;comparative genomics 6
TABLEOFCONTENTS Chapter1:Introduction......................................................................................... 9 Themainreasonsthatmotivatedthedevelopmentoftheproposedworkflow. 1.1Motivation................................................................................................... 9 1.2Projectobjectives........................................................................................ 11 Chapter2:Stateoftheart..................................................................................... 13 Background of the current stateoftheart methods that intent to address the difficultiesposedbythefieldofmetagenomics. 2.1Introduction................................................................................................. 13 2.2Standalonetools.......................................................................................... 14 2.2.1MEGAN.............................................................................................. 14 2.2.2MetaPhlAn.......................................................................................... 16 2.2.3MOCAT............................................................................................... 18 2.3Webbasedtools.......................................................................................... 19 2.3.1MGRAST........................................................................................... 19 2.3.2EBIMetagenomics.............................................................................. 20 Chapter3:Analysisanddesign............................................................................. 23 Analysis of the requisites of the proposed method and formalization of these to fullydefinethefunctionalityofthedevelopedworkflow. 3.1Analysisofgeneralaspects......................................................................... 23 3.2Pipelinedesignandspecifications.............................................................. 24 3.2.1Generalmappingpipeline................................................................... 24 3.2.2Specificfinegrainedpipeline............................................................. 27 3.3Specificationsoftheworkflowmanager..................................................... 29 Chapter4:Implementationandresults............................................................... 31 The system implementation is explained in detail, showing how each independent module solves a particular problem. Results are shown and neatly explained, illustratingtheadvantagesoftheusageoftheproposedmethod. 4.1Introduction................................................................................................. 31 4.2Implementationandresults......................................................................... 35 4.2.1Toolsforthemappingstep.................................................................. 35 4.2.1.1FormatREF........................................................................... 35 4.2.1.2MGRindex............................................................................ 36 4.2.1.3BlastToJoiner,GeckoToJoinerandReverseComplement.... 37 4.2.1.4TaxoMaker............................................................................ 39 4.2.1.5FastaFreqs,SearchSpaceandKLambda............................... 40 7
4.2.1.6FJoiner................................................................................... 41 4.2.1.7Gmap..................................................................................... 43 4.2.2Toolsforthefinegrainedlevelandfurtherprocessing....................... 47 4.2.2.1 QualityMapping................................................................... 47 4.2.2.2GeneParser............................................................................ 48 4.2.2.3CodingAbundances.............................................................. 49 4.2.2.4AnnotationalMapping.......................................................... 51 4.2.2.5GenomeProfileofAccumulatedReads............................... 53 4.2.2.6GenomeCutter,SpecificRegionsandSpecificReads......... 54 4.2.2.7CoverageIdentityMatrix..................................................... 55 4.2.3GalaxyImplementation....................................................................... 57 4.2.4Galaxy’stooldefinitionlanguage........................................................ 57 Chapter5:Usecasesvalidation............................................................................. 61 Verification of the correctness of the workflow is done via utilization of use cases, discussingandcomparingtheresults. 5.1Leanandobesesamples.............................................................................. 61 5.1.1Datasetavailability.............................................................................. 61 5.1.2ComparisonwithMEGAN................................................................. 62 5.1.3Statisticalvalidation............................................................................ 63 5.2Sickandhealthyswinesamples.................................................................. 64 5.2.1Datasetavailability.............................................................................. 64 5.2.2ComparisonwithMEGAN................................................................. 65 5.2.3Statisticalvalidation............................................................................ 66 Chapter6:Conclusions.......................................................................................... 69 Conclusions of the software, future improvements and general comments of the workdeveloped. 6.1Conclusionsofthedevelopedworkflow...................................................... 69 6.1.1.Conclusionesdeltrabajodesarrollado................................................ 70 6.2Furtherextensionsandimprovementsontheworkflow.............................. 71 6.3Acknowledgments........................................................................................ 71 Chapter7:Bibliography......................................................................................... 73 AppendixI:GuidedExerciseusingGalaxy. 8
CHAPTER1. INTRODUCTION 1.1.Motivation TheuseofInternetinthescientificcommunityhaspromotedanenormousquantityof availabledatainalmostanyexistingfield,enablingamorespecificallyorientedresearch. Suchcircumstances,however,areresultinginanincreasingnumberofsoftwareprojects beingdevelopedforparticularscientificandprocessingtasksthatoncefinishedare abandonedsinceitbecomesextremelydifficultfortheaverageusertoruntheparticulartools onhis/herownandthereforegivesuprapidly.Asasideeffect,manyofthesedevelopedtools arebeingincludedasWebServices(inadvanceWS),facilitatingremoteanduniversalaccess andremovingthesomewhattediousstepofdownloadingandinstallingtoolsonlocal machines.However,theseincreasingnumberofofferedWSaretriggeringascenariowhere processingtoolsexistfornearlyeverything,butitisbecomingharderandhardertoreach them,andhasgottentoapointwheretheterm“WebServicesFragmentation”isgaining popularity.Inaddition,suchWSdispersioniscausingworktoberepeatedoverandoverby differentresearchteams,forsearchingforanalreadyexistingWScanbecomparedtolooking foraneedleinahaystack,andthusscientificteamsaredemandingforregistersitesto facilitatethesocalledWSdiscovery.Furthermore,independentWScanbecombinedto producenewdataprocessingnodesthatifnotavailablewouldhaveneedcompletelynew softwaredevelopments.This“combination”oftoolsandservicesarecalledWorkflows, definedas“theorchestratedandrepeatablepatternofbusinessactivityenabledbythe systematicorganizationofresourcesintoprocessesthattransformmaterials,provideservices orprocessinformation ” .Workflowshavebeenusedsince1921byTheRailwayEngineer 1 journal(LawrenceSaunders;S.R.Blundstone,1921)butthetermwasnotusedinthesame senseastodaysinceFrederickTaylorandHenryGantt(knownforhispopularGanttcharts). AutomatedworkflowmanagerssuchasGalaxy,TavernaandMyExperiment[5](among others)wereborntoofferadefinitionlanguagetocombinetools,alongwiththeexecution componentsneededtocompletelyrunthepreviouslydefinedworkflows.Thewideuseof theseautomatedworkflowmanagersmainlymotivatedthedevelopmentofaseriesof independenttoolsandservices,thatputtogetherwouldcomposeafullpipelinewhich,in addition,couldbebuiltandregisteredandusedasaserviceviaautomatedworkflow managers. Additional reasons that motivated the development of this workflow is the exponential growth of available metagenomic (“beyond genomes”) data due to the improvements in the 1https://www.ftb.ca.gov/aboutFTB/Projects/ITSP/BPM_Glossary.pdf 9
Computing. In addition, the complete dependency of the BLAST suite to run the sequence comparison steps (up to MEGAN 4 [19]) makes the software monolithic and slow, since there are new methods to compute homology searches in less time and with almost the same results (e.g. KRAKEN [20], CLARK [21] and GECKO [22]). Furthermore, MEGAN does not offer support for lowabundant genomes and only considers as results of interest those that have a large number of mapped reads, hence discarding any taxonomic result that is not directly supported by large mapped readsabundance. A conceptual drawback of MEGAN is its lack of modularpipeline design. MEGAN software only allows for fixed processing, i.e. it is not possible to contrast results using different sources, it can not be expanded upon wish and necessity, and more important, MEGAN’s individual capabilities can not be combined to produceotherresultsandthereforecannotbereusedwithease. 2.2.2.MetaPhlAn MetaPhlAn[23]isacomputationaltoolforprofilingthecompositionofmicrobial communitiesfrommetagenomicshotgunsequencingdata.MetaPhlAncanuseeitherBLAST orBowtie2[24]toperformasearchofthespecificmarkergenesagainstareferencedatabase. MetaPhlAnisafastmetagenomicabundanceestimationtoolwhichmapsreadstoa setofselectedmarkersequencesthatareuniqueforeachorganisminthedatabase.The markersequencesarecarefullyselectedsuchthatareadcanonlymatchtoonemarkerand canthereforebeassignedtoadistinctorganism.Togetherwiththesmallsizeofthemarker sequencedatabase,thismakesMetaPhlAnveryfastwhiletheaccuracyiscomparabletoother referencebasedmethods.However,sincethegenomesarereducedtoshortmarkersequences, thereisnopossibilitytodetectorevenquantifydifferencesbetweenthesequencedorganism andthereference.Organismsthatcontainmarkersequencesofdifferentreferencestrains(e.g. byhorizontalgenetransfer)showupintheresultsmultipletimesanditisnotpossibleto detectsuchcases. 16
Figure3:MetaPhlAn’sprocessingpipelineoverview.(1)Inputstagetakesametagenomeproducedby shotgunsequencing(e.g.454pyrosequencing,Illumina),(2)theoutputisatableshowingthereads abundanceperspeciesand(3),theendingvisualizationstagethatenablesuserstointeractandexport heatmapsandcladograms,amongothers. IncontrasttoMEGAN,MetaPhlAnisnotdesignedtoworkwithpairedendreadsand thereforethesetypeofreadsshouldbeconcatenatedinonesinglefiletoworkwith. MetaPhlAnhasbeendevelopedunderUNIXenvironmentsinthePythonprogramming language,anditismostlyorientedtoitsuseonshell,thoughitalsoprovidesaGalaxy workflowinstancethatcanbedownloadedandinstalled. MetaPhlAnproducesatabulartextfilewiththereadsabundanceperspecies,along withaseriesofheatmapsandcladogramsthathelpdetectsignificantvariationbetween speciesandunderstandtheassignationofreadstotaxa,respectively. Figure4:CladogramproducedbyMetaPhlAnshowingtheassignationofreadstogenomesbyspecies, phylumandfamily.Althoughthesetypeofgraphsallowtotakeafirstglanceatthetaxonomical results,itbecomeshardtopreciselydeterminetheabundancedifferencesamongphylum. However, since MetaPhlAn only classifies a subset of reads that map to one of its marker genes, its results are only estimations of the species that populate the sample and a largeportionoftheDNAmaterialisunused,andthereforebecomesunusedinformation. In a field as effervescent as metagenomics, which is steadily expanding itself with new experiments and knowledge extracting methods, it becomes necessary to enable a handy, easy extension of the capabilities of the analysis software, thus avoiding developing functional 17
packages that will not be used due to new requested experiments. In this line, MetaPhlAn offers a fully functional pipeline design which, in contrast with MEGAN, can be extrapolated to reuse in newer software developments and integrated into larger webservers that offer metagenomicsanalysis. 2.2.3.MOCAT MOCAT [25] is a software pipeline for metagenomic sequence assembly and gene prediction for taxonomic and functional abundance profiling. Since it is a gene prediction tool, it searches for single copies of particular marked genes across different databases. The automated generation and efficient annotation of nonredundant reference catalogs by propagating precomputed assignments from 18 databases covering various functional categoriesallowsforfastandcomprehensivefunctionalcharacterizationofmetagenomes. The second version of MOCAT released supports Illumina single and pairedend reads in raw FastQ format. MOCAT2 [26] generates taxonomic and functional profiles, as wellasassemblereadsandpredictgenesinassembledsequences. Figure5:MOCAT’sprocessingpipeline.Thetopboxinorangerepresentsthedatagathering stepfromtheenvironment,whichmustbedonepriortorunningtheworkflow.Inlightpurple,the MOCATpipelineshowsthatsomestepscanbebothskippedorransequentially,e.g.onlytheread 18
trimming,assemblyandgenepredictionstepsaremandatorytoproduceanestimationofthe taxonomicalresults.Ontheotherhand,readalignmentandremoval,alongwiththeassemblyrevision, areoptionalandshouldbedoneforqualityimprovementifdesired. Sequencedrawreadsarequalitycontrolledbyremovinglowqualityreads,e.g.reads arealignedtothehumangenomeusingSOAP[27],andthenthesereadsareremovedsincein mostcasestheirpresenceisduetocontaminationofthesample.Thisaprioristepofaligning readscanbeusedasanestimationofthereadsabundances.Theremainingreadsare consideredtobeofhighquality,andareassembledintocontigswhicharelaterusedto performthegenepredictionstepusingProdigal[28]orMetaGeneMark[29]. MOCATisimplementedinPerl,andsincemanyofthetoolsusedinthepipelineare fromthirdpartydevelopers,itisrequiredtomanuallyobtaintheselicenses,forexamplein thecaseofMetaGeneMark.Furthermore,MOCATlacksasimpleinterfaceimplementationto helpusersruntheirexperiments,andisexecutedviacommandline,whichitmakesextremely hardforresearcherswhosefieldofknowledgeandapplicationisnotcomputersciences,e.g. biologistsandgeneticists.Inaddition,MOCATisanexampleofmetagenomicanalyser packagethatusesonlyportionsoftheoriginalmetagenomes(forexamplebylookingfor specifichighlyconservedgenes,i.e.genecalling)and/oragainstproteindatabasesbecause supposedlyonlycodingregionshaveafunctionalrole,andthatthereforenoncodingregions canbediscarded.Furthermore,suchprocedureremovesahighpercentageofdatatocompare andprocess,thusspeedinguptheanalysisprocessbutinexchangeoflosingDNAmaterial 2.3.Webbasedtools 2.3.1.MGRAST ThemetagenomicsRASTserver[30](MGRAST)isaSEEDbasedenvironmentthat allowsuserstouploadmetagenomesforautomatedanalyses.MGRASTprovidesthe annotationofsequencefragments,theirphylogeneticclassification,functionalclassification ofsamples,andcomparisonbetweenmultiplemetagenomes.Inaddition,italsocomputesan initialmetabolicreconstructionforthemetagenomeandallowscomparisonofmetabolic reconstructionsofmetagenomesandgenomes.Twoapproachesareusedtoperformthe taxonomicalclassificationofreads:(1)MGRASTbuildsclustersofproteinsatagiven percentageofidentitylevelusingQIIME[31],andthenthelongestsequenceofeachcluster iscomparedusingidentityasfilterinadatabasesearchviasBLAT,animplementationofthe BLAT algorithm;(2)performinggenepredictionbysearchingforribosomalRNAsagainsta 11 nonredundantintegrationoftheSILVA[32],Greengenes[33]andRDP[34]databases. 11http://genome.ucsc.edu/cgibin/hgBlat 19
MGRASTalsousestheNCBItaxonomytoperformthetaxonomicalclassification,in similitudetoMEGAN.Thefunctionalprofilesareavailablethroughcomparisonwithdata sourcesthatprovidehierarchicalinformation.Abundanceprofilesarethemainoutputfor displayinginformationonthedatasets.TheMGRASTannotationpipelineusuallydoesnot provideasingleannotationforeachsubmittedfragmentofDNA. Figure6:MGRAST’sprocessingpipelineshowsfromqualitycontrol(1,2,3and4)to taxonomicalclassificationviageneprediction(5and6)or/andRNAdetection(8,9).Thepipelinecan beexecutedusingoneortwoways:(1)toplayer,usinggenepredictionandproteinidentificationor (2)usingRNAdetection.Whilethefirstonelooksforannotatingreadsbyspecifichighlyconserved regions(genes)theotherlooksforannotationviatranscriptomesdetection. However,dataprivacyisoneoftheconcernsofthescientistsusingthistool:firstly, theyareworriedaboutuploadingtheirunpublishedand/orconfidentialdatatoapublic websiteandsecondly,thepriorityofjobanalysissubmittedtosuchwebsitearesubjectedto theconfidentialityleveloftheinputdata(withlowerpriorityandthereforelongerwaiting timesforprivatedata)andsizeofthesubmitteddata.Forexample,forshotgunmetagenomes themedianwaitingexpectedtimeis7–10days,whereasforamplicondatasets,thepipeline usuallyclearsdatasetswithin24h,sincetheseareofsmallersizeingeneral. 2.3.2.EBIMETAGENOMICS EBI Metagenomics [35] (EMG in advance) has been developed as a freetouse, largescale analysis platform for metagenomic sequence data. The resource is capable of processing sequenced metagenomic and metatranscriptomic reads, 16S rRNA amplicon data and usersubmitted sequence assemblies. Regardless of the data source (metagenomic, metatranscriptomic, amplicon or assembly), EMG provides a standardized analysis workflow, capable of producing taxonomic diversity and functional annotations. As a result, analyses 20
can be compared both within and across projects at a broad level, and across different data types(e.g.metagenomicversusmetatranscriptomic). The EMG pipeline is a combination of existing tools that are widely used in quality control, feature prediction and taxonomic and functional prediction, such as TIGRFAMs [36] or BioPython [37]. Rather than developing an entirely new repository for metagenomic data, EMG has partnered with the European Nucleotide Archive [38] to provide a permanent metagenomics data archive and data sharing/publication service. ENA's existing infrastructure and interfaces for data submission ensures that datasets are described with detail and contextual data, and are made available to the scientific community for data mining purposes andformetaanalysis. 21
Figure7:EMG’sprocessingpipelineshowtwoprocessingwaysafterthetraditionalquality controlsteps.Firstly,thequalitycontrolstepisperformedusingTrimmomatic[39].Twopossible processwaysareofferedbyEMG:(1)leftside,taxonomicreadsabundanceestimationbasedon16S rRNAsgenes,and(2)functionalassignmentbasedoncodingregions(CDS). Once reads (of either metagenomic or amplicon origin) are uploaded, they are queued for analysis on the pipeline. The waiting times depend mostly on size and if the uploaded data is correctly annotated, fulfilling with the standards, as one of the key points in EMG is the possibility of future users to run their own experiments with such data. The analysis pipeline consists of quality control, ribosomal RNA (using rRNAselector [40]) and protein prediction (via FragGeneScan [41]), function prediction (using InterProScan 5 [42]) and taxonomic prediction steps, which uses QIIME. The EMG pipeline is implemented under a Python frameworkthatmanagestheexecutionofstepsandjobswithinqueuesonaclustercomputer. 22
CHAPTER3. ANALYSISANDDESIGN In this section we will determine which requisites should be fulfilled by the proposed method. These include, among others, how should inputs be treated, the format of the output files, which times are reasonable for particular tasks, how should the workflow be accessed and executed, etc. Some aspects will be of obliged need due to the nature of the metagenomics analysis itself. We will reflect below the top level requisites and leave low level and detailed issues to be addressed later on the implementation chapter. In addition, specificationswillbegiventoprovideafirstgrasponthecapabilitiesoftheproposedmethod. 3.1.Analysisofgeneralaspects 1. Platformusage: a. Mostbioinformaticssoftware(andingeneral,DataSciencerelated)is designedtobeusedunderUNIXenvironments.Thisisduetothepowerful nativeshell,theallaspectscustomization,andthatUNIXisnotonlymeantto beusedasaclientoperatingsystem(suchasWindowsorMac)butalsoasa developerplatform.Therefore,theproposedmethodwillbedevelopedtouse underUNIXoperatingsystems,suchasUbuntu orDebian . 12 13 2. Programminglanguage(s): a. Metagenomicsprocessingdemandshighmemoryusage(e.g.loadingwhole genomedatabasesintomemory,suchasdonebyKRAKEN)aswellasCPU time(e.g.aligninglargesequencescanbeupto foreachpairof(n)O2 sequencesdependingonthealgorithmused)andthereforerequirealowlevel programminglanguagethatenablesthedevelopertocorrectlymanagehow memoryisallocatedandreleased,andhowCPUcyclesarespent.Thus,theC programminglanguagehasbeenchosen,sinceitisalowlevellanguage almostasfastandeffectiveaswritingcodedirectlyonmachinecode. b. Plotting,graphsandchartswillbeproducedusingtheRprogramming language,sinceitoffersahighlevelabstraction,acompletesuiteofstatistical packagestosimplifycodingtasksandamodularorientedimplementation,in termsofdirectshellexecutionunderUNIXsystems. 3. Typeofmetagenomicanalysissoftware: a. Acrossallpossiblemetagenomicanalysiscategories(e.g.functional annotation,comparativemetagenomics,assembly,etc)wewillorientthe 12http://www.ubuntu.com/ 13https://www.debian.org/ 23
proposedmethodtowardsasequencesimilaritysearch,thatis,performa comparisonofallreadsagainstareferencedatabaseandmapthereadsbased onthepropertiesofthealignedfragments(e.g.similarityandlengthofthe match).Asmentionedbefore,thisisanexpensivecomputationalprocessand thereforerequirespackagestoacceleratethepreprocessing,comparison, parsingandmappingprocedures. 4. Resultsflexibility: a. Toenableusingdifferenthomologysearchalgorithmsandinordernottoforce usingaparticularone,allresultsproducedbysuchalgorithmswillbeparsedto aconcreteformat.SeetheImplementationchapterforthespecificformat. 5. Workflowdivision: a. Thewholeworkflowwillbedividedintwo: i. Ageneralmappingworkflowthatwillbeonchargeofmappingreads totaxa. ii. Aspecificworkflowforpostprocessingresultsatafinegrainedlevel. 3.2.Pipelinedesignandspecifications Asstatedearlier,thesoftwarewillbedevelopedasaworkflow,i.e.apipelineof singlemodulesthattakeeitherrawdataorprocesseddatabyapreviousmoduleandproduce thenextlevelofprocesseddata.Thefollowinglistcontainsthedifferentmodulesthatwillbe included.LowerleveldefinitionsandproceduredetailsaregivenintheImplementation section. 3.2.1.Generalmappingpipeline Userinputsofthegeneralmappingpipeline: 1. AcustomdatabaseofnucleotidesinFASTAformat,concatenatedinasinglefile. 2. AmetagenomeinFASTAformat.Ifworkingwithpairedendreads,theseshallbe concatenatedinasinglefile,thusconsideringthemasindependentreads, 3. Acuadrangularscoringmatrixforeachnucleotidebyrowsandcolumnstohandle scoringalignments. Listedbelowarethemodulesthatwillcomposetheworkflowwithinputsandoutputs. Abriefdescriptionwillbegivenalongwithtwoseparatedlines,onelabeledwithIN for Inputs andOUT forOutputs . 24
1. SeqTrimNext[43]:Athirdpartytoolthatwillbeusedtotrimreadsproducedby 454pyrosequencing. IN: Theinputreadsfile. OUT: Thetrimmedreadsfile. 2. Trimmomatic:Anotherthirdpartyqualitycontroltoolthatwillbeusedtoprocess IlluminaHiSeq2000reads. IN: Theinputreadsfile. OUT: Thetrimmedreadsfile. 3. FormatREF:Amodulethatwilladdanidentifiertoeachsequenceofthedatabase usingspecialcharacterstodealwithdifferentFASTAheaders. IN: Theinputdatabase. OUT: Theformatteddatabase. 4. MGRindex:AtooltoproduceanindexforaFASTAfilecontaininguseful informationsuchasstartingpositions,lengthandmaskedresidues. IN: Readsfile. OUT: Twobinaryindexesthatcontainastructureperreadwiththereadparticular information,oneinthesameorderastheinputreadsfile,andtheotherinsortedorder toenablebinarysearch. 5. SearchSpace:Atooltocalculatethenumberofnucleotidesinafile. IN: Formattedinputdatabase. OUT: Atextfilecontainingthenumberofresidues. 6. FastaFreqs:Aprogramthatwillcalculatethefrequencypernucleotideofthe databasetosupportthecalculationofKarlinandLambdaparameters. IN: Theinputformatteddatabase. OUT: Afilecontainingthefrequenciesperresiduewiththesumofthefrequencies being1. 7. TaxoMaker:Aprogramthatproducesthetaxonomyfileofthedatabasealongwith otherinformation,suchasgenomeslengthandidentifier. IN: Formattedinputdatabase. OUT: Tabulatedtextfilecontainingdatapereachgenomeinthedatabase. 8. GECKO:AnexternalprogramtocomputehomologysearchesbetweentwoFASTA files. IN: Inputtrimmedmetagenome,Inputformatteddatabase OUT: Afragmentsbinaryfilewhichcontainsthealignmentsofeveryfragment reported. 9. ReverseComplement:AmodulethatproducesthereversecomplementofaFASTA file. IN: Formattedinputdatabase. OUT: Thesameinputdatabasebutwithitssequencesbeingcomplementedand reversed.Inaddition,sequenceswillbeinreverseorderaswell. 25
Figure8:Workflowdiagramdividedintwosubsystemsandapriorqualitycontrollayer.Thered enclosedareaisthegeneralmappingworkflow,whereasthegreenenclosedareaisthe specificgenomelevelworkflow.Thetopside,beigeenclosedareaisthepreprocessingstepwhich involvesqualitycontrolofmetagenomes. AlthoughdevelopedasapipelineofCmodulestobecalledsequentiallyusingUNIX’ shell,itseemedthattheaverageuserwouldfindtroubleinrecurrentlyexecutingscriptsone afteranother,andthatsuchdoingcouldpotentiallyintroduceerrorsinthepipelineprocessing. Furthermore,evenattestingstages,thedevelopermyself,wouldcommiterrorsduetothe largenumberofmodulestobecalledinordertoanalysesinglemetagenomes.Suchreasons promotedtheuseofaworkflowmanager,andforthatgoal,itwasdecidedtobetheGalaxy workflowmanager,foritnotonlyallowedthemanagementofjobsandqueueswithinthe 32
workflow,butalsoenabledtorunthesoftwarebothinlocalmachinesandonlineinstances. Seethesection4.4.0.fordetailsontheGalaxyimplementation,whichisalreadyinstalled withtheworkflowitselfintheprovidedvirtualmachine. InordertofitthepipelinerepresentationinA4sheetsandthereforeallowprintingthe document,thegeneralworkflowhasbeencutintwopiecesfromleft(startoftheworkflow) toright(endoftheworkflow).SeeFigure9andFigure10. Figure9:Leftside(inputsandprocessingstart)ofthegeneralworkflow.Theworkflowcanvasshown istheillustrationoftheconnectionbetweenmodulesinputsandoutputs.Itbecomesclearerthat manuallycallingeachmodulecouldproducepotentialhumanintroducedbugs. 33
Figure10:Rightside(processingendandvisualizationoutputs)ofthegeneralworkflow. Starmarkedtooloutputsareproducedresultsthatwillbeconserveduponfinishingtheexecution. On the other hand, the specific workflow (See Figure 11) does fit in a single A4 sheet. Whereas the general mapping workflow provides the tools needed to perform a basic taxonomic classification of metagenomic samples from its very beginning (i.e. preprocessing and homology search), the specific workflow is aimed at a posteriorfinegrained step. It is run with the results from the general workflow and two additional files, namely the genome for which to perform the specific execution in FASTA format and the annotation file from the GenBankdatabase. 34
Figure11:Fulloverviewofthespecificworkflow,whichproducesresultsatagenomelevelusingthe outputfromthegeneralmappingprocess.Theleftlayeristhedatainputmodules,thecenterlayer comprisetheprocessingtools,whereasthelastlayeraremainlyoutputresultsandvisualizationplots. 4.2.Implementationandresults Inthissectionwewillbrieflydiscussthemodulesthatcomposetheproposedmethod, pointingoutwhichissueswerefoundandhowtheywereaddressed.Inaddition,therunning timeandcomplexityintermsoftheinputwillbeprovided.Wewilldifferentiatebetweenthe generalprocessingtoolsthatinsomeaspectareresponsibleforthemapping,andthe postprocessingfinegrainedleveltoolsthatusethemappingfilestoproducefurtherresults. 4.2.1Toolsforthemappingstep 4.2.1.1.FormatREF ThissmallprogramwasincludedtointernallyassignanIDtoeachgenomeinthe databasethroughouttheanalysisprocess.Thereasonbehindthisadditionisthedisagreement 35
ontheidentifiersusedinFASTAfilestouniquelyidentifygenomes.Althoughdatabasessuch asGenBankprovidewithGIid’sandaccessionnumbers,itisnotpossibletoknowiftheuser whichhasmadeacustomselectionofspecieshas(bychanceorwill)modifiedthesenumbers orifheisusinganoldsetofgenomeswithoutdatedreferences.Thereforeitseemed appropriatetoincludeasimplenumbertouniquelyidentifyagenome. Theprogramitselfproducesacopyofthedatabaseaddinganumbertoeachgenome, whichdependsontheorderofthegenomesinthefile.TheFASTAheaderstartswitha“>” symbolandisfollowedbyanidentifierandoptionaldescriptions,uptoabreakline.The numberisaddedbetweentheopening“>”symbolandtherestoftheline,e.g: PriortoFormatREF AfterFormatREF >NC_004307.2Bifidobacteriumlongum NCC2705chromosome >|k|NC_004307.2Bifidobacterium longumNCC2705chromosome Table1:ChangesappliedbyFormatREFaretheinclusionofasimpleidentifierthatwillavoiddealing withdifferenttypesofFASTAheaders. The complexity of FormatREF is linear in the size of the input (the database) and therefore is , for it sweeps the database through copying characters and writing them to (n)O anotherfile. 4.2.1.2.MGRindex MetagenomestendtobelargeFASTAorFASTQfiles,ofnormallyhundredsand hundredsofmegabytes(ifnotmore).Alongtheexecutionoftheworkflow,itwillbeneeded toloadreadsand/orpropertiesaboutthem,e.g.thelengthofaparticularreadtocomputethe percentageofidentity,coverage,etc.Furthermore,linearlyscanningametagenomefileevery timeaparticularreadorfieldisneededwouldbeslow,redundantandabsolutelyintractable. Tofulfillthisrequirement,twoindexesarecreated:onethatisnotsortedandservesto holdproperties;andasortedonetoallowlogarithmictimeaccessusingbinarysearch.The storedpropertiesbytheseindexesareshowninthetablebelow: Field Description id Theidentifierthatcomesafterthe“>”FASTAheader startingsymbol. rNumber Readnumbercountingfrom0. 36
rLen Lengthofthereadinnucleotides. rLmasked Lengthofmaskednucleotides. nonACGT Numberoflettersthatarenotnucleotidese.g.the“N”is usedtorepresentanunknownnucleotide. pos Positionofthereadinthemetagenomefile. Lac Accumulatedlengthofallreadsthatappearpriorinthe file. Table2:Datathatwillbestoredinbothindexes.Thisinformationwillbeparticularlyusefulwhen translatingthecoordinatesoffragmentsofthesequencecomparisonsoftwareGECKO,whichyields suchcoordinatesasanaccumulatedsumofthelengthofallprevioussequences. Thefieldsarestoredasabinarystructure.Inthecaseofthesortedindex,aQuickSort [41]algorithmisusedtosortin time.Therefore,whenareadidentifierisfound,a(nlogn)O binarysearchisperformedtoobtainthebinarystructurethatbelongstoit.TheMGRindex programtakesonebyoneallreadsandfirstlystoresthem(unsortedbinaryfile),whichtakes being thenumberofreads.Afterwards,suchfileissortedwithQuickSortandstored(n)O n againinadifferentfile.Therefore,thecomplexityisthecompositionofthetwocomplexities independently,andmakeitupto whichisstilllinear.(nnlogn)O+ 4.2.1.3.BlastToJoiner,GeckoToJoinerandReverseComplement Forbothsequencecomparisonsoftware(GECKOandBLAST),theresultswillneed tobeparsedtoaformatthatisusablebytheworkflow.ThisisdoneusingGeckoToJoinerand BlastToJoiner,respectivelyofthesoftwareused.Thetwoarequicklinearparsersthatstore allrelevantinformationaboutthereportedfragments,suchasthematchinglength,the numberofidentitiesortheinsertedgaps.Alignments,however,arenotstoredsincetheyare nolongernecessarytoperformthemapping;onlytheirpropertiesareneededtodecidewhich arethebestcandidatesforeveryread. Theformatconsistsinatabbasedtextfilecomposedofaheaderwithoneormore 12tuples. t12 n,k(k,score,dentities,length,similarity,igaps,egaps,strand,rStart,rEnd,gStart,gEnd)= i Whereeachfieldcorrespondsto: Field Description Fragmentnumber Foreveryread,thepositioninthelistofreportedfragments 37
Score Reportedscoreforthefragment Identities Numberofidentities Matchinglength Lengthofthematch %Similarity Thenumberofidentitiesdividedbythematchinglength iGaps Numberofopeninggaps eGaps Numberofextensiongaps Strand Strandofthefragmentrepresentedwith++or+ Startinread Startingpositionofthematchinthereadcoordinates Endinread Endingpositionofthematchinthereadcoordinates Startingenome Startingpositionofthematchinthegenomecoordinates Endingenome Endingpositionofthematchinthegenomecoordinates Table3:Fieldsoftheparsingformatthatisusedintheworkflow.Inordertouseanysequence comparisonsoftware,itwillbeneededtoprovideaparsertotheformatshownabove. Multiplefragmentsindifferentgenomescanbereportedforthesameread.Therefore, onetomanyfragmentswillbeassociatedtoaheader,oneforeachuniquecombinationof read andgenome, alongwithlength : >ReadID >Genomeaccessnumberandname Genomelength Table4:Headerinformationthatrepresentsatuplereadgenome.Eachtupleisauniquecombination ofareadandagenomeandcanonlyappearonce. AnRscriptisusedtoplottheaccumulateddistributionoffragmentsbypercentageof identityandcoverageinlog10scale.SeeFigure12. Inaddition,inthecaseofusingGECKOasthesequencecomparisonsoftware,and thereforeGeckoToJoinerasparser,itbecomesnecessarytoproducethesortedandunsorted indexesforthereversecomplementofthedatabase,sincethecoordinatesoffragments locatedinthereversecomplementarysequencegivenbyGECKOneedtobetransformedtoa forwarddatabasesystem.Tocreatesuchindexes,firstwewillneedtocreatethereverse complementofthedatabaseandinreverseorder.TheprogramReverseComplementdoesthis byfirstlysweepingthedatabasefileandstoringthestartingpositionsofeachsequence(say )inavector.Thenthefilepointerispositionedoneachpositionstartinginthe andN Nth downtothe ,loadingthesequenceintomemoryandwritingitinreverseorderwhile1st complementingitatthesametime.Thusifthedatabaseisofsize ,afirstsweepisneeded,n 38
plusanothersweepforeverysequencetoloadintomemoryandalastonetowriteand complement,hencemakingthecomplexitytobe whichis .(n)O+n+n(3n)O Figure12:This3Dplotshowstheaccumulatednumberoffragments(Zaxis)inlogarithmicscaleby percentageofidentity(Yaxis)andpercentageofcoverage(Xaxis).Spiked(x,y)pointsrepresent regionsthathaveahighernumberofreportedfragments. BothGeckoToJoinerandBlastToJoinerruninlineartimedependingonthesizeofthe input.Theseprogramsarepreparedtorunextremelyfastsincetheinputfilescanbeofsizeup toseveralhundredgigabytes.Thecomplexityis forboth,being thenumberof(nlogn)O n fragmentspresentinthecomparisonfile.Theuseoftheindexesdiscussedinsection4.2.1 enablethebinarysearchtoretrievereadpropertiesinlogarithmictime. 4.2.1.4.TaxoMaker Theworkflowusesataxonomyfilethatmapsid’stogenomesandstoresthenames andlengthofthegenomesinthedatabase.Theseareusedforvisualreasons,i.e.showlabeled axisforgenomesinplotsorfiles,whichisalmostamusthavewhendealingwithlarge databases.Thelengthisusedtocalculatetheexpectedvalueofeachfragmentsinalaterstep intheworkflow.Thesetaxonomyfilesarealsousedtoincorporatecustomrelationships 39
betweenthegenomes,aswellascreatingphylumsofsimilargenomes,e.g.differentstrainsof MycoplasmaHyopneumoniaecouldbegroupedtogether.Seebelowanexampleoftaxonomy file: SpecieLevel 1 SpecieLevel 2 Accesscode SpecieName 1 0 ref|NC_021831.1| Mycoplasmahyopneumoniae J 1 1 ref|NC_017509.1| Mycoplasmahyopneumoniae 168 2 0 ref|NC_015144.1| WeeksellavirosaDSM 16922 3 0 ref|NZ_DS264342.1| RuminococcusobeumATCC 29174Scfld0253 3 1 ref|NZ_DS264341.1| RuminococcusobeumATCC 29174Scfld0254 3 2 ref|NZ_DS264340.1| RuminococcusobeumATCC 29174Scfld0255 Table5:Exampleoftaxonomicaldescriptionfilewheretwospeciesareboundedtogetheras substrainsusingthesamenumberforSpecieLevel1,asecondspeciehasnosubstrainattachedanda thirdspeciecomposedofthreecontigs. Thesefilescanbebuildwithanytexteditororspreadsheetsoftware.However,itcan becomeverytedioustobuildthesefileswhenworkingwithlargedatabaseswhereseveral genomesarepresent,andthereforethetoolTaxoMakerproducesthesetypeoffields automatically.Itrunsinlineartimesincetheprocedurefollowedistoparsethedatabasefile andaddnewentriestothetaxonomyfileforeachgenomefound. 4.2.1.5.FastaFreqs,SearchSpaceandKLambda Thesetoolsaretreatedasonce,sincethethreeareexecutedtogetherandare completelydependant.Thegoalofthesetoolsistoextractaseriesofparametersneededto latercomputetheexpectedprobabilityofreportingaparticularfragmentgiventhelengthand propertiesofboththequerymetagenomeandthedatabase.FastaFreqsproducesafile containingthefrequenciesperresidue,andshouldbeusedonthedatabase,alongwith searchSpace,whichcomputesthetotalnumberofresiduesinthedatabase.Inparticular,the searchSpaceprogramisnomorethanacombinationofwellknownunixscripts: 40
grepv“>”database.fasta|wc Whichproducesthetotalnumberofresidueswhencomputingthedifferencebetween thethirdcolumn(numberofcharacters)andthefirstcolumn(numberofnewlines),thus removing\ncharacters. TheKLambdatoolproducestheKandLambdaparameterasstudiedbyKarlinand Atschul[45]givenascoringmatrixofnucleotidesandthefrequenciesperresidue.These parametersarethenusedtocomputetheexpectedprobabilityofafragmentinthemapping moduletochoosethebestfragment. 4.1.2.6.Fjoiner Metagenomesareunculturedsamples,andthebacteriapresentinsuchenvironmental communitiessufferfromhighevolutionarypressureduetoallinteractionsbetweenspecies, causingmostlysmallmutations,insertions,deletions,translocations,etc.Suchisthemain reasonbehindthedevelopmentofthetoolFjoiner,whichtriestoextendfragmentsby performingaNeedlemanWunschalignmentbetweenthosethatbelongtothesamereadand genome.Alignedfragmentsthatsurpassidentityandcoveragethresholdsarekeptandstored. ThereforeFjoinerproduceslongerfragments,facilitatingthedecisionstepofthebest candidateforeveryread.GenomespresentinthedatabaseareloadedintoRAMtoaccelerate theprocess,sincesequenceretrievalfromfilesbecomessloweraslargerastheseare.Then thecomparisonfileisreadsequentiallyloadingallfragmentsbelongingtoatuple ReadGenome .Everypairofcandidatefragmentsaretriedtobealignedusingacustom tworowsNeedlemanWunschalgorithm.Tablexshowsthepseudoalgorithm: 1.Loadgenomesintomemory 2.Quicksortgenomesbyaccessionnumber 3.Loadreadsindex 4.Foreveryreadinthecomparisonfile a.Sortall fragmentsbelongingtothereadbyreadk coordinates b.Foralltuplesofconsecutivefragments( whilef)fnn+1 kn− 1 > i. If( distanceinbothgenomeandreadf)fnn+1 coordinatessatisfythethreshold 1.Loadthecompletereadusingbinarysearch 2.PerformNeedlemanWunschscoringmatrix between( and(.rf.r)fn1n+1 2 .gf.g)fn1n+1 2 3.Iftheobtainedalignmentsatisfiesminimum coverageandsimilaritythresholds: 41
Figure20:Optionscomparisonplot.TheXaxisarethegenomesbyspeciesandtheYaxisisreads abundance(topleft),averagepercentageofcoverage(topright),averagepercentageofidentity (bottomleft)andtheaveragelengthinnucleotides(bottomright). Figure 20 shows an optionscomparison plot between the first best candidate (changing color in the plot) for all genomes and the second candidate (blue line in all plots). For example, the topleft plot shows that only the genomes with ID 4 and 6 have significant distance in terms of mapping abundance respect to their second options, thus ensuring a better mapping. It can be observed in the rest of the plots that their average coverage (topright), average identities (bottomleft) and average length (bottomright) is superior in all cases to thatoftheirsecondbestoptions. The QualityMapper tool is linear on the size of the input, since its working principle is sweeping the whole binary mapping file adding up the values that conform the different statisticindicatorsshown. 4.2.2.2.GeneParser Althoughthistooldoesnotusethemappingfile,itseemedmoreconvenienttoputit underthepostmappingtools,sinceitisofnousewithouttheseinfurthersteps.GeneParser wasdevelopedtoparseGenBank’sannotationfilestokeeponlytheneededfields,sincethese 48
typeoffilesoftentendtobelarge,andfordatabasesofconsiderablesizeitcanbecomealot ofunusedinformation,alongwiththedifficultiesofcomputinghumanreadablefilesthat havenotbeentranslatedtosimpletabularformats.InthelineofGeckoToJoinerand BlastToJoiner,GeneParserproducesatabulartextfilewithasmanyrowsasannotations containedinthefileretrievedfromGenBank.Thefollowingtableshowsanexampleofwhat fieldsarekeptandtheirmeaning. Gene start Geneend Strand (reserved) Locustag Product 687 3158 f MHO_RS00005 membraneprotein insertaseYidC Table6:Extractedfieldsandparsedformatforanannotationfile.Thegenestartandgeneend fieldrefertothestartingandendingcoordinatesofthegeneinthegenome.Thestrandshowsan‘f’for forward,and‘r’forreverse.Areservedfieldisleftunusedforfutureimplementations.Thelocustag fieldshowsthetaggiventothegene,whereastheproductfieldreferstothedescriptionofthegene anditsfunctionalprofile. Table 6 shows an example of parsed annotation file containing one gene. More fields could be added in case of desire for particular and more advanced experiments, such as functional annotation. In particular, annotation files are used to obtain coding and noncoding readsabundance and toperformotherexperimentssuchasdifferentialexpressionorannotatedsearches. 4.2.2.3.CodingAbundances Thepresenceofspeciesinasamplecanbemeasuredusingthereadsabundance. However,thismeasuresmightnotbeenoughtoensurethecompositionofametagenome,and thereforemoreexperimentsarerequiredtobecarriedout.ThetoolCodingAbundances producesalistofabundancesofannotatedandunannotatedreads,alongwithafilepereach genomethatincludestheid’softhereadsassignedthatfallunderannotatedregions. Nonetheless,theconceptofannotatedreadisabitfuzzy:ReadsaresmallDNAsequences, andasinglereadmighthavebeenpartofagene,mightbethecompletegeneoralarger sequencethatincludesoneormoregenes,andatlast,itmighthavebeenthesocalledjunk DNA.Toaddressthistypeofsituation,mappedreadsaredividedinthefollowingcategories: (1) Matchedreads:Thosethatarenotannotated,notevenpartially. (2) Semiannotatedreads:Thosethatarepartiallyannotatedononesideor another. (3) Fullyannotatedreads:Thosethatarecompletelyinsideageneregion. (4) Extracoveredreads:Thosethatincludemorethanasinglegeneregion,and maycontainmore. Theproceduretoidentifythetypeofreadsfollowsthealgorithmdepictedbelow: 49
1.Sequentiallyopentheparsedannotationfiles. 2.Iftheparsedfileofannotationscontainsanyannotation: a.Loadannotationsintoafixedlengthbufferinmemory. b.Sweep the binary mapping file and extract the reads that maptothegenomewhoseannotationfilebelongsto. c.Sequentially sweep the annotations contained in the fixedbufferforeveryread. i.If a hit is found, mark the read and write it to theoutput. 3.Iftherearemoreannotationfiles: a.Gotostep1. 4.Elsestop. Figure21:PseudocodefortheCodingAbundancesprogram.Theusedapproachrunsincubictime, whichcouldbeslightlyreducedtoquadratictime.RAMbasedapproachedwouldefficientlyreduce computationtime,howevertheywouldbeunsuccessfulforlargedatabaseswithmanyannotations. It might be noticed by the reader that the running complexity could be reduced. In fact, the binary mapping file is swept through for every genome in the database, and the buffer containing the annotation files is linearly traversed for every read mapped to the genome. Therefore, given genomes, reads and a maximum of annotations per genome, g n k the running complexity is around , which is cubic. However, the number of (g n )O*3*k genomes is usually very small compared to the number of reads and same for the number g n of annotations . Therefore in the average case and thus the algorithm is k g≪k≪n tractable. A further modification to improve running times would be performing a binary searchovertheannotationscontainedinthebuffer,thusreducingitto .(g n ogk)O*3*l Thedifferentialexpressionexperimentisconductedonaparticulargenomebyonly consideringthereadsmappedtoannotatedregionsofthegenomeandcomparingabundance betweendifferentsamples(SeeFigure22)inthesamewayasRNAseqtranscriptome expressionanalysisisperformed. 50
Figure22:DNAseqdifferentialexpressionplot.Eachpointrepresentsanannotatedregionfora particulargenome.Inthexaxisandyaxis,thepercentageofreadsthataremappedtoeachannotated regiondividedbythetotalmappedreads. 4.2.2.4.Annotationalmapping Inthesamelineasthetooldiscussedabove,theSpecificDistributionofMatchestool aimstoofferacomparisonbetweenthequalityofthemappedreads,thedistributionof fragmentsreportedforthemetagenome(thefirstisasubsetofthesecond)andthetypeof matchesthatsaidparticulargenomehas.Thereforethistoolrequires(1)thebinarymapping fileoftheexperiment,(2)thereporteddistributionoffragmentsbythesequencecomparison software(SeeFigure23)and(3)theannotationfileforthetargetgenome. Thetoolproducesamatrixofthetypeofmatchesofsize where isthek003 × 1 k numberofrowsthatrepresentlength,and100thenumberofcolumns,thatrepresenta percentageandnormallythepercentageofidentityisused.Forexample,amatrixofsize impliesthattherearenoreadsoflengthlargerthan30bp.Sinceeverythreerows0 009 × 1 representonebasepair,row30willbeabasicmatchoflength10bp,row31willbea semiannotatedmatchoflength10bpandrow32willbeafullyannotatedreadoflength10 bp.Thenextthreerowsfollowthesameschemebutfor11bp.Thereforethetypeofmatch dependingontherowisextractedbyperformingthemodulusoperationandcomparingthe resttoeither0,1or2. 51
Figure23showsanexampleofannotationalmapping.Themainkeyofthisplotisto showthetypeofthematchedreadsforagenome(iftheseareannotatedornot)andto comparetheirproperties(identitiesandlength)withthefulldistributionoffragments. Figure23:Annotationalmappingplot.TheXaxisisthepercentageofidentityofthematches.The Yaxisisthelengthofthematches.TheZaxisistheaccumulationoffragmentsinlog10scale.The specifictypeofmatchesareplottedinpurple,greenandyellowtoillustratethedifferenttypeof matches. Thefollowingremarksshouldbetakenintoaccounttounderstandtheplot: 1. Inthebackground,the3Dgreydistributionisthedistributionoffragmentsbylength andpercentageofidentity. 2. Bycolors,andplottedontothegreydistribution,thetypeofmatchthatwasmappedto thegenome.Thetypeofmatchcanbeeitherunannotated(green),partiallyannotated (yellow)andfullyannotated(purple). 52
4.2.2.5.GenomeProfileofAccumulatedReads Animportanttermingenomicsandalsointhemetagenomicsfieldisthecoverage. Wehavespokenaboutcoverageinreferenceofameasurethatisthepercentageofamatching fragmentdividedbythetotallengthofaread.However,thistermisknownasanindicatorof howmanyreadsorDNAsequencespopulateorcover agenome.Inourcase,itisofinterest tocheckhowmanyreadsmaptoeachpositionofagenome,astodefineameasureof continuityalongagenome.Thehigherconcentrationofreadsalongthewholegenome,the better. Thebinarymappingfileisusedtofindallreadsmappedtoaparticulargenome,and theregionofthegenomewherethesereadsaremappedtoarestoredandaddedup.Intheend, aprofileoftheamountofaccumulatedreadsalongthegenomeisobtainedandplottedusing anRscript.SeeFigure24.However,bacterialgenomesareusuallyoflength1Mbpandupto severalmillions,andthereforeitbecomesdifficulttoplotmillionsandmillionsofpoints alongastraightlineinanimage.Forthismatter,theRscriptusesasmoothingwindowof parametrizedsize(saya1000bp)thatcomputestheaveragenumberofmappedreadsper 1000bp,andthusagenomeof1Mbpwouldonlyrequire pixels.Biggerwindows 103 106= 103 canbeusedtofitlargergenomes.Hencethecomplexityis sinceitisneededto(n)O+k sweepthebinarymappingfileonceandthelengthofthegenometoperformthesmoothing windowaverage. Figure24:Thegenomeprofileofaccumulatedreadsshowshowmanynucleotidesaremappedtoeach partofthegenome.Theplotissmoothedusingaparametrizedslidingwindowofsize1000basepairs, whichmeanseachpositionistheaverageofnucleotidesofthesurrounding1000basepairs. 53
4.2.2.6.GenomeCutter,SpecificRegionsandSpecificReads Wehaveproposedthequestion“Howmanyreadsmappedtoagenomeareneededto acceptorrejectthepresenceofaspeciesinasample?”earlierintheAnalysisandDesign section.Currentstateoftheartmetagenomicanalysissoftwaredonotofferanyprocedureto addressthisquestion.Thedevelopedworkflowoffersasetoftoolscapableofextracting thosereadsthataremappedtoregionsthatonlybelongtoonegenome,i.e.supportingwith strongevidencethataspeciesispresenteveninsituationswhereitsreadsabundanceislow. Figure25:Exampleofspecificregionsdetectedinadatabasewith .Eachhorizontalline 4N= representsagenome.Theydonotneedtobenecessarilyofthesamelength,theyareshowedinthe pictureinordertofacilitateitsunderstanding.Eachrectangle(butthedottedones)representsa fragmentsharedbetweenthegenomeofthelinewhereitispresentandthereferencegenome.The dottedrectanglesrepresentthespecificregionsofthereferencegenome,whereclearlythereareno fragmentsinanyofthegenomes. IfwetakealookatFigure25wewillseethatthistypeofproblemneedsgenomevs. genomesequencecomparisons.Forthatmatter,thefirststepofthistoolisrunningGECKO toperformacomparisonbetweenthetargetgenomeandtherestofthedatabase,thusleading to comparisonsforadatabaseofsize .Theoutputofthecomparisonistakenand 1N− N fedtotheSpecificRegionstool,whichwillextracttheregionsthatarefreeofmaps.The outputisproducedasalistofcoordinatesthatrepresentspecificregionstothetargetgenome. Then,theSpecificReadstoolusesthelistofregionsandthebinarymappingfileofthe experimenttocomputewhichreadsaremappedtothespecificregions. GECKOisathirdpartysoftwareandthereforeitscomplexitywillnotbeanalyzed.In ordertorunthecomparisonwithGECKO,wewillneedfirsttoremovethetargetgenome fromthereferencedatabase,sinceifitisincludeditwillproduceafullalignmentand fragmentsateverypositioninthegenomedependingonthekmersize.Toavoidthis,the GenomeCutterprogramremovesaparticulartargetgenomefromthereferencedatabase.Its complexityisstraightlinearsinceitonlyrequirestocopyandwritesequencesfromonefile toanother,skippingthetargetedone. 54
TheSpecificRegionstoolloadsanarrayofbytesoflengthequaltothatofthe genome,whereeachbyterepresentseithermappedregionorfreeregion.Thiscouldbe improvedusing ofthememoryneededbyimplementingitatabitlevel,however,eventhe 4 1 humangenomecouldbefittedintoRAMusingthebyteprocedure(~3GB).Finally,a sweepthroughisperformedonthearray,filteringthoseregionsthataresmallerthanagiven length(e.g.30basepairs).Thustherunningtimeis being thenumberoffragments(2n)On producedbytheGECKOcomparison.TheSpecificReadstoolrunsonthelistofregionsand thebinarymappingfile.First,theregionsareloadedintoramandthenforeverymappedread tothegenomethearrayissweptthroughcheckingifareadispartiallyorcompletelyinside, similarlytothetoolSpecificDistributionofMatches.Thereforeit’srunningtimeis dependentonthesizeofthebinarymappingfileandthenumberofregions.However,the numberofspecificregionsisinverselyproportionaltothelengthofthedatabase,andthe numberofmappedreadsisdirectlyproportionaltothesizeofthedatabase(amongothers), thereforeitthenumberofspecificregionsishigh,thenumberofmappedreadsissmaller, makingitinaverageoflessercomplexitythanquadratic. 4.2.2.7.CoverageIdentityMatrix Toprovideanothermeasureofqualityrelativelytothewholemappingdistribution, forexample,toshowpropertiesofallmappedreadsatonce,theCoverageidentitymatrix programwasdeveloped.Itproducesa tabulartextfilematrixwhere isthen×m×k n percentageofidentity(from0to100), isthepercentageofcoverage(from0to100asm well)and istheaccumulationofmappedreads.SuchmatrixisreadbyanRscriptplotsak 3Ddistributionofthematrix(seeFigure26)applyinglogarithmicscaletothe matrixsincek theaccumulationofmappedreadscanbeextremelyhighindispersedpointsandratherlow onothers. 55
Figure26:3DPlotofthedistributionmatrixproducedbyCoverageIdentityMatrix.Onthexaxis,the percentageofidentityofthemappedreads.Ontheyaxis,thepercentageofcoverageofthemapped reads.Onthezaxis,theaccumulationofmappedreadsforagivenpointofpercentageofidentityand coverageinlogarithmicscale.Noticethatthedistributionissegmentedat60%identityand30% coverage. Figure26showsthatthemostaccumulationoccursatthepoint100%identityand 100%coverage,whichimpliesthemappingmodulemappedmostreadswithgoodproperties. However,itcanalsobeseenthataconsiderableamountofmappedreadshadverylow coverage,henceonlyasmallpartofthereadwasmapped.Moreaccumulationtowardsthe backcornermeansabettermapping. ThecomplexityofCoverageIdentityMatrixislinearonthesizeofthebinary mappingfile,foritsweepssuchfileandcountsthenumberofmappedreadsina twodimensionalvector.Onceitfinishessweeping,thevectoriswrittentodiskasa matrix,whichisconstantandthereforenotaccountedforthecomplexity.00 001 × 1 56
4.2.3.Galaxyimplementation As mentioned in other sections of this manuscript, the whole workflow is available under a Galaxy workflow manager implementation, which makes it simple to run experiments without having to face the difficulties that arise from shellcalls that have many parameters and that have to be run in particular order. Installing a Galaxy instance in a local machine might be tedious in some cases and, furthermore, if we only wish to try out we would need to remove all installed plugins after checking the software out. Therefore we have provided with an Ubuntu virtual machine with all binaries installed along with the Galaxy instance to avoid bothering the final user with installing tutorials and derived issues. It is recommend that the virtual machine is run at least with 4GB of RAM, since metagenomics processing requires stillaconsiderableamountofmemory. In addition, a Galaxy Guided Exercise is included in the documentation which shows how to run the virtual machine, start the Galaxy Instance and also how to run different types ofexercises. 4.2.4.Galaxy’stooldefinitionlanguage GalaxyusestheXML[46]languagefortooldefinition.Thisis,eachdevelopedbinary thatcomposetheworkflowmusthaveanXMLfilewhichuniquelydeterminestheinputsof thebinaryexecutable(alongwithinformationsuchastheformatoftheinputs,thelabelsand helpsthatwillbedisplayed,etc.),theoutputs(similarlytotheinputs)andthesystemcallto thebinaryexecutable.SeedetailsattheGalaxyToolConfig tutorial. 15 TheXMLfilesmustbefirstcreatedforeachtoolandthenaddedtothetoolregistry (namelytothetool_conf.xml.samplefile)underasectionthatincludesthenameofthe workflowandthepathstoall.xmltools. <toolid="mgReadsIndex2"name="MGRindex"> <description>Createsindexesforthefastaccessoffasta files</description> <inputs> <paramname="input_fasta_file"type="data"format="fasta" label="Fastafile"help="Createanindexofafastafile"/> </inputs> <command>/home/galaxy/galaxy/tools/metagecko/binaries/mgReadsIndex2 $input_fasta_file$input_fasta_file;mv${input_fasta_file}.nonsorted $nosorted_index_output0;mv${input_fasta_file}.sorted 15https://wiki.galaxyproject.org/Admin/Tools/ToolConfigSyntax 57
σ H0: MEGAN = σMG σ = σ H1: MEGAN / MG Whichyieldedapvalueof0.69,whichimpliesthesuppositionthatstandard deviationsareunknownbutyetequalatanyregular .α Contrastedhypothesisonthedifferenceofmeans: μ H0: MEGAN − μMG = 0 μ = H1: MEGAN − μMG / 0 The resulting pvalues was of 0.99 and therefore the null hypothesis was accepted at any , concluding that the difference of means of the paired readabundance values was α nearlyzero,andthusthereexistnosignificantdifferencebetweentheresults. 5.2.Sickandhealthyswinesamples 5.1.1.Datasetavailability For the second comparison, we used two collections of metagenome samples from the respiratory tract of healthy and diseased Sus scrofa taken from the swine slaughterhouse. The samples had been sequenced using the 454pyrosequencing technology. The diseased samples contains 678,500 reads whereas the healthy sample contains 581,383 reads. The analysis is based on the understanding that the pigs are naturally infected, and therefore the point of observing significant variation in particular bacteria between the diseased and healthy sample suggest that such species are related to the condition of the swines. Unfortunately, the dataset used in this comparison has not been made public at the time of writing the manuscript and remainsyetprivate. The raw samples’ read average length was 1,200 bp before trimming and filtering. Filtering was performed using Replicates and trimming was performed using LUCY [48]. The average length of the reads dropped from 1,200 bp to 500 bp in average, which is expected for the type of reads after removing adapters, lowcomplexity regions, etc. The reference database was built using 50 custom genomes which are known to be related to the respiratory tract of swines. The included bacteria was comprised of mostly actinomycetales, bacteroidetes, firmicutes, proteobacteria and mycoplasmataceae, and accounted for a total of 283 DNA sequencesduetoscaffolds,contigsanddifferentstrains. 5.2.2.ComparisonwithMEGAN 64
Asstatedearlier,inordertoprovethattheresultsoftheproposedworkfloware consistentwiththoseofothermetagenomicanalysissoftwaresuites(intermsofabundancein thetaxonomicclassification),thefollowingtestwasperformedusingresultsfromBLASTN basedonmetagenomicsamplesfromtherespiratorytractofswines.Both,theworkflow(MG workflow)andMEGANwereexecutedusingthesameinputfromBLASTNandranwith defaultparameters. Oncomparisonofthesickandhealthymetagenomes(inadvance,M01andM02, respectively)basedonMEGAN,theabundanceplot(SeeFigure28)showssimilarresultsand tendencetoours.Readsabundanceweregroupedbystrainstoshowconcordance.However, whereasthemappingprocedureofametagenomeusingMEGANcanlastnearlyhalfanhour, theproposedworkflowtook1minuteand22.8secondstoanalyzethediseasedmetagenome and1minute15.178secondsforthehealthyonewhenBLASTthecomparisonhadbeendone withBLAST.RuntimeexecutionsweremeasuredusingastandardInteli5machinewith4GB ofRAM. Althoughsomesmalldissimilitudeappeartoexistbetweenthereadsabundance results,itisimportanttonoticethattheplotisshowingreadsabundanceinlogarithmicscale, andthereforeamplifyingthevariationsofsmallandlowabundantgenomes.Thisalso suggeststhat,asexpected,thedevelopedworkflowcanreportmoreevidenceintermsof lowabundantgenomesinpartduetothefragmentsextensionstepandtheadequate combinationoffilters. 65
Figure28:ReadsabundancecomparisonbetweenMEGANandtheproposedmethodperstrains.Top plotshowsthereadsabundanceforthediseasedmetagenome(M01)inlogarithmicscale.Bottomplot showsthereadsabundanceforthehealthymetagenome(M02)inlogarithmscale. 5.2.3.Statisticalvalidation Acontrastonthedifferencesofmeanswasperformedassumingbothsamples belongedtoanormaldistribution,inordertodeterminewhetherornotthepropositionthat 66
bothresultscamefromthesamedistributioncouldbeacceptedornot.Forthispurpose,the readsabundancepergenomewereconsideredaspaireddatapointsinthecontrasted hypothesis. ● Apriorihypothesisonthecontrastoftheequalityofstandarddeviationsforthe metagenomeM01: σ H0: MEGAN = σMG σ = σ H1: MEGAN / MG Whichyieldedapvalueof0.51inthediseasedmetagenome,whichimpliesthe suppositionthatstandarddeviationsareunknownbutyetequalatanyregular .α ● Apriorihypothesisonthecontrastoftheequalityofstandarddeviationsforthe metagenomeM02: σ H0: MEGAN = σMG σ = σ H1: MEGAN / MG Apvalueof0.58wasobtainedinthehealthymetagenome,whichimpliesthe suppositionthatstandarddeviationsareunknownbutyetequalatanyregular .α ● ContrastedhypothesisonthedifferenceofmeansforthemetagenomeM01: μ H0: MEGAN − μMG = 0 μ = H1: MEGAN − μMG / 0 The resulting pvalue of the M01 metagenome was of 0.29 and therefore the null hypothesis was accepted at any standard , concluding that the difference of means of the α paired readabundance values was nearly zero, and thus there exist no significant difference betweenthetaxonomicresults. ● ContrastedhypothesisonthedifferenceofmeansforthemetagenomeM02: μ H0: MEGAN − μMG = 0 μ = H1: MEGAN − μMG / 0 The resulting pvalue of the M02 metagenome was of 0.23 and therefore the null hypothesis was accepted at any standard , concluding that the difference of means of the α paired readabundance values was nearly zero, and thus there exist no significant difference betweenthetaxonomicresults. 67
68
CHAPTER6. CONCLUSIONS 6.1.Conclusionsofthedevelopedworkflow The proposed workflow is an example of simple lowcoupled modules that when put together compose a complete processing method capable of producing fast and reliable results in open environment with easytohandle formats. In the computational sense, our intention was to provide an alternative platform for metagenomics analysis, which, in essence, could be easily adapted and expanded to satisfy the growing demand of experiments that researchers are proposing in this novel field, while maintaining a lowlevel programming strategy that keptcomputationspeedasoneofthepriorities. Furthermore, metagenomics is an effervescent field and there are still a number of questions that need to be addressed before a stable version of a definitive data analysis software becomes available. Currently, metagenomic analysis tools generally represent a closed environment and offer few configuration options and limited extension possibilities. In this sense, our aim was to develop a software framework to which other modules could be added. An additional motivation to develop this software was the need for software sensitive enough to detect the presence of lowabundance species. Our intent was to provide data in standard and editable formats that facilitate further analysis with external software. The proposed workflow software offers several notable advantages over the software currently available in the market. Firstly, the use of GECKO enables this software to compute similarity searches in the samples against a collection of genomes in a reasonable time. Providing different mapping alternatives helps set up a sort of quality measures of the mapping process basedonabundancedifferencesacrossmappingalternatives. The proposed software is designed to provide evidence of the presence of lowabundance species by finding particular specific regions of genomes with mapped reads. These mapped reads provide strong evidence of the species present in samples. The methods developed for assessing and evaluating the quality of mapping also improve accuracy and reliability in terms of the identification of the species present in a sample. However, from our perspective, the most important contribution of this workflow software is that it offers the possibility of incorporating new modules to extend the analysis workflow by showing datafile specifications, which enables finegrained metagenomic data analysis via a simple Galaxy interfacewhichrunsexperimentsatthecostofafewclicks. 69
Inaddition,theproposedworkflowdepictedonthismanuscript,hasbeensubmittedto thehighlyrankedjournalBMCGenomics asaresearcharticle,andisundersecondround 17 inspectionatthetimeofwritingthisconclusions. 6.1.1.Conclusionesdeltrabajodesarrollado Elmétododesarrolladoesunejemplodesimplesmóduloscuyopocoacoplamiento permitelacomposiciónenunsistemacompletodeprocesamientocapazdeproducir resultadosrápidayfiablementeenunentornoabiertoyconformatosdedatossencillosde manejar.Entérminoscomputacionales,nuestraintenciónfueladeproveerunaplataforma alternativaparaelanálisisdemetagenomas,lacual,enesencia,pudieraserfácilmente adaptadayexpandidaparasatisfacerlagrandemandadeexperimentosquelosinvestigadores sugierenenestecampotannovedoso,alaparquemanteniendounenfoquedeprogramación abajonivelcuyaprioridadprincipaleselreducidocostecomputacional. Elcampodelametagenómicaestáenalzayaúnsiguenexistiendograncantidadde preguntasquedebeserrespondidasantesdequeunaversiónestándardepaquetedeanálisis demetagenómicaestédisponible.Actualmente,lasherramientasdeprocesamientode metagenomasestángeneralmentedesarrolladasbajounentornocerrado,ofreciendoporello pocasposibilidadesdeconfiguraciónyflexibilidad,alavezqueescasaampliación.Eneste sentido,hemospretendidodesarrollarunmarcodetrabajoalqueselepudieranañadirnuevos módulossindificultad.Ademásdeello,otramotivaciónhasidolanecesidadactualde softwaremásfino,capazdedetectarytratarespeciescuyaabundanciafuerabaja.Nuestra intenciónfue,además,proveerresultadosdemaneraestándaryabiertaparafacilitarun análisisposteriorconherramientasexternas.Elmétodopropuestoofrececiertasventajas sobreelsoftwareactualmenteexistente.Primariamente,elusodelpaqueteGECKOpermite computarcomparacióndesecuenciasentremuestrasmetagenómicasygenomaspatrónen tiemporazonable.Laflexibilidadenlasalternativasdelafasedemapping nospermite extrapolarmedidasdecalidaddelprocesobasándonosenlasdiferenciasdeabundancias. Elpaquetedesarrolladofuediseñadoconobjetodeproveerevidenciassobrela presenciadeespeciesconpocaabundanciabuscandozonasparticularesdegenomascon lecturasasignadas,pueséstasdansoportemásfuertesobrelapresenciadedichogenomas. Losmétodospropuestosparaevaluarlacalidaddelaasignacióntaxonómicapromueven mayorfiabilidadenlaidentificacióndeespecies.Sinembargo,desdenuestropuntodevista, lacontribuciónmásimportanteeslaposibilidaddeincorporarnuevosmódulosparaextender elflujodetrabajo,dadaslasespecificacionesdelosformatosdesalida,permitiendoasuvez unaexperienciaricadeprocesamientodedatosconsólounospocosclicksatravésdelgestor deflujosdetrabajoGalaxy. 17https://bmcgenomics.biomedcentral.com/ 70
Porúltimo,elflujodetrabajodesarrolladodescritoenestedocumentohasido presentadoalarevistadeimpactointernacionalBMCGenomicscomounartículode investigación,yseencuentraactualmente,alafechadeescrituradeestemanuscrito,bajola segundarondaderevisión. 6.2.Furtherextensionsandimprovements AlthoughwetriedtosqueezethelimitedtimethatisdedicatedtotheEndofDegree Projecttofitallthemodulesandresultswewantedintheproposedworkflow,itwas unavoidabletoleavesomethingsoutofthepanorama.Nevertheless,itisourintention,for boththestudentandthetutor,tokeepworkingonthisplatform,developingprocessingtools thatanswerkeyquestionsinmetagenomics,facilitatingthestilldifficultprocessingstepsofa wholeexperiment,reducingcomputationalrequirementsinbothtimeandspace,providing withnoveltytoolstoassertnewprocessingstandards,etc. Morespecifically,thereareafewthingsthatwedefinitelywishedtoimplement,but sadlyhadnottimeforit.Thefollowingproposalsfallundersuchcategory: ● Aphylumtaxonomictreeviewthatenablesuserstoaddupreadsabundanceon differentspecieslevels. ● Acompletesetofdifferentmappingalternatives,suchasincludingtheLCA algorithm. ● Apreselectionsystemtosignificantlyreducethesizeofthereferencedatabaseto speedupcomputationtime. ● Amorecompletedfunctionalprofileofthespeciescontainedinthesamples,e.g. settingupmetabolicpathways. 6.3.Acknowledgements Iwishtoexpressmysinceregratitudetomyadvisor,Prof.OswaldoTrelles,forthe invaluableopportunitythathegavemebyofferingmeaplaceinhisresearchteamandinthe Bioinformaticsworld,andforhavingfaithinmyabilitiesduringthedevelopmentofthis work.IalsowishtothankallofmycolleaguesattheBitlabteam,andspeciallytoÓscarfor allthesupportandprovidedguidance. Lastbutnottheleast,Iwishtothankmyfamilyforallwhattheyhavedoneforme throughtheyears,andtoAnnette,forwithoutherdaytodaysupportnoneofthiswouldhave beenpossible. 71
72
CHAPTER7. BIBLIOGRAPHY 1. EnisAfgan,DannonBaker,MariusvandenBeek,DanielBlankenberg,DaveBouvier, MartinČech,JohnChilton,DaveClements,NateCoraor,CarlEberhard,Björn Grüning,AysamGuerler,JenniferHillmanJackson,GregVonKuster,EricRasche, NicolaSoranzo,NiteshTuraga,JamesTaylor,AntonNekrutenko,andJeremy Goecks.TheGalaxyplatformforaccessible,reproducibleandcollaborative biomedicalanalyses:2016update.NucleicAcidsResearch (2016)doi: 10.1093/nar/gkw343. 2. Wilde,Michael,etal."Swift:Alanguagefordistributedparallelscripting."Parallel Computing 37.9(2011):633652. 3. Wolstencroft,Katherine,etal."TheTavernaworkflowsuite:designingandexecuting workflowsofWebServicesonthedesktop,weborinthecloud."Nucleicacids research (2013):gkt328. 4. Mardis,ElaineR."Theimpactofnextgenerationsequencingtechnologyongenetics." Trendsingenetics 24.3(2008):133141. 5. DeRoure,D.,Goble,C.andStevens,R.(2009):TheDesignandRealisationofthe myExperimentVirtualResearchEnvironmentforSocialSharingofWorkflows. FutureGenerationComputerSystems 25,pp.561567. 6. Benson,DennisA.,etal."GenBank."Nucleicacidsresearch 41.D1(2013):D36D42. 7. Kanehisa,Minoru,andSusumuGoto."KEGG:kyotoencyclopediaofgenesand genomes."Nucleicacidsresearch 28.1(2000):2730. 8. Mende,DanielR.;AlisonS.Waller;ShinichiSunagawa;AinoI.Järvelin;Michelle M.Chan;ManimozhiyanArumugam;JeroenRaes;PeerBork(20120223). "AssessmentofMetagenomicAssemblyUsingSimulatedNextGeneration SequencingData".PLoSONE . 9. RichterichP.(1998):Estimationoferrorsin"raw"DNAsequences:avalidation study.GenomeRes .8(3):251–259. 10. Huson,DanielH.,andNicoWeber."MicrobialcommunityanalysisusingMEGAN." Methodsinenzymology 531(2012):465485. 11. Altschul,S.F.,Gish,W.,Miller,W.,Myers,E.W.&Lipman,D.J.(1990)"Basiclocal alignmentsearchtool."J.Mol.Biol .215:403410. 12. Gish,W.&States,D.J.(1993)"Identificationofproteincodingregionsbydatabase similaritysearch."NatureGenet .3:266272. 13. MorgulisA.,CoulourisG.,RaytselisY.,MaddenT.L.,AgarwalaR.,&SchäfferA.A. (2008)"DatabaseindexingforproductionMegaBLASTsearches."Bioinformatics 15:17571764. 73