Este texto documenta el proceso asociado al estudio de la historia evolutiva de los datos de salmones provenientes del Laboratory of Genomics, Molecular Ecology and Evolutionary Studies, Universidad de Santiago de Chile.

1. Lectura de datos

Para leer las secuencias se listará cada secuencia de DNA que está en la carpeta Data, para luego abrirle usando la función read.fasta() de la biblioteca seqinr. Se usarán funciones que operan sobre cadenas de texto, biblioteca stringr, con el fin de estandarizar y reducir el nombre de las secuencias, dejándolas en el siguiente formato: secuencias$ALKBH5_Oncorhynchus_kisutch_LOC109886893.

#Bibliotecas
library("seqinr")
library("kmer")
library("phangorn")
library("stringr")
library("phytools")
library("treespace")
library("igraph")
library("ggplot2")
library("factoextra")
library("phangorn")

#--------------------------
#Apertura de archivos
#--------------------------
#Hacer listado de archivos en carpeta
archivos = list.files(pattern = "\\.fna$",recursive = T)

#Abrir archivo por archivo generando etiquetas para los nombres
secuencias=NULL
for (a in 1:length(archivos))
{
  secuencias=c(secuencias,read.fasta(archivos[a],whole.header = T))
  localizacion=strsplit(names(secuencias[a])," ")[[1]][2]
  especie=str_replace(strsplit(strsplit(names(secuencias[a]),"organism=")[[1]][2],"]")[[1]][1]," ","_")
  gen=strsplit(archivos[a],"/")[[1]][2]
  copia=""
  if (grepl("like",archivos[[a]],fixed=F)) {copia="like"}
  nuevo_nombre=paste(gen,copia,especie,localizacion,sep="_")
  names(secuencias)[a]= nuevo_nombre
}

2. Árbol de referencia

Para crear el árbol filogenético de referencia asociado a las secuencias de DNA enviadas por el equipo especialista, se usará la estrategia kmers (Gamage et al. 2020). Con ello se evita sesgar las características de cada gen al efectuar alineamiento múltiple de las mismas, creando directamente una matriz de distancia sobre las secuencias no alineadas al usar la función kdistance() de la biblioteca kmer. El parámetro \(k\) se variará desde 1 a un máximo de 8, limitado por memoria, comparándole mediante distancia Euclidiana.

El árbol respectivo se generará usando el método de Neighbour joining (Saitou and Nei 1987).

#-----------------------------
# Creación de árbol kmers
#-----------------------------

#Creación de árboles usando kmers
arboles=NULL
set.seed(1)
for (a in 1:8)
{
  #Kmers k:1-8
  distancia = kdistance(secuencias, k = a,method = "euclidean")
  
  #Neighbour joining
  arbol=NJ(distancia)
  
  #Eliminación de aristas negativas (si aparecen)
  arbol$edge.length[which(arbol$edge.length<0)]=0
  
  #Eliminación de multifurcaciones
  arbol=multi2di(arbol)
  
  arbol=midpoint(arbol)

  #Almacenamiento de árboles
  arboles[[a]]=arbol
}

#Resultados en estructura multiPhylo
class(arboles)="multiPhylo"

Posteriormente se calculará la distancia Robinson-Foulds (Briand et al. 2020) entre árboles usando la función multiDist() de la biblioteca treespace, con el fin de encontrar el árbol representativo (medoide). La matriz de distancia respectiva con y sin normalizar su amplitud, con respecto a la cantidad máxima de operaciones topológicas puede ser representada como sigue:

#-----------------------------
# Búsqueda de medoide
#-----------------------------

#Distancia sin normalizar
distancia_RF=as.matrix(multiDist(arboles))
print(ncol(distancia_RF))
## [1] 8
print(nrow(distancia_RF))
## [1] 8
for (a in 1:length(arboles))
{ 
  {for (b in 1:length(arboles))
  {distancia_RF[a,b]=RF.dist(arboles[[a]], arboles[[b]], normalize = F) }
  }
}


#Gráfico
colnames(distancia_RF)=paste("tree",colnames(distancia_RF),sep="_")
rownames(distancia_RF)=paste("tree",rownames(distancia_RF),sep="_")
fviz_dist(dist.obj = as.dist(distancia_RF), order = F, lab_size = 12,gradient = list(low = "white", mid = "gray", high = "black"))

#Distancia normalizada
distancia_RF_norm=distancia_RF

for (a in 1:length(arboles))
{ 
  {for (b in 1:length(arboles))
  {distancia_RF_norm[a,b]=RF.dist(arboles[[a]], arboles[[b]], normalize = T) }
  }
}

#Gráfico
fviz_dist(dist.obj = as.dist(distancia_RF_norm), order = F, lab_size = 12,gradient = list(low = "white", mid = "gray", high = "black"))

