Genòmica Computacional

1 GENÒMICA COMPUTACIONAL Autor: Sergi Alcaide i Portet Ponent: David Carrera Pérez Codirectors: David Torrents Arenales Friman Sánchez Castaño Treball de Fi de Grau 12 d’octubre de 2015 2 Contingut 1 Introducció ..................................................................................................................................................... 4 1.1 Part biològica ...................................................................................................................................... 4 1.2 Motivació del projecte ..................................................................................................................... 6 1.3 Objectius ................................................................................................................................................ 7 1.4 Estat de l’art ......................................................................................................................................... 8 1.4.1 Antecedents en mutacions somàtiques .......................................................................... 8 1.4.2 Estat actual d’estudis poblacionals ................................................................................... 8 1.4.3 Eines utilitzades: MPI, Suffix-Tree, C++ ........................................................................... 9 1.5 Competències .................................................................................................................................... 10 2 Planificació ................................................................................................................................................... 12 2.1 Descripció de les tasques ............................................................................................................. 12 2.2 Explicitació de l’ordre lògic de les tasques ........................................................................... 16 2.3 Diagrama de Gantt .......................................................................................................................... 17 2.4 Anàlisi de les modificacions en la planificació inicial ....................................................... 19 2.5 Valoració d’alternatives i pla d’acció ....................................................................................... 19 3 Pressupost .................................................................................................................................................... 21 3.1 Identificació i estimació dels costos ........................................................................................ 21 4 Cos del treball ............................................................................................................................................. 26 4.1 Flux de dades de l’aplicació ......................................................................................................... 26 4.1.1 Suffix Tree .................................................................................................................................. 26 4.1.2 Flux de dades ........................................................................................................................... 28 4.2 Creant la 1a versió........................................................................................................................... 28 4.2.1 Algoritme general de l’aplicació ...................................................................................... 28 4.2.2 Estructura de l’aplicació ...................................................................................................... 29 4.3 Anàlisi ................................................................................................................................................... 40 4.4 Millores ................................................................................................................................................ 48 4.4.1 Master amb Suffix Tree (Master Tree) ............................................................................ 48 4.4.2 Partitions ................................................................................................................................... 49 4.4.3 Reverse ....................................................................................................................................... 52 4.4.4 Lectura d’arxius comprimits ............................................................................................. 54 4.4.5 Comunicació Master-Slaves ................................................................................................ 54 4.5 Resultats .............................................................................................................................................. 55 4.6 Estudis poblacionals ....................................................................................................................... 57 5 Conclusions .................................................................................................................................................. 59 5.1 Conclusions generals...................................................................................................................... 59 3 5.2 Ètica del projecte i regulacions adherides ............................................................................ 60 5.3 Sostenibilitat i compromís social .............................................................................................. 60 5.3.1 Dimensió econòmica ............................................................................................................ 60 5.3.2 Dimensió social ....................................................................................................................... 60 5.3.3 Dimensió ambiental .............................................................................................................. 61 6 Referències .................................................................................................................................................. 62 4 1 Introducció 1.1 Part biològica Avui en dia una de les preocupacions en quant a salut que més incideixen en la població és la del càncer, això ha provocat l’augment de recursos destinats a la investigació en aquest camp. Quins canvis es produeixen en el genoma d’un individu que provoquin càncer, com es poden evitar aquests canvis... Moltes preguntes han sorgit en aquest camp i són molts els científics que intenten donar-hi resposta. Una de les estratègies que s’està emprant més en els estudis actuals de la genòmica del càncer es basa en el principi de comparar genomes de cèl·lules sanes amb cèl·lules tumorals del mateix pacient per tal d’identificar els canvis (mutacions somàtiques) que s’han produït en el genoma tumoral que puguin ser responsables de la malaltia. Per això, es prenen dues mostres del pacient, una d’una zona sana i l’altra d’una zona afectada pel càncer, i se’n extreu l’ADN (àcid desoxiribonucleic) per procedir a la seva lectura, és a dir, la seva seqüenciació 1 . El genoma humà està compost per quatre nucleòtids: adenina, guanina, citosina i timina, simbolitzats per les lletres A, G, C, i T, respectivament. El procés de seqüenciació converteix les molècules d’ADN en milers de milions de lectures (text) que es recullen en fitxers de format fast queue (.fq) on, a més de les subseqüències extretes també tenim informació que ens indica la qualitat/fiabilitat de cada mostra. Per tal de neutralitzar els errors de seqüenciació que puguem induir del procés (que no és perfecte), el procés de lectura cobreix tot el genoma moltes vegades (unes 30 o 60, també anomenat coverage 2 ), amb la conseqüència directa que el nombre de dades també es multiplica amb la mateixa constant. Tenint en compte que el genoma humà consta de 3,2·109 lletres (nucleòtids) aproximadament i que s’esperen no més de 10.000 mutacions en cada tumor, la identificació d’aquests canvis suposa un repte computacional important a nivell algorísmic i d’implementació. A més cal considerar que les subseqüències que llegim del genoma no estan ordenades i per tant tenim milers de milions de subseqüències desordenades i que estan replicades vàries vegades (degut al coverage). A més a més, si considerem que les millores en les tècniques de seqüenciació estan facilitant aquest procés i generant milers de genomes seqüenciats, els protocols d’anàlisis de seqüència també han de considerar, més enllà de la detecció fiable de mutacions, que aquest procés sigui escalable i que arribi a un nivell de paral·lelització que permeti processar centenars de mostres a la vegada. D’aquesta manera un problema que era inicialment de caràcter biològic passa a ser un problema de caràcter computacional. Per trobar les mutacions somàtiques es requereixen dos genomes per tant 6.400 milions de bases que acostumen a tenir coverage d’entre 30 i 60 (384 GBytes 3 ). Per poder fer l’anàlisi a més de la informació essencial (subseqüències) es requereix també informació addicional que 1 Les tècniques de seqüenciació llegeixen la cadena d’ADN de la mostra i mitjançant la reacció que ofereix cada punt de la mostra a un estímul lumínic es classifica en els diversos tipus de nucleòtids (adenina, citocina, guanina o timina). A més a més assigna a cada nucleòtid un caràcter ascii que indica la fiabilitat de la classificació. 2 Mot que s’empra per indicar el nombre de vegades que ha estat seqüenciat cada genoma. El procés de seqüenciació no és un procés exacte i produeixen errors durant l’extracció de les bases. Per aquesta raó normalment els genomes se seqüencien múltiples vegades per poder eliminar/filtrar aquests errors. 3 Comptant que cada caràcter s’emmagatzema en 8 bits (1 Byte) i cada genoma té 3.200 milions de bases i un coverage de 60. 5 relaciona, comptabilitza i guarda informació associada (com són un nom que identifica la subseqüència, nosaltres no la utilitzarem, i la qualitat de la mostra) . Com es pot comprovar el problema que es pretén tractar és un problema complex especialment pel que fa el volum de dades. Malgrat s’han invertit molts esforços en trobar solucions a aquest problema i s’han generat nombrosos programes per la detecció de mutacions somàtiques, aquest procés encara no esta resolt i segueix suposant un repte important per la biomedicina. De fet, amb la facilitat de seqüenciació, l’anàlisi d’aquestes seqüencies s’ha convertit en el coll d’ampolla dels estudis de les malalties a nivell genòmic. Pràcticament tots els programes existents es basen en la comparació de les seqüencies del pacient contra un genoma de referència. Aquest pas, malgrat facilita el procés, també comporta unes limitacions inherents al procés d’alineament, fent que no s’assoleixin nivells de fiabilitat suficients per la majoria d’estudis, o per l’aplicació d’aquesta metodologia a la clínica 4 . Altres alternatives, com el programa SMUFIN (Somatic Mutation Finder), es basen en un altre tipus de processament primari de les dades, el que resulta en una millor detecció de les mutacions. Degut a que aquests programes s’han realitzat en un entorn estrictament biològic, les implementacions existents encara deixen espai per moltes millores que requereixen solucions més computacionals, també en quant a trobar la millor combinació software-hardware. Aquest projecte es basa en el disseny i implementació de solucions alternatives per la identificació de mutacions somàtiques a partir de la comparació de dos genomes, un sa i un altre malalt. Per això hem escollit, com a base, parts de l’estratègia utilitzada pel programa SMUFIN, desenvolupat al grup de genòmica computacional del BSC (Barcelona Supercomputing Center). Malgrat aquest programa ja ha superat a altres existents en rapidesa i en fiabilitat, encara queden molts punts per millorar-lo. A més, gran part del codi d’SMUFIN no és interpretable ni avaluable, probablement degut a les pressions de temps entre les que es va generar, i té moltes dependències externes pel que fa a llibreries. Un altre camp que s’està potenciant molt últimament és el d’estudis poblacionals, on es busquen quines similituds i diferències tenim dins d’una mateixa població. Una gran mostra de la importància que està tenint aquest camp és el gran nombre de projectes de seqüenciació que existeixen actualment com UK10K [1], GoNL [2] o el famós projecte dels Estats Units “Precision Medicine” amb l’objectiu de seqüenciar un milió de voluntaris (amb un pressupost inicial de 215 milions de dòlars). Nosaltres també som conscients de la importància d’aquest camp i creiem que en el nostre cas també podem aplicar un algoritme semblant al de SMUFIN (en quan a l’ús del Suffix Tree), que pugui analitzar aquestes dades. És per això que també ens plantegem l’objectiu de que la nostra aplicació també serveixi per a fer estudis poblacionals, el que representa un augment del nombre de genomes que haurem d’introduir com a entrada i per tant un augment molt important pel que fa el nombre de dades que utilitzarà el programa. A més el canvi que es fa respecte les mutacions somàtiques és que haurem de contraposar tots els genomes contra els altres i que hi haurà moltes més divergències que en el cas d’un mateix individu, degut a la dimensió d’aquest altre objectiu i la limitació de temps ens centrarem en aconseguir fer una versió senzilla que validi l’ús del Suffix Tree en estudis poblacionals. 4 Sobretot degut a que al contraposar-lo amb el genoma de referència, tenim moltes altres variacions i per tant es complica encara més el trobar les diferències entre sa i tumoral. 6 1.2 Motivació del projecte L’SMUFIN és una aplicació desenvolupada al BSC (Barcelona Supercomputing Center), concretament al grup de Computational Genomics, que es basa en el principi de comparar els genomes sans i tumorals d’un mateix pacient per tal de trobar les mutacions somàtiques que puguin provocar el càncer i no utilitzar el genoma de referència per alinear-los. Per a cada mutació trobada n’indica de quin tipus és i la posició dins el genoma sencer. En el moment de la seva publicació [3], va demostrar ser un dels pioners en aquesta estratègia amb les bones estadístiques que s’indiquen en l’article, tant en especificitat, sensitivitat i rapidesa, especialment en les mutacions de medul·loblastomes en nens, leucèmia mieloide aguda, càncer de pròstata i limfomes. El flux de dades general de SMUFIN és el següent: Il·lustració 1: Flux de dades simplificat del programa SMUFIN La part estrictament biològica com hem comentat consisteix en agafar cèl·lules sanes i tumorals d’un pacient i seqüenciar-les, amb això obtenim fitxers en format fastq dels dos tipus de cèl·lules: El format fastq és el següent: Per a cada read ens trobem quatre línies consecutives: 1. La primera és un nom que indica el read per exemple: HW6217 2. La segona és el read pròpiament dit, i que conté els nucleòtids (identificats per les quatre lletres A, C ,G i T): ACGTAGTCAGCTA 7 3. La tercera és un signe de suma: + 4. Per acabar tenim una línia amb el mateix nombre de caràcters que la segona que ens indiquen la qualitat de cada nucleòtid segons una escala 5 . A partir d’aquí comença la part computacional, per començar es fa un filtratge de les dades per tal de reduir el volum de dades, filtrant així les dades de menor qualitat que no poden ser considerades com a correctes. Tot seguit es comparen les seqüències de cada genoma i s’alineen les parts conflictives que poden desencadenar una possible mutació, a continuació s’identifiquen el tipus i posició de cada mutació i es presenten en un format estàndard que en permet la validació. Un cop fet l’alineament també es poden guardar les dades en aquell punt en format BAM 6 . A més també es permet començar l’aplicació amb arxius d’aquest format, permeten d’aquesta manera continuar amb el flux de dades utilitzant dades que provinguin d’una execució anterior o inclús d’un altre programa. 1.3 Objectius L’objectiu principal d’aquest projecte és presentar una implementació vàlida utilitzant de rerefons la idea de SMUFIN de contraposar les seqüències sanes i tumorals d’un pacient de la primera part del flux de dades de SMUFIN en concret, llegir els arxius fastq de forma efectiva, i amb la informació d’aquests arxius crear una estructura de tipus Suffix Tree 7 , que ens permeti buscar-hi les mutacions. Aquesta estructura s’haurà de crear entre diferents nodes i per tant haurà d’estar correctament paral·lelitzada per tal d’aprofitar els màxims recursos dels que es disposi (sobretot pel què fa a memòria). A més intentarem en la mesura del possible reduir el nombre de llibreries utilitzades per no tenir un gran nombre de dependències externes. També s’intentarà utilitzar aquesta estructura per donar una alternativa vàlida per el problema dels estudis poblacionals. En aquest apartat buscarem demostrar que l’alternativa que presentem es pot aplicar en aquest camp. Nosaltres utilitzarem un conjunt reduït de dades i en validarem el resultat i demostrarem que és escalable a nivell poblacional. En ambdós casos el principal problema serà el gran ús de memòria i és per això que en els casos que haguem de prendre decisions entre reduir la memòria o reduir el temps d’execució, escollirem la primera. 5 Per a més informació consultar: https://en.wikipedia.org/wiki/FASTQ_format 6 Els fitxers .bam són fitxers binaris que normalment estan lligats a fitxers en format .sam que permeten mitjançant eines com SAMtools veure gràficament l’alineament fet i indexat amb tots els reads d’entrada. 7 Més endavant veurem com és aquesta estructura. 8 1.4 Estat de l’art 1.4.1 Antecedents en mutacions somàtiques L’augment que hi hagut darrerament pel que fa a la investigació i detecció de mutacions somàtiques ha provocat un creixement pel que fa a nombre de projectes relacionats amb aquest camp en tots els àmbits tant a nivell estatal, europeu i com a nivell mundial, provocant l’aparició de diferents aplicacions que són capaces d’analitzar genomes i detectar-ne mutacions somàtiques. La diferència que tenen aquestes aplicacions amb SMUFIN és que totes utilitzen sempre alguna eina per indexar els genomes amb l’anomenat Genoma de Referència, fet que els hi facilita la feina en el moment d’alinear els reads però que provoca l’augment de falsos positius degut a les diferències entre el genoma del pacient i el de referència. Com es veu en l’article de SMUFIN[3], totes les aplicacions que hi estan mencionades, tenen uns resultats pel que fa a especificitat i sensibilitat menors a aquest, degut entre d’altres al problema que genera la indexació amb el genoma de referència. Veiem d’aquesta manera que SMUFIN actualment té una força avantatge respecte els seus competidors i per tant consolida la idea d’utilitzar l’estructura de Suffix Tree, ja que no tenim constància que cap dels seus competidors l’utilitzi. Programes detectors de mutacions somàtiques semblants a SMUFIN:  Pindel [4]  BreakDancer [5]  GenomeSTRIP [6]  Genome Analysis Toolkit (GATK) [7]  CREST [8]  Delly[9]  MuTect [10]  ... Tot i tenir estratègies o pipelines diferents, en un punt del programa tots tenen en comú que acaben utilitzant alguna eina com són el MAQ (Mapping and Assembly Qualities) o el SSAHA2 (Sequence Search and Alignment by Hashing Algorithm) que els permeten alinear amb el genoma de referència. Un dels projectes que té més similituds és un que s’està duent a terme en paral·lel, també en el BSC, pel grup d’investigació de David Carrera, en què es pretén crear una estructura amb els reads normals i tumorals. En aquest cas però utilitzant en lloc de totes les subseqüències els anomenats k-mers 8 . Aquesta estratègia és més ràpida i no té un consum tant alt de memòria però implica guardar menys informació en l’estructura final. En el fons es tracta de la mateixa idea que el Suffix Tree però reduint-te la quantitat d’informació que aquest emmagatzema i sense connectar els diferents nodes. 1.4.2 Estat actual d’estudis poblacionals Els estudis poblacionals o genètica de poblacions, és una branca de la genètica que intenta descriure la variació i distribució de la freqüència al·lèlica per explicar els fenòmens evolutius. Per això es defineix a una població com un grup d’individus de la mateixa espècie que estan 8 Un k-mer és una subseqüència de mida k extreta d’una seqüència més gran. Una de les possibilitats per a l’anàlisi de seqüències d’ADN és generar totes les subseqüències de mida k per posteriorment recopilar dades de repeticions o divergències entre subseqüències. 9 aïllats, comparteixen el mateix hàbitat i es reprodueixen entre ells. Així doncs les poblacions al llarg del temps pateixen diverses mutacions, afectades per la selecció natural per exemple, i l’objectiu dels estudis poblacionals és la d’estudiar les similituds i diferències que es produeixen dins d’una mateixa població per poder arribar inclús a comparar diferents poblacions. Pel que fa al nostre projecte, el que farem és consolidar la idea del Suffix Tree tot aplicant-lo també en aquest tipus d’estudis. Per la manca de temps, no farem un anàlisi ni una aplicació dedicada a aquest objectiu sinó que modificarem algunes parts de la nostra aplicació per a que sigui viable l’ús de quasi la mateixa eina. La mateixa essència del títol ens fa indicar que el programa ha de suportar un gran volum de dades ja que ara no estarem parlant d’un sol individu sinó de varis i que per tant posarem a prova l’escalabilitat de l’aplicació. Com hem comentat anteriorment aquesta part de la genètica està en augment també gràcies a les millores en el camp de la informàtica en quant a volum de dades, i sobretot a la paral·lelització que permet reduir el temps d’execució significativament. No és d’estranyar doncs que el nombre d’aplicacions relacionades amb els estudis poblacionals no pari de créixer, aquí tenim una llista d’algunes d’elles:  Churchill [11]  PyPop [12]  PopGenome [13]  Stacks [14]  Genome STRiP 1000g [15]  ... Dia a dia la llista va en augment però de moment no tenim cap constància de cap altre projecte que hagi utilitzat una implementació de Suffix Tree o de Suffix Array. 1.4.3 Eines utilitzades: MPI, Suffix-Tree, C++ MPI MPI (Message Protocol Interface) és un sistema d’enviaments de missatges estàndard i portable dissenyat per un grup d’investigadors tant d’empresa com acadèmics per a funcionar en un gran nombre d’ordinadors paral·lels. Està orientat per a sistemes amb memòria distribuïda. El seu estàndard defineix una sintaxis i una semàntica que són el nucli de les rutines de la llibreria a més és portable a varis llenguatges com Fortran, C++, C o Java, tot i que els més comuns són els dos primers. La gran avantatge respecte a altres llibreries d’enviament de missatges és que aquesta ha set implementada per a quasi tots els sistemes d’arquitectura de memòria distribuïda gràcies a la seva velocitat i portabilitat l’han fet ser un dels més utilitzats en sistemes d’aquest tipus de memòria. Hem escollit aquesta tecnologia perquè és de les més utilitzades en el camp de memòria distribuïda, a més és compatible amb el llenguatge que utilitzem i ràpida en quant a velocitat. Una altra avantatge és que ja està instal·lat al MareNostrum III, en concret la versió Open MPI 1.8.1. 16 2.2 Explicitació de l’ordre lògic de les tasques Tasca/Subtasca Tasca/Subtasca Predecessora Abast del projecte - Planificació temporal Abast del projecte Gestió econòmica i sostenibilitat Planificació temporal Presentació preliminar Gestió econòmica i sostenibilitat Contextualització i bibliografia Presentació preliminar Plec de condicions Contextualització i bibliografia Presentació oral i document final Plec de condicions Anàlisi del model i algoritme Suffix Tree Presentació oral i document final Anàlisi del model MPI Anàlisi del model i algoritme Suffix Tree Familiarització amb l’entorn de treball i execució Anàlisi del model MPI Implementació Familiarització amb l’entorn de treball i execució Anàlisi i millores Implementació Redacció memòria Implementació Presentació Redacció memòria Taula 1: Explicitació de l'ordre de les Tasques/Subtasques 17 2.3 Diagrama de Gantt A continuació veiem el diagrama de Gantt de la planificació final on les tasques en verd il·lustren les tasques que engloben altres tasques: Il·lustració 2: Diagrama de Gantt de la Planificació final 18 L’ordre de les tasques anteriors és la següent, a més en la imatge s’il·lustren les durades, la data inicial i final i els predecessors de cada tasca: Il·lustració 3: Característiques de les tasques 19 2.4 Anàlisi de les modificacions en la planificació inicial Durant la realització s’han produït desviacions respecte a la planificació inicial, la més important de totes ha set que inicialment vàrem preveure de fer l’aplicació utilitzant el paradigma de MapReduce però al començar a investigar en el tema ens va sorgir la idea d’utilitzar MPI que ha set l’elecció final del projecte degut a que les complicacions que produïa la inicialització de MapReduce ens semblaven que contrarestaven forces respecte a MPI que no requereix de cap inicialització complexa. El fet és que la planificació inicial com hem dit es va basar en MapReduce que finalment no hem utilitzat, però que en un futur ens agradaria donar-li una oportunitat val a dir, per tant les dues planificacions divergeixen bastant en conceptes però en el fons es tracten de dues planificacions semblants però amb tecnologies diferents. Una altra desviació important és que inicialment volíem incloure la nostra aplicació dins el flux de dades de SMUFIN però de moment no serà així ja que l’objectiu ha canviat a crear una nova versió de SMUFIN i per tant aquest part també divergeix de la planificació inicial ja que ara no es requerirà ni fer l’estudi ni moure l’aplicació dins SMUFIN. Una de les desviacions (petita) que més hem patit ha set produïda degut a la demanda d’execucions al MareNostrum III. Degut a la demanda feta per a cada execució de nombre de processadors, ens hem trobat que en molts casos les execucions tardaven un dia en executar-se i per tant mentre es posava una execució a la cua d’execucions s’havia d’aprofitar el temps en que aquesta no s’executava en altres tasques. Pel què fa a les tasques de familiarització amb l’entorn de treball, les d’estudi dels models (canvi de MapReduce a MPI) utilitzats i l’etapa final de redacció i presentació del treball no han patit cap canvi. 2.5 Valoració d’alternatives i pla d’acció Gràcies a l’aplicació de la metodologia Scrum, que és molt dinàmic, ens ha permès que tot i que es produïssin aquests canvis abans mencionats no s’hagi vist alterat el treball de fons amb els mencionats Sprints 11 i les reunions setmanals amb un dels codirectors, ja que cada setmana es marcaven un objectiu per la setmana següent i per tant hem anat modelant el treball a fer de manera seqüencial i regular. El pla d’acció que ha sigut necessari és el següent: Si en una de les Reunions setmanals es detecta una desviació de la planificació caldrà:  Si una tasca s’ha finalitzat abans del que s’havia establert, no hi ha cap problema, s’avançarà a la següent tasca.  Si una tasca dura més del planificat, s’allarga i es comença la següent més tard. Si l’endarreriment és molt significatiu caldrà que l’estudiant dediqui més hores i que intenti tenir una solució menys efectiva en termes de temps i consum de recursos de la màquina ja que en una tasca posterior ja es mirarà d’optimitzar. 11 L’Sprint és el període en què es realitza l'increment del producte. 20 Degut a que els recursos entre les tasques són molt semblants, aquestes desviacions no haurien de produir cap problema en aquest apartat. 21 3 Pressupost 3.1 Identificació i estimació dels costos A continuació es detallaran els diferents elements a considerar en l’estimació del pressupost del projecte. Primer es calcularan els costos directes, seguidament els indirectes i per acabar es posarà tot en context i es calcularan contingències i imprevistos. Costos Directes Es farà un càlcul estimat per a cada tasca programada. En moltes de les tasques s’utilitzen recursos gratuïts i que no generen cap tipus de cost, en aquests casos només s’esmentarà els recursos mencionats. En els recursos que generin costos directes es calcularan els preus unitaris, la seva vida útil, la seva amortització estimada i per acabar el preu final en funció dels paràmetres anteriors. Es calcula una mitjana de 253 dies laborables per any en funció dels pròxims 4 anys que s’utilitzarà pel càlcul de l’amortització i es calcularà en base a 8 hores laborables. El cost del supercomputador MareNostrum III en les tasques que s’utilitza es posarà el preu unitari però no es calcularà l’amortització ni el preu final degut a que és un supercomputador d’ús públic. Els recursos humans del projecte només els constitueixen un estudiant d’enginyeria Informàtica especialitzat en computadors, no es tindran en compte altres actors. El cost per hora és de 20 euros. Gestió del Projecte Els recursos de cost zero que s’utilitzen en aquesta tasca són: l’aplicació Gantter, un visualitzador de PDF, el Racó de la FIB, l’Atenea de la UPC, el Dropbox, el Google Drive. La resta de recursos són els següents: Unitats Preu Unitari (€) Vida Útil Amortització estimada Preu (€) Ordinador Portàtil (inclou ratolí) 1 1000 4 9.26 9.26 Microsoft Office 1 100 3 1.24 1.24 Càmera 1 120 4 1.11 1.11 Microsoft Office Project 1 460 3 5.68 5.68 Recursos humans (Enginyer Informàtic) 75 20 - - 1500 Total 1517.28 Taula 2: Costos directes GEP 22 Estudi previ a la implementació L’estudi previ conté les 3 subtasques següents:  Anàlisi del model i algoritme Suffix Tree (taula2)  Anàlisi del model MPI (taula 3)  Familiarització amb l’entorn de treball i execució(taula 4) Totes tres tasques són molt semblants amb el context de que són tasques d’estudi per a l’estudiant i que per tant no tenen costos molt significatius, les tres comparteixen el següent recurs gratuït: visualitzador de PDF. Els recursos que comporten costos directes es detallen a continuació: Anàlisi del model i algoritme Suffix Tree Unitats Preu Unitari (€) Vida Útil Amortització estimada Preu (€) Ordinador Portàtil (inclou ratolí) 1 1000 4 1.85 1.85 Recursos humans (Enginyer Informàtic) 15 20 - - 300 Total 301.85 Taula 3: Costos directes Anàlisi del model i algoritme Suffix Tree Anàlisi del model MPI Unitats Preu Unitari (€) Vida Útil Amortització estimada Preu (€) Ordinador Portàtil (inclou ratolí) 1 1000 4 1.85 1.85 Recursos humans (Enginyer Informàtic) 20 20 - - 400 Total 401.85 Taula 4:Costos directes Anàlisi del model MPI Familiarització amb l’entorn de treball i execució Unitats Preu Unitari (€) Vida Útil Amortització estimada Preu (€) Ordinador Portàtil (inclou ratolí) 1 1000 4 1.24 1.24 MareNostrum III 1 22700000 3 - - Recursos humans (Enginyer Informàtic) 10 20 - - 200 Total 201.24 Taula 5:Costos directes Familiarització amb l’entorn de treball i execució 23 Implementació de l’aplicació utilitzant MPI Els recursos sense costos d’aquesta tasca són un compilador, el servidor de repositoris (que s’utilitzarà un de gratuït), un compte al MareNostrum III, un editor de text gratuït com el Sublima Text o l’Emacs i un comparador d’arxius gratuït com el diff. Degut a la durada i la dificultat d’aquesta tasca és possible que puguin sorgir imprevistos que comportin que la durada es vegi augmentada i en conseqüència els costos, com per exemple que l’estudiant no aconsegueixi fer que l’algoritme funcioni correctament per alguns casos o que es tardi més del planificat en la comprovació dels resultats de l’algoritme. En tot es pot preveure un 10% de risc en aquesta tasca, i aquest mateix percentatge serà calculat a continuació i mostrat a la següent taula. Unitats Preu Unitari (€) Vida Útil Amortització estimada Preu (€) Ordinador Portàtil (inclou ratolí) 1 1000 4 32.11 32.11 MareNostrum III 1 22700000 3 - - Recursos humans (Enginyer Informàtic) 686 20 - - 13720 Total 13752.11 Imprevistos - 523.21(10%) 1375.21 Total + Imprevistos 15127.32 Taula 6: Costos directes Implementació de l’aplicació utilitzant MPI Anàlisi i millores Els recursos amb cost zero d’aquesta tasca són els següents: un compilador, un servidor de repositoris, un compte al MareNostrum III, un editor de text gratuït com el Sublime Text o l’Emacs, jocs de proves (proporcionades pel BSC i de descàrrega gratuïta), eines gratuïtes que analitzen traces d’un programa com el Tareador i una eina per visualitzar aquestes traces com el Paraver i Extrae, Doxygen per tal de generar la documentació del codi (gratuït) i un comparador d’arxius gratuït com el diff o el cmp. No es preveuen riscos en aquesta tasca degut a que encara no en sabem l’abast. Unitats Preu Unitari (€) Vida Útil Amortització estimada Preu (€) Ordinador Portàtil (inclou ratolí) 1 1000 4 7.35 7.35 MareNostrum III 1 22700000 3 - - Recursos humans (Enginyer Informàtic) 60 20 - - 1200 Total 1207.35 Taula 7: Costos directes Anàlisi i millores 24 Redacció de la memòria Aquesta tasca tindrà els següents recursos a cost zero: un visualitzador de PDF, el Racó de la FIB, el Dropbox i l’editor de text Microsoft Office. Unitats Preu Unitari (€) Vida Útil Amortització estimada Preu (€) Ordinador Portàtil (inclou ratolí) 1 1000 4 6.18 6.18 Microsoft Office 1 100 3 0.82 0.82 Recursos humans (Enginyer Informàtic) 70 20 - - 1400 Total 1407 Taula 8: Costos directes Redacció de la memòria Presentació Per acabar l’última tasca tindrà els següents recursos a cost zero: un visualitzador de PDF, el Racó de la FIB i el Dropbox. Unitats Preu Unitari (€) Vida Útil Amortització estimada Preu (€) Ordinador Portàtil (inclou ratolí) 1 1000 4 6.18 2.18 Microsoft Office 1 100 3 0.82 0.82 Recursos humans (Enginyer Informàtic) 20 20 - - 400 Total 403 Taula 9: Costos directes Presentació Costos Indirectes Els costos indirectes d’aquest projecte són comuns en totes les tasques i es calcularan per a totes les tasques a la vegada, són els següents: Unitats Preu Unitari (€/h) Vida Útil Amortització estimada Preu (€) Llum 956 0.05 - - 47.8 Quota ADSL 956 0.09 - - 86.04 Total 133.84 Taula 10: Costos indirectes de tot el projecte 25 Pressupost Total del projecte Unitats Preu Unitari (€/h) Vida Útil Amortització estimada Preu (€) Tasca1(directes) 1517.28 Tasca2(directes) 301.85 Tasca3(directes) 401.85 Tasca4(directes) 201.24 Tasca5(directes) 15127.32 Tasca6(directes) 1207.35 Tasca7(directes) 1407 Tasca8(directes) 403 Indirectes 133.84 Total acumulat 20700.73 Contingència 5% Total acumulat 1035.04 Total sense IVA 21735.77 Total amb IVA (21%) 26300.28 Taula 11: Pressupost Total del projecte El factor de risc es calcula en funció de la dificultat i les hores destinades a la tasca, en moltes d’elles és zero degut a que la tasca en sí no té molta dificultat. Per altra banda, s’ha calculat el percentatge de la contingència. Per aquest projecte s’ha estimat un 5%, ja que el pressupost que està a continuació està molt detallat i que molts dels recursos usats no tenen cap cost. Finalment s’aplica l’IVA (21%) sobre el cost acumulat del pressupost. 32 viable de tenir uns arxius que s’utilitzaran una vegada (almenys en l’aplicació que estem construint) replicats en molts nodes i per tant hem de ser realistes. Degut a aquest entrebanc, tindríem que tots els nodes estaran competint per a llegir un mateix arxiu i que per tant estaríem seqüenciant la lectura i llegiríem els arxius tantes vegades com nodes utilitzem, a més caldria aplicar els mateixos filtres a cada node. De manera que al final no tindríem un nivell tant alt de paral·lelització com esperàvem. Hauríem serialitzat tant la lectura dels fitxers com el filtratge de qualitat, o sigui que a sobre és llegirien tantes vegades com nodes tenim i que en tots s’hi haurien d’aplicar els mateixos filtres. Estratègia 3 En aquesta tercera estratègia torna a aparèixer el concepte de master (veure il·lustració següent), podríem dir que és molt semblant a la primera estratègia pel que fa a la lectura d’arxius però la diferència la trobem en la part de la comunicació. Aquí alliberem càrrega al master que només ha d’enviar els reads a un dels altres nodes, anomenats slaves. Com a punts forts doncs tenim que només llegim un cop l’entrada i a la part de la comunicació ja només enviem la informació necessària, havent filtrat la resta al master, i que la comunicació del master no es veu tant saturada com en la primera estratègia. En contraposició ens trobem que si tenim n slaves, si definim el processat dels reads i enviament al següent node com una iteració, llavors l’slave n no començarà a treballar fins a la iteració n, això provoca que inicialment molts nodes estiguin inactius i que perdem eficiència en termes de paral·lelització, sobretot si utilitzem un gran nombre de nodes. Una correcció que podria millorar el problema és prioritzar en primer terme el transport de reads abans que el processat en cada node. Il·lustració 8: Estratègia Inicial 2 33 Un cop vistes les tres estratègies vàrem descartar la segona estratègia ja que la disposició que requeria dels arxius no era la que s’utilitzava al Mare Nostrum III, llavors ens vam quedar entre les dues que tenien master. Teníem la intuïció que la lectura seria més ràpida que el processat dels reads, és per això que ens vam decantar per sobrecarregar al master amb més feina i alliberar càrrega de treball als slaves i per tant vam apostar per la primera estratègia. MPI Pel que fa a MPI, durant l’execució es crea el que s’anomena un MPI_Datatype que no és més que la creació del paquet (com a estructura) que volem enviar de tal manera que es pugui enviar amb MPI. Per a això construíem el paquet molt simplement, es tracta de x paquets (on “x” és la mida del buffer que utilitzem d’enviament) que estant formats del read de mida fixa i de l’identificador. D'aquesta manera s'aprofita al màxim l'enviament ja que l'estructura que creem conté els camps justos que necessitem. Com hem comentat abans l’enviament es fa de manera bloquejant, entre les dues opcions principals de MPI: MPI_Send i MPI_SSend, vam escollir la segona. La diferència entre les dues rau que en la segona no es fa realment l’enviament fins que no rebem el senyal que el node objectiu ha fet primer el MPI_Recv (rebuda bloquejant) corresponent, i per tant ens assegurem que es rebrà correctament. Amb MPI_Send podríen haver-hi complicacions amb tantes comunicació degut als buffers interns que utilitza MPI. Anteriorment també enviàvem una altra senyal que era el nombre de parelles que contenia l'enviament, però vam descobrir que utilitzant l'estructura MPI_Status podíem obtenir aquest valor directament. Per a fer-ho cal utilitzar la comanda MPI_GET_COUNT, que ens indica el nombre d’elements enviats. Un dels paràmetres de MPI_Recv és el de màxim número de paquets de tipus datatype que esperem rebre, en aquest camp doncs només cal posar-hi la mida màxima del buffer que serà el nombre màxim que rebrem, ja que tan el buffer del master i el dels slaves tenen la mateixa mida. Il·lustració 9: Estratègia inicial 3 34 Una possible millora del programa seria la utilització de MPI_Bcast, que en principi està pensat per a reduir el temps que es tarda en comunicar un node amb la resta, com és el nostre cas. Creació d'Arbres Degut a termes biològics 15 els punts que busquem es troben a partir dels 30 (més o menys) nivells de profunditat el que comporta que si trobem una manera d’identificar els 30 nivells que tenim per sobre podem estalviar-nos d’inserir-los en l’arbre. La idea que se’ns va acudir és la de convertir els 30 primers nivells, que correspondrien als 30 primers nucleòtids, a un número, de manera que per qualsevol seqüència es pogués traduir a un nombre i donat aquest nombre també fóssim capaços d’obtenir la seqüència. Aquests 30 nucleòtids els vam anomenar “prefix”, i el fet de guardar el prefix en format d’un número 16 ens permetia reduir la memòria utilitzada. La idea doncs de no inserir els primers 30 nivells ens comporta, que de cop tenim molts “subarbres” i que són identificables gràcies a aquest nombre. Llavors un cop el grup de reads ha arribat a un slave aquest busca en tots els sufixes dins cada read i en el cas que algun d’ells estigui en el rang que ha de processar els insereix dins de l’arbre corresponent, per fer això analitzem per a cada sufix, el seu prefix (els x primers caràcters del sufix, per defecte 30, però es configurable amb el paràmetre --prefix_size) i el converteix amb un nombre. La fórmula és la següent: ∑ 𝑥4∗𝑖 𝑖=𝑚𝑖𝑑𝑎 𝑝𝑟𝑒𝑓𝑖𝑥 𝑖=0 Inicialment vam fer una versió que dividia els prefixes de manera consecutiva: el primer slave és queda els j primers prefixes, el segon els j següents. On j és igual a: 𝑗 =𝑛ú𝑚𝑒𝑟𝑜 𝑑𝑒 𝑝𝑟𝑒𝑓𝑖𝑥𝑒𝑠 𝑛ú𝑚𝑒𝑟𝑜 𝑑𝑒 𝑠𝑙𝑎𝑣𝑒𝑠 Però de seguida ens vam adonar que hi havia un desequilibri pel que fa a la càrrega dels prefixes, i això és degut a que en el nostre ADN apareixen més sovint alguns nucleòtids que d’altres, en concret A i T apareixen amb més freqüència, com es pot veure en la següent il·lustració, que és la distribució dels nucleòtids en el cromosoma vint-i-dos que nosaltres utilitzem: 15 Bàsicament a partir d’entre 25 i 35 nivells es considera que identifiquem un punt dins el genoma (unicitat), ja que si agaféssim per exemple un valor menor, com pot ser 4, amb 4 nucleòtids no seríem capaços de definir la zona on estem ja que el nombre de vegades que 4 nucleòtids en concret (per exemple ATGA) apareixen dins el genoma és molt elevat (poca unicitat). 16 En concret un nombre de 64 bits, anomenat “unsigned long long” en C++, molt inferior a crear 30 nivells dins l’arbre. On x és:  0 si A  1 si C  2 si G  3 si T 35 Il·lustració 10: Repartiment de nucleòtids en les dades d'entrada del cromosoma 22 Vegeu il·lustració següent que reflecteix el repartiment dels sufixes pel cromosoma vint-i-dos, executat en 16 nodes (15 slaves i 1 master): Il·lustració 11: Distribució consecutiva 3,81% 4,39% 5,66%5,79% 5,13% 5,91% 6,54%6,97%6,73% 7,33%7,84%8,24%8,34%8,28% 9,03% 12345678910 11 12 13 14 15 0,00% 1,00% 2,00% 3,00% 4,00% 5,00% 6,00% 7,00% 8,00% 9,00% 10,00% Slaves % prefixes Repartiment consecutiu 36 Com es pot veure la càrrega de treball no està ben repartida ja que tenim que el primer slave té prop del 4% dels reads mentre que el 15 en té un 9%. Per aquest motiu vam voler canviar la manera en què distribuíem els prefixes ja que això té una implicació directe amb l’ús de memòria que es fa en el node i amb la computació que aquest fa. La solució que se’ns va acudir va ser la de distribuir consecutivament però amb granularitat 1, això vol dir anar repartint un per un els prefixes: el primer pel primer slave, el segon pel segon slave ... i quan arribem al final tornem a començar. De manera que el què fem és fer el mòdul del prefix entre el nombre de slaves. El resultat amb les mateixes condicions que l’anterior va ser el següent: Com es pot veure la distribució va ser un èxit total permeten així distribuir la càrrega de treball entre els nodes de manera simple i eficaç amb diferències màximes de 0’05%. Tot i que les proves només es van fer amb un cromosoma, en la posteritat s’ha comprovat que es continua mantenint aquesta distribució correctament, si s’utilitza un nombre no potència de 2 com a número de slaves. Per a distribuir els arbres dins els slaves utilitzem un diccionari (o taula de hash). Per a designar aquest estructura vam pensar primer, en com serien els accessos als elements d'aquesta estructura. La resposta és fàcil són accessos aleatoris, això s’explica veient que mentre construïm l’arbre, l’ordre en què van arribant els reads és tal i com estant al fitxer d’entrada i per tant no es garanteix que vinguin de forma ordenada (s’entén ordenada com que els prefixes siguin consecutius), i a més per a cada read busquem tots els sufixes i els prefixes d’aquests tampoc estan ordenats, per tant com el seu nom indica farem accessos aleatoris al diccionari. L’altre punt que buscàvem a l’estructura era que no utilitzés molta memòria, que sempre és un dels nostres objectius, i per tant destinar el màxim de memòria en els arbres. En definitiva necessitàvem una estructura que en termes de memòria no fos molt complexa, no volíem sobrecarregar la memòria amb l’estructura externa i per altra banda que potenciés els accessos 6,67%6,66%6,66%6,68%6,67%6,65%6,66%6,70%6,66%6,66%6,67%6,67%6,67%6,67%6,66% 12345678910 11 12 13 14 15 0,00% 1,00% 2,00% 3,00% 4,00% 5,00% 6,00% 7,00% 8,00% 9,00% 10,00% slaves %prefixes Repartiment consecutiu, granularitat 1 Il·lustració 12: Distribució consecutiva, granularitat 1 37 aleatoris en termes de velocitat d’accés i d’inserció. Amb aquest últim paràmetre com a principal punt a favor ens vam decidir en utilitzar el diccionari o taula de hash 17 . Per a definir l'estructura final vam escollir que la clau, fos el nombre que utilitzàvem per definir el prefix, ja que és una propietat única per a cada arbre i que necessitàvem guardar sí o sí. Per al camp valor que és el camp que realment té la informació que busquem vam posar-hi, no l'estructura de l'arbre directament, sinó simplement un punter a aquesta estructura ja que així evitem sobrecarregar l'estructura. Durant l’elecció de l’estructura que utilitzaríem primer ens vam declinar cap a la taula de hash per defecte de gcc que s'anomena unordered_map però ens vam adonar que utilitzava bastanta memòria i buscant alternatives ens vam trobar amb un anàlisis [20] que bàsicament era un benchmark de taules de hash en què es comparaven diferents aspectes de l'estructura, com per exemple temps de lectura en accessos seqüencials, temps d’insercions, memòria utilitzada... Així que ens vam centrar en les gràfiques de temps d'inserció aleatòria i memòria utilitzada i en les dues teníem una vencedora clara diferent: en termes de insercions aleatòria: dense_hash_map (de Google) i per l'altra banda en memòria utilitzada: sparse_hash_map (també de Google). Com hem fet durant tot el projecte escollirem basant-nos primer amb la memòria utilitzada i és per això que finalment vam escollir la sparse_hash_map. Tot i que en l'anàlisi s'explica que és dues o tres vegades més lenta que l'altra alternativa, també es comenta que és 5 o 6 vegades menor en l’ús de memòria. Així que resumint el què fa un slave al rebre un conjunt de reads és: mirar per a cada read tots els sufixes i si un d'ells l'hem d’incloure mirem si existeix l'arbre amb el mateix prefix dins l'estructura de hash, si no existeix el creem, i inserim el sufix dins l'arbre. Obtenció de resultats i validació de les mutacions dins l'arbre Un cop vam arribar a aquest punt vam adonar-nos que ens faltava una manera de validar l'Arbre (el conjunt de tots els arbres dels slaves), i comprovar que contenia totes les dades que havíem inserit. Com que estàvem treballant amb el cromosoma vint d'un «pacient irreal» ja que el material genètic amb el que treballàvem s'havia creat artificialment, sabíem en tot moment on hi havien les mutacions i de quin tipus es tractaven, de manera que ens va ajudar molt en la validació de l'aplicació. Així que gràcies a que també teníem els arxius BAM del cromosoma, que ens permetien veure quins reads estaven involucrats en cada mutació, se'ns va acudir de fer l'estratègia que nosaltres vam anomenar Read-Oriented. Read-Oriented L’estratègia Read-Oriented es basa en la idea com hem comentat de que estem provant l’aplicació amb un joc de proves, el genoma “In silico”, del qual ja coneixem les solucions i per tant ens permet validar el resultat fàcilment, l’objectiu doncs d’aquesta estratègia és només la de validar la creació de l’arbre i extreure dades del coverage en uns determinats punts. 17 Una taula de hash és una estructura que associa dos conceptes: claus i valors. Per a fer-ho el què fa és: per a cada clau li aplica una “funció de hash”, que no és més que una transformació en un número únic que identifica una posició, i per tant amb aquest número accedim el valor que li correspon a la clau. 38 El procés consisteix en buscar determinats sufixes que intervenen en una mutació i en concret aquells que no tenen la mutació en el seu prefix. Un cop trobats els punts dins l’arbre ens centrem en veure quin nombre de sufixes conservem, així que podrem comprovar si encara mantenim el mateix nombre de coverage i el més important si els reads associats a als punts són els correctes i per tant estan correctament inserits dins l’arbre, fet que vol dir que l’hem generat correctament. Per a realitzar el procés necessitem crear uns arxius d’entrada, nosaltres hem decidit dedicar un arxiu per a cada mutació i per tant tindrem tants arxius com punts de mutació. A l’interior dels arxius tindrem els reads tumorals associats a la mutació, i per cada sufix la posició de la mutació en el seu interior. L’encarregat de llegir els arxius serà com és natural el master, llavors per cada arxiu el master el llegirà, i enviarà als slaves: el sufix a buscar i el punt de mutació. Els slaves per a cada read que rebin hauran de buscar tots els sufixes fins que el prefix arribi al punt de mutació, això es tradueix en que buscarem en tot l’arbre els punts en que tant els reads normals com els tumorals coincideixen i es divergeixen ja que és el punt exacte de la mutació. Imaginem que tenim la següent seqüència i que el punt vermell és el punt on es produeix la mutació [veure il·lustració següent]: Llavors en l’arbre buscarem tots els primers sufixes, en què com a molt el prefix arribi a la “T” abans de la mutació, no té sentit fer-ne més ja que com que a partir d’aquest la mutació estarà dins el prefix farà que els sufixes del tumoral i el normal no caiguin a la mateixa branca, per exemple amb un prefix de 5: el normal seria GCATT mentre que el tumoral GCATG, prefixes diferents, branques diferents i nosaltres el què busquem és el punt just on divergeixen normal i tumoral. Suposant un cas idoni el resultat esperat seria: Com es pot veure suposant un coverage de 30, significaria que no hem perdut cap read, per suposat aquest és un cas ideal i és complicat que aparegui a la realitat. Per tal de visualitzar el Read normal: AGCAACGACGACGCATTCCAT Read tumoral: AGCAACGACGACGCATGCCAT Il·lustració 13: Read normal i tumoral en un punt de mutació tipus Point Mutation Il·lustració 14: Cas idoni d'un punt de mutació 39 resultat de la cerca el què és, primer com que hem inserit tots els sufixes que anem a buscar ja que els reads que busquem provenen de la mateixa entrada, sabem que en tots tindrem algun resultat i que per tant el master espera rebre resposta de tots els sufixes. Per això el què fem és que el master fa un MPI_Recv per a cada sufix i per a identificar-los utilitzem el MPI_Tag ja que ens és independent quin slave ens ho enviï. La forma en que rebem la resposta és: primer rebem 10 comptadors que identifiquen els comptadors del node “pare”, en el nostre cas la T (30,30) i els altres 8 corresponents a cada “fill”. Seguidament rebem el nombre i el valor dels id, que identifiquen tots els reads involucrats en aquell punt. La forma en que presentem els resultats és la següent: A la part de dalt de l’arxiu hi posem la mateixa capçalera que a l’arxiu d’entrada, primer el punt de mutació segons el genoma de referència, després els coverage originals (obtinguts del BAM) després el sufix original que hem buscat i la posició de la mutació dins seu. Seguidament ja venen tots els sufixes ordenats de forma descendent en funció de la seva llargària. Per a cada sufix tindrem en aquest ordre: la mida del sufix, el sufix buscat, els 10 comptadors que corresponen com hem dit als comptadors del pare i dels fills, i els ids involucrats, on per a cada id tindrem l’arxiu i la línia dins l’arxiu on està el read que el conté: Tree-Oriented Com hem vist, l’estratègia del Read-Oriented només ens serveix pel cas del genoma “In silico”, ja que sabem els sufixes exactes que hem de buscar per a trobar els punts de mutació. A la vida real doncs aquesta estratègia no és viable, i per això n’hem hagut d’idear una altra. La idea que hi ha darrera és molt senzilla, això si, que sigui senzilla no vol dir que sigui fàcil d’implementar ja que com veurem s’han de buscar uns filtres que permetin trobar tots els punts de mutació però intentant sempre de reduir el màxim possible els falsos positius. La idea doncs, no és més que fer una cerca en tot l’arbre i aplicar uns determinats filtres en cada punt per determinar si aquest és un punt de mutació o no, per a fer-ho doncs ens fixarem principalment amb els comptadors del node en qüestió i els dels seus fills. Pel que fa a la obtenció de resultats, tenim un problema i és que hem de trobar uns filtres adequats que ens treguin els falsos positius però que no filtri cap dels punts correctes. Per Il·lustració 15: Capçalera de l'arxiu de sortida Il·lustració 16: Format d’exemple d’arxiu de sortida 40 aquesta idea, també ens hem servit que tenim els resultats del cromosoma amb el que treballem i per tant és intentar idear uns filtres partint dels punts trobats abans amb l’estratègia de ReadOriented amb la gràcia que ara hem utilitzat aquesta mateixa estratègia en reads aleatoris, que sabem que no contenen cap mutació. D’aquesta manera podem provar diferents filtres mirant en tot moment el nombre de mutacions que aquests filtren i intentar ajustar-los al màxim per a filtrar els punts aleatoris, que serien falsos positius. Així doncs ens guiem sobretot amb la hipòtesis que el genoma In silico pot simular un genoma natural i per tant que els passos que anem seguint i corroborant amb aquest són aplicables a la realitat. Val a dir però que hem de tenir en compte que no ho estem provant amb l’arbre degut a la facilitat que ens aporta treballar només amb aquests punts, però tenim en compte que els falsos positius poden ser més nombrosos en l’arbre ja que només n’hem agafat un subconjunt. De moment encara estem treballant en buscar uns bons filtres, en aquesta part hem tingut l’ajuda d’alguns dels biòlegs i biotecnòlegs del grup on ens trobem ja que aquí podríem dir que es tornen a encreuar la feina d’informàtic i d’altres disciplines. 4.3 Anàlisi Un cop acabada la primera versió vam instrumentar-la per a poder avaluar-ne el rendiment i veure els seus punts , que serien els que prioritzaríem alhora de buscar optimitzacions. Degut al gran detall que volíem arribar amb l’anàlisi de l’execució per tal de no crear arxius per analitzar immensos hem decidit ajustar els paràmetres del programa per que s’adeqüin millor. En concret hem reduït el nombre de nodes utilitzats i els arxius d’entrada han estat substituïts per uns de més petits. Les traces han estat generades utilitzant l’aplicació Extrae i visualitzades gràcies a l’eina Paraver, les característiques de l’execució es mostren a continuació:  Nombre de nodes: 6 (5 slaves i 1 master)  Total dades d’entrada més reduïda de 15,6 GB (cromosoma 20) a 1790 KB  Readblock : 50 MB  Mida buffer d’enviament: 6,48 MB (60 reads)  A més només entrem fastq1 tant de normal com tumoral i per tant ni fastq2 de normal ni de tumoral A continuació veurem la traça general de l’aplicació, tot il·lustrada amb les funcions i les comunicacions entre els diferents nodes: 41 Cada color identifica una funció dins els codi com podem veure a la llegenda, a continuació explicarem què es fa en cada funció per a tal de poder comentar la gràfica:  End: Té varis significats, la primera part ens indica el pas de paràmetres inicials, mentre que la part del final indica primer que els slaves esperen la senyal de comunicació que indica que s’han acabat els arxius d’entrada (les línies grogues primes representen la comunicació), i un cop rebuda aquesta senyal representa que els nodes (tant master com slaves) eliminen les estructures i s’esperen a acabar tots junts.  Send_reads(master): La funció send_reads és l’encarregada de copiar els reads llegits al buffer d’enviament, i si aquest s’omple enviar-los als slaves  Insert_normal(slaves): Funció encarregada d’inserir un sufix dins l’arbre que li correspon.  Insert_tumoral(slaves): Ídem que l’anterior però pel cas dels sufixes tumorals.  Read_file(master): Funció encarregada com el seu nom indica de llegir un arxiu d’entrada, en concret de llegir blocs d’aquests de mida fixada amb el paràmetre Il·lustració 17: Traça General de l'execució Il·lustració 18: Llegenda de la traça 48 Les mides del buffer han estat seleccionades en funció del nombre d’elements (read + id) que hi caben que és el segon eix que hem indicat a sota del de MB, que va des de 100 fins a 107 elements. En la gràfica s’il·lustra perfectament els dos efectes descrits prèviament veient que s’estabilitza a partir dels 5*103 elements i que comença a augmentar de nou als 5*106. Gràcies a aquest anàlisi ara ja sabem quins són els valors indicats pels nostres buffers de manera que sigui més eficient tant la lectura com l’enviament de dades. 4.4 Millores 4.4.1 Master amb Suffix Tree (Master Tree) Com hem comentat anteriorment un cop vist les estadístiques i les traces de l’execució de la primera versió, sobretot la de l’eficiència, ens hem plantejat la idea que el node master també insereix-hi reads, reduint així la càrrega de treball dels slaves i de retruc augmentar l’eficiència del master, que té un valor bastant baix. L’estratègia que hem seguit és que un cop al master envia els reads als slaves, aquest també insereix els mateixos reads que acaba d’enviar (evidentment inserirà els sufixes que li pertoquin). La taula següent, realitzada en el joc de proves reduït, utilitzat en el cas de la versió inicial, mostra els resultats obtinguts d’eficiència: 0,0108 0,108 0,54 1,08 4,32 6,48 8,64 10,8 54 108 540 1080 0 1000 2000 3000 4000 MB TIME(S) Impacte de la mida del buffer d'enviament (Cro. 20) #Elements 100 103 5*103 104 4*104 6*104 8*104 105 5*105 106 5*106 107 Il·lustració 23: Temps d'execució en funció de la mida del buffer d'enviament 49 Veiem que la taula ha canviat bastant, la gran millora la veiem, com era previsible en el node master que augmenta d’un 2,2% a un 78,7% l’estat de running, tot i que quasi tots els slaves empitjoren el resultat final és d’un 81,97% d’eficiència superant així el 76,9% de la versió inicial. Un cop vista l’estadística vam decidir de provar aquesta nova millora amb el joc de proves gran (cromosoma 20) i malauradament vam veure com els resultats que aquesta obtenia eren inferiors en quant a temps ja que teníem un increment del 81%. És per això que tot i que en els jocs de prova més petits tingui millors resultats (un 5% més ràpid), nosaltres hem de tenir en consideració que l’aplicació s’utilitza per a grans volums de dades i el resultat per aquests és negatiu. En definitiva no hem inclòs aquesta optimització en l’aplicació final. 4.4.2 Partitions Una altra estratègia que se’ns va acudir va ser l’anomenada partitions. Aquesta es basa en la idea que quan analitzem l’arbre per a buscar possibles punts de mutació, mirem un punt en concret de l’arbre i la resta de l’arbre no té cap implicació en aquell moment o sigui que és independent de la resta de l’arbre. Gràcies a aquesta característica, podem construir per exemple la meitat de l’arbre, analitzar-la i un cop analitzada, eliminar aquesta part i fer el mateix per l’altra part. O anar més endavant i dividir l’arbre en tantes parts com vulguem. Aquesta és la part que nosaltres anomenem partitions, el nombre de particions que aplicarem a l’arbre i que és un paràmetre configurable. D’aquesta manera si per exemple posem un partitions de 5, estem dividint la memòria que necessitem en 5 (cas ideal) i per tant estem eliminant la limitació de memòria que poden tenir algunes màquines per a executar l’aplicació, val a dir però que implica un increment en el temps que tarda l’aplicació degut entre altres que s’han de llegir les entrades x vegades (5 en l’exemple) i que s’han d’eliminar les estructures creades cada x vegades també. Així podríem dir que aquesta optimització ens permet alternar Taula 17: Estat dels processadors durant l'execució amb Master Tree 50 la memòria necessitada per l’aplicació segons el requeriment de les màquines on s’executa tot saben que una reducció d’aquesta implica un increment pel que fa al temps. A més aquesta característica podem aplicar-la als estudis poblacionals, tot i l’impacte en el temps d’execució, la memòria ja no serà el problema principal a tenir en consideració, que fins ara era el coll d’ampolla principal que podíem trobar-hi. Per a poder realitzar aquesta nova característica dividim del rang total de prefixes que existeixen entre el nombre de partitions que s’ha configurat, d’aquesta part n’anomenen un bloc de treball. En cada bloc de treball tenim un prefix mínim i un de màxim que identifiquen el conjunt de prefixes a tractar en aquella iteració, tenint en compte que l’última pot tenir un pèl més de càrrega de treball. Llavors quan els slaves reben un read i en revisen els prefixes el primer que comproven és si està dins els rang d’aquell bloc de treball. En definitiva aquesta nova característica de l’aplicació farà que sigui més portable, al poder modificar el consum de memòria utilitzat (més partitions -> menys memòria requerida) permeten així d’executar l’aplicació en clústers que no tinguin tanta memòria disponible (podríem dir que treu la limitació de memòria), el que suposa un gran avanç pel que fa a la part dedicada als estudis poblacionals. A continuació tenim un gràfic amb els resultats d’executar el cromosoma 20 del genoma “In silico” utilitzant 75 nodes. Com a guia hem posat també els valors de la versió anterior per veure que en cas d’utilitzar una sola partició no tenim cap pèrdua de rendiment. Com podem veure el comportament és l’esperat, a mesura que augmentem les particions es redueix la memòria utilitzada i augmenta el temps d’execució. Tot i això ens esperàvem una reducció pel que fa a memòria més elevada. Encara que sembli que el punt òptim és entre 11 i 12 particions, cal recordar que cada una de les sèries de valors té la seva escala diferent al costat i si ens hi fixem veiem que un possible punt òptim seria el 4, ja que té una reducció en el temps comparat amb el seu anterior a part de la de memòria. Creiem que això és producte d’una distribució incorrecte pel que fa a càrrega de treball entre els diferents blocs de treball o particions. A més és important esmentar que l’arbre que utilitzem té un comportament peculiar 0 2000 4000 6000 8000 10000 12000 0 200 400 600 800 1000 1200 1400 1600 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Temps (s) Memòria (GB) Número de particions Partitions Cromosoma 20: 75 nodes Memòria (GB) Memòria versió antiga Temps(s) Temps versió antiga Il·lustració 24: Gràfic d'execució amb partitions 51 quan inserim sufixes, i és que sempre intentem tenir els sufixes comprimits i per tant a mesura que augmentem el nombre de sufixes inserits, no només augmentem la memòria degut a aquesta descompressió que cal fer sinó que també augmentem el temps d’inserció ja que cal anar descomprimint per tal de poder inserir. Així que a l’aplicar l’optimització de partitions, també estem reduint en certa mesura aquest comportament i per tant reduint el temps d’inserció en l’arbre gràcies a l’efecte que acabem de descriure. Una altra comparativa que ens agradaria mostrar no només amb la versió inicial sinó també amb la primera optimització Master Tree, és la de l’eficiència, els resultats extrets per 2, 4, 8 i 16 són els següents: Taula 19: Estat dels processadors, partitions 2 Taula 21: Estat dels processadors, partitions 8 Taula 20:Estat dels processadors, partitions 16 Taula 18: Estat dels processadors, partitions 4 52 El que veiem és que l’eficiència respecte la versió inicial ha augmentat força amb valors semblants a la de l’optimització de Master Tree arribant gairebé als 82%, sobretot quan augmentem el nombre de particions, que no és més que dividir l’arbre en més parts i de retruc fer treballar el master ja que ha de llegir més vegades els arxius. Veiem també que a l’augmentar el nombre de particions estem accentuant el desbalanceig de càrrega que hi ha entre els nodes slaves i que tot i que augmentem el nivell d’eficiència general de l’aplicació és possible que si augmentem en gran mesura les particions al final tindrem que el desbalanceig serà tant alt que provocarà un pèrdua de rendiment general. Pel que fa als temps d’execució, veiem que per aquest joc de proves petit no hi ha gaire diferències, des dels 4,57 segons el més ràpid (1 partició) al 4,85 segons de l’últim (16 particions). O sigui que pel que fa a temps ja tenim la gràfica anterior que ens il·lustra el comportament que aquest té en funció del nombre de particions. En definitiva aquesta optimització creiem que ha reduït en gran importància l’ús de memòria, eliminant teòricament la limitació de memòria dels equips en compensació d’un augment del temps d’execució i com hem comentat abans pensem que és un gran avenç pel que fa als estudis poblacionals, permeten teòricament un numero il·limitat de pacients 20 , i la possibilitat de executar-se en qualsevol màquina sense importar la quantitat de memòria que aquesta tingui. A més també hem aconseguit millorar l’eficiència de l’aplicació sobretot gràcies a augmentar la càrrega de treball que fa el node master, tot i la pèrdua de rendiment d’alguns slaves. 4.4.3 Reverse Una de les dades extretes de la primera versió i que ens preocupaven era la pèrdua de coverage que teníem en els punts investigats, tant els llocs de mutació que havíem buscat com altres punts a l’atzar, tot mantenint el mateix nivell de profunditat que els punts de mutació per poderlos comparar. Per pèrdua de coverage ens referim a tenir uns comptadors molt menors en comparació el coverage inicial, en el cas de l’”In silico” de 30. Nosaltres vam atribuir aquesta pèrdua de coverage a quatre factors, principalment:  Mutació dins prefix  Reads amb contaminació  Mutació heterozigota  2 hèlixs La primera és deguda a que quan un read conté la mutació dins els primers 30 nucleòtids (o sigui dins el prefix), llavors aquest read, no anirà cap al mateix arbre que els normals ja que és diferent en el prefix i per tant no estarà en el punt on busquem. Això provocarà que estarem inserint el sufix en una altra part de l’arbre però que no l’estarem utilitzant, aquests es podrien intentar trobar un cop hem analitzat una mutació i tenim la zona reconstruïda, a l’identificar el punt de mutació buscar tot els sufixes que puguin estar relacionats o sigui que siguin pròxims i buscarlos dins l’arbre amb i sense mutació. 20 Cal mencionar que en els estudis poblacionals no s’acostuma a tenir un gran nombre de coverage i per tant els arxius d’un individu tenen una mida molt més petita. 53 La segona és quan tenim una contaminació en el procés de seqüenciació, que és un procés físic i que per tant té un marge d’error que provoca que el read en qüestió sigui incorrecte, tot i que treballem amb el genoma “In silico”, aquest també té uns errors induïts per tal de que s’assembli a un genoma real seqüenciat. L’altre possible cas és que sigui una mutació heterozigota, que és aquella mutació que només és present en un dels dos al·lels 21 del pacient. Per tant pot ser que tinguem reads que cobreixen una mutació d’un al·lel (que realment té la mutació) i reads de l’altre. Aquest punt només es produït en el cas de les zones tumorals i pot passar tant als comptadors tumorals com els normals. Malauradament no és possible corregir aquest efecte ja que es provocat pel mètode de seqüenciament i per tant inherent en els arxius d’entrada, però com que sabem que existeix podem mirar de suavitzar els filtres de comptadors per tal que no es filtrin aquest tipus de mutacions. Per acabar a més a més, tenim el defecte que durant la seqüenciació s’agafa una de les dues hèlix que formen l’ADN indistintament i aquesta hèlix estan unides d’una manera determinada. I és que tot i ser la mateixa informació en les dues, està invertida amb el sentit que quan en una tenim la guanina (G) a l’altra tenim la citosina (C) i recíprocament, també tenim aquesta relació amb la adenina (A) i la timina(T). I a més quan la màquina selecciona l’altra hèlix el read es llegeix al revés. O sigui que a part d’estar capgirada a la vegada també esta invertida, si entenem inversió per el canvi de A a T (i viceversa) i de C a G (i viceversa). Per això nosaltres vam pensar en inserir el que anomenem el reverse, que no és més que per cada read, crear un read revers que vol dir capgirar-lo i canviar A per T, T per A, C per G i G per C (o sigui invertint-lo). D’aquesta manera augmenten en gran quantitat el volum de dades que introduïm al programa però recuperem part del coverage. Els resultats van ser els següents: Cromosoma 20 Arxius BAM (originals) Suffix Tree Normal Tumor Normal Tumor Sense Revers 29.79 29.98 10.89 10.32 Amb Revers 58.79 60.02 21.22 21.12 Taula 22: Valors de coverage amb revers i sense dins l’arbre Veiem doncs que hem complert l’objectiu que teníem d’augmentar el coverage, ja que l’hem duplicat en tots els casos. El problema és que per això hem hagut de duplicar l’entrada del programa que ens ha fet duplicar la memòria requerida i quasi triplicar el temps d’execució. És per això que tot i assolir l’objectiu, deixarem l’opció dins el programa però estarà inhabilitada per defecte. 21 Cadascuna de les formes alternatives que pot presentar un gen que ocupa el mateix lloc en un cromosoma determinat o en dos cromosomes homòlegs, i que expressa diferentment un mateix caràcter. 54 4.4.4 Lectura d’arxius comprimits És molt comú en aquest camp, degut a que els arxius d’entrada són de grans dimensions, que estiguin comprimits és per això que vàrem creure necessari la característica de poder llegir arxius d’aquest tipus, en concret comprimits amb Zip (amb extensió per exemple .gz). Així que vam crear una nova funció de lectura de fitxers fastq totalment renovada que ens permetia llegir els arxius comprimits directament sense necessitat de descomprimir-los tots sencers, sinó utilitzant l’estratègia que havíem emprat anteriorment d’utilitzar el buffer de lectura. Com que aquest algoritme era més lent vàrem decidir de mantenir l’algoritme antic i utilitzar aquest nou en el cas només dels arxius comprimits. A continuació veiem una gràfica de temps de la lectura dels fitxers en funció de les mides de readblock, tant l’algoritme inicial com el nou, això si, el nou utilitza arxius comprimits: Com hem comentat el rendiment del nou algoritme és pitjor, val a dir però que tracta amb arxius comprimits i que si ho volguéssim fer amb la versió inicial hauríem de descomprimir les dades prèviament a l’execució i per tant tot i que sembli pitjor de moment és l’única opció que tenim. 4.4.5 Comunicació Master-Slaves La versió inicial té una versió molt bàsica pel que fa a comunicació entre nodes és per això que vam decidir de millorar aquest apartat, per a fer això vam pensar en dues crides diferents: MPI_Ibcast Vam estar desenvolupant un versió que utilitzava la crida MPI_Ibcast, que és fer un Broadcast, però de manera asíncrona, la complexitat alhora de programar es bastant alta perquè és necessari crear i administrar un buffer de sortida que s’actualitzi a mesura que o bé el master acaba un bloc de lectura o bé tots els nodes han informat que s’ha llegit un bloc i per tant es pot eliminar. Tot i que vam arribar a implementar una versió inicial que funcionava amb arxius de mida reduïda, no vam ser capaços de poder-la utilitzar en arxius més grans i a més la poca informació que vàrem trobar de la funció ens va fer sospitar que no devia ser tant bona respecte a la seva homòloga síncrona. Per altra banda actualment estem treballant amb una versió que utilitza la funció MPI_Bcast per a realitzar les comunicacions entre els nodes ja que és molt més eficient en quant a temps i ens permetria millorar per tant el rendiment de l’aplicació. La dificultat que aporta aquesta funció és que tots els nodes tenen que realitzar la crida en el mateix punt i per tant cal fer una fusió de funcions de master i slaves per tal que puguin arribar a fer-la. 0,00030,00050,0010,0030,005 0 , 01 0, 0 3 0,05 0,1 0,3 0,6 1 2 5 10 20 50 60 0 500 1000 1500 SIZE(MB) TIME (SEC) TEMPS DE LECTURA DE FITXERS versió1 versió2_gz Il·lustració 25: Gràfica de temps de la lectura de fitxers en les diferents versions 55 4.5 Resultats Com hem comentat abans per tal de capturar resultats hem hagut d’implementar uns filtres amb l’ajuda dels companys del grup multidisciplinari, per això a continuació presentarem alguns resultats amb alguns d’aquests filtres i combinacions d’alguns d’ells, recordem que estem utilitzant l’estratègia anomenada Tree Oriented. Per començar presentarem les dades del cromosoma amb el què hem treballat: Cromosoma 20 Point Mutations: 168 Inversions 6 Deletions 12 Insertions 18 Taula 23: Mutacions en el cromosoma 20 Degut al gran nombre de les Point Mutations ens hem decidit d’extreure els resultats d’aquestes ja que en els altres casos el nombre és massa petit per tal d’intentar extreure algun patró que els diferenciï. Les Point Mutations és la mutació més simple, es tracta del canvi d’un sol nucleòtid, ja sigui perquè ha canviat a un altre nucleòtid, perquè se n’ha eliminat un o perquè se n’ha afegit un de més. El primer filtre que vam dissenyar és el següent: valor ≤ 100*𝑇𝑖 𝑇𝐴+𝑇𝐶+𝑇𝐺+𝑇𝑇 On Ti és el comptador tumoral en un punt i TA, TC ... són els comptadors tumorals dels fills. El resultat que teníem en funció del valor del filtratge és el següent: Il·lustració 26: Resultat d'aplicar el 1r filtre Com veiem el nombre de mutacions que filtrem augmenta a mesura que el filtre és més rigorós, com per altra banda era previsible, d’aquesta manera veiem que el valor òptim estaria entre 2 i 0 3 6 9 12 15 18 21 24 27 30 33 36 39 42 45 48 51 54 57 60 63 66 69 72 75 78 81 84 87 90 93 96 99 0 10 20 30 40 50 60 70 80 90 100 Percentatge de mutacions no filtrades 56 5. El problema que se’ns presenta és que si escollim un nivell tant baix és molt probable que el nombre de nodes que passaran el filtre serà tant gran que realment no haurà tingut molt impacte l’aplicació d’aquest filtre. Per això vam voler repetir aquest filtre amb aquell conjunt de punts aleatoris que no contenien mutació per tal de veure quin era el “gruix” de punts no tumorals que realment filtràvem amb aquest filtre. El comportament és el següent: Il·lustració 27: Resultat d'aplicar el 1r filtre en els dos tipus de posicions Veiem doncs que en realitat si que estem filtrant els nodes naturals, per exemple si agafem el valor de 20, estem mantenint un 92,16 % de les mutacions i en canvi dels punts de no mutacions els estem reduint al 30%. Tot i així com hem vist abans si establim el filtre a 20 estem perdent un nombre de mutacions massa elevat pels objectius que ens vam marcar. És per això que en comptes d’aplicar aquest filtre amb el valor de 20 vam decidir d’aplicar-lo a 4 i buscar algun altre filtre que solucionés el volum de punts de no mutació. El filtre que vam pensar per eliminar els falsos punts de mutació va ser el següent: 𝑣𝑎𝑙𝑜𝑟 ≤ 100∗𝑇𝑖 𝑁𝑖+𝑇𝑖 Degut a que en els punts de mutació un dels seus fills conté un gran nombre de tumorals en comparació als normals. El filtre no elimina els punts en què un dels seus fills compleix la desigualtat anterior. D’aquesta manera obtenim el següent resultat: 04812 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84 88 92 96 100 0 10 20 30 40 50 60 70 80 90 100 Valor del filtratge % no filtrat no_mutations point_mutations 57 Il·lustració 28: Resultat d'aplicar el 2n filtre en els dos tipus de punts Com veiem el resultat del 2n filtre és un èxit ja que com veiem en la il·lustració anterior ens permet filtrar en gran mesura els punts de no mutació i no perdem un gran nombre de mutacions, veiem doncs per exemple el punt de 60, en què mantenim un 96,41% de les mutacions i reduïm els punts de no mutació fins a un 21’5%. 4.6 Estudis poblacionals Per tal d’elaborar la versió d’estudis poblacionals vam partir del codi de la versió somàtica amb les optimitzacions abans descrites. La principal restricció que ens vam trobar va ser la limitació de temps de manera que les solucions proposades queden lluny del nivell d’optimització que hauríem desitjat que tinguessin. El principal problema que se’ns planteja en aquest apartat és que ara no construirem un arbre a partir de dos genomes i comparant diferències un amb l’altre sinó que a part de que tindrem un nombre més elevat de genomes, aquests hauran de comparar-se tots amb tots dins l’arbre. La manera en que vam idear per solucionar aquest problema va ser la de no utilitzar els comptadors tumorals dels nodes i només utilitzar els comptadors normals. Per tal d’identificar en cada node a quin pacient pertanyien els sufixes vam idear una taula que ens indiqués els identificadors (ids) per a cada pacient com la que veiem a continuació: 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84 88 92 96 100 0 10 20 30 40 50 60 70 80 90 100 VALOR DEL FILTRATGE % NO FILTRAT point_mutations no_mutations Il·lustració 29: Taula d'identificadors dels pacients