minimap2 vs BLASTN

Hola,
recientemente Heng Li publicó un trabajo (https://doi.org/10.1093/bioinformatics/bty191) describiendo un nuevo alineador genérico de nucleótidos que se llama minimap2, que podéis descargar en https://github.com/lh3/minimap2.

Figura tomada de https://doi.org/10.1093/bioinformatics/bty191

En el artículo se compara minimap2 en diferentes escenarios contra otros softwares alternativos, incluyendo su antecesor BWA mem y se destaca su velocidad y su versatilidad, ya que es capaz de alinear lecturas cortas, secuencias largas e incluso también puede alinear saltando intrones.

Yo lo que he hecho ha sido una prueba rápida para compararlo con BLASTN en el escenario habitual de GET_HOMOLOGUES-EST, donde se comparan por ejemplo todos los genes de una planta (Brachypodium distachyon) contra todos los genes de otra especie cercana (Oryza sativa). Esto es lo que he hecho:

# how many sequences
$ grep -c "^>" *fna
Bdistachyon.fna:36647
Osativa.fna:42189

# index and BLASTN search
$ ncbi-blast-2.6.0+/bin/makeblastdb -in Osativa.fna -dbtype nucl
$ ncbi-blast-2.6.0+/bin/blastn -query Bdistachyon.fna -db Osativa.fna \
   -out Bdistachyon.Osativa.blastn.tsv -dbsize 100000000 -evalue 1e-5 -outfmt 6
real	0m40.937s
user	0m40.280s
sys	0m0.636s

# index [assuming up 80% sequence identity] and minimap2 search
$ minimap2/minimap2 -x asm20 -d Oryza.mmi Osativa.fna
$ time minimap2/minimap2 Oryza.mmi Bdistachyon.fna > Bdistachyon.Osativa.minimap.paf
real	0m2.084s
user	0m3.360s
sys	0m0.300s

Ahora echemos un ojo a los alineamientos resultantes. Selecciono un par de secuencias, primero de BLASTN:

BdiBd21-3.2G0760100.1   LOC_Os01g70090.1        87.839  847     95      5       31      876     37      876     0.0     987
BdiBd21-3.2G0521100.1   LOC_Os01g37510.1        85.652  683     92      3       91      773     103     779     0.0     713

y ahora de minimap2, en formato PAF:

BdiBd21-3.2G0760100.1	876	155	776	+	LOC_Os01g70090.1	876	161	776	181	621	60	tp:A:P	cm:i:16	s
1:i:179	s2:i:0	dv:f:0.0980
BdiBd21-3.2G0521100.1	777	110	653	+	LOC_Os01g37510.1	783	122	659	87	543	60	tp:A:P	cm:i:10	s
1:i:85	s2:i:0	dv:f:0.1196

Al maneo para estos dos ejemplos podemos observar que:
i) el mejor hit de BLASTN y minimap coinciden
ii) los alineamiento de BLASTN son más largos


Hasta pronto, buenas vacaciones,
Bruno

Plant Genome Evolution 2017 (y III)

Hola,
termino esta serie con las del último día de este congreso. El próximo en dos años.

Edit09102017
Our poster "Pan-genomes: estimating the true genomic diversity of plant species" is available at https://digital.csic.es/handle/10261/156147

