¿Por qué es problemático el duplicado de PCR?

En la edición anterior, mencionamos brevemente por qué MarkDuplicates se encuentra en la tercera posición del flujo de trabajo. En esta edición, explicaremos esa razón desde una perspectiva estadística.

Durante la preparación de la librería, se realizan entre 20 y 30 ciclos de PCR para amplificar las fragmentos originales de ADN (inserts) hasta niveles suficientes para su secuenciación. Esta amplificación es inherentemente sesgada.

• Sesgo de GC: Las regiones con extremos de contenido de GC se amplifican menos.
• Sesgo de longitud: Los inserts más cortos se amplifican mejor.
• Sesgo del número inicial de copias: Si ciertos originales se amplifican por casualidad en los primeros ciclos, ese seso se amplifica exponencialmente.

Como resultado, las lecturas de secuenciación contienen una mezcla de copias de un solo insert original que ha sido replicado cientos de veces. Contar estas copias como observaciones independientes distorsiona las estadísticas de cobertura y lleva a llamadas de variantes con una confianza errónea. MarkDuplicates actúa como la primera barrera para controlar este sesgo.

¿Qué lecturas se consideran duplicadas?

Las reglas de MarkDuplicates son sorprendentemente simples.

Las lecturas que comparten la misma coordenada de referencia (misma posición, mismo sentido y mismo punto de inicio del CIGAR) se consideran duplicados derivados de un solo original.

En el caso de paired-end, los puntos de inicio de R1 y R2 deben ser idénticos. Si solo las coordenadas de R1 coinciden pero las de R2 difieren, provienen de orígenes diferentes.

Selección del representante

De cada grupo de coordenadas idénticas, se conserva una sola lectura como original y se marcan las demás como duplicadas. El criterio para elegir al representante es la lectura con la suma más alta de calidades de base totales. Es decir, se conserva la copia de mayor calidad y se asigna la marca 1024 a las demás.

• Regla en el código: flag |= 0x400 (=1024).
• Las lecturas no se eliminan del archivo; permanecen intactas.
• Las herramientas posteriores, como HaplotypeCaller, ignoran estas lecturas en sus cálculos estadísticos al detectar esta marca.

¿Por qué solo marcar y no eliminar?

Existen tres razones principales:

1. Reversibilidad: Si más adelante se desea reevaluar con otro algoritmo, basta con quitar la marca.
2. Cálculo de indicadores de control de calidad (QC): Para calcular la tasa de duplicación, las lecturas duplicadas deben permanecer en el archivo.
3. Algunos análisis utilizan los duplicados: En RNA-seq o en análisis de cobertura muy baja, los duplicados pueden utilizarse como señal biológica.

Esta convención de "solo marcar" es la base de la flexibilidad de los flujos de trabajo prácticos.

Duplicados ópticos — La identidad del duplicado óptico

Entre las lecturas que se alinean en la misma coordenada, hay aquellas que están físicamente cercanas en el flujo de célula (flowcell). Esto ocurre cuando un solo clúster es interpretado erróneamente por el secuenciador como dos clústers distintos. No son duplicados de PCR, sino duplicados ópticos (optical duplicates).

MarkDuplicates extrae las coordenadas x e y del nombre de la lectura (por ejemplo, HWI-ST177:290:C0TECACXX:1:1101:1225:2130) y clasifica por separado como duplicados ópticos a aquellos pares de lecturas que se encuentran dentro de un radio de 100 píxeles.

riD_TILE : 1101, x: 1225, y: 2130
riD_TILE : 1101, x: 1230, y: 2135   # proximidad → posible duplicado óptico

La opción OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 es el umbral. Plataformas recientes como NovaSeq X deben aumentarse a aproximadamente 2500. Si se establece este valor incorrectamente, los conteos de duplicados ópticos resultarán anómalos.

Comando de ejecución

gatk MarkDuplicates \
    -I S1.sorted.bam \
    -O S1.dedup.bam \
    -M S1.metrics.txt \
    --VALIDATION_STRINGENCY LENIENT \
    --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 \
    --TMP_DIR /data/tmp
samtools index S1.dedup.bam

Solo necesitas seleccionar unas tres opciones.

• -M metrics.txt: Archivo de informe (por ejemplo, la tasa de duplicación por biblioteca).
• --VALIDATION_STRINGENCY LENIENT: No detenerse aunque haya lecturas que infrinjan el estándar SAM.
• --OPTICAL_DUPLICATE_PIXEL_DISTANCE: Valores específicos de la plataforma.
• --TMP_DIR: Especificar particiones grandes.

Cómo leer metrics.txt

LIBRARY   UNPAIRED_READS  READ_PAIRS  UNMAPPED  UNPAIRED_DUP  READ_PAIR_DUP  READ_PAIR_OPT_DUP  PERCENT_DUPLICATION  ESTIMATED_LIBRARY_SIZE
lib1      100             49,999,900  200,000   50            2,500,000      15,000             0.05                 950,000,000

Revisamos siempre estos dos indicadores.

