Este documento guía a través de los pasos para alinear secuencias y construir redes de haplotipos utilizando los paquetes msa y haplotypes, con apoyo de ape, Biostrings, seqinr y pegas.

⚠️ Recuerde establecer su carpeta de trabajo antes de iniciar. Esta acción también puede ser definida por medio de las pestañas en la parte superior Session > Set Working Directory > Choose Directory.

Instalación y llamado de paquetes

Instalación de paquetes

#install.packages("BiocManager")
#BiocManager::install("msa")
#BiocManager::install("Biostrings")

#Nota: Si está ejecutando el código por primera vez o no ha usado las bibliotecas mencionadas anteriormente, instalar las bibliotecas eliminando el símbolo de anotación "#" al inicio de cada comando. 
#install.packages("ape")
#install.packages("haplotypes")
#install.packages("seqinr")
#install.packages("pegas")

#Nota: Si está ejecutando el código por primera vez o no ha usado las bibliotecas mencionadas anteriormente, instalar las bibliotecas eliminando el símbolo de anotación "#" al inicio de cada comando.

Para cargar los paquetes

#Notar que algunas bibliotecas también se usan para reconstruir árboles filogenéticos. 
library(Biostrings)
library(ape)         # Distancias genéticas y NJ
library(msa)         # Alineamiento múltiple
library(haplotypes)  # Manipulación secuencias genéticas y redes de haplotipo
library(seqinr)      # Análisis exploratorio y visualización de secuencias biológicas
library(pegas)       # Agrupamiento en árboles

Cargar y alinear las secuencias

Para leer el archivo .fasta (aún no alineado)

secuencias_lobster <- readDNAStringSet("lobster.fasta", format = "fasta")

#Visualizar nuestras secuencias
secuencias_lobster
## DNAStringSet object of length 187:
##       width seq                                             names               
##   [1]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR01_Pg025
##   [2]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR02_Pg026
##   [3]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR03_Pg027
##   [4]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR03_Pg028
##   [5]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR04_Pg029
##   ...   ... ...
## [183]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR03_Pg106
## [184]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR04_Pg107
## [185]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR02_Pg151
## [186]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR02_Pg167
## [187]   400 TCCTAAGCAAGAATCAAACTTC...TCTACACTTAAAATATTAACCA AR01_Pg187

Para alinear con MUSCLE (puede cambiarse por “ClustalW”)

aligned_lobster <- msa(secuencias_lobster, method = "Muscle")

aligned_lobster
## MUSCLE 3.8.31   
## 
## Call:
##    msa(secuencias_lobster, method = "Muscle")
## 
## MsaDNAMultipleAlignment with 187 rows and 400 columns
##       aln                                                  names
##   [1] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR03_Pg173
##   [2] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR03_Pg174
##   [3] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR03_Pg137
##   [4] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR01_Pg115
##   [5] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR02_Pg130
##   [6] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR02_Pg131
##   [7] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR02_Pg132
##   [8] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR03_Pg133
##   [9] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR04_Pg138 
##   ... ...
## [180] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR01_Pg038
## [181] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR01_Pg037
## [182] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR01_Pg036
## [183] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR01_Pg035
## [184] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR01_Pg034
## [185] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR01_Pg033
## [186] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR01_Pg032
## [187] TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA AR01_Pg031
##   Con TCCTAAGCAAGAATCAAACTTCTTC...TGTCTACACTTAAAATATTAACCA Consensus
aligned_lobster.DNAbin <- as.DNAbin(aligned_lobster)
#¿Qué clase de objeto tenemos aquí?
class(aligned_lobster.DNAbin)
## [1] "DNAbin"

PASO OPCIONAL: Exportar el alineamiento a su computador

Si deseamos exportar este alineamiento en archivo .fasta, para su visualización o edición en otro software.

aligned_lobster.DNAStringSet <- DNAStringSet(aligned_lobster)

writeXStringSet(aligned_lobster.DNAStringSet, filepath = "lobster_alineado.fasta")

Luego, puede cargar su archivo de alientamiento desde el computador, en caso deba cortar y editar flancos.

lobster.aligned <- read.dna("lobster_alineado.fasta", format = "fasta")
#Nota: El archivo cargada se actualizará a formato DNAbin, compatible con los siguientes pasos.

Estimar diversidad haplotípica y nucleotídica

Obteniendo valores de diversidad haplotípica con comando hap.div() de pegas

hap.div(aligned_lobster.DNAbin, variance = TRUE)
## [1] 0.8579150135 0.0004894735

Obteniendo valores de diversidad nucleotídica con comando nuc.div() de pegas

