Cómo funciona

Alcance y compatibilidad

De qué se responsabiliza get_MNV al anotar indels y MNV, qué deja en manos de tu llamador de variantes, y dónde están sus límites.

get_MNV es un anotador de variantes y un generador de resúmenes de haplotipos para llamadas de variantes ya existentes. No reemplaza a un llamador de variantes: el VCF o el TSV de iVar de entrada sigue siendo el responsable de decidir qué alelos existen. Cuando se proporciona un BAM, get_MNV usa la evidencia de las lecturas para contar el soporte y para emitir haplotipos combinados exactos de indel/SNV/MNV observados en las mismas lecturas.

Para cómo se lee un indel a partir de los alineamientos, qué cuenta cada número de la salida, y cómo se descubren y reportan los haplotipos locales, mira Indels y haplotipos locales. Para si las variantes de un haplotipo van juntas de verdad, mira Ligamiento.

Compatibilidad con llamadores#

Distintos llamadores pueden representar el mismo haplotipo local de maneras diferentes. Por ejemplo, FreeBayes se basa en haplotipos y puede emitir SNPs, indels, MNPs y eventos complejos como llamadas de haplotipo en bruto. Su propia documentación recomienda descomponer esas llamadas tras el llamado si se necesitan registros primitivos de SNP/indel.

bcftools norm se ocupa de la otra cara del problema: puede alinear a la izquierda y normalizar indels con una referencia FASTA, separar sitios multialélicos y atomizar variantes complejas o MNV en registros más simples. Ese paso de normalización resulta útil antes de comparar conjuntos de llamadas procedentes de distintos llamadores.

GATK HaplotypeCaller es un llamador, no un anotador. Recurre al ensamblaje local y a parámetros de genotipado como la ploidía de la muestra, el número máximo de haplotipos y la distancia de fusión de MNP. Por ahora, get_MNV no realiza este descubrimiento de haplotipos de novo ni el modelado de verosimilitud de genotipos.

Alineamiento a la izquierda#

Una inserción o una deleción dentro de un homopolímero o de una repetición en tándem puede escribirse en varias posiciones que describen el mismo cambio: borrar un GCT de GCTGCTGCT quita la primera copia, la segunda o la tercera. Los alineadores ponen el hueco en la de más a la izquierda, y get_MNV cuenta el soporte de un indel a partir de las lecturas cuyo CIGAR lleva ese hueco exacto, así que un alelo escrito más a la derecha no encuentra ninguna lectura que lo respalde. Un VCF escrito a partir de una descripción HGVS desplazada hacia 3', o por un llamador que coloca los indels a su manera, hace justo eso.

--left-align lleva cada indel a su posición equivalente más a la izquierda sobre la referencia de --fasta antes de que nada más lo lea, que es la normalización que hace bcftools norm -f ref.fa:

• Una inserción o una deleción se desplaza a la izquierda mientras el resultado describa el mismo cambio, y se escribe con la base anterior como ancla. Si llega a la primera base del contig, donde no hay base anterior, se ancla en la base siguiente, como pide la especificación VCF.
• Se recortan las bases que REF y ALT comparten en cualquiera de los dos extremos, como las recorta --normalize-alleles, así que --left-align incluye ese flag.
• Las sustituciones (SNV y MNV) y los registros complejos (una deleción y una inserción en un mismo registro, o un indel y una sustitución) se recortan pero nunca se mueven: una vez recortados, sus dos alelos acaban en bases distintas, así que ninguna otra posición los describe. Un registro complejo no se parte en registros más simples, que es lo que hace bcftools norm -a.
• Se quedan como están escritos: los alelos simbólicos que get_MNV no escribe de forma explícita (<INV>, <INS>, eventos de más de 10.000 bases), los breakends, *, los alelos con un código de ambigüedad como R, que bcftools también rechaza, y un REF que no coincide con la referencia, que detiene la ejecución igual que detiene a bcftools.
• Un <DEL> o un <DUP> que get_MNV escribe a partir de la referencia es un indel explícito desde ese momento y se alinea a la izquierda como cualquier otro. bcftools mueve también los simbólicos y reescribe su END, así que un <DEL> de más de 10.000 bases, que get_MNV deja como está escrito, puede quedar aquí más a la derecha que en la salida de bcftools.
• La entrada TSV de iVar también se alinea a la izquierda. iVar lee sus indels del alineamiento, así que normalmente ya están donde los tienen las lecturas.
• A diferencia de bcftools, que lee cada registro por separado, un indel se detiene en una base que cambia otro registro de la muestra en las mismas moléculas. iVar y los llamadores de haplotipos escriben una deleción detrás de una sustitución anclada en la base sustituida (108 A>T junto a 108 AA>A), y en las moléculas que llevan las dos la repetición acaba en la sustitución, como acaba el recorrido de un alineador. Un registro en otras moléculas no lo detiene: otro alelo indel del mismo sitio, uno que el llamador fasó aparte, o uno con menos de la mitad de la frecuencia del indel.
• El END de un registro movido, que se conserva con --keep-original-info, se reescribe para el registro tal como queda escrito. Los registros que se convierten en un mismo cambio al moverse se conservan, como los conserva bcftools sin -d, y un aviso los nombra, porque cada uno se cuenta con todas las lecturas de ese cambio.

