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.
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
## Descripción del conjunto de datos
El conjunto de datos para la inferencia filogenética está compuesto por 20 secuencias de nucleótidos
# Librerías
library("msa")
## Loading required package: Biostrings
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## Loading required package: stats4
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:utils':
##
## findMatches
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: XVector
## Loading required package: GenomeInfoDb
##
## Attaching package: 'Biostrings'
## The following object is masked from 'package:base':
##
## strsplit
library("BiocManager")
## Bioconductor version '3.17' is out-of-date; the current release version '3.18'
## is available with R version '4.3'; see https://bioconductor.org/install
##
## Attaching package: 'BiocManager'
## The following object is masked from 'package:msa':
##
## version
library("muscle")
library(ggmsa)
## Registered S3 methods overwritten by 'ggalt':
## method from
## grid.draw.absoluteGrob ggplot2
## grobHeight.absoluteGrob ggplot2
## grobWidth.absoluteGrob ggplot2
## grobX.absoluteGrob ggplot2
## grobY.absoluteGrob ggplot2
## Registered S3 methods overwritten by 'treeio':
## method from
## MRCA.phylo tidytree
## MRCA.treedata tidytree
## Nnode.treedata tidytree
## Ntip.treedata tidytree
## ancestor.phylo tidytree
## ancestor.treedata tidytree
## child.phylo tidytree
## child.treedata tidytree
## full_join.phylo tidytree
## full_join.treedata tidytree
## groupClade.phylo tidytree
## groupClade.treedata tidytree
## groupOTU.phylo tidytree
## groupOTU.treedata tidytree
## inner_join.phylo tidytree
## inner_join.treedata tidytree
## is.rooted.treedata tidytree
## nodeid.phylo tidytree
## nodeid.treedata tidytree
## nodelab.phylo tidytree
## nodelab.treedata tidytree
## offspring.phylo tidytree
## offspring.treedata tidytree
## parent.phylo tidytree
## parent.treedata tidytree
## root.treedata tidytree
## rootnode.phylo tidytree
## sibling.phylo tidytree
## ggmsa v1.6.0 Document: http://yulab-smu.top/ggmsa/
##
## If you use ggmsa in published research, please cite:
## L Zhou, T Feng, S Xu, F Gao, TT Lam, Q Wang, T Wu, H Huang, L Zhan, L Li, Y Guan, Z Dai*, G Yu* ggmsa: a visual exploration tool for multiple sequence alignment and associated data. Briefings in Bioinformatics. DOI:10.1093/bib/bbac222
library(seqinr)
##
## Attaching package: 'seqinr'
## The following object is masked from 'package:Biostrings':
##
## translate
library("ape")
##
## Attaching package: 'ape'
## The following objects are masked from 'package:seqinr':
##
## as.alignment, consensus
## The following object is masked from 'package:muscle':
##
## muscle
## The following object is masked from 'package:Biostrings':
##
## complement
library (rentrez)
library("phangorn")
library("gridExtra")
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:BiocGenerics':
##
## combine
library("ggplot2")
library("ggtree")
## ggtree v3.8.2 For help: https://yulab-smu.top/treedata-book/
##
## If you use the ggtree package suite in published research, please cite
## the appropriate paper(s):
##
## Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam.
## ggtree: an R package for visualization and annotation of phylogenetic
## trees with their covariates and other associated data. Methods in
## Ecology and Evolution. 2017, 8(1):28-36. doi:10.1111/2041-210X.12628
##
## G Yu. Data Integration, Manipulation and Visualization of Phylogenetic
## Trees (1st ed.). Chapman and Hall/CRC. 2022. ISBN: 9781032233574
##
## Guangchuang Yu. Data Integration, Manipulation and Visualization of
## Phylogenetic Trees (1st edition). Chapman and Hall/CRC. 2022,
## doi:10.1201/9781003279242
##
## Attaching package: 'ggtree'
## The following object is masked from 'package:ape':
##
## rotate
## The following object is masked from 'package:Biostrings':
##
## collapse
## The following object is masked from 'package:IRanges':
##
## collapse
## The following object is masked from 'package:S4Vectors':
##
## expand
library(cowplot)
library(bios2mds)
## Loading required package: amap
## Loading required package: e1071
## Loading required package: scales
## Loading required package: cluster
## Loading required package: rgl
## This build of rgl does not include OpenGL functions. Use
## rglwidget() to display results, e.g. via options(rgl.printRglwidget = TRUE).
# 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)
# 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",
)
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))
nj_tree=bionj(dist_tn93)
# Quitar multifurcaciones
nj_tree = multi2di(nj_tree)
# Enraizar árbol en
nj_tree =root(nj_tree,"kf514433")
f = function(x) bionj(dist.dna(x))
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", ladderize = T,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 = "top")
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))
# 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 = 15
)
# Comparación 1
t1 = 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);")
t1= groupOTU(t1,branches)
t1=root(t1,"kf514433")
t2=reorder(t1,order="postorder")
t1 = bootstrapped_conf_tree
T1 <- ggtree(t1 ) +
theme_tree2(legend.position='none', plot.margin = unit(c(0,0,0,0),"cm")) +
geom_tiplab() + xlim(0,5.2)
## ! The tree contained negative edge lengths. If you want to ignore the edges,
## you can set `options(ignore.negative.edge=TRUE)`, then re-run ggtree.
T2 <- ggtree(t2, branch.length = 'none') +
theme_tree2(legend.position='none', plot.margin = unit(c(0,0,0,0),"cm")) +
geom_tiplab(hjust =1) + xlim(0,5.2)
d1 = T1$data[T1$data$isTip,]
d1$x[] = 1
d2 = T2$data[T2$data$isTip,]
d2$x[] = 2
TTcon = rbind(d1, d2)
T1 = ggtree(t1,branch.length = 'none' ) +
theme_tree2(
legend.position='none',
plot.margin = unit(
c(0,0,0,0),
"cm")
) +
geom_tiplab(
size = 4,
aes(color=group),
align=T
)+
xlim(0,10.5)
T2 = ggtree(t2,branch.length = 'none') +
theme_tree2(
legend.position='none',
) +
geom_tiplab(
size = 4,
aes(color=group),
offset = -3.6
)+
scale_x_reverse( limits = c(14, 0) )# if required add within scale_x_reverse() -> limits = c(5.2, 0))
L1 = ggplot(TTcon, aes(x = x, y = y, colour = label, group = label)) + geom_line() +
theme_void() + theme(legend.position="none", plot.margin = unit(c(1,0,1,0),"cm"))
cowplot::plot_grid(T1, L1 ,T2, nrow = 1, align = "hv")
# Comparación 2
comparePhylo(t1, t2, plot = T, force.rooted = FALSE,
use.edge.length = FALSE, commons = TRUE,
location = "bottomleft",type="cladogram")
## => Comparing t1 with t2.
## Both trees have the same number of tips: 20.
## Both trees have the same tip labels.
## Trees have different numbers of nodes: 18 and 21.
## Both trees are unrooted.
## Both trees are not ultrametric.
## 6 splits in common.
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 examinar 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, pero no sucede de la misma manera en todo el árbol filogenético. Esta discrepancia podría deberse a varios factores. 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.
La reproducción del experimento revela que, en términos generales, 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.