• PERCENT_DUPLICATION: 0.05 = 5%. En WGS, el rango normal es del 3 al 10%. Si supera el 20%, se recomienda preparar una nueva biblioteca.
• ESTIMATED_LIBRARY_SIZE: Estimación de inserciones únicas. Refleja la complejidad de la biblioteca; cuanto mayor, mejor.

Estas dos cifras son el informe práctico de calidad de la biblioteca.

Cálculo manual — Dos escenarios

Supongamos que 10 lecturas se alinean en la misma coordenada.

Escenario A: WGS con cobertura de 30x

• Genoma humano de 3×10⁹ pb × 30x / lectura de 150 pb = 600 millones de pares de lecturas.
• La probabilidad estadística de que 10 lecturas se alineen en la misma coordenada es baja (cobertura media 30, desviación estándar 5).
• En la práctica, esto suele deberse a señales de duplicación óptica o por PCR.

Escenario B: Secuenciación dirigida (amplicones)

• Biblioteca amplificada mediante PCR solo para regiones específicas de ciertos genes.
• Es normal que haya miles de lecturas en la misma coordenada; las lecturas originales provienen realmente de esa coordenada y no son artefactos de PCR.
• En este caso, se omite MarkDuplicates o se utiliza una evaluación de duplicados basada en UMI.

Es decir, las reglas cambian completamente según el diseño experimental.

UMI — La solución fundamental al problema de la complejidad de la biblioteca

Determinar las duplicaciones por PCR basándose únicamente en la coordenada conduce a falsos positivos. Esto es especialmente cierto en muestras de bajo volumen o secuenciación dirigida. La solución fundamental es UMI (Identificador Molecular Único).

Se añade un código aleatorio de 8 a 12 pb a cada inserción original durante la etapa de preparación de la biblioteca. Si el encabezado de la lectura contiene el UMI, se determina la identidad del original mediante la combinación coordenada + UMI. Incluso si hay 10 lecturas en la misma coordenada, si los UMIs son diferentes, se trata de originales distintos.

• fgbio: Herramienta de procesamiento de UMI de Fulcrum Genomics.
• umi_tools: Herramienta de Python desarrollada por cgat.
• GATK4 también incluye MarkDuplicatesWithMateCigar y UmiAwareMarkDuplicatesWithMateCigar.

En la cuantificación de genes de baja expresión en RNA-seq o en la detección de VAF bajo en biopsias líquidas, el uso de UMI se ha convertido en el estándar. Cabe tener en cuenta que el flujo de trabajo básico de las Buenas Prácticas de GATK no está diseñado para trabajar con UMIs.

Práctica pequeña en Colab

!apt-get install -y samtools default-jre >/dev/null
!wget -q https://github.com/broadinstitute/picard/releases/download/3.1.0/picard.jar
!wget -q https://github.com/samtools/samtools/raw/master/examples/toy.bam
!java -jar picard.jar MarkDuplicates \
    I=toy.bam \
    O=toy.dedup.bam \
    M=toy.metrics.txt \
    VALIDATION_STRINGENCY=LENIENT
!cat toy.metrics.txt | head -15

Como toy.bam es un ejemplo muy pequeño, sus métricas estadísticas son artificiales. Aun así, permite familiarizarse con la lectura de metrics.txt.

Mapeo CS

• Agrupación mediante hash: agrupa las lecturas en un mapa hash usando (coordenada de referencia, hebra) como clave y selecciona una representante por grupo.
• Selección del máximo: toma el argmax de la suma de calidades de base dentro del grupo.
• Marcado con flag frente a eliminación: convención reversible equivalente a un borrado lógico.

Consulte el episodio de DryBench «Eliminación de duplicados mediante hash» (hash-based-deduplication).

Conclusión

Aunque MarkDuplicates parece una herramienta sencilla, es una etapa esencial cuyo objetivo estadístico es controlar el sesgo de PCR. En el siguiente episodio (S08) pasaremos a BQSR y estudiaremos con detalle estadístico por qué deben recalibrarse las puntuaciones Phred y cómo se aprende esa recalibración.

Para profundizar

• Documentación oficial de Picard MarkDuplicates (https://broadinstitute.github.io/picard/): fuente primaria de las opciones.
• Ebbert M. et al. (2016), Evaluating the necessity of PCR duplicates removal from next-generation sequencing data and a comparison of approaches. BMC Bioinformatics 17:239. Estudio empírico sobre la necesidad de eliminar duplicados.
• Smith T. et al. (2017), UMI-tools: modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy. Genome Research 27:491. Artículo original de UMI-tools.
• Broad Institute BroadE — MarkDuplicates deep dive, clase en YouTube que muestra visualmente el origen óptico de los duplicados.
• Documentación oficial de fgbio (https://github.com/fulcrumgenomics/fgbio): estándar práctico para pipelines con UMI.

Una vez entendido intuitivamente por qué el flag 1024 no implica eliminación, las opciones de las etapas posteriores del pipeline resultan mucho más claras. Continuemos con BQSR en el siguiente episodio.