Iniciando "FASTQ to Paper"

A partir de aquí comienza el punto de partida práctico de la serie "FASTQ to Paper". En los niveles anteriores de Micro (M01～M30), hemos derivado manualmente el alfabeto de los algoritmos: alineación, indexación, HMM y filogenia. A lo largo de las siguientes 25 entregas, seguiremos cómo estos algoritmos se transportan a través de formatos de archivos prácticos, pasan por diversas herramientas y culminan en las figuras de un artículo científico.

Y el primer archivo de este viaje es precisamente FASTQ. Cada lectura cruda proveniente del secuenciador se encuentra contenida en este formato; si este archivo está defectuoso, toda la tubería de procesamiento posterior se distorsionará. Por lo tanto, no nos apresuremos.

Cuatro líneas constituyen una lectura

Cada lectura en un archivo FASTQ consta exactamente de cuatro líneas.

@SRR1039508.1 HWI-ST177:290:C0TECACXX:1:1101:1225:2130/1
CATTGCTGATACCAATTAGGTCGTGAGGACTGGGCAAGTGTAAAGCT
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHGGFHHHHHHHGGGGHHFB

• 1ª línea (cabecera): Identificador de lectura que comienza con @. La cadena posterior codifica el secuenciador, el flujo de celda, la baldosa y las coordenadas.
• 2ª línea (secuencia): Cadena de caracteres A/C/G/T/N. Es la secuencia de bases real de esta lectura.
• 3ª línea (separador): +. Aunque a veces se repite la cabecera, en la práctica suele ser solo un carácter +.
• 4ª línea (cadena de calidad): Cadena ASCII que tiene exactamente la misma longitud que la secuencia. Cada carácter representa el puntaje de calidad de la base en esa posición.

Estas cuatro líneas se repiten millones o incluso cientos de de millones de veces dentro del archivo. Para un genoma humano con una cobertura de 30x, esto implica unos 100 millones de lecturas; el tamaño del archivo es de entre 30 y 50 GB, incluso comprimido con gzip.

Puntajes Phred: probabilidad de error expresada como logaritmo

¿Por qué la cadena de calidad consiste simplemente en letras como HHHHHH...? Cada carácter en realidad es un número entero codificado en ASCII llamado puntaje Phred.

El puntaje Phred Q se define de la siguiente manera:

• Q10 → probabilidad de error 0,1 (1 de cada 10 lecturas es incorrecta)
• Q20 → probabilidad de error 0,01
• Q30 → probabilidad de error 0,001
• Q40 → probabilidad de error 0,0001

Existe un valor que conviene memorizar por experiencia práctica: Q30 es el "límite inferior para una reanálisis genómico confiable". La mayoría de las lecturas generadas por Illumina NovaSeq superan Q30, y PacBio HiFi oscila entre Q20 y Q40. Si la mayoría de las lecturas tienen un valor Q10 o inferior, esto debe interpretarse como señal de restos de adaptadores o agotamiento de los reactivos.

La razón para usar logaritmos es la conveniencia computacional. Si se almacenara directamente en el archivo una probabilidad de 0,0001, ocuparían cuatro dígitos; con Q40, basta con dos dígitos enteros. Además, al multiplicar probabilidades de múltiples posiciones, en escala logarítmica simplemente se suman. Esto es exactamente la misma técnica de probabilidades logarítmicas utilizada en informática.

Phred33 vs Phred64 — historia de dos codificaciones

Ahora viene la verdadera trampa: hay dos estándares para mapear el número entero Q a caracteres ASCII.

Es decir, si en el archivo se observa el carácter H,

• Si es Phred33, entonces Q = ord('H') − 33 = 72 − 33 = Q39 (casi perfecto)
• Si es Phred64, Q = ord('H') − 64 = 72 − 64 = Q8 (casi un error)

Esto significa que el mismo carácter tiene significados completamente diferentes en dos mundos distintos. El problema surge aquí: descargar datos antiguos de SRA e insertarlos directamente en la tubería puede provocar una situación en la que se pierdan dos días preguntándose "¿por qué es esta la calidad?".

Principios de detección automática

¿Cómo sabemos qué codificación tiene un archivo? La regla práctica es la siguiente:

1. Examine el valor mínimo del carácter de calidad de las primeras decenas de miles de lecturas.
2. Si el valor mínimo está entre los códigos ASCII 33 y 58 → Phred33 (porque los caracteres con valores ≤ 58 no existen en Phred64).
3. Si solo aparecen valores mínimos ASCII ≥ 64 → sospeche de Phred64.
4. En caso de duda, decida según el año de publicación del conjunto de datos (después de 2012 = casi con certeza Phred33).

El mensaje Encoding: Sanger / Illumina 1.9 que aparece en la consola cuando se ejecuta FastQC por primera vez es exactamente esta determinación.

Analizador de transmisión en Python de 20 líneas

Como el archivo pesa 30 GB, no se puede cargar completamente en la memoria. La práctica estándar es procesarlo en modo streaming, leyendo cuatro líneas a la vez.

import gzip

def iter_fastq(path):
    opener = gzip.open if path.endswith(".gz") else open
    with opener(path, "rt") as f:
        while True:
            header = f.readline()
            if not header:
                return
            seq = f.readline().rstrip()
            plus = f.readline()
            qual = f.readline().rstrip()
            yield header.rstrip(), seq, qual

def phred33_to_q(qual_str):
    return [ord(c) - 33 for c in qual_str]

