Dos versiones del mismo genoma: el caso de la cebada

Hola,

una de las plantas con las que hemos trabajado más en los últimos años es la cebada (Hordeum vulgare), un cultivo de enorme importancia en el secano de la península ibérica (ver por ejemplo https://www.eead.csic.es/barley). Puedes leer aquí sobre otros genomas de plantas y para qué los usamos.

La cebada es una especie diploide con un genoma haploide de 5.3 Gbp repartidos en 7 cromosomas, con un %GC=44.5. Es un genoma de gran complejidad que ha requerido un esfuerzo de al menos un década por parte de un consorcio internacional. Las diferentes etapas de la secuenciación y ensamblaje del cultivar Morex se pueden recorrer en:

• https://doi.org/10.1105/tpc.110.082537 (2011)
  
• https://www.nature.com/articles/nature11543 (2012)
  
• https://www.nature.com/articles/nature22043 (2017, IBSC_v2)
• https://doi.org/10.1186/s13059-019-1899-5 (2019)
  
• https://doi.org/10.1093/plcell/koab077 (2021, MorexV3)

Mi intención es mostrar cuánto pueden cambiar dos versiones del mismo genoma a medida que se incorporan nuevas tecnologías, en este caso lecturas largas PacBio (circular consensus sequencing, CCS). 

En este ejemplo comparamos la versión más reciente (MorexV3) respecto a una anterior (IBSC_v2). Si consultamos el European Nucleotide Archive, podemos ver sus estadísticas generales (verás que he puesto en negrita los cambios más notables):
 

IBSC_v2 (GCA_901482405.1)

Total Length:           4,833,791,107
Ungapped Length:        4,446,895,020
Chromosomes & Plasmids:             0
Spanned Gaps:               1,030,196
Scaffolds:                          8
Scaffold N50:             657,224,000
Contigs:                    1,030,204
Contig N50:                    19,388

MorexV3 (GCA_904849725.1)

Total Length:           4,225,577,519
Ungapped Length:        4,224,251,725
Chromosomes & Plasmids:             7
Spanned Gaps:                     162
Scaffolds:                        290
Scaffold N50:             610,333,535
Contigs:                          452  
Contig N50:                69,630,691 

Se puede ver que el nuevo ensamblaje es más pequeño, contiene muchos menos huecos (gaps) y tiene 7 cromosomas construidos a partir de un conjunto de apenas 452 contigs enormes.

En cuanto a la anotación, tomando datos de Ensembl Plants, si la versión IBSC_v2 tenía 39,841 genes codificantes de proteínas de "alta confianza, HC", la más reciente tiene 35,825. Finalmente, si mapeamos con bwa SNPs de IBSC_v2 en ventanas de 150 nucleótidos contra el nuevo genoma, de un total de 16.5M logramos 13.4M mapeos únicos y perdemos 2.7M SNPs por el camino (un 16%).

En la versión de 2021 Mascher et al hacen incapié en la significante ganancia en la calidad de las anotaciones de elementos repetidos, como se ve en la figura para los retrotransposones del tipo BARE1: 

 

                               Fuente: https://doi.org/10.1093/plcell/koab077

En resumen, el nuevo genoma es sustancialmente diferente al anterior y visto lo visto, seguramente no sea el último, aunque todos ellos han sido muy útiles sin duda,

Bruno

Apuntes sobre ensamblaje de genomas de plantas

Buenas, ayer asistimos Ernesto Igartua y yo al 6th CNAG Symposium on Genome Research: Agrigenomics, organizado por el Centro Nacional de Análisis Genómico en Barcelona, donde a menudo contratamos servicios de secuenciación.

Allí presentamos nuestro trabajo con cebada, junto a otros colegas que trabajan en ganadería, piscicultura y agricultura y utilizan herramientas de la genómica contemporánea.

Como curiosidades me apunté que André Eggen, de Illumina, mencionó que comparando razas bovinas habían imputado SNPs mezclando genotipos de baja densidad (chips de ~10K SNPs), con genomas completos, alcanzando millones de SNPs. Por cierto, habían usado el software propietario DeNovoMAGIC para ensamblar genomas bovinos.

Otra cosa fue que los peces que estudian Franscesc Piferrer y su grupo tienen un mecanismo de metilación en función de la temperatura para controlar la producción de hormonas sexuales, algo que me recordó mucho a la memoria de vernalización en las plantas.

Pero además de estas charlas, y de visitar las salas de secuenciación y de servidores del CNAG, tuvimos dos sesiones casi seguidas donde repasamos los últimos métodos de ensamblaje y validación de genomas de plantas de la mano de Tyler Alioto y Gareth Linsmith. Éstas son mis notas:

Detección de contaminantes en las lecturas/reads: 
kraken : https://ccb.jhu.edu/software/kraken

Ensamblajes híbridos y diploides, combinando lecturas cortas y largas y estrategias más complejas para genomas de individuos heterocigotos.

• reads cortos, generalmente Illumina, de entre 100 y 300b, para alcanzar profunidades de al menos 30X en cada tipo de librería: 
• paired-end (PE) con insertos de por ejemplo 400 y 730pb 
• mate-pair (MP) con insertos de 4 y 8Kb para superar la longitud de la mayoría de secuencias repetidas

• reads largos, generalmente PacBio o de Oxford Nanopore. EN CNAG usan secuenciadores minIon para producir lecturas de 11.5Kb de media, alcanzando longitudes máximas > 100kb. Gareth comentó que en manzano necesitaron 60x, y eso que era material doble haploide. Este tipo de reads requieren consensos calculados con software como Sparc, Racon o Nanopolish.

En cuanto a ensambladores, Tyler destacó DISCOVAR de novo y Platanus, más adecuado para individuos con moderadas tasas de sitios heterocigotos. Pero advirtió del efecto negativo que tiene la heterocigosis sobre N50. En cambio, Gareth mencionó que primero ensambla las lecturas cortas con SOAPdenovo sin resolver las burbujas de Bruijn para luego luego combinar los reads largos con DBG2OLC y CANU.

Estrategias complementarias de ensamblaje
Datos de RNAseq para scaffolding con AGOUTI y Rascaf.

Pools de fósmidos como los empleados en el genoma de la ostra.
Mapas ópticos con enzimas nickasas que cortan cada 10Kb, con Bionano.
Dovetail genomics, aproximación basada en Hi-C.

Herramientas para corregir y finalizar genomas
PILON : https://github.com/broadinstitute/pilon/wiki
BESST : https://github.com/ksahlin/BESST

Estrategias para evaluar y validar genomas
Aparte del criterio clásico de sintenia respecto a especies cercanas, ambos mencionaron los problemas de evaluar un ensamblaje solamente por su N50 sin mirar por ejemplo los genes core anotados, por ejemplo con BUSCO, el sucesor de CEGMA. Gareth mencionó ALE para calcular la verosimilitud de un ensamblaje dadas las librerías de secuencias y KAT para comparar los k-meros originales de los reads con los del ensamblaje, que deberían coincidir, o para determinar la fracción de sitios heterocigotos:

Frecuencias de k-meros de los genotipos B73 y Mo17 de maíz, tomada de http://www.nature.com/articles/srep42444.

Casi se me olvida mencionar la comparación entre el mapa físico y el genético como criterio de calidad, muy útil en el genoma de manzano o en el de la cebada:

Comparación entre las posiciones de marcadores en una población de mapeo en cebada y sus posiciones en los mapas físico IBSC y POPSEQ de cebada, tomada de http://link.springer.com/article/10.1007%2Fs11032-015-0253-1.

Hasta  pronto,
Bruno

One-liner para calcular el N50 de un ensamblaje

Hola,
el estadístico N50 se usa a menudo para resumir a grosso modo un ensamblaje de secuencias, que generalmente es un fichero FASTA con una serie de contigs.
Modificando la definición de la wikipedia, N50 es la longitud de contigs tal que usando contigs de igual o mayor tamaño recuperamos la mitad de las bases de ese ensamblaje.

Figura tomada de http://bmpvieira.github.io/assembly14.

Hoy solamente quería compartir un comando one-liner de Perl que nos permite calcularlo:

$ perl -lne 'if(/^>(\S+)/){$h=$1} else {$l=length($_);$TL+=$l;$L{$h}+=$l} \

END{ foreach $s (sort {$L{$b}<=>$L{$a}} keys(%L)){ $t+=$L{$s}; \

if($t>$TL/2){ print "N50=$L{$s}";exit }}}' assembly.fasta

Cuidado al copiarlo al terminal desde el blog, mejor hacerlo línea a línea,
hasta pronto, Bruno

Apuntes de genómica microbiana

Hola,
esta semana hemos participado en un taller de 6 días organizado por Mirko Rossi (Universidad de Helsinki, Finlandia) donde hemos repasado los principales aspectos de la genómica comparativa de bacterias. Todo esto desde un punto de vista muy aplicado, dado que la mayor parte de los alumnos eran investigadores de la esfera veterinaria. Nuestra aportación fue una sesión de tres horas y media dedicada a la determinación del pangenoma de un grupo de bacterias de interés, con la ayuda de nuestro software GET_HOMOLOGUES.

Estimación de las dimensiones del core-genoma (izq) y pan-genoma (der)
de una colección de 12 plásmidos plncA/C, genomas de ejemplo en el capítulo
"Robust identification of orthologues and paralogues for microbial pan-genomics
using GET_HOMOLOGUES: a case study of pIncA/C plasmids", parte del libro
http://www.springer.com/biomed/medical+microbiology/book/978-1-4939-1719-8 .

A parte de esto, en esta entrada quería sobre todo contar algunas de las cosas que otros profesores a los que pude escuchar contaron, por si os pueden ayudar en vuestros proyectos:

• Con Andrey Prjibelski, del Algorithmic Biology Lab de San Petersburgo, aprendimos a ensamblar genomas bacterianos con el software SPAdes, y luego a comparar y a evaluar distintos ensamblajes con ayuda de QUAST. Nos recomendó http://nucleotid.es para comparar de una manera objetiva las virtudes y defectos de distintos ensambladores disponibles ahora mismo, aunque comentó que dada la variedad de escenarios de los usuarios y de los tipos de lecturas/reads existentes no hay ahora mismo ninguno claramente superior a los demás. Sí remarcó que algunos están mal documentados y son por tanto más difíciles de utilizar. Me quedé un poco con la idea de que Velvet y SPAdes son los mejores para genomas bacterianos. Del primero menciona que es muy estable, de calidad moderada, pero muy limitado con lecturas del tipo mate pair. El segundo es el que desarrollan ellos y fue el que usamos en la sesión. SPAdes trabaja siempre sobre una base de lecturas Illumina/IonTorrent, que por defecto se corrigen antes de ser utilizadas, pero puede hacer ensamblajes híbridos con secuencias de otras plataformas, incluso Sanger, útiles para resolver secuencias repetidas, como mate pairs convencionales o Nextera (hqmp, de alta calidad). El programa tiene además un módulo (--careful) para corregir por mapeo errores (indels,mismatches) en los contigs generados al ensamblar. Su mayor ventaja respecto a Velvet es probablemente que ensambla iterativamente con K-meros diferentes para superar las limitaciones de K pequeños (ciclos) con los grandes y viceversa (falta de solapamientos cortos, ver este blog y este otro para una comparación más amplia).  
• Bastien Chevreux, de DSM Nutritional Products AG, filosofó acerca de distintos problemas y estrategias para secuenciar y ensamblar genomas bacterianos con su software MIRA. La diferencia fundamental de este ensamblador respecto a los anteriores es que no se basa en grafos de De Bruijn, como Velvet o SPAdes, si no en la acumulación y alineamiento de lecturas solapantes, razón por la que es menos escalable y por tanto más lento que ellos. Este software a mi me parece sin duda más complejo de usar que SPAdes, pero tiene a su favor una lista de correo muy activa y documentación extensa. MIRA se ejecuta simplemente sobre un fichero manifiesto, editado por el usuario, que puede tener una sintaxis bastante compleja, pero que para un mapeo es algo sencillo, como por ejemplo:

  project= nombre_de_proyecto
  job= mapping,genome,accurate
  
  #input data comes next
  readgroup
  is_reference
  data= /path/to/file.gbk
  strain= Ecoli_strainA
  
  readgroup
  technology= solexa
  autopairing # busca PEs por ti
  data= /path/to/reads?.fastq
  strain= Ecoli_strainK12

• Rauni Kivistö, de la Universidad de Helsinki, nosexplicó el uso de los programas ACT y BRIG para visualizar datos de comparación de genomas basados en BLAST.
  
•  Veronika Vonstein y Ross  Overbeek, de la Fellowship for the Integration of Genomes, presentaron su plataforma para la anotación automática de genomas RAST que se construye sobre SEED, una colección de componentes, algo parecido a operones o trozos de rutas metabólicas, que llevan curando a mano desde hace dos décadas y que les permiten proyectar en tiempo real, según se van secuenciando nuevos genomas, anotaciones de alta calidad, que tienen en cuenta el contexto genómico y no se basan en BLAST, sino en eficientes búsquedas de K-meros específicos de familias de proteínas prealineadas que llaman FIGFams. Me sorprendió que un trabajo de esta envergadura no tuviera más repercusión, de hecho yo no lo conocía, pero se explica porque durante mucho tiempo trabajaron en el sector privado. Ahora mismo tienen financiación pública y se apoyan en hardware del Argonne National Laboratory.
  

Creo que esto es todo, un saludo y hasta la próxima,
Bruno

grafos de De Bruijn

Hola,
un artículo educativo publicado recientemente en Nature Biotechnology me ha vuelto a recordar los grafos de De Bruijn, que ya había mencionado de pasada en una entrada anterior sobre vectores de sufijos. Hoy los veremos un poco más de cerca, por su importancia para la reconstrucción de secuencias genómicas a partir de fragmentos más pequeños que han sido previamente secuenciados.

Los modernos ensambladores de secuencias emplean grafos de De Bruijn para guardar en memoria los prefijos (nodos) y sufijos (aristas) de las lecturas o reads obtenidas en experimentos de secuenciación. Este tipo de representación permite reducir en parte la redundancia natural de este tipo de datos a la vez que facilita su ensamblaje posterior. La figura  compara los métodos de ensamblaje tradicionales (a) con los actuales (d). Los tradicionales se usan para pegar entre si secuencias obtenidas por el método de Sanger, que suelen ser largas, de buena calidad y en números pequeños. Los métodos que acompañan a los secuenciadores actuales se adaptan al nuevo escenario de secuencias cortas y en números varios órdenes de magnitud más elevados.

fuente: http://www.nature.com/nbt/journal/v29/n11/full/nbt.2023.html

Como se muestra en el panel d de la figura, un grafo de De Bruijn facilita la reconstrucción del genoma de partida, que se obtiene buscando un ciclo euleriano que recorra todas las aristas (sufijos) una sola vez. Esta estrategia es computacionalmente mucho más barata que la alternativa de buscar ciclos hamiltonianos, que en la práctica es imposible al ser un problema de tipo NP completo. Sin embargo, como pasa a menudo, la realidad se resiste a ser capturada por una sola teoría, e incluso los grafos de De Bruijn son todavía incapaces de resolver el problema de las secuencias repetidas, que aparecen en un genoma o región genómica no una sinó muchas veces.

En el Russell's blog encontraréis código fuente en Python para aprender cómo se construyen estos grafos. Espero haberos despertado la curiosidad, un saludo,
Bruno