Alineamiento

Configuración del Entorno y Dependencias

Antes de comenzar con los análisis bioinformáticos, este bloque de código evaluará si tu entorno de R cuenta con todas las dependencias necesarias. Si te falta algún paquete, el código lo descargará e instalará automáticamente (esto puede tardar un par de minutos la primera vez).

# 1. Definir los paquetes requeridos
paquetes_cran <- c("seqinr", "rentrez","ape","ggplot2")
paquetes_bioc <- c("Biostrings","pwalign","msa","ggmsa")

# 2. Verificar e instalar paquetes de CRAN (repositorio estándar)
for (pkg in paquetes_cran) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
  }
}

# 3. Verificar e instalar herramientas de Bioconductor (repositorio biológico)
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
for (pkg in paquetes_bioc) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    BiocManager::install(pkg, update = FALSE)
  }
}

# 4. Cargar todas las librerías al entorno
library(Biostrings)

## Warning: package 'S4Vectors' was built under R version 4.5.3

library(seqinr)
library(rentrez) # Herramienta para conectar con la base de datos del NCBI
library(pwalign)
library(msa)
library(ape)
library(ggmsa)
library(ggplot2)

1. Alineamiento Global y Local de Secuencias de ADN

Ejercicio

Encuentre el alineamiento global óptimo y los alineamientos locales de las siguientes dos secuencias utilizando el paquete Biostrings en el entorno de R:

Secuencia 1: "ATGAAATATACAAGTTATATCTTGGCTTTTCAGCTCTGCATCGTTTTGGGTTCTCTTGGC"
Secuencia 2: "ATGAACGCTACACACTGCATCTTGGCTTTGCAGCTCTTCCTCATGGCTGTTTCTGGCTGT"

(Nota: Para propósitos de demostración en el código, se utilizarán fragmentos representativos de estas secuencias).

Solución

Utilizaremos la función pairwiseAlignment del paquete Biostrings para realizar tanto el alineamiento global óptimo como el local. La puntuación para coincidencias (match) es de +1 y para discordancias (mismatch) es de -1. La penalización por apertura de brecha (gap opening) es de 5 y la penalización por extensión de brecha (gap extension) es de 2.

# Cargar la librería necesaria
library(Biostrings)

# Definir las secuencias de ADN
s1 <- DNAString("ATGAAATATACAAGTTATAT")
s2 <- DNAString("ATGAACGCTACACACTGCAT")

# Crear la matriz de sustitución de nucleótidos
mat <- pwalign::nucleotideSubstitutionMatrix(match = 1, mismatch = -1, baseOnly = TRUE)

# 1. Alineamiento Global
globalAlign <- pairwiseAlignment(s1, s2, substitutionMatrix = mat, gapOpening = 5, gapExtension = 2)

# 2. Alineamiento Local
localAlign <- pairwiseAlignment(s1, s2, type = "local", substitutionMatrix = mat, gapOpening = 5, gapExtension = 2)

Resultados del Alineamiento

A continuación se muestran los resultados obtenidos para el alineamiento global:

globalAlign

## Global PairwiseAlignmentsSingleSubject (1 of 1)
## pattern: ATGAAATATACAAGTTATAT
## subject: ATGAACGCTACACACTGCAT
## score: 4

A continuación se muestran los resultados obtenidos para el alineamiento local:

localAlign

## Local PairwiseAlignmentsSingleSubject (1 of 1)
## pattern: [1] ATGAAATATACA
## subject: [1] ATGAACGCTACA
## score: 6

2. Descarga de Secuencias Reales y Análisis de Homología Evolutiva

Ejercicio

En esta sección, ilustraremos el concepto de homología evolutiva comparando la secuencia de la proteína GTP asa entre tres especies: Humano (Homo sapiens), Ratón (Mus musculus) y Guajolote (Meleagris gallopavo).

Utilizaremos el paquete rentrez para descargar las secuencias reales del NCBI y generaremos despliegues gráficos (Dot Plots) con parámetros relajados para observar las diferencias en la conservación estructural.

ID Humano (Homo sapiens): CAG38816.1
ID Ratón (Mus musculus): NP_001399410.1
ID Chimpancé (Meleagris gallopavo): XP_085080516

Solución

Primero, utilizaremos el paquete rentrez para consultar el servidor del NCBI, descargar la información biológica en formato FASTA y cargarla directamente a nuestro entorno de R.