Pamela Soltis habla de los genomas de los helechos, que tienen muchos cromosomas y pueden haber sufrido varias rondas de poliploidización y por tanto experimentan silenciamiento genómico a gran escala. Encuentran que los individuos estudiados han perdido al menos un alelo, pero no son pérdidas fijadas en la población. Han estudiado la expresión específica de genes homeólogos en tetraploides y ven que aproximadamente la mitad muestran un sesgo hacia uno de los parentales.  Después habla de la aneuploidía compensada, donde los individuos de una población conservan el número de cromosomas pero no el patrón aditivo de los parentales, con trisomías y monosomías por ejemplo. Luego pasa a hablar de que no siempre coinciden en el tiempo la producción de duplicaciones genómicas (WGD) con la radiación de especies. Mientras las Asteraceas sí coinciden, en muchos otros ejemplos hay un retraso (http://www.sciencedirect.com/science/article/pii/S1369526612000465).

Jeffrey Chen continúa la sesión de poliploides, que normalmente son más grandes y vigorosos y experimentan un efecto de dominancia parental epigenética, que en muchos casos es heredable y reversible. Ellos trabajan con Arabidopsis suecica y hacen híbridos inter-específicos con A. thaliana y A. arenosa para estudiar cómo se modifican los ritmos circadianos y la fotosíntesis. También estudian el algodón tetraploide, porque hasta ahora no tienen el genoma de arenosa y en Gossypium hay más recursos genómicos. En estos materiales están estudiando como los subgenomas A y D tienen diferentes marcas en histonas que explican la dominancia de uno sobre el otro (https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1229-8). Entre los genes con silenciamiento específico hay genes de domesticación, floración y dormancia. Termina con un repaso de los efectos de la poliploidización en diferentes espeicies, donde hay cambios genéticos (Brassica), epigenético (A. thaliana, algodón) o ambos (trigo).

Michael Barker habla precisamente de poliploides del género Brassica. Empeiza recordando que un tercio de las plantas son poliploides, y que incluso muchas diploides son derivadas de eventos de poliploidización antiguos. Recuerda también que la poliploidización ha sido en general previa a la domesticación. Su pregunta es si los paleólogos, genes de poliploidizaciones ancestrales, están enriquecidos entre los genes de domesticación y sus datos en Brassica apuntan en esa dirección. En general sus datos indican que la edad del evento de duplicación más reciente de una especie determina su variabilidad genética.

Arp Schnittger explica como la poliploidización inducida con colchicina se usa en mejora y como, en general, para generar gametos aneuploides debe fallar algo en el control del huso en la meiosis (splindle checkpoint).  Encuentran que a diferencia de los animales, la formación del huso en plantas se aborta muchos antes en caso de estrés y eso podría explicar la facilidad de formación de poliploides (https://www.ncbi.nlm.nih.gov/pubmed/27816818).

Toni Gabaldón explica las herramientas desarrolladas en su grupo que en este caso usaron para el estudio filogenómico del olivo (Lamiales). Muestra datos de profundidad de sintenia con respecto al café. Muestra como definen eventos de duplicación sobre árboles de genes con el algoritmo de solapamiento de especies (https://academic.oup.com/bioinformatics/article/27/1/38/201693/Assigning-duplication-events-to-relative-temporal). Explica en detalle como genes de un híbrido aparecerán en un árbol como parálogos, pero se pueden distinguir porque la topología resultante tendrá menos copias ancestrales que las esperadas para una duplicación. Le preguntan si ha cruzado sus datos de árboles génicos con datos de sintenia, y dice que todavía no, porque su ensamblaje y el de las especies vecinas están fragmentados todavía.

Steve Maere también habla de duplicaciones genómicas y de cómo se retienen o pierden genes después. Destaca que los TFs parecen preferir multiplicarse por duplicaciones completas antes que por cambios a menor escala (http://www.pnas.org/content/102/15/5454.full). Relaciona estos patrones con el balance de dosis, que se preserva en el primer caso, pero no en el segundo. Se pregunta si balance de dosis realmente está detrás de estos patrones y quieren responder mirando familias de proteínas en diferentes especies. Encuentran que ninguna familia se retiene por completo, y que hay un rango amplio de pérdida, cuando mapean los genes sobre bloques sinténicos. Las familias más conservadas tienen anotaciones relacionadas con regulación y señalización, y destacan los TFs. Cuando miran a familias de TFs, los WRKY y F-box aparecen mucho más conservados que los MADS-box. Además, las familias más retenidas divergen menos en secuencia. Finalmente muestra que las familias con más retención, candidatas a ser más sensibles al balance de dosis, contienen genes que tienen fenotipos sensibles a dosis.

Pat Edger habla de dominancia en los genomas híbridos y alopoliploides estudiando dos especies parentales del género Mimulus y su híbrido, comparando después la expresión de genes de ambos subgenomas, observando una dominancia clara de M. luteus en la mayor parte de los casos. (http://www.plantcell.org/content/early/2017/08/16/tpc.17.00010). Luego introduce su trabajo sobre Fragaria vesca, todavía sin publicar.

Olivier Panaud empieza hablando sobre la variabilidad de tamaños genómicos de angiospermas, con una distribución con la moda en 600Mb pero llegando a superar los 5Gb (valores C). Sin embargo, los genes ocupan generalmente entre 100 y 200 Mb, lo que muestra que la diferencia es el espacio que ocupan los transposones (TEs). El modelo actual supone que el tamaño aumenta a medida que se acumulan TEs y se reduce por deleción. Por ejemplo, Oryza australiensis es 2x O. sativa a causa solamente de 3 familias de TEs. Han mirado un montón de especies del género con esta perspectiva (por ejemplo; http://www.nature.com/ng/journal/v46/n9/full/ng.3044.html). Han estimado que la vida media de transposones LTR-RT en arroz es mucho más corta que en animales, de aproximadamente 1.7Myr. Actualmente están analizando TEs en 3000 variedades de arroz y para hacerlo rápidamente indexan secuencias de TEs con la transformada BW y mapean los reads de esos genomas. Termina con una diapo donde hacen GWAS con el fenotipo de CNV de familias de TEs y encuentran un pico en el transposón en sí, lo que sugiere que son factores ambientales los que lo hacen saltar. En principio se puede GWAS usando TEs en vez de SNPs con cualquier carácter fenotípico.

Ezrha Mizrachi habla de la regulación de crecimiento secundario en leñosas, sobre todo Eucaliptus. Su trabajo reciente es sobre la identificación de genes y rutas metabólicas implicadas en la producción de xilema, sobre todo en cloroplastos y mitocondrias   (http://www.pnas.org/content/114/5/1195.short). Menciona que encuentran copias casi completas de más de 100 genes cp en el genoma nuclear, pero con sus datos de RNAseq descubren que no se expresan. Tiene un artículo reciente donde ponen a punto métodos de ensamblaje de transcriptomas de novo, pero no lo encuentro todavía, creo que en PLoS ONE. Solo menciona que no le gusta SOPAdenovo.

Jen Wisecaver centra su charla en el estudio de redes de coexpresión en metabolismo secundario, en rutas especializadas, combinando datos de expresión de diferentes condiciones y especies (http://www.plantcell.org/content/early/2017/04/13/tpc.17.00009). Dice que cada genoma tiene decenas de módulos coexpresados de 10 a 40 genes, y que capturan todos los que hay descritos que se agrupen en el mismo cromosoma.

Klass Vandepoele habla de inferencia de redes de TFs en plantas. Usa datos de Y1H,ChIPseq, DNAase-seq, coexpression data (GENIE3), PBMs, phylogenetic –profile sites (http://www.plantphysiol.org/content/early/2016/06/03/pp.16.00821). Los combinan y los comparan con el conjunto de datos de AtRegNet y calculan especificidad y sensibilidad (los mejores son Y1H y ChIPseq). Entonces entrenan un predictor supervisado que es mucho mejor que todos los métodos por separado. Usan enriquecimiento como estadístico. Usando este predictor anotan TFs y sus dianas con gran acierto en A. thaliana (https://www.biorxiv.org/content/early/2017/08/09/173559).

 

O Tzfadia habla de TranSeq, un protocolo barato de secuenciación de extremos 3’ de transcritos, que en sus manos usaron para mejor de manera significativa los modelos de genes en el genoma de referencia de tomate. Además, observan que este tipo de lecturas capturan con precisión los patrones de expresión de genes que obtuvieron por TruSeq. Finalmente cuenta que este tipo de secuencias capturan CNV con gran precisión. 

Jornada de agrigenómica 2017

Hola,
hoy he estado por la mañana en la 1ª jornada de agrigenómica organizada por QUAES en la UPV, en Valencia. Los 6 ponentes que hemos participado hemos hablado sobre las distintas aplicaciones de las tecnologías genómicas y de secuenciación en la mejora de cultivos. Éstas son mis notas.

Iraida Amaya nos habló de sus proyectos en torno a la mejora de caracteres en la fresa, usando marcadores SNP obtenidos con chips y por genotipado por secuenciación que clasifican por subgenoma. Estos proyectos tienen, además de genómica, metabolómica para el seguimiento de los 20 compuestos volátiles más importantes para el aroma de la fresa. Como hitos destaca la obtención de dos marcadores PCR para genotipar con un 91% de precisión fresas con aroma aceptable, así como su trabajo en curso para caracterizar por agroinfiltración FvMYB10 como regulador del color rojo de la fresa, y del color blanco de sus mutantes, que portan un cambio de marco de lectura. Para ello usaron la metodología QTLseq propuesta por Takagi, 2013.

Guillermo Marco, bioinformático de Sistemas Genómicos, nos dió un curso abreviado sobre diseño y análisis de experimentos de RNAseq, apoyándose en la revisión de Conesa, 2016. Como en un trabajo reciente nuestro (Cantalapiedra, 2017),  recomienda siempre validar las muestras por PCA para descartar réplicas y probar diferentes tuberías de análisis y estrategias de ensamblaje para elegir la más adecuada en cada caso. Su mensaje general es que el diseño del experimento RNAseq determina su éxito.

Patricia Pascual, directora de I+D de Agrícola El Bosque SL, nos presentó el diseño de su incipiente programa de mejora en mora y lo justificó como estrategia empresarial.

David Calvache, del INIA, nos habló sobre el uso de marcadores moleculares en los exámenes técnicos de nuevas variedades desde el punto de vista de la propiedad intelectual (título de obtención) y las licencias de producción y comercialización (registro de variedades). El exámen técnico permite comprobar la distinción, homogeneidad y estabilidad de una variedad candidata. Discute el uso de marcadores molecular para determinar caracteres y la idea de distancia mínima respecto a variedades existentes. Esta distancia puede ser fenotípica y/o molecular, en base a una batería de marcadores. Explica, si lo entiendo bien, que la distinción de variedades esencialmente derivadas la hace la justicia ordinaria apoyándose en distancias moleculares y pedigrís.

Esther Esteban, de la oficina española de variedades vegetales,nos recuerda los retos para la agricultura del futuro, en torno a la seguridad alimentaria, la sostenibilidad, la adaptación al cambio climático, las enfermedades emergentes, los estreses abióticos o las mejoras nutricionales. Explica las diferencias de normativa entre Europa y otros países sobre las nuevas tecnologías de mejora, incluyendo cisgénesis, transgénesis o las ZFN. Cita el trabajo de Lusser, 2011. Desde entonces han surgido nuevas tecnologías como la biología sintética y la edición mediante CRISPR/Cas9. Hace una reflexión sobre las políticas en torno a la biotecnología en Europa frente al resto del mundo.

Yo hablé sobre genómica de plantas y la aplicación de las herramientas de ultrasecuenciación en la mejora de cebada, apoyándonos en la colección nuclear de cebadas españolas, y usando ejemplos de trabajos recientes nuestros [ get_homologues-est , oídio , RNAseq en sequía , adaptación climática ] y figuras de la estupenda revisión de  Bevan, 2017.

Un saludo,
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

Seminario de introducción a la metagenómica

Os quiero anunciar que el 19 de Octubre impartiré el workshop titulado "Introduction to Bioinformatics applied to Metagenomics and Community Ecology" como parte de la conferencia ​Community Ecology for the 21st Century (Évora, Portugal).

Si estáis interesados, podéis contactar a los organizadores de la conferencia en el siguiente enlace, todavía hay plazas disponibles en el workshop.

Durante el curso presentaré la nueva herramienta AmpliTAXO para el análisis sencillo y online de datos de NGS de RNA ribosomal y otros marcadores.

El curso consistirá en 2 partes, la primera teórica donde se expondrán los retos de la metagenómica, las posibilidades de las nuevas técnicas de secuenciación y el funcionamiento de las herramientas de análisis más habituales (UPARSE, QIIME, MOTHUR). La segunda parte será práctica y consistirá en el análisis de datos metagenómicos reales obtenidos por NGS.

Podéis encontrar más información en inglés en mi nuevo blog y próximamamente en la página de la conferencia (pendiente de actualizar).

Os dejo un pequeño adelanto de los contenidos...

Cómo elegir software para simular reads

Hola,
si trabajas con cierta frecuencia con datos de ultrasecuenciación, es decir, con ficheros en formato FASTQ, a menudo te habrás encontrado con la situación de que te gustaría tener más datos para validar un programa publicado, o para diseñar tu propia tubería de análisis. Pues bien, una posibilidad real desde hace varios años es simular tus propios reads o lecturas a partir de secuencias genómicas del organismo en cuestión y parámetros como la tasa de error o la plataforma de secuenciación empleada. Como ocurre frecuentemente en Biología Computacional, hay una gran variedad de software publicado para esta tarea y elegir no es trivial. El diagrama a continuación es un árbol de decisión que os guiará para esta tarea:

Figura prestada de Escalona, Rocha & Posada (2016) doi:10.1038/nrg.2016.57 http://www.nature.com/nrg/journal/v17/n8/abs/nrg.2016.57.html

Buen finde,
Bruno

pseudoalineamientos y conteos de tránscritos

Hola,
durante el análisis de unos experimentos de RNAseq mi colega Carlos Cantalapiedra y yo nos desesperábamos al calcular los valores de expresión de cada transcrito en una serie de condiciones, dado que para cada una de ellas era necesario alinear las lecturas/reads originales (en formato FASTQ) contra los transcritos ensamblados. Al menos éste era el protocolo habitual antes de aparecer el algoritmo genérico de pseudoalineamiento, que se describe en la siguiente figura:

Pseudoalineamiento de reads a partir de su composición en k-meros (k=3 en el ejemplo) y un  vector de sufijos de un transcriptoma. Figura tomada de http://bioinformatics.oxfordjournals.org/content/32/12/i192.full.  

Este tipo de algoritmos, implementados en software como kalllisto o RapMap (éste último con licencia GPL), permiten estimar de manera muy eficiente y precisa con qué transcritos son compatibles las lecturas de un archivo FASTQ, es decir, a qué secuencias de un transcriptoma mapean, sin necesidad de alinear base a base los reads. Los alineamientos se pueden calcular, si fuera preciso, después del mapeo. 

Hasta la fecha se han propuesto al menos dos implementaciones del algoritmo genérico de pseudoalineamiento, que se diferencian fundamentalmente en las estructuras de datos utilizadas, como se explica en estos blogs (1 , 2) [en el segundo se hace un repaso a la evolución de este tipo de algoritmos y a su nomenclatura].

En los próximos párrafos muestro dos aplicaciones usando ejecutables de RapMap y como datos un fichero (tr.fna) con 67K transcritos de Arabidopsis thaliana, ensamblados de novo con trinity, y otro (1.f1) con 89M de lecturas Illumina SE en formato FASTQ.

1. Mapeo de reads (20 cores, fichero de salida SAM)

          operación   tiempo   comando      

BWAmem       índice    1m28s   bwa-0.7.12/bwa index tr.fna    

BWAmem        mapeo    9m57s   bwa-0.7.12/bwa mem -t 20 tr.fna 1.fq > 1.sam 

RapMap       índice      58s   RapMap-0.2.2/bin/rapmap quasiindex -t tr.fna -i ./ 

RapMap        mapeo    5m56s   RapMap-0.2.2/bin/rapmap quasimap -t 20 -i ./ -r 1.fq \

                                 -o 1.sam 

Por comparación con BWA queda claro que incluso generando alineamientos SAM estos algoritmos son mucho más rápidos. 


2. Conteo de transcritos (20 cores, implementación Sailfish)

          operación   tiempo   comando      

Sailfish     índice    1m06s   SailfishBeta-0.10.0/bin/sailfish index -t tr.fna -o qindex \

                                  -p 20   

Sailfish      mapeo    1m09s   SailfishBeta-0.10.0/bin/sailfish quant -i qindex/ -r 1.fq \

                                  -p 20 -o 1.bur-0.sf --auxDir tmpsf --dumpEq --libType SF
  

Como se muestra en la tabla, la estimación de valores de expresión de transcritos por pseudoalineamiento y conteo de reads es muy eficiente. En este ejemplo se genera un archivo de salida (quant.sf) con este contenido:

Name      Length EffectiveLength TPM NumReads
TR1|c0_g1_i1 517 316.623      5.3984 214.782
TR2|c0_g1_i1 390 191.501     3.36608     81
TR3|c0_g1_i1 699 498.611      10.502 658
...

Entre las opciones del program está la posibilidad de calcular los conteos con bootstraping, como hace kallisto, aunque en la versión que yo he probado es experimental.

En la siguiente entrada veremos más aplicaciones,
un saludo,
Bruno

mapeo fino de genes por NGS

Buenas,
esta semana copio aquí una reseña de un trabajo recientemente publicado de Carlos P Cantalapiedra, autor habitual de este blog y próximo doctor del grupo, donde se explica el proceso para localizar un loci responsable de una resistencia a infección por parte de hongos, combinando genética clásica y secuenciación de nueva generación: http://www.eead.csic.es/spreading/showspreading?Id=416

Pongo aquí una de las figuras del artículo:

Genotipo de varias líneas de cebada en torno al locus que confiere resistencia. En naranja, genotipos como los del parental resistente. En verde, genotipos como los del parental susceptible. La captura de exoma permite reducir la zona de búsqueda al punto donde se unen ambos genotipos (punto 211721 dentro del recuadro). Adaptada de https://dl.sciencesocieties.org/publications/tpg/first-look/pdf/plantgenome2015.10.0101.pdf.

La referencia del artículo completo, en inglés, es:

Cantalapiedra CP, Contreras-Moreira B, Silvar C, Perovic D, Ordon F, Gracia MP, Igartua E, Casas A. (2016) A cluster of NBS-LRR genes resides in a barley powdery mildew resistance QTL on 7HL. The Plant Genome. Early access. DOI: 10.3835/plantgenome2015.10.0101. URL.

Hasta luego,

Bruno

Calculando experimentos de secuenciación

Buenas,
hoy necesitábamos calcular cuántos individuos (de una especie monocotiledónea) podríamos secuenciar con cierta profundidad en un secuenciador Illumina, pensando en el HiSeq2500 en concreto. Al final decidimos apostar por una profundidad promedio de 80x, para is sobre seguro, usando parejas de lecturas de 2x125b. Buscando en Internet encontré rápidamente una calculadora del propio fabricante que igual algunos no conocéis y puede ayudar a hacer esto rápidamente.

Figura tomada de http://www.danielecook.com/calculate-depth-coverage-bam-file.

Vayamos con un ejemplo con la calculadora
[ http://support.illumina.com/downloads/sequencing_coverage_calculator.html ]:

0. tipo de secuenciación: DNA             [se puede elegir RNA también]
1. protocolo: whole-genome sequencing  [otras: Nextera, Truseq, custom]
2. tamaño del genoma: 320Mbp
3. profundidad deseada: 80x
4. % de duplicados: 2%                [valor por defecto]
5. instrumento: HiSeq 1500/2500

Volumen total de secuenciación requerido: 26,1Gb   [26.122.448.980b]

En mi ejemplo, usando el protocolo v4, esto equivale a 0.42 líneas o lanes, lo que significa que podría poner hasta 2 muestras por línea.

Hasta luego,
Bruno

PD Me dicen mis colegas Dave Des Marais y Pat Edger que la longitud de un genoma (de plantas en este caso)  puede estimarse aproximadamente a partir del contenido en DNA del núclo usando la fórmula long = 1C * 980.

Genotipado HLA y software disponible

En la presente entrada intentaré hacer una recopilación de software para el genotipado de HLA (MHC humano), en la sección de comentarios podéis añadir otras herramientas que intentaré incorporar al texto. Para escribir la entrada me resultó muy útil el siguiente post en Biostar forum y Omictools.

Antes de empezar es necesario mencionar que existe una base de datos muncial de alelos HLA denominada IMGT-HLA que recopila las miles de secuencias conocidas (y públicas) de genes y transcritos para esta familia. Todas las herramientas de genotipado emplearán las secuencias de esta base de datos para realizar sus predicciones.

Número de alelos registrados para cada tipo de HLA en la base de datos IMGT-HLA (Septiembre 2015).

Las herramientas de genotipado de HLA alinean (map) las reads de NGS a las secuencias de referencia de los principales loci HLA humanos presentes en la mencionada base de datos (clase I: HLA-A, HLA-B y HLA-C, clase II: HLA-DR y HLA-DQ).

Dependiendo del software, se puede procesar reads de secuenciación genómica, exómica o transcriptómica, aunque cuando se quieren analizar cientos/miles de individuos en un único experimento se suelen enriquecer por la técnica de secuenciación de amplicones (diseñando primers que amplifican las regiones menos conservadas, ver entrada anterior) o captura con sondas específicas para HLA.

El mapeo de reads a las secuencias de referencia puede realizarse directamente o tras realizar un ensamblaje de novo previo. Realizando un ensamblaje previo será más fácil encontrar alelos únicos de HLA puesto que los contigs resultantes darán menos mapeos ambiguos. Sin embargo el ensamblaje de esta familia de genes parálogos generará también un gran número de ensamblajes erróneos o quimeras (falsos contigs mezcla de dos secuencias análogas). A su vez el mapeo directo de reads puede generar ambigüedades puesto que numerosas reads alinearán con múltiples referencias a la vez.

Típica estrategia de genotipado por mapeado (alineamiento) de reads a sequencias de referencia. A la izquierda las reads son ensambladas de novo antes de ser alineadas. A la derecha las reads son directamente alineadas. Imagen modificada de Warren et al. (2012).

Listado de software para el genotipado de HLA

Sólo se listan herramientas libres para uso académico ordenadas por orden cronológico de la última versión del software:

• seq2HLA (Jun 2015): diseñado para procesar reads de RNA-Seq, mapea las mismas a las secuencias alélicas de referencia (IMGT-HLA) generando genotipos con una puntuación de probabilidad para los mismos y los niveles de expresión de los alelos predichos.
• HLAreporter (May 2015): primero filtra las reads que mapean a los diversos alelos de un único gen, las ensambla de novo y los contigs resultantes son de nuevo alineados a los alelos de referencia iniciales para asignar genotipos.
• HLAminer (Feb 2015): realiza un ensamblaje de novo de las reads (de casi cualquier procedencia) para después alinear los contigs resultantes contra los alelos de referencia.
• Optitype (Apr 2014): otro método que acepta diversos tipos de datos y también se basa en el mapeo a secuencias exónicas de referencia. Los resultados del mapeo son representados en forma matricial, las reads en filas y los alelos en columnas. En la matriz se identifican como máximo 2 alelos que explican el mayor número de reads mapeadas.
• PHLAT (Feb 2014): además de analizar datos genómicos, transcriptómicos y exómicos, ha sido también testado con datos de amplicones. Mapea reads a las secuencias de referencia seleccionando múltiples alelos candidatos y selecciona la pareja de alelos con la mayor probabilidad de acontecer juntos.
• HLAforest (Jan 2013): similar a seq2HLA, analiza reads de RNA-Seq, aunque puede ser usada con otro tipo de datos reduciendo su precisión.
• ATHLATES (Jun 2012): similar a HLAminer, filtra y ensambla las reads para después identificar exones de IMGT-HLA en los contigs ensamblados. Está diseñado para reads de sequenciación de exoma.
• GATK-HLA Caller (Dec 2011): similar a seq2HLA, alinea, filtra y calcula probabilidades para cada genotipo.

Por último explicaré el software diseñado en mi laboratorio...

• AmpliHLA (Sep 2015), no es el mejor, simplemente es diferente. Está únicamente enfocado al análisis de datos de secuenciación de amplicones usando primers que amplifiquen diferentes regiones de los genes HLA de interés y etiquetas de DNA que diferencien las muestras.

AmpliHLA requiere un pre-procesado online de los datos de NGS con la herramienta AmpliSAS. AmpliSAS clasifica las reads por muestra/individuo, corrige errores de secuenciación y filtra artefactos de secuenciación y PCR. AmpliSAS está diseñado para el genotipado de cualquier tipo de gen, especialmente si no tenemos alelos de referencia previos (como generalmente ocurre con muchos organismos cuyo genoma no ha sido secuenciado o regiones complejas del genoma como los genes que codifican las moléculas de MHC).

Un archivo Excel generado tras el análisis con AmpliSAS ha de ser introducido en el formulario de AmpliHLA y el programa automáticamente unificará marcadores (diversas regiones amplificadas de un mismo gen) y buscará sus secuencias en la base de datos humana para genotipar con la máxima precisión posible cada individuo. El principal inconveniente es el requerimiento de múltiples PCRs y diversos marcadores por gen para conseguir un genotipado de calidad. La principal ventaja es la obtención de un genotipado de-novo que permite descubrir alelos incluso si no están presentes en la base de datos humana.

Esquema de funcionamiento de la herramienta de genotipado de novo mediante secuenciación de amplicones: AmpliSAS. Primero las reads son separadas por muestras y marcadores. Después se realiza un clustering de los errores de secuenciación con sus secuencias de origen. Por último se filtran las reads minoritarias y se asignan alelos por cada muestra y marcador.

Secuenciación de amplicones y genotipado de alto rendimiento

Secuenciación de amplicones (SA)  es una traducción aproximada al español de la técnica de "Amplicon sequencing" que junto con las tecnologías de secuenciación masiva (del inglés new generation sequencing, NGS) permite genotipar cientos/miles de individuos en un único experimento.

La secuenciación de amplicones (SA) consiste en secuenciar los productos de múltiples PCRs. Un amplicón se define como el conjunto de secuencias obtenidas de cada PCR individual.

Antiguamente se realizaban PCRs individuales y se secuenciaban uno a uno los productos. Con las nuevas técnicas de NGS, podemos incluir etiquetas de DNA (por ej. una secuencia única de 6 nucleótidos) diferentes para cientos de muestras o individuos y clasificar más tarde las secuencias o reads resultantes de una única secuenciación (Binladen et al. 2007; Meyer et al. 2007).

Esquema de etiquetado y amplificación para la secuenciación de amplicones.

Mediante esta técnica podremos genotipar individuos y distinguir los diferentes alelos (con la secuenciación tradicional a veces es complicado separar alelos de un mismo gen). El principal problema de las técnicas de NGS es su alta tasa de error, que a su vez puede ser compensada incrementando la profundidad de secuenciación (el número de reads). Otros problemas pueden ser los errores de la polimerasa o la generación de quimeras (una secuencia mezcla de otras).

En el siguiente enlace podemos ver un vídeo explicativo del proceso de secuenciación de amplicones:
http://www.jove.com/video/51709/la-secuenciacin-de-prxima-generacin-de-16s-arn-ribosomal-genes?language=Spanish

Básicamente existen 4 etapas en el análisis por AS con NGS:

1. Diseño experimental de los primers usados para amplificar los genes de interés (marcadores) y las etiquetas a usar para distinguir los diferentes individuos o muestras.
2. Amplificación por PCR de los marcadores en el laboratorio, generalmente se realiza una PCR por cada muestra.
3. Secuenciación de los productos de amplificación. Las tecnologías de NGS más usadas para SA son: Illumina, 454 e Ion Torrent.
4. Análisis bioinformático de los datos de secuenciación. El análisis incluye separación de las reads en amplicones, corrección de errores de secuenciación, filtrado de reads minoritarias/contaminantes y generación de genotipados.

Etapas de la técnica de secuenciación de amplicones mediante NGS.

Aplicaciones 

SA es utilizado para realizar clasificaciones taxonómicas usando genes como: cytochrome c oxidase subunit 1 (CO1), genes rRNA (16S/18S/28S), genes específicos de plantas (rbcL, matK, and trnH-psbA) y espaciadores internos nucleares (ITSs) (Kress et al. 2014; Joly et al. 2014). Los genes anteriores se distinguen por una tasa de mutación suficientemente rápida como para distinguir especies cercanas y a la vez suficientemente estables como para distinguir congéneres.

Lista de genes habitualmente usados como marcadores taxonómicos (fuente: Kress et al. 2014)

Un experimento pionero de SA fue la determinación de la diversidad microbiana en aguas marinas profundas (Sogin et al. 2006), usando primers flanquenado la región V6 hipervariable de la subunidad 16S rRNA bacteriana. Dicho estudio descubrió miles de poblaciones minoritarias de organismos no conocidos con anterioridad.

Otro gran campo de aplicación es el genotipado de familias de genes de alta complejidad, como el complejo mayor de histocompatibilidad, que poseen múltiples loci y diferente número de copias entre individuos, incluso de la misma especie (Babik et al. 2010; Lighten et al. 2014). El complejo mayor de histocompatibilidad (MHC) de clase I y II codifica receptores celulares que presentan antígenos a las células del sistema inmune y son los genes más polimórficos conocidos en vertebrados . El MHC humano también se conoce como HLA (Human Leukocyte Antigen) y juega un papel clave en la compatibilidad en el transplante de órganos. Los loci del MHC son tan polimórficos que no hay dos individuos en una población no endogámica que posean el mismo conjunto de alelos (excepto gemelos).

Estadísticas del número de alelos conocidos para la familia de genes del HLA (fuente: base de datos IMGT-HLA)

Hasta hace poco era necesario clonar y secuenciar uno por uno los diferentes alelos de este tipo de genes para conseguir una secuencia fiable. Actualmente tan tediosa tarea puede ser simplificada mediante un único experimento de NGS que incluya múltiples individuos y múliples genes. El secuenciador de nueva generación (ej.: Illumina, 454 o Ion Torrent) leerá las sequencias individuales de cada uno de los alelos. Actualmente existen incluso kits comerciales para simplificar el proceso: Illumina TruSeq Custom Amplicon, Roche 454 Fluidigm Access Array or Life Technologies Ion Torrent Ion AmpliSeq.