El cambio de proteína y las descripciones HGVS no dependen de dónde se escribió el registro, porque HGVS coloca un indel en su posición más 3' diga lo que diga la entrada. Lo que se mueve es Positions, Reference Bases, Base Changes, Event Components, el registro del VCF de salida y, con ellos, el soporte de lecturas. El alineamiento a la izquierda ayuda a ese soporte solo cuando el alineador colocó los huecos a la izquierda, que es lo que hace la mayoría; cuando las lecturas ponen un hueco en otro sitio, lo dice la advertencia sobre lecturas que cubren el sitio sin soporte exacto.

Está desactivado por defecto, en la línea de comandos y en la app de escritorio. La única diferencia deliberada con bcftools es un código de ambigüedad en la referencia: bcftools lee R como A, mueve un indel a través de ella y escribe esa A en REF, que entonces no coincide con el FASTA, mientras que get_MNV detiene el indel en la R. La suite tests/leftalign pasa las dos herramientas por miles de indels y registros complejos colocados al azar en las repeticiones del genoma incluido, escritos de varias maneras, y exige el mismo registro de las dos.

Qué aporta get_MNV#

• Agrupación de SNVs por codón en filas SNP, MNV y mixtas SNP/MNV.
• Normalización de inserciones/deleciones de iVar en alelos anclados compatibles con VCF.
• Alineamiento a la izquierda de indels sobre la referencia con --left-align, como hace bcftools norm -f (mira Alineamiento a la izquierda).
• Efectos proteicos de indels y alelos complejos a partir del feature seleccionado, la hebra, la fase del GFF y el offset proteico del transcrito.
• Soporte derivado del BAM para observaciones individuales de SNP/MNV y para eventos exactos de indel o indel complejo.
• Filas complex_indel adicionales solo cuando el haplotipo combinado completo de indel/SNV/MNV se observa en las lecturas del BAM.
• Haplotipos locales leídos de las moléculas. Los eventos cercanos se agrupan en una ventana por proximidad, pero la proximidad solo propone: son las lecturas que cruzan esa ventana las que deciden qué combinaciones se reportan, y cada una sale con el número de moléculas que la llevan. Una combinación que ninguna lectura lleva no se emite, incluida la que es subconjunto de otra que sí se observa, porque una molécula con A, B y C no es evidencia de una molécula con solo A y B. Dos combinaciones que conviven de verdad salen las dos, cada una con su recuento. No hay tope al número de variantes que una ventana admite. Cada candidata se juzga sobre su propio tramo de referencia, y en el caso de una inserción ese tramo llega una base más allá del ancla, porque una inserción vive entre dos bases de referencia y una lectura que termina en el ancla no ha visto nada de la unión.

Reglas de límites#

Las inserciones de VCF son eventos interbase anclados en la base de referencia anterior. Una inserción se considera que solapa un feature solo cuando la secuencia insertada cae entre dos bases de referencia dentro de ese feature. Una inserción anclada en la última base del feature queda fuera de él.

Las deleciones de VCF están ancladas, pero su efecto biológico es el tramo de referencia eliminado. Una deleción anclada justo antes de un CDS/gen afecta a ese feature siempre que las bases eliminadas lo solapen.

Una deleción que va desde las bases de antes del codón de inicio hasta dentro del gen se lee en su posición más 3' a lo largo del transcrito, dentro del primer exón, que es donde la coloca HGVS. Una deleción en una repetición puede escribirse en varios sitios, y alineada a la izquierda, como la escriben --left-align y bcftools norm, borrar el segundo ATG de un gen que empieza por ATGATG justo después de una G empieza en esa G; leída donde está escrita, solo conservaría dos bases codificantes y saldría como corrimiento de marco. Colocada en su posición más 3' queda dentro del gen y se lee Met2del, con los mismos codones y el mismo c. que la misma deleción escrita ahí. La que en esa posición sigue llegando a la base de antes del codón de inicio pierde el inicio. La fila conserva el registro tal como está escrito.

