Mapeo Genómico e Identificación de Variantes Patogénicas
Niko Alain Cruz Sancén
06 de August de 2026
1 Objetivo de Aprendizaje
En esta sesión vamos a integrar un flujo de trabajo bioinformático completo usando secuencias cortas (FASTQ) reales. El objetivo es que aprendas a evaluar la calidad de las lecturas, limpiarlas, alinearlas contra el genoma humano y, lo más importante, detectar variantes genéticas con potencial patogénico.
Al final, entenderás cómo estas herramientas nos permiten identificar mutaciones (como las del gen PAH en la fenilcetonuria) para tomar decisiones en la práctica clínico-nutricional.
2 1. Carga de Librerías
Primero, traemos nuestras herramientas de trabajo. Cada librería tiene un propósito específico en nuestro flujo:
QuasR y Rbowtie: Son el motor de alineamiento; nos ayudan a pegar nuestras lecturas contra el genoma.
Rsamtools y GenomicRanges: Nos permiten manipular los archivos genómicos pesados y buscar por coordenadas específicas.
ShortRead: Es genial para leer y manipular archivos FASTQ.
3 2. Preparación de los Datos
Para esta sesión práctica, nuestro profesor ya preparó un subset de datos clínicos con las lecturas necesarias para el análisis (250,000 secuencias). Esto nos ahorrará tiempo de cómputo y ancho de banda.
Instrucciones previas:
Entra a este enlace de Google Drive: Datos Proyecto PKU
Descarga los archivos
subset_PKU_1.fastq.gzysubset_PKU_2.fastq.gz.Crea una carpeta llamada
proyecto_PKUen el mismo lugar donde guardaste este archivo de R y mueve los archivos descargados ahí.
Una vez que tengas los archivos en tu carpeta, corre el siguiente bloque para verificar que R los detecte correctamente:
# Definimos las rutas donde deben estar los archivos que descargaste
subset_fwd <- "proyecto_PKU/subset_PKU_1.fastq.gz"
subset_rev <- "proyecto_PKU/subset_PKU_2.fastq.gz"
# Verificamos que los archivos existan
if(file.exists(subset_fwd) & file.exists(subset_rev)) {
cat("¡Excelente! Los archivos de secuenciación están listos para analizarse.\n")
} else {
stop("R no encuentra los archivos. Revisa que la carpeta se llame 'proyecto_PKU' y esté en el mismo directorio que tu script.")
}## ¡Excelente! Los archivos de secuenciación están listos para analizarse.4 3. Preparación del Genoma de Referencia
Para buscar la fenilcetonuria, no necesitamos todo el genoma. El gen responsable (PAH) vive en el Cromosoma 12. Vamos a descargar exclusivamente este cromosoma desde la base de datos de UCSC y lo descomprimiremos para que nuestro alineador pueda leerlo fácilmente.
dir.create("referencia", showWarnings = FALSE)
fasta_referencia_gz <- "referencia/chr12.fa.gz"
fasta_referencia <- "referencia/chr12.fa"
# Descargamos desde UCSC
if(!file.exists(fasta_referencia) && !file.exists(fasta_referencia_gz)){
download.file("[https://hgdownload.soe.ucsc.edu/goldenPath/hg38/chromosomes/chr12.fa.gz](https://hgdownload.soe.ucsc.edu/goldenPath/hg38/chromosomes/chr12.fa.gz)",
destfile = fasta_referencia_gz, method = "wget")
}
# Descomprimimos el archivo para que QuasR pueda leerlo
if(file.exists(fasta_referencia_gz) && !file.exists(fasta_referencia)){
system("gunzip -f referencia/chr12.fa.gz")
}5 4. Control de Calidad (QC) de Lecturas Crudas
Antes de alinear, tenemos que revisar qué tan “sucias” vienen
nuestras secuencias. Usamos fastqcr para generar un reporte
visual súper amigable en formato HTML.
Presta atención a las gráficas de calidad (Phred scores): todo lo que esté en la zona verde es bueno, lo rojo hay que filtrarlo.
# Nota: Este bloque no se ejecuta al tejer el PDF/HTML.
# Córrelo directamente en tu consola de R para ver los resultados en tu navegador.
library(fastqcr)
#fastqc_install() # Descomentar la primera vez que se use en la computadora
dir.create("proyecto_PKU/fastqc_results", showWarnings = FALSE)
fastqc(fq.dir = "proyecto_PKU", qc.dir = "proyecto_PKU/fastqc_results")Antes de cualquier alineamiento, debemos conocer el estado de nuestros datos crudos. Generaremos un reporte de calidad en formato PDF.
Revisa el archivo QC_Report_Raw.pdf en tu directorio
para observar la distribución de calidad (Phred scores) y el contenido
de GC.
6 5. Preprocesamiento y Filtrado
Basándonos en el QC previo, vamos a limpiar nuestros datos. Le diremos al programa dos cosas clave:
minLength = 50: Si la secuencia es muy cortita, bórrala. No nos sirve.nBases = 0: Si la máquina de secuenciación dudó y puso una letra ‘N’ (ambigua), descarta esa lectura. Queremos datos 100% seguros.
# Definimos los archivos de entrada y salida
archivos_raw <- c("proyecto_PKU/subset_PKU_1.fastq.gz",
"proyecto_PKU/subset_PKU_2.fastq.gz")
archivos_filtrados <- c("proyecto_PKU/filtrados/subset_PKU_1_filtered.fastq.gz",
"proyecto_PKU/filtrados/subset_PKU_2_filtered.fastq.gz")
dir.create("proyecto_PKU/filtrados", showWarnings = FALSE)
# Ejecutamos la limpieza SOLO si los archivos no existen ya
if(!file.exists(archivos_filtrados[1]) | !file.exists(archivos_filtrados[2])) {
cat("Filtrando secuencias...\n")
preprocessReads(filename = archivos_raw,
outputFilename = archivos_filtrados,
truncateStartBases = 0,
truncateEndBases = 0,
minLength = 50,
nBases = 0)
} else {
cat("Los archivos ya estaban filtrados. ¡Saltando este paso para ahorrar tiempo de cómputo!\n")
}## Los archivos ya estaban filtrados. ¡Saltando este paso para ahorrar tiempo de cómputo!# Creamos un archivo de texto (índice) que apunte a nuestras secuencias ya limpias
archivo_muestras_filtered <- "proyecto_PKU/samples_filtered.txt"
muestras_filtradas <- data.frame(
FileName1 = archivos_filtrados[1],
FileName2 = archivos_filtrados[2],
SampleName = "Paciente_PKU"
)
write.table(muestras_filtradas, archivo_muestras_filtered, sep = "\t", row.names = FALSE, quote = FALSE)7 6. Control de Calidad Post-Filtrado
Es una excelente práctica clínica volver a correr el análisis de calidad sobre los datos limpios para asegurar que el filtro hizo su trabajo.
# Cargamos la librería específicamente para este paso
library(fastqcr)
# Creamos la carpeta
dir.create("proyecto_PKU/fastqc_results_post", showWarnings = FALSE)
# Ejecutamos el reporte
cat("Generando reporte de calidad post-filtrado...\n")
fastqc(fq.dir = "proyecto_PKU/filtrados", qc.dir = "proyecto_PKU/fastqc_results_post")8 7. Alineamiento y Mapeo Genómico
Aquí ocurre la magia. Vamos a empatar nuestros “pedacitos” de ADN (lecturas) contra el mapa completo del Cromosoma 12.
Detalle técnico importante: Como sacamos un muestreo
aleatorio en el Paso 2, perdimos la sincronía de nuestras secuencias
pareadas (Forward y Reverse). Para evitar errores, le diremos a
qAlign que trate estas secuencias como independientes
(Single-End).
# 1. Creamos el índice en modo Single-End
archivo_indice <- "samples_raw_single.txt"
muestras_single <- data.frame(
FileName = archivos_filtrados,
SampleName = c("Paciente_PKU", "Paciente_PKU")
)
write.table(muestras_single, archivo_indice, sep = "\t", row.names = FALSE, quote = FALSE)
# 2. Alineamos contra el cromosoma 12
cat("Iniciando el alineamiento (esto puede tomar unos momentos)...\n")## Iniciando el alineamiento (esto puede tomar unos momentos)...proyecto <- qAlign(sampleFile = archivo_indice,
genome = fasta_referencia,
aligner = "Rbowtie")
# 3. Revisamos cuántas lecturas "pegaron" con éxito
print(alignmentStats(proyecto))## seqlength mapped unmapped
## Paciente_PKU:genome 133275309 26170 4726089 8. Identificación de Variantes (Variant Calling)
Ya tenemos un archivo .bam (Binary Alignment Map).
Ahora, le pediremos a R que se asome específicamente en las coordenadas
del gen PAH, y nos diga si la secuencia de nuestro
paciente tiene letras distintas a las del genoma de referencia.
# Extraemos la ruta correcta del archivo BAM
bam_file <- alignments(proyecto)$genome$FileName[1]
cat("Analizando el archivo BAM:", bam_file, "\n")## Analizando el archivo BAM: /home/alainsancen/dipplomado_bioinformatica/modulo2/sesion4/proyecto_PKU/filtrados/subset_PKU_1_filtered_cf9c465fc84b.bam# Configuramos la exigencia de calidad (Phred score mínimo de 20)
p_param <- PileupParam(max_depth = 250,
min_base_quality = 20,
min_mapq = 20)
# Le damos a R las coordenadas exactas del gen PAH (Cromosoma 12)
region_interes <- GRanges(seqnames = "chr12",
ranges = IRanges(start = 102836000, end = 102917000))
# Extraemos las variantes encontradas
variantes <- pileup(bam_file,
scanBamParam = ScanBamParam(which = region_interes),
pileupParam = p_param)
# Mostramos el resultado
cat("\n¡Variantes detectadas en la región de interés!\n")##
## ¡Variantes detectadas en la región de interés!# Cargamos la librería de visualización
library(ggplot2)
# Verificamos que la tabla no esté vacía antes de graficar
if(nrow(variantes) > 0) {
# Creamos el gráfico estilo lollipop (puntos y líneas)
grafica <- ggplot(variantes, aes(x = pos, y = count, color = nucleotide)) +
# Línea base para saber dónde cae exactamente en el eje X
geom_segment(aes(x = pos, xend = pos, y = 0, yend = count), alpha = 0.5, show.legend = FALSE) +
# El punto grande que resalta el color y la altura
geom_point(size = 4, alpha = 0.9) +
scale_color_brewer(palette = "Set1") +
theme_minimal() +
labs(
title = "Frecuencia de Nucleótidos por Posición (Gen PAH)",
subtitle = "Cobertura de lecturas en el Cromosoma 12",
x = "Posición Genómica Exacta",
y = "Número de Lecturas (Count / Cobertura)",
color = "Nucleótido detectado"
) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold", size = 14)
)
print(grafica)
} else {
cat("No hay suficientes variantes detectadas en esta muestra para graficar.\n")
}
# Cargamos la librería para manipular datos
library(dplyr)
# Transformamos la tabla original al formato que buscas
tabla_resumen <- variantes %>%
# Agrupamos por cromosoma y posición exacta
group_by(seqnames, pos) %>%
summarise(
# Juntamos todas las letras encontradas en esa posición separadas por comas
alelos_detectados = paste(unique(nucleotide), collapse = ","),
# Sumamos las lecturas totales para saber la profundidad real
cobertura_total = sum(count),
.groups = 'drop'
) %>%
# Creamos tu columna con el formato específico para PolyPhen (agregando REF/)
mutate(
variante_formato = paste0(seqnames, ":", pos, " REF/", alelos_detectados)
) %>%
# Le decimos a R que use estrictamente el select de dplyr
dplyr::select(variante_formato, seqnames, pos, alelos_detectados, cobertura_total)
# Mostramos el resultado a los alumnos
cat("\nResumen de variantes listas para análisis en bases de datos clínicas:\n")##
## Resumen de variantes listas para análisis en bases de datos clínicas:10 Conclusión
Los pasos que harían en clase son:
Entrar a la página oficial de ClinVar (NCBI).
En la barra de búsqueda, se puede escribir el nombre del gen y la posición para filtrar. Por ejemplo:
PAH[gene] AND 102840374. Asegurándonos de que esté en la versión del genoma humano GRCh38, que fue la que descargamos.ClinVar les arrojará una tabla con los reportes de otros hospitales en el mundo.
¿Qué pueden encontrar ahí?
Pathogenic / Likely Pathogenic: Confirmado, esta mutación rompe la enzima. Aquí es donde tú, como nutriólogo, entras a diseñar la dieta libre de fenilalanina.
Benign / Likely Benign: Es un simple polimorfismo. El paciente tiene una letra diferente, pero su cuerpo procesa la fenilalanina sin problemas. No hay que cambiarle la dieta.
VUS (Variant of Uncertain Significance): No es claro el efecto de portar esa mutación. Se mantiene bajo vigilancia.