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.

Es primavera: lidiando con Flower

Hola,

vamos a hablar un poco sobre el trabajo con secuencias obtenidas de experimentos de NGS (Next Generation Sequencing; para despistados) con el secuenciador 454. En concreto, nos vamos a centrar en una herramienta que puede servir para tener un primer vistazo de los datos obtenidos, antes de realizar el ensamblaje, mapeo, ...

Flower (Bioinformatics, 2011) es un programa desarrollado por Ketil Malde, también implicado en el desarrollo de FlowSim. Está escrito con Haskell, el lenguaje de programación funcional, y anteriormente se distribuía el paquete como tal (última versión flower_v0.7) pero hoy día viene en la librería biosff (v0.2). Se trata de una alternativa GPL a las SFF Tools que distribuye Roche con el paquete Data Analysis de su software propietario.

Para instalarlo lo más sencillo es utilizar cabal, gestor de paquetes de Haskell, y bajar además ghc, entorno para desarrollo y compilación en ese lenguaje.

sudo apt-get install ghc cabal-install
cabal install biosff

Vamos a desarrollar algunos ejemplos con los datos de la entrada SRR000001 del NCBI SRA, run que se llevó a cabo en un equipo Roche 454 GS FLX en 2007. Quizás la opción más simple es la que nos ofrece un vistazo rápido del experimento, dándonos información de índice generado (creo que es un índice a la manera de un suffix array), el número de reads, número de flows y la key usada para el run: secuencia de bases que se añade a toda la muestra durante la preparación de las librerías.

flower -i SRR000001.sff
Index: (782054672,9419708)
Num_reads: 470985
Num_flows: 400
Key: TCAG

En éste caso vemos que hay un index (se regeneró mediante sfffile -o, ya que las secuencias de SRA no suelen llevar el índice). Con aquel run obtuvieron 470K reads en 400 flows (100 ciclos: 100 nts/read de media esperada). Por último nos indica que la key es TCAG (aunque al parecer siempre ha sido la misma).

Flower (algo así como FLOW ExtracteR), puede generar FASTA y FASTQ con scores en Phred+33 o basados en la codificación de Illumina. La carencia de un comando para obtener un FASTQ es una de las cosas que más parece echar de menos la comunidad que utiliza SFF Tools. Para obtener el FASTQ en Phred+33 con Flower:

 flower -Q SRR000001.sff > SRR000001.fastq

Con la opción -s se obtienen datos sobre cada uno de los reads del experimento. La salida del programa da los campos separados convenientemente por tabuladores.

 flower -s SRR000001.sff > SRR000001.reads
# name........     date......     time....     reg     trim_l     trim_r     x_loc     y_loc     len     K2     trimK2     ncount     avgQ     travgQ
EM7LVYS01C1LWG     2 7-04-10      15:13:00     01      1          235        1131      1422      255     0.74   0.76       0          26.38    26.33

K2 es una métrica de calidad "K-square", que es la suma de cuadrados de las distancias desde el entero más cercano para cada flow value. Como el sistema de scoring de 454 se basa en estimar la longitud de un homopolímero en base a la señal registrada en un flow, un entero supondría la definición de una longitud de homopolímero con exactitud. ncount es el número de flow cycles sin valor determinado o basecalls ambiguos, que también llevan asociado un error si el flow value no es exactamente 0. avgQ es la calidad promedio del read. Todos los campos que comienzan con tr son equivalentes a los anteriores, pero tras realizar el trimming indicado en el SFF.

Además, para el análisis del experimento también se pueden generar datos para representación gráfica de los flow values. Con la opción -h se obtiene la distribución de flow values por nucleótido.

 flower -h SRR000001.sff > SRR000001histogram.txt
Score     A        C       G        T
0.00     361230    586168  620918   409514
0.01     136705    284992  305021   168203
0.02     182110    407799  437550   228386
0.03     239742    571509  614918   306284
0.04     311296    779612  839794   401568
0.05     397364   1029955 1110715   513386
0.06     503853   1318643 1412014   648700
....
así hasta 99.99

Con la opción -H se obtiene la distribución de flow values para los flow cycle. Si representamos el resultado de esta salida, tenemos un gráfico como el siguiente, donde se muestra claramente que el sistema está optimizado para obtener la menor tasa de error posible. Como ya se ha comentado, la tasa de error está asociada a los puntos entre cada dos enteros. Para poder dar el mayor quality score posible, la mayor confianza posible de que la longitud del homopolímero es la predicha por el basecall, se debe maximizar el número de flow values lo más cercanos que sea posible a un entero.

Imagen tomada de Bioinformatics (2010) 26 (18): i420-i425

En un principio, Flower incluia una herramienta para seleccionar reads del SFF (flowselect). Sin embargo, en la versión de biosff esta utilidad ha desaparecido. Podría venir a suplir las opciones -i y -e de sfffile, que generalmente se usan ya durante el proceso de análisis, para generar subconjuntos de los datos originales una vez se tienen datos de mapeo o ensamblaje.

Como conclusión, Flower es más que una alternativa a SFF Tools para el pre-análisis de las secuencias de 454, en vista de la buena representación de los datos, el buen desempeño que tiene y las opciones extra que nos aporta. Por otro lado, sigue sin cubrir algunas de las utilidades que presenta SFF Tools, sfffile en especial, pero es código libre así que todos estamos invitados a mejorarlo.

saludos!

Carlos