Una inserción o deleción que puede colocarse entera fuera de la proteína no cambia ninguna base codificante y no es del gen: se da como se dan las bases de fuera del gen, intergénica, esté escrito como esté su registro. En el codón de parada ese sitio es su posición más 3' a lo largo del transcrito, como la coloca HGVS. El caso habitual es una serie de As que empieza en un TAA: una A de menos, alineada a la izquierda sobre el codón de parada, se leería si no como un corrimiento de marco en el último residuo de la proteína, cuando el codón de parada sigue ahí de cualquier forma. Una deleción que en su posición más 3' sigue quitando bases del codón de parada se lee tal como está escrita, y una secuencia codificante que no acaba en un límite de codón, como una fila CDS que acaba a mitad de codón en el borde de un exón, no tiene codón de parada que rebasar. En el codón de inicio el sitio es su posición más 5' a lo largo del transcrito: una A de menos o de más en una serie de As que acaba en la A de un ATG deja un ATG donde lo encuentra el ribosoma, y la misma proteína detrás, aunque HGVS escriba el cambio sobre el codón de inicio (c.1del, c.1dup). La que deja un codón de iniciación donde empezaba el codón de inicio, como hace una copia de las primeras bases de la proteína, sigue siendo del gen: el ribosoma lee la copia.

En ningún extremo se mueve un indel más allá de una base que cambia otro registro en las mismas moléculas, porque describiría otro haplotipo. Se da por hecho que los registros comparten las moléculas del indel salvo que digan lo contrario: dos alelos indel de un mismo sitio, como los deja un registro multialélico partido, son alternativas; un llamador puede fasar dos registros en haplotipos distintos; y un registro con menos de la mitad de la frecuencia del indel falta en la mayoría de las moléculas que lo llevan. Un haplotipo fasado que contiene un indel que queda fuera de la proteína no se da como haplotipo propio, porque las filas de sus otras variantes ya dicen lo que hace a la proteína. Un registro complejo que escribe ese indel junto a una sustitución en una sola línea se sigue leyendo como un solo cambio.

Cuando una fila MNV solapa un indel, get_MNV conserva la fila MNV como contexto posicional y marca su efecto en el aminoácido como Unknown con Change Type = Indel overlap. Si el BAM respalda el evento combinado completo, get_MNV emite por separado una fila exacta complex_indel.

El soporte exacto de indels complejos tiene en cuenta el CIGAR. Una lectura debe producir la misma secuencia ALT local y contener los componentes de inserción/deleción esperados. Esto importa en las combinaciones netamente neutras, donde una inserción más una deleción pueden producir la misma secuencia que un MNV simple bajo otro alineamiento.

Notas sobre CDS eucariota#

Para anotaciones eucariotas GFF/GTF, usa --gff-features CDS de modo que get_MNV pueda aprovechar la fase del GFF y los identificadores de transcrito. Cuando las filas CDS llevan transcript_id o Parent, get_MNV construye un modelo de CDS empalmado para cada transcrito y evalúa la agrupación por codones, los efectos del MNV en el aminoácido y el contexto de frameshift del indel sobre esa secuencia codificante completa.

Esto implica que los codones repartidos entre uniones de exones se pueden anotar como un único MNV a nivel de transcrito, y los SNP aguas abajo se marcan como frameshift solo cuando el desplazamiento neto del indel codificante aguas arriba permanece fuera de marco. Las filas CDS sin un identificador de transcrito utilizable mantienen el comportamiento anterior por característica.

Limitaciones actuales#

• Los indels solo se alinean a la izquierda cuando lo pide --left-align; sin él, un indel se empareja con las lecturas donde lo escribió la entrada. Los registros complejos se recortan, no se parten en registros más simples, y no se unen registros de líneas distintas; bcftools norm hace las dos cosas cuando una comparación las necesita.
• Los genotipos, la ploidía, los conjuntos de fase y las verosimilitudes de genotipo no se vuelven a estimar. En datos eucariotas diploides/poliploides, interpreta con cautela los MNV heterocigotos sin fase a nivel de exón salvo que la evidencia del VCF o del BAM confirme que los alelos están en el mismo haplotipo.
• El ensamblaje local de novo queda fuera del alcance. get_MNV solo combina alelos ya presentes en la entrada y, cuando es posible, confirmados por el soporte de lecturas del BAM.
• El descubrimiento de haplotipos locales se limita a ventanas de eventos cercanos; la reconstrucción de haplotipos muy grandes debería seguir a cargo de un llamador dedicado.
• El ensamblaje local completo sigue fuera del alcance: get_MNV reanota los alelos presentes en el archivo VCF/iVar de entrada en lugar de descubrir nuevas variantes candidatas.

Parámetros de ajuste#

Estos flags modifican el comportamiento que tiene en cuenta los indels. Dos de ellos ya no adoptan por defecto el comportamiento original de la herramienta, porque ese comportamiento acertaba menos veces de las que fallaba: el cambio de marco solo se propaga desde un indel aguas arriba mayoritario, y la profundidad del locus del indel se cuenta desde la base ancla.