#Medoide
ind_medoide=which.min(colSums(distancia_RF)/nrow(distancia_RF))
print(paste("El árbol medoide es el número",ind_medoide,sep=" "))
## [1] "El árbol medoide es el número 6"
#Guardado de árboles en formato newick
write.tree(arboles,"arboles_referencia.newick")

Finalmente, los árboles representativos se almacenarán en formato newick en caso de que se requiera liberar recursos del procesador, siendo el árbol númerp 6 el medoide.

La siguiente figura corresponde a una previsualización del árbol medoide. Para una mejor visualización se recomienda usar las aplicaciones Figtree o MEGA.

plot(arboles[[ind_medoide]],cex=0.6,type = "phylogram",use.edge.length = T)
edgelabels(round(arboles[[ind_medoide]]$edge.length,4), col = "black",frame = "none", cex=0.6)

3. Árbol MSA

Con el fin de comprobar la idoneidad del árbol de referencia obtenido en la sección anterior respecto a una metodología incorrecta, como podría ser el uso de alineamiento múltiple sobre todas las secuencias (MSA, por su sigla en inglés), se inferirá el árbol respectivo. Para ello se guardarán las secuencias en un archivo, transformándolas desde el formato SeqFastadna al DNAStringSet usando funciones de la biblioteca Biostrings (código comentado). Posteriormente este archivo se abrió en MEGA, realizando el correspondiente MSA, almacenándolo en el archivo secuencia_alineada.fasta.

El árbol asociado al MSA se inferirá usando el método de Neighbour joining.

#-----------------------------
# Inferencia de árbol MSA
#-----------------------------

#library(Biostrings)
#library(seqinr)
#FUN = function(x)
#    paste(getSequence(x), collapse = "")
#sec_fasta=as(vapply(secuencias, FUN, character(1)), "DNAStringSet")
#writeXStringSet(sec_fasta, 'secuencia.fasta')

#Lectura de secuencias alineadas
secuencia_alineada=read.phyDat(file = "secuencia_alineada.fasta",type = "DNA",format = "fasta")

#Cálculo de árboles 
distancia_Hamming=dist.hamming(secuencia_alineada)
arbol_msa=NJ(as.matrix(distancia_Hamming))

#Eliminación de aristas negativas (si aparecen)
arbol_msa$edge.length[which(arbol_msa$edge.length<0)]=0

#Puntos
arbol=midpoint(arbol)

#Eliminación de multifurcaciones
arbol_msa=multi2di(arbol_msa)

Como se puede apreciar, la mayoría de las hojas del árbol teóricamente relacionadas por los genes estudiados no quedan asociadas, validando al árbol de referencia basado en kmers.

#-----------------------------
# Comparación árbol kmers vs MSA
#-----------------------------

#Graficar árbol
plot(arbol_msa,cex=0.6,type = "phylogram",use.edge.length =T)
edgelabels(round(arboles[[ind_medoide]]$edge.length,4), col = "black",frame = "none", cex=0.6)

#Distancia y homologación de nombres
arbol1=arboles[[ind_medoide]]
arbol1$tip.label=str_replace_all(arbol1$tip.label, fixed("_"), "")
arbol2=arbol_msa
arbol2$tip.label=str_replace_all(arbol2$tip.label, fixed(" "), "")
distancia_nnorm=RF.dist(arbol1,arbol2,normalize = F)
distancia_norm=RF.dist(arbol1,arbol2,normalize = T)

print(paste("La distancia no normalizada entre ambos árboles es",distancia_nnorm,"operaciones",sep=" "))
## [1] "La distancia no normalizada entre ambos árboles es 256 operaciones"
print(paste("La distancia normalizada entre ambos árboles es",round(distancia_norm,3),sep=" "))
## [1] "La distancia normalizada entre ambos árboles es 0.992"

5. Referencias

Briand, Samuel, Christophe Dessimoz, Nadia El-Mabrouk, Manuel Lafond, and Gabriela Lobinska. 2020. “A Generalized Robinson-Foulds Distance for Labeled Trees.” BMC Genomics 21 (S10). https://doi.org/10.1186/s12864-020-07011-0.
Gamage, Gihan, Nadeeshan Gimhana, Indika Perera, Shanaka Bandara, Thilina Pathirana, Anuradha Wickramarachchi, and Vijini Mallawaarachchi. 2020. “Phylogenetic Tree Construction Using k-Mer Forest- Based Distance Calculation.” International Journal of Online and Biomedical Engineering (iJOE) 16 (07): pp. 4–20. https://doi.org/10.3991/ijoe.v16i07.13807.
Saitou, N, and M Nei. 1987. “The neighbor-joining method: a new method for reconstructing phylogenetic trees.” Molecular Biology and Evolution 4 (4): 406–25. https://doi.org/10.1093/oxfordjournals.molbev.a040454.