Laboratorio Filogenia
Claudia Rivera
Adolfo González
Juan Giglio
Deparamento de Ingeniería Informática, Universidad de Santiago de Chile
Adolfo González
Juan Giglio
Deparamento de Ingeniería Informática, Universidad de Santiago de Chile
1 Introducción
El virus zoonótico SARS-CoV-2, responsable de la pandemia de COVID-19, emergió en Wuhan, China, en noviembre de 2019. La pandemia ha afectado a más de 107 millones de personas globalmente y ha sido fatal para más de 2.3 millones de personas. En Hawai, se han reportado más de 26,000 casos y 400 muertes. Para entender la evolución del virus, se examina su genoma, especialmente el gen S que codifica la proteína de espiga. Las interacciones entre esta proteína y el receptor humano ACE2 permiten al virus ingresar a las células. Mutaciones en el gen S pueden afectar la transmisión y la virulencia del virus. En el paper de Maison et al. (2021) se analiza un segmento de 969 pb del gen S de SARS-CoV-2 en pacientes de Hawai para comprender las adaptaciones en la proteína de espiga, crucial para el desarrollo de vacunas.
En el siguiente trabajo se replica el análisis de de alineamiento e inferencia filogenética del artículo de Maison et al. (2021) en R. En el trabajo de Maison et al. (2021) se caracterizan segmentos de 969bp del gen S del SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) correspondientes a dos muestras de pacientes en Hawaii infectados en el año 2020.
1.1 Método de Maison et al. (2021)
En el artículo de Maison et al. (2021) el proceso de alineamiento múltiple de secuencias de nucleótidos se lleva a cabo utilizando el programa MUSCLE. La construcción del árbol filogenético se realiza utilizando el software MEGAX 19.20, para obtener un árbol inicial se aplican los algoritmos Neighbor-Join y BioNJ a una matriz de distancia estimada mediante el modelo de Tamura-Nei. Se utiliza el método de Máxima Verosimilitud con 1,000 réplicas de remuestreo (bootstraps) para evaluar la robustez del árbol. Finalmente, el árbol filogenético generado por MEGAX se enraíza utilizando FigTree versión 1.4.4. La elección del punto de raíz se basa en la secuencia de referencia del coronavirus humano alfa 299E (KF514433), proporcionando un marco de referencia para la interpretación de la evolución del SARS-CoV-2.
Para construir la filogenia se utilizan 20 secuencias de coronavirus, 11 SARS-CoV-2 y 9 no SARS-CoV-2. Ver Figura 1.1.
knitr::include_graphics("hjhsw8003_0052_fig003.jpg")
Figura 1.1: Árbol filogenético, se utiliza el texto amarillo para el linaje A de betacoronavirus, rojo para el betacoronavirus del linaje B, negro para el betacoronavirus libaje C, azul para el betacoronavirus linaje D y púrpura para para el alfacoronavirus. Fuente: Maison et al. (2021)
1.2 Descripción del conjunto de datos
El conjunto de datos para la inferencia filogenética está compuesto por 20 secuencias de nucleótidos
1.3 Reproducción del análisis filogenético
1.3.1 Obtención de datos
# Librerías
library("msa")
library("BiocManager")
library("muscle")
library(ggmsa)
library(seqinr)
library("ape")
library (rentrez)
library("phangorn")
library("gridExtra")
library("ggplot2")
library("ggtree")
library(cowplot)
library(bios2mds)
library(phytools)
set.seed(0)# Se descargan los mismos conjuntos de datos, para los genomas completos se filtra el gen "S"
ids=c(
"mw237664", # <1..>972
"mw237663", # <1..>972
"mt627421", # 21563..25384
"mt407659", # 21515..25336
"mt994395", # 21550..25371
"mw035565", # 21509..25330
"mw035511", # 21509..25330
"mw064483", # 21509..25330
"nc_045512", # 21563..25384
"mt344949", # 21563..25384
"mt093571", # 21563..25384
"nc_014470", # 21391..25170
"ay274119", # 21492..25259
"ky417146", # 21493..25260
"kc881006", # 21492..25262
"nc_019843", # 21456..25517
"ef065513", # 20974..24798
"kf530099", # 23621..27706
"ay391777", # 23644..27729
"kf514433" # 20539..24051
)
start=c(
1,
1,
21563,
21515,
21550,
21509,
21509,
21509,
21563,
21563,
21563,
21391,
21492,
21493,
21492,
21456,
20974,
23621,
23644,
20539
)
end=c(
972,
972,
25384,
25336,
25371,
25330,
25330,
25330,
25384,
25384,
25384,
25170,
25259,
25260,
25262,
25517,
24798,
27706,
27729,
24051
)
# Groups
branches=list(
BetacoronavirusB=c("mw237664", # <1..>972
"mw237663" , # <1..>972
"mt627421" , # 21563..25384
"mt407659" , # 21515..25336
"mt994395" , # 21550..25371
"mw035565" , # 21509..25330
"mw035511" , # 21509..25330
"mw064483" , # 21509..25330
"nc_045512", # 21563..25384
"mt344949" , # 21563..25384
"mt093571" , # 21563..25384
"nc_014470", # 21391..25170
"ay274119" , # 21492..25259
"ky417146" , # 21493..25260
"kc881006") , # 21492..25262
BetacoronavirusC=c("nc_019843"), # 21456..25517
BetacoronavirusD=c("ef065513") , # 20974..24798
BetacoronavirusA=c("kf530099" , # 23621..27706
"ay391777") , # 23644..27729
Alfacoronavirus="kf514433" # 20539..24051
)
seqs<-rentrez::entrez_fetch(db = "nucleotide",
id = ids,
rettype = "fasta",
)
file_conn=file("covid.fa")
writeLines(seqs,file_conn)
close(file_conn)1.3.2 Alineamiento múltiple
# Al igual que en el estudio de Maison se utiliza MUSCLE, esta vez utilizando su implementación en R
seqs_vec=readDNAStringSet('covid.fa')
seqs_ss=DNAStringSet(seqs_vec,start=start,end=end)
names(seqs_ss)=ids
seqs_align<-msa(inputSeqs = seqs_ss,method="Muscle",type="dna")
saveRDS(seqs_align, "seqs_align.RDS")
msaPrettyPrint(
seqs_align,
output="pdf",
askForOverwrite=FALSE,
showLogo = "top",
showLegend=F,
consensusColor="ColdHot",
)knitr::include_graphics("seqs_align.pdf")Figura 1.2: Resultado del alineamiento múltiple. Fuente: Elaboración propia
1.3.3 Cálculo de la matriz de distancia
Se construye el árbol filogenético inicial aplicando el algoritmo bionj sobre la matriz de distancia TN93 calculada con el método de Tamura y Nei Tamura and Nei (1993). Tamura y Nei desarrollaron un método matemático para estimar el número sustituciones transicionales y tranversionales por sitio y en total, considerando exceso de transiciones, frecuencias diferentes de nucleótidos, y la variación del ratio de entre diferentes sitios Tamura and Nei (1993).
# Cálculo de distancia de Tamura Nei entre pares de secuencias de dna, el modelo de Tamura y Nei
seqs_align<-readRDS("seqs_align.RDS")
seqs_phydat=as.phyDat(seqs_align)
dist_tn93=dist.dna(
as.DNAbin( seqs_phydat),
model="TN93"
)
heatmap(as.matrix(dist_tn93))
Figura 1.3: Matriz de distancia
1.3.4 Construcción del árbol filogenético
nj_tree=bionj(dist_tn93)
# Quitar multifurcaciones
nj_tree = multi2di(nj_tree)
# Enraizar árbol en
nj_tree =root(nj_tree,"kf514433")
# No se utiliza TN93 porque genera valores faltantes
f = function(x) bionj(
dist.dna(
x,
#model="TN93"
))
bootstrapped_trees = boot.phylo(nj_tree, as.DNAbin(seqs_phydat), f, trees = T,B=1000)##
Running bootstraps: 100 / 1000
Running bootstraps: 200 / 1000
Running bootstraps: 300 / 1000
Running bootstraps: 400 / 1000
Running bootstraps: 500 / 1000
Running bootstraps: 600 / 1000
Running bootstraps: 700 / 1000
Running bootstraps: 800 / 1000
Running bootstraps: 900 / 1000
Running bootstraps: 1000 / 1000
## Calculating bootstrap values... done.# Add bootstrap confidence values to the nodes of the original tree
bootstrapped_conf_tree = addConfidences(nj_tree, bootstrapped_trees$trees)
bootstrapped_conf_tree = groupOTU(bootstrapped_conf_tree,branches)
# Define a color palette for visualizing bootstrap support
color_palette = colorRampPalette(c("red3", "mediumseagreen"))(100)
color_positions = round(bootstrapped_conf_tree$node.label * 97, 0)
colors = color_palette[color_positions]
# Create a plot of the bootstrapped tree using ggtree
g = ggtree(bootstrapped_conf_tree, layout = "dendrogram",branch.length = 'none')
g = g + geom_tiplab(size = 4,aes(color=group)) # Add labels to terminal nodes
g = g + theme(plot.margin = margin(t = 0, r = 0, b = 150, l = 0, unit = "pt"), legend.position = "right")
g = g + geom_nodepoint(size = 4, color = colors) + theme(legend.position='top')
g = g + labs(title = "BioNJ tree Bootstrapping (1000 trees)")
#g=g+ geom_text(aes(color=place, label=value), hjust=1, vjust=1.4, size=3)
# Create a data frame with values from 1 to 100 for color mapping
color_data = data.frame(x = 1:100, y = 1)
color_data$colors = color_palette
# Create a color bar plot using ggplot2
g1 = ggplot(color_data, aes(x = x, y = y, fill = colors)) +
geom_bar(stat = "identity", width = 1, position = "identity") + scale_fill_manual(values = rev(color_data$colors)) +
labs(title = NULL,
x = "% of node support",
y = NULL) +
theme_minimal() + theme(axis.text.y = element_blank(), # Remove y-axis labels
axis.title.y = element_blank(),
panel.grid = element_blank()) # Remove y-axis title
g1 = g1 + guides(fill = "none")
# Arrange the tree plot and the color bar plot
grid.arrange(g, g1, nrow = 2, ncol = 1, heights = c(0.8, 0.2))
Figura 1.4: Árbol filogenético.
# Gráfico con alineamiento entre pos 1500 y 3900
seqs_bios2mds=msaConvert(seqs_align, "bios2mds::align")
seqs_phydat=as.phyDat(seqs_align)
export.fasta(seqs_bios2mds)
# https://www.biostars.org/p/132909/
msaplot(
p=g,#ggtree(bootstrapped_conf_tree,branch.length = 'none'),
fasta="alignment.fa",
window=c(1500, 3900),
width=5,
height=0.8,
offset = 50
)
Figura 1.5: Árbol filogenético con secuencia entre las posiciones 1500 y 3900.
# Se construye el árbol filogenético de Maison et al.
tmaison = read.tree(text="((((((((((mw237664:1,mw237663:1):1,(mt627421:1,mt407659:1):1):1,mt994395:1):1):1,(mw035565:1,mw035511:1):1):1,mw064483:1):1):1):1,((nc_045512:1,(mt344949:1,mt093571:1):1):1,(nc_014470:1,(ay274119:1,(ky417146:1,kc881006:1):1):1):1):1):1,(nc_019843:1,(ef065513:1,(kf530099:1,ay391777:1):1):1):1,kf514433:1);")
tmaison= groupOTU(tmaison,branches)
tmaison=root(tmaison,"kf514433")
tpropuesto = bootstrapped_conf_tree # Comparación
comparePhylo2(
tmaison,
tpropuesto,
plot = T,
use.edge.length = FALSE,
#commons = TRUE,
location = "bottomright",
type="phylogram",
)
Figura 1.6: Ramas compartidas entre los árboles filogenéticos.
## => Comparing tmaison with tpropuesto.
## Both trees have the same number of tips: 20.
## Both trees have the same tip labels.
## Trees have different numbers of nodes: 21 and 18.
## Both trees are unrooted.
## Both trees are not ultrametric.
## 5 splits in common.tangleTree<-cophylo(tmaison,tpropuesto,rotate=FALSE,use.edge.length=F)
plot(tangleTree)
Figura 1.7: Comparación entre hojas de los árboles filogenéticos.
1.3.5 Discusión
La interpretación del árbol filogenético en el marco de la investigación sobre las características genéticas y la filogenia de la secuencia del gen S de SARS-CoV-2 en Hawái, con un énfasis en la mutación P681H, implica el examen de la disposición de las ramas y nodos en el árbol. Las ramas cercanas indican una mayor similitud genética, mientras que los nodos representan antepasados comunes y las ramas divergentes sugieren eventos evolutivos. Es crucial localizar las muestras de Hawái en el árbol para entender cómo se relacionan genéticamente con cepas de otras regiones. La longitud de las ramas puede ofrecer información sobre la antigüedad de eventos evolutivos, y es necesario analizar la presencia y distribución de la mutación P681H en las diferentes ramas del árbol. Al replicar el experimento, se ha observado que, en general, las relaciones evolutivas se mantienen con un mayor soporte, de igual forma en todo el árbol filogenético. Para el árbol de la Figura 1.4, se aplica la librería ggtree para la representación visual de un árbol filogenético derivado de bootstrapping, con una atención especial a la personalización estética y la inclusión de información específica. Este enfoque se centra en la presentación visual de un árbol bootstrap generado a partir de 1000 repeticiones, ofreciendo insights sobre la robustez de las ramificaciones. En dicho árbol se puede observar una diferencia al obtenido en el experimento inicial, esta discrepancia podría deberse a varios factores. El bootstrapping es un método estadístico que implica la generación de múltiples árboles a partir de subconjuntos aleatorios de los datos. Cada árbol bootstrap representa una variante plausible de la estructura filogenética dada la variabilidad inherente en los datos. En primer lugar, la variabilidad en las técnicas de muestreo podría influir, especialmente si el muestreo de secuencias no se llevó a cabo de manera representativa en términos de regiones geográficas o momentos temporales. Además, los errores en la alineación de secuencias, como asignaciones incorrectas de bases o la omisión de regiones críticas, podrían afectar la precisión del árbol. La elección de un modelo de evolución molecular inadecuado también puede jugar un papel crucial, ya que diferentes regiones genéticas pueden evolucionar bajo modelos distintos. La longitud de la secuencia utilizada en el análisis también es un factor importante; secuencias más cortas pueden no contener suficiente información para resolver algunas relaciones evolutivas. Problemas de convergencia del algoritmo, especialmente en métodos de máxima verosimilitud o bayesianos, podrían requerir ajustes de parámetros para obtener resultados más precisos. Las ramas del árbol con bajo soporte podrían indicar áreas de incertidumbre, y estas regiones podrían ser el resultado de datos ruidosos, variabilidad genética real o insuficiente información. Además, la variabilidad natural en las secuencias genéticas o eventos evolutivos recientes podría introducir complejidades en la biología evolutiva real y afectar la construcción precisa del árbol. Es fundamental revisar cada etapa del análisis, desde la adquisición de datos hasta la interpretación de resultados, ajustando parámetros según sea necesario y comparando con estudios similares para abordar cualquier problema potencial. Además, la utilización de versiones actualizadas del software y la adhesión a las mejores prácticas son cruciales para garantizar la robustez y fiabilidad de los resultados. La variabilidad en la reproducción de árboles filogenéticos entre programas como MEGAX 19.20 y R puede deberse a diversas razones fundamentales. En primer lugar, la disparidad en los algoritmos y métodos utilizados para construir estos árboles puede conducir a resultados divergentes. Cada programa puede emplear enfoques diferentes, como métodos basados en distancias, máxima parsimonia o máxima verosimilitud, lo que influye en la topología final del árbol. Otra fuente potencial de diferencias radica en la configuración de parámetros. Los valores predeterminados y la interpretación de ciertos parámetros pueden variar entre programas, lo que afecta directamente a la construcción del árbol filogenético. Además, los modelos de evolución molecular utilizados para describir la dinámica evolutiva de las secuencias biológicas pueden diferir, introduciendo variaciones en los resultados. La forma en que se manejan y procesan los datos de secuencias también puede contribuir a discrepancias. Desde la limpieza de datos hasta la alineación de secuencias y el tratamiento de gaps, las variaciones en el proceso de manipulación de datos pueden influir en la generación del árbol. Además, las diferencias en la implementación del software, incluyendo aspectos técnicos como la precisión numérica y la gestión de memoria, pueden afectar los resultados finales. Es importante considerar también la posibilidad de actualizaciones de software, ya que distintas versiones pueden introducir cambios en los algoritmos o en la forma en que se procesan los datos. Por último, la interpretación y visualización de resultados puede variar entre programas, con opciones de visualización diferentes que influyen en la representación gráfica del árbol filogenético. En resumen, la variabilidad entre árboles filogenéticos generados en diferentes programas es multifacética y puede ser resultado de diversas influencias desde los métodos algorítmicos hasta la manipulación de datos y la interpretación de resultados.
1.3.6 Conclusiones
La investigación genética del virus SARS-CoV-2 en Hawái se orienta a comprender la evolución y las características genéticas del virus en esta región. En este sentido, el análisis de un árbol filogenético permite analizar la relación genética entre diferentes cepas del virus y la mutación P681H, arrojando luz sobre su evolución y dispersión en Hawái y otras regiones. En este sentido, la disposición de las ramas y nodos en el árbol filogenético entrega información acerca de las similitudes genéticas, antepasados comunes y eventos evolutivos. Se puede inferir que la investigación genética enfrenta desafíos, como la variabilidad en las técnicas de muestreo, errores en la alineación de secuencias, la elección de modelos de evolución molecular adecuados, la longitud de la secuencia utilizada y problemas de convergencia del algoritmo. La variabilidad en la reproducción de árboles filogenéticos entre diferentes programas de software puede deberse a diferencias en los algoritmos que utilizan, configuración de parámetros, modelos de evolución molecular, manejo de datos de secuencias, implementación de software y actualizaciones de software. En este sentido se observa que, para obtener resultados más precisos en la investigación genética y la interpretación de árboles filogenéticos, es esencial revisar cuidadosamente cada etapa del análisis, ajustar parámetros según sea necesario, comparar con estudios similares y seguir las mejores prácticas en el campo. La reproducción del experimento revela en términos generales que las relaciones evolutivas se mantienen con un soporte más robusto; sin embargo, se observa una variación en la estabilidad de todo el árbol filogenético. Este fenómeno sugiere que ciertas ramas del árbol son más resilientes, posiblemente menos susceptibles a ruido o variaciones en los datos, mientras que otras regiones pueden ser más sensibles a factores específicos. La consistencia en las relaciones evolutivas con mayor soporte respalda la confianza en la estructura filogenética general de esas áreas. Sin embargo, la discrepancia en la estabilidad del árbol también destaca la necesidad de realizar ajustes adicionales o validaciones específicas, tales como la revisión de parámetros, la optimización de modelos de evolución y la evaluación de la calidad de las secuencias utilizadas. La variabilidad observada podría indicar la presencia de eventos evolutivos recientes o puntos de divergencia genética que requieren una investigación más detallada. A pesar de estas variaciones, la capacidad de mantener relaciones evolutivas sólidas en gran parte del árbol sugiere que la estructura filogenética general sigue siendo informativa para comprender la evolución de SARS-CoV-2. Sin embargo, es crucial abordar estas variaciones con un enfoque crítico, considerando cuidadosamente las limitaciones del estudio y reconociendo las áreas de incertidumbre como oportunidades para futuras investigaciones y refinamientos metodológicos.