Aportación Bioinformática a la Biología de Especies Marinas

AUTOR: Hicham Benzekri http://orcid.org/0000-0003-1064-5075 EDITA: Publicaciones y Divulgación Científica. Universidad de Málaga Esta obra está bajo una licencia de Creative Commons Reconocimiento-NoComercialSinObraDerivada 4.0 Internacional: http://creativecommons.org/licenses/by-nc-nd/4.0/legalcode 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 APORTACÓN BIOINFORMÁTICA A LA BIOLOGÍA DE ESPECIES MARINAS HICHAM BENZEKRI UNIVERSIDAD DE MÁLAGA Departamento de Biología Molecular y Bioquímica Facultad de Ciencias y Plataforma Andaluza de Bioinformática Edificio de Bioinnovación Málaga. junio de 2016 APORTACIÓN BIOINFORMÁTICA A LA BIOLOGÍA DE ESPECIES MARINAS Memoria presentada por: Hicham Benzekri Para optar al grado de Doctor por la Universidad de Málaga Tesis realizada bajo la dirección del Dr. M. Gonzalo Claros Díaz en la Plataforma Andaluza de Bioinformática y el Departamento de Biología Molecular y Bioquímica de la Universidad de Málaga Fdo. Hicham Benzekri Vº.Bº. DIRECTOR DE LA TESIS DOCTORAL: Fdo. M. Gonzalo Claros Díaz Málaga, 11 de noviembre de 2015 D. M. Gonzalo Claros Díaz, Investigador de la Plataforma Andaluza de Bioinformática y el Departamento de Biología Molecular y Bioquímica de la Universidad de Málaga. CERTIFICA: Que Don Hicham Benzekri, Ingeniero agronomo, ha realizado bajo mi dirección en el Departamento de Biología Molecular y Bioquímica y en la Plataforma Andaluza de Bioinformática de la Universidad de Málaga, el trabajo de investigación recogido en la presente memoria de Tesis Doctoral que lleva por título: “Aportación bioinformática a la biología de especies marinas”. Tras la revisión de la presente Memoria se ha estimado oportuna su presentación ante la Comisión de Evaluación correspondiente, por lo que autorizo su exposición 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álaga, 11 de noviembre de 2015 El Director de la Tesis, Dr. D. M. Gonzalo Claros Díaz 5 AGRADECIMIENTOS Agradecimientos En primer lugar quiero agradecer de manera especial y sincera al Dr. M. Gonzalo Claros Díaz por aceptarme para realizar esta tesis bajo su dirección y por haberme facilitado los medios para poder llevar a cabo este trabajo de investigación. Le agradezco sus atentas y rápidas respuestas a las diferentes dudas e inquietudes surgidas durante el desarrollo de este trabajo, su paciencia, confianza, y todo el apoyo que me ha ofrecido. Quiero expresar también mi más sincero agradecimiento al Dr. Manuel Manchado Campaña, del grupo de investigación de IFAPA El Toruño (Cádiz), por haber financiado el trabajo de tesis, por su apoyo y confianza, y por haber sido tan atento y dispuesto a resolver todas mis dudas durante el periodo de investigación. Agradezco al Dr. Oswaldo Trelles Salazar por haberme acogido en su departamento durante un periodo de esta tesis, y a los miembros de su grupo de investigación: Nono, Javier, Maxi, Johan, Alfredo, José Manuel y Vicki, gracias por serme tan útiles en el trabajo y por todos esos momentos agradables que pasamos juntos. A los compañeros de la plataforma de bioinformática: Rocío, Darío, Rafa, Noé, Rosario, Isa, Pedro y David, muchas gracias por estar presentes y ser tan atentos y amistosos. Nunca olvidaré esos momentos que pasé en este ambiente tan agradable y lleno de aprendizaje y de buen humor. A Rocío gracias por iniciarme en este mundo de la bioinformática y por tu inestimable ayuda y apoyo durante todos estos años. A Noé, gracias por haberme introducido en el trabajo de la tesis y por todo los que me has enseñado; a Darío gracias por ser tan amable ofreciéndome tu ayuda en cuanto la necesitaba y por todo lo que me has enseñado; a Rafa gracias por tu cordialidad y tu ayuda. A los demás compañeros del Edificio de bioinnovación: Pepi, Diego, Carlos, Xhanti y Rocío gracias por los ratos agradables que pasamos juntos durante las comidas y por vuestra simpatía y afecto. A mis grandes amigos que conocí en España: Said, Issam, Victor y su novia Luli, y Paco, les agradezco muchísimo su apoyo y el estar tan presentes en mis momentos difíciles. A mi familia, mis hermanos, mis cuñados, mis sobrinos y en especial a mis padres, ya que sin ellos nada de esto habría sido posible. A todos los que confiaron en mí y me apoyaron. 7 ACRÓNIMOS Y ABREVIATURAS Acrónimos y Abreviaturas ADN: ácido desoxirribonucleico ARN: ácido ribonucleico ADNc: copia en ácido ribonucleico ARNm: ácido ribonucleico mensajero ARNnc: ARN no codificante ARNr: ácido ribonucleico ribosómico ARNt: ARN de transferencia API: interfaz de programación de aplicaciones BAC: cromosoma artificial bacteriano BLAST: Basic Local Alignment Search Tool CPU: Unidad central de procesamiento CEG : Core Eukaryotic Genes dATP: trifosfato de desoxiadenosina dCTP: trifosfato de desoxicitidina ddNTP: didesoxinucleótidos trifosfato DDBJ: Banco de Datos de ADN de Japón dGTP: trifosfato de desoxiguanosina dTTP: trifosfato de desoxitimidina EC: Enzyme Commision EBI: European Bioinformatics Institute EMBL: Laboratorio Europeo de Biología Molecular ENA: European Nucleotide Archive EST: etiquetas de secuencias expresadas (del inglés Expressed Sequence Tag ) GB: gigabyte GO: Gene ontology HSP: High-scoring Segment Pairs HER: High-entropy region ID: identificador INDEL: inserciones y deleciones IPR: Interpro JSON: formato ligero para el intercambio de datos (del inglés JavaScript object notation ) KEGG: Kyoto Encyclopedia of Genes and Genomes MID: Molecular Identifier miRNA: microARN MP: lecturas pareadas de tipo mate-pair NCBI: National Center for Biotechnology Information NGS: secuenciación de nueva generación (del inglés next generation sequencing ) NHGRI: National Human Genome Research Institute Nt: nuleótido OLC: overlap-layout-consensus ORF: Marco abierto de lectura, (del inglés open reading frame ) pb: pares de bases 14 INDICE GENERAL III.3.1.2. Array-Jobs ..................................................................................................... 56 III.3.1.3. Gemas de Ruby ............................................................................................. 58 III.3.1.4. Aptana Studio ..................................................... ¡Error! Marcador no definido. III.3.2. Para el ensamblaje de secuencias ..................................................................... 60 III.3.2.1. CAP3 ............................................................................................................. 60 III.3.2.2. Mira ............................................................................................................... 61 III.3.2.3. EULER-SR .................................................................................................... 62 III.3.2.4. Velvet ............................................................................................................ 62 III.3.2.5. Oases ............................................................................................................. 63 III.3.2.6. SOAPdenovo-trans ........................................................................................ 64 III.3.2.7. CABOG ......................................................................................................... 65 III.3.2.8. Ray ................................................................................................................ 66 III.3.3. Para el tratamiento y mejora de los ensamblajes ............................................. 67 III.3.3.1. GAM-NGS .................................................................................................... 67 III.3.3.2. SOAPdenovo Scaffolder ................................................................................ 69 III.3.3.3. SSPACE ........................................................................................................ 70 III.3.3.4. SOAPdenovo GapCloser ............................................................................... 70 III.3.4. Para analizar secuencias .................................................................................. 71 III.3.4.1. SeqTrimNext ................................................................................................. 71 III.3.4.2. Full-LengtherNext ......................................................................................... 72 III.3.4.3. AutoFact ....................................................................................................... 72 III.3.4.4. Basic Local Alignment Search Tool (BLAST) .............................................. 73 III.3.4.5. Sma3s ............................................................................................................ 74 III.3.4.6. RAST ............................................................................................................ 75 III.3.4.7. Bowtie2 ......................................................................................................... 75 III.3.4.8. Samtools ........................................................................................................ 76 III.3.4.9. CD-HIT ......................................................................................................... 77 III.3.4.10. MREPS ....................................................................................................... 77 III.3.5.11. Gigabayes .................................................................................................... 78 III.3.4.12. Tablet .......................................................................................................... 79 III.3.4.13. GEvo ........................................................................................................... 79 III.3.4.14. NUCMER .................................................................................................... 79 III.3.5. Scripts escritos ................................................................................................. 79 III.3.5.1. Para generar lecturas artificiales ................................................................... 80 15 INDICE GENERAL III.3.5.2. Para combinar dos ficheros ACE .................................................................. 80 III.3.5.3. Para calcular el tamaño de inserto real en lecturas pareadas. ...................... 80 III.3.5.4. Para eliminar ensamblajes erróneos .............................................................. 81 III.3.5.5. Para seleccionar los transcritos que representan las sondas de microarray. .. 81 III.3.5.6. Para calcular las estadísticas sobre los marcadores SSR ............................... 81 III.3.5.7. Para solapar los contigs genómicos ............................................................... 82 III.3.5.8. Para la corrección automática de los ensamblajes (Cominer) ....................... 82 III.3.5.9. Para ordenamiento de los contigs o scaffolds y acabado de los genomas (ICMapper) ................................................................................................................. 83 III.4. Datos biológicos ...................................................................................................... 83 III.4.1. Secuencias transcriptómicas ............................................................................. 83 III.4.1.1. Para el transcriptoma de Solea senegalensis ................................................. 83 III.4.1.2. Para el transcriptoma de Solea solea ............................................................ 84 III.4.1.3. Para el transcriptoma de Tisochrysis lutea ................................................... 84 III.4.2. Secuencias genómicas ....................................................................................... 85 III.4.2.1. Secuencias genómicas de Photobacterium damselae subsp. piscicida ............ 85 III.4.2.2. Secuencias genómicas de Solea senegalensis .................................................. 85 IV. Resultados y discusión ................................................................................................. 89 IV.1. Algoritmos preparatorios. ....................................................................................... 89 IV.1.1. Definición de nuevas especies contaminantes ................................................... 89 IV.1.2. Selección de las condiciones de ensamblaje de novo ......................................... 90 IV.1.2.1. Ensamblaje de transcriptomas ...................................................................... 90 IV.1.2.2. Ensamblaje de genomas ................................................................................ 94 IV.1.2.3. Combinación de ensamblajes de varios programas ....................................... 99 IV.1.2.4. Tamaño real del inserto en las lecturas pareadas ....................................... 101 IV.1.3. CoMiner: validación de ensamblajes ............................................................... 102 IV.1.4. Ordenamiento de los scaffolds y acabado de los genomas ............................... 108 IV.2. Ensamblajes de transcriptomas .............................................................................. 115 IV.2.1. Estrategia tomada como modelo ..................................................................... 115 IV.2.2. Transcriptomas de Solea senegalensis y Solea solea ........................................ 119 IV.2.2.1. Preprocesamiento de las librerías ................................................................ 119 IV.2.2.2. Ensamblaje ................................................................................................. 120 IV.2.2.3. Anotación de los transcriptomas. ................................................................ 125 IV.2.2.4. SoleaDB, una base de datos para explorar los transcriptomas de Solea ..... 126 16 INDICE GENERAL IV.2.2.5. Calidad del transcriptoma de los dos lenguados .......................................... 130 IV.2.2.6. Análisis comparativo entre las dos especies de lenguado ............................. 133 IV.2.2.7. Los transcriptomas de lenguado como fuente de marcadores moleculares .. 138 IV.2.2.8. Selección de transcritos para un estudio de microarray de oligonucleótidos a realizar sobre S. senegalensis ..................................................................................... 141 IV.2.3. Transcriptoma de la microalga Tisochrysis lutea ............................................ 142 IV.2.3.1. Preprocesamiento de las librerías ................................................................ 142 IV.2.3.2. Ensamblaje .................................................................................................. 144 IV.2.3.3. Anotación del transcriptoma e inclusión en base de datos ......................... 145 IV.2.3.4. Análisis del transcriptoma de Tisochrysis lutea .......................................... 146 IV.2.4. Transcriptoma de Ruditapes decussatus (almeja fina) .................................... 152 IV.2.4.1. Preprocesamiento de las librerías ................................................................ 152 IV.2.4.2. Ensamblaje .................................................................................................. 153 IV.2.4.3. Anotación e inclusión en base de datos ....................................................... 154 IV.2.4.4. Análisis del transcriptoma de R. decussatus ............................................... 155 IV.3. Ensamblajes de genoma ......................................................................................... 156 IV.3.1Genoma de Photobacterium damselae subsp. piscicida ................................... 156 IV.3.1.1. Preprocesamiento de las lecturas brutas ..................................................... 157 IV.3.1.2. Ensamblaje .................................................................................................. 159 IV.3.1.3. Mejora del ensamblaje ................................................................................. 159 IV.3.1.4. Anotación de los dos borradores de genomas .............................................. 161 IV.3.1.5. Análisis comparativo entre las dos cepas de P. damselae subsp. piscicida .. 164 IV.3.1.6. Descubrimiento de los plásmidos de L091106-03H ...................................... 168 IV.3.2. Genoma del lenguado senegalés ( Solea senegalensis ). ...................................... 170 IV.3.2.1. Preprocesamiento de las lecturas genómicas brutas .................................... 170 IV.3.2.2. Ensamblaje .................................................................................................. 174 IV.3.2.3. Evaluación del genoma ensamblado ............................................................ 177 IV.3.2.3.1. Prueba de compleción de los genes en el ensamblaje ........................ 177 IV.3.2.3.2. Confirmación de la sintenia entre los lenguados ............................... 178 IV.3.2.4. Los super-scaffolds como posibles cromosomas del lenguado senegalés ....... 180 IV.3.2.5. Validación de los super-scaffolds ................................................................. 183 IV.3.2.5.1. Localización de los transcritos de lenguado. ..................................... 183 IV.3.2.5.2.Validación de los marcadores moleculares del lenguado senegalés ...... 184 V. Conclusiones ................................................................................................................ 199 17 INDICE GENERAL VI. Bibliografía ................................................................................................................ 203 VII. Apéndices ................................................................................................................. 223 Apendice A: Script para crear lecturas artificiales .......................................................... 223 Apendice B: Script para calcular el tamaño de inserto en lecturas Illumina a partir de un fichero de mapeo (SAM) ................................................................................................. 226 Apendice C: Script para calcular estadísticas sobre marcadores SSR a partir del fichero de salida MREPS ................................................................................................................ 228 Apendice D: Software utilizado en las diferentes etapas del análisis de los transcriptomas y genomas estudiados ........................................................................................................ 235 INDICE GENERAL 18 19 INTRODUCCIÓN Parte I Introducción 20 INTRODUCCIÓN 21 INTRODUCCIÓN I. Introducción I.1. La importancia de la secuenciación a día de hoy La secuenciación del ADN es el proceso de determinación del orden preciso de nucleótidos en una molécula de ADN. Incluye a todas los métodos y tecnologías utilizadas para determinar el orden de cuatro bases: adenina, guanina, citosina y tiamina. Determinar la secuencia de ADN es útil en el estudio de la investigación básica de los procesos biológicos fundamentales. La genética molecular ha acelerado significativamente la investigación y los descubrimientos en biología y causaron una revolución científica en varios campos:  En el campo de la investigación médica, el diagnostico de muchas enfermedades es mucho más eficaz. Se puede evaluar la susceptibilidad de las personas a enfermedades específicas y suministrar los medicamentos más adecuados contra estas enfermedades.  En el campo legal, las pruebas de paternidad y forenses son más precisas con las tecnologías de la PCR y secuenciación.  En ganadería y agricultura está facilitando la mejora de los animales y vegetales en relación con su capacidad productiva o resistencia a las enfermedades. Las técnicas de secuenciación de ADN han evolucionado mucho a lo largo de los últimos años y han generado un gran impacto sobre la investigación científica aplicada a la biología y a la medicina, donde se ha pasado de secuenciar los genes uno a uno de modo manual en los años 80 [1, 2] a proyectos de secuenciación de mil o más genomas en la actualidad [3-5]. En los siguientes apartados, vamos a ver cómo se ha pasado de la secuenciación manual de varios cientos de nucleótidos a la secuenciación automática con capilares y posteriormente al desarrollo de la secuenciación de nueva generación (NGS, del inglés Next Generation Sequencing). Como la NGS genera secuencias de 75 a 800 nt de uno o varios millones de moléculas a la vez en reacciones a escala nanotecnológica, producen un volumen de secuencias tal que su procesamiento necesita la ayuda de la bioinformática. La secuenciación de los organismos representa la primera etapa del proceso de obtención de sus transcriptomas o genomas anotados (figura I.1) 22 INTRODUCCIÓN Figura I.1: Esquema general del flujo de trabajo para obtener genomas o transcriptomas anotados I.1.1. Secuenciación de tipo Sanger, hoy considerado “clásico”. En 1977 se publicaron a la vez el método de secuenciación química de Maxan y Gilbert [6] y el método de secuenciación de terminación de cadena de Sanger [7], también conocido como método enzimático o secuenciación por didesoxinucleótidos. Este método, esquematizado en la figura I.2, se basa en la formación de un extremo 3’ terminador mediante la incorporación de didesoxinucleótidos trifosfato (ddNTP) por lo que la ADN polimerasa no puede seguir extendiendo al carecer de grupo hidroxilo en 3’. La ADN polimerasa sintetiza hebras complementarias a la que se quiere secuenciar, por lo tanto necesita un cebador, un oligonucleótido diseñado para que se hibride con el extremo 3’ de ésta, además de una mezcla de dNTP (uno de ellos radiactivo) y ddNTP. Para llevar a cabo la secuenciación se realizan cuatro reacciones por separado, cada una con un ddNTP diferente. De este modo se obtienen secuencias truncadas en diferentes posiciones de 3’ del nucleótido correspondiente a cada tubo. Estas moléculas se separaran por tamaño en una electroforesis, utilizando un carril para cada tubo, para visualizar el patrón de bandas del que se puede deducir la secuencia [7]. Frederick Sanger fue premiado con su segundo Nobel de química en 1980 por su contribución en la determinación de las secuencias de los ácidos nucleicos, compartiendo el premio con Walter Gilbert y Paul Berg [8], lo que puede dar una idea de lo que su técnica Secuenciaci o n Preprocesamiento Ensamblaje Sanger Roche/454 Illumina Anotaci o n Algoritmos voraces Algoritmos OLC Algoritmos basados en De Bruijn 23 INTRODUCCIÓN aportó a la investigación científica, posibilitando la secuenciación, entre otros, del genoma humano [9, 10]. Figura I.2: (tomada de http://mol-biol4masters.masters.grkraj.org) Método de secuenciación por didesoxinucleótidos de Sanger El método de Sanger se adaptó posteriormente para la secuenciación automática en capilares, sustituyendo el marcaje con isótopos radiactivos por un marcaje con 4 fluorocromos diferentes. Con ello se consiguió realizar la reacción en un único tubo y migrarla en un solo capilar (figura I.3-a), años más tarde las máquinas pasaron a incorporar varios capilares (de 4 a 384) para paralelizar este proceso. Después de tres décadas de continuo perfeccionamiento, es posible obtener secuencias de unas 1000 pb con una fiabilidad del 99,999 % y con un coste de unos 0,5 $ por kilobase (Figura I.4). 30 INTRODUCCIÓN Figura I.8: (tomada de [34]) Comparación entre las técnicas de secuenciación (a) extremos emparejados (paired-end) y (b) parejas conjugadas (mate-pair) I.2. Reconstrucción de fragmentos de secuencias I.2.1. Preprocesamiento. Después de obtener las lecturas de una secuenciación es indispensable preprocesarlas para eliminar los elementos que provienen de la manipulación del ADN o del ARN en el laboratorio. Así se obtendrá únicamente la parte útil y se descartan los fragmentos de baja calidad, contaminantes, vectores, adaptadores [35], y otros artefactos. En las nuevas tecnologías, además, es posible encontrar otros elementos nuevos que se deben considerar en el preprocesamiento. Por ejemplo, el etiquetado con MID (Roche) u otras etiquetas implica que hay que retirarlas de la secuencia final, y además hay que usarlas para agrupar las lecturas por experimentos. En resumen, la calidad y la fiabilidad del posterior ensamblaje 31 INTRODUCCIÓN dependerá de un buen preprocesamiento [36] ya que, de otra manera, los consensos obtenidos podrían contener numerosos errores. Ya existen programas para preprocesar secuencias de tipo Sanger, como SeqClean ( http://compbio.dfci.harvard.edu/tgi/software /), Lucy [37], ESTPrep [38] y SeqTrim [36]. Algunas herramientas bioinformáticas como TagCleaner [39] y NovoBarCode (Novocraft) se han diseñado para identificar los MID y otras etiquetas, y clasificar las secuencias procedentes de los diferentes experimentos. En Galaxy [40] ( https://usegalaxy.org ) se puede acceder a distintas herramientas, como por ejemplo FastQC, Groomer, Splitter, etc., para montar un flujo de trabajo de limpieza. En cuanto a los contaminantes, los más comunes en el laboratorio suelen ser hongos y bacterias, sobre todo procedentes de la rizosfera [41]. También se consideran contaminantes los restos de tejidos humanos procedentes de los investigadores y microorganismos utilizados con frecuencia en el laboratorio, como E.coli [42] y Agrobacterium tumefaciens [43], este último sobre todo cuando se trabaja con plantas. Cuando se trata de la secuenciación organismos animales, también se consideran contaminantes los organismos que se utilizan como alimentación y/o las especies de bacterias que normalmente infectan al organismo estudiado. Si no se eliminan las secuencias procedentes de organismos contaminantes se corre el riesgo de considerarlas parte del organismo que se está estudiando [44]. Para la detección y eliminación de secuencias contaminantes existen programas como DeconSeq [45]. Ninguno de los programas anteriores se basta por sí solo para un buen preprocesamiento. Por eso, para la NGS se desarrollaron nuevas herramientas bioinformáticas especializadas como SeqTrim-Next [46] o el paquete Galaxy [40]. Después del preprocesamiento de las lecturas viene la etapa de su ensamblaje en secuencias más largas. Existen varios algoritmos de ensamblaje: I.2.2. Conceptos básicos sobre ensamblaje de secuencias Un ensamblaje (o montaje) de secuencias se refiere al alineamiento y mezcla de multiples fragmentos de una secuencia de ADN mucho mayor para reconstruir la secuencia original. En el ensamblaje las lecturas se agrupan en contigs y los contigs en scaffolds. Los contigs representan el alineamiento multiple de lecturas (o fragmentos de lecturas). Los scaffolds (a veces llamados supercontigs o metacontigs) definen el orden y la orientación de los contigs y las longitudes de los huecos entre los contigs (figura I.9). 32 INTRODUCCIÓN Figura I.9: (tomada de [47]) Estrategia de secuenciación utilizando secuencias pareadas. A partir de las lecturas generadas en un proyecto de secuenciación de cualquier tecnología o estrategia, se obtienen las secuencias consenso de los contigs mediante el ensamblaje con un programa informático. Las secuencias de los contigs se organizan entonces formando un scaffold basándose en la información de las secuencias pareadas para ordenar contigs no solapantes. Finalmente los scaffold pueden ordenarse igual que el genoma identificando en ellos marcadores que se conozcan en el genoma, como por ejemplo STS, marcadores moleculares y genes (círculos rojos) Existen dos tipos principales de ensamblaje: Ensamblaje comparativo (por mapeo) Basándose en un genoma secuenciado previamente y que suponemos sea similar al que se quiere ensamblar. El procedimiento básico tratará de colocar cada una de las lecturas en la posición adecuada utilizando el genoma de referencia como guía. Ensamblaje de novo Intenta reconstruir la secuencia de ADN completa a partir de las lecturas sin ningún tipo de conocimiento previo acerca del genoma a ensamblar. Busca lecturas cuyo final coincida con el principio de otra de forma que se puedan unir para formar fragmentos mayores hasta completar el genoma. I.2.3. Algoritmos voraces ( greedy ) Los primeros ensambladores de NGS utilizaban algoritmos voraces [48, 49] que aplican heurísticas en cada fase del algoritmo mediante la cuales pretenden buscar soluciones parciales óptimas y una operación básica: dada una lectura o contig, añade una lectura o 33 INTRODUCCIÓN contig más. Esta operación se repite hasta que no sea posible. Cada operación utiliza el solapamiento con más puntuación para generar la siguiente unión. La función de puntuación mide, entre otros indicadores, el número de bases coincidentes en el solapamiento. De esta manera, los contigs crecen por extensión tomando la lectura que se encuentra por el siguiente solapamiento con más puntuación. Los algoritmos voraces pueden atascarse cuando un contig puntual toma lecturas que deben ser de otros contigs. Los algoritmos voraces son implícitamente algoritmos de grafos ya que simplifican drásticamente el grafo al considerar sólo los enlaces con más puntuación, como optimización pueden instanciar sólo un solapamiento por cada lectura que examinan y pueden también descartar cada solapamiento inmediatamente después de la extensión del contig. Al igual que todos los ensambladores, los algoritmos voraces necesitan mecanismos para evitar incorporar solapamientos falsos en los contigs, un ensamblador que se base en solapamientos falsos unirá secuencias no relacionadas a cualquier lado de una repetición. El primer ensamblador descrito para las lecturas cortas fue SSAKE [50] que se diseñó para lecturas no pareadas de longitud uniforme. SHARCGS [51] también trabaja con longitudes uniformes, alta cobertura y lecturas cortas no pareadas, pero añade una funcionalidad de prey posprocesado a SSAKE. El preprocesado filtra lecturas erróneas mediante comprobación de un número mínimo de coincidencias exactas de máxima longitud en otras lecturas. VCAKE [52] es otro algoritmo de extensión iterativo distinto a SSAKE y a SHARCGS, puede incorporar coincidencias imperfectas durante la ampliación de los contigs. VCAKE se utilizó en combinación con NEWBLER [12] en una integración para datos híbridos de Solexa y 454 [53]. Otra integración combinó NEWBLER y el CELERA ASSEMBLER [10, 54] para datos híbridos de 454 y Sanger [55]. Ambas integraciones rompen los contigs del primer ensamblador para producir pseudolecturas adecuadas para el segundo ensamblador. Esta última integración ajusta la cobertura de lectura y los indicadores de calidad en las pseudolecturas que genera, lo cual ayuda al segundo ensamblador a dar un peso a los contigs con gran cobertura desde el primer ensamblaje. I.2.4. Algoritmos de solapamiento-composición-consenso (OLC). La estrategia OLC (del inglés overlap-layout-consensus ) fue ampliamente utilizada en los ensambladores para datos Sanger y fue optimizada para genomas grandes en diversos tipos de software incluyendo el CELERA ASSEMBLER [10, 54], ARACHNE [56] y CAP y PCAP [57], EDENA, CABOG, TIGR, ATLAS, PHRAP, PHUSION y SHORTY, ha sido estudiada ampliamente [58]. Los ensambladores que utilizan esta estrategia están muy orientados a ensamblajes de novo , utilizan un grafo de solapamientos y operan ejecutando tres fases: (i). Búsqueda de solapamientos mediante la comparación de lecturas pareadas todas contra todas. El software precalcula los k-meros contenidos en todas las lecturas, selecciona 34 INTRODUCCIÓN los solapamientos candidatos que comparten dichos K-meros y calcula las alineaciones utilizando un K-mero como semilla. El descubrimiento de solapamientos es sensible al tamaño de los K-meros, a la longitud mínima de solapamiento y al porcentaje mínimo de identidad requerido en dichos solapamientos. El comportamiento de estos parámetros se ve afectado por los errores de secuenciación y una baja cobertura. Valores mayores conducen a una mayor precisión pero con contigs más cortos. (ii). La construcción y manipulación de un grafo de solapamientos conduce a un diseño del grafo basado en un conjunto de lecturas representativo, es decir, el grafo no necesita incluir todas las secuencias base, por lo que estos grafos pueden utilizarse en grandes genomas ya que pueden ajustar el tamaño del grafo a la cantidad de memoria que utilizan. (iii). Generación de la secuencia consenso mediante un alineamiento múltiple de lecturas, aunque no haya un método eficiente para calcularla [59]. Esta fase se puede ejecutar en paralelo, por contigs. Existen dos ensambladores que aplican la estrategia OLC a lecturas cortas de las plataformas Illumina y SOLiD: el software EDENA [60], que fue desarrollado para lecturas no pareadas de longitud uniforme, descarta lecturas duplicadas y busca todos los solapamientos perfectos libres de errores. Elimina solapamientos individuales que son redundantes con otros solapamientos mediante una aplicación del algoritmo de reducción de solapamientos transitivos [61]. El otro software es SHORTY [62] que trata el caso especial donde unas pocas lecturas pueden actuar como semillas para conseguir lecturas cortas y sus extremos emparejados, y mediante iteraciones utiliza contigs como semillas para generar nuevos contigs. I.2.5. Algoritmos basados en grafos de De Bruijn Esta tercera aproximación para ensamblaje de secuencias se utiliza principalmente para lecturas cortas de las plataformas de SOLiD e Illumina. Se basa en grafos de frecuencias de K-meros para tratar grandes cantidades de lecturas cortas, ya que los grafos de K-meros no requieren el análisis de todos los solapamientos mediante la comprobación de todos contra todos, ni tampoco es necesario almacenar las lecturas individuales en sus solapamientos, ni tampoco comprime las secuencias repetitivas. Por el contrario, los grafos de K-meros no mantienen secuencias originales y son grandes consumidores de memoria para grandes genomas, aunque con los sistemas de memoria distribuida se comportan bastante bien [63]. Al final, el grafo de De Bruijn hay que resolverlo con una estrategia euleriana (figura I.10) que se basa en las siguientes etapas: (i). Fragmentación de las lecturas secuenciadas en una colección de oligomeros (K-meros), donde todas las lecturas tienen la misma longitud (k). Los cebadores redundantes se comprimen en uno solo para reducir la complejidad del análisis, aunque se conserva la información del número de veces que se repiten. Los valores de k suelen asignarse entre 19 y la longitud de las lecturas; en la práctica, estos valores suelen escogerse basándose en un 35 INTRODUCCIÓN experimento previo con datos similares [64]. Esta fragmentación de lecturas evita el paso limitante de tener que comparar todas las secuencias entre sí. Figura I.10: (tomada de http://www.homolog.us ) División de una lectura en k-meros y construcción de un grafo de De Bruijn con el conjunto de k-meros para reconstruir la lectura original (ii). Construcción de un grafo de De Bruijn con el conjunto de k-meros, de modo que cada nodo contiene un k-mero de longitud k – 1 que se solapa exactamente en k – 2 nucleótidos con otros nodos. (iii). Búsqueda de los caminos eulerianos que reconstruyen la secuencia original de la que proceden todas las lecturas relacionadas. Esta estrategia se ve muy afectada por las repeticiones y por los errores de secuenciación que los ensamblajes por solapamiento [63]. Además, como el ADN es de doble cadena, a la hora de reconstruir los caminos eulerianos puede darse el caso que algún k-mero aparezca en el sentido de la transcripción y otro en antisentido, provocando secuencias consenso artefactuales. Algunos ejemplos de los ensambladores más populares que utilizan esta estrategia son Euler-SR [65], Velvet [11], SOAPdenovo [66] y ABySS [67] y el Ray [68] (este último permite la paralelización del ensamblaje). I.2.6. Estrategias para la secuenciación de transcriptomas. La secuenciación de pools de ADNc a menudo se utiliza para caracterizar el transcriptoma de un organismo de un modo rápido y barato. El transcriptoma de un organismo engloba al conjunto de genes que se expresan (EST) en una célula o conjunto de células, que incluyen tanto a los ARN que codifican proteínas como los no codificantes (ARNnc). Este tipo de secuenciación proporciona información sobre los genes de un organismo a bajo coste en comparación con la secuenciación de genomas, ya que solo se 36 INTRODUCCIÓN investigan las regiones que se transcriben, en lugar del genoma completo. La generación de EST a partir de ARNm se considera la estrategia más frecuente y más útil para descubrir genes [69]. En los casos en los que el análisis transcriptómico está enfocado a estudiar los genes que codifican proteínas, es importante utilizar secuencias de ARN enriquecidas en poli - (A)+, ya que de este modo se reduce en gran parte los ARN no deseados, como los ARN pequeños y los abundantes ARN ribosómicos (ARNr). También es deseable realizar la normalización de las muestras de ARN, ya que así se reduce la aparición de los transcritos más abundantes y se aumenta la de los poco abundantes [70]. I.2.7. Estrategias para la secuenciación de genomas En la actualidad se utilizan principalmente dos estrategias para abordar la secuenciación de genomas, BAC a BAC (cromosomas artificiales de bacterias), o mediante la secuenciación de todo el genoma (WGS, del inglés whole genome sequencing). La elección de la estrategia más apropiada dependerá del tamaño y de la complejidad del genoma que se quiere secuenciar, así como de las tecnologías utilizadas y del presupuesto disponible. Otra opción más barata y rápida consiste en la secuenciación del transcriptoma para el estudio de los genes que se expresan (EST, etiquetas de secuencias expresadas, del inglés expressed sequence tags) en determinadas condiciones experimentales, útil para obtener información en especies no modelo con genomas complejos. BAC a BAC: El genoma o un cromosoma que se quiere secuenciar se trocea en fragmentos solapantes de menor tamaño (100 a 200 kpb aproximadamente) que se clonan en vectores BAC (figura I.11-A.1). Después se seleccionan los BAC y se ordenan en un mapa físico (figura I.11-A.2). Al final, cada BAC se secuencia por separado tras una digestión parcial con una enzima de restricción que producirá fragmentos solapantes al azar, en lo que se conoce por el término inglés shotgun (figura I.11-A.3). Al secuenciar los BAC de uno en uno se consigue reducir la complejidad del ensamblaje, de modo que esta estrategia puede resultar útil para genomas con gran cantidad de repeticiones. Sin embargo, crear la genoteca de BAC y mapear cada clon sobre el genoma supone una enorme cantidad de esfuerzo, además de conllevar mucho tiempo y un gran coste económico. Por eso, esta técnica se está utilizando cada vez menos (figura I.12). WGS: En esta estrategia se evita la generación de genotecas BAC fragmentando directamente el genoma al azar en elementos solapantes de menor tamaño (figura I.11-B.1). Los fragmentos se secuencian directamente con tecnologías de NGS (figura I.11-B.2), o se clonan primero y se secuencian después en el caso de las tecnologías de secuenciación basadas en capilares. Debido a su mayor sencillez, esta es la estrategia de secuenciación de genomas más utilizada hoy en día (figura I.12). A diferencia del caso anterior, aquí toda la complejidad de la secuencia del genoma tendrá que resolverse con aplicaciones bioinformáticas, ya que no se dispone de ninguna información orientadora. 37 INTRODUCCIÓN Figura I.11: Comparación de las estrategias de secuenciación (A) BAC a BAC, y (B) WGS Figura I.12: (tomada de [71]) Número de genomas secuenciados cada año desde 1995. (a) Número de genomas secuenciados utilizando WGS y con otras estrategias. (b) Cada columna indica el tamaño acumulado (en miles de millones de pares de bases) de los genomas secuenciados por WGS (columna de la izquierda) y por otros métodos (columna de la derecha) 38 INTRODUCCIÓN I.3. Anotación de los ensamblajes I.3.1. Se anota por similitud Una vez se han ensamblado los unigenes de un transcriptoma o se han reconstruido las regiones genómicas, se procede a anotar las secuencias ya que el interés de tener las secuencias de un transcriptoma o de un genoma reside en la obtención de información biológica que permita seguir profundizando en el funcionamiento de los organismos y su relación con el entorno. Este incremento de conocimiento científico podrá posteriormente traducirse en aplicaciones prácticas que conlleven beneficios económicos y mejoras en la calidad de vida. La anotación de genes se lleva a cabo habitualmente de un modo automático e informatizado. Existen muchas herramientas bioinformáticas libres, como Blast [72], que resulta de gran utilidad para encontrar parecidos significativos con genes conocidos que están almacenados en las bases de datos. La anotación por similitud de secuencia tiene algunas limitaciones, en especial en los organismos no modelo, donde el número de genes con anotaciones es menor, sobre todo porque se basan principalmente en la información que se conoce de especies modelo como Arabidopsis thaliana [31] y el pez cebra ( Danio rerio ). La anotación basada en similitud será tan exacta como la anotación de las secuencias almacenadas en las bases de datos con las que se compara. Para comparar las secuencias de nucleótidos de los unigenes obtenidos durante el ensamblaje de un transcriptoma suele emplearse con frecuencia Blastx con bases de datos que contienen proteínas conocidas, como RefSeq_protein de GenBank [73] o Swiss-Prot y TrEMBL, ambas de UniProt [74]. Estas dos bases de datos de secuencias de proteínas, tienen diferentes niveles de anotación. Por ejemplo, Swiss-Prot está revisada manualmente, y TrEMBL, utiliza anotaciones automáticas procedentes de referencias cruzadas [74]. Suele ser frecuente, cuando se trabaja con especies no modelo, encontrarse una gran cantidad de ortólogos de función desconocida. Además hay que tener en cuenta que si no se toman las precauciones adecuadas que consisten en la eliminación de los contaminantes en las lecturas originales y de los contigs quimericos, a veces se incorporan a las bases de datos secuencias con errores de anotación o secuencias que realmente son contaminantes o artefactos [44]. Si se considera como buena una secuencia mal anotada, este error también se mantendrá en nuestras secuencias anotadas por similitud. También existe la posibilidad de anotar con Blastn para confirmar que los transcritos reconstruidos con el ensamblaje fueron secuenciados también en otros experimentos de EST. I.3.2. Anotaciones bioinformáticamente útiles Gracias a la anotación por similitud se puede enriquecer la información de los transcritos reconstruidos con una definición y con otros elementos de anotación, como la ontología de genes (GO, del inglés Gene Ontology ) [75], las rutas metabólicas en las que intervienen los transcritos, recogidas en los mapas KEGG [76], los dominios proteicos registrados en InterPro [77], y en caso de tener actividad enzimática, el código de la EC (del 39 INTRODUCCIÓN inglés, Enzyme Commission ), además de otras anotaciones no funcionales como SNP, SSR y miRNA. Gene Ontology: El objetivo del consorcio de la GO es producir un vocabulario controlado y dinámico, aplicable a todos los seres vivos, incluso aunque el conocimiento acerca de la función de los genes y de las proteínas en las células esté en continuo cambio [75]. Con este fin, se desarrollaron tres ontologías independientes: procesos biológicos, funciones moleculares y componentes celulares. La GO, como todas las ontologías, mantiene un lenguaje controlado que es útil para las personas y para los ordenadores, ya que cada elemento de la ontología tiene un código numérico (por ejemplo, GO:1234567) práctico para los análisis bioinformáticos, y una descripción informativa para los usuarios. Además, la GO mantiene una relación jerárquica entre sus elementos. La jerarquización de elementos permite la anotación a distintos niveles, debido al colapsamiento de anotaciones diferentes de la misma rama. Estas ontologías están disponibles en http://www.geneontology.org/ Mapas KEGG: Los mapas de rutas metabólicas de KEGG son diagramas gráficos que representan interacciones moleculares y redes metabólicas, procesos con información genética, ambiental, procesos celulares, sistemas orgánicos y enfermedades humanas [76]. Se pueden consultar en http://www.genome.jp/kegg/ , y son de gran utilidad para relacionar entre sí los genes que participan en una misma ruta, algo muy útil cuando se realizan análisis funcionales. Código EC: La Comisión Internacional para Enzimas se creó en 1956 en la Unión Internacional de Bioquímica y Biología Molecular para evitar que la misma enzima recibiera diferentes descripciones o nombres. Su función fue aportar un nombre descriptivo sobre la reacción que cataliza la enzima, y un código único para cada función enzimática. Las enzimas se nombran con cuatro números separados por puntos, EC 1.2.3.4, donde el primer número da la característica más genérica y los siguientes son cada vez más específicos. Así, el primer número puede tomar 6 valores; 1 para oxidoreductasas, 2 para las transferasas, 3 para las hidrolasas, 4 para las liasas, 5 para las isomerasas y 6 para las ligasas. Todo lo referente con estos códigos se puede encontrar en http://www.chem.qmul.ac.uk/iubmb/enzyme/ . Eso sí, solo las proteínas con actividad enzimática tendrán un código EC. InterPro: La base de datos InterPro integra modelos predictivos de varios repositorios (Pfam, PRINTS, PROSITE, SMART, ProDom, PIRSF, SUPERFAMILY, PANTHER, CATHGene3D, TIGRFAMs y HAMAP). Cada uno se centra en diferentes aspectos biológicos o utiliza una metodología distinta para encontrar el denominador común de las secuencias. El propósito de InterPro es combinar los puntos fuertes de cada uno de los repositorios para poner a disposición de la comunidad científica una única fuente con información estructurada sobre familias de proteínas, dominios y regiones funcionales [77]. InterPro está disponible en http://www.ebi.ac.uk/interpro/ . Existe una herramienta que se dedica expresamente a encontrar los dominios InterPro para una colección de proteínas concretas: InterProScan ( http://www.ebi.ac.uk/Tools/pfa/iprscan/ ). 46 INTRODUCCIÓN Sistemas operativos: UNIX OSX 10.6: Sus comandos se utilizaron para la ejecución de programas y para la visualización y manipulación de datos en general. UNIX SLES 10.2.5: Suse Linux Entrerprise Server se utilizó para la ejecución de programas en el sistema de colas de los supercomputadores Picasso SuperDome y Picasso - cluster, y para la consulta y manipulación de datos en general. I.5. Interés económico de las especies marinas Los recursos marinos son vitales dada la creciente demanda mundial de proteína animal, derivada del aumento demográfico (8 mil millones de seres humanos para 2030) y de los hábitos dietéticos actuales, que relacionan la salud y el consumo de pescado [95, 96]. El pescado es considerado un alimento fundamental en la dieta, ya que aporta proteínas de fácil digestión, ácidos grasos esenciales de la serie n-3 e importantes minerales y vitaminas, sobre todo A y D [97]. Su consumo ayuda a un adecuado equilibrio de ácidos grasos poliinsaturados n-3/n-6, relacionado con la expresión génica, lo que previene de modo natural algunas de las actuales enfermedades crónicas o degenerativas [98]. La producción pesquera mundial ha aumentado de forma constante en las últimas cinco décadas (figura I.14). El suministro de peces comestibles se ha incrementado a una tasa media anual del 3,2 %, superando así la tasa de crecimiento de la población mundial del 1,6 %. El consumo aparente mundial de pescado per cápita aumentó de un promedio de 9,9 kg en la década de 1960 a 19,2 kg en 2012, según las estimaciones preliminares [99]. Este incremento notable se ha debido a una combinación de crecimiento demográfico, aumento de los ingresos y urbanización, y se ha visto propiciado por la fuerte expansión de la producción pesquera y la mayor eficacia de los canales de distribución. Una gran parte de los recursos marinos son obtenidos a partir de productos de la pesca extractiva (figura I.14). Mientras la producción mundial de este tipo de pesca conoce una estabilidad en las últimas décadas, la producción de la acuicultura sigue creciendo y está alcanzando mundialmente un desarrollo espectacular (figura I.14), lo que lo convierte en el sector alimentario de más rápido crecimiento. Esta tendencia tiene que mantenerse o subir, para sustentar los actuales niveles de consumo de productos pesqueros, pues la población mundial sigue su aumento geométrico. 47 INTRODUCCIÓN Figura I.14: (tomada de [99]) Producción mundial de la pesca de captura y la acuicultura Aunque sea de relativamente reciente aparición, la acuicultura se ha caracterizado por demandar desde el principio una diversificación de especies que se puedan producir industrialmente. Fruto de esta precocidad, se emprendieron de forma generalizada numerosas actuaciones dirigidas a conocer las posibilidades de cultivo de un notablemente elevado número de peces. Durante los últimos tres lustros, son más de treinta las especies de peces que han sido objetivo de estudio en el área mediterránea [100]. Sin embargo, este esfuerzo se ha visto poco reflejado en la producción acuícola, que continúa monopolizada por las especies inicialmente desarrolladas, es decir: dorada ( Sparus aurata ) y lubina ( Dicentrarchus labrax ). Las razones para este desfase se podrían resumir en la dificultad tecnológica y la inadecuada posición de mercado de las numerosas especies exploradas. En España, estas dos especies junto con los mejillones forman las especies de acuicultura pioneras y ampliamente comercializadas. En la figura I.15 se puede ver el valor de la producción de éstas en comparación con otras menos comercializadas. Sin embargo, la saturación alcanzada en el mercado por su producción ha traído consigo la diversificación de especies acuícolas como el cultivo del besugo ( Pagellus bogaraveo ), el pargo ( Pagrus pagrus ), el abadejo ( Pollachius pollachius ), la corvina ( Argyrosomus regius ), el atún ( Thunnus thynus ) y el lenguado senegalés ( Solea senegalensis , Kaup 1858). Éste, junto con el lenguado común ( Solea solea ), desde los últimos treinta años aparecen en el sur de Europa como unos buenos candidatos para la diversificación de los mercados europeos. Una vez superados los problemas iniciales, el lenguado ha demostrado ser una prometedora especie para la acuicultura marina [101]. Existen estudios relacionados con su acuicultura en el litoral gaditano y en estuarios portugueses de hace ya más de veinte años [102]. Uno de los principales problemas con los que se encontró el cultivo del lenguado en sus inicios fueron las bajas tasas de crecimiento en juveniles y la aparente gran susceptibilidad a enfermedades, sobre todo pasteurelosis, vibriosis, mixobacteriosis y enfermedades virales. Estos problemas han sido en parte solventados a lo largo del tiempo mediante la mejora de las dietas y el estudio del estrés debido a las grandes densidades de cultivo. Sin embargo 48 INTRODUCCIÓN existen otros problemas relacionados con la reproducción en cautividad, el cultivo larvario, la nutrición y la respuesta immunitaria. Las cuales podrían ser solventados más fácilmente mediante un conocimiento en profundidad de los mecanismos moleculares y fisiológicos implicados. Figura I.15: (tomada de [103]) Valor de la producción (primera venta) de las especies de acuicultura marina (por grupos) en España en 2012, en millones de euros (JACUMAR) 49 OBJETIVOS Parte II Objetivos 50 OBJETIVOS OBJETIVOS 51 II. Objetivos El objetivo general de esta tesis es de mejorar el conocimiento a nivel genómico del lenguado senegalés ( Solea senegalensis ), el lenguado común ( Solea solea ), y de algunas especies afines a través de un análisis bioinformático y del desarrollo de protocolos de almacenamiento de los datos genómicos de estas especies, por el fin de suministrar elementos útiles que puedan abrir vías en la mejora genética de las mismas. Para ello se seleccionarán las mejores herramientas bioinformáticas disponibles y se optimizará su uso para sacar un rendimiento óptimo de los datos brutos disponibles. Así mismo, se desarrollarán aquellas nuevas herramientas y flujos de trabajo que sean necesarias para facilitar el tratamiento de los datos y mejorar la calidad de los resultados. Tanto los datos originales como los procesados se integrarán finalmente en una base de datos para ofrecer a la comunidad científica un transcriptoma y genoma depurado con las anotaciones más útiles. Los objetivos generales presentados anteriormente pueden desglosarse en los siguientes subobjetivos más específicos: 1Selección y comparación de los programas de ensamblaje de las secuencias disponibles de una parte los programas específicos de lecturas largas, de otra parte aquellos específicos de lecturas cortas. Además, establecer un modelo para el testeo, selección y optimización de los programas de ensamblaje que se van a utilizar y métodos para su uso combinado. 2Desarrollo de flujos de trabajo para el ensamblaje del transcriptoma y del genoma de la especies problema de este trabajo. Estos flujos incluyen tanto las herramientas disponibles en la bibliografía como las herramientas propias con el fin de utilizar los dos tipos de datos (lecturas largas o cortas) o una combinación de ambos, y proporcionar el resultado más optimizado para la posterior anotación. 3Inclusión la información en una base de datos web previamente desarrollada que permita almacenar los datos de modo estructurado y extensible para facilitar el acceso a ellos por parte de la comunidad científica. Aportar unas mejoras a la base de datos tanto en los tipos de información disponible como en la interfaz, jerarquía y organización de los datos así que la forma de acceder a ellos. 52 OBJETIVOS Figura II.1: Esquema de los objetivos. Creación de flujo completo de ensamblaje de las lecturas de transcriptomica. Elaboración de herramientas para el testeo y selección previa de los ensambladores usados y otras para la post evaluación y la mejora de los resultados. Posteriormente la información se integra en bases de datos disponibles para la comunidad científica 53 Parte III Materiales y métodos 54 MATERIALES Y MÉTODOS MATERIALES Y MÉTODOS 55 III. Materiales y métodos III.1. Equipos informáticos Durante el desarrollo de esta tesis se han utilizado varios ordenadores personales y supercomputadores, todos ellos con sistema UNIX. Los supercomputadores forman parte de la infraestructura del Centro de Supercomputación y Bioinformática que la Universidad de Málaga tiene en el edificio de Bioinnovación del Parque Tecnológico de Andalucía. Hasta 2012, el trabajo de esta tesis se realizó en dos supercumputadores: Picasso SuperDome y Picasso-Cluster. Picasso-cluster: Se trata de un cluster de 80 núcleos Intel Xeon E5450 a 3,00 GHz (con arquitectura x86) y con 160 GB de RAM, divididos en 10 blades con 8 núcleos y 2 GB de RAM para cada uno. Los blades están interconectados por red InfiniBand. El gestor de colas es PBS Pro y el sistema operativo es el Suse Linux Entrerprise Server v. 10.2.5. Picasso: se trata de un HP SuperDome con 128 núcleos Intel Itanium2 a 1,6 GHz y con 400 GB de memoria compartida. Ordenador y memoria se encuentran en dos racks conectados entre sí. Un tercer rack contiene el sistema operativo y el scratch (espacio del disco duro para el almacenamiento de datos temporales). El gestor de colas es PBS Pro y el sistema operativo es el Suse Linux Entrerprise Server v. 10.2.5. A principios del año 2013, se actualizaron todos los recursos de supercomputación, y todas las maquinas comparten ahora la misma arquitectura y están unificados con un único sistema de cola que envía los trabajos a los nodos más adecuados. Las nuevas instalaciones se componen en varios tipos de máquinas:  Clúster Intel E5-2670 con 768 cores a 2.60 GHz por core y con un RAM total de 3 TB y dispone de una red infiniband.  Máquinas de memoria compartida con 2 TB de RAM cada una. Tiene 560 cores a 2.40 GHz por core. Dispone de una red infiniband.  Clúster AMD Opteron 6176 con 984 cores, 4 TB de RAM y 246 TB de scratch.  Clúster de visualización ESX para máquinas virtuales : o 2 máquinas a 2.66 GHz y 32 GB de RAM cada una. o 3 máquinas a 2.66 GHz y 16 GB de RAM cada una. o 4 máquinas a 2.13 GHz y 124 GB de RAM cada una. o Almacenamiento compartido con un total de 750 TB (brutos). 62 MATERIALES Y MÉTODOS Para más información consúltese el manual del programa, mencionado al comienzo de este apartado. MIRA3 fue uno de los programas que utilizamos para ensamblar las lecturas largas de los transcriptomas tratados en este trabajo. III.3.2.3. EULER-SR EULER [110] es un algoritmo para ensamblaje de secuencias cortas y de los primeros que utiliza un camino euleriano implementado mediante un grafo de Bruijn. Según sus autores tiene la ventaja con respecto a los mecanismos de diseño por solapamiento de que es mucho más eficaz con las repeticiones ya que en vez de enmascararlas genera un grafo que permite crear una estructura de repeticiones del genoma para poder tratarlas. El software EULER fue desarrollado para lecturas Sanger y fue posteriormente modificado y denominado EULER-SR [111, 112] para lecturas cortas de 454 GS20 [9], lecturas cortas no pareadas de Illumina/Solexa [113]. EULER se basa en tres principios fundamentales:  Representa las lecturas como enlaces y los solapamientos como nodos en un grafo de Bruijn.  Realiza un ensamblaje eficiente mediante la reducción a un problema de camino euleriano: cada enlace debe visitarse una sola vez.  Las repeticiones se tratan mediante el uso de múltiples enlaces para una lectura repetida. Es fácil de ejecutar, basta con indicar el fichero de entrada y el valor de k-mero que se quiere utilizar para obtener el resultado en el fichero my_results_file.txt: Assemble.pl my_file.fasta kmer_value > my_results_file.txt Para las secuencias de transcriptómica 454, Euler fue el ensamblador de tipo De Bruijn que mejores resultados devolvía, por lo que lo utilizamos como parte del flujo de trabajo del ensamblaje de este tipo de lecturas. III.3.2.4. Velvet Velvet [11, 114] es un conjunto de algoritmos desarrollados por Daniel R. Zerbino y Ewan Birney para el tratamiento de grafos de de Bruijn para ensamblaje de secuencias genómicas. Velvet hace un uso bastante extenso de herramientas de simplificación de grafos para reducir caminos no intersectantes con nodos simples. Esta simplificación comprime el grafo sin pérdida de información. El programa ejecuta la fase de simplificación durante la 63 MATERIALES Y MÉTODOS construcción del grafo y, de nuevo, en varias ocasiones durante el proceso de ensamblaje. Esta técnica, introducida como «eliminación de únicos» para grafos de k-mero es similar a la formación de unitigs en los grafos de solapamiento [61] y en los ensambladores OLC [61]. La ejecución de Velvet se hace en dos etapas de la forma siguiente: srun velveth out_folder my_kmer -fastq -longPaired –separate long_paired_reads1.fastq long_paired_reads2.fastq -long long_single_reads.fastq -shortPaired –separate short_paired_reads1.fastq short_paired_reads2.fastq > salida1 srun velvetg Velveth_29 -exp_cov auto -cov_cutoff 10 > salida2 En la primera línea de comando, my_kmer se reemplaza con al valor del k-mero, en la última versión de Velvet este valor puede ir hasta 101. –fastq es para indicar que los ficheros introducidos son de tipo .fastq. –longPaired indica que las lecturas son largas y pareadas, -separate indica que las parejas de lecturas vienen en dos ficheros independientes, -long indica que las lecturas son largas y single, -shortPaired indica que las lecturas son cortas y pareadas y por ultimo out_folder es el nombre de carpeta donde se guardarán los resultados de la ejecución. En la segunda línea de comando, out_folder representa la carpeta creada anteriormente, -exp_cov indica la cobertura esperada y -cov_cutoff indica la cobertura mínima. Velvet formó parte del estudio que hicimos sobre la comparación de los programas de ensamblaje para las lecturas largas transcriptomicas y genómicas. III.3.2.5. Oases Oases [115] es un extensión de Velvet desarrollada por Marcel Shulz y Daniel Zerbino, destinada a ensamblar lecturas de transcriptómica . Oases importa un ensamblaje preliminar producido por Velvet [11, 114], y agrupa los contigs en grupos pequeños llamados loci. Luego aprovecha de las lecturas pareadas y las lecturas largas, para construir isoformas utilizando la similitud entre los grafos de Bruijn y los grafos de ayuste (splicing). La ejecución de Oases se compone de tres líneas de comando. Las dos primeras corresponden a una ejecución normal del Velvet tal como se describió más arriba con la diferencia de que en el segundo comando se añade la opción -read_trkg yes que permite el seguimiento de las posiciones de lecturas en el ensamblaje. En tercera línea es donde interviene el Oases recuperando los datos preliminares de Velvet contenidos en out_fol der. srun velveth out_folder my_kmer -fastq -longPaired –separate long_paired_reads1.fastq long_paired_reads2.fastq -long long_single_reads.fastq -shortPaired –separate short_paired_reads1.fastq short_paired_reads2.fastq > salida1 srun velvetg out_folder -read_trkg yes > salida2 srun oases out_folder -ins_length 169 -ins_length_sd 47 -min_trans_lgth 100 -cov_cutoff 3 -min_pair_count 20 > salida3 64 MATERIALES Y MÉTODOS En la tercera línea donde interviene el programa Oases, la opción -min_pair_count indica el número de pares de lecturas necesarios para considerar una conexión entre dos contigs montar un scaffold. Oases se utilizó en el ensamblaje de las lecturas de Solea senegalensis y Solea solea con los parámetros que dieron el mejor resultado en la fase de pruebas. III.3.2.6. SOAPdenovo-trans Es un ensamblador de novo de transcriptomas basado en el SOAPdenovo2 [116] y utiliza un algoritmo de ensamblaje paralelizable basado en grafos De Bruijn. Está preparado para manejar el ayuste alternativo y a los diferentes niveles de expresión en los transcritos. El paquete el programa consiste en dos ejecutables cuyos la utilización de uno u otro depende del k-mero usado.  SOAPdenovo-Trans-31kmer: para k-meros hasta 31  SOAPdenovo-Trans-127mer: para k-meros hasta 127 La ejecución del programa se efectúa en cuatro etapas de la forma siguiente: SOAPdenovo-Trans-31kmer pregraph -s config_file -K 29 -p 16 -d 1 -o my_assembly SOAPdenovo-Trans-31kmer contig -g my_assembly SOAPdenovo-Trans-31kmer map -s config_file -g my_assembly -K 29 -p 16 SOAPdenovo-Trans-31kmer scaff -g my_assembly -F -p 16  En la primera línea donde se utiliza el comando pregraph , –s indica el fichero de configuración donde viene la información relativa a las lecturas que se van a utilizar (ver descripción más abajo), -K indica el valor de k-mero, -p indica el número de CPU, -d indica la frecuencia mínima de k-mero a considerar, -o indica el nombre de carpeta del proyecto.  En la segunda línea donde se utiliza el comando contig se indica el nombre de la misma carpeta del proyecto -g .  En la tercera línea donde se utiliza el comando map con –s se indica el fichero de configuración, con –g la carpeta del proyecto con –K el valor de k-mero, con –p el número de CPU.  En la cuarta línea donde se utiliza el comando scaff, -g indica la carpeta del proyecto, -F permite rellenar los huecos en los scaffolds y –p el número de CPU. También se puede ejecutar todas las etapas de forma seguida utilizando el comando all SOAPdenovo-Trans-31kmer all -s config_file -K 29 -p 16 -d 1 -D 4 -F -o output 65 MATERIALES Y MÉTODOS El fichero de configuración tiene el siguiente formato: #longitude máxima de las lecturas max_rd_len=150 [LIB] #tamaño medio de inserto avg_ins=180 #Si la secuencia se tiene que revertir, poner 0 para indicar que no reverse_seq=0 #En que parte(s) utilizar la lectura (escoger 3 para utilizarla tanto en contigs como en scaffolds) asm_flags=3 #En qué orden las lecturas se utilizarán en el scaffolding (se considera en caso que haya más de una librería de lecturas). rank=1 # Número mínimo de pares de lecturas pareadas para realizar una conexión. pair_num_cutoff=10 #Lecturas pareadas en formato fastq (utilizar q1 y q2) q1=paired_reads1.fastq q2=paired_reads2.fastq #lecturas simples en formato fastq (utilizar q) q=sequences_illum_454.fastq Se pueden añadir más librerías escribiendo la etiqueta [LIB] seguida con los nombres de lecturas y sus características siendo el mismo modelo descrito arriba. SOAPdenovo-trans fue utilizado para ensamblar los datos cortos de Illumina de Tisochrysis lutea y Ruditapes decussatus III.3.2.7. CABOG CABOG (Celera Assembler with Best Overlap Graph) [10, 54] es un ensamblador de novo de tipo OLC de secuencias largas (por ejemplo Sanger y 454) de ADN para genomas completos. CABOG ha contribuido de manera importante al avance de la genómica, incluyendo el primer ensamblaje completo un genoma de un organismo multicelular y el primer genoma diploide de un individuo humano. El software de CABOG es un sistema modular compuesto de múltiples programas que interactúan a través de interfaces bien definidas de tal manera que se puede modificar el orden secuencias de los programas y sus parámetros para cubrir necesidades específicas. La ejecución de CABOG se realiza con las dos líneas de comando siguientes con este mismo orden: fastqToCA -insertsize 5000 1500 -libraryname my_library -technology '454' -type 'sanger' -innie -reads my_single_reads.fastq -mates my_paired_reads.fastq > out_file.frg runCA -d my_work_folder -p my_prefix -unitigger=utg out_file.frg 66 MATERIALES Y MÉTODOS En la primera línea de comando –insertsize indica la información relativa a los tamaños de inserto a saber la media y la desviación típica (los dos valores separados por espacio), -libraryname indica el nombre de la librería de lecturas, -technology indica la tecnología de secuenciación de las lecturas, -innie indica que la orientación de las lecturas pareadas es “Forward-Reverse” (en caso que la orientación es “Reverse-Forward”, se especifica –outtie ), -reads indica el nombre del fichero .fastq de lecturas single, -mates indica el nombre del fichero de lecturas pareadas, -out_file.frg indica el nombre de fichero de salida .frg. Las lecturas pareadas vienen intercaladas en el mismo fichero .fastq. En la segunda línea de comando –d indica el nombre de la carpeta del trabajo, -p indica el prefijo de los ficheros generados, –unitigger indica el modulo utilizado para construir unitigs. CABOG ofrece tres módulos, de los cuales se utiliza solo uno a la vez. El modulo original (utg) es el más adecuado para datos de Sanger, el módulo “best overlap graph unitigger” (bog) es el más compatible con los datos de 454 y el módulo “bogart unitigger” (bogart) es el más compatible para los datos de Illumina. Al final de la línea de comando se añade el nombre del fichero .frg generado anteriormente. CABOG fue utilizado para ensamblar las lecturas 454 de Photobacterium damselae . III.3.2.8. Ray Ray [68], es un ensamblador que utiliza el grafo de De Bruijn para construir su estructura de datos. La principal característica que podemos destacar de este ensamblador es que es un ensamblador que utiliza memoria distribuida para realizar sus operaciones. Los otros ensambladores disponibles, utilizan memoria compartida. Por lo tanto el ensamblador Ray utiliza MPI [117] como interfaz de paso de mensajes entre todos los nodos que estén funcionando de manera simultánea. Gracias a esta característica podemos aprovechar aún más el potencial del superordenador Picasso, y no tendremos que limitarnos a las 16 CPU que contienen los nodos de cálculo o las 64 CPU que contienen los nodos bigmem. La ejecución del programa Ray se hace de la forma siguiente: mpiexec -np 256 Ray -k 51 -p my_paired_reads1.fastq my_paired_reads2.fastq -s my_single_reads.fastq -route-messages -connection-type debruijn -routing-graph-degree 4 -o illumina_01 donde –np indica el número de CPU utilizadas, –k indica el valor de k-mero, -p indica los nombres de los dos ficheros de lecturas pareadas separados por espacio, -s indica el nombre del fichero de lecturas simples. El programa Ray tiene una amplia gama de configuraciones que se pueden consultar en http://denovoassembler.sourceforge.net/manual.html. En caso que se trate de ensamblajes donde se utilizan grandes cantidades de lecturas, convendría activar los routing messages , añadiendo a la línea de comando lo siguiente: -route-messages -connection-type debruijn -routing-graph-degree 4 67 MATERIALES Y MÉTODOS El grado de los nodos de grafos utilizando la opción -routing-graph-degree según el número de núcleos que se van utilizar tal como viene explicado en la tabla III.1 Número de CPUs Grado de nodos Diámetro de nodos Configuración 256 4 4 (4*4*4*4=256) -routing-graph-degree 4 512 8 3 (8*8*8) -routing-graph-degree 8 1024 4 5 (4*4*4*4*4) -routing-graph-degree 4 1024 32 2 (32*32) -routing-graph-degree 32 1024 2 10 (2^10) -routing-graph-degree 2 Tabla III.1: Configuración de los Grados de nodos De Bruijn en el programa Ray, según el número de CPUs utilizado Ray fue utilizado para ensamblar las secuencias de Illumina de Solea seneganensis III.3.3. Para el tratamiento y mejora de los ensamblajes III.3.3.1. GAM-NGS En el ensamblaje del genoma de Solea senelgalensis (Apartado IV.3.2.2) era necesario reconciliar los contigs procedentes de las diferentes librerías de lecturas Illumina. Para ello se utilizó GAM-NGS [118] que es un programa que se utiliza para combinar dos o más ensamblajes genómicos con el fin de mejorar la contigüidad y la exactitud, y así obtener un ensamblaje único (de reconciliación) con mejores características que los originales. Las regiones que representan el mismo locus (llamados bloques) en los ensamblajes se identifican a través del mapeo de las lecturas y se guardan en un grafo ponderado. La fase de la combinación se lleva a cabo con la ayuda este grafo, lo que permite también la resolución de las regiones problemáticas. Para reconciliar dos ensamblajes con GAM-NGS es necesario elegir uno de ellos como «maestro» (ensamblaje de base). El otro se definirá como «esclavo» y que se utilizará para complementar el primero. Para reconciliar ensamblajes con GAM-NGS se requiere una etapa previa de mapeo de las lecturas sobre cada uno de ellos utilizando Bowtie2 [119], posteriormente se generará un fichero de alineamiento indexado en formato binario, BAM, utilizando samtools [120]. Estos dos últimos programas se describirán de forma más detallada en el apartado III.3.4. La etapa de mapeo de lecturas se realiza mediante las siguientes líneas de comando: module load bowtie/2.2.4 bowtie2-build -f assembly_1.fasta my_indexes 68 MATERIALES Y MÉTODOS bowtie2 my_indexes -q -1 paired_reads1.fastq -2 paired_reads2.fastq -U single_reads.fastq --very-sensitive -p 8 -S assembly_1.sam module load samtools #convertir el fichero de format SAM a formato BAM samtools view -bS assembly_1.sam > assembly_1.bam #ordenar el fichero BAM samtools sort assembly_1.bam sorted.assembly_1 #filtrar el fichero BAM (eliminar las lecturas no mapeadas) samtools view -b -F 4 sorted.assembly_1.bam > filtered.sorted.assembly_1.bam # crear un index al fichero BAM samtools index filtered.sorted.assembly_1.bam La ejecución del programa GAM-NGS se realiza en varias etapas: PATH=`pwd -P` THREADS_NUM=4 GAM_CREATE=gam-create GAM_MERGE=gam-merge # Etapa 1: Creacion de una carpeta donde guardar los ficheros de salida mkdir -p ${PATH}/gam-ngs_merge # Etapa 2: Preparacion de los ficheros de entrada echo -e "${PATH}/Alignments/assembly_1/sorted.assembly_1.bam\n300 500" >${PATH}/gam-ngs_merge/assembly_1.pe.list.txt echo -e "${PATH}/Alignments/assembly_2/sorted.assembly_2.bam\n300 500” >${PATH}/gam-ngs_merge/assembly_2.pe.list.txt # Etapa 3: Construccion de los bloques ${GAM_CREATE} --master-bam ${PATH}/gam-ngs_merge/assembly_1.pe.list.txt --slave-bam ${PATH}/gam-ngs_merge/assembly_2.pe.list.txt --min-block-size 10 --output ${PATH}/gam-ngs_merge/out >${PATH}/gam-ngs_merge/gam-create.log.out 2>${PATH}/gamngs_merge/gam-create.log.err # Etapa 4: Reconciliacion ${GAM_MERGE} --blocks-file ${PATH}/gam-ngs_merge/out.blocks --master-bam ${PATH}/gamngs_merge/assembly_1.pe.list.txt --master-fasta ${PATH}/Assembly/assembly_1/assembly_1.fasta --slave-bam ${PATH}/gam-ngs_merge/assembly_2.pe.list.txt --slave-fasta ${PATH}/Assembly/assembly_2/assembly_2.fasta --min-block-size 10 --output ${PATH}/gam-ngs_merge/out --threads ${THREADS_NUM} >${PATH}/gam-ngs_merge/gam-merge.log.out 2>${PATH}/gam-ngs_merge/gam-merge.log.err En la etapa 2 se genera en un archivo de configuración donde se imprimen los nombres de ficheros del mapeo con sus rutas así que el tamaño de inserto mínimo y máximo. 69 MATERIALES Y MÉTODOS En la etapa 3, --master-bam indica el nombre del fichero BAM de mapeo de las lecturas sobre el ensamblaje considerado como maestro, --slave-bam indica el nombre del fichero BAM de mapeo de las lecturas sobre el ensamblaje considerado como esclavo. En la etapa 4, --master-fasta indica el fichero de contigs o scaffolds del ensamblaje maestro, --slave-fasta indica indica el fichero de contigs o scaffolds del ensamblaje esclavo, --min-block-size indica el número de lecturas mínimo para considerar un bloque , --output indica el prefijo de los ficheros de salida, y con --threads se indica el número de cores a utilizar. III.3.3.2. SOAPdenovo Scaffolder En el ensamblaje del genoma de Solea senegalensis (Apartado IV.3.2.2) teníamos disponibles muchas librerías de lecturas que fueron aprovechadas para unificar los contigs. Uno de los programas que utilizamos para este propósito fue SOAPdenovo Scaffolder que realmente es un módulo del programa de ensamblaje genómico SOAPdenovo [116] que se utiliza para unir los contigs preensamblados con SOAPdenovo, sin embargo también puede usarse para unir contigs generados por otros ensambladores. La ejecución del programa requiere una preparación previa de los contigs utilizando la herramienta finalFusion ( https://sourceforge.net/projects/soapdenovo2/files/Prepare/ ), lo cual se hace con la siguiente línea de comando: finalFusion -g Scaff -K 31 -c contigs.fasta -D donde –g indica el prefijo de para la salida, -K indica el tamaño de k-mero, -c indica el fichero de contigs y con -D se activa el modo de preparación de los contigs. Se generaran varios ficheros que tienen el prefijo Scaff. La unificación de los configs ( scaffolding ) se realiza con las siguientes líneas de comando: # mapeo de las lecturas sobre los contigs SOAPdenovo-63mer map -p 16 -s config -g Scaff 1>map.log 2>map.err # unificacion de los contigs SOAPdenovo-63mer scaff -p 16 -g Scaff -N 600000000 -F 1>scaff.log 2>scaff.err donde –p indica el número de cores a utilizar, -s indica el fichero de configuración (véase el detalle de este fichero en el apartado III.3.2.6) –g indica el prefijo de los ficheros generados en la etapa de preparación, -N indica el tamaño esperado del genoma y con –F se activa el relleno de los huecos en los scaffolds 70 MATERIALES Y MÉTODOS III.3.3.3. SSPACE En los flujos de trabajo de los ensamblajes genómicos de Photobacterium damselae (Apartado IV.3.1.2) y Solea senegalensis (Apartado IV.3.2.2) era necesario unificar los scaffolds aprovechando de la información de las lecturas disponibles. Para ello utilizamos SSPACE [121] que es un programa independiente que sirve para el scaffolding de los contigs preensamblados mediante el uso de lecturas pareadas y tiene la ventaja de ser utilizado también para unir scaffolds . Este programa utiliza la información del mapeo de las lecturas pareadas sobre los contigs o scaffolds para evaluar el orden, la distancia entre ellos y la orientación por el fin de establecer conexiones entre los mimos. La ejecución de SSPACE se realiza con una sola línea de comando de la siguiente forma: SSPACE_Standard_v3.0.pl -l libraries.txt -s contigs.fasta -T 16 donde –l indica el nombre de un fichero de texto incluyendo la información sobre las lecturas pareadas que se utilizarán en la unificación de los contigs o scaffolds , –s indica el nombre del fichero de contigs o scaffolds de entrada, -T indica el número de CPU a utilizar. El fichero libreries.txt tiene el siguiente formato: lib1 bowtie lib1_paired_reads1.fastq lib1_paired_reads2.fastq 3000 0.25 FR lib2 bowtie lib2_paired_reads1.fastq lib2_paired_reads2.fastq 8000 0.25 FR donde la columna 1 corresponde el nombre de la librería, la columna 2 corresponde al programa utilizado para mapear las lecturas sobre los contigs o scaffolds en el cual se puede elegir Bowtie [119] o BWA [122] las columnas 3 y 4 corresponden a los nombres de ficheros de lecturas pareadas, las columnas 5 y 6 corresponden al tamaño de inserto de las lecturas pareadas y el error mínimo permitido en el tamaño de inserto respectivamente. Por ejemplo para un tamaño de inserto 3000 y un error de 0,25, la distancia entre las lecturas puede tener un error de 3000 * 0,25 = 750 en ambos sentidos. La última columna se refiere a la orientación de las lecturas pareadas (FR: directa-inversa, RF: inversa-directa, FF: directa - directa, RR: inversa-inversa) III.3.3.4. SOAPdenovo GapCloser SOAPdenovo GapCloser ( http://sourceforge.net/projects/soapdenovo2/files/ GapCloser ) es un programa incluido en el paquete de SOAPdenovo que se utiliza para cerrar los huecos (indicado con N yuxtapuestas) que se crean durante la etapa de scaffolding con SOAPdenovo Scaffolder o con otro ensamblador . Está destinado a genomas grandes de plantas y animales aunque también funciona bien para genomas de bacterias y hongos. La ejecución de este programa se realiza con la siguiente orden: GapCloser -b config_file -a scaffolds.fasta -l 155 -t 64 -o gap_closed_scaffolds.fasta 71 MATERIALES Y MÉTODOS donde –b indica el nombre del fichero de configuración específico a SOAPdenovo (verse el detalle de este fichero apartado III.3.2.6), -a indica el nombre del fichero de scaffolds de entrada –l indica la longitud máxima de las lecturas (este valor puede tomar un valor máximo de 155), -t indica el número de cores a utilizar, y -o el nombre del fichero de salida. SOAPdenovo GapCloser se utilizó en el flujo de trabajo de los ensamblajes genómicos de Photobacterium damselae y Solea senegalensis , para rellenar los huecos con nucleótidos para mejorar la calidad de los scaffolds. III.3.4. Para analizar secuencias III.3.4.1. SeqTrimNext Antes de proceder al ensamblaje de las lecturas Roche/454 y Illumina tanto transcriptomicas como genómicas fue necesario preprocesarlas para descartar los fragmentos de baja calidad, contaminantes, vectores, adaptadores, y otros artefactos, para ello utilizamos SeqTrimNext. Se trata de una herramienta para el preprocesamiento de secuencias de nueva generación desarrollada en la Plataforma Andaluza de Bioinformática. Se puede utilizar por línea de comandos, como servicio web basado en REST y como herramienta web. Para más información véase su portal http://www.scbi.uma.es/ seqtrimnext. A continuación se muestra un ejemplo de la ejecución de SeqTrimNext a través del sistema de colas de Picasso. # se inicializa SeqTrimNext en el Cluster . ~seqtrimnext/init_env # se recogen las CPU asignadas. srun hostname -s > workers seqtrimnext -t my_template.txt -Q paired_seq1.fastq,paired_seq2.fastq -w workers en la que -t indica que el fichero my_template.txt es una plantilla en la que se detallan los parámetros que se desean utilizar. En la versión web del programa se ofrecen plantillas con valores por omisión para varias situaciones (genómica, transcriptómica, amplicones, secuencias de plantas, etc.). El parámetro -w indica que el programa se ejecutará en los núcleos con las ID que se recogen en el fichero workers . Con -Q se indican los nombres de los ficheros con las secuencias (conviene consultar la ayuda del programa para ver más opciones y formatos de entrada). 78 MATERIALES Y MÉTODOS III.3.5.11. Gigabayes Para detectar las variaciones genéticas SNP en los ensamblajes de las lecturas Roche/454 se utilizó Gigabayes [131] que es un programa optimizado para el análisis de millones de lectura NGS alineados a una secuencia consenso que puede ser un transcrito o un ADN genómico. El descubrimiento de los SNP con Gigabayes se hace a través del análisis de un fichero .ACE que es un fichero de alineamiento de las lecturas sobre la secuencia del consenso, el cual lo generan la mayoría de los ensambladores OLC. El análisis de Gigabayes se hace en dos etapas. En la primera etapa se crea el fichero binario a partir de los datos de alineamiento de las lecturas y en la segunda se analiza la información contenida en este fichero binario. gigaBuild --fd my_file.fasta --fq my_file.qual --ace ./my_file.ace --gig my_file.ace.gig gigaBayes --gig my_file.ace.gig --gff newfile.ace.gff --anchor --ploidy diploid -- indel --debug --CRL 4 --PSL 0.9 --log my_file.ace.log  En la primera línea de comando --fd indica el fichero de lecturas utilizado en el ensamblaje , --fq indica el fichero de calidades, --ace indica el fichero de alineamiento y --gig indica fichero binario de salida.  En la segunda línea de comando --gig indica fichero binario generado anteriormente, --gff indica un fichero de tipo .gff que se general durante el análisis, --anchor indica la utilización de secuencia de anclaje para buscar el polimorfismo. --ploidy se refiere a la ploidía del organismo (“haploid” o ”diploid”), la opción --debug permite que se muestren mensajes de depuración del programa. - -CRL indica la cobertura mínima en una posición para ser considerada. --PSL indica el valor mínimo de probabilidad de un SNP para que se muestre en el informe, -- log indica el fichero de salida donde se imprime toda la información relativa a los SNPs encontrados. Para incluir en el análisis la búsqueda de los Indels se añade la opción --indel . Gigabayes se utilizó para descubrir los SNP en los proyectos de ensamblaje de Solea senegalensis donde se utilizaron únicamente las lecturas de 454/Roche. 79 MATERIALES Y MÉTODOS III.3.4.12. Tablet Se trata de una herramienta gráfica interactiva para la visualización de ensamblajes [132]. Admite multitud de formatos diferentes: ACE, AFG, MAQ, SOAP2, SAM, BAM, FASTA, FASTQ y GFF3. Puede descargarse en http://bioinf.scri.ac.uk/tablet/ y es muy sencillo de utilizar, tan solo hay que abrir el programa y elegir el fichero de ensamblaje que se quiere visualizar. Tablet se utilizó para la visualización de los ensamblajes de lecturas Roche/454 y para observar las corrección aportadas por el programa CoMiner (Apartado IV.1.3) III.3.4.13. GEvo Para poder comprobar la colinealidad de los genomas de Solea senegalensis y Cynoglossus semilaevis (Apartado V.3.2.3.2), usamos GEvo ( https://genomevolution.org /CoGe/GEvo.pl ) [133] que es una herramienta web diseñada para comparar regiones genómicas entre múltiples organismos utilizando diferentes algoritmos de comparación, lo cual permite de identificar patrones en la evolución de los genomas. III.3.4.14. NUCMER En el ensamblaje del genoma de Solea senegalensis ( Apartado IV.3.2.2), era necesario alinear los contigs por el fin de combinar aquellos que solapan entre ellos. Para ello se utilizó NUCMER que es un programa del paquete MUMmer [134] muy utilizado el alineamiento global de genomas enteros que sean completos o en estado de borrador y que tiene la ventaja de ser muy rápido en sus tareas de alineamiento. Se ejecuta con los dos siguientes ordenes sucesivos: nucmer -p nucmer reference_genome.fasta query_genome.fasta show-coords -T nucmer.delta > nucmer.coords En la primera línea, reference_genome.fasta indica el nombre del fichero fasta que incluye al genoma escogido como referencia, query_genome.fasta indica el nombre del fichero fasta que incluye al genoma en el cual se buscarán alineamientos con el primer genoma. Con esta línea se genera el fichero nucmer.delta III.3.5. Scripts escritos En este apartado se describen los scripts que fueron escritos a lo largo de este trabajo para tratar y procesar los datos genómicos de las especies estudiadas. 80 MATERIALES Y MÉTODOS III.3.5.1. Para generar lecturas artificiales Este script escrito en perl (apéndice A), se utilizará en el apartado IV.1.2 para generar lecturas artificiales por el fin de testear los programas de ensamblaje. Se ejecuta con la orden: generate_artif_rd.pl sequences.fasta 300 500 30 donde el primer parámetro representa el nombre del fichero de secuencias de donde se extraerán las lecturas artificiales, el segundo parámetro representa la longitud mínima de las lecturas, el tercer parámetro la longitud máxima y el cuarto parámetro la cobertura de las lecturas. Este script se utilizó para crear secuencias artificiales para utilizarlas en la verificación de los programas de ensamblaje de lecturas largas. III.3.5.2. Para combinar dos ficheros ACE Es un script en Ruby que combina dos ficheros ACE de MIRA y CAP3 (ver apartado IV.1.2.3 más adelante). Utiliza la gema scbi_ace (apartado III.3.1.2) y se ejecuta de la siguiente manera: combine_ace_files.rb assembly_cap3.ace assembly_mira.ace mira_reads_index donde el primer parámetro representa el nombre del fichero ACE de CAP3 el segundo parámetro representa el nombre del fichero ACE de MIRA, y el tercer parámetro el índice de los identificadores de secuencias de MIRA Este script se utilizó para combinar los ficheros ACE de MIRA y CAP3 en los ensamblajes de algunas versiones del transcriptoma Solea senegalensis donde solo se utilizarón lecturas 454/Roche. III.3.5.3. Para calcular el tamaño de inserto real en lecturas pareadas. Este script escrito en Ruby (apéndice B) utiliza el fichero SAM de mapeo de una librería de lecturas sobre una referencia (que puede ser un ensamblaje preliminar de estas lecturas formado de contigs) para calcular el tamaño de inserto real de esta librería. El script se utiliza de la siguiente forma: insert_size_from_sam.rb mapping_file.sam en el que mapping_file.sam representa el nombre del fichero de mapeo previamente mencionado. Este script fue utilizado para calcular el tamaño de inserto real en las lecturas pareadas (Roche/454 o Illumina) utilizadas en este trabajo. 81 MATERIALES Y MÉTODOS III.3.5.4. Para eliminar ensamblajes erróneos Este script escrito en Ruby recorre el resultado de alineamiento Blast entre los transcritos procedentes de lecturas cortas con los procedentes de las lecturas largas (considerados como más fiables) y entre los transcritos procedentes de lecturas cortas detecta aquellos que tienen una estructura sentido/antisentido para poder eliminarlos. Se ejecuta de la siguiente forma: misassembled_trascripts.rb blast_alignment 95 80 donde el primer parámetro indica el fichero Blast de alineamiento, el segundo indica el porcentaje de identidad mínimo para considerar un alineamiento, y el tercero la longitud de alineamiento mínimo. Este script se utilizó en los flujos de trabajo de los ensamblajes de transcriptomica para eliminar los ensamblajes erróneos en los contigs procedentes de los ensamblajes illumina. III.3.5.5. Para seleccionar los transcritos que representan las sondas de microarray. Es un script en Ruby que utiliza el resultado de Full-LengtherNext y el de dos predictores de secuencias codificantes: testcode [135] y ORF predictor [136] ( http://proteomics.ysu.edu/tools/OrfPredictor.html ), para seleccionar automáticamente los transcritos que representarán las sondas para un experimento de microarrays. La ejecución se realiza de la siguiente forma: fln2microarray.rb fln_annotation_file test_code_file orf_predictor_file transcripts.fasta 44000 donde fln_annotation_file indica el nombre del fichero de anotación de Full-LengtherNext, test_code_file indica el nombre del fichero la salida de testcode, orf_predictor_file indica el nombre del fichero de salida del programa ORF predictor, transcripts.fasta indica el nombre del fichero fasta de transcritos y con el ultimo parámetro se indica el número de transcritos requeridos para el microarray. Este script se utilizó para seleccionar los transcritos de Solea senegalensis que se enviaron para realizar un experimento de microarrays III.3.5.6. Para calcular las estadísticas sobre los marcadores SSR Es un script escrito en Ruby (apéndice C) que se utiliza para recorrer el fichero de salida de MREPS (apartado III.3.4.10) que sirve para de detección de marcadores SSR en un transcriptoma, y desde ello calcula estadísticas sobre estos marcadores tanto los que están dentro del ORF como aquellos que están en la UTR. 82 MATERIALES Y MÉTODOS La ejecución de este script se realiza de la siguiente forma: ssrs_stats.rb mreps_file transcripts_ids fln_pt_seqs donde mreps_file indica el nombre del fichero de salida de MREPS, transcripts_ids indica el nombre de un fichero que incluye a los identificadores de transcritos (un id por línea) de los cuales se obtendrán las estadísticas de los SSR incluidos en los mismos. fln_pt_seqs indica el nombre del fichero de anotación de Full-LengtherNext Este script se utilizó para extraer las estadísticas sobre los marcadores SSR detectados en los transcriptomas estudiados. III.3.5.7. Para solapar los contigs genómicos En los ensamblajes genómicos, los contigs pueden compartir zonas comunes en sus bordes. Éstas se pueden detectar con un programa de alineamiento global como NUCMER (apartado III.3.4.14). Con este script en Ruby se utiliza el resultado de NUCMER para solapar los contigs según unos criterios bien definidos y utilizando la información de sintenia con ICMapper (ver apartado III.3.5.9 a continuación) de dos especies cercanas diferentes. Globalmente, hay dos criterios para solapar los contigs: sea que tengan unos alineamientos muy estrictos (longitud y porcentaje de identidad del alineamiento muy altos), o que estos parámetros sean mínimos y que los dos contigs se encuentran contiguos sobre la base de la sintenia con especies cercanas. Después de obtener el fichero de coordenadas del alineamiento del genoma consigo mismo utilizando NUCMER, la ejecución del script se realiza de la siguiente forma: contigs_ovlp.rb –n nucmer.coords –l 100 –i 80 –fr species_1.ordering –sr species_2.ordering –c contigs.fasta donde –n indica el nombre del fichero de coordenadas que representan el alineamiento de NUCMER. –l indica la longitud mínima inicial de alineamiento (filtro preliminar), -i indica el porcentaje de identidad inicial mínimo (filtro preliminar), -fr y –sr indican la información de ordenamiento sobre la base de la sintenia con dos especies cercanas. Con –c se indica el nombre del fichero de contigs en los cuales se realizaran los solapamientos. Este script fue utilizado para combinar los contigs solapantes del genoma de Solea senegalensis. III.3.5.8. Para la corrección automática de los ensamblajes (Cominer) Este programa esta descrito de forma detallada más adelante en el apartado IV.1.3. Su ejecución se realiza de esta forma: 83 MATERIALES Y MÉTODOS cominer.rb –a assembly.ace –t 40 –e 2 –l 40 -p donde –a indica el fichero ACE a analizar –t indica la longitud máxima a recortar en porcentaje de la lectura, -e indica el porcentaje mínimo de los errores en las lecturas para esta sea considerada. –l indica la longitud mínima de lecturas a considerar en el ensamblaje nuevo, -p es para generar un plot sobre la zonas de distribución de los errores. III.3.5.9. Para ordenamiento de los contigs o scaffolds y acabado de los genomas (ICMapper) Este programa escrito en Perl ordena los contigs o scaffolds de un borrador de genoma basándose en una referencia; esta descrito de forma detallada más adelante en el apartado IV.1.4; ICMapper en se ejecuta en dos etapas: Alineamiento Blast entre los contigs o scaffolds y un genoma de referencia acabado: blastn -task 'blastn' –db reference.fasta -query scaffolds.fasta –dust no -culling_limit 1 -outfmt '6 qseqid sacc pident length qstart qend sstart send qlen slen sstrand' > result.blast Ordenamiento de los contigs o scaffolds: icmapper.pl –f result.blast –m 100 –d 7000 –g 10000 donde –f indica el nombre del fichero de salida del alineamiento Blast entre el borrador de genoma y la referencia, -m indica el alineamiento mínimo a considerar, -d y –g respectivamente indican la separación máxima de las diagonales donde se ubican los HSP y la longitud máxima de los huecos entre los HSP (debajo de estos dos parámetros, los HSP se pueden considerar dentro del mismo cluster). Este programa se utilizó en el ensamblaje del genoma de Solea senegalensis donde nos fue útil para confirmar la contigüidad de los contigs solapantes y los scaffolds unificados tomando como referencia el genoma de Cynoglossus semilaevis . También se utilizó para construir los super-scaffolds de Solea senegalensis sobre la base de los cromosomas de Cynoglossus semilaevis III.4. Datos biológicos III.4.1. Secuencias transcriptómicas III.4.1.1. Para el transcriptoma de Solea senegalensis Los datos utilizados consistieron en lecturas largas de 454/Roche y lecturas cortas de Illumina. 84 MATERIALES Y MÉTODOS Lecturas 454/Roche: 14 librerías simples, entre ellas 12 corresponden a 3 zonas del organismo que son testículo/ovario, hipófisis y el hipotálamo, una al sistema inmunitario y una al sistema de osmoregulación. El número total de lecturas fue 5 663 225 con una longitud media de 757 nt Lecturas Illumina: Esta lecturas provienen de experimentos de RNAseq que incluyen a un experimento temporal donde se obtuvieron muestras en los diferentes estadios de desarrollo de S. senegalensis (huevos, embriones, larvas, y estadio metamórfico) y otros experimentos de incubación en diferentes sustancias: AR (ácido retinoico), ATRA (all-trans AR), DEAB (4-diethyl amino-benzaldehyde), DMSO (dimethyl sulfoxide) y TTNPB (agonista de rars) . Se seleccionó una muestra de cada condición experimental para formar un total de 37 muestras que corresponden a lecturas cortas pareadas de Illumina de longitud inicial 2 x 76 nt. III.4.1.2. Para el transcriptoma de Solea solea Para este experimento las lecturas utilizadas son exclusivamente de Illumina obtenidas de un experimento de RNAseq donde se usaron tratamientos análogos a los del experimento de RNAseq de S. senegalensis (apartado anterior) y consisten en 43 librerías de lecturas tipo pareadas de longitud inicial 2 x 100 nt. III.4.1.3. Para el transcriptoma de Tisochrysis lutea Se utilizaron tanto datos de 454/Roche como de Illumina, y las características de las lecturas por cada librería se muestran en tabla III.2. Tabla III.2: Características de las lecturas brutas utilizadas en el ensamblaje del transcriptoma de Tisochrysis lutea Tecnología Índice de la librería Nº de librerías Normalización Tipo Tamaño de inserto (nt) Longitud de lecturas (nt) Nº de lecturas Illumina 1 3 Si PE 150 2 x 100 925 682 496 2 7 No PE 150 2 x 76 612 903 376 454 1 2 Si SE - 547 1 151 067 2 1 No SE - 505 781 III.4.1.4. Para el transcriptoma de Ruditapes decussatus Se utilizó una sola librería de 127 327 692 lecturas Illumina pareadas de tipo PE y de longitud inicial 2 x 75 nt. 85 MATERIALES Y MÉTODOS III.4.2. Secuencias genómicas III.4.2.1. Secuencias genómicas de Photobacterium damselae subsp. piscicida Se utilizaron exclusivamente secuencias de 454/Roche (tabla III.3) Tabla III.3: Características de las lecturas brutas utilizadas en el ensamblaje de Photobacterium damselae subsp. piscicida Cepa Tecnología Índice de las librerías Nº de librerías Tipo Tamaño de inserto Longitud de lecturas (nt) Nº de lecturas L091106-03H 454 1 1 MP 8 kb 509 148 622 2 1 SE - 1 195 297 269 DI21 454 1 2 MP 8 kb 445 433 717 2 1 SE - 550 187 433 III.4.2.2. Secuencias genómicas de Solea senegalensis En este ensamblaje se utilizaron tanto lecturas de 454/Roche como de Illumina, las características de las diferentes librerías utilizadas se ilustran en la tabla III.4 Tabla III.4: Características de las lecturas brutas utilizadas en el ensamblaje del genoma de Solea senegalensis Tecnología Índice de la s librerías Nº de librerías Tipo Sexo Tamaño de inserto Longitud de lecturas (nt) Nº de lecturas Illumina 1 27 PE Hembras 350 2 x 100 2 x 152 869 031 542 2 6 PE Hembras 300 2 x 100 1 005 527 592 3 12 MP Hembras 3 kb 2 x 100 1 109 021 114 454 1 4 MP Hembras 3 kb 551 2 103 984 2 9 MP Hembras 8 kb 546 5 236 315 3 1 MP Hembras 20 kb 579 1 018 187 86 RESULTADOS Y DISCUSIÓN 87 RESULTADOS Y DISCUSIÓN Parte IV Resultados y discusión 94 RESULTADOS Y DISCUSIÓN IV.1.2.2. Ensamblaje de genomas A continuación se evaluó la capacidad de los 4 programas del apartado anterior para ensamblar lecturas artificiales genómicas, y también se añadió el programa CABOG, específico para ensamblar genomas de más de 100 Kb. La evaluación se llevó a cabo con el protocolo experimental ilustrado de la figura IV.2, muy similar al utilizado en la parte de transcriptómica, con la diferencia que la valoración de las características del ensamblaje utilizaba diferentes criterios que son el N50, N90 y la coherencia entre los contigs y la secuencia original. El N50 es un valor definido como la longitud del contig (siempre partiendo desde el contig de mayor tamaño) a partir del cual está contenido el 50% de los nucleótidos ensamblados [140], en el caso de N90 es 90%. N50 y N90 son parámetros muy utilizados para evaluar la calidad de los borradores de genomas ya que dan una idea sobre la longitud de los contigs respecto a la longitud total. Las lecturas artificiales simples se crearon en puntos aleatorios del genoma y con longitudes aleatorias entre 300 y 500 pb para que sean lo más parecidas a las lecturas reales de 454. Las lecturas extraídas se dejaron como exactas (sin errores) para simplificar los ensamblajes y dejar como única variante la complejidad del genoma del organismo, y para permitir una comprobación precisa de los contigs resultantes. En la figura IV.3 se muestra una comparación entre la distribución de las frecuencias de las coberturas de lecturas reales (SRR097627) de Saccharomyces cerevisiae con una cobertura aparente de 4X mapeadas sobre el cromosoma III de esta especie y lecturas artificiales extraídas desde el cromosoma III con la misma cobertura. Según vemos, las lecturas artificiales tienen una distribución normal de frecuencias de coberturas muy parecida a aquella en las lecturas reales. Por tanto, el criterio con el que se generan las lecturas artificiales simula en buena medida a un experimento de secuenciación. Las muestras del genoma utilizadas incluyen todo o parte del genoma de 3 organismos con longitud y grado de complejidad creciente (tabla IV.4) para ver el comportamiento de los programas de ensamblaje frente a diferentes tipos de genomas y evaluar su capacidad a ensamblar las zonas repetitivas. 95 RESULTADOS Y DISCUSIÓN Figura IV.2: Las diferentes etapas de comprobación de un programa de ensamblaje de lecturas genómicas. Figura IV.3: Distribución de la cobertura de lecturas reales y artificiales cuando la cobertura es de 4X sobre el cromosoma III de levadura Tabla IV.4: Muestras de genoma utilizadas para la evaluación de los programas Organismo Secuencia utilizada Longitud Enterobacteria phage λ Todo el genoma 48 502 pb Saccharomyces cerevisiae Cromosoma 3 316 617 pb Homo sapiens Cromosoma 14 (extremo) 500 000 pb 96 RESULTADOS Y DISCUSIÓN En la tabla IV.5 se recogen los resultados de contigs resultantes de los ensamblajes para las secuencias genómicas indicadas en la tabla IV.4. También se indican las características básicas de esos contigs y su calidad (basada en las identidades encontradas por Blast entre el contig y la secuencia original). Así, un contig se considera “bueno” si tiene sea un alineamiento único o varios alineamientos que están en el mismo orden en el contig y en la secuencia original. Para los contig “buenos” se incluye otro índice de calidad, que es el número de huecos ( gaps ) al ser indicadores de la discontinuidad del alineamiento entre el contig y la secuencia original; cuantos menos huecos, mejor es la reconstrucción. Se muestran solo los resultados correspondientes a la cobertura 30X por el mismo motivo que con las pruebas para transcriptómica. Respecto al bacteriófago λ de las Enterobacteria, cuyo genoma es relativamente simple, los ensamblajes fueron rápidos y todos los programas produjeron un único contig cuyo tamaño estaba muy cerca del tamaño del genoma. En todos los casos el contig se alineó con la secuencia original, sin huecos (tabla IV.5-A). CAP3, MIRA y EULER produjeron el resultado más cercano al esperado. En el caso del cromosoma III de Saccharomyces cerevisiae (tabla IV.5-B), cuya secuencia es más larga y que lleva más repeticiones que el bacteriófago, los programas produjeron un borrador de genoma con varios contigs, pero con cobertura alta respecto a la secuencia original (Entre 97,8% en EULER y 99,9% en MIRA). El mal resultado de EULER se debe a que generó el borrador de genoma más fragmentado (49 contigs) con alineamientos discontinuos con la secuencia original (7 huecos, mientras que los otros ensambladores no presentan ninguno). CAP3, MIRA y CABOG generaron un número de contigs más reducido con el menor número por parte de CAP3 (2 contigs). Por tanto, parece que EULER no va a ser un buen ensamblador de secuencias genómicas. Respecto al cromosoma humano (tabla IV.5-C), aunque la cobertura de lecturas fue la misma que en los casos anteriores, los programas se mostraron más conservadores y produjeron más contigs en general. CAP3 generó menos contigs que los demás programas, pero solo 8 de sus 9 contigs fueron “buenos”, y uno de los “buenos” tenía un hueco. En la figura IV.4 se ilustra la estructura del Contig9 generado por CAP3 en comparación con la secuencia original. El orden invertido de los alineamientos entre el contig y la secuencia original y los errores que se observan en las lecturas ensambladas indican que los extremos 3’ y 5’ de esta zona del genoma comparten una secuencia repetitiva que hizo equivocarse al ensamblador. En el caso de VELVET y EULER el borrador de genoma generado fue mucho más fragmentado que los demás (368 y 2 564 contigs respectivamente). Solo 94 (3,7%) de los contigs de EULER fueron superiores a 500 pb. De estos 94 contigs, 2 fueron malos y los demás, a pesar de ser “buenos”, incluían 20 huecos, mientras que VELVET generó 368 contigs, 35 (9,5%) de ellos superiores a 500 pb, y todos fueron “buenos” y sin huecos. MIRA tuvo claramente el mejor resultado con menor número de contigs, y mayor N50 (298 365), además todos sus contigs fueron “buenos”, carecían de huecos y cubrían toda la secuencia del genoma original (100,4%). CABOG tuvo un resultado parecido al de MIRA, aunque con menor N50 (108 539) y menor cobertura de la secuencia original (94,2%). Por tanto, parece 97 RESULTADOS Y DISCUSIÓN que MIRA y CABOG son los ensambladores más idóneos para reconstruir secuencias genómicas de alta complejidad. Tabla IV.5: Resultados de la evaluación de los programas de ensamblaje de novo de secuencias genómicas con una cobertura 30X (A) Enterobacteria phage λ Parte del genoma: Todo Nº de lecturas artificiales (30X): 3 643 Tiempo de ejecución (s) Contigs Contigs “buenos” Nº > 500 nt N50 N90 Nº Suma (nt) Cob. (%)* Huecos CAP3 73 1 1 48 415 - 1 48 415 99,8 MIRA 75 1 1 48 415 1 48 415 99,8 0 CABOG 78 1 1 48 347 - 1 48 347 99,7 0 VELVET 39 1 1 48 409 - 1 48 409 99,8 0 EULER 10 1 1 48 415 - 1 48 415 99,8 0 (B) Saccharomyces cerevisiae Parte del genoma: Cromosoma 3 Nº de lecturas artificiales (30X): 23 789 Tiempo de ejecución (s) Contigs Contigs “buenos” Nº > 500 nt N50 N90 Nº Suma (nt) Cob. (%)* Huecos CAP3 640 2 2 302 907 11 907 2 314 814 99,4 0 MIRA 496 5 5 186 396 23 565 5 316 298 99,9 0 CABOG 528 4 4 185 448 23 205 4 313 165 98,9 0 VELVET 300 21 6 114 172 70 111 6 312 719 98,8 0 EULER 82 49 9 69 034 21 783 9 309 615 97,8 7 (C) Homo sapiens Parte del genoma: últimos 500kb del cromosoma 14 Nº de lecturas artificiales (30X): 37 581 Tiempo de ejecución (s) Contigs Contigs “buenos” Nº > 500 nt N50 N90 Nº Suma (nt) Cov. (%)* Huecos CAP3 1 586 9 9 149 272 33 754 8 416 534 83,3 1 MIRA 1 430 7 7 298 365 193 585 7 502 085 100,4 0 CABOG 677 13 13 108 539 29 820 13 470 800 94,2 0 VELVET 596 368 35 49 352 1 554 35 468 406 93,7 0 EULER 152 2 564 94 1 642 72 92 275 654 55,1 20 (*) Porcentaje del genoma original cubierto por los contigs ensamblados 98 RESULTADOS Y DISCUSIÓN Figura IV.4: Ejemplo de error de ensamblaje detectado en el Contig9 del ensamblaje CAP3 para el fragmento del cromosoma 14 del genoma humano Para comprobar si el comportamiento de los ensambladores podría ser dependiente de la cobertura, analizamos el número de huecos que presentaban los contigs “buenos” producidos en el caso del genoma humano (tabla IV.5-C) para coberturas que van de 5X a 60X (tabla IV.6). En el caso de CAP3, MIRA y EULER se nota un aumento de huecos cuando la cobertura era baja, aunque EULER siempre produjo huecos en todas las condiciones analizadas. Tabla IV.6: Número de huecos en los contigs “buenos” según la cobertura en el caso del fragmento de cromosoma 14 humano Cobertura 5x 10x 20x 30x 40x 50x 60x CAP3 4 1 3 1 2 0 1 MIRA 8 1 1 0 0 0 0 CABOG 0 0 0 0 0 0 0 VELVET 1 1 0 0 0 0 0 EULER 13 25 28 20 12 13 9 En conclusión, los programas de ensamblaje estudiados se mostraron eficientes para ensamblar genomas simples, mientras que en genomas más complejos se ponen de manifiesto mejor sus diferencias. Entre ellos destacan MIRA y CABOG al producir el mejor resultado en cuanto a número de contigs y su calidad. Como recomendación, conviene utilizar MIRA con una cobertura ≥30X. Para ensamblajes con poca cobertura o que ésta sea poco homogénea, es más recomendable utilizar CABOG ya que genera contigs más fiables en estas condiciones. 99 RESULTADOS Y DISCUSIÓN IV.1.2.3. Combinación de ensamblajes de varios programas Puesto que no hay programa bioinformático óptimo para ningún problema biológico, un buen criterio consiste en combinar los resultados de distintos algoritmos con los mismos datos. En el caso de los ensamblajes, parece conveniente combinar los contigs de dos o más programas de ensamblaje para obtener el mejor ensamblaje posible. Ensamblajes de transcriptoma En el apartado IV.1.2.1 hemos visto que los programas MIRA (algoritmo de tipo OLC) y EULER (algoritmo de De Bruijn) generaron los mejores transcriptomas. Así, suponiendo que hay lecturas que se ensamblarán mejor con un algoritmo que con otro, definimos un método para reconciliar los contigs de estos dos ensambladores para sacar un solo transcriptoma. Para este propósito utilizamos CAP3, un programa inicialmente destinado a ensamblar lecturas de tipo Sanger y apto para ensamblar secuencias largas como los contigs resultantes de un ensamblaje de transcritos. Como resultado tendremos supercontigs donde habrá contigs de MIRA y EULER, y singulones que son contigs que no encontraron con quién alinearse. La combinación de estos dos grupos de secuencias esperamos que devuelva un transcriptoma final mejor que el de cualquiera de los ensambladores por separado (figura IV.5). Figura IV.5: Reconciliación entre los contigs de MIRA y EULER utilizando CAP3 La validación de esta estrategia fue realizada en nuestro laboratorio por Noé Fernández Pozo [87] para ensamblar el transcriptoma de Pinus pinaster y presentarlo a la comunidad científica en la base de datos EuroPineDB [141]. En la tabla IV.7 se observa cómo el ensamblaje en supercontigs con CAP3 devuelve un número de singulones más supercontigs menor que la suma de contigs de EULER y MIRA por separado (primera fila de la tabla), y hay mayor número y porcentaje de transcritos completos y únicos. También se ha validado posteriormente en nuestro laboratorio con otros transcriptomas, tanto de especies animales como vegetales [140, 142, 143]. 100 RESULTADOS Y DISCUSIÓN Tabla IV.7: Comparación del ensamblaje de lecturas de 454/Roche para generar el transcriptoma de Pinus pinaster. La tabla está sacada del manuscrito en preparación del programa Full-LengtherNext (Seoane et al, en preparación) Los programas CAP3 y MIRA utilizan el formato ACE para incluir a los datos de las lecturas que han permitido formar la secuencia consenso del contig, mientras que EULER no utiliza este formato. Pensando en que el ensamblaje en formato ACE podría servir para inspeccionar la calidad de ensamblaje de un transcrito (véase más adelante, apartado IV.1.3), y que también podría servir para detectar SNP con GigaBayes [131], decidimos crear un ACE de los supercontigs generados por CAP3, como se ilustra en la figura IV.6. Para ello, se desarrolló un script en Ruby (véase materiales y métodos; apartado III.3.5.2) que importa la información de los dos ficheros ACE y en el fichero ACE del CAP3 reemplaza los contigs del MIRA con las lecturas ensambladas que les corresponden (figura IV.6). Figura IV.6: Combinación de los ficheros ACE de MIRA y CAP3. Los contigs de MIRA (A) se reemplazan con las lecturas que les corresponden (B) para generar un nuevo fichero ACE (C) 101 RESULTADOS Y DISCUSIÓN Ensamblajes de genoma En los ensamblajes de genoma, no es tan fácil hacer una fusión de ensamblajes, ya que puede causar un borrador de genoma con duplicaciones. Además, tiene que hacerse con un programa diferente, GAM-NGS [118], en el que hay que definir un ensamblaje de referencia o maestro, que luego se puede extender con otro ensamblaje que denominan “esclavo”. Comparativemente a otros programas probados como GARM [144], GAM-NGS tiene la ventaja de utilizar la información de mapeo las secuencias de lecturas sobre los contigs o scaffolds para definir regiones llamadas “bloques” que representan el mismo locus entre los diferentes contigs o scaffods , lo cual permite una fusión más fiable de estos. Dado que cada vez se secuencian menos genomas con 454/Roche y más con secuencias cortas de Illumina (cuyos ensambladores incluyen diferentes versiones de algoritmos basados en los grafos de De Bruijn), otra manera de combinar distintos ensamblajes es usar diferentes k -meros y combinarlos. Algunos ensambladores, como SOAPdenovo [66] lo hacen por sí solos, mientras que otros, como Trinity, solo son capaces de usar un k -mer. En los próximos capítulos veremos con más detalle la estrategia seguida para ensamblar el genoma de una bacteria con 454/Roche, y el del lenguado combinando 454/Roche e Illumina. IV.1.2.4. Tamaño real del inserto en las lecturas pareadas Para poder ordenar los contigs obtenidos se utilizan lecturas pareadas, con las que definimos los scaffolds . El tamaño de inserto en las lecturas pareadas (PE o MP) varía de una librería a otra según las plataformas de secuenciación y el protocolo experimental utilizado. La mayoría de los programas de ensamblaje y de scaffolding requieren incluir la media y la desviación típica del tamaño real de los insertos. Para obtener este valor, desarrollamos un script en Ruby (véase materiales y métodos; apartado III.3.5.3) con el que ensamblamos lecturas pareadas que sigue el protocolo ilustrado en la figura IV.7. Las secuencias pareadas usadas para el ensamblaje se mapean con Bowtie sobre los contigs resultantes para medir a qué distancia caen los extremos de las secuencias pareadas y así estimar el tamaño medio y la desviación típica. 102 RESULTADOS Y DISCUSIÓN Figura IV.7: Protocolo utilizado para calcular la media y la desviación típica del tamaño de inserto en lecturas pareadas. El programa usado para el ensamblaje es, en cada caso, el que requiere los valores de tamaño de inserto, como por ejemplo, CABOG IV.1.3. CoMiner: validación de ensamblajes En cualquier tipo de ensamblaje de novo existen errores. La identificación de los ensamblajes erróneos constituye una tarea difícil debido a la gran cantidad de datos y a la variación de la calidad de estos datos por culpa de complicaciones mecánicas y bioquímicas de los secuenciadores. Para la reconstrucción de un transcriptoma o genoma fiable, la corrección de errores en los ensamblajes requiere esfuerzos adicionales de validación manual [145]. La manera habitual de medir la calidad de los ensamblajes se basa en estadísticas simples, como el contig más largo, la longitud media de los contig, la suma de las longitudes y el N50. De hecho, se suele dar preferencia a los contigs largos [146], sin reparar muchas veces en que puedan ser quiméricos o resultantes de un ensamblaje erróneo. Se han desarrollado validadores de ensamblajes basados en parámetros más complejos, como es Hawkeye [147]. Sirve para facilitar la inspección visual de los ensamblajes en el que se muestran los puntos de conflicto para reducir al mínimo el tiempo dedicado a la identificación, edición y corrección. A diferencia de otros editores de contigs como GAP5 [148], Hawkeye combina predictores computacionales con visualización interactiva. Su principal problema está en que es manual, por lo que pierde utilidad a medida que aumenta la cantidad de contigs de un ensamblaje. Por eso, aunque la inspección visual y edición manual sean métodos eficaces, estos pueden constituir unas tareas engorrosas, particularmente para ensamblajes con los datos de secuenciación de nueva generación (NGS), por esta razón se desarrolló un programa de validación automática de los contigs basado en criterios independes, llamado amosvalidate [146]. No obstante, esta herramienta solo marcaba las regiones que aparecían mal ensambladas, y la corrección requiere una nueva visualización y una edición manual utilizando Hawkeye . 103 RESULTADOS Y DISCUSIÓN En este trabajo nos planteamos la detección de errores de ensamblaje y su corrección automática, para lo que desarrollamos un programa llamado CoMiner (verse el apartado III.3.5.8 de materiales y métodos) que fue programado en Ruby y comprobado en un iMac dual core a 3,06 GHz con 4 Gb de RAM. El programa lee y escribe ensamblajes en el formato ACE [149], que es un formato generado por varios programas basados en el algoritmo OLC como los tradicionales Phrap, CAP3 y GAP4-5, y los desarrollados para la NGS de lecturas largas (Newbler, CABOG, Arachne, Minimus y TIGR Assembler). Con el objetivo de incrementar la fiabilidad del ensamblaje sin la intervención humana, el algoritmo de CoMiner se divide en 4 etapas principales (figura IV.8): (i) descubrimiento de las regiones de alta entropía (HER, por High-entropy region ); (ii) identificación de las lecturas conflictivas; (iii) edición de las lecturas (recorte o eliminación); (iv) verificación del contig y escritura del nuevo fichero ACE. Figura IV.8: (tomada de [150]) Flujo de trabajo del algoritmo de CoMiner. A. Descubrimiento de las regiones de alta entropía e identificación de las lecturas conflictivas. B: Edición o eliminación de las lecturas conflictivas en un contig según la distribución de errores. C: Una vez editadas las lecturas conflictivas, se verifica la cobertura del contig. Este se divide en dos o más contigs si es necesario, luego se guarda en el fichero ACE Veamos cada una de estas etapas con más detalle: (i) Descubrimiento de HER: El objetivo de esta etapa es localizar los fragmentos donde las lecturas no alinean perfectamente. Se eligió una medida para definir la calidad de un alineamiento, la entropía de la secuencia consenso en la posición i [- H(i)] tal como se describió en [151]. Se calcula entonces la frecuencia de cada uno de los 4 nucleótidos en cada posición de la secuencia consenso, lo que se utiliza para calcular la entropía en cada posición del consenso. Los SNP y discordancias 110 RESULTADOS Y DISCUSIÓN Figura IV.13: Esquema del flujo de trabajo del programa ICMapper 111 RESULTADOS Y DISCUSIÓN Figura IV.14: Esquema de las diferentes etapas del tratamiento de los HSPs de Blast con ICMapper. A: clusterización de los HSP por diagonal y por longitud de huecos (gaps). B: Formación de alineamientos resumidos por cada clúster. C: Eliminación de las repeticiones y ordenamiento de las secuencias del borrador de genoma 112 RESULTADOS Y DISCUSIÓN Figura IV.15: Dot plot ejemplo de un alineamiento entre contigs genómicos reales de una cepa de la bacteria Lactococcus garvieae y una cepa de Lactococcus lactis (referencia). Antes del tratamiento con ICMapper (A) los fragmentos de alineamiento son pequeños y dispersos además que se nota la existencia de muchos alineamientos repetidos. Después del tratamiento ICMapper (B), los alineamientos resumidos aparecen más largos, de hecho son más representativos del contig ya que cubren zonas de la referencia más grandes, lo que permite un posicionamiento más fiable del contig. De otra parte se eliminaron los alineamientos repetitivos, los cuales probablemente fueron un producto de zonas repetitivas entre los dos genomas Para comprobar las posibilidades de ICMapper, se compararon sus resultados con los de 3 programas del mismo tipo, que son Oslay [154], Projector2 [155] y Mauve Aligner [156]. Los datos utilizados consistieron en borradores de genomas artificiales (contigs) creados a partir del genoma acabado de varias bacterias. En la creación de los contigs artificiales se incluyeron separaciones o solapamientos con longitudes aleatorias para tener una estructura de borrador de genoma parecida a los casos reales. Como referencia se utilizaron cepas de la misma especie o de especies diferentes del mismo género. La evaluación de la precisión del ordenamiento consistió en comparar el orden de los contigs establecido por los programas con su orden real ya conocido. En este estudio, para alinear los contigs artificiales contra el genoma de referencia se utilizó el la versión 2.2.20 de Blast con la siguiente configuración: -F F para desactivar el programa “DUST” de reconocimiento de secuencias de baja complejidad y evitar que se interrumpe el alineamiento Blast en estas zonas; -G 5, para afectar el valor 5 al coste por abrir huecos lo que permite reducir los huecos en los alineamiento, y -E 10 para afectar el valor 10 al coste por extender los huecos lo que permite reducir la longitud de los huecos abiertos. El resto de parámetros se dejaron por defecto. Tal como se observa en la tabla IV.9, ICMapper y Mauve aligner generan mejores resultados que los dos demás programas, En caso de la especie Mycobacterium tuberculosis CDC1551, a la diferencia de los demás programas, ICMapper tuvo un porcentaje de ordenamiento de 100% en los dos casos donde la referencia era una cepa de la misma especie. En general ICMapper tuvo una precisión de ordenamiento superior o igual que Mauve aligner 113 RESULTADOS Y DISCUSIÓN y muy superior que Projector2 y Oslay, excepto el caso de la bacteria Lactococcus lactis subsp, cremoris MG1363, donde ICMapper tuvo una precisión de ordenamiento inferior a Mauve aligner (95.58% y 98.57% respectivamente). Tabla IV.9: Porcentaje de pares de bases totales correctamente ordenadas utilizando ICMapper y tres otros programas. (B: Borrador de genoma, R: Referencia) Suma de longitudes de los contigs correctamente ordenados por cada programa Nº contigs artificiales Suma de bases ICMapper Mauve aligner Projector2 Oslay B: Mycobacterium tuberculosis CDC1551 144 4 311 730 4 311 730 4 307 728 4 292 255 4 244 863 R: Mycobacterium tuberculosis H37Rv 100% 99,9% 99,54% 98,44% B: Mycobacterium tuberculosis CDC1551 144 4 311 730 4 311 730 4 295 761 4 290 321 4 252 754 R: Mycobacterium tuberculosis KZN 100% 99,62% 99,50% 98,63% B: Mycobacterium tuberculosis CDC1551 144 4 311 730 3 942 363 3 919 573 3 028 947 3 549 569 R: Mycobacterium avium 104 91,43% 90,90% 70,24% 82,32% B: Lactococcus lactis subsp, cremoris MG1363 88 2 497 652 2 387 292 2 461 938 2 087 813 2 319 609 R: Lactococcus lactis subsp. lactis Il1403 95,58% 98,57% 83,59% 92,87% B: Streptococcus pneumoniae 670-6B 77 2 197 707 2 122 247 2 122 247 1 993 147 2 111 903 R: Streptococcus pneumoniae G54 96,56% 96,56% 90,69% 96,09% B: Bacillus amyloliquefaciens DSM7 140 3 920 736 3 888 585 3 880 362 3 234 464 3 878 842 R: Bacillus amyloliquefaciens FZB42 99,17% 98,97% 82,49% 98,93% B: Bacillus thuringiensis BMB171 183 5 254 098 5 241 972 5 074 150 4 924 055 5 139 721 R: Bacillus thuringiensis serovar konkukian str. 97-27 99,76% 96,57% 93,71% 97,82% B: Bacillus thuringiensis BMB171 NC_014171 183 5 254 098 5 160 653 5 176 617 4 889 963 5 007 363 R: Bacillus thuringiensis str. Al Hakam 98,22% 98,52% 93,06% 95,3% B: Escherichia coli 536 154 4 795 594 4 667 799 4 616 347 4 023 816 4 474 579 R: Escherichia coli 55989 97,33% 96,26% 83,9% 93,3% Valor medio del porcentaje de los contigs correctamente ordenados 97,6% 97,3% 88,5% 94,9% Por otra parte, se comparó el tiempo de ejecución de ICMapper con Mauve aligner (el programa que generó el resultado más parecido). La tabla IV.10 muestra el tiempo transcurrido durante todas las etapas del proceso (en caso de ICMapper, se incluye a tiempo transcurrido durante la etapa del alineamiento Blast). En todos los casos se observa una diferencia grande entre los tiempos de ejecución de los dos programas. ICMapper llegó a ser 20 veces más rápido que Mauve aligner. Esta diferencia en el tiempo de ejecución reside en el hecho de que Mauve aligner utiliza su propio algoritmo de alineamiento y ejecuta varios ciclos de ordenamiento después de los cuales escoge el mejor resultado basándose en criterios preestablecidos, mientras que ICMapper de una parte aprovecha de la rapidez del algoritmo Blast y de otra parte realiza la tarea de ordenamiento en un ciclo único y eficaz. 114 RESULTADOS Y DISCUSIÓN Tabla IV.10: Comparación entre los tiempos de todas las etapas de ejecución de los programas ICMapper y Mauve aligner Tamaño del genoma Tiempo global de ejecución (s) ICMapper Mauve aligner Mycobacterium tuberculosis CDC1551 4 311 730 27 105 Lactococcus lactis subsp, cremoris MG1363 2 497 652 5 98 Streptococcus pneumoniae 670-6B 2 197 707 5 57 Bacillus amyloliquefaciens DSM7 3 920 736 7 61 Valor medio del tiempo de ejecución (s) 11 80 Contrariamente a Mauve aligner que esta destinado a genomas de tamaño pequeño, especialmente los genomas microbianos [156], ICMapper fue capaz de ordenar tanto borradores de genomas de microorganismos como borradores de genomas grandes de eucariotas. Ej. Genoma de Solea senegalensis (Apartado IV.3.2 a continuación) 115 RESULTADOS Y DISCUSIÓN IV.2. Ensamblajes de transcriptomas IV.2.1. Estrategia tomada como modelo Los ensamblajes y anotaciones del transcriptoma realizados en este trabajo se basaron sobre una flujo de trabajo que se desarrolló en nuestro grupo de investigación [87] para el transcriptoma del pino ( Pinus pinaster; figura IV.16) Preprocesamiento: Las lecturas se preprocesan con SeqTrimNext, desarrollado también en nuestro laboratorio a partir de SeqTrim [152] ( http://www.scbi.uma.es/seqtrimnext/ ; [46]) Ensamblaje: A continuación se describe la estrategia del procesamiento y ensamblaje que se utilizó para el transcriptoma del pino (figura IV.17). Se pueden distinguir dos partes: ensamblaje, propiamente dicho, y reconciliación. El ensamblaje de lecturas largas se realiza con el método descrito en el apartado IV.1.2.3 basado en la combinación de los ensamblajes de MIRA [109] y EULER [110], el cual fue optimizado para extraer la máxima información de las lecturas y que ésta sea fiable. El ensamblaje de lecturas cortas se realiza con el programa ABySS [67] utilizando una serie de k-meros de 43 a 73. Los contigs de este ensamblaje se agrupan con CD-HIT [129] y después TGICL [157] a fin de eliminar las redundancias generadas principalmente por la utilización de múltiples k-meros. La reconciliación final consiste en un ensamblaje CAP3 [105] de los grupos de secuencias útiles procedentes de los ensamblajes anteriores que son los siguientes: (i) Contigs de MIRA (ii) Debris de MIRA que son codificantes (el resto se consideran artefactos de secuenciación) (iii) Contigs de EULER sobre los que mapean lecturas (iv) Contigs de EULER sobre los que no mapean lecturas, pero son codificantes (el resto se consideran artefactos de ensamblaje) (v) Contigs agrupados de ABySS. El ensamblaje da lugar a supercontigs y singulones que se juntarán para formar el transcriptoma final tal como se explicó más detalladamente en el apartado IV.1.2.3. Verificación: Este paso se realiza con Full-LengtherNext [87] que sirve para evaluar de un modo rápido si el ensamblaje tiene una buena proporción de genes completos, si está muy fragmentado o si posiblemente contiene un alto número de secuencias sin información 116 RESULTADOS Y DISCUSIÓN biológica. Esta verificación puede utilizarse tanto para evaluar la calidad del transcriptoma final como para evaluar la calidad de los ensamblajes de lecturas para poder optimizarlos con el ajuste de los parámetros de ensamblaje. Full-LengtherNext, después de las ultimas mejoras (Seoane et al., en preparación) tiene la posibilidad de crear un transcriptoma representativo que incluye los genes anotados únicos y los codificantes putativos. De este modo se reduce la fragmentación del transcriptoma que después resulta más manejable y se mejora la calidad de los resultados en los estudios posteriores (por ejemplo de RNAseq). Figura IV.16: (tomada de [87]): Flujo de trabajo propuesto para obtener un transcriptoma de pino anotado a partir de lecturas de NGS y de tipo Sanger. Las secuencias de tipo Sanger se preprocesan con SeqTrim y las de NGS con SeqTrimNext. Todas las lecturas limpias se ensamblan juntas con MIRA, y las secuencias limpias de NGS se ensamblan además con EulerSR, pero sin mezclar los diferentes experimentos. Luego se re-ensamblan los contigs de Euler-SR y MIRA utilizando CAP3, posteriormente se verifican los unigenes del ensamblaje con FullLengtherNext y si se considera que el ensamblaje es correcto, se anotan los unigenes obtenidos con AutoFact, Blast2GO, Gigabayes y MREPS Anotación: La anotación que se propuso en el flujo (figura IV.16) se completa con los programas AutoFact, Blast2GO, Gigabayes y MREPS. De este modo, los transcritos quedan anotados con tres descripciones del producto del gen, que permiten a los investigadores confirmar si 117 RESULTADOS Y DISCUSIÓN tres programas que usan distintos criterios para anotación coinciden en la información encontrada. También se añaden términos GO, códigos EC, InterPro, rutas KEGG, SNP y SSR que permitirán realizar estudios computerizados de los resultados con programas de enriquecimiento biológico o detección de posibles marcadores moleculares para mejora genética. Figura IV.17: (tomada de [158]) Esquema del flujo de trabajo del procesamiento, ensamblaje y reconciliación de las lecturas 454 y Illumina de Pinus pinaster Generación de una base de datos: El primer prototipo de base de datos transcriptómica desarrollada por nuestro laboratorio fue EuroPineDB [141], que poco después fue mejorada y ampliada en SustainPine [158] (figura IV.18). Están basadas en Ruby On Rails tal y como se describe en material y métodos (apartado I.3.4.1). La principal mejora entre las dos está en que la web, la base de datos propiamente dicha, y los cálculos se alojan en una estructura de 3 máquinas virtuales [46] para agilizar su uso. También se mejoró el sistema de paginación interna y la forma de navegar por el contenido a base de tabuladores interactivos. 118 RESULTADOS Y DISCUSIÓN Figura IV.18: Base de datos de SustainPineDB en su estado inicial. A: Arquitectura de la base de datos. B: Página de inicio. C: página de anotaciones 119 RESULTADOS Y DISCUSIÓN IV.2.2. Transcriptomas de Solea senegalensis y Solea solea El lenguado senegalés ( Solea senegalensis ) y el lenguado común ( Solea solea ) son dos especies de peces planos importantes desde el punto de vista económico y evolutivo, tanto en la pesca como en la acuicultura. A pesar de la disponibilidad de algunos recursos genómicos que se describieron recientemente, se necesitan aún más experimentos de secuenciación para establecer un transcriptoma completo y representativo. Por otra parte, un análisis comparativo del transcriptoma de ambas especies puede ayudar a entender la evolución de pez plano. IV.2.2.1. Preprocesamiento de las librerías Un total 37 y 43 librerías de Illumina de Solea senegalensis y Solea solea respectivamente (1 800 millones de lecturas) y 14 librerías adicionales de Roche/454 para Solea senegalensis (5,6 millones de lecturas) fueron preparadas en el laboratorio del grupo de investigación de IFAPA El Toruño (Cádiz). En la tabla IV.11 se presenta un resumen del preprocesamiento con SeqTrimNext de las lecturas Roche/454 e Illumina de todas las librerías que se pusieron a nuestra disposición. Se observa que la mayoría de las lecturas pareadas fueron útiles (83,3% y 79,5% en S. senegalensis y S. solea respectivamente). La fuente de contaminaciones más importantes fue de secuencias ribosómicas y mitocondriales. Otros contaminantes menos representados en los datos de Illumina fueron secuencias de los organismos utilizados para la nutrición de las larvas (artemias (8-21%) y rotiferos (2-4%), principalmente) y otros microorganismos diferentes, principalmente de hongos. Tabla IV.11: Resumen del pre-procesamiento de las lecturas originales de los transcriptomas de S. senegalensis y S. solea 126 RESULTADOS Y DISCUSIÓN Figura IV.24: Flujo de trabajo de la anotación de los transcriptomas de lenguado IV.2.2.4. SoleaDB, una base de datos para explorar los transcriptomas de Solea SoleaDB fue construida sobre la base de SustainPineDB para almacenar los distintos ensamblajes de los transcriptomas de S. senegalensis y S. solea , y los diferentes tipos de anotaciones. En la figura IV.25 se puede observar el esquema de las tablas de la base de datos en SoleaDB y en la figura IV.26 se muestran unas capturas de pantalla de la interfaz web. Las principales mejoras de SustainPineDB aplicadas en SoleaDB fueron las siguientes: (i) Mejora en el modo de acceso a los datos de anotaciones. En efecto, ahora se puede acceder a las anotaciones específicas a cada versión del transcriptoma en una estructura con pestañas horizontales que permiten visualizar de forma independiente la información relativa a cada versión del transcriptoma (figura IV.26-B). (ii) Se añadieron estadísticas completas sobre cada ensamblaje y sus anotaciones, la cual se puede consultar en la página “Assembly info”. (iii) En la lista de transcritos, se añadió la posibilidad de buscar por nombre de transcritos o por anotaciones y de otra parte se implementó una función que permite resaltar las palabras claves comunes entre las descripciones de los 3 programas de anotación, lo que permite al usuario de localizar fácilmente las descripciones comunes. (iv) En la página donde se presenta la información sobre el transcrito, se añadió un campo de curación donde se pueden introducir notas o información sobre el transcrito en cuestión. (v) En la página de la información sobre el ensamblaje, se implementó una herramienta que permite extraer toda la 127 RESULTADOS Y DISCUSIÓN Figura IV.25: Esquema de las tablas de base de datos SoleaDB que guardan las misma estructura utilizada en SustaPineDB con algunos elementos adicionales 128 RESULTADOS Y DISCUSIÓN información útil sobre una lista de transcritos y exportarla a un fichero en formato tabulado, el cual se puede abrir con programas como Excel (figura IV.26-C). (vi) En la página de los KEGG se incorporó la posibilidad de visualizar las mapas de las rutas metabólicas en las cuales se resaltan en color verde los códigos ECs propios de la especie estudiada. En la figura IV.27 se puede observar un ejemplo de los mapas KEGG a los que se puede acceder desde la web. 129 RESULTADOS Y DISCUSIÓN Figura IV.26: Capturas de pantalla de la interfaz de la base de datos SoleaDB. A: Portada de la web. B: Ilustración de la pestaña “Assemblies” en la que se muestra toda la información sobre todas las versiones de transcriptomas y subversiones. C: captura de la parte de pantalla que corresponde a la pestaña “Assembly info” donde hay información general sobre el transcriptoma, así como ficheros descargables y herramientas 130 RESULTADOS Y DISCUSIÓN Figura IV.27: Mapa de KEGG de la ruta del metabolismo de alanina, aspartato y glutamato, en el que las enzimas anotadas en S. senegalensis aparecen resaltadas en verde. IV.2.2.5. Calidad del transcriptoma de los dos lenguados Un análisis detallado del resultado de anotación Full-LengtherNext (tabla IV.14) revela primero que, a pesar del número muy elevado de los transcritos, el número de artefactos fue muy bajo. Segundo, que el transcriptoma de S. senegalensis v4 mejoró con respecto a la v3 en términos longitud de los transcritos y número de los ORF completos, aunque ambas versiones presentaron un número similar de ID diferentes de ortólogos, lo que muestra la importancia de utilizar muchos tejidos para la preparación de las librerías Roche/454 a la hora de obtener una mayor representación de los genes. Tercero, el número de los ORFs completos diferentes fue alto y similar entre los transcriptomas de las dos especies (tabla IV.14 “Diferent complete ORFs”), poniendo de manifiesto que los transcriptomas de lenguado son igual de fiables. Finalmente, S. senegalensis v4 tuvo menor número de ID de ortólogos y menor número de ORF completos diferentes que S. solea v1 lo que puede ser debido a una mayor 131 RESULTADOS Y DISCUSIÓN fragmentación en S. senegalensis ya que, según parece, los datos de Roche/454 solo contribuyeron en incrementar la ortología en el ensamblaje final. Tabla IV.14: (tomada de [143]) : Características de los transcriptomas ensamblados de S. senegalensis y S. solea El alto número de lecturas utilizadas en el ensamblaje de los dos transcriptomas de lenguado puede haber favorecido la acumulación de errores [162, 165]. La evaluación de la fiabilidad del transcrito se basó inicialmente en el mapeo de las lecturas útiles de dos librerías aleatorias de cada transcriptoma utilizando Bowtie2 [119]. Como el 96,7%-98% de las lecturas mapearon sobre los transcritos, los errores de ensamblaje pueden considerarse insignificantes. Curiosamente, el transcrito más largo de los transcriptomas de S. senegalensis v4 y S. solea v1 no es un artefacto, ya que en ambas especies corresponde a una proteína parecida a la titina cuyo ARNm es muy largo (94 446 bp) y que anteriormente ya se ha ensamblado en la tilapia (Número de acceso XM_005460929). El hecho de que este transcrito es 6 veces más largo en S. senegalensis v4 que en S. senegalensis v3 supone la contribución importante de las lecturas cortas en el ensamblaje final. Los transcritos sin ortólogo como fuente de nuevos transcritos de lenguado. El alto número de transcritos sin ortólogo (tabla IV.14) mereció un análisis más profundo. En función del índice de testcode [135], más del 91% de los transcritos son desconocidos (“Unknown”), lo que puede explicar en buena parte la falta de ortología. Para 132 RESULTADOS Y DISCUSIÓN verificar el origen de estos transcritos desconocidos, se mapearon lecturas genómicas (desde varias librerías genómicas de S. senegalensis disponibles en nuestro laboratorio) sobre los transcritos desconocidos de S. senegalensis v4, lo que dio lugar a 462 568 (91,25%) transcritos desconocidos mapeados (440 385 con más de 10 lecturas). Este alto porcentaje de mapeo indica que estas secuencias no fueron artefactos de ensamblaje y que pueden corresponder a fragmentos genómicos o transcritos inmaduros que copurifican con los transcritos maduros. Por otra parte, entre 7,30% y 8,37% de los transcritos sin ortólogo en S. senegalensis v4 y S. solea v1, respectivamente, algunos mostraron un índice Testcode >0,94 (tabla IV.14, “Putative new transcripts”), por lo que tienen alta probabilidad de ser transcritos codificantes. La fracción no redundante de estos transcritos (14 451 y 15 503 transcritos en S. senegalensis v4 y S. solea v1, respectivamente; tabla IV.14), sobre la base de la ausencia de ortólogos en la base de datos UniProtKB, salvo que sean contaminantes de especies no secuenciadas aún, pueden considerarse proteínas (o fragmentos de proteínas) “nuevas” de lenguado. Transcriptoma de referencia en cada especie de lenguado El número alto de transcritos ensamblados indicó una sobreestimación de los transcriptomas de lenguado cuando se comparan con otros teleósteos [159, 160]. Probablemente, la sobrevaloración venga de los transcritos que pueden representar alelos, parálogos, transcritos fragmentados, isoformas o secuencias alternativas de ayuste, ARNm inmaduros o incluso una combinación de ellos. Por lo tanto, para los estudios de estos transcriptomas conviene utilizar el transcriptoma de referencia generado por FullLengtherNext que incluye un número de transcritos mucho más reducido y cercano al teórico. El transcriptoma de referencia para S. senegalensis ( S. senegalensis v4.1) consistió en 59 514 transcritos y para S. solea ( S. solea v1.1) consistió en 54 005 transcritos (tabla IV.14, última fila) Para verificar la representatividad de los transcritos que constituyen los transcriptomas de referencia, se mapearon las lecturas útiles (tabla IV.11) sobre los dos transcriptomas de referencia. Como resultado, el 82,3-87,5% de las lecturas mapearon sobre los transcritos, mientras que el 76,5-93,3% de los transcritos recibieron el mapeo de más de una lectura, lo que muestra que estos transcritos representan adecuadamente el transcriptoma. Otra verificación del transcriptoma de referencia (v4.1 y v1.1) se basó en un análisis de ortología con el pez cebra (43 132 entradas disponibles en RefSeq y 42 555 en ENSEMBL) utilizando Top-blast. En S. senegalensis v4.1, 39 851 transcritos de referencia fueron ortólogos a 21 542 entradas RefSeq y 20 753 entradas Ensembl de pez cebra (figura IV.28, izquierda). En S. solea v1.1, 34 949 transcritos de referencia fueron ortólogos a 20 594 entradas RefSeq y 19 632 entradas Ensembl de pez cebra (figura IV.28, derecha). Estos números muestran que un cierto número de alelos, ARNm inmaduros y genes específicos de linaje (incluso algunos transcritos quiméricos no detectados) pudieron haber sido incluidos en el transcriptoma de referencia. Por otra parte, como el número de ID de RefSeq y Ensembl corresponde casi a la mitad de los transcritos de lenguado, es probable que se hayan incluido ambos alelos de cada gen en el transcriptoma de referencia. Esta hipótesis también está soportada por el hecho de que las muestras analizadas corresponden a animales silvestres y a larvas, que son 133 RESULTADOS Y DISCUSIÓN mayoritariamente heterocigotos. Es importante señalar que el número de los ID diferentes de pez cebra en RefSeq o Ensembl es cerca de los 21 000 reportados recientemente en otro lenguado ( Cynoglossus semilaevis ) y no son tan diferentes de los 26 206 genes que fueron recientemente reportados en pez cebra [159, 166]. Entonces, se puede sugerir que los transcriptomas de referencia de lenguado desarrollados en este trabajo cubren la mayoría de los genes del lenguado. Figura IV.28: Anotación de los transcriptomas de referencia con los ortólogos de pez cebra utilizando los ID de RefSeq y Ensembl Sse: número de transcritos de Solea senegalensis con un ID; Sso: número de transcritos de Solea solea con un ID. R: número de ID únicos de RefSeq; E: número de ID únicos de Ensembl IV.2.2.6. Análisis comparativo entre las dos especies de lenguado Un análisis comparativo entre los transcriptomas de las dos especies de lenguado puede revelar nuevas pistas sobre su biología y evolución y proporcionar otros elementos para confirmar la representatividad de su transcriptoma de referencia. S. solea y S. senegalensis muestran una clara similitud funcional Para comprobar la similitud de los dos lenguados desde el punto de vista funcional, comparamos la abundancia de los términos GO de los transcritos con ortólogo de la tabla IV.14. La distribución de los términos GO por categorías entre las dos especies de lenguado reveló que eran similares (figura IV.29). El número más alto de transcritos anotados con procesos biológicos fue asociado con procesos metabólicos (15,2%) y celulares (22,2%) (figura IV.29-A). En componentes celulares, las categorías más representadas fueron células (36,3%) y orgánulos (22,1%) (figura IV.29-B). En función molecular, el mayor número de transcritos anotados fue en la categoría de actividad catalítica (30,4%) (figura IV.29-C). Curiosamente, las actividades reguladoras de canales y antioxidantes fueron representadas solo en S. senegalensis. La ausencia de estas dos actividades en S. solea seguramente se origina de la secuenciación de las muestras o del proceso del ensamblaje, sobre todo que están muy poco representadas en S. senegalensis , aunque no podemos descartar la existencia de una 134 RESULTADOS Y DISCUSIÓN Figura IV.29: (tomada de [143]) distribución de los Gos en los dos transcriptomas según los procesos biológicos (A), componentes celulares (B) y función molecular (C) 135 RESULTADOS Y DISCUSIÓN razón biológica. En conclusión, los transcriptomas de lenguado son aparentemente similares del punto de vista biológico y funcional. S. solea y S. senegalensis muestran una alta similitud de los genes Otra comparación de los transcriptomas de lenguado se basó en su parecido con una tercera especie que hace de referencia, o sea, en la ortología con pez cebra. La figura IV.30 muestra que el 78,4% de los transcritos ortólogos en los lenguados con un ortólogo en pez cebra tuvieron una similitud ≥95% a nivel de nucleótido, y 1 437 de ellos eran idénticos a 100%. Como se esperaba, los transcritos sin ortología con pez cebra tuvieron un grado de identidad menor (92-96%; figura IV.30). Hay que señalar que 41 transcritos sin ortólogo tuvieron la misma secuencia y 49 transcritos mostraron 99% de identidad entre las dos especies de lenguado. Todos estos datos demuestran el alto nivel de similitud entre los dos transcriptomas de lenguado. Suponiendo que los transcritos de lenguado que comparten el mismo ortólogo de pez cebra también pueden considerarse como ortólogos de lenguado, un total de 39 851 transcritos de referencia en S. senegalensis y 34 949 transcritos de referencia en S. solea comparten 17 562 IDs para RefSeq y 17 031 IDs para Ensembl, lo cual refleja claramente el alto nivel de ortología entre estas dos especies de lenguado, y de nuevo sugiere que en el transcriptoma de referencia están incluidos los dos alelos de cada gen (al haber aproximadamente el doble de transcritos de lenguado que ID de pez cebra). Figura IV.30: (tomada de [143]): Distribución de los niveles de similitud entre los transcriptomas de referencia de los dos lenguados para los transcritos con (barras negras) o sin (barras gris) ortólogo de pez cebra S. solea y S. senegalensis contienen transcritos específicos del lenguado y de peces planos. Los ortólogos verdaderos entre ambas especies de lenguado se obtuvieron después de realizar un Blast reciproco entre los transcriptomas de referencia. En este análisis, dos 142 RESULTADOS Y DISCUSIÓN IV.2.3. Transcriptoma de la microalga Tisochrysis lutea El potencial económico de las microalgas está en crecimiento continuo [176]. Sus aplicaciones son múltiples y van desde la purificación de las aguas hasta la alimentación animal y humana [177]. Isochrysis galbana recientemente llamada Tisochrysis lutea (Tiso) [178] es una microalga muy utilizada en el sector de acuacultura como fuente de proteínas y vitaminas. Por eso se utiliza en la dieta alimentaria de las larvas y la fase juvenil de los peces. Estudios sobre el transcriptoma permitirían ampliar la investigación sobre varios aspectos como el metabolismo y el ciclo de vida de esta especie. IV.2.3.1. Preprocesamiento de las librerías Las muestras de Tisochrysis lutea secuenciadas se prepararon en el laboratorio del grupo de investigación de IFAPA El Toruño (Cádiz). El resumen del preprocesamiento con SeqTrimNext de las lecturas normalizadas y sin normalizar obtenidas con las tecnologías Roche/454 e Illumina (figura IV.33) se presenta en la tabla IV.17. Los contaminantes más representados fueron el plastidio, el RNA ribosómico y otros microorganismos, especialmente hongos como Aspergillus niger y Schizosaccharomyces pombe. El porcentaje de lecturas no normalizadas útiles fue alto para Roche/454 y muy alto para Illumina (69,9% y 93,0% respectivamente) mientras que en las lecturas normalizadas estos porcentajes fueron bastante bajos tanto para Roche/454 como para Illumina (32,4% y 31,8% respectivamente). Estos bajos porcentajes de lecturas normalizadas útiles se deben al exceso de secuencias contaminadas (52,3% en Roche/454 y 51% en Illumina) de las cuales, la mayor parte fueron lecturas del plastidio. 143 RESULTADOS Y DISCUSIÓN Figura IV.33: Estrategia de pre-procesamiento, ensamblaje y la reconciliación del transcriptoma de Tisochrysis lutea Tabla IV.17: Resumen del preprocesamiento de las lecturas originales del transcriptoma de T. lutea Referencia a la figura IV.33 Roche/454 Illumina Total lecturas #1 Normalizadas 1 151 067 925 682 496 No normalizadas 505 781 612 903 376 Longitud media Normalizadas 547 100 No normalizadas 76 Lecturas rechazadas #2 Normalizadas 602 104 (66,2%) 628 141 506 (68%) No normalizadas 48 828 (25,7%) 34 565 829 (5,6%) Contaminación Normalizadas 602 104 (52,3%) 472 474 387 (51%) No normalizadas 48 828 (9,7%) 21 926 541 (3,6%) Total de lecturas útiles #3 Normalizadas 373 149 (32,4%) 294 921 874 (31,8%) No normalizadas 353 784 (69,9%) 569 737 072 (93,0%) Lecturas pareadas Normalizadas - 205 762 484 (22,2%) No normalizadas - 555 513 714 (90,9%) Lecturas simples Normalizadas - 89 159 390 (9,6%) No normalizadas - 14 223 358 (2,3%) Longitud media Normalizadas 233 84 No normalizadas 67 144 RESULTADOS Y DISCUSIÓN IV.2.3.2. Ensamblaje El transcriptoma de T. lutea fue inicialmente ensamblado siendo el mismo protocolo utilizado para el transcriptoma de Solea senegalensis (apartado IV.2.2.2) donde se emplearon las lecturas Roche/454 y las lecturas Illumina normalizadas. Con esto se obtuvo la versión 1 del transcriptoma de la microalga. Para obtener más cobertura del transcriptoma y más genes completos, se realizó otro ensamblaje donde se añadieron las lecturas no normalizadas (RNA-seq). En este segundo ensamblaje se aplicó el mismo protocolo utilizado para ensamblar el transcriptoma de S. senegalensis (Apartado IV.2.2.2), excepto en la parte del ensamblaje de las lecturas cortas donde Oases fue reemplazado por SOAPdenovo-Trans [179] ya que este era más rápido para ensamblar grandes cantidades de lecturas [161], como en este caso. De hecho SOAPdenovo-Trans fue capaz de ensamblar todas las librerías Illumina (normalizadas y no normalizadas) juntas en menos de 1 día, lo cual no fue posible con Oases utilizando el k-mero 23 debido a que el proceso de ensamblaje fue interrumpido al consumir el tiempo máximo asignado a los procesos encolados (una semana) (figura IV.33). En la tabla IV.18 se muestra el número de secuencias en cada etapa del protocolo del ensamblaje así como las características del transcriptoma generado. El número total de transcritos fue 149 517 con una longitud media de 499. Tabla IV.18: Número de secuencias en cada etapa del ensamblaje de T. lutea y las características de la versión 1 del transcriptoma Referencia a la figura IV.33 454 Contigs de MIRA #4 27 668 Contigs de EULER Mapeados #6 12 067 No mapeados codificantes #7 71 Contigs totales 454 #8 39 806 Illumina Contigs de SOAPdenovo #10 191 682 Secuencias útiles #11 190 473 Reconciliación CAP3 #9 149 517 145 RESULTADOS Y DISCUSIÓN Limpieza adicional aplicada sobre de los transcritos. Aunque las muestras secuenciadas procedían de un cultivo homogéneo (poca variabilidad genética), el número de transcritos fue muy superior al número de proteínas de la haptofita más cercana Emiliania huxleyi (39 125). Esto nos hizo sospechar que los transcritos podrían incluir secuencias que no eran de T. lutea . Por lo tanto, realizamos una nueva búsqueda de contaminantes que consistió, en un primer lugar, en alinear los transcritos contra el genoma humano ya que este por su tamaño muy grande no estaba incluido en la base de datos de contaminantes de SeqTrimNext. Un total de 37 130 secuencias se alinearon con el genoma de humano y por lo tanto fueron descartadas. En un segundo lugar, el resto de los transcritos fueron comparados contra la colección de nucleótidos ( Nucleotide collection nr/nt ) de NCBI utilizando blastn ( http://blast.ncbi.nlm.nih.gov/Blast.cgi ), lo cual permitió descartar otras 42 082 secuencias de las cuales algunas fueron de la bacteria Kordia algicida, pero la mayoría fueron de plantas . Después de esta segunda “limpieza” de contaminantes que, según nos constó, procedía de las librerías de Illumina sin normalizar, el número final de transcritos se quedó en 70 306, los cuales formaron la versión 2 del transcriptoma. En la tabla IV.19 se muestra una comparación entre las dos versiones 1 y 2 de T. lueta. Se nota que aunque el número de transcritos se redujo a más de la mitad, el número de transcritos con tamaño superiores a 500 nt permaneció casi invariable, por lo que la mayoría de los transcritos eliminados como contaminantes desde la versión 1 en su totalidad eran pequeños (<500 nt). La eliminación de estos elementos contaminantes hizo que la longitud media de los transcritos pasara de 499 pb a 904 pb y el N50 de 1447 a 1646. Tabla IV.19: Características de las dos versiones del transcriptoma de T. lutea T. lutea v1 T. lutea v2 Transcritos 149 517 70 306 Nº transcritos > 500 nt 37 417 36 385 Longitud media 499 904 N50 1447 1646 Transcrito más largo 10 309 10 309 IV.2.3.3. Anotación del transcriptoma e inclusión en base de datos En la anotación del transcriptoma de T. lutea (v1 y v2) se aplicó el mismo flujo de trabajo utilizado en Solea senegalensis (apartado IV.2.2.3). Los datos de anotación resultantes fueron incluidos en la base de datos IsochrysisDB ( http://www.scbi.uma.es/isochrysisdb/ ) cuya portada se muestra en la figura IV.34. La estructura y la interfaz de la base de datos se basó en SoleaDB, la de Solea senegalensis (apartado IV.2.2.4). 146 RESULTADOS Y DISCUSIÓN Figura IV.34: Portada de la base de datos IsochrysisDB IV.2.3.4. Análisis del transcriptoma de Tisochrysis lutea En la figura IV.35 se muestra la distribución de longitudes de los transcritos en la versión 2 del transcriptoma. Se nota una abundancia de los transcritos con longitudes entre 100 y 300 nt, En el margen de longitudes superiores a 300 y hasta 2 000 nt, el número de transcritos fue bastante uniforme. Figura IV.35: distribución de las longitudes de los transcritos de T. lutea v2 147 RESULTADOS Y DISCUSIÓN En la tabla IV.20 se muestra el resumen del análisis con Full-LengtherNext del transcriptoma de T. lutea v2. Se nota que la gran mayoría de transcritos fueron validos, ya que los artefactos formaron solo 4,7% del transcriptoma. El porcentaje de los transcritos anotados con un ID de ortólogo fue 34%, de estos IDs de ortólogos 47,6% fueron únicos y 23,2% fueron completos únicos. La parte más grande de los transcritos fueron sin ortólogo (65,2%), aunque una buena parte de ellos (42,9%) fueron codificantes por lo que representaron los nuevos transcritos o bien formas inmaduras y truncadas de otros transcritos conocidos (resultados no mostrados). Después de la eliminación los artefactos y la resolución de las quimeras, el transcriptoma global se quedó con 69 800 transcritos. De estos, FullLengtherNext seleccionó 23 671 transcritos que formaron el transcriptoma representativo. Tabla IV.20: Características del transcriptoma ensamblado de T. lutea v2 Transcritos % Nº de transcritos 70 306 Artefactos 3 339 4,7 Transcritos validos 70 306 100 >500pb 35 554 50,6 >200pb 50 035 71,2 Transcrito más largo 10 309 - Transcritos con ortólogo 23 939 34,0 Diferentes ortólogos IDs 11 392 47,6 ORFs completos 9 199 38,4 ORFs completos diferentes 5 546 23,2 C-terminus 6 661 27,8 N-terminus 2 764 11,5 Internal 5 315 22,2 ARNs no codificantes putativos 27 0,0 Transcritos sin ortólogo 45 834 65,2 Nuevos transcritos putativos 19 645 42,9 Desconocidos 26 189 57,1 Transcriptoma de referencia 23 671 - A partir del transcriptoma global (versión 2), 6 889 secuencias fueron anotadas con al menos un GO. El análisis de las funciones procedentes de la anotación con los términos GO (figura IV.36) muestra que el número más alto de transcritos anotados con procesos biológicos fue asociado con procesos celulares (26%) y metabólicos (24%) (figura IV.36-A). En lo que se refiere a los componentes celulares, las categorías más representadas fueron células (35%) y orgánulos delimitados con membranas (19%) (figura IV.36-B). Respecto a la función molecular, el mayor número de transcritos anotados fue en la categoría de actividad catalítica (52%) y fijación a receptores (36%) (figura IV.36-C). Por otro lado, 790 transcritos fueron asignados a un mapa de KEGG universal. Globalmente, las enzimas que intervienen en las rutas más importantes como aquellas del metabolismo del nitrógeno, de la fotosíntesis y del metabolismo de azúcar (figura IV.37) estaban bien representadas en los mapas de KEGG 148 RESULTADOS Y DISCUSIÓN Figura IV.36 : Distribución de los GO del transcriptoma de T. lutea según los procesos biológicos (A), componentes celulares (B) y función molecular (C) 149 RESULTADOS Y DISCUSIÓN Figura IV.37: Rutas metabólicas relacionadas con el metabolismo del nitrógeno (A), fijación de carbono (B) y el metabolismo de azúcar (C) según KEGG, en donde aparecen las enzimas de todos los organismos. Se resaltan en verde las enzimas de T. lutea 150 RESULTADOS Y DISCUSIÓN Marcadores SSR contenidos en el transcriptoma de T. lutea Por la importancia que tienen los marcadores moleculares en la exploración del genoma y en la mejora genética, en la tabla IV.21 se presentan algunas estadistas sobre los marcadores SSR encontrados en la versión 2 del transcriptoma global y en el de referencia. Se observa que en ambos casos las repeticiones trinucleotídicas fueron los SSR más abundantes seguidas con las repeticiones dinucleotídicas y superiores a tetra-nucleótidos («SSR, otros»). Tabla IV.21: Estadísticas de los SSR para el transcriptoma global y de referencia de T. lutea Transcriptoma global Transcriptoma de referencia Transcritos con algún SSR 9 500 4 650 SSR totales 12 277 5 937 SSR di-nucleotídicos 2 677 1 068 SSR trinucleotídicos 5 280 2 865 SSR tetranucleotídicos 1 763 8 40 SSR, otros 2 557 1 164 Motivo más abundante GA (644) GA (254) Comparación de T. lueta con otras especies. El transcriptoma de referencia de T. lutea v2 fue comparado con las proteínas de 6 microalgas cercanas (2 clorófitas, 2 glaucofitas, una rodófitas y una haptofita). El número de proteínas para cada microalga y el número de les genes homólogos entre T. lutea y estas microalgas de referencia se muestran en la figura IV.38. Como se esperaba, al ser Emiliana huxleyi una especie del mismo linaje y muy cercana de T. lutea, tiene la proporción más alta de genes homólogos con T. lutea (43%). Esta proporción se redujo al 18% con las microalgas de referencia de otros linajes. Por otra parte, el transcriptoma de referencia de T. lutea fue comparado con las proteínas RefSeq y Ensembl de E. huxleyi. (figura IV.39), mostrando que 10 216 transcritos fueron ortólogos a 8 170 entradas RefSeq y 8 109 entradas de Ensembl. Muchas haplofitas, incluida E. huxleyi, se sabe que tienen un ciclo de vida haplodiplonte en el cual ambas fases haploide y diploide son capaces de una reproducción asexual [180]. En T. lutea , el nivel de ploidía no está descrito [181]. En cultivo, los Isochrysidaceae tienen un único morfotipo sin calcificación y se cree que la calcificación de la fase diploide se implicó solo una vez en el origen de los cocolitóforos y aparentemente desapareció temprano en la historia evolutiva de los Isochrysidaceae [182]. Así, los genes encontrados en T. lutea son probablemente de la etapa haploide aunque no podemos ser seguros de ello. La realización de más estudios sobre el ciclo de vida de esta especie pueden ofrecer más conocimientos sobre la interacción de los alelos y la función de sus genes. 151 RESULTADOS Y DISCUSIÓN Figura IV.38: Análisis comparativo del transcriptoma de referencia de T. lutea con las proteínas de microalgas de referencia. Los tamaños de los círculos son relativos al número de transcritos o proteínas. Códigos de colores: verde: clorófitas, amarillo: rodófitas, naranja: glaucofitas, azul: haptofitas. Los números en los círculos periféricos forman la proporción (%) de los genes de los organismos de referencia homólogos con T. lutea. Los números en los círculos interiores forman la proporción (%) de los genes de T. lutea homólogos con los organismos de referencia Figura IV.39: Anotación del transcriptma de referencia con los ortólogos de Emiliania huxleyi utilizando los IDs de RefSeq y Ensembl