Plataforma de supercomputación para bioinformática

AUTOR: Diego Darío Guerrero Fernández http://orcid.org/0000-0001-6749-0962 EDITA: Publicaciones y Divulgación Científica. Universidad de Málaga Esta obra está sujeta a una licencia Creative Commons: Reconocimiento - No comercial - SinObraDerivada (cc-by-nc-nd): Http://creativecommons.org/licences/by-nc-nd/3.0/es Cualquier parte de esta obra se puede reproducir sin autorización pero con el reconocimiento y atribución de los autores. No se puede hacer uso comercial de la obra y no se puede alterar, transformar o hacer obras derivadas. Esta Tesis Doctoral está depositada en el Repositorio Institucional de la Universidad de Málaga (RIUMA): riuma.uma.es tesis doctoral Plataforma de supercomputaci´on para bioinform´atica Diego Dar ´ ıo Guerrero Fern´ andez Universidad de M´ alaga Departamento de Biolog´ıa Molecular y Bioqu´ımica Facultad de Ciencias Plataforma Andaluza de Bioinform´atica Edificio de Bioinnovaci´on M´alaga. Mayo de 2015 Plataforma de supercomputaci´on para bioinform´atica Memoria presentada por: Diego Dar´ıo Guerrero Fern´andez Para optar al grado de Doctor por la Universidad de M´alaga Tesis realizada ba jo la direcci´on del Dr. M. Gonzalo Claros D´ıaz en la Plataforma Andaluza de Bioinform´atica y el Departamento de Biolog´ıa Molecular y Bioqu´ımica de la Universidad de M´alaga Fdo. Diego Dar´ıo Guerrero Fern´andez Vo.Bo.DIRECTOR DE LA TESIS DOCTORAL: Fdo. M. Gonzalo Claros D´ıaz M´alaga, mayo de 2015 2 D. M. Gonzalo Claros D ´ ıaz, Investigador de la Plataforma Andaluza de Bioinform´atica y el Departamento de Biolog´ıa Molecular y Bioqu´ımica de la Universidad de M´alaga. CERTIFICA: Que Don Diego Dar ´ ıo Guerrero Fern´ andez, Ingeniero en Inform´atica, ha realizado bajo mi direcci´on en el Departamento de Biolog´ıa Molecular y Bioqu´ımica y la Plataforma Andaluza de Bioinform´atica de la Universidad de M´alaga, el trabajo de investigaci´on recogido en la presente memoria de Tesis Doctoral que lleva por t´ıtulo: “Plataforma de supercomputaci´on para bioinform´atica”. Tras la revisi´on de la presente Memoria se ha estimado oportuna su presentaci´on ante la Comisi´on de Evaluaci´on correspondiente, por lo que autorizo su exposici´on y defensa para optar al grado de Doctor. Y para que as´ı conste, en cumplimiento de las disposiciones legales vigentes, firmo el presente certificado. M´alaga, mayo de 2015 El Director de la Tesis, Dr. D. M. Gonzalo Claros D´ıaz 3 D. M. Gonzalo Claros D ´ ıaz, Investigador del Departamento de Biolog´ıa Molecular y Bioqu´ımica de la Universidad de M´alaga. INFORMA: Que Don Diego Dar ´ ıo Guerrero Fern´ andez, Ingeniero en Inform´atica, ha realizado bajo mi direcci´on en el Departamento de Biolog´ıa Molecular y Bioqu´ımica y la Plataforma Andaluza de Bioinform´atica de la Universidad de M´alaga, el trabajo de investigaci´on recogido en la presente memoria de Tesis Doctoral que lleva por t´ıtulo: “Plataforma de supercomputaci´on para bioinform´atica”, y que la misma cumple los requisitos de idoneidad necesarios para ser presentada por compendio de publicaciones. La memoria est´a avalada por los siguientes art´ıculos: - AlignMiner: a Web-based tool for detection of divergent regions in multiple sequence alignments of conserved sequences. - Guerrero, D., Bautista, R., Villalobos, D. P., Cant´on, F. R., and Claros, M. G. - Algorithms for Molecular Biology. - SCBI MapReduce, a New Ruby Task-Farm Skeleton for Automated Parallelisation and Distribution in Chunks of Sequences: The Implementation of a Boosted Blast+ - Dar´ıo Guerrero-Fern´andez, Juan Falgueras, and M. Gonzalo Claros - Computational Biology Journal - GENote v.: A Web Tool Prototype for Annotation of Unfinished Sequences in Nonmodel Eukaryotes - Bioinformatics for Personalized Medicine - No´e Fern´andez-Pozo, Dar´ıo Guerrero-Fern´andez, Roc´ıo Bautista, Josefa G´omez-Maldonado, Concepci´on Avila, Francisco M. C´anovas, M. Gonzalo Claros - Lecture Notes in Computer Science Y para que as´ı conste, en cumplimiento de las disposiciones legales vigentes, firmo el presente informe. M´alaga, mayo de 2015 El Director de la Tesis, Dr. D. M. Gonzalo Claros D´ıaz 4 Proyectos de investigaci´on Incremento de la eficiencia en el uso del agua en Vitis vin´ıfera L.: bases gen´eticas y fisiol´ogicas para una mejor adaptaci´on al cambio clim´atico (2014-2017, RTA2013-00068-C03-02, MINECO-INIA). IP: M.G. Claros Implementaci´on de TECnolog´ıas INNOvadoras de mejora gen´etica en lenguado senegal´es (Solea senegalensis) y dorada (Sparus aurata) para la optimizaci´on de su producci´on industrial (INNOTECSS) (2014-2017, RTA2013-00023-C02-01, MINECO-INIA). IP: Manuel Manchado Campa˜na Desarrollo de herramientas bioinform´aticas para los estudios gen´omicos y transcript´omicos a partir de datos de secuenciaci´on de lecturas cortas de alto rendimiento para las especies que no tienen un organismo modelo de referencia (NEOGEN) (1-4-11 a 30-4-2016; Proyecto de Excelencia de la Junta de Andaluc´ıa, P10-CVI-6075). IP: M. G. Claros Arquitecturas, compiladores y aplicaciones en multiprocesadores - TIN2010-16144 [20112013] Ip: E. L´opez Zapata y Oscar Plata Genomic tools in maritime PINE for enhanced biomass production and SUSTAINable forest Management (SUSTAINPINE). (2010-2013; MICINN and FP7-PLANT-KBBE Scientific Advisory Board PLE2009-0016). IP: Francisco C´anovas Art´ıculos y cap´ıtulos de libros Canales, J., Bautista, R., Label, P., G´omez-Maldonado, J., Lesur, I., Fern´andez-Pozo, N., . . . C´anovas, F. M. (2014). De novo assembly of maritime pine transcriptome: Implications for forest breeding and biotechnology. Plant Biotechnology Journal, 12(3), 286–299. http://doi.org/10.1111/pbi.12136 Benzekri, H., Armesto, P., Cousin, X., Rovira, M., Crespo, D., Merlo, M., . . . Manchado, M. (2014). De novo assembly, characterization and functional annotation of Senegalese sole (Solea senegalensis) and common sole (Solea solea) transcriptomes: integration in a database and design of a microarray. BMC Genomics, 15(1), 952. http://doi.org/10.1186/14712164-15-952 Dar´ıo Guerrero-Fern´andez, Juan Falgueras, M. Gonzalo Claros (2013). SCBI MapReduce, a New Ruby Task-Farm Skeleton for Automated Parallelisation and Distribution in Chunks of Sequences: The Implementation of a Boosted Blast+. Computational Biology Journal: 10/2013; 2013. DOI:10.1155/2013/707540 Dar´ıo Guerrero-Fern´andez, Rafael Larrosa, and M. Gonzalo Claros (2013). FQbin a compatible and optimized format for storing and managing sequence data. IWBBIO 2013, page 337-344. 5 No´e Fern´andez-Pozo, Dar´ıo Guerrero-Fern´andez, Roc´ıo Bautista, Josefa G´omezMaldonado, Concepci´on Avila, Francisco M. C´anovas, M. Gonzalo Claros (2012). GENote .1: A Web Tool Prototype for Annotation of Unfinished Sequences in Non-model Eukaryotes. Bioinformatics for Personalized Medicine - 10.1007/978-3-642-28062-7 7 Claros, M. G., Bautista, R., Guerrero-Fern´andez, D., Benzerki, H., Seoane, P., and Fern´andez-Pozo, N. (2012). Why Assembling Plant Genome Sequences Is So Challenging. Biology. http://doi.org/10.3390/biology1020439 Fern´andez-Pozo, N., Canales, J., Guerrero-Fern´andez, D., Villalobos, D. P., D´ıaz-Moreno, S. M., Bautista, R., ... and Claros, M. G. (2011). EuroPineDB: a high-coverage web database for maritime pine transcriptome. BMC genomics, 12(1), 366. Guerrero, D., Bautista, R., Villalobos, D. P., Cant´on, F. R., and Claros, M. G. (2010). AlignMiner: a Web-based tool for detection of divergent regions in multiple sequence alignments of conserved sequences. Algorithms for Molecular Biology : AMB, 5, 24. http://doi.org/10.1186/1748-7188-5-24 Art´ıculos en preparaci´on o revisi´on Rosario Carmona, A. Zafra , Pedro Seoane, A. Castro, Dar´ıo Guerrero-Fern´andez,Trinidad Castillo, Ana Medina-Garc´ıa, Francisco M. C´anovas, Jos´e F. Aldana-Montes, Ismael NavasDelgado, Juan D. Alch´e , M. Gonzalo Claros - ReprOlive: a Database with Linked-Data for the Olive Tree (Olea europaea L.) Reproductive Transcriptome - Frontiers in Journal (2015). Guerrero-Fern´andez, D., Bocinos, A., Bautista, R.Fern´andez-Pozo, Juan Falgueras and Claros, M. G. SeqTrimNext: pre-processing sequence reads for next-generation sequencing projects. Guerrero-Fern´andez and Claros, M. G. InGeBIOL: A web interface generator for command line tools. Fern´andez-Pozo, N., Guerrero-Fern´andez, D., Bautista, R. and Claros, M. G. FULLLENGTHERNEXT: A tool for fine-tuning de novo assembled transcriptomes of non-model organisms. Comunicaciones orales en congresos Dar´ıo Guerrero-Fern´andez, No´e Fern´andez-Pozo, Almudena Bocinos, Roc´ıo Bautista and M. Gonzalo Claros. “Highly efficient pre-processing of NGS reads and identification of full-length genes” - JBI2012, Barcelona. Enero 2012. 6 D. Guerrero, A. Bocinos, R. Bautista, J. Falgueras, M.G. Claros. “SeqTrimNext: preprocessing for NGS”. RES Scientific Seminar of Supercomputing and Next Generation Sequencing, M´alaga. Marzo 2011. D. Guerrero, “Infraestructuras del DataCenter” - IIR Datacenter design. Madrid. Noviembre 2011. D. Guerrero, “Introducci´on al uso de las herramientas bioinform´aticas de la PAB” - I Curso PAB de an´alisis de micromatrices - M´alaga. 2010. Otras colaboraciones y comunicaciones a congresos P. Seoane, R. Carmona, R. Bautista, D. Guerrero-Fern´andez, M.G. Claros. AutoFlow: an easy way to build workflows. International Work-Conference on Bioinformatics and Biomedical Engineering IWBBIO14. Granada, abril 2014. Pedro Seoane, Rosario Carmona, Roc´ıo Bautista, Dar´ıo Guerrero-Fern´andez y M.G. Claros. (( Using Autoflow, a workframe to resolve workflows, to build a de novo plant transcriptome)) . Plant Genomics Congress, Londres (UK), 12-13 de mayo de 2014 Hicham Benzekri, Dar´ıo Guerrero-Fern´andez, Roc´ıo Bautista, and M.G. Claros. “Detecting and correcting mis-assembled reads in contigs”. International Work-Conference on Bioinformatics and Biomedical Engineering IWBBIO13. Granada, marzo 2013 Hicham Benzekri, Roc´ıo Bautista, Dar´ıo Guerrero-Fern´andez, No´e Fern´andez-Pozo, M. G Claros. (( A reliable pipeline for a transcriptome reference in non-model species)) . International Conference The next NGS Challenge: data processing and integration. Valencia, mayo 2013 H. Benzekri, N. Fern´andez-Pozo, D. Guerrero-Fern´andez, R. Bautista, M.G. Claros “Aproximaci´on bioinform´atica al transcriptoma de Solea y su disponibilidad en SoleaDB”. Biotecnolog´ıa y recursos gen´omicos aplicados a la acuicultura. Avances logrados en AQUAGENET. Puerto Real, mayo 2012. N. Fern´andez-Pozo, D. Guerrero-Fern´andez, R. Bautista, J. G´omez-Maldonado, C. Avila, F.M. Canovas, M.G. Claros. (( GeNOTE: a web tool for annotation of non-model eukaryotic, unfinished sequences)) . Workshop on Bioinformatics for Personalized Medicine (X Jornadas de Bioinform´atica). M´alaga, 27-29/05/10 INT D. Guerrero-Fern´andez, R. Bautista, D.P. Villalobos, F.R. Cant´on, M.G. Claros. “Detection of divergent regions in aligned conserved sequences with AlignMiner”. Workshop on Bioinformatics for Personalized Medicine (X Jornadas de Bioinform´atica). M´alaga, 27-29/05/10 INT 7 ´ INDICE GENERAL ´ INDICE GENERAL VI DISCUSI ´ ON 197 14.Discusi´on final 199 14.1. Facilitar acceso a los recursos de supercomputaci´on . . . . . . . . . . . . . . . . . 199 14.2. Aprovechamiento de recursos inform´aticos . . . . . . . . . . . . . . . . . . . . . . 200 14.3. Nuevos algoritmos para biolog´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 14.4. Difusi´on y almacenamiento ordenado de los resultados de investigaci´on . . . . . . 201 VII BIBLIOGRAF´ IA 203 VIII AP´ ENDICES 223 A. Otras publicaciones 225 14 ´ INDICE GENERAL ´ INDICE GENERAL Acr´onimos y Abreviaturas Aa amino´acido ADN ´acido desoxirribonucleico AFP Apple File Protocol API interfaz de programaci´on de aplicaciones (del ingl´es Application Programming Interface) ARN ´acido ribonucleico BAM Binary Alignment Map CIFS Common Internet File System CPD Centro de Procesamiento de Datos CPU unidad central de procesamiento DNS Domain Name Server EBI European Bioinformatics Institute EC Enzyme Commission EMBL Laboratorio Europeo de Biolog´ıa Molecular EST etiquetas de secuencias expresadas (del ingl´es Expressed Sequence Tag) FDR Fourteen Data Rate FIFO First In First Out FTP protocolo de transferencia de archivos (del ingl´es File Transfer Protocol) GO Gene Ontology GPFS General Parallel File System GPU unidad de procesamiento gr´afico HPC High Performance Computing IP Internet Protocol KEGG Kyoto Encyclopedia of Genes and Genomes KVM Kernel-based Virtual Machine LDAP Lightweight Directory Access Protocol MID identificadores multiplexados MPI Message Passing Interface NCBI National Center for Biotechnology Information NFS Network File System NGS secuenciaci´on de nueva generaci´on (del ingl´es next generation sequencing) NIS Network Information Service nt nucle´otido OLC Overlap/Layout/Consensus PAB Plataforma Andaluza de Bioinform´atica PCR reacci´on en cadena de la polimerasa QDR Quadruple Data Rate RAM memoria de acceso aleatorio RDP Remote Desktop Protocol REST Representational State Transfer RoR Ruby-on-Rails SAM Sequence Alignment Map SBS Sequencing by Synthesis SCBI Centro de Supercomputaci´on y Bioinform´atica de la UMA SMB Server Message Block 15 16 ´ INDICE GENERAL SMFS Single Molecule Fluorescent Sequencing SMRT Single Molecule Real-Time sequencing SMS Single Molecule Sequencing SNP polimorfismos mononucleot´ıdicos SRA sequence Read Archive SSH Secure SHell SSR repeticiones de secuencias simples TB terabyte, unidad de almacenamiento equivalente a 1000 o 1024 gygabytes UMA Universidad de M´alaga URL Uniform Resource Locator VNC Virtual Network Computing WGA secuenciaci´on de todo el genoma (del ingl´es whole genome sequencing) X11 Interfaz gr´afica para sistemas Linux/UNIX 17 Parte I INTRODUCCI ´ ON CAP´ ITULO 1. BIOINFORM´ ATICA 19 Cap´ıtulo 1 Bioinform´atica 1.1. Visi´on general La bioinform´atica se puede considerar una de las ´areas cient´ıficas que m´as r´apido est´a avanzando en los ´ultimos a˜nos. Continuamente se producen nuevos avances tecnol´ogicos que permiten realizar experimentos m´as complejos y con un coste m´as reducido. Esto ha propiciado que la cantidad de resultados obtenidos est´e aumentando considerablemente [143]. Para que el enorme volumen de informaci´on que se est´a acumulando sea lo m´as ´util posible, hay que desarrollar repositorios complejos para estos conocimientos que es necesario mantener vivos y accesibles para futuros proyectos. Un ejemplo claro de esta tendencia podemos encontrarlo en la gr´afica de la figura 1.1, donde se observa la evoluci´on de la base de datos sequence Read Archive (SRA) [102], que almacena datos en bruto de proyectos de secuenciaci´on masiva. Es f´acil intuir que esta cantidad de datos ser´a cada vez mayor, y esto ya est´a provocando problemas de procesamiento y almacenamiento considerables, no s´olo desde el punto de vista de que la cantidad de datos generada sobrepase el l´ımite manipulable por una persona por medios manuales, sino que tambi´en sobrepasa el l´ımite de los datos que una persona puede analizar con un sistema inform´atico medio y es necesario utilizar recursos de supercomputaci´on. Recientemente se ha Figura 1.1: Evoluci´on de la cantidad de secuencias brutas almacenadas en la base de datos SRA. Fuente: NCBI-http://www.ncbi. nlm.nih.gov/Traces/sra/ acu˜nado el t´ermino BigData [8] para referirse a los problemas relativos a la manipulaci´on de ingentes cantidades de datos, que no son ´unicos de la bioinform´atica sino que se encuentran en ´areas tan diversas como la mercadotecnia, las redes sociales, y el an´alisis de mercados o tendencias [135]. Hoy en d´ıa se proponen t´ecnicas y entornos de trabajo como el MapReduce o Hadoop en grandes granjas de servidores para procesar tal cantidad de datos en tiempos razonables [107, 159], para lo que se necesita de manera indiscutible la bioinform´atica. La inform´atica adem´as puede aplicarse a la biolog´ıa para resolver problemas, sobre todo 20 CAP´ ITULO 1. BIOINFORM´ ATICA cuando la cantidad de datos disponibles hace inviable su manipulaci´on manual o tradicional, lo que abri´o el campo de la bioinform´atica. Algunas aplicaciones de la bioinform´atica son: An´alisis estad´ısticos e interpretaci´on autom´atica de resultados. Estudio de redes funcionales mediante la interacci´on entre genes. An´alisis de la estructura, funci´on e interacciones de las prote´ınas. Este fue uno de los primeros problemas que se abord´o con la bioinform´atica y sigue siendo un problema que exige grandes recursos inform´aticos. Ensamblaje de secuencias y anotaci´on. An´alisis de expresi´on. Miner´ıa de datos y de textos. Gen´omica comparativa, filogenia, estudios evolutivos y clasificaci´on de especies. Determinaci´on de mutaciones, reconocimiento de marcadores y polimorfismos mononucleot´ıdicos (SNP). 1.2. Algunas limitaciones De un modo u otro, todos los campos de la bioinform´atica presentan una serie de problemas que limitan su m´aximo aprovechamiento. Entre ellos nos encontramos con: Enorme variedad de peque˜nos programas que realizan funciones individuales. Conexi´on insatisfactoria entre dichos programa porque se necesita la traducci´on de formatos de archivos entre un programa y otro, as´ı como la manipulaci´on de los datos para adaptar su estructura. Mucho trabajo repetitivo que se puede simplificar con el desarrollo de flujos de trabajo autom´aticos para poder aplicar un mismo procedimiento a diferentes conjuntos de datos de forma eficiente. Gran cantidad de datos. En sus dos vertientes m´as problem´aticas: muchos archivos peque˜nos y archivos extremadamente grandes. Almacenamiento ordenado y difusi´on de grandes cantidades de datos a los usuarios finales de forma satisfactoria. Tiempos de ejecuci´on demasiado elevados. Poca paralelizaci´on de las aplicaciones o de los algoritmos. Altos requisitos de memoria RAM. Los algoritmos se pensaron para peque˜nas cantidades de datos, al ampliar los datos de entrada los requisitos de memoria se multiplican. Complejidad en la instalaci´on y manejo de los programas. Casi todos los programas se utilizan en un terminal mediante l´ınea de comandos lo que provoca que no sean c´omodos para la mayor´ıa de investigadores, que aunque disponen de los conocimientos necesarios para su ´area, no son expertos en inform´atica. 1.3. Ordenadores virtuales Dada la complejidad de instalaci´on de algunas herramientas bioinform´aticas y la enorme cantidad de herramientras disponibles, una de las posibilidades que han tenido los usuarios menos expertos en bioinform´atica es usar distribuciones de sistemas operativos con este tipo de programas ya instalados [60], ya sea en una CAP´ ITULO 1. BIOINFORM´ ATICA 21 m´aquina real o en forma de m´aquina virtual. De esta forma, la comunidad cient´ıfica no tiene que dedicar un esfuerzo adicional a preparar un entorno de trabajo con las herramientas bioinform´aticas que necesitan. Las primeras distribuciones las encontramos de mano de los laboratorios del CERN y Fermilab, cada uno con sus propias distribuciones Linux que usaban internamente (CERN Linux y Fermi Linux respectivamente). A principios del 2004 decidieron unir esfuerzos y lanzar Scientific Linux [52] que es una de las distribuciones m´as utilizadas en la actualidad. Tambi´en existen otras distribuciones que siguen la misma filosof´ıa, como Poseidon Linux [59]. Las alternativas m´as usadas [177] son BioLinux [60], BioPuppy [96], DNALinux [20] o LXtoo [239] que se centran en el software ´util para distintos aspectos de la bioinform´atica. Cubren la instalaci´on unificada y simple de paquetes de herramientas bioinform´aticas en los ordenadores de los usuarios y tambi´en la instalaci´on de m´aquinas virtuales que proporcionan el entorno totalmente configurado. Esta es una buena idea para facilitar el acceso de los investigadores a las herramientas bioinform´aticas necesarias, pero no ofrece el acceso a los recursos de supercomputaci´on, por lo que no se pueden usar para tareas computacionales intensivas. Ser´ıa ideal disponer de alg´un tipo de acceso de escritorio remoto que sea compatible con el uso de recursos de supercomputaci´on para las situaciones en las que los an´alisis no se pueden realizar en una m´aquina personal. 1.4. Infraestructuras centralizadas Por motivos de aprovechar mejor la financiaci´on, la organizaci´on interna, el montaje, soporte de los programas, la formaci´on y mantenimiento del sistema, en todas partes del mundo se est´an desarrollando infraestructuras centrales para apoyar los an´alisis bioinform´aticos de los grupos de investigaci´on [114] Aunque existen muchas variantes en funci´on de su evoluci´on hist´orica y el contexto cient´ıfico en el que se generaron, una de las grandes ventajas es que las personas que la componen son capaces de aprender nuevas metodolog´ıas y expandirse hacia nuevas ´areas de conocimiento de la bioinform´atica con m´as facilidad que los componentes de los grupos de investigaci´on de laboratorio, pero no deja de ser importante contratar personal que tenga conocimientos afianzados en alguna de las ´areas en las que dar´a servicio el centro, ya que eso proporcionar´a una garant´ıa de su capacidad de adaptarse a otras, o al menos aportar un punto de vista alternativo a los problemas. [115]. A´un as´ı, siguen siendo necesarios m´as bioinform´aticos para procesar todos los datos que se generan hoy en d´ıa [41], por lo que no parece que sean las infraestructuras el paso limitante, sino el personal capaz de usarlas. Quiz´a uno de los problemas sea que la carrera profesional de bioinform´atico no est´a bien definida, y se suele basar en colaboraciones multidisciplinarias que no obtienen una recompensa personal o profesional a la medida del esfuerzo invertido, pues se les considera m´as un apoyo t´ecnico que un colaborador indispensable en la investigaci´on, seg´un comenta el codirector del Bioinformatics Service Center de la Universidad de Texas (EE. UU.) [41]. De hecho, uno de los principales problemas que aparecen al ofrecer servicios de bioinform´atica de forma centralizada es que en muy pocos proyectos de investigaci´on se pueden aplicar an´alisis automatizados (por ejemplo, solo el 20 % en la Universidad de Texas), mientras que la mayor´ıa requieren una adaptaci´on de los an´alisis al problema biol´ogico concreto, o sea, necesitan investigaci´on bio- 22 CAP´ ITULO 1. BIOINFORM´ ATICA inform´atica. Por tanto, hay que tender hacia unos an´alisis bioinform´aticos en los que, aunque haya una infraestructura centralizada que facilite la parte inform´atica y bioinform´aticos que conozcan muy bien el problema biol´ogico que hay que tratar, lea los art´ıculos cient´ıficos relevantes del ´area, y est´e completamente inmerso en la investigaci´on, tambi´en se impliquen de lleno los cient´ıficos promotores de la investigaci´on y colaboren en el an´alisis de los datos y en la toma de decisiones, pues ellos son realmente los que tienen mayor experiencia en el ´area [80]. La mayor´ıa de estos servicios centrales para bioinform´atica necesitan el apoyo de recursos de supercomputaci´on. Algunos centros optan por un modelo h´ıbrido [116] en el que alojan sus propios servidores en centros de c´alculo que les proporcionan un valor a˜nadido en cuanto a mantenimiento y experiencia en administraci´on de equipos del personal, mientras que otros delegan toda la carga computacional en servicios de supercomputaci´on externos, y se limitan exclusivamente a usar los programas all´ı instalados, liber´andose totalmente de la gesti´on de los servidores y la infraestructura. En el caso de supercomputaci´on, una infraestructura centralizada implica menos problemas que en el caso de la bioinform´atica, pues los ordenadores desde el punto de vista bioinform´atico son s´olo una herramienta y no una variable de la investigaci´on como puede ser las t´ecnicas de secuenciaci´on o las aplicaciones utilizadas para los an´alisis. CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES 23 Cap´ıtulo 2 Supercomputaci´on para principiantes 2.1. Computaci´on de altas prestaciones Cuando los problemas bioinform´aticos que estamos abarcando no pueden solucionarse con los equipos inform´aticos que puede tener disponibles un laboratorio normal, llega la hora de echar mano de la supercomputaci´on. La supercomputaci´on consiste b´asicamente en coordinar la utilizaci´on de grandes cantidades de recursos de c´omputo para resolver problemas espec´ıficos, como puede ser el ensamblaje o an´alisis de datos biol´ogicos [170]. Dichos recursos pueden estar formados por grandes servidores de memoria compartida o por una agrupaci´on de gran cantidad de equipos m´as peque˜nos llamado cluster. Cada tipo de sistema presenta unas ventajas respecto al otro, y cada uno de ellos es apropiado para resolver un tipo de problemas, por lo cual un centro de supercomputaci´on que tenga la intenci´on de ser vers´atil, debe disponer de todos ellos para adaptarse a problemas de todo tipo. Aparte de los sistemas de c´alculo (unidad central de procesamiento (CPU) y memoria), un supercomputador necesita una serie de elementos no menos importantes, como pueden ser un alojamiento apropiado, almacenamiento y copias de seguridad, redes de alto rendimiento y programas que faciliten el uso y gesti´on de los recursos. El conjunto de elementos f´ısicos necesarios suele llamarse CPD [17]. En los siguientes apartados se explicar´a con m´as detalles cada uno de sus elementos. 2.2. Infraestructura del CPD Un supercomputador no puede alojarse en cualquier habitaci´on con un aire acondicionado dom´estico, es necesario contar con una infraestructura correctamente dimensionada [18] para el equipamiento que va a alojar. En la figura 2.1 podemos ver el CPD de la Universidad de M´alaga (Universidad de M´alaga (UMA)), cuyo dise˜no hemos esquematizado en la figura 2.2 y que nos servir´a de apoyo para identificar todas las partes que componen un CPD. La sala destinada a un CPD debe proporcionar un espacio amplio, con el menor n´umero posible de elementos que entorpezcan la distribuci´on de los sistemas y adem´as tener una altura suficiente para poder alojar los armarios (racks). A esa altura habr´a que descontar la necesaria para la posible instalaci´on de un suelo flotante capaz de soportar el peso de los equipos que van a instalarse encima (un ´unico ar- 30 CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES Herramientras de monitorizaci´on: Estas herramientas sirven para poder vigilar en todo momento lo que est´a pasando [9] tanto en los equipos como en el propio CPD, (por ejemplo, fallo del aire acondicionado, aver´ıas de alg´un componente de los equipos, ca´ıdas en el rendimiento, problemas en el cableado de la red de comunicaciones, etc.). Nagios [19] y Ganglia [136] son dos de las herramientas m´as utilizadas, aunque existen otras variantes como ICINGA [66], que es la utilizada en el CPD de la UMA. Aplicaciones espec´ıficas para afrontar las investigaciones que quieren realizar los usuarios finales. Las hay de muchos tipos y a menudo requiere conocimientos avanzados para poder instalarlos correctamente en un entorno distribuido como el que estamos tratando. En el CPD de la UMA se utiliza Modules [65], un sistema de carga por m´odulos que ayuda a mantener los programas organizados en diferentes categor´ıas. Cada usuario puede as´ı decidir qu´e programa y versi´on desea utilizar en funci´on de sus necesidades. Sistema de colas: Se encarga del funcionamiento correcto de un cluster de computaci´on al orquestar las ejecuciones de los diferentes programas que lanzan los usuarios y garantizar el uso correcto de todos los recursos del mismo. Entre los sistemas de colas m´as utilizados se encuentran Slurm [92], PBS, SGE [68], Torque/Maui [90], Moab [47] o Condor/HTCondor [218]. A pesar de tanta variedad, todos funcionan de un modo muy parecido. En el CPD de la UMA se utiliza Slurm como sistema de colas, de manera que el usuario accede a la m´aquina de entrada, crea un fichero (script) donde define los recursos que necesita (memoria, duraci´on del trabajo, n´umero de CPU o unidad de procesamiento gr´afico (GPU), etc.), e indica el programa que desea ejecutar junto con sus argumentos de entrada. Luego env´ıa el script al sistema de colas con el comando correspondiente y ´este organiza la ejecuci´on en base a los recursos solicitados, prioridades asignadas al usuario, otros trabajos en ejecuci´on, y otros par´ametros (m´as adelante, en en la figura 2.5 y en el apartado 2.5.3, se explicar´a esto con ficheros reales). Cuando en el cluster quedan libres suficientes recursos para comenzar la ejecuci´on, el sistema de colas ejecuta el trabajo en los servidores de c´alculo apropiados sin intervenci´on del usuario. Durante el proceso, adem´as de los resultados del programa que ha decidido lanzar el usuario, se suelen obtener otros ficheros con los resultados que el usuario hubiese visto en la pantalla si hubiese estado delante durante la ejecuci´on. Visto de forma simplificada, el sistema de colas act´ua como un delegado al que encargamos la ejecuci´on del programa y luego nos cuenta qu´e ha pasado. 2.5. Supercomputaci´on en la pr´actica El principal obst´aculo a la hora de aprovechar los recursos de un supercomputador es que se quieran ejecutar principalmente programa que solo sean capaces de aprovechar un n´ucleo de una CPU, con lo que no se har´a un uso ´optimo de los recursos. Peor a´un es que haya usuarios que crean que s´ı se puede, con lo que reservan varias CPU, de las que solo se usar´a una y las dem´as estar´an inutilizadas mientras dure la ejecuci´on, con lo que se est´an desaprovechando de forma grave los recursos. Los usuarios, que muchas veces tienen que actuar como programadores, tendr´an que saber o aprender las t´ecnicas de paralelizaci´on y creaci´on de flujos de trabajos que exponemos a continuaci´on, adem´as de utilizar correctamente el sistema de colas, si desean obtener los resul- CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES 31 tados necesarios para sus investigaciones en el menor tiempo posible. 2.5.1. Paralelizaci´on Un supercomputador no tiene porqu´e ser m´as r´apido que un ordenador convencional si no se es capaz de resolver el problema en paralelo. Muchos usuarios entran en el sistema y ejecutan en vivo alg´un programa para ver lo r´apido que es, llev´andose la sorpresa de que es s´olo un poco m´as r´apido que su ordenador personal. La verdadera potencia de un supercomputador se debe a que se pueden ejecutar problemas en paralelo para reducir proporcionalmente el tiempo usado respecto a una ejecuci´on de forma secuencial. Conseguir esto implica un esfuerzo importante de programaci´on, aunque existen otros paradigmas m´as complejos como las redes de actores [76] que se utilizan para obtener paralelismo, as´ı como lenguajes con aproximaciones puramente paralelos, como Chappel [192], Scala [158] o Swift [233] que pueden simplificar en gran medida la paralelizaci´on. Existen muchas formas de clasificar el paralelismo, pero vamos a centrarnos en las clasificaciones basadas en (1) qu´e parte del problema se paraleliza, (2) d´onde se ejecutan los trabajos y (3) c´omo se comunican entre s´ı los procesos paralelos. Respecto a la la parte del algoritmo que se paraliza, cuando la ejecuci´on paralela se obtiene porque un mismo programa, algoritmo, o conjunto de instrucciones se ejecuta de forma simult´anea sobre diferentes conjuntos de datos de entrada, hablamos de paralelismo orientado a datos. Cuando el paralelismo se obtiene al aplicar diferentes programas, algoritmos o tareas de forma simult´anea a un mismo conjunto de datos, se dice que est´a orientado a tareas [15]. La clasificaci´on basada en el lugar que se ejecutan los trabajos [164] los divide en: (a) programas paralelos, cuando su ejecuci´on no sale de un s´olo ordenador y utiliza s´olo los recursos disponibles en una m´aquina, ya sea mediante la creaci´on de hilos de ejecuci´on o de procesos adicionales, y (b) programas distribuidos, cuando un programa es capaz de utilizar varios ordenadores conectados entre s´ı. Si atendemos a la forma de comunicaci´on entre los procesos que se ejecutan en paralelo, los paradigmas m´as frecuentes son: Paralelismo compulsivo: Se trata del tipo de procesamiento paralelo m´as f´acil de conseguir y m´as eficiente, aunque eso no significa que sea el m´as habitual. Se da cuando el problema que est´as resolviendo permite partir los datos de entrada en conjuntos m´as peque˜nos y ejecutarlos por separado sin influir en el resultado final. O bien cuando diferentes programas se pueden aplicar a la vez de forma independiente al mismo conjunto de datos. De este modo, un programa de ejecuci´on lineal que cumpla estos requisitos puede ejecutarse en paralelo con el n´umero de equipos que deseemos. Los procesos (ordenadores) que trabajan en este tipo de ejecuci´on lo hacen de forma totalmente independiente, ni siquiera tienen la informaci´on de que existe otro ordenador trabajando en lo mismo que ´el. La mayor´ıa de las veces, este tipo de problemas se puede abordar sin necesidad de modificar los algoritmos implicados, desde una perspectiva de alto nivel que utiliza los propios sistemas de colas para lanzar bloques de trabajos con caracter´ısticas comunes. Todos los sistemas de colas los facilitan, ya sea de forma nativa o mediante extensiones el concepto de arrayjob, que consiste en a˜nadir de forma autom´atica un n´umero nde trabajos id´enticos al sistema de colas en los que ´unicamente var´ıa el valor de una variable de entorno. En funci´on de esa variable de entorno, cada ejecuci´on del pro- 32 CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES grama podr´a decidir el paquete de datos sobre el que desea realizar la ejecuci´on. Estos trabajos, al ejecutarse de forma independiente y posiblemente simult´anea en diferentes servidores de c´alculo del supercomputador, proporcionan el paralelismo deseado. MapReduce: Se utiliza mucho en el an´alisis de BigData, y comparte similitudes con el paralelismo compulsivo. La idea principal es b´asicamente la misma: una fase de mapeo (para preparar en paquetes los datos de entrada), la ejecuci´on parcial de cada uno de los paquetes en diferentes CPU, y finalmente una fase de reducci´on (para reorganizar o combinar los resultados parciales obtenidos de cada paquete de entrada para obtener el resultado final). Este tipo de paralelismo requiere una fase de sincronizaci´on entre los ordenadores que est´an trabajando en conjunto antes de comenzar la etapa de reducci´on que en muchas ocasiones provoca retrasos de ejecuci´on y hace disminuir el rendimiento de la paralelizaci´on. Al tratarse de un algoritmo b´asico, puede implementarse en cualquier lenguaje, aunque adem´as existen librer´ıas (frameworks) que facilitan su aplicaci´on, como son el propio MapReduce que utiliza Google [107], o Hadoop [25, 6], que es una implementaci´on de c´odigo abierto desarrollada en Java que inici´o Yahoo y que hoy d´ıa es parte del proyecto Apache. Implementaciones como Pydoop [111] proporcionan acceso a Hadoop desde Python, y comienzan a implantarse nuevas evoluciones como YARN [227] que persiguen un mayor rendimiento bajo grandes cargas realizando una gesti´on distribuida de los trabajos frente a la gesti´on centralizada que se realizaba anteriormente. MapReduce y Hadoop son excelentes candidatos para muchos an´alisis bioinform´aticos [217]. Granjas de tareas (Task-farm): Puede considerarse una variante din´amica de MapReduce en la que un proceso principal coordina un conjunto de tareas independientes, asigna su ejecuci´on a procesadores distintos (granja de procesadores) y realiza una recolecci´on de los resultados a medida que los diferentes procesadores van finalizando su ejecuci´on [50]. Las granjas de tareas suelen usar un reparto de tareas din´amico, mediante el cual se van asignando nuevas tareas a los procesadores que van quedando libres, un comportamiento que las hace especialmente ´utiles para entornos con equipos de rendimiento heterog´eneos, y evita que se queden procesadores sin utilizar porque hayan procesado su conjunto de datos m´as r´apidamente que otro. Paso de mensajes con Message Passing Interface (MPI) [220]: Hay casos en los que ejecutar un programa en paralelo/distribuido no suele ser tan sencillo como trocear los datos de entrada y aplicar un paralelismo compulsivo o MapReduce. Los problemas m´as complejos requieren una programaci´on m´as minuciosa donde se alternan fases de ejecuci´on en paralelo con otras fases de sincronizaci´on o ejecuci´on lineal. En esos casos se debe utilizar alg´un sistema de paso de mensajes entre los equipos involucrados para organizar los diferentes trabajos que est´an realizando cada uno [43]. Este tipo de ejecuci´on se beneficia del uso de una red de baja latencia, como puede ser InfiniBand o Myrinet, para disminuir al m´ınimo el tiempo en que el proceso inicial est´a parado a la espera de una respuesta de otro servidor. Cuanto m´as r´apida sea la comunicaci´on, m´as eficiente ser´a este tipo de ejecuci´on y se malgastar´a menos capacidad de c´alculo. Existen m´ultiples implementaciones de MPI, entre las que OpenMPI[71], MPICH[34] y MVAPICH[83] son las m´as utilizadas en supercomputaci´on. CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES 33 Paralelizaci´on autom´atica: Consiste en procesar autom´aticamente el c´odigo de un programa de ejecuci´on lineal y obtener cierto grado de paralelismo. Existen compiladores comerciales, como los de Intel o The Portland Group, que proporcionan la opci´on de paralelizaci´on autom´atica, aunque lo m´as extendido es utilizar herramientas como OpenMP [42], que buscan zonas en el c´odigo fuente del programa que se puedan ejecutar en paralelo (como pueden ser bucles o condicionales m´ultiples que no tengan dependencias internas) y produce una versi´on modificada del c´odigo fuente que puede ser compilada con los compiladores tradicionales para obtener una versi´on paralelizada. El programador puede ayudar a la herramienta d´andole pistas en forma de comentarios dentro del programa de qu´e partes ´el considera que son aptas para paralelizar. A veces da resultados espectaculares, y en otros casos no tanto [120], pero para bases de c´odigo que ya est´an escritas y no merece la pena reescribir desde cero de forma paralela constituye la mejor forma de obtener m´as rendimiento. Otro uso habitual es combinar OpenMP con MPI para dise˜nar un programa paralelo distribuido [175]. OpenMP se centra en programas escritos en C [98] y Fortran [13], pero la misma idea se puede encontrar en otros lenguages como R [89] de la mano de pR [119], Java [31] o Python [212] con Pydron [146]. 2.5.2. Flujos de trabajo y automatizaci´on Para realizar un trabajo de investigaci´on, lo normal es que haya que utilizar varios algoritmos para obtener el resultado buscado. Es m´as, los distintos algoritmos suelen estar en paquetes de programas diferentes que habr´a que hacer compatibles entre s´ı. Para ello, lo m´as recomendable es crear un flujo de trabajo (workflow) que automatice dicha combinaci´on Datos entrada Programa 1 Decisión/ bifurcación Programa 2 Programa 4 Programa 3 Datos salida Figura 2.4: Un flujo de trabajo consiste en encadenar ejecuciones de trabajos. y compatibilizaci´on, que de otra forma ser´ıa muy tediosa. Un flujo de trabajo permitir´a aplicar los mismos pasos a diferentes conjuntos de datos de entrada, sin tener que ejecutar manualmente cada uno de los programas [183]. De hecho, un flujo de trabajo no es m´as que un programa simplificado y resumido que hace llamadas a otros programas, ya sean instalados en la propia m´aquina o remotamente v´ıa web services, lo que permite aprovechar todas las opciones de control de flujo que se utilizan en programaci´on: puede ser totalmente lineal, incluir bifurcaciones bas´andose en los resultados previos (figura 2.4), realizar bucles de un determinado programa, realizar ejecuciones en paralelo, etc. Hay muchas opciones para crear flujos de trabajo, como por ejemplo: que el usuario haga un programa compi- 34 CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES lado que llame a otros programas instalados en la m´aquina en el orden deseado. Es una pr´actica habitual de muchos programas que realizan gran parte de sus funciones realizando llamadas a otros programas externos, como es el caso de TopHat [99] —realiza llamadas a Bowtie [108], a scripts externos de Python [171] y a Samtools [117]— para realizar un flujo de trabajo complejo con el objetivo de mapear lecturas de RNA-Seq. La ganancia de rendimiento obtenida por usar un lenguaje compilado en la creaci´on de flujos de trabajos no es significativa, puesto que el c´odigo necesario para lanzar y organizar la ejecuci´on de programas externos s´olo representa un peque˜no porcentaje del tiempo total de ejecuci´on del flujo de trabajo. usar un lenguaje de scripting como Ruby, Python, Bash, Perl, etc., para crear un programa que llame a otras aplicaciones [39]. Esta opci´on es mucho m´as flexible que usar un programa compilado, puesto que permite realizar modificaciones m´as r´apidamente y es la opci´on utilizada en muchos programas como SeqTrim [58] o GATK y Dintel [142], as´ı como en proyectos como SoleaDB [23] donde se utilizan pipelines para realizar ensamblajes e importaciones autom´aticas, o en MspJI-seq [84] para estudiar la metilaci´on de genomas; usar herramientas como Conveyor [124] o AutoFlow [198] que facilitan construir flujos de trabajo en l´ınea de comando; usar alguna herramienta que proporcione colecciones de servicios web preparados para ser enlazados, como son el MowServ [176] o BioExtract [132]; utilizar herramientas visuales, algunas tambi´en basadas en repositorios de servicios web, como Taverna [223], Galaxy [70], Tavaxy [1], Yabi [88], Pipeline Pilot [191], Ergatis [162], Kepler [128], Triana [216], Dyscovery Net [184] o YAWL [225]; a pesar de la sencillez de dise˜no y ejecuci´on que proporcionan gracias a sus elaboradas interfaces, se encuentran limitadas en cuanto a flexibilidad, o n´umero de herramientas que pueden utilizar, y algunas de ellas exigen que su ejecuci´on se realice en la misma m´aquina donde se ha dise˜nado el flujo y de una forma visual, por lo que no se integran bien con los sistemas de colas utilizados en supercomputaci´on. Para grandes cantidades de datos y uso en supercomputaci´on es mucho m´as eficiente echar mano de frameworks que de herramientas visuales porque permitan el uso de la l´ınea de comandos. Las ´ultimas versiones de Taverna y Galaxy est´an dando pasos para facilitar su ejecuci´on en sistemas de supercomputaci´on y en nuestro grupo de investigaci´on hemos desarrollado AutoFlow [198] basado en Ruby y Bash y sin interfaz gr´afica espec´ıficamente adaptado para crear, ejecutar y compartir flujos de trabajo con paralelizaci´on. Independientemente del sistema utilizado para crear el flujo de trabajo, ´este consistir´a en ejecutar un programa con los datos de entrada originales, comprobar los resultados para ver si son correctos y dichos resultados pasarlos a otro programa despu´es de haber convertido su formato si es necesario. Esta cadena de programas puede hacerse el n´umero de veces que se necesite hasta obtener el resultado final. 2.5.3. Utilidad del sistema de colas En un sistema de supercomputaci´on, el modo de trabajar difiere mucho de c´omo se trabaja en un ordenador de sobremesa. Cuando se CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES 35 utiliza un ordenador de escritorio, todo lo que ocurre aparece de una manera u otra en la pantalla. En supercomputaci´on esto no es posible puesto que el usuario final no tiene acceso al ordenador donde se est´a ejecutando su trabajo, un equipo alojado en un CPD que puede estar muy lejos del usuario y que adem´as ¡no tiene pantalla!. En supercomputaci´on se trabaja con un modelo desconectado como el resumido en la figura 2.5: el usuario accede con su nombre y contrase˜na a un servidor concreto (m´aquina de entrada) con un protocolo de comunicaci´on remoto que suele ser Secure SHell (SSH) por motivos de seguridad. En esa m´aquina el usuario preparar´a su trabajo (o un flujo de trabajo), subir´a al servidor los ficheros de entrada que quiere analizar desde su ordenador para almacenarlos en su espacio privado, y cuando lo tiene todo preparado, lo ejecuta. No tiene sentido que el trabajo se ejecute en directo porque se malgastar´ıan los recursos de supercomputaci´on, sino que el sistema de colas decidir´a cu´ando y d´onde ejecutarlo. El usuario podr´a desconectarse e incluso apagar su ordenador (que funciona solo como terminal) si quiere, porque el sistema de colas ya se encargar´a de gestionar su trabajo. Como se puede intuir, todos los programas no son aptos para ejecutarlos en un supercomputador. Solo tendr´an sentido los que tengan un modo de l´ınea de comandos que permita indicarle al programa qu´e hacer sin necesidad de utilizar ventanas gr´aficas, y tambi´en debe ser capaz de proporcionar los resultados sin necesidad de mostrarlos en una pantalla, es decir, debe escribirlos en ficheros de salida, ya sea texto, im´agenes, v´ıdeo o cualquier otro formato que el usuario pueda inspeccionar posteriormente. El sistema de colas no deja de ser un intermediario entre el usuario y los recursos de c´omputo, y se encarga de gestionar que el uso de dichos recursos sea ´optimo. Para ello requieren una puesta a punto continua por parte de los administradores de sistemas para adaptarse a las necesidades de los usuarios sin llegar a perjudicarlos ni malgastar los recursos. 2.5.4. Acceso al supercomputador de la UMA Al supercomputador de la UMA se accede mediante una conexi´on SSH al servidor picasso.scbi.uma.es con el comando ssh [email protected] desde un terminal de OS X o Linux, o con alg´un cliente de SSH de terceros como puede ser Putty [161] o WinSSH [27] si se utiliza Windows R (figura 2.6). Al acceder a su cuenta, el usuario se encuentra en el directorio de inicio alojado en el almacenamiento compartido. Este espacio dispone de una cuota de disco limitada, y podr´a usarse para almacenar los archivos que se consideren m´as importantes. Para trabajar con ficheros de datos que necesitan una cuota m´as amplia y con un acceso a disco m´as r´apido, deber´a cambiarse a su directorio de SCRATCH simplemente con cd \$SCRATCH Las cuotas de disco asignadas a los usuarios son variables y se establecen seg´un las necesidades particulares. Para intercambiar ficheros con Picasso, se debe usar un programa que soporte el protocolo SFTP, ya sea por l´ınea de comandos o alguno con interfaz gr´afica como FileZilla [61] o CyberDuck [49] (figura 2.7). 36 CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES Parte pública del sistema Parte oculta del sistema Usuario 1 - Conectar 2 - Preparar script 3 - Enviar trabajo 8 - Ejecutar 6 - Ok Cluster servidores Servidor entrada (login) Servidor sistema colas 9 - Finalizado 14 - Finalizado 7 - Desconectar 10 - Conectar 11 - ¿Finalizó ya? 12 - ¿Finalizó ya? 13 - Finalizado 16 - Desconectar Sesión 1 Sesión 2 4 - Añadir a la cola 5 - Encolado 15 - Ver o descargar resultados Figura 2.5: Modelo de trabajo desconectado de un sistema de colas La preparaci´on del trabajo consiste en crear un archivo de texto con un formato concreto para el int´erprete Bash [35], con una serie de comentarios que indican al sistema de colas Slurm qu´e recursos necesita. A continuaci´on se muestra un ejemplo concreto donde se indica que necesitaremos 10 procesadores y 2 Gb de memoria RAM durante 10 horas: #!/usr/bin/env bash CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES 37 Figura 2.6: Acceso por SSH a los ficheros almacenados en Picasso con el programa de terceros SSH Secure Shell #Number of desired cpus: #SBATCH --cpus=10 #Amount of RAM needed for this job: #SBATCH --mem=2gb #The time the job will be running: #SBATCH --time=10:00:00 #Load desired software module load blast_plus/2.2.28+ #the program to execute with its parameters: blastp -db swissprot -query prot.fasta -out test1.blast -num_threads 10 Este script se guardar´a en un fichero, por ejemplo ejemplo.sh, para poder enviarlo al sistema de colas con el comando sbatch ejemplo.sh Cuando el usuario estime que deber´ıa estar terminado o reciba una notificaci´on de que ha finalizado, entrar´a de nuevo en la m´aquina de entrada, y preguntar´a al sistema de colas qu´e ha pasado con su trabajo. En Slurm se usa el comando squeue 38 CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES Figura 2.7: Acceso remoto por SFTP a los ficheros almacenados en Picasso. De esta forma se pueden subir o bajar ficheros desde el ordenador personal al supercomputador. para un resultado similar al del siguiente listado: picasso$> squeue -a JOBID NAME USER STATE CORES ====================================== 54327 ejemplo.sh dario [R] 4 54425 gaussian.sh almu [R] 64 Si el trabajo ha terminado, encontrar´a los ficheros de resultado en su almacenamiento privado, junto con los ficheros de salida que el usuario podr´ıa haber visto por la pantalla de su ordenador, o los posibles mensajes de los errores que se hayan producido. Si tuviera que ejecutar un trabajo repetidas veces con diferentes datos de entrada, se puede echar mano de grupos de trabajos oarrayjobs. Gracias a algunas extensiones en Slurm, es tan sencillo como a˜nadir una l´ınea indicando el rango de trabajos a ejecutar y se lanzar´a un trabajo independiente por cada valor del rango. Cada trabajo recibir´a dicho valor en una variable de entorno para saber qu´e paquete de datos tiene que utilizar. En el siguiente ejemplo de c´odigo se ilustra un arrayjob en el que se solicita el env´ıo de 50 trabajos independientes que ser´an ejecutados por el sistema de colas a medida que vayan quedando huecos libres. Si existe sitio para ejecutarlos todos a la vez y no hay ninguna restricci´on en la configuraci´on del sistema de colas, todos los trabajos podr´ıan ejecutarse en paralelo: #!/usr/bin/env bash #Number of desired cpus : #SBATCH  cpus=10 #Amount of RAM needed for this job : #SBATCH mem=2gb #The time the job will be running : #SBATCH  time=10:00:00 #MAKE AN ARRAY JOB,SLURM ARRAYID will take values from 1to 50 #SARRAY  range=150 CAP´ ITULO 2. SUPERCOMPUTACI´ ON PARA PRINCIPIANTES 39 #Load desired software module load blast plus /2.2.28+ #the program to execute with its parameters : blastp db swissprot query prot${SLURM ARRAYID }.fasta out test${SLURM ARRAYID}.blast num threads 10 El fichero que define el arrayjob se env´ıa al sistema de colas con el comando sarray fichero.sh con lo que se a˜nadir´an ntrabajos al sistema de colas que aparecer´an como trabajos independientes cuando se consulte su estado con squeue: picasso$> squeue -a JOBID NAME USER STATE CORES ====================================== 54327 array[1] dario [R] 4 54328 array[2] dario [R] 4 54329 array[3] dario [R] 4 54330 array[4] dario [R] 4 54331 array[5] dario [R] 4 54332 array[6] dario [R] 4 .. ... .. ... .. ... 54376 array[50] dario [R] 4 46 CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES cada estrategia experimental, y que esta elecci´on influir´a directamente en el procesamiento bioinform´atico de los resultados. 3.5. Aplicaciones Seg´un la tecnolog´ıa y las necesidades experimentales, las aplicaciones de la ultrasecuenciaci´on se pueden dividir en tres grandes grupos. Secuenciaci´on de novo:una secuenciaci´on de novo tiene sentido cuando se pretende descubrir el genoma o transcriptoma de un organismo sobre el que no existe este tipo de informaci´on. En estos casos obtendremos mejores resultados y de forma m´as r´apida si se emplean dos t´ecnicas de secuenciaci´on diferentes, ya que cada tecnolog´ıa suplir´a sus defectos con la otra [236, 12]. Por ejemplo las secuencias obtenidas con el 454 tienen una mayor longitud, lo que facilita el ensamblaje inicial y produce un conjunto desordenado de contigs de semilla. Al complementar esta informaci´on con secuencias pareadas (ya sean de 454 o Illumina), conseguiremos determinar el orden lineal en que se sit´uan dichos contigs que conforman el boceto del genoma y, posteriormente, podemos ampliar su cobertura y rellenar los huecos. Utilizar las secuencias de un HiSeq para el ensamblaje inicial en estos casos puede ser perjudicial debido a que poseen una longitud menor que dificulta el proceso de ensamblaje cuando existen grandes zonas repetitivas en el genoma que se est´a estudiando, aunque los genomas relativamente simples puedan ser ensamblados sin mucha dificultad con secuencias cortas [166]. Resecuenciaci´on: se utiliza para secuenciar organismos del que ya disponemos de un genoma o transcriptoma fiable, con el objetivo de estudiar cambios entre individuos concretos, como SNP o microvariaciones [157, 179], realizar an´alisis de gen´omica funcional para determinar c´omo interact´uan los genes y las prote´ınas, estudiar las modificaciones evolutivas que ha sufrido una especie [73], realizaci´on de genotipado para diferenciar organismos biol´ogicos [148], o incluso estudios de gen´omica de poblaciones que proporcionan informaci´on sobre c´omo afectan las variaciones gen´eticas a grupos de poblaci´on [152]. Por norma general, la resecuenciaci´on es mucho m´as econ´omica que una secuenciaci´on de novo y proporciona un reto computacional bastante menor [36]. Si bien este tipo de secuenciaci´on podr´ıa hacerse sin problemas con cualquiera de las tecnolog´ıas descritas en la Secci´on 3.4, las secuenciaciones de HiSeq nos proporcionan mayor cantidad de datos a menor coste. Expresi´on g´enica (RNA-seq): los avances en secuenciaci´on permiten aplicarla a ´areas que antes ser´ıa impensable. Por ejemplo, antes de aparecer la ultrasecuenciaci´on, para estudiar los genes que se expresaban ante un est´ımulo concreto se utilizaban micromatrices [113, 51, 69], una t´ecnica muy delicada que permit´ıa conocer la expresi´on simult´anea de un conjunto finito de posibles genes candidatos sobre los que se ten´ıan indicios de que podr´ıan expresarse. Hoy d´ıa la ultrasecuenciaci´on permite de una forma mucho m´as inmediata y con menos manipulaci´on manual secuenciar las muestras sometidas a un est´ımulo concreto y despu´es cuantificar qu´e genes (de forma global, sin necesidad de seleccionar candidatos) se estaban expresando m´as en el momento de la secuenciaci´on en funci´on de la cantidad de lecturas que se obtienen de cada transcrito. El concepto general parte de que si no se realiza ninguna manipulaci´on de la muestra para evitarlo, la probabilidad de que las lecturas pertenecientes a cada gen aparezcan es la misma CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES 47 Segunda generación Segunda generación Segunda generación Segunda generación Segunda generación Segunda generación Tercera generación Tercera generación Fabricante Roche 454 Illumina Illumina SOLiD Life Technologies Life Technologies Pacific Biosciences Oxford Nanopore Secuenciador 454 GS FLX HiSeq 2500 HiSeq X Ten 5500 XL Ion DPM Ion Proton PacBio RS II MinION Tecnología Pirosecuenciación Secuenciación por síntesis (SBS) Secuenciación por síntesis (SBS) Ligación, codificación de 2 bases Secuenciación por síntesis en chip Secuenciación por síntesis en chip Molécula única en tiempo real - SMRT Molécula única (SMS) nanoporo Nº lecturas 1 M 4000 M 6000 M 400 M 6 M 80 M 0.05 M (50 K) 0.06 G (60 K) Longitud de las lecturas 700 - 1000 nt 125 nt 150 nt 75 nt 400 nt 200 nt 8500 nt (up to 40000 nt) 5000 nt (up to 90000 nt) Total nucleótidos 1 Gnt 1000 Gnt 1800 Gnt 30 Gnt 2 Gnt 14 Gnt 0.4 Gnt (400 Mnt) 0.2 Gnt (200 Mnt) Duración de la reacción 24 horas 6 días 3 días 2 - 7 días 5 horas 5 horas Minutos Desde minutos a horas Figura 3.3: Resumen de las t´ecnicas de secuenciaci´on de segunda y tercera generaci´on (resumido de [125, 174, 194, 126]). en todos los casos. Por ese motivo, si hacemos una secuenciaci´on y las secuencias pertenecientes a un gen aparecen relativamente m´as veces, es porque ese gen est´a m´as expresado. Respecto a la t´ecnica de secuenciaci´on a utilizar, nos encontramos en un caso similar al caso de la resecuenciaci´on, donde lo importante es la cantidad de secuencias, por lo que a mayor cantidad de datos, m´as fiable ser´a el resultado. Tiene por tanto sentido utilizar t´ecnicas que produzcan gran cantidad de secuencias peque˜nas como pueden ser HiSeq o SOLiDTM. 3.6. An´alisis bioinform´aticos Una vez obtenidas las lecturas, todo el an´alisis posterior recae en la bioinform´atica. Con respecto a las secuencias, lo podemos dividir en las siguientes grandes etapas: preprocesamiento, ensamblaje y anotaci´on (o mapeo), y an´alisis estad´ıstico o funcional. A continuaci´on se describen con m´as detalle. 3.6.1. Preprocesamiento Todas las lecturas de ultrasecuenciaci´on necesitan un preprocesamiento por m´as que los fabricantes digan que lo que proporcionan ya es adecuado para su uso. Gracias al preprocesamiento, se incrementar´a la calidad de los resultados, lo que facilitar´a el posterior tratamiento de los mismos. En la etapa de preparaci´on de muestras se a˜nadieron diferentes elementos artificiales (figura 3.2) que habr´an de desaparecer de la parte ´util de la secuencia para no proporcionar resultados indeseables o artefactuales. En esta etapa tambi´en hay que eliminar cualquier secuencia incluida de forma no intencionada, como pueden ser contaminantes de otros organismos que no nos interesan para el experimento [58]. 3.6.2. Ensamblaje El ensamblaje consiste en obtener cadenas de ADN/ARN relativamente grandes (contigs) que idealmente estar´ıan ordenadas (gracias a las secuencias pareadas) en forma de sca↵olds, a partir de los peque˜nos trozos de secuencias que hemos obtenido de uno o varios experimentos de secuenciaci´on (figura 3.4; [149]). Existen varios tipos de algoritmos para realizar este proceso, pero b´asicamente consisten en la comparaci´on de todas las secuencias obtenidas en el experimento y la creaci´on de diferentes 48 CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES CATGCATGCTAGCTGCGATGACTGATG CATGCATGCTAGCTGCAGATGACTGATG AGATGACTGATGCATGCATGCTAG GACTGATGACGTGAAT CGTGAATCATGCATG Ensamblaje CATGCATGCTAGCTGCAGATGACTGATG CATGCATGCTAGCTGCAGATGACTGATGCATGCATGCTAGCTGCAGATGACTGATGACGTGAATCATGCATG Lecturas obtenidas del secuenciador Reordenación de contigs CONTIG 2 CONTIG 1 ... Contig1: CATGCATGCTAGCTGCAGATGACTGATGCATGCATGCTAGCTGCAGATGACTGATGACGTGAATCATGCATG Formación de contigs Se busca las zonas de coincidencia entre las propias secuencias. Se obtienen diferentes contigs (agrupaciones de secuencias) inconexos entre sí. Los contigs se pueden reordenar usando experimentos adicionales o secuencias pareadas CGTGAATCATGCATGCGTGAACTGACGTACACTACGTATCATGCATG Contig2: CGTGAATCATGCATGCGTGAACTGACGTACACTACGTATCATGCATG ... Figura 3.4: Distintas etapas del ensamblaje, desde las lecturas originales (arriba) a los contigs ordenados en sca↵olds (abajo). grafos o tablas de relaciones ponderadas para determinar qu´e secuencias est´an solapadas con otras y en qu´e medida lo est´an. Siguiendo las relaciones de solapamiento y algunos heur´ısticos para acelerar las decisiones, se consigue formar una cadena mayor. Como es de suponer, los algoritmos m´as antiguos no est´an preparados para utilizar m´ultiples CPU y necesitan mantener todas las secuencias en memoria durante la fase de c´alculo de los solapamientos. Esto implica que deben utilizarse m´aquinas con mucha memoria RAM (incluso del orden de 1 o 2 terabyte, unidad de almacenamiento equivalente a 1000 o 1024 gygabytes (TB)) para realizar ensamblajes de novo complejos. El resultado de este paso ser´a un conjunto de contigs que se podr´a usar para otros an´alisis comparativos, e incluso para realizar ensamblajes incrementales o mapeos de experimentos de resecuenciaci´on posteriores. El problema del ensambaje est´a lejos de solucionarse, al menos hasta que la tecnolog´ıa de secuenciaci´on SMRT baje del 13 % de error que presenta actualmente [173]. Adem´as, no todos los organismos se ensamblan con la misma facilidad [45] ni todos los ensambladores funcionan igual en todos los organismos [55, 32]. CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES 49 Para ensamblar dos lecturas hay que detectar las regiones en las que solapan, que deben tener una longitud m´ınima para ser fiables, siempre teniendo en cuenta el tama˜no de las lecturas que tenemos disponibles. Si son lecturas de 50 nt de longitud no ser´ıa razonable exigir un solapamiento de 40 nt, mientras que ese solapamiento s´ı es perfectamente v´alido para ensamblar lecturas largas. Por la forma de detectar estas regiones solapantes, podemos encontrar diferentes tipos de algoritmos [140]: Algoritmos voraces En los algoritmos de ensamblaje voraces (figura 3.5) se parte de una lectura inicial y despu´es se selecciona el mejor candidato entre el resto de lecturas para extender la cadena. El mejor candidato ser´a aquella lectura que solape m´as nucle´otidos con la lectura inicial. El proceso se repite de forma iterativa tomando como secuencia inicial el resultado de la iteraci´on anterior, hasta que ya no queden m´as candidatos posibles para extender la cadena. Aunque no se construya una estructura de grafos propiamente dicha, s´ı que se est´a recorriendo un grafo virtual al realizar en cada paso la b´usqueda del mejor candidato. Este sistema de ensamblaje es el m´as primitivo de todos y presenta claros problemas con las zonas repetitivas, que quedar´ıan reducidas a una sola repetici´on. Ejemplos de algoritmos que encajan en esta categor´ıa son SSAKE [229], VCAKE [91] o SHARCGS [54]. Algoritmos Overlap/Layout/Consensus (OLC) La aproximaci´on OLC (figura 3.6) se realiza en tres etapas (Overlap/Layout/Consensus): a) En la primera etapa se buscan las regiones solapantes entre secuencias, para lo cual se realiza una comparaci´on de cada secuencia contra todas las dem´as. Esta comparaci´on se realiza en funcion de los porcentajes de identidad y del tama˜no m´ınimo del solapamiento junto con un par´ametro muy importante para los ensamblajes denominado tama˜no de k-mero, que normalmente elige el usuario, pero otras veces se encuentra fijado en el c´odigo fuente del algoritmo. El tama˜no de k-mero lo podemos considerar como la unidad b´asica de comparaci´on que se utilizar´a en la b´usqueda de solapamientos [44]. Al utilizar un k-mero de longitud n, se calculan todas las posibles combinaciones de nucle´otidos de dicha longitud, y despu´es se procesan todas las lecturas determinando los k-meros que contienen en com´un entre ellas. Esta abstracci´on acelera los procesos de ensamblaje y permite un ahorro considerable de memoria. b) En la segunda etapa se procesa el grafo de las relaciones establecidas entre las lecturas para establecer la forma m´as probable del ensamblaje atendiendo al n´umero de k-meros que comparten los posibles solapamientos. c) Finalmente se alinean las secuencias que solapan entre s´ı para elegir las secuencias consenso que ser´an los contigs resultantes del ensamblaje. Algunos ensambladores que utilizan esta aproximaci´on son CAP3 [85], PCAP [86], Celera Assembler/CABOG [139], Arachne [21] o Newbler [106]. Algoritmos basados en grafos de De Bruijn Un grafo de De Bruijn se utiliza para representar solapamientos entre cadenas de s´ımbolos, lo cual encaja perfectamente con el problema del ensamblaje [46]. Los ensambladores basados en grafos de De Bruijn suelen utilizarse para lecturas cortas. En la figura 3.7 se representa de forma simplificada c´omo se a˜nade cada lectura al grafo, primero se divide en oligomeros (k-meros) de un tama˜no k elegido por el usuario, y de cada k-mero se extraen los extremos izquierdo y derecho de tama˜no k1, extremos que se almacenan como nuevos nodos en el grafo de De Bruijn unidos entre s´ı con un arco dirigido desde el extremo izquierdo al derecho. Si los nodos ya existen en 50 CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES 3Solapar y repetir el proceso 2Elegir candidato que mejor solapa 1Lectura inicial ACTGTGACTGACACTACGACGACGTACGTTACG... GTACGTTACGACTGTGACTGACACTACGACTACGTTACG... GTACGACGTACGTCACTGACTGACGTAGCCGTACGCAGTCAG GTAACGTACGTACGTACGTACGTCACACGTACTGACGTA ... ACTGTGACTGACACTACGACGACGTACGTTACG.ACTGTGACTGACACTACGACTACGTTACG..... Algoritmo de ensamblaje voraz Figura 3.5: Procedimiento de ensamblaje con un algoritmo voraz 1Detección de solapamientos 3Recorrer grafo y obtener consenso 2Crear grafo relaciones ...ACGTACGTTACGACTGTGACTGACACTACGACTACGTTACG.... Algoritmo de ensamblaje OLC L1: TACGTTACGACTGTG L2: GTTACGACTGTGACTG L3: CTGTGACTGACACTAC L1: TACGTTACGACTGTG L2: GTTACGACTGTGACTG L3: CTGTGACTGACACTAC L4: ACTACGACTACGTTACG L4: ACTACGACTACGTTACG L1: TACGTTACGACTGTG L2: GTTACGACTGTGACTG L1: TACGTTACGACTGTG L3: CTGTGACTGACACTAC L1: TACGTTACGACTGTG L4: ACTACGACTACGTTACG Lecturas desordenadas Comparación de lecturas ... GTTACGACTGTG CTGTG CTGTGACTG L1 L2 L3 L4 ACTAC Los nodos representan las lecturas, los arcos las zonas que solapan Figura 3.6: Procedimiento de ensamblaje con un algoritmo OLC CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES 51 2Crear grafo con los k-mers obtenidos 1Partir lecturas en k-meros 3Recorrido del grafo Algoritmo de ensamblaje con grafos de De Bruijn Para k=6 ACTGA CTGAC TGACT GACTG ACTGAC De cada k-mero, se añaden dos nodos al grafo, uno del extremo izquierdo y otro del derecho, ambos de longitud k-1. Junto a un arco que los una. ACTGA CTGAC Se repite el proceso mientras queden k-meros, añadiendo los nodos que falten, pero siempre creando las arcos que los unen. k-mero k-1 izquierda k-1 derecha ACTGA CTGAC TGACT GACTG Se comienza a recorrer el grafo por el nodo que tiene una flecha más de salida que de entrada. Terminará en el nodo que tiene una flecha más de entrada que de salida. ACTGA CTGAC TGACT GACTG 1.Inicio 234 5 678 9.Fin La cadena resultado se formará tomando la primera letra de cada nodo por el que vamos pasando, excepto el nodo final que se añade completamente. L1: ACTGACTGACTG K1: ACTGAC K2: CTGACT K3: TGACTG K4: GACTGA K5: ACTGAC K6: CTGACT K7: TGACTG Figura 3.7: Procedimiento de ensamblaje con grafos de De Bruijn el grafo, ´unicamente se a˜nade el arco que los une. La secuencia ensamblada se obtiene recorriendo el grafo desde el extremo inicial hasta el final, concatenando la primera base del kmero almacenado en cada uno de los nodos recorridos, excepto el k-mero del nodo final que se une por completo al resultado. Al formar el grafo mediante la adici´on de las lecturas de forma iterativa, se evita hacer comparaciones de todas las secuencias entre s´ı. Adem´as es habitual que las secuencias no se almacenen en el grafo y que las repeticiones sean anotadas con un contador num´erico en vez de utilizar arcos adicionales, estrategias que ayudan a ahorrar memoria durante la ejecuci´on del algoritmo [122]. Ensambladores que usan esta t´ecnica son Velvet [240], ALLPATHS2 [133], Euler-SR [40], ABySS [203], SOAPdenovo [130] o PASHA [127]. 3.6.3. Mapeo Realizar un mapeo es mucho m´as simple que un ensamblaje (figura 3.8; [185]). En este caso se dispone de un genoma o transcriptoma de referencia ya ensamblado que ayuda a realizar el alineamiento. As´ı se consigue saber a qu´e zona de ese genoma corresponde cada una de las secuencias que hemos obtenido del experimento. Para ello se compara cada secuencia con el genoma de referencia para as´ı determinar a qu´e zona del genoma de referencia se parece m´as (de forma ideal ser´a una coincidencia exacta, aunque pueden existir peque˜nas variaciones, una veces debidas a errores de secuenciaci´on y otras a polimorfismos). Puesto que la comparaci´on es uno a uno, s´olo es necesario mantener un conjunto m´ınimo de secuencias en memoria (la referencia y la secuencia activa), y el consumo de recursos por lo tanto es mucho menor. En cuanto al uso de procesadores, este tipo de problemas es f´acilmente paralelizable ya que las comprobaciones se pueden ha- 52 CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES cer varias a la vez sin que el resultado de una comparaci´on tenga influencia alguna sobre las dem´as. Cada tecnolog´ıa de secuenciaci´on dispone de sus propios programas de mapeo comerciales, como son Eland de Illumina, Corona de SOLiDTM, o Reference Mapper de 454/Roche, pero existe una gran cantidad de programas de c´odigo abierto que suponen alternativas viables. Igual que ocurr´ıa con los algoritmos para ensamblaje, existen diferentes aproximaciones en funci´on de la forma en que se realizan los alineamientos de las lecturas con los genomas o transcriptomas de referencia: Indexaci´on de las lecturas Todas las lecturas se indizan en un hash contra el cual se alinea el genoma o transcriptoma de referencia. Ejemplos de programas de mapeo que encajan en esta categor´ıa son MAQ [118], RMAP [205] o SHRiMP [186]. Indexaci´on del genoma Se crea el hash al indizar el genoma o transcriptoma usado como referencia, justo al rev´es que en la primera aproximaci´on, y son las lecturas las que se alinean de forma iterativa contra el hash del genoma. PASS [38] y GASSST [180] son algoritmos que siguen esta estrategia. B´usqueda inversa Utiliza la transformada de Burrows–Wheeler [189] para encontrar todas las coincidencias exactas entre las lecturas y la referencia, que posteriormente se complementa con un algoritmo de b´usqueda inversa (backtracking) para encontrar coincidencias, aunque exista una peque˜na cantidad de errores en las lecturas. Este tipo de algoritmos engloba herramientas como Bowtie [108] [109] y SOAPaligner [121]. H´ıbridos Algoritmos como GenomeMapper [196] o Stampy [129] combinan diferentes estrategias de b´usqueda junto a modelos estad´ısticos para mejorar la precisi´on y velocidad de los algoritmos tracionales. GenomeMapper incluso permite el mapeo simult´aneo con varias referencias, por lo que el tiempo de mapeo se reduce considerablemente. Algunos algoritmos, como SOAP3-dp [131], incorporan el uso de GPU para acelerar los resultados de b´usquedas. En esta etapa podemos obtener diferentes resultados: un archivo de contigs oscaffolds ya ensamblados, o datos de mapeo que pueden servir para cuantificar expresi´on o determinar variaciones (SNP o mutaciones) de la muestra procesada respecto al genoma de referencia. Los resultados de mapeo tradicionalmente se proporcionaban en un fichero de texto tabulado en un formato llamado Sequence Alignment Map (SAM) [190], pero a causa del incremento en el volumen de datos producido por las t´ecnicas de ultrasecuenciaci´on se hizo necesario el desarrollo de su equivalente binario Binary Alignment Map (BAM) que ocupa menos espacio. Para manipular ficheros en ambos formatos existen las librer´ıas SAMtools [117] y BAMtools [16]. 3.6.4. Anotaci´on y an´alisis Los ensamblajes, sean de genomas o transcriptomas, necesitan una anotaci´on [208], que consiste en comparar los contigs obtenidos con otras secuencias, principalmente de bases de datos p´ublicas, para establecer si el parecido entre la secuencia nueva y la conocida es suficiente para permitirnos asignarle las mismas funciones que tiene la secuencia conocida [187]. Las anotaciones tambi´en son bastante costosas computacionalmente hablando, ya que las bases de datos con las que se suelen comparar CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES 53 CATGCATGCTAGCTGCAGATGACTGATG Mapeo/Alineamiento CATGCATGCTAGCTGCAGATGACTGATG Lecturas obtenidas del secuenciador Cuantificación del mapeo ...ACACGCCATGCATGCATGCTAGCTGCAGATGACTGATGCATGCATGCTAGCTGCGATGACTGATGACGTGAATCATGC - GENOMA DE REFERENCIA Se busca la zona de coincidencia entre cada una de las secuencias y el genoma de referencia, no las coincidencias entre las propias secuencias como ocurre en el ensamblaje. Las secuencias se van apilando, así podemos ampliar la cobertura, buscar mutaciones, o incluso cuantificar la expresión relativa de los genes dependiendo del número de secuencias que se han apilado en su dominio. GEN A GEN B ... CATGCATGCTAGCTGCAGATGACTGATG CTGCGATGACTGATGACGTGAATC ...ACACGCCATGCATGCATGCTAGCTGCAGATGACTGATGCATGCATGCTAGCTGCGATGACTGATGACGTGAATCATGC - GENOMA DE REFERENCIA CATGCATGCTAGCTGCAGATGACTG ATGCATGCTAGCTGCAGATGACTGATGCA ATGCATGCATGCTAGCTGCAGATGACTGA CGATGACTGATGACGTGAATC CGATGACTGATGACGTGAATC Figura 3.8: Mapeo de lecturas sobre una secuencia de referencia son bastante grandes y dependiendo del grado de exactitud que se le exija a la comparaci´on, podemos obtener m´as o menos resultados. Podemos agrupar el tipo de anotaciones en dos grandes grupos (figura 3.9). Por un lado, las anotaciones m´as elementales que se realizan con programas tipo Blast+ [37], hmmer [93] o exonerate [204], que realizan comparaciones de secuencias contra bases de datos conocidas [231] como son GenBank [22] o RefSeq [168] de NCBI, o UniProtKB [10, 14]. Estas anotaciones dan una idea inicial en lenguaje natural a los investigadores sobre qu´e contienen sus secuencias, pero el lenguaje humano es muy dif´ıcil de interpretar por parte de los programas de an´alisis (muchas de estas anotaciones est´an escritas en lenguaje natural y cada persona las escribe de la forma que le parece m´as correcta). Por otro lado, se pueden realizar anotaciones que siguen una ontolog´ıa determinada [193] como pueden ser Gene Ontology (GO) [75] o InterPro [145], o determinar la ruta metab´olica para obtener los mapas Kyoto Encyclopedia of Genes and Genomes (KEGG) [95] y c´odigos Enzyme Commission (EC) de la actividad enzim´atica asociada. Con este tipo de anotaciones tenemos la ventaja de que cada anota- 54 CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES a) Lenguaje natural Diferencias entre anotaciones en lenguaje natural y codificadas >sequence1 - DNA replication and chromosome cycle and related with the L-ortinithine transmembrane transporter activity >sequence2 - Activity in L-ornithine transmembrane transporter, chromosome cycle and DNA replication ... b) Codificada o usando Ontología >sequence1 - GO:0000064;GO:0000067 >sequence2 - GO:0000067;GO:0000064 ... ... GO:0000064 => L-ornithine transmembrane transporter activity GO:0000067 => DNA replication and chromosome cycle ... identificador secuencia Anotación en lenguaje natural, no tiene porqué estar escrita de la misma forma Los códigos de las anotaciones tienen una descripción textual única: identificador secuencia Las anotaciones codificadas aunque estén desordenadas son equivalentes. Figura 3.9: Representaci´on de las diferencias entre una anotaci´on textual y anotaciones que siguen una codificaci´on u ontolog´ıa. ci´on realizada tiene un c´odigo ´unico en vez de una descripci´on de texto, y por lo tanto son ´utiles para posteriormente realizar un an´alisis inform´atico de los resultados. Los programas habituales para anotar con texto natural y con c´odigos son Maker [79], Sma3s [144], Blast2GO [48] y Autofact [104]. De esta etapa idealmente obtendr´ıamos un listado de contigs con una serie de anotaciones que nos indicar´an para qu´e sirven, incluso la prote´ına que codifican, si hemos conseguido secuenciar el gen completo o s´olo una parte, etc. Lo m´as normal es que se anoten una serie de contigs y otros queden sin anotaciones, por lo que si ´estos se consideran de inter´es ser´ıa necesario realizar estudios y experimentos adicionales para determinar en qu´e est´an implicados [23]. Con los contigs ya anotados se puede empezar a interpretar y hacer an´alisis de los resultados, ya sean an´alisis estad´ısticos o b´usquedas en bases de datos especializadas (David, Cytoscape, String, Ingenuity Pathways, Gene Investigator, etc.) que pueden a˜nadir informaci´on m´as detallada, mostrar las rutas metab´olicas en las que est´an implicados dichos contigs (los genes representados por ellos), o incluso aportarte bibliograf´ıa relacionada con los mismos. 3.7. Difusi´on de resultados Una labor muy importante dentro de la investigaci´on (y que tambi´en est´a presente en las solicitudes de financiaci´on de proyectos de investigaci´on) es la presentaci´on y difusi´on de los resultados finales, as´ı como de todos los datos que se han obtenido durante el estudio. Para presentar las conclusiones y resultados cient´ıficos se utilizan las v´ıas habituales: art´ıculos en revistas especializadas, aportaciones a congresos y presentaciones de tesis o proyectos. La difusi´on de los datos obtenidos, por ejemplo los datos de secuenciaci´on suele realizarse por medio de organismos especializados (NCBI, Laboratorio Europeo de Biolog´ıa Molecular (EMBL), Japan Data Bank, GenBank, etc...), que centralizan todos los datos sobre genes y prote´ınas y los ponen disponibles para la comunidad cient´ıfica, de modo que otros investigadores puedan apoyarse en tu trabajo para avanzar. CAP´ ITULO 3. ULTRASECUENCIACI´ ON PARA PRINCIPIANTES 55 Tambi´en existen repositorios de datos brutos de secuenciaci´on o de experimentos con micromatrices, como son SRA [110], ArrayExpress [165] y GeneOmnibus [56], para que el investigador los descargue y manipule seg´un su criterio para el mismo u otros objetivos que el que llev´o a su creaci´on. No es menos importante que revistas como Nucleic Acids Research yPlant Cell Physiology dediquen un n´umero anual a las bases de datos biol´ogicas, o que Oxford University Press haya dado salida a una revista dedicada exclusivamente a bases de datos, Database. Como parte de esta necesidad de difusi´on, cada vez son m´as los proyectos en los que desean realizar bases de datos con una interfaz sencilla para exponer sus resultados y promocionar su investigaci´on lo m´aximo posible. Estos proyectos pueden beneficiarse de herramientas que han emergido en los ´ultimos a˜nos para facilitar la creaci´on de aplicaciones web que presenten los datos almacenados en bases de datos. Entre dichas herramientas nos encontramos con Ruby on Rails [74], que es una combinaci´on del lenguaje de programaci´on Ruby [62] junto a unas librer´ıas de programaci´on (Rails) que siguen un paradigma de programaci´on llamado ”convenio mejor que configuraci´on” (Convention over configuration [232]) y una arquitectura de programaci´on basada en modelo/vista/controlador [53] que ayudan en la organizaci´on de c´odigo y proporcionan un aumento de la productividad del programador. A partir de los cimientos plantados por Ruby on Rails, han aparecido librer´ıas que siguen los mismos paradigmas para otros lenguajes de programaci´on, y ya podemos encontrar alternativas maduras como Django [78] para Python [171], o CakePHP [26] para PHP [112]. Para ahorrarse desarrollos completos, nos encontramos con herramientas como GBrowse [209], que permiten publicar una base de datos gen´omicos con anotaciones gracias a una interfaz autoadaptable. 62 CAP´ ITULO 5. AUTOSERVICIO DE BIOINFORM´ ATICA ser desesperantemente lentas debido a que X11 transfiere todos los datos relacionados con la aplicaci´on para su ejecuci´on remota, aunque la ejecuci´on suele ser fluida, pues la comunicaci´on con el servidor una vez terminada la carga inicial es m´ınima. El usuario final se conecta por ssh a una m´aquina linux remota, ejecuta la aplicaci´on gr´afica con un comando en el terminal y en ese momento empieza la descarga de la mayor´ıa de datos necesarios para ejecutar el programa en el ordenador cliente y realizar la ejecuci´on directamente all´ı. El ordenador cliente debe tener instalado un entorno X11 para poder ejecutar la aplicaci´on, que suele venir instalado en Linux/UNIX, pero no en el resto de sistemas operativos. Esta soluci´on no se adapta a lo que necesitamos al no permitir el acceso a servidores con sistemas Windows R . Virtual Network Computing (VNC): Puede instalarse para todos los sistemas operativos actuales, incluso m´oviles o tablets [214]. Consta de dos partes: un servidor que se instala en la m´aquina que se desea controlar remotamente, y un cliente que se instala en la m´aquina que se va a usar para controlar la m´aquina remota. En el cliente, el usuario abre un programa, escribe la direcci´on Internet Protocol (IP) o nombre del servidor al que quiere conectarse, se identifica con su usuario y clave de VNC que le garantiza el acceso al equipo, y entonces se le abrir´a una ventana con la pantalla completa del equipo remoto. Esta soluci´on da una experiencia de usuario mejorada respecto al X11 en lo relativo a la velocidad de arranque inicial y a la interacci´on con la m´aquina remota, puesto que el usuario ve el escritorio completo, tal como lo estar´ıa viendo si estuviese sentado delante del monitor del ordenador remoto. El funcionamiento de VNC se basa en el env´ıo del bitmap completo de la pantalla a trav´es de la red (y no los datos), por lo cual solo parecer´a lento en las redes de poca velocidad. Para mejorar este aspecto, actualmente existen implementaciones (como la utilizada por Apple Inc. en su OS X), que disponen de algoritmos de calidad adaptativa y que adem´as detectan qu´e zonas de la pantalla es necesario redibujar para que se env´ıen ´unicamente los trozos de pantallas que deben refrescarse. Esta opci´on ser´ıa perfectamente v´alida para nuestro cometido, pero su integraci´on con dominios Active Directory R de Microsoft R para permitir que varios usuarios se conecten a la vez a un mismo equipo no es totalmente transparente. Exigir´ıa adem´as la instalaci´on de un programa cliente en todas las m´aquinas de usuario. Remote Desktop Protocol (RDP): Implementado por Microsoft R , el Remote Desktop Protocol [222] viene instalado de serie en los sistemas operativos Windows actuales R , por lo que se evitar´ıa la instalaci´on de un programa cliente en la mayor´ıa de las m´aquinas. Los equipos con OS X o Linux s´ı que necesitar´ıan la instalaci´on de un peque˜no programa gratuito. El RDP se integra correctamente con Active Directory [155], lo que permite usar las mismas cuentas de usuario que est´en dadas de alta en el supercomputador para acceder a la propia m´aquina. Siempre que sea posible, RDP no enviar´a una imagen de la pantalla completa a trav´es de la red, sino una primitiva que indica al programa cliente qu´e hacer, como crear una ventana de tales dimensiones, moverla, cerrarla, poner un bot´on o texto en tal sitio o una barra de desplazamiento en el otro. Esto hace que su uso en redes lentas sea muy eficaz. Adem´as, para mejorar el rendimiento, incluye opciones para limitar el n´umero de colores a mostrar, ignorar adornos de ventanas y fotos de fondo de escritorio del equipo remoto (no ocupa lo mismo una foto que enviar una primitiva que dice: fondo azul), e incluso pueden traspasar el canal de sonido del equipo remoto al cliente, integrar los dispositivos de almace- CAP´ ITULO 5. AUTOSERVICIO DE BIOINFORM´ ATICA 63 namiento local del cliente en el equipo servidor, o utilizar m´ultiples monitores. Esta es sin duda la opci´on elegida, aunque s´olo sea posible utilizarla para controlar remotamente m´aquinas con Windows R . Otras opciones comerciales: Existen otros programas que podr´ıan suplir perfectamente las necesidades de un protocolo de acceso remoto, pero a costa de esquemas de licencias y pagos peri´odicos basados en el n´umero de clientes conectados y en los servidores disponibles. Esto los hace menos atractivos para entornos de investigaci´on donde los recursos econ´omicos siempre son muy escasos. As´ı, la compa˜n´ıa Citrix Systems, Inc., tiene el programa Citrix Workspace Suite para que se pueda acceder de forma segura a aplicaciones, escritorios, datos y servicios desde cualquier dispositivo. Los usuarios (todos) tendr´ıan que instalarse el programa gratuito Citrix Receiver. Otra opci´on hubiera sido elegir TeamViewer [72], que aunque permite un uso espor´adico gratuito, requiere un pago peri´odico si se pretende hacer un uso regular. 5.3. Elecci´on del directorio de usuarios centralizado Cuando hay un conjunto de m´aquinas diferentes, lo deseable es permitir que cualquier usuario pueda utilizar cualquier m´aquina con las mismas credenciales, tanto porque facilita la gesti´on de los usuarios y las m´aquinas, como porque simplifica los accesos a los usuarios. Esto conlleva la existencia de un directorio de usuarios centralizado. El administrador realiza el alta del usuario en el sistema central y todos los equipos cliente identifican a sus usuarios en este servidor, por lo que no es necesario dar el alta de forma individual en cada equipo. Existen diferentes alternativas en funci´on de la plataforma que se utilice, entre las que cabe destacar: Network Information Service (NIS): Antes se le denominaba p´aginas amarillas [94] y era muy utilizado. Todav´ıa hoy d´ıa existe una gran cantidad de servidores que lo usan, a pesar de que se trata de una organizaci´on plana, y la tendencia general actualmente es instalar LDAP que proporciona una configuraci´on mucho m´as flexible. Se puede utilizar para identificar usuarios en redes de equipos Linux/UNIX. Active Directory R :Se trata de la soluci´on comercial de Microsoft R [155] para los servidores con Windows NT o Windows Server, y puede utilizarse para controlar la identificaci´on de los usuarios y m´aquinas cliente que usen un sistema operativo Windows R que se encuentren asociadas al dominio del servidor. Se trata, sin embargo, de un sistema con licencia variable respecto al n´umero de m´aquinas o usuarios y que s´olo permite identificar m´aquinas con Windows R instalado. LDAP: Organiza el directorio de usuarios y sus datos en un ´arbol de entradas jerarquizadas reconfigurable [199]. Por ejemplo, en el nivel m´as alto de la jerarqu´ıa podr´ıa encontrarse un centro de investigaci´on, bajo el cual podr´ıan definirse diferentes grupos de investigaci´on, que a su vez pueden agrupar personas, m´aquinas, impresoras, servicios DNS o DHCP, o cualquier entrada que se nos ocurra definir. Tambi´en se pueden definir campos con informaci´on personalizada para cada objeto que se almacena en un directorio LDAP. A la hora de realizar la autentificaci´on, las m´aquinas cliente pueden definir en qu´e rama del ´arbol va a comenzar la b´usqueda, lo que permite establecer con facilidad diferentes niveles de acceso a di- 64 CAP´ ITULO 5. AUTOSERVICIO DE BIOINFORM´ ATICA ferentes grupos de m´aquinas. Por ejemplo las m´aquinas de un laboratorio estar´ıan configuradas para comenzar la b´usqueda a partir del nodo de ese laboratorio, por lo que no tendr´ıan acceso a los datos de un laboratorio distinto. LDAP puede integrarse con un servidor Samba para ofrecer servicios de identificaci´on de usuarios simulando una m´aquina Active Directory, pero adem´as lo puede utilizar el resto de m´aquinas Linux/UNIX u OS X para identificar sus usuarios o incluso ofrecer servicios de list´ın telef´onico. La elecci´on por tanto para el sistema de directorio centralizado de usuarios ser´a LDAP (en concreto su implementaci´on OpenLDAP) y un controlador de dominio Samba. 5.4. Elecci´on del almacenamiento compartido Acceder a diferentes m´aquinas con el mismo usuario sin disponer de un sistema de almacenamiento compartido es algo que no tiene mucho sentido porque los archivos se quedar´ıan repartidos por las diferentes m´aquinas que el usuario hubiese utilizado. Para ello existen protocolos que permiten que distintas m´aquinas compartan los mismos archivos. A su vez, cada usuario debe tener su espacio de almacenamiento privado en un servicio de almacenamiento compartido. Los servicios de almacenamiento compartido m´as habituales para los entornos de escritorio son los siguientes: NFS: Se trata del protocolo m´as utilizado en entornos Linux/UNIX y consiste en un servidor NFS al que se conectan los diferentes usuarios desde sus m´aquinas cliente mediante un cliente NFS que suele venir instalado en el sistema [202]. Apple File Protocol (AFP): Se trata del protocolo que los equipos con sistemas operativos de Apple utilizan para compartir archivos entre ellos. Desde la llegada de OS X, ´estos tambi´en pueden acceder a NFS y Samba sin ning´un problema. Server Message Block (SMB)/Common Internet File System (CIFS): El protocolo Server Message Block (SMB) lo utilizan las redes de Windows R para compartir archivos. Tambi´en es el protocolo utilizado en los dominios Active Directory R para proporcionar a los usuarios el acceso a sus datos y directorios compartidos [153]. Samba: Adem´as de actuar como un controlador de dominio equivalente a Active Directory R , proporciona los mismos servicios que SMB/CIFS para entornos Linux/UNIX y OS X [172]. Como podemos ver en la figura 5.1, el servidor permite definir los directorios y permisos que cada usuario tendr´a disponible, pero adem´as, al integrarlo con LDAP, nos permitir´a que un usuario obtenga autom´aticamente un punto de montaje con sus datos cuando se conecta a cualquiera de las m´aquinas. En realidad, los datos est´an alojados en el mismo sistema de almacenamiento centralizado usado en la parte de supercomputaci´on, situaci´on que permite que el usuario tenga a su disposici´on en una m´aquina con entorno visual y herramientas gr´aficas los mismos datos que en el entorno de sistema de colas y l´ınea de comandos del supercomputador. Esta es por tanto la opci´on que hemos elegido para nuestro sistema de almacenamiento compartido. CAP´ ITULO 5. AUTOSERVICIO DE BIOINFORM´ ATICA 65 Almacenamiento centralizado x 60 x 10 Servidor Samba Servidor LDAP Espacio usuario Espacio usuario Espacio usuario Identificación Acceso a disco compartido Figura 5.1: Almacenamiento compartido para clientes basado en Samba 5.5. Elecci´on del sistema de virtualizaci´on La virtualizaci´on [167] consiste en dividir un ordenador real, con su memoria, almacenamiento, procesadores y red, en un conjunto de m´aquinas ficticias o virtuales, e instalar en cada una de esas subm´aquinas otro sistema operativo (que se dice que est´a virtualizado) que pueden ser iguales o distintos entre s´ı, con sus aplicaciones y programas correspondientes. Podemos por lo tanto considerar esta t´ecnica como un modo de consolidar recursos que permite reducir los costes de adquisici´on de equipos y de gesti´on de su funcionamiento a cambio de un impacto m´ınimo en el rendimiento de la m´aquina virtualizada. Para usar virtualizaci´on necesitamos usar un hipervisor, que es un programa que se sit´ua debajo del sistema operativo virtualizado y que se encarga de proporcionarle los distintos servicios que ofrecer´ıa una m´aquina real: acceso a disco, acceso a memoria, ejecuci´on de c´odigo, acceso a la red, etc. El hipervisor ser´a el encargado de interceptar las instrucciones que ejecute la m´aquina virtual y realizar las adaptaciones necesarias antes de pasarlas al hardware subyacente. La intervenci´on del hipervisor es tan liviana que el rendimiento de una m´aquina virtual es muy cercano al de una m´aquina real [87]. El buen rendimiento tambi´en se debe al uso de procesadores y chipsets modernos en los que los fabricantes incorporan extensiones para mejorar aun m´as el rendimiento de las m´aquinas virtualizadas [4], como son las tecnolog´ıas VT-X [151] y VT-d [2] de Intel R o AMD-V [5] del fabricante AMD R . En funci´on de la capa en la que se sit´ue el hipervisor entre el sistema virtualizado y el hardware real distinguimos dos tipos principales [11]: Hipervisor de tipo 1 o baremetal:el hipervisor se instala directamente sobre el hardware de la m´aquina (figura 5.2), sin necesidad de tener un sistema operativo anfitri´on que le d´e soporte. El hipervisor captura las peticiones del sistema virtualizado y las pasa directamente a los elementos correspondientes de su hardware, realizando las tareas de planificaci´on y encolado necesarias. Los baremetals son los m´as eficientes al no requerir un sistema operativo intermedio. Como ejemplos tenemos VMware ESX o VMware ESXi [154], XEN [169], MS Hiper-V [97] o Kernel-based Virtual Machine (KVM)[101]. Hay quien considera que KVM es un hipervisor tipo 2 (v´ease el p´arrafo siguiente), pero en realidad es tipo 1 ya que 66 CAP´ ITULO 5. AUTOSERVICIO DE BIOINFORM´ ATICA Hipervisor de tipo 2Hipervisor de tipo 1 SERVIDOR (HARDWARE) SERVIDOR (HARDWARE) MV 1 MV NMV 3 HIPERVISOR (Software) SISTEMA OPERATIVO ANFITRION ... HIPERVISOR (Software) MV 1 MV NMV 3 ... Figura 5.2: Comparaci´on entre los hipervisores de tipo 1 y 2 se ejecuta en el mismo nivel que el kernel del sistema operativo y no encima de ´este. Este sistema es el m´as utilizado para virtualizaci´on de servidores y es el que utilizaremos en nuestro sistema porque proporciona mayor rendimiento y aprovechamiento de los recursos, lo que permite alojar m´as m´aquinas virtuales por cada servidor real utilizado. Hipervisor de tipo 2 o hosted:En este caso, el hipervisor se instala sobre un sistema operativo (figura 5.2) ya existente (ya sea Linux/UNIX, OS X, Windows,... ). El hipervisor pasa las peticiones del sistema operativo virtualizado hacia el sistema operativo real que las pasa a su vez al hardware de la m´aquina. Este tipo de sistema tiene un rendimiendo algo inferior al tipo 1 debido a que hay recursos que ya se est´an utilizando para el sistema operativo anfitri´on y adem´as existe una capa adicional que puede ralentizar las peticiones. Por otro lado, tienen la ventaja de que el hipervisor puede convivir con los dem´as programas de usuario que est´en en dicha m´aquina y esto proporciona mayor comodidad para virtualizaci´on en equipos de escritorio. Ejemplos son VMware Workstation, VMware Player, VMware Fusion o VMware Server, QEMU o VirtualBox. ¿Virtualizaci´on completa o paravirtualizaci´on? El sistema operativo virtualizado puede ser un sistema operativo sin modificaciones, en cuyo caso hablamos de virtualizaci´on completa, o puede ser un sistema operativo especialmente modificado para comportarse mejor en entornos virtualizados, en cuyo caso hablamos de paravirtualizaci´on. Existen tambi´en drivers paravirtualizados para dispositivos espec´ıficos (como la tarjeta de red) que ayudan a incrementar el rendimiento cuando se utiliza virtualizaci´on completa. En el sistema que vamos a desarrollar en CAP´ ITULO 5. AUTOSERVICIO DE BIOINFORM´ ATICA 67 este trabajo usaremos la virtualizaci´on completa con hipervisores de tipo 1, en concreto con el sistema ESX de VMware instalado en un cluster de dos servidores HP DL380-G6, con 8 n´ucleos y 32GB de RAM cada uno, que proporcionan las ayudas de Intel R para la virtualizaci´on. 5.6. Broker de m´aquinas virtuales Servidores virtualización ESX Servidor web Usuario Gestor de reservas Servidor del broker Máquina virtualN Máquina virtual2 Máquina virtual1 Cliente del broker Etiquetado de máquinas Figura 5.3: Arquitectura del broker de m´aquinas virtuales Una vez definido que utilizaremos m´aquinas virtuales con un sistema operativo Windows R y un dominio Samba con autentificaci´on LDAP para el acceso de los usuarios y los datos, surge el problema principal: son m´aquinas virtuales a las que se va a acceder con un protocolo de acceso remoto (RDP), lo que requiere la constante vigilancia de qui´en est´a utilizando cada m´aquina, qu´e m´aquinas est´an libres y cu´ales est´an en uso. Para que este problema no repercuta en la facilidad del uso del sistema ni en la comodidad de los usuarios, hemos desarrollado lo que llamamos un broker de m´aquinas virtuales con una interfaz web que tiene la siguiente arquitectura (figura 5.3): Servidor del broker:Se trata de un servidor escrito en perl que recibe peri´odicamente el estado de cada una de las m´aquinas virtuales. El estado incluye informaci´on sobre el posible usuario que est´a registrado en la m´aquina, y sirve para mantener una base de datos con los usuarios activos y m´aquinas ocupadas. Las m´aquinas tendr´an etiquetas para determinar su tipo, pues no todas tendr´an las mismas caracter´ısticas ni las mismas aplicaciones instaladas. Cliente del broker:Se trata de un programa en Delphi que se ejecuta como un servicio en cada una de las m´aquinas Windows R ,yse conecta con el servidor del broker para facilitarle los datos del usuario que est´a usando la m´aquina en ese momento. Si no hay usuario, el broker entender´a que la m´aquina est´a libre. Servidor web: Muestra a trav´es de la web un listado de los tipos de m´aquinas disponibles correspondientes a cada etiqueta, bas´andose en la informaci´on que mantiene el servidor del broker. Desde esta web podr´a seleccionar el tipo de m´aquina deseado, lo que desencadena la ejecuci´on del gestor de reservas (apartado 5.6), y se descarga un archivo con extensi´on .rdp con 68 CAP´ ITULO 5. AUTOSERVICIO DE BIOINFORM´ ATICA los datos de acceso necesarios para conectarse a una m´aquina virtual libre. Este fichero se abrir´a autom´aticamente (en Windows R ) o con un doble clic en el programa cliente del protocolo RDP y le pedir´a las credenciales de acceso para abrir la m´aquina. En la figura 5.4 se muestra el proceso de conexi´on a al escritorio remoto. El usuario podr´a elegir entre: (1) una m´aquina B´ ASICA que contiene todas las aplicaciones que tenemos en la Plataforma Andaluza de Bioinform´atica para gen´omica, transcript´omica, expresi´on g´enica, etc., que se puede consultar en la p´agina web http://www.scbi. uma.es/site/scbi/software; (2) una m´aquina GEDECYDER, que consiste en una m´aquina b´asica m´as la aplicaci´on GE DeCyder-2D para an´alisis diferencial de geles de prote´ınas; y (3) la m´aquina DISCOVERY que est´a tambi´en basada en la m´aquina b´asica y contiene los programas Discovery Studio, HyperChem y PDQuest. Gestor de reservas: Se trata de un programa escrito en Perl que hace de intermediario entre la p´agina web que presenta los datos al usuario y el servidor que recibe la informaci´on de estado de las m´aquinas. Este programa se comunica con el servidor del broker para obtener el n´umero de m´aquinas disponibles de cada etiqueta y, posteriormente, cuando el usuario realiza la petici´on de una m´aquina con una etiqueta concreta, se encarga de solicitar que esa m´aquina sea retirada del conjunto de m´aquinas disponibles y reservada para que el usuario pueda conectarse. El control proporcionado por el broker permite adem´as monitorizar cuantas m´aquinas se est´an utilizando en cada momento. Adem´as, se le pueden a˜nadir scripts al sistema para arrancar o apagar m´aquinas adicionales cuando sea necesario. El broker dar´a de baja autom´aticamente las m´aquinas apagadas o incorPortal web Selección de tipo máquina Conexión remota por RDP Figura 5.4: Etapas de la conexi´on a los escritorios remotos a trav´es del portal web. porar´a las m´aquinas reci´en arrancadas a la reserva de m´aquinas disponibles. CAP´ ITULO 5. AUTOSERVICIO DE BIOINFORM´ ATICA 69 Portal web Usuario Servidores virtualización ESX Máquina virtualN Máquina virtual2 Máquina virtual1 Cliente broker Broker Almacenamiento centralizado x 60 x 10 1Reservar máquina Escritorio remoto RDP 2Dar conexión Samba + LDAP 3 - Conectar Figura 5.5: Esquema de funcionamiento del broker de m´aquinas virtuales 5.7. Conclusiones Hemos creado un sistema de escritorios virtuales basado en el uso de hipervisores de tipo 1 (secci´on 5.5) con acceso remoto mediante protocolo RDP (secci´on 5.2) en combinaci´on con un sistema de autentificaci´on LDAP (secci´on 5.3) y almacenamiento compartido basado en Samba (secci´on 5.4) que permite solicitar los recursos a trav´es de la web gracias al broker descrito en la secci´on 5.6, y proporciona a los usuarios la comodidad necesaria para acceder a programas instalados, ya sea desde sus despachos, laboratorios o incluso desde su propia casa. En la figura 5.5 podemos ver c´omo se conectan entre s´ı los diferentes elementos que hemos ido introduciendo en las secciones anteriores: el usuario inicialmente reserva una m´aquina a trav´es de un portal web, el portal consulta al broker e indica al usuario los datos de conexi´on a la m´aquina. Cuando el usuario se conecta mediante su programa de escritorio remoto (RDP) a la m´aquina virtual indicada por el broker, se encuentra con que su espacio reservado en el almacenamiento centralizado ya ha sido montado a trav´es de Samba y por tanto tiene sus datos disponibles en la unidad Z: de la m´aquina (figura 5.6). Para intercambiar archivos deber´a usar un programa de SFTP conect´andose al almacenamiento centralizado del supercomputador, como explicamos en el apartado 2.5.4. Adem´as, al disponer de m´aquinas con distintos conjuntos de aplicaciones, los propios usuarios decidir´an qu´e tipo de m´aquina y herramientas usar´an para realizar su trabajo. Como ya hemos indicado, los datos de usuario que se utilizan en todas las m´aquinas est´an almacenados en un almacenamiento compartido, con lo que no tienen que estar copiando o moviendo datos de una m´aquina a otra. 70 CAP´ ITULO 5. AUTOSERVICIO DE BIOINFORM´ ATICA 1 - Acceso a archivos compartidos por SFTP 2Acceso a los archivos desde el escritorio remoto Figura 5.6: Acceso a los archivos por SFTP y a trav´es de la conexi´on a escritorio remoto de la m´aquina virtual. Este trabajo planta la semilla del “Autoservicio de supercomputaci´on para bioinform´atica”, que seguir´a evolucionando con el resto de trabajos aportados en esta tesis. El sistema implementado actualmente permite el acceso a tres tipos de m´aquinas virtuales distintas, aunque este n´umero puede configurarse simplemente a˜nadiendo una nueva etiqueta para la IP de la m´aquina en el servidor del broker. CAP´ ITULO 6. GENERADOR DE INTERFACES WEB 71 Cap´ıtulo 6 Generador de interfaces web para dar acceso a programas desarrollados en l´ınea de comandos Mi contribuci´on: El desarrollo completo de la herramienta InGeBIOL. , { tag : ”ANNOTATION TRANSCRIPTOMICS” , title : ”Transcriptomics” } ] }, #second category { tag : ”ASSEMBLY”, title : ”Assembly tools” } ] A more complex hierarchy with three levels of categories is exemplified in Appendix A.4.4. Every command defined in InGeBIOL is expected to be tagged with at least one category. Therefore, the command will appear on the branches corresponding to those tags. If no category (tag) is defined for a command, it is then not shown in the dock but it is still accessible by using its direct URL such as https://your_ingebiol_site/sessions/new/ Command_folder_name (to know about command folder names, see subsection 3.3). 3.2.4 Data storage configuration The file globals.rb within the folder “globals” (Figure 2) requires modifications to indicate where jobs (DATA PATH) and private repositories (PRIVATE DATA PATH) are going to be stored. For example, the installation in the PAB environment requires that jobs and the private repository are saved in the SCRATCH disk; thus, the corresponding environment variables are modified as follows: #file :globals/globals .rb #Path where jobs are going to be saved DATA PATH = /mnt/scratch/ingebiol data #Where to save private repository files PRIVATE DATA PATH = /mnt/scratch/ ingebiol priv 3.2.5 Queue system The job sent from the web interface of InGeBIOL is directed to the computer CPUs. If InGeBIOL is installed in a location with access to a supercomputer, as is the case in PAB, the job must use the corresponding queue system syntax. Therefore, it must be specified in the QSUB CMD environment variable declared in queue system.rb file within the “globals” folder (Figure 2). The minimum customisation requires to uncomment the corresponding environment variable (PBS and Slurm), while other queue systems require defining the command. In the PAB, the queue system Slurm is used, therefore, the modification will be the following: #file :globals/queue system .rb #command used to send job to slurm QSUB CMD = /usr/bin/sbatch #command used to send job to PBS #QSUB CMD = / usr/pbs/bin/qsub #Command used to launch local jobs LOCAL CMD = bash Additionally, the shell used to execute local jobs must be declared in the environment variable LOCAL CMD; by default, the most common shell for linux, bash, is declared. If any other shell is used, it can be specified in this variable. A complete description of parameters within this file is given in Appendix A.4.3. 3.3 Defining new commands Each command handled by InGeBIOL gathers all configuration files within a folder inside the “commands” folder (Figure 2). As a quick example, we are going to describe the changes necessary to define a simple web interface for GENote .1 [11]. The most important files are common.json, to control the name and appearance of the command on the presentation page (Figure 5A); one or more numbered stageXX.json files to define the input files and parameters that will be requested to the user in the submit form (Figure 5B) as well as the command line tools to be executed with submitted data; and a final results.json file, where it is defined how to show on the web page the results obtained after the execution (Figure 6). As a result, adding a new tool to a running InGeBIOL installation is very easy and fast, as it only requires duplication of an already existing folder (for example the “demo” folder (Figure 2) corresponding to a generic command, see Appendix subsection A.3) and the modification of the desired parameters taking advantage that JSON files are human-readable. 3.3.1 Customising the command presentation page The command line GENote .1 [11] is defined by the files within the “genote” folder within the “commands” folder in Figure 2. The common.json configuration file serves to define the presentation page, 6 Navigation dock Login panel Command name and description REST API definition Navigation dock Job submission formJobs queue Text input widget } } } } } B A Job results File selector input widget Numeric input widgets Fig. 5: Interface of the GENote .1 command implemented in InGeBIOL: A) Login and presentation page, B) Job submission form shown after user login, composed of di↵erent input widgets. 7 File output widgets . } } } } } . Image output widgets Fig. 6: Display of results for GENote .1 using output widgets 8 as shown in Figure 5A, including name (title), descriptions (short description), logotypes, etc., as well as the category or categories (category tags) of the navigation dock into which this command will be included. Here it is the important part of the file whose modification customises the web page: { #Write here the desired title for the command title : ”GENote v.&beta ; ” , #Provide ashort description short description : ”A WEB INTERFACE FOR GENote”, #Command version version : ”v &beta ;0.1”, #Show command in the dock in dock : true , #Under these categories : category tags : [ ANNOTATION GENOMICS ] } The complete set of parameters that can be defined in this file is detailed in Appendix subsection A.3. 3.3.2 Input data form A command can be presented as a single-page HTML form (such as in Figure 5B) or as multiple forms shown in consecutive pages called ‘stages’. Every stage is defined in its corresponding stageXX.json file (where XX are numbers used to sort the stage files and determine their presentation order). Each stage is comprised of one or more input widgets to capture input data from the user and one or more commands with the calls to the command line tools that will be sent for execution to the queue system. A definition in JSON code requires that a name (title) is given to the stage, the list of widgets is declared (input params), the calls to command line tools is defined (command list), and the use of a queue system and necessary resources are configured (submit file header). InGeBIOL can also submit jobs for execution to the local machine, but this is not a good practice since it could collapse the server and should only be used in exceptional situations (only with very fast commands). A generic example is as follows: #stage template { #title for this stage title : SUBMIT FILES , #list of input widgets input params : [ INPUT WIDGET 1, INPUT WIDGET 2, ... INPUT WIDGET N ], #commands to be executed with input widgets values command list :[ COMMAND LINE CALL 1, COMMAND LINE CALL 2, ... COMMAND LINE CALL M ], #send this command list to queue system use queue system : true , #request 6cpus and 10gb for 10 hours to queue system submit file header :[ #SBATCH  cpus=6 , #SBATCH mem=10gb , #SBATCH  time=10:00:00 ] } When users gain access to this command via web, they are required to fill at least the compulsory fields (in red) and then press “Send” button. If the error validation implemented in each widget encounters any issue, the user is notified to fix it. When no validation error is found, the command line calls, input data and parameters are converted to a queue system script and sent for execution. If two or more stages are defined for the command, they will appear successively in lexicographical order. Defining input widgets: We call ‘input widget’ to a small piece of JSON that defines the widget properties using “parameter:value” tuples in this way: #INPUT WIDGET 1definition { ”input type”:text , ”id”:field1 , ”title”:Text field } And renders this widget on screen: These three properties in the listing are mandatory for all widgets and their meanings are the following: •input type: is the field type, selected from text,integer,float,file,file popup, checkbox,radio,popup or separator. 9 •id: a unique identifier to distinguish between all input parameters of a stage. This identifier will be used to place each input widget value in the command line calls, by replacing any widget id found in the command line call definition with the corresponding value. •title: this text string is shown on screen to inform the user about the widget content. Widgets may also contain the following optional properties that present a default behaviour if they are not specified: •default value: The default value for the input widget that is shown in the submit form so that the user can modify it before submission. The use of default values can speedup user interaction with the submission form, and give them a hint about what is a valid value for this input widget. •required: Widgets can be mandatory or optional. When required is true, the field content is checked by the validation engine, and the submission process is aborted if empty. If it is false or not defined, the widget is optional and can have an empty value and no error will be reported. •size: It serves to indicate the width of the field in the form. If no size is specified, then the default value stablished by the HTML specification is used. •tooltip: When specified, a small info icon is placed near the field, and the tooltip text is shown when the mouse is moved over it. •command switch: It defines how the value of this widget will replace its respective widget id placeholder when placed in the definition of a command line call. By default it will take the value %s, that indicates to replace the whole value of the widget without modifications. A value of ‘‘-option %s’’ will be replaced by ‘‘-option widget value’’, this way any command option or flag can be directly generated by the widget. There are a few optional parameters with the purpose of input data validation in the case of text, integer or float input widgets: •validation regex: A Ruby regular expression that will be used to validate input data on text input widgets. If no value is specified, then no validation will be reported. •validation lower limit and validation upper limit: For integer or float input widgets, these two parameters define the accepted range of input values. If no value is set for the lower limit, then no restriction is imposed for a minimum value. If no value is set for the upper limit, then no maximum value is forced. •validation error msg: The error message that will be reported to the user if the validation fails. A complete list of input widgets already defined in InGeBIOL can be found at Appendix subsection A.1. Of course, other new types of widgets can be defined. The interface defined for GENote .1 (Figure 5B) uses one stage with five input widgets: one text widget for the job name, two integer widgets for minimum identity and hit length, a float widget for minimum E-value and a file popup widget to present a file selector to receive the input file. Executing command-line tools: Usually, only one command line will be executed after a stage, but it could be the case that some selection requires the execution of several equivalent or di↵erent tools. Therefore, a list of tools to be executed is defined with a JSON snippet that contains the list of required files (the execution will be halted if any file listed on required files is missing) and the command necessary to launch the tool or application. In the listing below, the definition of the unique command line call necessary to launch GENote .1 is shown, indicating in required files that input file.fasta is required for execution: #COMMAND LINE CALL 1definition { required files :[ input file .fasta ], command : . ˜ genote/init env ;genote .rb input file .fasta evalue min ident min hit length 6>output .txt } The defined command line call will be parsed and references to input widgets ( evalue,min ident and min hit length) will be replaced by their values or command switch expressions if present. For example, if input widgets for GENote .1weredefined like this: { input type : integer , id : min ident , title : Minimum identity (%) , command switch : i %s , 10 ... } , { input type : float , id : evalue , title : Minimum evalue (eg .: 1e 6) , ... } , { input type : integer , id : min hit length , title : Minimum hit length (nt), command switch : L %s , ... } If the user provides the value 90 for min ident, 1e-6 for evalue and 15 for min hit length in their respective input widgets, the executed command line will be: .˜genote/initenv ; genote.rb input file . fasta 1e6i90L156>output . txt Submitting to a queue system: Execution of jobs in any queue system requires the definition of several parameters, such as number of CPUs, memory and execution time. Although they could be introduced by the stage form, we decided to fix these values as a submit file header variable on the stage configuration file. In the use case of GENote .1, we are requesting six cores and 10 GB of RAM with an estimated execution time of 10 hours: submit file header :[ #SBATCH  cpus=6 , #SBATCH mem=10gb , #SBATCH  time=10:00:00 ] 3.3.3 Output definition: Alongside stage files, there is only one final results.json file, and it is used by InGeBIOL to define the output of the command and render the results page (Figure 6) that is shown when the user clicks on the job name in the jobs list. This file contains a list with the necessary output widgets selected from the predefined ones (Appendix A.2), or new, customly-defined ones. A generic example of final results.json can be: { #title of the final results title : RESULTS , #define alist of output widgets params : [ OUTPUT WIDGET 1, OUTPUT WIDGET 2 ]} where a title is given, and the params variable contains a list of the required output widgets. An output widget is defined as follows: #OUTPUT WIDGET 1defining an output file { type : file , file : output files/annotations .txt , title : Annotations , tooltip : } where the items have the same meaning than in the input widgets, and is rendered on screen with this appearance: The disk icon near the file name allows the user to download the file (the same as clicking on the name), while the bookmark icon will add the file to the user file repository. The execution of GENote .1 provides images and text files Figure 6. Therefore, the output stage contains two image ouput widgets that directly show the images on the web page, and five file ouput widgets that show the name of each file. 3.3.4 REST web service Since web services are now widely used as a way to support interoperable machine-to-machine interactions and REST (Representational State Transfer) has gained widespread acceptance across the web and has been adopted by mainstream web 2.0 service providers [4], InGeBIOL automatically adds REST compatibility [3] to every configured command. Therefore, all commands accessible via web can also be integrated with other tools that use REST web services to launch jobs, query job status, download result files, or any of the predefined actions. REST web service protocol, directly based on HTTP requests, is faster and easier to use than the overloaded SOAP based ones [23], but InGeBIOL makes it even easier because, for each command, users can find a specific help page with the parameters and calling examples for a simple REST client like the curl command line application (Figure 7). 11 Fig. 7: Generated REST API documentation for the customised command of GENote .1 12 3.4 Interactive user interface Web interfaces used to be static, but the arrival of AJAX has pave the way to produce interactive web interfaces that can resemble in some way the menubased software. Here it is explained how to exploit interactivity in InGeBIOL simply adding some lines to the input widget definition in the stage file. Self-validation: As stated previously (section 3.3.2), InGeBIOL comes with automatic input data validation out of the box. This self-validation can be based on regular expressions (text), or value-range definition (integer or float). Its main advantage is to avoid errors in the execution pipeline due to user mistakes. Defining a range can also give users a clue about valid values for the input widget because the valid range is shown on the interface. For example, if we define a range of values for a integer to be checked adding the following lines in the widged definition: #Input widget with range validation { type : integer , ... validation lower limit : 20 , validation upper limit : } the user will see this limits in the interface: Regular expression validation can force the user to specify a controlled string. For example, if a field is expected to contain only nucleotides, it can be stated as follows: #Input widget with RE validation { type : text , ... validation regex : [ACTGactg], validation error msg : You must provide a valid nucleotide sequence (ACTG) } Show and hide parameters: When the interface contains many optional parameters, having all of them in view can be annoying. To prevent this, InGeBIOL can group input widgets into multiple optional groups that will be shown grouped to the user with a button to dynamically toggle their visibility as in the following image: This can be done by adding a separator input widget, with a title to be shown on screen and a folding name that will serve to tag other widgets. The default state of this separator is declared on folding closed. #Define afolding separator initially closed { input type : separator , title : Advanced parameters , folding name : advanced , folding closed : true } As a result, the title with the toggle button is shown on screen (as in the previous image). To control the widgets that will toggle their presence with this separator, the widget must include the field folding name in its definition: #input widget with folding name { input type : text , id : species , folding name : advanced , title : Species name , ... } Additionally, individual input widgets can be hidden/shown (without being tagged with any folding name) based on any other input widget value (for example a popup selector). This is done using the visible if parameter into the widget definition: #input widget visible only when using the fasta option { input type : file popup , id : quality file , visible if : ”file type”==” fasta” ... } 13 Execution templates: A command definition can have a conditional execution expression that will be evaluated after replacing the widget placeholders (widget id) with their values. If the result of evaluating the expression is true, the command line call protected by it will be executed, but if the result is false, its execution will be skipped. This feature enables the definition of execution templates (sets of predefined parameters) by setting up a popup input widget with the di↵erent available templates, and then using the selected template to conditionally execute the command line application with di↵erent predefined parameters. In the listing below, a popup input widget with two assembly options (Unitigger or Best Over Graph) is defined: #popup input widget with two execution templates { input type : popup , id : assembly method , title : Assembly methods , default value : bog , values : {Unitigger :utg , Best Over Graph :bog }, ... } This code will show on screen a popup menu to choose one of the available options: where the option selected on the popup will refer to one of the possibilities of the command line call definitions. In other words, the exec if can use its value to conditionally execute the commands. The analysis of this widget requires the definition of two command line calls, one to be be executed when assembly method input widget have the value “bog” (that means that the user has selected the “Best Over Graph” method on the popup), and the other will be executed when assembly method contains the value “utg” corresponding to the selection of the “Unitigger” assembly method. #this command line call will be executed when bog is selected on the popup { exec if : ”assembly method”==”bog” required files :[ pname .frg ], command : runCAOBT .pl unitiger=bog createACE=1 pname .frg }, #this command line call will be executed when utg is selected on the popup { exec if : ”assembly method”==”utg” required files :[ pname .frg ], command : runCAOBT .pl unitiger=utg createACE=0 pname .frg } Interative help: One optional help file in HTML or text format can be used to give the user information about the usage of the command. It will be shown in an independent panel side by side to the submit form. The user can dismiss or recall it when necessary. The help file must be placed in the same folder as the stage definition file, and its name specified on the stage definition file by using the help file variable, as stated in the listing below: #stage template with help file { #title for this stage title : SUBMIT FILES , help file : help .html , #list of input widgets input params : [ INPUT WIDGET 1, INPUT WIDGET 2, ] ... } In this example, the HTML content of the help.html file that is shown in the following listing: 1.Select an input format <br/> <br/> 2.Provide a sequences file in the selected input format <br/> <br/> 3.Select a template file <br/> <br/> 4.You can provide advanced params , l i k e species or contaminants database <br/> <br/> 5.Submit your job by clicking on ”Send” button will render a help screen with five steps indicating how to fill the submit form: 14 3.5 Input and output layouts The layout of input and output widgets is automatically organised without any user intervention. However, it has been implemented an extra level of interface customisation making use of layout templates that serve to define how input or output widgets should be placed on screen. To use a custom template, a template variable pointing to the template file must be specified on the stage.json or final results.json files, depending if the template will be used for input or output widgets respectively. In the listing below: #using an RoR HTML template file in a stage configuration { title : SUBMIT FILES , stage type : submit , enabled : true , template : stage1 .html .erb , ... atemplate file stage1.html.erb is defined for input widgets to declare a table with two columns to place two widgets side by side. It must be in the “commands” folder (Figure 2), and can contain any valid HTML code. Each <%= html[’widget id’] %> snippet inside the template will be replaced by the code necessary to show the widget referenced by widget id: <table border=”0” width=”100%”> <tr> <td> <%=h t m l [ widget id 1]%> </td> <td <%=h t m l [ widget id 2]%> </td> </tr> </table> 3.6 Workflows can be implemented The current implementation of InGeBIOL is not only able to provide web interfaces to simple command-line tools and gather them as a portal: workflows and pipelines can also be given web interfaces. It is recommended to describe complex workflows in dedicated tools like Conveyor [20] or Autoflow [26], and then use their command line execution capabilities to call them from InGeBIOL using workflows as if they were a black box. We mentioned (section 3.4) that the command list is an array of command elements defined in the following way: { exec if : conditional expression , required files :[ list of required files ], command : command to be executed } where the optional exec if field is a piece of Ruby code that will be evaluated (it can contain references to values of any input widget of the stage), and the associated command line will be executed only when the evaluation returns true. This feature can be used to control the execution path in simple workflows, for example the simple workflow outlined in this schema: Input file Selected file type seqtrim.pl command Convert SFF to fasta Extract ZIP fasta file zip filesff file where the command line to “convert SFF” will only be executed if the input type is a SFF file, and the command line to “extract ZIP” will be executed if the selected input file is a zip file. The result of both alternative branches is a file of sequences in fasta format to feed the seqtrim.pl 15 4. Image: it serves to display a single result image in a determined size. It is defined by the following code: { ”id”:”image1”, ”type”:”image”, ”file name”:”image file .png”, ”content type”:”image/png”, ”title”:”Some plot graphic”, ”options”:{”height”:”200”,”width”:”250”}, ”linked”:true , ”tooltip”:”” } and produces the following output widget: 5. Image browser: it displays a set of images in a single browser panel, useful when an undetermined number of images need to be shown. It is defined by the following code: { ”id”:”stat images”, ”type”:”image browser”, ”title”:”Stats graphs”, ”file”:”graphs/ ” , ”content type”:”image/png”, ”options”:{”height”:”100”,”width”:”150”}, 22 ”max images”:0 , ”columns”:4, ”linked”:true , ”show names”:true , ”ignored images”:[”any file name .png”] } and produces the following output widget: 6. Command output: it can display the output text obtained after executing a command line. This is not designed to execute full commands here, but to show excerpts of result files using fast commands such as head,tail,wc,grep, etc.. It is defined by the following code: { ”id”:”output1”, ”type”:”command output”, ”command”:”du sh ”, ”title”:”Generated files”, ”tooltip”:”” } and produces the following output widget: A.3 Implementing a new command Here it is described a step-by-step implementation of a new command within InGeBIOL. Every command has everything within the configuration folder config/commands/ that contains a few JSON configuration files defining the input widgets, command-ine calls to be executed and output widgets. When InGeBIOL 23 is installed, it will come with a demo command that will not appear by default in the dock, but that can be accessed at http://your_ingebiol_site/ingebiol/session/new/demo). Installation of the very first command can be done simply by duplicating this “demo” folder within the “commands” folder, changing the folder name and editing the JSON files contained in it. For example, the customisation of a command that converts all text found on a file to uppercase using the tr linux command is illustrated. The tr command simply translates one character into another and prints the result on the screen. Hence, the following command: echo abcabc |tr a A will output the word AbcAbc, because it changes all the lowercase “a” chars with the uppercase version: “A”. To convert all downcase chars to uppercase, we could specify all chars, the range of chars, or use the special lowercase class denoted by [:lower:] , and its uppercase equivalent [:upper:]: echo abcabc |tr [: lower :] [: upper :] that will produce ABCABC, but since we need to convert whole files, the full command will have to receive an input file.txt and output an output file.txt taking advantage of standard input and output redirections: tr [: lower :] [: upper :] <input file .txt >output file .txt The third command line call above will be used in the example outlined in the following steps: 1. Go to “commands” folder: cd config/commands 2. Duplicate the “demo” folder changing its name to ‘uppercase’: cp -r demo uppercase 3. Go to “uppercase” folder: cd uppercase 4. Edit common.json file and change the values of the title,short description and long description variables to change the appearance of the presentation page used for the command, and set in dock to true so that it is shown on the dock (subsubsection 3.3.1): ######################################################## #Common parameters of command # #You only can modify the values of each key . # #Comments must be in their own line preceded with a”#” #symbol . # ######################################################## { #Write here the desired title for the command ”title”:”Uppercase converter”, #Provide ashort description ”short description”:”A demo interface for UPPERCASE using InGeBiol”, #And along one ”long description”:”Long description for command”, #Command version ”version”:”v 0.9a1”, #Contact url 24 ”contact url”:”http://www .scbi .uma .es”, #Corporate logo name ,file must be in public/images ”corporate logo”:”logogrande .png”, #Command logo name ,file must be in public/images ”command logo”:”ingebiolBig .png”, #Small corporate logo name ,file must be in public/images ”small corporate logo”:”logoscbi .png”, #Small command logo name ,file must be in public/images ”small command logo”:”ingebiolSmall .png”, #Copyright ”copyright”:”Copyright 2009 SCBI”, ”in dock”:true , ”category tags”:[ASSEMBLY ], #Download links shown under login box links :[ { text : Download , url : http://www .scbi .uma .es/downloads } ], #Articles shown at bottom of the page ”articles title”:If you find this command useful ,please cite :, ”articles”:[ { text : Name of the article , url : http:// url of the article } ] } 5. Edit the input params variable of the stage0001.json file to add a text input widget for the job name, and a file popup input widget to receive the input file.txt that will be converted to uppercase (subsubsection 3.3.2). Add the previously discussed tr command line call to the command list variable in the same stage file to convert the input file.txt to uppercase and generate the output file.txt (section 3.3.2): ######################################################## #Common parameters of command # #You only can modify the values of each key .Be careful #with punctuation symbols ,since this is ajson file and #needs to keep its structure . # #Comments must be in their own line preceded with a”#” #symbol . # #You can add as many stages as you want ,they are shown #in alphabetical order . # ######################################################## { ”title”:”SUBMIT FILES”, 25 ”stage type”:”submit”, ”enabled”:true , ”help file”:”help .html”, ”submit button title”:”Next”, ”input params”:[ #text input widget for job name { ”input type”:”text”, ”id”:”job name”, ”title”:”Text field”, ”default value”:”” , ”required”:true , ”size”:”30”, ”tooltip”:”Provide acustom name for the job”, ”validation regex”:”” , ”validation error msg”:”” , ”standard attribute”:true , ”save in presets”:true }, #file popup input widget for input file .txt { ”input type”:”file popup”, ”id”:”input file”, ”title”:”File with popup”, ”default value”:”” , ”command switch”:”%s ” , ”filter”:”templates/ . txt”, ”global filter”:[”” ], ”required”:true , ”save as”:”input file .txt”, ”max file size”:1, ”tooltip”:”Any hint” } ], ”command list”:[ #command line call used to convert input file .txt to uppercase and store it on output file .txt { ”required files”:[] , ”command”:”tr [:lower :] [:upper :] <input file .txt >output file .txt” } ] , ”use queue system”:false } 6. Edit final results.json adding to the params variable a file output widget to let the user download the output file.txt, and a command output widget to show on screen the first 20 lines of the same output file.txt. ######################################################## #Final results of command # #You only can modify the values of each key . # #Comments must be in their own line preceded with a”#” #symbol . # ######################################################## { ”title”:”SELECTED RESULTS”, ”stage type”:”final results”, 26 ”enabled”:true , ”template”:”” , ”params”:[ #file output widget to show adownload link to output file .txt { ”type”:”file”, ”file”:”output file .txt”, ”title”:”Output file in uppercase”, ”tooltip”:”” }, #command output to show 20 first lines of output file .txt { ”type”:”command output”, ”command”:”head n20 output file .txt”, ”title”:”Head of output file”, ”tooltip”:”” } ] } 7. It can be tested by pointing the browser to http://your_ingebiol_site/ingebiol/session/new/ demo. If there is any error in configuration files, an error message will appear instead of the web interface of the command. If configuration files are correct, a presentation page similar to the Figure 5A will be shown. There, the user can login as a guest introducing any valid e-mail, write a job name, upload a file, and send the job. It will appear on the job list section, where its name can be clicked to view the results page including a download link with the output file.txt file and the twenty first lines of it. A.4 Global configuration InGeBIOL have a set of features that a↵ect to the integration with existing environments. Global configuration files are located at /home/rails/ingebiol/config/global folder (this may change if you installed it on another location). A.4.1 File: globals.rb This is a ruby file used to define global configuration parameters that a↵ect to all interfaces. Example file with comments: #This file defines global configuration parameters #Path where jobs are going to be saved DATA PATH = F i l e . e x p a n d path( File . join (BASE PATH, ../ ,guidata )) #Path where config files are saved CONFIG PATH = F i l e . e x p a n d path( File . join (BASE PATH, ../ ,config )) #Where to save repository files PRIVATE DATA PATH = F i l e . e x p a n d path( File . join (BASE PATH, ../ ,guidata ,private )) #Path of commands configurations COMMAND CONFIG = F i l e . e x p a n d path( File . join (CONFIG PATH, commands )) #Path of globals configuration GLOBALS CONFIG = F i l e . e x p a n d path( File . join (CONFIG PATH, global )) #Path to user scripts USER SCRIPTS PATH = F i l e . expand path( File . join (CONFIG PATH, scripts )) #How are stages numbered STAGES PATTERN = stage .json 27 #Standard attributes are saved in afile with this name STANDARD ATTR JSON = std attr .json #This script is used to get additional info about ajob GET JOBINFO SCRIPT = get job info .rb #These are the titles shown in the job list in the web interface JOBLIST TITLES JSON = joblist titles .json #This tag is replaced by all command swithes of every input param #when it is found in acommand description COMMAND SWITCHES TAG = COMMAND SWITCHES #This is the default command to show when no one is defined DEFAULT COMMAND = login A.4.2 ldap.rb This is a ruby file used to define the values of the LDAP authentication server and related configuration: Example file: #This file is used to define the values of the LDAP #authentication server and related configuration : #set to true to use LDAP authentication USE LDAP AUTH = false #set to true to use LDAP authentication while allowing GUEST access USE MIXED AUTH = true #LDAP authentication server LDAP HOST = 10.200.190.1 #How to query users to LDAP server (may change between OpenDirectory and OpenLDAP servers) LDAP USERS UID= uid=%s , ou=people ,dc=mylab ,dc=uma ,dc=es A.4.3 queue system.rb This is a ruby file used to define the commands used to submit jobs to queue systems: Example file: RUNNING FILE = RUNNING QUEUED FILE = QUEUED #command used to send job to slurm QSUB CMD = /usr/bin/sbatch #command used to send job to PBS #QSUB CMD = / usr/pbs/bin/qsub QSUB CMD = /usr/pbs/bin/qsub qroutex86 #sudo command to be used with QSUB CMD if needed QSUB SUDO = /usr/bin/sudo ubioperl #Command used to launch local jobs LOCAL CMD = bash LOCAL SUDO = 28 A.4.4 categories.json #This file defines an hierarchy of categories .Then each program/interface can be tagged to belong to any of these categories . # #Programs not tagged with avalid category are found under the ”BAD CATEGORY”title . [ { ”tag”:”PREPROCESS”, ”title”:”Seq preprocessing”, ”tooltip”:”Software to preprocess raw sequences prior to assembly” } , { ”tag”:”ASSEMBLY”, ”title”:”Seq assembly”, ”tooltip”:”Software to assembly already preprocessed genomics or transcriptomics data”, ”children”:[ { ”tag”:”ASSEMBLY GENOMICS”, ”title”:”Genomics”, ”children”:[ { ”tag”:”SHORT SEQUENCES”, ”title”:”Using short sequences” }, { ”tag”:”LONG SEQUENCES”, ”title”:”Using long sequences” } ] } , { ”tag”:”ASSEMBLY TRANSCRIPTOMICS”, ”title”:”Transcriptomics” } ] } , { ”tag”:”ANNOTATION”, ”title”:”Seq annotation”, ”tooltip”:”Software to annotate genomics or transcriptomics data after assembly”, ”children”:[ { ”tag”:”ANNOTATION GENOMICS”, ”title”:”Genomics” } , { ”tag”:”ANNOTATION TRANSCRIPTOMICS”, ”title”:”Transcriptomics” } ] } ] References [1] Stephen Altschul, Barry Demchak, Richard Durbin, Robert Gentleman, Martin Krzywinski, Heng Li, Anton Nekrutenko, James Robinson, Wayne Rasband, James Taylor, and Cole Trapnell. The anatomy of successful computational biology software. Nature Biotechnology, 31:894–897, 2013. [2] Michael B¨achle and Paul Kirchberg. Ruby on rails. IEEE Software, 24:105–108, 2007. 29 [3] Thomas Bayer. REST Web Services. . . . 2002) http://www. oio. de/public/xml/rest-webservices. ..., pages 3–5, 2007. [4] Djamal Benslimane, Schahram Dustdar, and Amit Sheth. Services mashups: The new generation of web applications. IEEE Internet Computing, 12(5):13–15, 2008. [5] Daniel Blankenberg, Gregory Von Kuster, Nathaniel Coraor, Guruprasad Ananda, Ross Lazarus, Mary Mangan, Anton Nekrutenko, and James Taylor. Galaxy: A web-based genome analysis tool for experimentalists, 2010. [6] Brandi L. Cantarel, Ian Korf, Sofia M C Robb, Genis Parra, Eric Ross, Barry Moore, Carson Holt, Alejandro S´anchez Alvarado, and Mark Yandell. MAKER: An easy-to-use annotation pipeline designed for emerging model organism genomes. Genome Research, 18(1):188–196, 2008. [7] Tim Carver and Alan Bleasby. The design of Jemboss: A graphical user interface to EMBOSS. Bioinformatics, 19:1837–1843, 2003. [8] B Chevreux, T Pfisterer, B Drescher, A J Driesel, W E G M¨uller, T Wetter, and S Suhai. Using the miraEST Assembler for Reliable and Automated mRNA Transcript Assembly and SNP Detection in Sequenced ESTs. Genome Research, 14:1147–1159, 2004. [9] J Falgueras, A J Lara, N Fern´andez-Pozo, F R Cant´on, G P´erez-Trabado, and M G Claros. SeqTrim: a high-throughput pipeline for pre-processing any type of sequence read. BMC Bioinformatics, 11:38, 2010. [10] Juan Falgueras, Antonio J. Lara, Francisco R. Cant´on, Guillermo P´erez-Trabado, and M. Gonzalo Claros. SeqTrim - A validation and trimming tool for all purpose sequence reads. In Advances in Soft Computing, volume 44, pages 353–360, 2007. [11] M.Gonzalo Fern´andez-Pozo, No´e and Guerrero-Fern´andez, Dar´ıo and Bautista, Roc´ıo and G´omezMaldonado, Josefa and Avila, Concepci´on and C´anovas, FranciscoM. and Claros. GENote v.:A Web Tool Prototype for Annotation of Unfinished Sequences in Non-model Eukaryotes. In Bioinformatics for Personalized Medicine, pages 66–71. 2012. [12] David Flanagan and Yukihiro Matsumoto. The Ruby Programming Language, volume 159. 2008. [13] Don Gilbert. Pise: software for building bioinformatics webs. Briefings in bioinformatics, 3:405–409, 2002. [14] Dar´ıo Guerrero, Roc´ıo Bautista, David P Villalobos, Francisco R Cant´on, and M Gonzalo Claros. AlignMiner: a Web-based tool for detection of divergent regions in multiple sequence alignments of conserved sequences. Algorithms for molecular biology : AMB, 5:24, 2010. [15] T Ho↵mann and N Laird. fgui: A Method for Automatically Creating Graphical User Interfaces for Command-Line R . . . . Journal of Statistical Software, 2009. [16] X Huang and A Madan. CAP3: A DNA Sequence Assembly Program. Genome Research, 9:868–877, 1999. [17] JSON.org. Introducing JSON, 2014. [18] Sudhir Kumar and Joel Dudley. Bioinformatics software for biologists in the genomics era, 2007. [19] A Lara, G P´erez-Trabado, D Villalobos, S D´ıaz-Moreno, F Cant´on, and M G Claros. A Web Tool to Discover Full-Length Sequences: Full-Lengther. In E Corchado, J M Corchado, and A Abraham, editors, Innovations in Hybrid Intelligent Systems, pages 361–368. Springer, 2007. [20] Burkhard Linke, Robert Giegerich, and Alexander Goesmann. Conveyor: A workflow engine for bioinformatic analyses. Bioinformatics, 27(7):903–911, 2011. [21] J R Miller, A L Delcher, S Koren, E Venter, B P Walenz, A Brownley, J Johnson, K Li, C Mobarry, and G Sutton. Aggressive assembly of pyrosequencing reads with mates. Bioinformatics, 24(24):2818–2824, 2008. 30 [22] Brad A. Myers and Mary Beth Rosson. Survey on user interface programming. In CHI, pages 195–202, 1992. [23] Kumar Potti Pavan, Ahuja Sanjay, and Prodano↵Zornitza. Comparing Performance of Web Service Interaction Styles : SOAP vs . REST. In 2012 Proceedings of the Conference on Information Systems Applied Research, pages 1–24, 2012. [24] Paolo Romano, Ezio Bartocci, Guglielmo Bertolini, Flavio De Paoli, Domenico Marra, Giancarlo Mauri, Emanuela Merelli, and Luciano Milanesi. Biowep: a workflow enactment portal for bioinformatics applications. BMC bioinformatics, 8 Suppl 1:S19, 2007. [25] Torsten Seemann. Ten recommendations for creating usable bioinformatics command line software. GigaScience, 2:15, 2013. [26] Pedro Seoane, Rosario Carmona, Roc´ıo Bautista, and Others. AutoFlow: an easy way to build workflows. [27] Ilhami Visne, Erkan Dilaveroglu, Klemens Vierlinger, Martin Lauss, Ahmet Yildiz, Andreas Weinhaeusel, Christa Noehammer, Friedrich Leisch, and Albert Kriegner. RGG: a general GUI Framework for R scripts. BMC bioinformatics, 10:74, 2009. [28] Erik Wilde and Robert J. Glushko. XML Fever, 2008. [29] J. Sergio Zepeda and Sergio V. Chapa. From desktop applications towards ajax web applications. In 2007 4th International Conference on Electrical and Electronics Engineering, ICEEE 2007, pages 193–196, 2007. 31 110 CAP´ ITULO 7. EJECUCI´ ON DISTRIBUIDA DE TAREAS CAP´ ITULO 8. COMPRESI´ ON DE SECUENCIAS 111 Cap´ıtulo 8 Formato comprimido para almacenar y manejar las lecturas de los ultrasecuenciadores en un supercomputador Mi contribuci´on: En este trabajo desarroll´e inicialmente la librer´ıa de compresi´on en Ruby a modo de prototipo para evaluar la viabilidad del proyecto. Posteriormente entre Rafael Larrosa y yo dise˜namos el formato comprimido final y realizamos la implementaci´on en C. Una vez terminada la liber´ıa en C, realic´e las mediciones de rendimiento con los distintos tipos de secuencias y una gema scbi fqbin para facilitar su uso en Ruby. 112 CAP´ ITULO 8. COMPRESI´ ON DE SECUENCIAS CAP´ ITULO 8. COMPRESI´ ON DE SECUENCIAS 113 T´ıtulo del art´ıculo: FQbin a compatible and optimized format for storing and managing sequence data Autores: Guerrero-Fern´andez, D.; Larrosa, R. and Claros, M. G. Publicaci´on: IWBBIO 2013 - International Work-Conference on Bioinformatics and Biomedical Engineering Proceedings Volumen:2013 pp. 337-344 A˜no: 2013 DOI/ISBN:ISBN978-84-15814-13-9 Abstract Existing hardware environments may be stressed when storing and processing the enormous amount of data generated by next-generation sequencing technology. Here, we propose FQbin, a novel and versatile tool in C for compressing, storing and reading such sequencing data in a new and Fasta/FastQ-compatible format that outperforms the existing proposals. It is based on the general-purpose zLib library and o↵ers up to 10X compression. The compressed file is read and decompressed up to 3X faster than a FastQ file is read, and a nearly ´ınstant’random access to every entry in the FQbin container is allowed. Fast file reading is maintained even in shared storage environments, where dif-ferent processes are simultaneously accessing the same FQbin file. Slow networks can take even more advantage from FQbin. 114 CAP´ ITULO 8. COMPRESI´ ON DE SECUENCIAS 115 Parte IV DESARROLLO DE HERRAMIENTAS BIOINFORM ´ ATICAS CAP´ ITULO 9. DETECCI´ ON DE REGIONES DIVERGENTES 117 Cap´ıtulo 9 Herramienta web para detectar las regiones m´as divergentes en un alineamiento de secuencias Mi contribuci´on: He desarrollado el algoritmo de AlignMiner as´ı como la interfaz web interactiva en la que se incluye tanto la detecci´on de las regiones divergentes, como una herramienta para dise˜nar los cebadores ´optimos que distingan cada secuencia. 118 CAP´ ITULO 9. DETECCI´ ON DE REGIONES DIVERGENTES CAP´ ITULO 9. DETECCI´ ON DE REGIONES DIVERGENTES 119 T´ıtulo del art´ıculo: AlignMiner: a Web-based tool for detection of divergent regions in multiple sequence alignments of conserved sequences. Autores: Guerrero, D., Bautista, R., Villalobos, D. P., Cant´on, F. R., and Claros, M. G. Publicaci´on: Algorithms for Molecular Biology. Volumen:5 A˜no: 2010 DOI:http://dx.doi.org/10.1186/1748-7188-5-24 Abstract Background Multiple sequence alignments are used to study gene or protein function, phylogenetic relations, genome evolution hypotheses and even gene polymorphisms. Virtually without exception, all available tools focus on conserved segments or residues. Small divergent regions, however, are biologically important for specific quantitative polymerase chain reaction, genotyping, molecular markers and preparation of specific antibodies, and yet have received little attention. As a consequence, they must be selected empirically by the researcher. AlignMiner has been developed to fill this gap in bioinformatic analyses. Results AlignMiner is a Web-based application for detection of conserved and divergent regions in alignments of conserved sequences, focusing particularly on divergence. It accepts alignments (protein or nucleic acid) obtained using any of a variety of algorithms, which does not appear to have a significant impact on the final results. AlignMiner uses di↵erent scoring methods for assessing conserved/divergent regions, Entropy being the method that provides the highest number of regions with the greatest length, and Weighted being the most restrictive. Conserved/divergent regions can be generated either with respect to the consensus sequence or to one master sequence. The resulting data are presented in a graphical interface developed in AJAX, which provides remarkable user interaction capabilities. Users do not need to wait until execution is complete and can.even inspect their results on a di↵erent computer. Data can be downloaded onto a user disk, in standard formats. In silico and experimental proof-of-concept cases have shown that AlignMiner can be successfully used to designing specific polymerase chain reaction primers as well as potential epitopes for antibodies. Primer design is assisted by a module that deploys several oligonucleotide parameters for designing primers .on the fly”. Conclusions AlignMiner can be used to reliably detect divergent regions via several scoring methods that provide di↵erent levels of selectivity. Its predictions have been verified by experimental means. Hence, it is expected that its usage will save researchers’time and ensure an objective selection of the best-possible divergent region when closely related sequences are analysed. AlignMiner is freely available at http://www.scbi.uma.es/alignminer webcite. 0 7,5 15 22,5 30 BAC1 BAC2 BAC3 Number of contigs 0 12500 25000 37500 50000 BAC1 BAC2 BAC3 N50 length (nt) Suppl. Figure 1: Assembly of 3 different Pinus pinaster BAC clones using 454 paired-end reads and Newbler 2.3 with and without pre-processing using SeqTrimNext. A) Number of contigs obtained after assembling; clearly, SeqTrimNext pre-processed reads provide fewer contigs than the original reads preprocessed by Newbler. B) N50 value for the same BAC clones, reflecting that SeqTrimNext preprocessed reads can be assembled into longer contigs. Scaffolds of the three BAC clones after SeqTrimNext pre-processing are devoid of linkers, adaptors and E coli sequences. Only Newbler SeqTrimNext + Newbler A B CAP´ ITULO 11. AN´ ALISIS DE TRANSCRIPTOMAS 127 Cap´ıtulo 11 Herramienta para analizar transcriptomas de especies no modelo que se han ensamblado de novo Mi contribuci´on: He contribuido a paralelizar y distribuir la ejecuci´on del algoritmo en todas las etapas posibles y a resolver los problemas computacionales relacionados con el manejo de las bases de datos subyacentes. Tambi´en he desarrollado gemas de Ruby que han servido para compatibilizar las entradas y salidas del programa, que se pueden consultar en https:// rubygems.org/search?query=scbi. El art´ıculo est´a en el formato de la revista Bioinformatics porque se envi´o a ella, pero fue rechazado y lo estamos revisando para volverlo a enviar a otra revista. 128 CAP´ ITULO 11. AN´ ALISIS DE TRANSCRIPTOMAS BIOINFORMATICS Vol. 00 no. 00 2012 Pages 1–17 FULL-LENGTHERNEXT: A tool for fine-tuning de novo assembled transcriptomes of non-model organisms No´ e Fern´ andez-Pozo 1, Dar´ ıo Guerrero-Fern´ andez 2, Roc´ ıo Bautista 2and M. Gonzalo Claros 1,2⇤ 1Departamento de Biolog´ ıa Molecular y Bioqu´ ımica, Facultad de Ciencias, Universidad de M´ alaga, 29071, Spain. 2Plataforma Andaluza de Bioinform´ atica, Centro de Supercomputaci´ on y Bioinform´ atica, Edificio de Bioinnovaci´ on, Universidad de M´ alaga, 29590 M´ alaga, Spain. Received on XXXXX; revised on XXXXX; accepted on XXXXX Associate Editor: XXXXXXX ABSTRACT Motivation: De novo transcriptome assemblies devoid of any genomics or transcriptomics reference occur commonly when working with non-model species. There is no easy way to distinguish between artefacts and species-specific transcripts, or to identify the best possible assembly before deep annotation or complementary laboratory research. Results: FULL-LENGTHERNEXT has been designed as a parallelisable and distributable command-line, web-tool and REST web-service pipeline adapted to high-throughput transcriptomics. It has applications in (1) the classification of unigenes among coding for complete or incomplete proteins, (2) the construction of an open reading frame fixing the likely frame-shifts in unigene assemblies; (3) the discovery of putative species-specific unigenes; and (4) the extraction of putative non-coding RNAs. Combining the previous information, it provides a quick overview for the selection of the best de novo transcriptome assembly. Particularly, it has been shown that independent de novo assemblies using MIRA3 and Euler-SR can be reconciled using CAP3 in a more successful transcriptome. Also, unigenes receiving an Unknown status usually are short sequences that can be discarded for subsequent processings. Availability: FULL-LENGTHERNEXT can be freely downloaded or executed via web at http://www.scbi.uma.es/full lengther2 Contact: [email protected] Supplementary Information: Suppl. Files 1, 2, 3 and 4. 1 INTRODUCTION Next-generation sequencing (NGS) platforms can sequence a particular transcriptome in a fast and cost-effective way (see for example Angeloni et al., 2011; Malik et al., 2011; Vera et al., 2008; Kristiansson et al., 2009). Assemblers for transcriptomics reads (e.g. Chevreux et al., 2004; Pavel A. Pevzner and Waterman, 2001; Li et al., 2009) will produce contigs (usually considered unigenes) that, in the ideal case, should be close to the number of real genes expressed under the studied conditions. When analysing de novo transcriptomes, it is quite difficult to distinguish between artefacts, ⇤to whom correspondence should be addressed true new transcripts, and the estimated number of identified transcripts. Annotations cannot be relied on for distinguishing between them since they are based on sequence similarity, and most of the problematic cases will not have any orthologue in databases. In fact, in transcriptomics projects, 40-60% of unigenes do not match any similar sequences in databases (Paterson et al., 2010; Fern´ andez-Pozo et al., 2011; Parchman et al., 2010). Most of the times, researchers hypothesise that they derive from lineageor species-specific genes, but unfortunately a lot of them might simply be the result of misassembly or assembly of inaccurate reads (Paterson et al., 2010). For that reason, efforts should be invested in distinguishing real new genes from misassembled sequences. In de novo transcriptome assemblies, the identification of reconstructed transcripts coding for a complete protein has not received enough attention, even though they are extremely useful for (i) inferring the protein coded by the transcript, (ii) determining the genomic structure of genes in non-model organisms, (iii) facilitating gene identification efforts, and (iv) catalysing experimental research (Team, 2002). Identification of full-length cDNAs is a milestone in non-model species (Ralph et al., 2008; 1KP Project http: //www.onekp.com). No bioinformatics tool is available for identification of full-length transcripts in NGS projects. Tools such as TargetIdentifier (Min et al., 2005) or Full-Lengther (Lara et al., 2007), did exist for full-length transcript prediction in classical EST (expressed sequence tag) experiments. These softwares were prepared for low-throughput, cDNA-cloning based ESTs. In spite of their widespread use (Wang et al., 2010; Koop et al., 2008; Kim et al., 2008; Semova et al., 2006), both algorithms are ineffectual for working with unigenes from NGS for two reasons: (1) there is no clone supporting any sequence (which is the main rationale of these algorithms), and (2) they are not prepared for high-throughput, and cannot cope with new sequencing technologies. Another aspect of de novo transcriptome assembly is the constant development and improvement of NGS assemblers. But, there is no easy way to determine which is the best possible assembly when working with a non-model species. It is proposed here that this task could be assisted by determining the amount of long reconstructed transcripts and the number of contigs coding for different, complete proteins. c Oxford University Press 2012. 1 Sample et al We present FULL-LENGTHERNEXT, a tool designed to meet the need for de novo transcriptome assembly of non-model organisms described above. FULL-LENGTHERNEXT is adapted to NGS technologies and can work in parallel and in a distributed way to optimise computing resources. It can classify unigenes depending on their coding content, extract the putative protein from unigene sequences, suggest which unannotated unigenes could be speciesspecific coding sequences and is able to help in deciding which de novo assembly is the best for a non-model organism. Data supporting these capabilities are given. 2 APPROACH The overall algorithm underlying FULL-LENGTHERNEXT is depicted in Fig. 1. It divides the input file in chunks of chunkSize that are distributed in as many workers (cores) as were selected (see Suppl. File 1 for details about parameter meaning). When a worker finishes the analysis of a chunk, it receives a new chunk, and so on until all chunks have been analysed (Fig. 1A). Within each worker, sequences are compared successively against modules shown in Fig. 1A until a significant homology is found. Finally data and annotations (including the assigned status) are saved on disk. The presence of frame-shifts in unigenes is revealed by the partition of a reference sequence in several hits. In order to obtain a contiguous protein sequence for the unigene, FULLLENGTHERNEXT selects the correct frames (BLASTx hits) and joins them separated by the necessary Xs (representing unknown amino acids) to fill the gap between consecutive hits or to cover the overlapping residues in consecutive hits (‘Fix Frame Errors’ in Fig. 1B). Unigenes showing similarity in both strands with the same reference do not allow for a ORF reconstruction. They are tagged as misassembled (Fig. 1B) and are consequently discarded for further processing. Start and stop codons for the reconstructed ORF are then located in the unigene sequence based on the reference sequence. Depending on the cases found, FULL-LENGTHERNEXT will assign one of the following statuses (see “Suppl. File 2” for more detailed description): N-terminal for unigenes lacking the 3’-end coding region, C-terminal for unigenes lacking the 5’-end coding region, Complete for coding for a complete protein, Internal for unigenes lacking the 5’-end and 3’-end coding region; and Misassembled for unigenes with a clear missassembly. Unigenes classified as Complete are annotated with information from the same reference that enabled its analysis. Sequences that have passed through the three databases without retrieving any similarity could be putative species-specific genes or misassemblies (Paterson et al., 2010). TestCode algorithm (Fickett, 1982) implemented as a Ruby class was used for obtaining the probability that a sequence is coding. Since TestCode only renders reliable results for sequences longer than 200 bp, it was used for the analysis of unigenes whose predicted open reading frame was longer than 200 bp after the analysis of both strands (Fig. 1C). Unigenes can then be given the status of Coding if the TestCode score >0.74 (see “Suppl. File 2” for more detailed description). Unigenes with scores lower than 0.74 are candidates for non-coding sequences or misassemblies (including assemblies of artefactual reads). They, together with unigenes containing ORFs shorter than 200 bp, are then compared with the ncRNA database using BLASTn with an E-value of 103, to discern true non-coding RNAs from other putatively useless sequences. This analysis can provide the status of Putative ncRNA. Finally, unigenes without any similarity in any database and with a TestCode <0.74 receive the status Unknown. 3 METHODS 3.1 Architecture FULL-LENGTHERNEXT has been developed in Ruby on OSX and SLES Linux. It can be used as a command line, as a Web tool (http://www. scbi.uma.es/full_lengther2), or as a REST-based Web service. It can be installed as a Ruby gem using the command gem install full lengther next, or can also be downloaded from http://www. scbi.uma.es/downloads. Installation as a gem is preferred since other gem dependencies (such as SCBI MapReduce, scbi blast and scbi fasta) are automatically installed. Execution requires previous installation of Ruby 1.9.2 or higher, and BLAST+ (Camacho et al., 2009). FULLLENGTHERNEXT has been designed to use multiple parallel or distributed cores by using the Ruby gem SCBI MapReduce (Guerrero-Fernandez et al, submitted), based on one manager and a variable number of independent workers. The number of workers can be optimally made equal to the number of available cores in each machine or network. FULL-LENGTHERNEXT has been developed and tested using a cluster of 10 blades with 8 Xeon cores (x86 architecture) at 3.00 GHz and 16 GB of RAM per blade, and using PBS as queuing system. Details about launching FULL-LENGTHERNEXT and the resulting output is described in “Suppl. File 1”. 3.2 Reference databases FULL-LENGTHERNEXT provides one status to any sequence based on comparisons using curated protein databases and a heuristic algorithm. Protein databases contain only complete proteins: incomplete proteins were discarded on the basis of the feature field (incomplete proteins contain NON TER or NON CONS tags), description field (incomplete proteins contain Flags: Fragment), and the protein sequence (complete proteins must start by M). The first protein database is optional and must be constructed by the user. Its rationale is to gather the sequences closest to the studied organism in order to obtain the highest number of unigenes with an orthologue in this small protein database and reduce the time spent on data analysis. The FULLLENGTHERNEXT gem provides the script make user db.rb to build a user database; it only need specification of the database division (fungi, plants, invertebrates, etc.) and a specific taxonomic group. For example, the conifer-specific user database used below was constructed using the command: make user db.rb plants Coniferales. The following protein database consists of complete proteins derived from the UniProtKB/Swiss-Prot (Apweiler et al., 2004) database, including the different isoforms, from the same division in DBgroup mentioned above. The last complete protein database is derived from the non-redundant, automatically annotated UniProtKB/TrEMBL from the same division as above. Being the biggest database, it will be queried only by the small dataset of proteins that have not found a reasonably homologous hit in any previous database. Both customised databases can be automatically generated by means of the script download fln dbs.rb available with the FULL-LENGTHERNEXT gem. Sequences without similarity in any database and predicted as non-coding by FULL-LENGTHERNEXT are compared with a database of non-coding RNAs, constructed from the junction of Rfam (Griffiths-Jones et al. (2005)) and NONCODE (Bu et al. (2012)) databases and filtered with CDHIT (Li and Godzik (2006)) to avoid redundancy. The curated ncRNA database is downloaded when the above mentioned script download fln dbs.rb is 2 Fine-tuning de novo transcriptomes One Sequence >200bp? Find & Order ORFs Sense StrandAntisense Strand Separate SEQs By Strand Test Code Algorithm Get The Longest ORF Get The Longest ORF Get The Longest ORF NO YES Save Status >200bp? NO YES Find & Order ORFs ncRNA finding Coding Putative Coding Unknown Each worker Input file User DB SwissProt DB Trembl DB Coding gene prediction Print Saved Data In Files worker Split into chunks & distribution to workers worker worker worker ncRNA DB BLASTX Sequence chunk Reverse complement if necessary Fix Frame Errors Next DB Match? NO YES Misassembly? NO Find Start YES Find End Determine Status Complete?NO Save error Last DB? NO Save Annotations YES YES Coding gene prediction Any match? NOYES ncRNA finding A B C Fig. 1. Flowcharts setting out the FULL-LENGTHERNEXT algorithm. A: general scheme of the algorithm. B: details of the process involving the five modules that are executed within each worker in order to obtain one status for each sequence. C: a detailed scheme of the strategy for finding putative species-specific coding genes and non-coding RNAs; this is applied only to sequences without any similarity in protein databases. executed, and can also be downloaded from http://www.scbi.uma. es/downloads/FLNDB/ncrna\_fln\_100.fasta.zip. 3.3 Sequence datasets for validation Long-unigene dataset (LUD): this consists of a set of 10 000 sequences of 1202 bp in length (range: 6016-95 bp) on average from the putative transcriptome of Pinus pinaster (Fern´ andez-Pozo et al., 2011). Short-unigene dataset (SUD): this consists of a set of 10 000 sequences of 338 bp in length on average (range: 2145-40 bp) from the putative transcriptome of Pinus pinaster (Fern´ andez-Pozo et al., 2011). Complete gene dataset from GenBank (CGD-GB): this are the first 1000 non-redundant sequences from 5000 GenBank Homo sapiens sequences containing the description complete cds filtered with CD-HIT (Li and Godzik, 2006) with a similarity score of 80%. They were used as true positives (TP) for Complete status. The accession number of all sequences in CGD-GB can be found in the “Suppl. File 2”. Complete gene dataset from the mammalian gene collection (CGD-MGC): this are the first 1000 non-redundant sequences from 5000 Homo sapiens full-length sequences ftp1.nci.nih.gov/fasta/hs_mgc_mrna. fasta.gz filtered with CD-HIT as above. They were also used as another set of TP for Complete status. The accession number of all sequences in CGD-MGC can be found in the “Suppl. File 2”. Incomplete gene dataset (CG-Trimmed): all sequences in CGD-MGC dataset were modified by trimming 100 bp from each one, 50 bp from the start codon and 50 bp from the end codon. Resulting sequences are expected to be devoid of their start and stop codons and should therefore be considered an internal region of the original coding genes. This dataset served as true negatives (TN) for any status of incomplete sequence. Complete gene dataset simulating an assembly of paired-end reads (CG paired): assembly of paired-end reads providing transcript scaffolds containing internal blocks of Ns. To simulate such a result, 981 complete sequences qualified as sure Complete from CGD-MGC (Table 1) were disrupted at a random place in the second third of the sequence, with a random number of Ns (5, 10, 15, 20, 25, 30, 35, 40, 45 or 50). Ns substitute individual nucleotides so that the sequence length is preserved and a true gap is introduced. They were used as TP for Complete when sequences contain blocks of Ns. Complete gene datasets with indels: the 981 sequences qualified as sure Complete from CGD-MGC (Table 1) were manipulated to create 16 new datasets of 981 sequences each by introducing artificial insertions, deletions, and insertions and deletions (see details in Suppl. File 3). They were used as TP for Complete status for tolerance to indels. Hence, Non-coding sequence dataset (NCSD): it consists of 1000 intergenic sequences from Arabidopsis downloaded from an intergenic dataset from TAIR site. NCSD was used as TN for Coding status. 3 Sample et al Log CPU Number 0 5 10 15 20 Time (Hours) SUD SUD+UserDB LUD LUD+UserDB 8 CPUs 4 CPUs 2 CPUs Fig. 2. Efficiency of time usage by FULL-LENGTHERNEXT illustrated as the time (in hours) spent using 2, 4 and 8 CPUs for the analysis of SUD and LUD datasets, depending on the inclusion of the optative user database. Pyrosequencing reads from a transcriptome (NGS-T): this consists of 746 105 useful reads from a pyrosequencing reaction of the Pinus pinaster transcriptome described in (Fern´ andez-Pozo et al., 2011). 4 RESULTS AND DISCUSSION 4.1 Time efficiency Computation time spent by FULL-LENGTHERNEXT algorithm was assessed with the LUD and SUD datasets of real-world sequences with 2, 4 and 8 CPUs and with and without the conifer-specific user database, UniProt databases being in any case from the plant division. It can be shown (Fig. 2) that the shorter the sequences, the faster the analysis, and that the use of the user database decreases significantly the time spent on analysing datasets of large sequences. The inclusion of a species-specific database allowed 5496 LUD sequences (55%) and 2676 SUD sequences (27%), which could account for the speed-up observed with long sequences. The use of a user database is also providing an additional advantage: it offers more accurate results due to a comparison with closer orthologues (results not shown). In our laboratory, we have built transcriptomes for different non-model species, and the analyses took from 14 h 32 min using 8 cores for 89 544 unigenes (average length: 622 ±563 nt) to 8 h 25 min using 5 cores for 133 175 unigenes (average length: 256 ±149 nt), taking into account that execution time will depend on the lengths of unigenes, the amount of similarity matches in databases, and the load of the queuing system. In conclusion, in only a few hours FULL-LENGTHERNEXT analysis can provide a general picture of a transcriptome assembly before the heavy investment of time that is needed for a deep annotation of a transcriptome with a dedicated software. 4.2 Reliability of Complete status When GenBank sequences (CGD-GB) were analysed with FULLLENGTHERNEXT, 27 of the 1000 sequences were not classified as Complete (Table 1). A detailed inspection of these sequences (see Suppl. File 4 for a specific description) revealed that 15 were correctly qualified by FULL-LENGTHERNEXT since the testing sequences were devoid of the true start or stop codon; 6 cases were mispredicted because long isoforms usually give better scores than short isoforms; and the 2 Coding statuses were correct, even though no reference sequence was found in the databases. Depending on the criteria, a minimum of 4 and a maximum of 12 could be considered FULL-LENGTHERNEXT failures. However, the accumulation of sequence errors in CGD-GB prompted for another full-length reference dataset source with more accurate metadata. This alternative reference was found in The Mammalian Gene Collection (CGD-MGC), where only 9 of 1000 sequences were apparently misclassified, all of them as C-terminal (Table 1). Again, 6 were originated by the preference of long isoforms over short isoforms, 2 were missing the 5’-end of the gene, and the last was a chimeric sequence (see Suppl. File 4 for s detailed description). Since mispredictions could range from 0 to 6 and there are less misannotations, the CGD-MGC dataset seems to be more appropriate for FULL-LENGTHERNEXT tests. This is also supported by the fact that, during the analysis, the warning messages (1) ”frame errors” appeared 13 times in CGD-MGC and 48 times in CGD-GB; (2) ”unexpected stop codons” appeared 3 times in CGDMGC and 18 times in CGD-GB; and (3) ”no M at the beginning” appeared 7 times in CGD-MGC and 18 times in CGD-GB. Incomplete sequence status was tested using a collection of true incomplete proteins (CG Trimmed dataset). The analysis revealed that 983 were distributed amongst different incomplete statuses and 7 were Unknown (Table 1). Four cases were misclassified as sure Complete because they were similar to unknown TrEMBL proteins of small size (116 amino acids in length on average). Six cases were classified as Putative complete because the 5’ and 3’ deletions removed only a few amino acids from the predicted protein; in fact, FULL-LENGTHERNEXT gives warning messages indicating that these testing sequences are devoid of start and stop codons but could code for a protein of a size close to the reference. Overall reliability of FULL-LENGTHERNEXT was then assessed using the results for the CGD-MGC and CG Trimmed datasets (Table 1). Their corresponding sensitivity, specificity, precision and accuracy (Table 2, first row) were over 0.99, close to the ideal value of 1.0, suggesting that the heuristic algorithm of FULL-LENGTHERNEXT produces reliable positive and negative results. The low relative error quotient (REQ = 0.00807; Martin et al., 2004) supports that FULL-LENGTHERNEXT is a useful tool for the automatic and fast identification of sequences coding for complete and incomplete proteins. Further testing of the capabilities of FULL-LENGTHERNEXT was performed with new datasets obtained from the 981 sequences in CGD-MGC that received the sure Complete status, excluding the 9 entries classified as C-terminal and the 10 entries classified as Putative complete. 4.3 Robustness against sequence mistakes The main source of mistakes in datasets were sequence indels, which provoke frame-shifts giving rise to false start and/or stop codons. Since FULL-LENGTHERNEXT should be able to cope with these errors, the CG Ins, CG Dels and CG Indels datasets (see Suppl. File 3) have been created for testing this possibility. Most sequences in these datasets received a Complete status (Table 1), with sensitivity, specificity, precision and accuracy values over 4 Fine-tuning de novo transcriptomes Table 1. Status prediction for all sequence datasets used for FULL-LENGTHERNEXT validation. Testing dataset Complete N-Terminal C-Terminal Internal Coding ncRNA Unknown SPSP SP SP Complete transcripts CGD-GB 940 33 3120 0 0 20 0 1 CGD-MGC 981 10 00 90 000 0 0 Incomplete transcripts CG Trimmed 465 20 2 23 933 00 0 7 Full-length artificial paired-ends sequences CG paired 968 11 00 10 100 0 0 Complete transcripts with insertions CG Ins1 966 12 00 30 000 0 0 CG Ins2 961 16 00 40 000 0 0 CG Ins3 954 19 40 40 000 0 0 CG Ins4 969 8 00 40 000 0 0 CG Ins5 966 10 00 50 000 0 0 CG Ins6 920 53 20 42 000 0 0 Complete transcripts with deletions CG Dels1 957 17 10 51 000 0 0 CG Dels2 951 25 10 40 000 0 0 CG Dels3 956 20 00 50 000 0 0 CG Dels4 949 23 40 50 000 0 0 CG Dels5 963 14 10 30 000 0 0 CG Dels6 959 17 10 40 000 0 0 Complete transcripts with indels CG Indels1 951 23 20 32 000 0 0 CG Indels2 968 12 00 10 000 0 0 CG Indels3 957 20 00 40 000 0 0 CG Indels4 930 46 10 40 000 0 0 Coding status validation1 CGD-MGC ---- -- - 725 211 29 35 NCSD ---- -- - 10 58 4 928 S: sure; P: putative 1CGD-MGC (coding sequences) and NCSD (non-coding sequences) datasets were used for testing the reliability of the Coding status based on TestCode score. FULL-LENGTHERNEXT was slightly modified to ignore the comparison with any protein database to simulate the situation of none of the testing sequences matching any hit in databases, meaning that only the C part of the algorithm depicted in Fig. 1 was used. 0.99 (Table 2), and a very low REQ (<0.0039). Values were statistically equal (P>0.05) with respect to the original sequences, demonstrating that up to 3 indels events in a sequence did not significantly alter FULL-LENGTHERNEXT predictions. Values were also statistically equal (P>0.05) when all modifications were computed as a single case (Table 2, CG all indels row). Another source of confusion could be the presence of blocks of Ns within sequences. This usually occurs when the transcriptome has been constructed using paired-end reads. The capability of FULLLENGTHERNEXT to deal with this type of sequence was tested with the artificial dataset CG paired. Again, most sequences in this dataset were assigned a Complete status, with sensitivity, specificity, precision and accuracy values over 0.989 (Table 2), and a very low REQ (0.00613), also providing statistically significant equality (P>0.05) to the case of sequences without those N blocks (Table 2). The success of FULL-LENGTHERNEXT predictions regarding indels and N blocks indicates that it is is useful for working with real results from tentative unigenes obtained by single-reads and 5 Sample et al Table 2. Summary of positives and negatives observed in FULL-LENGTHERNEXT tests and its use in algorithm reliability validation Testing dataset TP TN FP FN Sensitivity Specificity Precision Accuracy REQ CGD-MGC 991 2 - 6 0.99398 0.99002 0.99001 0.99199 0.00807 CG Trimmed1990 10 - ----- CG ins 5854 - - 32 0.99456 0.99002 0.99829 0.99390 0.00359 CG dels 5851 - - 35 0.99405 0.99002 0.99829 0.99347 0.00385 CG indels 3907 - - 17 0.99567 0.99002 0.99745 0.99452 0.00346 CG all indels215612 - - 84 0.99465 0.99002 0.99936 0.99437 0.00301 CG paired 979 - - 2 0.99796 0.99002 0.98989 0.99394 0.00613 coding 936 - - 64 0.93600 0.92800 0.92857 0.93200 0.07265 non-coding - 928 72 - ----- TP: True Positives; TN: True Negatives; FP: False Positives; FN: False negatives. Sensitivity: the proportion of actual positives, which are correctly identified, as TP/(TP+FN). Specificity: the proportion of negatives, which are correctly identified, as TN/(TN+FP). Precision: the degree to which repeated measurements under unchanged conditions show the same results, that is, repeatability or reproducibility of the measurement, as TP/(TP+FP). Accuracy: proximity of results to the true value, as (TP+TN)/(FP+TP+FN+TN). Relative error quotient (REQ): a value for assessing prediction methods as (FN*W+FP)/TP(1+W), where low REQ represents a low proportion of errors and higher REQ indicates a higher proportion of errors (Martin et al., 2004). 1CG Trimmed dataset was used as source of TN and FP for the CGD-MGC dataset as well as for the derived datasets containing insertions and/or deletions, which is why specificity is the same for the first six values. 2CG all indels values are the sum of data of all insertion, deletion and indel datasets. paired-reads assembly, while coping with the presence of sequence mistakes. 4.4 Reliability of Coding status It is in the aim of FULL-LENGTHERNEXT to determine which could be species-specific genes amongst the sequences without similarity in databases. This goal was achieved by the piece of algorithm represented in Fig. 1C. The CGD-MGC dataset was used as source of TP and the NCSD as source of TN. The results of this test can be found in Table 1 (rows below “Coding status validation”). Sensitivity, specificity, precision and accuracy values were over 0.928 (Table 2), values that were statistically lower (P<0.05) than for the Complete status, but still close to the theoretical value of 1.0. As REQ is still low (0.07265; values of REQ <0.4are considered good predictors [Martin et al., 2004]), it can be suggested that Coding status is a successful indication that the transcript is likely to be species-specific and not a result of sequencing or assembling artefacts, laboratory manipulation, sequencing errors, inaccurate pre-processing (Falgueras et al., 2010) or misassembling. The Coding status does not refer to the completeness of the sequence, since this deep analysis lacked reliability (results not shown). Specificity and accuracy values (Table 2) suggest that most unigenes that do not receive the Coding status could be discarded for the analysis. However, the TestCode score is based on protein coding genes, although it is now known that there is an increasing number of biologically relevant long and short non-coding RNAs (ncRNAs) in deeply sequenced transcriptomes (Huang et al., 2011). In order to recover likely ncRNAs, unigenes without a status after TestCode analysis were compared with a non-redundant ncRNA database (see Reference databases for details). Unigenes without a ncRNA homologue will receive a final Unknown status. The reliability of Coding and Unknown statuses have been further studied with unigenes in the “CAP3” column in Table 3, corresponding to a FULL-LENGTHERNEXT analysis of realworld data in NGS-T (see below). Interestingly, the mean length of unigenes having an Unknown status is 315 nt (mode = 31 nt), while the mean length of unigenes with any other status is 756 nt (mode = 341 nt). This is evidence that unigenes with Unknown status might correspond to short sequences resulting from artefacts or misassemblies, or sequences so short that BLAST was not able to detect any significant similarity. Additionally, unigenes with the status of Coding (1495 + 1195) and Unknown (9619) were analysed using AutoFact (Koski et al., 2005) with the same stringency than FULL-LENGTHERNEXT (E-value <1025). Although AutoFact is not maintained and is not designed for highthroughput analysis, it was useful for this purpose since it uses up to 10 databases, including EST, KEGG, Pfam, UniRef90 or NR, for annotation. It must be mentioned that these databases were not filtered for complete proteins, so they may contain sequences that are absent from FULL-LENGTHERNEXT reference databases. A total of 63.9% and 61.5% of unigenes with the sure Coding or Putative coding status, respectively, were annotated with AutoFact. Annotations were mainly derived from EST databases [96.5% (1119 unigenes) in sure Coding and 99.4% (911 unigenes) in Putative coding], suggesting that these unigenes could derive from species-specific genes. As for unigenes with the Unknown status, 43.1% (4146 unigenes) were annotated by AutoFact, all of them but 1 from EST databases, suggesting that the remaining 5473 unigenes (56.9%) are strong candidates for being artefacts or misassemblies. Due to their short length, they can be discarded for further annotation in order to accelerate the process without a significant loss of information. However, users are invited to perform a relaxed annotation of Unknown unigenes using dedicated annotation tools in order to recover the unigenes that might contain useful information. 6 Fine-tuning de novo transcriptomes Table 3. Comparison of NGS-T sequences from a pine transcriptome assembled with Euler-SR, MIRA and CAP3 assemblers Euler-SR MIRA3 CAP31 #seqs % #seqs % #seqs % Unigenes 24147 100% 50813 100% 41246 100% Unigenes >500 bp 11729 48.57% 18453 36.32% 20115 48.77% Longest unigene (bp) 4483 7450 7450 With ortologue216779 69.49% 35334 69.54% 28314 68.65% Different ortologue IDs 9997 59.58% 14334 40.57% 13452 47.51% Complete transcripts 2766 16.49% 6328 17.91% 6164 21.77% Different complete transcripts 2606 15.53% 4055 11.48% 4640 16.39% Misassembled 11 0.07% 139 0.39% 30 0.11% Without ortologue27362 30.49% 15479 30.46% 12932 31.35% Coding 785 10.66% 2191 14.15% 1815 14.03% Putative coding 492 6.68% 1772 11.45% 1490 11.52% Putative ncRNA 6 0.08% 5 0.03% 8 0.06% Unknown 6085 82.65% 11511 74.37% 9619 74.38% Mapped reads3443 976 59.50% 637321 85.41% 651 643 87.33% 1Due to its overlap-layout-consensus design for Sanger sequencing, CAP3 cannot be used with the huge amount of reads provided by any NGS method. It has been therefore used for reconciliation of the assemblies obtained from Euler-SR and MIRA3. 2Percents for subclassifications of this category were calculated using this line as 100% reference. 3Mapping was performed with Bowtie 2.0 using the default parameters (Langmead et al., 2009) and the useful reads as input. 4.5 Identification of the best de novo assembly of a transcriptome Transcriptomics projects based on de novo assemblies usually lack external criteria for discerning which is the best possible assembly of sequencing reads. In this context, FULL-LENGTHERNEXT could be used as a first approach for identifying unigenes derived from true transcripts by means of the information contained in the summary stats.html file (see “Suppl. File 1”). As mentioned above, sequences qualified as Misassembled can be discarded. Unigenes with an orthologue are good candidates for well reconstructed transcripts. Depending on the assembler or the coverage, real transcripts can be reconstructed as non-contiguous contigs, or as several overlapping unigenes that correspond to the same transcript. It is deduced that the sum of unigenes with an orthologue is giving the incorrect idea that more transcripts have been reconstructed than really exists. Consequently, a better approach might be to pay attention to the number of unique, orthologous IDs, since this will give a better idea of how many different transcripts have been reconstructed. The number of unigenes receiving the Coding status could also be a good approach for determining the real number of species-specific genes that contain the studied transcriptome. We propose that the sum of the number of Coding unigenes plus the number of different orthologue IDs is the best estimate of the amount of transcriptome that an experiment has revealed, since the higher the sum, the better the assembly appears to be. Unfortunately, this simple sum does not take into account the contiguity of the assembly. We propose that the amount of different Complete transcripts reconstructed can give an indication of the contiguity of the assembly when no transcript reference is available. In other words, the greater the number of different Complete unigenes, the better the assembly is expected to be, the more reliable the assembly can be considered to be. The preceding rationale can be illustrated with the real-world data of NGS-T obtained from Fern´ andez-Pozo et al., 2011, where MIRA3 was used for the assembly. In Table 3, the ‘MIRA3’ column is an assembly equivalent to this one published (Fern´ andezPozo et al., 2011), where the final number of unigenes with orthologues (50 813) might suggest that the complete transcriptome has been covered. However, the number of different orthologue IDs (14 334) reflects that there could still be genes to be discovered. The number of Coding unigenes (2191 + 1772) does not account for the missing genes if one assumes that the pine genome, such as Arabidopsis, could have ⇠25 000 genes. Since the objective of MIRA3 is to avoid the co-assembly of different alleles or paralogues, the number of revealed transcripts could be overestimated in the ‘MIRA3’ column of Table 3 and the cited report (Fern´ andez-Pozo et al., 2011). Therefore, a different assembly was performed using the Eulerian assembler Euler-SR (Pavel A. Pevzner and Waterman (2001)), the algorithm of which is completely different from the overlap-layout-consensus developed for MIRA3. As expected, this assembly (‘Euler-SR’ column, Table 3) provides lower values than the MIRA3 column for all parameters, since Euler-SR is more sensitive to sequencing errors or allelic variations and stops transcript assembly when data do not clearly favour a sequence instance. However, percentages of unigenes > 500 bp (48.57%), of different orthologue IDs (59.58%) and different complete transcripts (15.53%) are higher, suggesting that this assembly reduces the multiplicity of unigenes for the same original transcript. Consequently, Euler-SR seems to provide the most reliable unigenes contained in NGS-T. The suggested interpretation is that Euler-SR reveals the minimum number of unigenes and MIRA3 reveals and overestimated number of unigenes, while the real number of revealed unigenes lies between both values. A reconciliation of both assemblies with CAP3 has been carried out, since CAP3 has been described as a highly reliable de novo 7 Sample et al BC002469, BC003607, BC000680, BC001040, BC002472, BC000709, BC007313, BC001140, BC009746, BC011596, BC006527, BC002497, BC001471, BC002492, BC000623, BC001586, BC001421, BC001727, BC007314, BC000683, BC000096, BC002950, BC000107, BC001863, BC001141, BC002560, BC007315, BC000024, BC004402, BC001446, BC000711, BC000126, BC001042, BC003547, BC008036, BC009752, BC001595, BC014390, BC019252, BC000051, BC001627, BC004119, BC003087, BC001658, BC007660, BC000694, BC001659, BC017365, BC003088, BC001613, BC001142, BC000025, BC001451, BC003548, BC002446, BC001866, BC001626, BC009763, BC006194, BC001398, BC017175, BC000635, BC001621, BC000566, BC001144, BC014636, BC003549, BC001145, BC000104, BC004349, BC002442, BC001044, BC003550, BC001109, BC001110, BC007672, BC001013, BC025414, BC017366, BC003551, BC004420, BC003552, BC002464, BC008767, BC001146, BC000567, BC001870, BC000101, BC000713, BC001689, BC012316, BC008740, BC002506, BC001147, BC008037, BC002444, BC001022, BC001871, BC001033, BC002448, BC000706, BC007656, BC000036, BC018823, BC017367, BC001387, BC000714, BC001690, BC000715, BC000601, BC001440, BC001734, BC003608, BC001148, BC001625, BC009175, BC022845, BC002520, BC008747, BC012530, BC008764, BC002609, BC002443, BC021965, BC019253, BC001873, BC001874, BC007317, BC007318, BC007319, BC001765, BC008039, BC001376, BC000716, BC000717, BC021192, BC003553, BC001046, BC003090, BC006493, BC007320, BC002447, BC001482, BC002453, BC004129, BC001355, BC003060, BC002955, BC003554, BC002487, BC002439, BC017369, BC001047, BC001766, BC014392, BC002586, BC007655, BC001444, BC002532, BC000091, BC000718, BC003609, BC008763, BC007321, BC001152, BC003091, BC001435, BC001767, BC000719, BC011599, BC017188, BC002488, BC004963, BC001660, BC000649, BC001441, BC006534, BC000632, BC001472, BC009806, BC003610, BC001484, BC001878, BC009808, BC001731, BC002956, BC000670, BC002489, BC001485, BC000062, BC001025, BC016155, BC003092, BC001423, BC003061, BC017178, BC004130, BC004132, BC001880, BC021967, BC014638, BC002578, BC008770, BC000102, BC000682, BC001438, BC001111, BC000590, BC001154, BC002959, BC005136, BC001881, BC001486, BC001661, BC010878, BC007323, BC001036, BC001768, BC006494, BC016758, BC001691, BC006535, BC002429, BC005137, BC009177, BC003613, BC000639, BC003555, BC004421, BC003556, BC021193, BC001155, BC000544, BC005115, BC000019, BC000720, BC001113, BC001156, BC001769, BC001386, BC002436, BC001466, BC001380, BC003093, BC001157, BC013348, BC001114, BC001454, BC001158, BC001883, BC001050, BC001884, BC001662, BC000235, BC001588, BC001771, BC000127, BC002960, BC008734, BC006495, BC007333, BC014301, BC002513, BC008906, BC000028, BC016027, BC000283, BC002479, BC000187, BC007340, BC005138, BC002461, BC001596, BC003062, BC001159, BC004136, BC003094, BC003615, BC001031, BC000069, BC000721, BC002573, BC001886, BC000722, BC002544, BC004137, BC002527, BC000674, BC000093, BC000629, BC000723, BC015557, BC003616, BC002515, BC000724, BC003557, BC008726, BC004352, BC014397, BC002549, BC001161, BC001887, BC001359, BC000010, BC013903, BC001634, BC005139, BC007348, BC001738, BC002567, BC001633, BC002485, BC000725, BC000097, BC020170, BC007349, BC002962, BC000627, BC005116, BC000652, BC021085, BC000070, BC002505, BC000038, BC001692, BC000581, BC008765, BC001163, BC001744, BC002503, BC001708, BC001392, BC000023, BC004138, BC001019, BC000603, BC001773, BC003096, BC017452, BC001052, BC001720, BC004892, BC008727, BC002585, BC003617, BC008929, BC000576, BC000661, BC023975, BC001164, BC017198, BC002577, BC005140, BC001888, BC000057, BC001664, BC023503, BC017199, BC014561, BC013968, BC001756, BC001422, BC002426, BC006496, BC000249, BC004964, BC000250, BC008750, BC001889, BC001165, BC000011, BC001115, BC002591, BC000114, BC000594, BC002965, BC000654, BC002539, BC000726, BC008777, BC004139, BC000660, BC003097, BC003098, BC000727, BC000578, BC001491, BC008730, BC017378, BC004140, BC017564, BC009235, BC005827, BC009503, BC000251, BC000631, BC002454, BC000117, BC001890, BC002579, BC004141, BC005141, BC001166, BC011601, BC021981, BC008930, BC017194, BC000054, BC003064, BC003558, BC000252, BC002599, BC002445, BC000185, BC001891, BC001492, BC025372, BC002555, BC000671, BC006458, BC006459, BC000644, BC003065, BC014563, BC007401, BC000728, BC000278, BC016759, BC004143, BC003100, BC000006, BC002530, BC000043, BC017565, BC007402, BC000260, BC002967, BC000729, BC001116, BC015961, BC006195, BC003559, BC000588, BC000213, BC001732, BC003560, BC001693, BC019256, BC001777, BC000039, BC001493, BC021190, BC001892, BC000574, BC001167, BC001712, BC019257, BC001665, BC000012, BC002968, BC017197, BC000657, BC001752, BC000253, BC000094, BC002970, BC001666, BC000550, BC016760, BC001778, BC008768, BC001055, BC001894, BC002432, BC002441, BC001465, BC003101, BC001169, BC003102, BC001056, BC000642, BC006460, BC000658, BC001416, BC000730, BC003103, BC001391, BC001779, BC011604, BC001808, BC000609, BC013382, BC003104, BC001709, BC001117, BC004965, BC017193, BC001599, BC007491, BC001494, BC017187, BC003619, BC009894, BC005143, BC002511, BC000029, BC000692, BC002468, BC005830, BC008745, BC000645, BC003561, BC001057, BC023976, BC006537, BC001600, BC001118, BC000650, BC000731, BC001639, BC000098, BC002525, BC003105, BC000732, BC002618, BC000617, BC001119, BC001364, BC003106, BC002604, BC001898, BC003563, BC024043, BC002558, BC002564, BC000044, BC009470, BC001496, BC001395, BC001172, BC003067, BC001430, BC000598, BC003107, BC003620, BC002592, BC002434, BC000055, BC000597, BC000548, BC000075, BC000543, BC009895, BC013383, BC013910, BC012798, BC008741, BC023504, BC006538, BC002973, BC001780, BC008775, BC005145, BC002541, BC001453, BC000595, BC000045, BC009896, BC007662, BC001058, BC001449, BC014431, BC000614, BC001641, BC001358, BC001602, BC003565, BC001383, BC006539, BC001173, BC000733, BC000053, BC003566, BC001120, BC000734, BC003108, BC014564, BC000266, BC000571, BC001497, BC002600, BC009898, BC007659, BC003567, BC007403, BC000060, BC001900, BC002517, BC012333, BC003568, BC002582, BC001450, BC003569, BC000183, BC004893, BC006462, BC001901, BC001782, BC006463, BC000046, BC001439, BC001642, BC000736, BC018648, BC003622, BC014565, BC000234, BC014566, BC000737, BC000738, BC001174, BC002538, BC001462, BC002430, BC018649, BC000739, BC001060, BC001123, BC001632, BC001620, BC003110, BC007404, BC001410, BC008061, BC007405, BC003623, BC036703, BC014433, BC006497, BC003068, BC005147, BC001414, BC005831, BC006196, BC002975, BC000254, BC002614, BC001394, BC002427, 14 Fine-tuning de novo transcriptomes BC002455, BC000092, BC001431, BC000033, BC002976, BC002594, BC002536, BC003624, BC002977, BC001432, BC000740, BC001500, BC008062, BC001903, BC002433, BC002476, BC002546, BC001176, BC002978, BC001721, BC001177, BC001125, BC000741, BC001630, BC001619, BC000233, BC001501, BC000665, BC001502, BC000690, BC003112, BC001178, BC001459, BC000086, BC008063, BC005832, BC000545, BC001670, BC000562, BC001503, BC015558, BC001605, BC001388, BC001746, BC001741, BC000232, BC000080, BC001904, BC012799, BC000182, BC001504, BC001372, BC001606, BC007407, BC006465, BC000095, BC017455, BC002584, BC001023, BC001754, BC008749, BC008064, BC002601, BC001785, BC001786, BC015969, BC003113, BC001179, BC003575, BC000013, BC002557, BC001061, BC003625, BC016320, BC002576, BC007493, BC001643, BC001505, BC011606, BC001180, BC001034, BC006498, BC006540, BC001360, BC002523, BC002979, BC000742, BC004966, BC004147, BC006541, BC019260, BC000211, BC017174, BC004967, BC000261, BC001506, BC008933, BC002466, BC000693, BC003576, BC006499, BC001644, BC008065, BC006197, BC021985, BC018652, BC006542, BC008934, BC001062, BC006501, BC023985, BC004151, BC000745, BC019236, BC004153, BC004154, BC006543, BC001063, BC006544, BC006504, BC008937, BC000747, BC006505, BC004155, BC001064, BC013923, BC000214, BC004156, BC018654, BC008938, BC007408 15 Sample et al SUPPLEMENTARY FILE 4: PATTERN OF INSERTION AND/OR DELETIONS TO OBTAIN THE COMPLETE DATASET WITH INDELS Manipulation of the 981 sequences qualified as sure Complete from CGD-MGC (Table 1) to create 16 new datasets of 981 sequences each by introducing artificial insertions and/or deletions as depicted in the following figure: Fig. 3. Pattern of artificially-introduced insertions and deletions in every sequence of CGD-MGC dataset, generating the 16 datasets of 981 sequences each indicated at the left. Note that in no case is the initial reading frame recovered after the indel combination. The obtained sequence datasets were the following: •CG Ins1 to CG Ins6: datasets including insertions of 1-2 nt in the first third, middle and/or last third of the sequence, shifting the reading frame one or more times. In some cases, there are up to three insertion events per sequence. •CG Dels1 to CG Dels6: datasets equivalent to the preceding ones, but with deletions of 1-2 nt instead of insertions. •CG Indels1 to CG Indels4: datasets including combinations of insertions and deletions of 1-2 nt in the first third and last third of the sequence for shifting twice the reading frame. 16 Fine-tuning de novo transcriptomes SUPPLEMENTARY FILE 5: DETAILED DESCRIPTION OF OF MISCLASSIFICATIONS OF FULL-LENGTHERNEXT In CGD-GB: A total of 27 sequences in CGD-GB were misclassified (Table 1), 20 of them being C-terminal,4N-Terminal, 2 were simply Coding and 1 was Unknown. Four (BC111790, BC103905, BC080189, BC015371) of the 20 C-terminal sequences were misclassified because the top-BLASTx subject (the subject with the highest score and E-value) corresponds to a confounding reference of the protein (Suppl. Fig. 4A) because it has an E-value and/or score slightly higher than the right reference sequence. For example, BC111790 gives an E-value = 0.0 for Q5SW24 (long reference) and Q5SW24-2 (short reference), but their scores differ (1068 for Q5SW24 and 1061 for Q5SW24-2). BC111790 should be analysed using the shorter reference Q5SW24-2 because identity spans the complete length of both sequences. However, BC111790 is also identical to the long reference Q5SW24, but only downstream of the M corresponding to its initial M, while the sequence upstream of this M has no significant similarity (only 7 similarities in 21 positions). The lack of similarity among those 21 amino acids while the rest are identical is a remarkable behaviour that can be explained only if the long reference Q5SW24 is considered to contain a 21 N-terminal amino acid-long extra sequence. The spurious similarity in the extra region accounts for the higher score in Q5SW24 with respect to Q5SW24-2 (Suppl. Fig. 4A.1). Moreover, this dual behaviour of similarity also prompts to select reference Q5SW24-2, and not the reference Q5SW24, as this is the right orthologue of BC111790. Similarly, 2 (BC104884, BC063037) of the misclassified 4 N-terminal sequences did have an equivalent explanation because the top-BLASTx subject was a confounding reference longer than the queries at the 3’ region (Suppl. Fig. 4A.2). The 6 sequences treated as cases of Suppl. Fig. 4A were mispredicted because the query sequence seems to be incomplete since it is shorter than a top-BLASTx confounding reference although, in fact, it contains the correct reference and should be classified as Complete. Other cases of apparent misprediction correspond to the other 2 sequences as N-terminal (Table 1), where the query sequence contained a frame-shift that misled to a premature stop codon (Suppl. Fig. 4B.1); therefore, the query sequence did not code for a complete protein since the C-terminal part is absent. As a result, the FULL-LENGTHERNEXT status was right while the query sequence metadata of complete cds were wrong. In 13 of the 20 cases where FULL-LENGTHERNEXT status was C-terminal (Table 1), the status was right because the query sequence does not contain the beginning of the gene: a frameshift mistakenly rises a false ATG as the start codon (Suppl. Fig. 4.B.2). Again, the FULL-LENGTHERNEXT was right while the sequence metadata were wrong. Finally, 2 of the 3 sequences that did not retrieve any reliable similarity in UniProtKB protein databases, were identified as coding sequences, suggesting that this part of the algorithm is reliable. For the remaining mispredicted sequences, it is difficult to assert whether the status or the metadata were right, since they contain too many sequencing errors and several subjects reveal different classes of status. 2FALSE STOP? FALSE ATG MEND STOP 5' 3' A M ATG Correct reference END STOP M 5' 3' FALSE STOP MEND 1ATG 5' 3' BRight prediction but wrong metadata M ATG STOP 5' 3' 1 2 END END N-terminal C-terminal N-terminal C-terminal Correct reference Confounding reference Wrong prediction by a confounding reference Confounding reference Fig. 4. Cases of complete sequences from CGD-GB classified as incomplete. The query (testing) sequence is in white and the reference sequence in database is in dark grey. A: cases where a true complete sequence was classified as (1) C-terminal because the query cannot contain an ATG located where required by the first M of the confounding reference, or (2) N-terminal because the query does not contain an stop codon at the position indicated by the confounding reference. B: cases where metadata of the query sequence indicated that it was complete but FULLLENGTHERNEXT predicts that it is an incomplete protein. In some cases (1), the query sequence finishes in a premature stop codon due to a frameshift. In other cases (2), the query sequence starts at an internal ATG, lacking the 5’-end of the gene; since it is also due to a frame-shift, it usually has an in-frame stop codon upstream of the false ATG. In summary, the 27 seemingly misclassified sequences from CGD-GB cannot be considered FULL-LENGTHERNEXT failures since 15 of them were misannotated and the prediction of FULLLENGTHERNEXT was the real status. Another 6 cases were analysed using a confounding reference and 2 out of the 3 sequences without UniProt orthologue were correctly classified as Coding. In CGD-MGC: In this dataset, only 9 sequences were apparently misclassified, all of them as C-terminal (Table 1). Two (BC008063 and BC023985) were truncated at the 5’ end, their metadata therefore being wrong and the FULL-LENGTHERNEXT status being right (as in Fig. 4B.2). Another one (BC001863) seems to be a chimeric sequence because it shows similarity with two different subjects, one of which is the sequence BC020999, which has the description WARNING: chimeric clone. The remaining 6 sequences receive a wrong status due to the use of a confounding reference (as in Fig. 4A). It is interesting to note that of the 10 sequences having a Putative complete status (Table 1), 7 (BC008039, BC003106, BC008367, BC017174, BC001223, BC002429, BC002461) did not start at a methionine, suggesting that the 5’-end of the sequence may not contain the true ATG. 17 146 CAP´ ITULO 11. AN´ ALISIS DE TRANSCRIPTOMAS CAP´ ITULO 12. ANOTACI´ ON DE SECUENCIAS GEN´ OMICAS 147 Cap´ıtulo 12 Prototipo para anotar secuencias gen´omicas incompletas Mi contribuci´on: Apoyo en la organizaci´on de la orientaci´on a objetos del c´odigo, paralelizaci´on del an´alisis y desarrollo de la interfaz web. 148 CAP´ ITULO 12. ANOTACI´ ON DE SECUENCIAS GEN´ OMICAS CAP´ ITULO 12. ANOTACI´ ON DE SECUENCIAS GEN´ OMICAS 149 T´ıtulo del art´ıculo: GENote v.: A Web Tool Prototype for Annotation of Unfinished Sequences in Non-model Eukaryotes Autores: No´e Fern´andez-Pozo, Dar´ıo Guerrero-Fern´andez, Roc´ıo Bautista, Josefa G´omezMaldonado, Concepci´on Avila, Francisco M. C´anovas, M. Gonzalo Claros Publicaci´on: Bioinformatics for Personalized Medicine - Lecture Notes in Computer Science. Springer Berlin Heidelberg Volumen:6620 A˜no: 2012 DOI:http://dx.doi.org/10.1007/978-3-642-28062-7_7 Abstract De novo identification of genes in newly-sequenced eukaryotic genomes is based on sensors, which are not available in non-model organisms. Many annotation tools have been developed and most of them require sequence training, computer skills and accessibility to sufficient computational power. The main need of non-model organisms is finding genes, transposable elements, repetitions, etc., in reliable assemblies. GENote v.is intended to cope with these aspects as a web tool for researchers without bioinformatics skills. It facilitates the annotation of new, unfinished sequences with descriptions, GO terms, EC numbers and KEEG pathways. It currently localises genes and transposons, which enable the sorting of contigs or sca↵olds from a BAC clone, and reveals some putative assembly inconsistencies. Results are provided in GFF3 format and in tab-delimited text readable in viewers; a summary of findings is provided also as a PNG file. 150 CAP´ ITULO 12. ANOTACI´ ON DE SECUENCIAS GEN´ OMICAS 151 Parte V DIFUSI ´ ON DE RESULTADOS MEDIANTE BASES DE DATOS Biology 2012, 1 458 83. Falgueras, J.; Lara, A.J.; Fernandez-Pozo, N.; Canton, F.R.; Perez-Trabado, G.; Claros, M.G. SeqTrim: A high-throughput pipeline for pre-processing any type of sequence read. BMC Bioinformatics 2010, 11, 38. 84. Guerrero-Fernaández, D.; Falgueras, J.; Claros, M.G. SCBI_MAPREDUCE: A task-farm, practical solution in Ruby for distribution of new and legacy bioinformatics software. IEEE Trans. Parallel. Distr. Syst. 2012, submitted for publication. 85. Paszkiewicz, K.; Studholme, D.J. De novo assembly of short sequence reads. Brief. Bioinform. 2010, 11, 457±472. 86. Nakamura, K.; Oshima, T.; Morimoto, T.; Ikeda, S.; Yoshikawa, H.; Shiwa, Y.; Ishikawa, S.; Linak, M.C.; Hirai, A.; Takahashi, H.; et al. Sequence-specific error profile of Illumina sequencers. Nucleic Acids Res. 2011, 39, e90. 87. Minoche, A.E.; Dohm, J.C.; Himmelbauer, H. Evaluation of genomic high-throughput sequencing data generated on Illumina HiSeq and genome analyzer systems. Genome Biol. 2011, 12, R112. 88. Hoffmann, S.; Otto, C.; Kurtz, S.; Sharma, C.M.; Khaitovich, P.; Vogel, J.; Stadler, P.F.; Hackermuller, J. Fast mapping of short sequences with mismatches, insertions and deletions using index structures. PLoS Comput. Biol. 2009, 5, e1000502. 89. Gilles, A.; Meglecz, E.; Pech, N.; Ferreira, S.; Malausa, T.; Martin, J.F. Accuracy and quality assessment of 454 GS-FLX Titanium pyrosequencing. BMC Genomics 2011, 12, 245. 90. Rasko, D.A.; Webster, D.R.; Sahl, J.W.; Bashir, A.; Boisen, N.; Scheutz, F.; Paxinos, E.E.; Sebra, R.; Chin, C.S.; Iliopoulos, D.; et al. Origins of the E. coli strain causing an outbreak of hemolytic-uremic syndrome in Germany. N. Engl. J. Med. 2011, 365, 709±717. 91. Balzer, S.; Malde, K.; Jonassen, I. Systematic exploration of error sources in pyrosequencing flowgram data. Bioinformatics 2011, 27, i304±309. 92. Miller, J.R.; Koren, S.; Sutton, G. Assembly algorithms for next-generation sequencing data. Genomics 2010, 95, 315±327. 93. Medvedev, P.; Pham, S.; Chaisson, M.; Tesler, G.; Pevzner, P. Paired de bruijn graphs: A novel approach for incorporating mate pair information into genome assemblers. J. Comput. Biol. 2011, 18, 1625±1634. 94. Compeau, P.E.; Pevzner, P.A.; Tesler, G. How to apply de Bruijn graphs to genome assembly. Nat. Biotechnol. 2011, 29, 987±991. 95. Earl, D.; Bradnam, K.; St. John, J.; Darling, A.; Lin, D.; Fass, J.; Yu, H.O.; Buffalo, V.; Zerbino, D.R.; Diekhans, M.; et al. Assemblathon 1: A competitive assessment of de novo short read assembly methods. Genome Res. 2011, 21, 2224±2241. 96. Huang, X.; Madan, A. CAP3: A DNA sequence assembly program. Genome Res. 1999, 9, 868±877. 97. Benzekri, H.; Bautista, R.; Guerrero-Fernández, D.; Claros, M.G. Departamento de Biología Molecular y Bioquímica, Facultad de Ciencias, Universidad de Málaga, 29071 Málaga, Spain, and Plataforma Andaluza de Bioinformática, Centro de Supercomputación y Bioinformática, Edificio de Bioinnovación, Universidad de Málaga, 29590 Málaga, Spain. Unpublished work, 2012. Biology 2012, 1 459 98. Lander, E.S.; Waterman, M.S. Genomic mapping by fingerprinting random clones: A mathematical analysis. Genomics 1988, 2, 231±239. 99. Aird, D.; Ross, M.G.; Chen, W.S.; Danielsson, M.; Fennell, T.; Russ, C.; Jaffe, D.B.; Nusbaum, C.; Gnirke, A. Analyzing and minimizing PCR amplification bias in Illumina sequencing libraries. Genome Biol. 2011, 12, R18. 100. Li, Z.; Chen, Y.; Mu, D.; Yuan, J.; Shi, Y.; Zhang, H.; Gan, J.; Li, N.; Hu, X.; Liu, B.; et al. Comparison of the two major classes of assembly algorithms: Overlap-layout-consensus and de Bruijn-graph. Brief. Funct. Genomics 2012, 11, 25±37. 101. FullLengtherNext. Available online: http://www.scbi.uma.es/fulllengthernext (accessed on 14 September 2012). 102. Loblolly Pine Genome Project. Available online: http://dendrome.ucdavis.edu/NealeLab/lpgp/ (accessed on 14 September 2012). 103. Díaz-Sala, C.; Cervera, M. Promoting a functional and comparative understanding of the conifer genome-implementing applied aspects for more productive and adapted forests (ProCoGen). BCM Proceedings 2011, 5, P158. 104. Kumar, S.; Blaxter, M.L. Comparing de novo assemblers for 454 transcriptome data. BMC Genomics 2010, 11, 571. 105. Sommer, D.D.; Delcher, A.L.; Salzberg, S.L.; Pop, M. Minimus: A fast, lightweight genome assembler. BMC Bioinformatics 2007, 8, 64. 106. Zheng, Y.; Zhao, L.; Gao, J.; Fei, Z. iAssembler: A package for de novo assembly of Roche-454/Sanger transcriptome sequences. BMC Bioinformatics 2011, 12, 453. 107. Iorizzo, M.; Senalik, D.A.; Grzebelus, D.; Bowman, M.; Cavagnaro, P.F.; Matvienko, M.; Ashrafi, H.; van Deynze, A.; Simon, P.W. De novo assembly and characterization of the carrot transcriptome reveals novel genes, new markers, and genetic diversity. BMC Genomics 2011, 12, 389. 108. Martin, J.; Bruno, V.M.; Fang, Z.; Meng, X.; Blow, M.; Zhang, T.; Sherlock, G.; Snyder, M.; Wang, Z. Rnnotator: An automated de novo transcriptome assembly pipeline from stranded RNA-Seq reads. BMC Genomics 2010, 11, 663. 109. Gnerre, S.; Maccallum, I.; Przybylski, D.; Ribeiro, F.J.; Burton, J.N.; Walker, B.J.; Sharpe, T.; Hall, G.; Shea, T.P.; Sykes, S.; et al. High-quality draft assemblies of mammalian genomes from massively parallel sequence data. Proc. Natl. Acad. Sci. USA 2011, 108, 1513±1518. 110. Simpson, J.T.; Wong, K.; Jackman, S.D.; Schein, J.E.; Jones, S.J.; Birol, I. ABySS: A parallel assembler for short read sequence data. Genome Res. 2009, 19, 1117±1123. © 2012 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/). Detecting and correcting mis-assembled reads in contigs Hicham Benzekri1, Darío Guerrero-Fernández1, Rocío Bautista1, and M. Gonzalo Claros1,2,* 1 Plataforma Andaluza de Bioinformática-SCBI, University of Málaga, C/ Severo Ochoa 34, 29590 Málaga, Spain {bhicham,dariogf,rociobm,claros}@uma.es http://www.scbi.uma.es 2 Molecular Biology and Biochemistry Department, University of Málaga, Campus de Teatinos s/n, 29071 Málaga, Spain http://www.bmbq.uma.es/fmp Abstract. De novo assemblies do not have the possibility of quality control with an external sequence. In fact, accuracy and reliability of these assemblies is highly affected by sequencing errors and mis-assemblies. Here, a frequencybased algorithm is developed in Ruby and intended to discern assembly errors from polymorphisms/read errors and then edit or remove the misassembled read(s) to provide more but highly reliable contigs. The software reads and writes the ACE assembly format. Transcriptome and genome assemblies were tested. Keywords: contigs, OLC, de novo, assembly. 1 Introduction Sequence assembly errors exist in any de novo assembly. Identification of misassemblies is a difficult issue due to the high amount of data and its error-prone quality because of biochemical and mechanical complications in sequencers. This usually requires additional efforts for manual validation of the most accurate reconstruction of the analyzed genome or transcriptome [1]. Too often, assembly quality is judged only by contig size or N50, with larger contigs being preferred [2], even though large contigs can be chimeric as a result of mis-assembling. A widely-used contig-testing tool is Hawkeye [3]. It can be used with assemblies of all sizes to facilitate the visual inspection of large-scale assembly data while minimizing the time needed to detect mis-assemblies and make accurate judgments for assembly quality. In fact, it guides users to the most likely areas of mis-assembly, allowing its manual edition and correction. In contrast to other contig editors such as GAP5 [4], Hawkeye combines computational predictors with interactive visualizations to decrease verification costs. However, visual inspection and manual edition are cumbersome tasks, particularly for assemblies from next-generation sequencing (NGS) data. This is the reason why amosvalidate [2], an automated * Corresponding author IWBBIO 2013. Proceedings Granada, 18-20 March, 2013 345 validation pipeline for contigs based on several independent criteria, was developed. But this tool only tagged regions that appear mis-assembled, and the correction requires visualization again and manual edition with Hawkeye. The aim of our work is to develop a fully automated algorithm called CoMiner with the aim of editing and correcting prominent mismatches in contigs, reducing the manual intervention dedicated to increase the quality of assemblies, based only on the contig assembly per se. 2 Implementation CoMiner was programmed in Ruby and tested in a dual core iMac at 3.06 GHz with 4 GB of RAM. Contig data can be read and written in ACE format [5], which is generated by various assembly programs, such as Phrap, CAP3, GAP4-5, Newbler, Arachne, Minimus and TIGR Assembler, all of them based on overlay-layoutconsensus algorithms (CoMiner is not ready for analyzing De Bruijn contigs, but could be adapted for mapping alignments in a near future). Fig. 1. Flowgram of CoMiner algorithm. A: Discovery of high entropy regions and identification of conflictive reads. B: Conflictive-read edition/removal within a contig depending on mismatch distribution. C: Once all conflictive reads in a contig have been edited, the contig coverage is verified, split in two or more contigs if necessary, and then saved into a new ACE file. Since the final goal of CoMiner is to increase the assembly accuracy without human intervention, the algorithm (Fig. 1) can be divided in the following main steps: (i) discovery of high-entropy regions (HERs); (ii) identification of conflictive read(s); (iii) read edition (trimming or removal); (iv) contig verification and saving in a new ACE file. (i) HER discovery: The aim of this step is to retrieve assembly fragments where reads do not align perfectly. We have elected the entropy of consensus nucleotide at position i [–H(i)] as a measure of the alignment goodness as described in [6]. Therefore, the frequency of each of the four nucleotides at each position of the Identification of conflictive/ candidate reads Mismatches at sequence ends Read trimming Read removal Mismatches spread over the read Regions with no coverage Regions with low coverage Write contigs & subcontigs in new ACE file Split contig into two or more subcontigs High entropy regions Frequency to entropy ACE file yes yes yes yes no no no no Reading of each contig Frequency table Read edition Read edition Next read Edited contig verification A B C IWBBIO 2013. Proceedings Granada, 18-20 March, 2013 346 consensus sequence is assessed and then used to calculate entropy at each consensus position. SNPs and point sequence are considered equivalent events in this rationale, and do not significantly affect assemblies unless they are closely located. This is the reason why entropy data were sieved by a fast Fourier transform as described [6] for converting contiguous sharp peaks into high entropy regions. Sequence ranges whose Fourier-transformed entropy is over a cutoff value that corresponds to the median entropy of the contig will be considered HERs and will focus subsequent analyses. (ii) Identification of conflictive-read(s): Several calculations are performed to determine whether a HER was caused by one or more mis-assembled reads or whether it was caused by mismatches scattered over all involved reads. These possibilities are discerned calculating the mismatch frequency of one read r [Fread(r)] as follows: for a contig containing m number of reads for which a k number of HERs have been defined, n(j) being the length of one HER, Fread(r) is calculated dividing the total number of mismatches of the read r within all HERs against the consensus by the total number of nucleotides involved in all HERs: The total mismatch frequency of the contig (Fcontig) is calculated dividing the total number of mismatches of every contig read within all HERs by the total length of all HER regions in each read, as follows: Both Fread(r) and Fcontig will define a robust Fcutoff value as: F read (r)=err(i,j) i=1 n ∑ j=1 k ∑ n(j) j=1 k ∑ F contig =err(i,j,r) i=1 n ∑ j=1 k ∑ r=1 m ∑ n(j,r) j=1 k ∑ r=1 m ∑ Fcutoff =Fcontig +K×Fread (r)−Fcontig r=1 m ∑ m $ % &' ( ) Fig. 2. Different instances of mismatch distribution in a conflictive read. A: All mismatches are located at one end; therefore, nucleotides within the distance d were trimmed from the read, provided that d < 40% of the read length, and the remaining read is longer than 40 nt. B: There are mismatches at both ends; again, nucleotides within distances d1 and d2 are trimmed provided that d1 + d2 < 40% of the read length, and the remaining read is longer than 40 nt. C: Mismatches are spread over the whole read; when the number of mismatches is over 2% of the read length, the complete read is removed. IWBBIO 2013. Proceedings Granada, 18-20 March, 2013 347 where K = 1.4826 to consider outliers only those values beyond the third quartile. Therefore, reads with Fread(r) > Fcutoff are considered conflictive and candidate for edition. (iii) Edition of candidate-read(s): The distribution of mismatches in candidate reads is analyzed as detailed in Fig. 2, driving to the trimming or removal of conflictive read(s) depending on mismatch distribution. (iv) Contig verification and saving: Read edition may modify the contig coverage and some region(s) may be now devoid of any read. The algorithm looks for this type of situations and splits the contig in two new, independent subcontigs, each one with a new, independent consensus sequence. There is a special case where a contig can contain an overlapping region of two reads flanked by a coverage of only one read (Fig. 3). This contig will be split into two independent subcontigs when the overlapping fragment is below 40 nt or the identity is below 90%. Unedited contigs, edited contigs and new subcontigs are then written into a new ACE file. Table 1: Results of two assemblies before (–) and after (+) CoMiner treatment Transcriptome Genome CoMiner – + – + Contig # 76 824 76 894 35 777 35 823 Mean contig size (nt) 429 428 471 470 N50 (nt) 484 482 506 505 N90 (nt) 251 252 307 307 Edited contigs 21 002 1465 Split contigs 146 74 Mapped contigs 35 412 35 458 Mapped nt 3 362 691 3 362 923 3 Results and Discussion CoMiner performance was tested for transcriptome and genome assemblies (Table 1). A total of 1 110 923 454/FLX reads from Solea senegalensis transcriptome were trimmed using SeqTrimNext (http://www.scbi.uma.es/seqtrimnext) [7] and then Fig. 3. Example of contig after read edition, where only two reads slightly connect two putative subcontigs. CoMiner will divide it in two new subcontigs. IWBBIO 2013. Proceedings Granada, 18-20 March, 2013 348 assembled using MIRA3 (http://www.scbi.uma.es/mira) with the standard parameters for 454/FLX data. In the resulting assembly (Table 1, Transcriptome columns), CoMiner detected HERs in 33 912 contigs (44.1%), but only edited 21 002 (27.3%), corresponding to contigs where at least one HER was caused by mismatches concentrated in at least one read. An example of this is shown in Fig. 4A, in which the algorithm identified a read containing several mismatches that were responsible for the wide, central HER. The right (3’) end of this read contained all mismatches and was consequently trimmed. A new entropy analysis after CoMiner automatic edition showed that the HER had disappeared (Fig. 4B), suggesting that the edited contig was more reliable than the initial one. Integrity of most contigs was unaffected by CoMiner edition, but 146 (0.7%) were split into two subcontigs and other 2 contigs were split into 3 different subcontigs each. It should be noted that when a subcontig consisted of only one read, it is not considered a contig and the read is removed from the final contig count. An example of contig splitting is shown in Fig. 5, where a 1570 bp contig was divided into two smaller subcontigs. Another example of this situation can be the chimeric contig group1_solea_c8241 of 1500 nt, since it was divided into a 5’ subcontig of 887 nt Fig. 4. Example of a 571 bp contig with several HERs before (A) and after (B) self-edition using CoMiner. After sieved entropy analysis, HERs spanning only one nucleotide are considered point errors or SNPs (in blue), while true HERs (in red) span two or more nucleotides. The wider HER (nt 273-296) is marked by a double arrow and magnified below, showing that all mismatches come from one single read. After edition (B), this wide HER disappeared. The other three smaller HERs were still present, suggesting that their mismatches were not concentrated in a single read and, therefore, will not be edited. IWBBIO 2013. Proceedings Granada, 18-20 March, 2013 349 with similarity to an ORM1 like protein (B0V340; E = 10–107) of 153 amino acids, and a 3’ subcontig of 601 nt without similarity in databases. As a result, automatic edition of 27.3% of contigs did not significantly increase the number of contigs and did not significantly change the general parameters of the assembly (Table 1, Transcriptome columns), while contig reliability was presumably improved. Genome DNA assembly was tested using 317 692 genomic reads of Arabitopsis thaliana (SRX105465). They were pre-processed with SeqTrimNext and assembled with CAP3 (http://www.scbi.uma.es/cap3) to obtain 35 777 contigs (Table 1, Genome columns). CoMiner detected 3259 contigs (9.1%) with one or more HERs, but only edited 1465 (4,1%), 74 of them being split into two or more contigs. Contigs before and after CoMiner treatment were mapped to A. thaliana genome using an in-house algorithm (H. Benzekri, unpublished results) to test the putative increase of contig reliability. A total of 35 412 (98.97%) and 35 458 (98.98%) contigs were mapped, respectively, providing a total of 3 362 691 and 3 362 923 mapped nucleotides, respectively (Table 1, Genome columns). When mapping was performed with a more restrictive mapper, such as Bowtie2 [8], 8453 original contigs (23.62%) and 8494 CoMiner-edited contigs (23.71%) were mapped. Edition slightly increased (1.0011.004 fold) the amount of mapped contigs and nucleotides in any case. Unfortunately, this increase is in the same range as the total contig number, indicating that more analyses are required to provide statistical significance for this weak increase. Even though CoMiner is currently only able to manage mismatches in overlaylayout-consensus assemblies, it seems to be a promising tool for automatic editing of mis-assembled reads. CoMiner performance was tested with transcriptome and genome data, and quality of edited contigs presumably seems improved. However, more real-world assemblies should be performed to give statistical significance to the qualitative results presented here. Finally, CoMiner edited contigs can always be Fig. 5. Example of a 1570 bp contig with several real HERs in red. The arrow is signaling a couple of HERs that was resolved by left-trimming one read and removing two other reads. Edition caused a gap within the contig that CoMiner resolved splitting it into two subcontigs, the left subcontig of 551 nt, and the right subcontig of 1018 nt. IWBBIO 2013. Proceedings Granada, 18-20 March, 2013 350 analyzed and visualized to search for more, different tentative errors by means of amosvalidate and Hawkeye, with the aim of spending less manual edition efforts. Acknowledgement The authors gratefully acknowledge Josep Planas (University of Barcelona, Spain) for kindly providing 454/FLX data from the AQUAGENET project (funded by InterregSudoe). We would also like to acknowledge the computer resources and technical support provided by the Plataforma Andaluza de Bioinformática of the University of Málaga, Spain. This study was supported by grants from the Spanish MICINN (BIO2009-07490) and Junta de Andalucía (P10-CVI-6075), as well as institutional funding to the research group BIO-114 and an agreement with AQUAGENET project. References 1. Istrail, S.; Sutton, G.G.; Florea, L.; Halpern, A.L.; Mobarry, C.M.; Lippert, R.; Walenz, B.; Shatkay, H.; Dew, I.; Miller, J.R.; Flanigan, M.J.; Edwards, N.J.; Bolanos, R.; Fasulo, D.; Halldorsson, B.V.; Hannenhalli, S.; Turner, R.; Yooseph, S.; Lu, F.; Nusskern, D.R.; Shue, B.C.; Zheng, X.H.; Zhong, F.; Delcher, A.L.; Huson, D.H.; Kravitz, S.A.; Mouchard, L.; Reinert, K.; Remington, K.A.; Clark, A.G.; Waterman, M.S.; Eichler, E.E.; Adams, M.D.; Hunkapiller, M.W.; Myers, E.W.; Venter, J.C. Whole-genome shotgun assembly and comparison of human genome assemblies. Proc Natl Acad Sci U S A 2004, 101, 1916-1921. 2. Phillippy, A.M.; Schatz, M.C.; Pop, M. Genome assembly forensics: finding the elusive mis-assembly. Genome Biol 2008, 9, R55. 3. Schatz, M.C.; Phillippy, A.M.; Shneiderman, B.; Salzberg, S.L. Hawkeye: an interactive visual analytics tool for genome assemblies. Genome Biol 2007, 8, R34. 4. Bonfield, J.K.; Whitwham, A. Gap5--editing the billion fragment sequence assembly. Bioinformatics 2010, 26, 1699-1703. 5. Gordon, D.; Abajian, C.; Green, P. Consed: a graphical tool for sequence finishing. Genome Res 1998, 8, 195-202. 6. Guerrero, D.; Bautista, R.; Villalobos, D.P.; Canton, F.R.; Claros, M.G. AlignMiner: a Web-based tool for detection of divergent regions in multiple sequence alignments of conserved sequences. Algorithms Mol Biol 2010, 5, 24. 7. Falgueras, J.; Lara, A.J.; Fernandez-Pozo, N.; Canton, F.R.; Perez-Trabado, G.; Claros, M.G. SeqTrim: a high-throughput pipeline for pre-processing any type of sequence read. BMC Bioinformatics 2010, 11, 38. 8. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2. Nat Methods 2012, 9, 357-359. IWBBIO 2013. Proceedings Granada, 18-20 March, 2013 351