# Descargar secuencias FASTA desde el NCBI
fasta_data <- entrez_fetch(db = "protein", 
                           id = c("CAG38816.1", "NP_001399410.1", "XP_085080516"), 
                           rettype = "fasta")

# Guardar temporalmente en un archivo
writeLines(fasta_data, "HRas.fasta")

# Leer el archivo FASTA
secuencias <- read.fasta("HRas.fasta", seqtype = "AA")

# Extraer las secuencias individuales
humano <- secuencias[[1]]
raton <- secuencias[[2]]
guajolote <- secuencias[[3]]

cat("Descarga exitosa.\nLongitud Humano:", length(humano), "AA\nLongitud Ratón:", length(raton), "AA\nLongitud Guajolote:", length(guajolote), "AA\n")

## Descarga exitosa.
## Longitud Humano: 189 AA
## Longitud Ratón: 189 AA
## Longitud Guajolote: 189 AA

A continuación, generamos los Dot Plots. Utilizaremos los mismos parámetros relajados (wsize = 1, nmatch = 1) para ambas comparaciones, de modo que el contraste visual sea evidente.

Comparación 1: Humano vs. Ratón

seqinr::dotPlot(humano, raton, 
        wsize = 1,    
        wstep = 1,    
        nmatch = 1,   
        pch = 15,     
        cex = 0.4,    
        main = "Homología: HRas (Humano vs Ratón)",
        xlab = "Secuencia Humana (CAG38816.1)", 
        ylab = "Secuencia de Ratón (NP_001399410.1)")

Comparación 2: Humano vs. Guajolote

seqinr::dotPlot(humano, guajolote, 
        wsize = 1,    
        wstep = 1,    
        nmatch = 1,   
        pch = 15,     
        cex = 0.4,    
        main = "Homología: HRas (Humano vs Guajolote)",
        xlab = "Secuencia Humana (CAG38816.1)", 
        ylab = "Secuencia de Guajolote (XP_085080516)")

Para dar el siguiente paso y realizar un Alineamiento Múltiple de Secuencias (MSA), vamos a utilizar el paquete msa. Este paquete es fantástico porque te da acceso directo a los algoritmos más famosos del mundo para esta tarea, como ClustalW, ClustalOmega y MUSCLE, directamente desde tu consola.

# La función c2s() de seqinr une los caracteres ('M', 'T', 'E'...) en una sola palabra ("MTE...")
seq_humano <- seqinr::c2s(humano)
seq_raton <- seqinr::c2s(raton)
seq_guajolote <- seqinr::c2s(guajolote)

# Creamos el set de secuencias con sus respectivos nombres para que la tabla sea legible
secuencias_hras <- AAStringSet(c(
  "Humano" = seq_humano,
  "Raton" = seq_raton,
  "Guajolote" = seq_guajolote
))

Los algoritmos de alineamiento en Bioconductor no leen letras sueltas (como las tiene seqinr), sino que necesitan un bloque de texto continuo empaquetado en un formato especial llamado AAStringSet (Conjunto de Cadenas de Aminoácidos).

Vamos a unir las letras y crear ese objeto:

# La función c2s() de seqinr une los caracteres ('M', 'T', 'E'...) en una sola palabra ("MTE...")
seq_humano <- seqinr::c2s(humano)
seq_raton <- seqinr::c2s(raton)
seq_guajolote <- seqinr::c2s(guajolote)

# Creamos el set de secuencias con sus respectivos nombres para que la tabla sea legible
secuencias_hras <- AAStringSet(c(
  "Humano" = seq_humano,
  "Raton" = seq_raton,
  "Guajolote" = seq_guajolote
))

Ahora lanzamos el algoritmo. En este caso, usaremos MUSCLE, que es excelente y muy rápido para proteínas de esta longitud.

# Ejecutar el alineamiento
alineamiento_multiple <- msa(secuencias_hras, method = "Muscle")

# Imprimir el resultado completo en la consola
print(alineamiento_multiple, show = "complete")