nuc.div(aligned_lobster.DNAbin, variance = TRUE)
## [1] 4.124941e-03 7.296537e-06

Construcción de red de haplotipos con haplotypes

Para obtener haplotipos

lobster_haps <- haplotype(aligned_lobster.DNAbin)

lobster_haps
## 
## Haplotypes extracted from: aligned_lobster.DNAbin 
## 
##     Number of haplotypes: 57 
##          Sequence length: 400 
## 
## Haplotype labels and frequencies:
## 
##       I      II     III      IV       V      VI     VII    VIII      IX       X 
##       2       1       1       4       1       1       1       5       1       1 
##      XI     XII    XIII     XIV      XV     XVI    XVII   XVIII     XIX      XX 
##       1       1       1       1       1       2       1       1       1       1 
##     XXI    XXII   XXIII    XXIV     XXV    XXVI   XXVII  XXVIII    XXIX     XXX 
##       1       8       1       2       1       6       1       7       2       1 
##    XXXI   XXXII  XXXIII   XXXIV    XXXV   XXXVI  XXXVII XXXVIII   XXXIX      XL 
##       1       1       1       1       1      23       4       1       2       1 
## ...
## (use summary() to print all)
summary(lobster_haps)
##       I      II     III      IV       V      VI     VII    VIII      IX       X 
##       2       1       1       4       1       1       1       5       1       1 
##      XI     XII    XIII     XIV      XV     XVI    XVII   XVIII     XIX      XX 
##       1       1       1       1       1       2       1       1       1       1 
##     XXI    XXII   XXIII    XXIV     XXV    XXVI   XXVII  XXVIII    XXIX     XXX 
##       1       8       1       2       1       6       1       7       2       1 
##    XXXI   XXXII  XXXIII   XXXIV    XXXV   XXXVI  XXXVII XXXVIII   XXXIX      XL 
##       1       1       1       1       1      23       4       1       2       1 
##     XLI    XLII   XLIII    XLIV     XLV    XLVI   XLVII  XLVIII    XLIX       L 
##       1       2       1       1       2       1       1       1       1       1 
##      LI     LII    LIII     LIV      LV     LVI    LVII 
##       1       1       1       1      10       1      65

Para construir la red de haplotipos con el método TCS

lobster_net <- haploNet(lobster_haps)

#Visualizamos
lobster_net
## Haplotype network with:
##   57 haplotypes
##   1210 links
##   link lengths between 1 and 3 step(s)
## 
## Use print.default() to display all elements.

Para obtener la frecuencia de cada haplotipo

hap_freq <- sapply(attr(lobster_haps, "index"), length)

#Visualizamos
hap_freq
##  [1]  2  1  1  4  1  1  1  5  1  1  1  1  1  1  1  2  1  1  1  1  1  8  1  2  1
## [26]  6  1  7  2  1  1  1  1  1  1 23  4  1  2  1  1  2  1  1  2  1  1  1  1  1
## [51]  1  1  1  1 10  1 65

Para graficar la red de haplotipos

plot(lobster_net,
     size = hap_freq,
     scale.ratio = 10,  # más espacio entre nodos
     labels = TRUE,
     show.mutation = FALSE,
     col = "skyblue",
     cex = 0.5,   # tamaño del texto
     font = 1)    # estilo del texto

¿Cómo podríamos mejorar el gráfico resultante? Por ejemplo, cambiar la etiqueta del haplotipo o colorearlo de acuerdo a su distribución geográfica.

OPCIONAL: Podemos exportar la red de haplotipos en SVG o PDF para editarla en otro programa. Para ello usaremos los comandos pdf() o svg()

svg("red_haplotipos.svg", width = 10, height = 10)
plot(lobster_net,
     size = hap_freq,
     scale.ratio = 10,
     labels = TRUE,
     show.mutation = FALSE,
     col = "skyblue",
     cex = 0.5,
     font = 1)
dev.off()

Aclaración: Cuando el objetivo es estudiar lo que ocurre dentro de una población —por ejemplo, en análisis de diversidad genética o en la reconstrucción de redes de haplotipos—, no es recomendable incluir un grupo externo (como una especie distinta). Esto se debe a que las diferencias genéticas acumuladas entre la población de estudio y el grupo externo son generalmente muy amplias, lo que provocaría una separación artificial en dos grupos bien definidos. Este efecto enmascara la variación interna de la población, que es precisamente lo que se busca analizar. Por ello, en estudios de estructura o relaciones internas, el uso de grupos externos puede distorsionar la interpretación biológica de los resultados.