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
## Descripción del conjunto de datos
El conjunto de datos para la inferencia filogenética está compuesto por 20 secuencias de nucleótidos
1.2 Reproducción del análisis filogenético
1.2.1 Obtención de datos
# 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':
##
## strsplitlibrary("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':
##
## versionlibrary("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/bbac222library(seqinr)##
## Attaching package: 'seqinr'## The following object is masked from 'package:Biostrings':
##
## translatelibrary("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':
##
## complementlibrary (rentrez)
library("phangorn")
library("gridExtra")##
## Attaching package: 'gridExtra'## The following object is masked from 'package:BiocGenerics':
##
## combinelibrary("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
##
## Guangchuang Yu. Data Integration, Manipulation and Visualization of
## Phylogenetic Trees (1st edition). Chapman and Hall/CRC. 2022,
## doi:10.1201/9781003279242
##
## S Xu, Z Dai, P Guo, X Fu, S Liu, L Zhou, W Tang, T Feng, M Chen, L
## Zhan, T Wu, E Hu, Y Jiang, X Bo, G Yu. ggtreeExtra: Compact
## visualization of richly annotated phylogenetic data. Molecular Biology
## and Evolution. 2021, 38(9):4039-4042. doi: 10.1093/molbev/msab166##
## 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# 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
)
seqs<-rentrez::entrez_fetch(db = "nucleotide",
id = ids,
rettype = "fasta",
)
file_conn=file("covid.fa")
writeLines(seqs,file_conn)
close(file_conn)1.2.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",
#y=c(1,),
#output='pdf',
askForOverwrite=FALSE,
showLogo = "top",
showLegend=F,
consensusColor="ColdHot",
)1.2.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))
1.2.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")
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)
# 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)## ! 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.g = g + geom_tiplab(size = 4, color = "black") # 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 = "NJ tree Bootstrapping (1000 trees)")
# 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))
1.2.5 Discusión
Se observa que en general se mantienen las relaciones evolutivas con mayor soporte.