## 
## MsaAAMultipleAlignment with 3 rows and 189 columns
##     aln (1..54)                                            names
## [1] MTEYKLVVVGAGGVGKSALTIQLIQNHFVDEYDPTIEDSYRKQVVIDGETCLLD Humano
## [2] MTEYKLVVVGAGGVGKSALTIQLIQNHFVDEYDPTIEDSYRKQVVIDGETCLLD Raton
## [3] MTEYKLVVVGAGGVGKSALTIQLIQNHFVDEYDPTIEDSYRKQVVIDGETCLLD Guajolote
## Con MTEYKLVVVGAGGVGKSALTIQLIQNHFVDEYDPTIEDSYRKQVVIDGETCLLD Consensus 
## 
##     aln (55..108)                                          names
## [1] ILDTAGQEEYSAMRDQYMRTGEGFLCVFAINNTKSFEDIHQYREQIKRVKDSDD Humano
## [2] ILDTAGQEEYSAMRDQYMRTGEGFLCVFAINNTKSFEDIHQYREQIKRVKDSDD Raton
## [3] ILDTAGQEEYSAMRDQYMRTGEGFLCVFAINNTKSFEDIHQYREQIKRVKDSDD Guajolote
## Con ILDTAGQEEYSAMRDQYMRTGEGFLCVFAINNTKSFEDIHQYREQIKRVKDSDD Consensus 
## 
##     aln (109..162)                                         names
## [1] VPMVLVGNKCDLAARTVESRQAQDLARSYGIPYIETSAKTRQGVEDAFYTLVRE Humano
## [2] VPMVLVGNKCDLAARTVESRQAQDLARSYGIPYIETSAKTRQGVEDAFYTLVRE Raton
## [3] VPMVLVGNKCDLPARTVETRQAQDLARSYGIPYIETSAKTRQGVEDAFYTLVRE Guajolote
## Con VPMVLVGNKCDLAARTVESRQAQDLARSYGIPYIETSAKTRQGVEDAFYTLVRE Consensus 
## 
##     aln (163..189)              names
## [1] IRQHKLRKLNPPDESGPGCMSCKCVLS Humano
## [2] IRQHKLRKLNPPDESGPGCMSCKCVLS Raton
## [3] IRQHKLRKLNPPDESGPGCMNCKCVIS Guajolote
## Con IRQHKLRKLNPPDESGPGCMSCKCVLS Consensus

Para hacer el árbol, necesitamos convertir nuestro alineamiento a un formato matemático, calcular una “matriz de distancia” (que mide qué tan diferente es cada secuencia de las demás) y luego construir el árbol usando el clásico algoritmo Neighbor-Joining.

# ==========================================
# PARTE 1: ÁRBOL FILOGENÉTICO
# ==========================================

# 1. Convertir el objeto msa a un formato compatible con cálculos de distancia
aln_seqinr <- msaConvert(alineamiento_multiple, type = "seqinr::alignment")

# 2. Calcular la matriz de distancia evolutiva (basada en similitud de identidad)
matriz_distancia <- dist.alignment(aln_seqinr, "identity")

# 3. Construir el árbol usando el algoritmo Neighbor-Joining (NJ)
arbol_filogenetico <- nj(matriz_distancia)

# 4. Dibujar el árbol en la ventana de gráficos (Plots)
plot(arbol_filogenetico, 
     main = "Árbol Filogenético de la Proteína HRas", 
     edge.width = 1,     # Grosor de las ramas
     font = 1,           # Letra en negrita
     cex = 1)          # Tamaño del texto

# ==========================================
# PARTE 2: GRÁFICO DE ALINEAMIENTO (CORREGIDO)
# ==========================================

# 1. Extraemos el alineamiento como un objeto estándar
secuencias_alineadas <- as(alineamiento_multiple, "AAStringSet")

# 2. PASO CLAVE: Guardamos el alineamiento en el disco duro
writeXStringSet(secuencias_alineadas, "alineamiento_HRas_final.fasta")

# 3. Le decimos a ggmsa que lea directamente desde el archivo
grafico_color <- suppressWarnings(
  ggmsa("alineamiento_HRas_final.fasta", 
        start = 120,       
        end = 189,         
        color = "Clustal", # Esquema de colores universal y robusto
        seq_name = TRUE) + 
    ggtitle("Región Hipervariable de HRas (Aminoácidos 120-189)") +
    theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
)

# 4. Imprimir el gráfico
print(grafico_color)

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

Alineamiento

Niko Alain Cruz Sancen

2026-07-09

Configuración del Entorno y Dependencias

1. Alineamiento Global y Local de Secuencias de ADN

Ejercicio

Solución

Resultados del Alineamiento

2. Descarga de Secuencias Reales y Análisis de Homología Evolutiva

Ejercicio

Solución