• --left-align (desactivado por defecto): lleva cada indel a su posición equivalente más a la izquierda sobre la referencia antes de anotarlo, como hace bcftools norm -f, para que coincida con el hueco que reporta el alineador. Mira Alineamiento a la izquierda.

• --frameshift-min-freq <0.0-1.0> (por defecto 0.5): frecuencia alélica mínima que un indel aguas arriba debe alcanzar para desplazar el marco de lectura de los codones SNV/MNV aguas abajo (el marcador (fs) y los tipos de cambio de frameshift). El valor por defecto propaga solo desde un indel mayoritario; pon 0.0 para propagar desde cualquiera. En poblaciones intrahospedador o mixtas, esto evita reetiquetar como frameshifted una sustitución de alta frecuencia aguas abajo a causa de un indel de baja frecuencia aguas arriba que casi con seguridad reside en otra molécula.

  La frecuencia que se compara con el umbral es la que dan las lecturas cuando hay un BAM, que es el mismo número que se reporta como EFREQ. Solo sin BAM ocupa su lugar el AF declarado por el llamador. Importa porque muchos llamadores no escriben AF: mientras la puerta consultaba solo el valor declarado, una entrada así dejaba pasar todos los indels a cualquier umbral. Cuando se proporciona un BAM, la propagación además se fasea por lecturas: para cada indel aguas arriba y un SNV aguas abajo dentro del alcance de un read, se inspeccionan las lecturas que abarcan ambos loci, y el corrimiento de marco no se propaga a ese codón cuando las lecturas que llevan el SNV apenas llevan el indel (están en trans, en moléculas distintas). Es un refinamiento conservador que solo suprime (nunca añade propagación que el gate de frecuencia no haría), y solo aplica donde hay lecturas que abarcan ambos loci; más allá del alcance de un read sigue gobernando el gate de frecuencia.

• --legacy-indel-depth (desactivado por defecto): restringe la profundidad del locus del indel (el denominador de EDP/EFREQ) a las lecturas que abarcan por completo el alelo REF. Por defecto se cuenta desde toda lectura que observe la base de anclaje, lo que evita el subconteo de profundidad y el sesgo de EFREQ que sufre una deleción de varias bases, donde las lecturas con solapamiento parcial quedan excluidas del denominador. Este flag recupera el denominador antiguo, más estrecho.

• --count-mates-separately (desactivado por defecto): cuenta los dos mates de un fragmento como dos observaciones en vez de una molécula. Por defecto son una. Un fragmento es una sola molécula de ADN secuenciada por los dos extremos, así que contar los mates por separado cuenta dos veces el solape, y trata una variante en cada mate como si no tuvieran relación cuando de hecho es la prueba de que van juntas: un codón partido por un intrón puede quedar sin respuesta lectura a lectura y resolverse con el par. Donde los mates solapan y discrepan sobre una base, uno de los dos se equivoca y no hay forma de saber cuál, así que se considera que la molécula no observó esa posición. Los datos single-end no se ven afectados.

  El descubrimiento de haplotipos locales también es a nivel de molécula: cada variante candidata se juzga sobre su propio tramo de referencia, así que un fragmento puede responder por una variante que cubre su primer mate y por otra que cubre el segundo. El conteo exacto del evento no lo es: una lectura tiene que reconstruir el alelo compuesto entero en una observación contigua, y eso no lo puede hacer ningún mate por separado. Un haplotipo local cuyas variantes caen en mates distintos se encuentra pero no se puede respaldar, y el umbral de soporte lo descarta en vez de reportarlo sin evidencia.

• --phased-indel-min-reads <N> (por defecto 2) y --phased-indel-min-freq <0.0-1.0> (por defecto 0.0): soporte mínimo del BAM que una fila de haplotipo indel/complejo con fase debe reunir para ser emitida. Una ventana densa en datos intra-host puede albergar varias combinaciones genuinas a baja frecuencia; elevar estos valores deja solo las bien soportadas.

Cuando un indel con cobertura de lecturas no produce ningún soporte de CIGAR exacto, get_MNV emite una advertencia: por lo general, esto significa que el indel de entrada no está alineado a la izquierda del mismo modo que el BAM (algo habitual en homopolímeros y repeticiones en tándem). Vuelve a ejecutar con --left-align, o normaliza antes la entrada con bcftools norm -f ref.fa, para que el alelo coincida con el alineamiento de las lecturas. Si la advertencia sigue con --left-align, las lecturas ponen el hueco en otro sitio o no llevan el indel.

Referencias#

• FreeBayes: https://github.com/freebayes/freebayes
• bcftools norm: https://samtools.github.io/bcftools/bcftools.html#norm
• GATK HaplotypeCaller: https://gatk.broadinstitute.org/hc/en-us/articles/30332006386459-HaplotypeCaller

Editar esta página en GitHub