# Ejemplo: calidad media de las primeras 100 lecturas
for i, (hdr, seq, qual) in enumerate(iter_fastq("reads.fastq.gz")):
    if i >= 100:
        break
    qs = phred33_to_q(qual)
    print(f"read {i}: len={len(seq)}, meanQ={sum(qs)/len(qs):.1f}")

Esto es todo. Detección automática de gzip, transmisión en streaming, resta del valor ASCII 33. Incluso si el archivo tiene 3 TB, este bucle consume memoria constante.

Validación de que la longitud de la secuencia sea igual a la longitud de la cadena de calidad

El código defensivo que siempre se incluye en la práctica al escribir un analizador FASTQ por primera vez es el siguiente:

for hdr, seq, qual in iter_fastq(path):
    if len(seq) != len(qual):
        raise ValueError(f"length mismatch: {hdr}")

Estos eventos —errores del secuenciador, pérdida de transmisión o fragmentación parcial de gzip— son capturados por esta aserción. Es mucho mejor que falle visiblemente aquí, en lugar de morir silenciosamente en las etapas posteriores del pipeline (BWA·GATK).

Sensibilidad práctica — Tres indicadores

Al abrir el informe de FastQC aparecen más de 20 elementos. Al principio no es necesario revisar todos; basta con observar tres.

1. Calidad de la secuencia por base: Avanzando de izquierda a derecha, observe cuánto tiempo se mantiene verde la caja Q30. Si la parte posterior del read colapsa por debajo de Q20, será necesario realizar un recorte 3′.
2. Secuencias sobrerepresentadas: Detección de restos de adaptadores. Si una secuencia específica se repite más del 5% de las veces en los reads, es una señal de contaminación por adaptadores o cebadores de PCR.
3. Niveles de duplicación de secuencias: Señal de sesgo de amplificación por PCR. En RNA-seq la duplicación natural es esperada, pero en WGS una duplicación superior al 30% indica que probablemente sea necesario volver a preparar la librería.

En el siguiente episodio (S02) resolveremos estos tres puntos en la práctica utilizando FastQC y herramientas de recorte automático (fastp · Trim Galore).

Abrir archivos reales en Colab

La forma más rápida de familiarizarse con los datos FASTQ reales en un entorno gratuito es utilizar muestras pequeñas de SRA. En Colab, ejecute solo las dos celdas siguientes.

!pip install pyfastx -q
!wget -q https://sra-pub-src-1.s3.amazonaws.com/SRR1039508/SRR1039508_1.fastq.gz -O r1.fq.gz

import pyfastx
fq = pyfastx.Fastq("r1.fq.gz")
print(f"Lecturas totales: {len(fq):,}, longitud media: {fq.avglen:.1f} bp")
print(f"GC%: {fq.gc_content:.1f}, codificación: {fq.phred}")

# Resumen de calidad de las primeras tres lecturas
for name, seq, qual in fq[:3]:
    qs = [ord(c) - fq.phred for c in qual]
    print(name, "meanQ=", sum(qs) / len(qs))

Cuando fq.phred muestra 33 o 64, la detección automática de la codificación ha funcionado.

Nota (entorno práctico): el nivel gratuito T4 de Colab es suficiente. Si una descarga de SRA es lenta, también puede utilizar el espejo de EBI ENA (https://www.ebi.ac.uk/ena/browser/).

Correspondencia de conceptos de informática

FASTQ combina en la práctica tres conceptos de informática.

• Codificación con desplazamiento ASCII: desplaza el punto de inicio para representar enteros como caracteres. Phred33/64 son simplemente convenciones con desplazamientos de 33 y 64.
• Probabilidad logarítmica: como las probabilidades se multiplican con frecuencia, se almacenan en escala logarítmica. Q = −10log₁₀P.
• Analizador de flujo: procesa el archivo mediante un iterador línea por línea porque este puede superar la memoria disponible. Es la misma filosofía que las tuberías de Unix.

El episodio de DryBench «Codificación ASCII y análisis de texto» (ascii-encoding-and-text-parsing) presenta estos conceptos desde la informática pura. FASTQ puede considerarse uno de los primeros casos en los que se materializan en la práctica biológica.

Conclusión: hasta la próxima entrega

Ahora podemos abrir el archivo, extraer el vector de calidad de cada lectura y detectar automáticamente la codificación. En el siguiente episodio, S02, eliminaremos los adaptadores y recortaremos las colas 3′ de baja calidad para preparar la entrada del pipeline. Después, en S03, alinearemos las lecturas contra el genoma de referencia con BWA-MEM.

Si desea profundizar

El texto principal es una reconstrucción propia de BPD. Para afianzar los principios, recomendamos los siguientes recursos.

• Illumina Knowledge Base — Understanding Illumina Quality Scores (documento técnico oficial). Destacan sus gráficas empíricas de calibración de puntuaciones Q.
• Sección FASTQ del NCBI SRA Handbook: cronología oficial de la evolución de la codificación.
• Ewing & Green (1998), Base-calling of automated sequencer traces using phred. II. Genome Research 8:186. Artículo original del algoritmo Phred y punto de partida de las puntuaciones Q actuales.
• Broad Institute BroadE YouTube — Clase de 15 minutos sobre la química de Illumina. Se visualizan las causas ópticas y químicas por las que disminuye la calidad en el extremo 3′.

Solo quienes han abierto archivos y escrito manualmente sus propios parsers pueden sentirse seguros al tomar decisiones de recorte en la siguiente entrega. Enciende Colab y descarga una secuencia SRR.