Obtencion de datos para la alineación

Antes que nada. Cada genoma puede contener ID’s correspondientes a plasmidos. En este tutorial solo usare los cromosomas de cada especie. Es decir, omitiremos los plásmidos.

for f in *.fna; do seqkit split --by-id $f;done

Luego creo una carpeta (only_chromosomes) y paso ahi los Fasta correspondientes a los cromosomas.

CACTUS

Vamos a usar otro programa para la alineación ya que el archivo maf de salida es mas informativo y facil de reformatear.

Datos de entrada

Primero necestiamos un archivo con los genomas de entrada: Este debe traer el arbol filogenetico en una linea y la etiqueta de los genomas con su respectiva ubicación.

Es muy importante que los nodos del el arbol no estén etiquetados o si lo están, que esta etoqueta sea distinta para cada nodo. De lo contrario habrá errores a la hora de correr cactus

Con el nombre de los nodos me refiero a la etiqueta que esta entre ) y :

((PCC_6605:0.327138,(Z-M001:0.118673,(PCC_8305:0.0945433,PCC_7418:0.0879512)0.652174:0.0310702)0.965638:0.141524)0.432679:0.022499,((PCC_10605:0.134131,((NIES-3709:0.0631663,PCC_6308:0.0700922)0.739832:0.037213,NIES-3708:0.0878375)0.87798:0.0442458)0.973352:0.136561,(NIES-4102:0.189317,(PCC_7437:0.0304549,NIES-3757:0.0301495)0.986676:0.123867)0.779102:0.0610787)0.4326790:0.0224990);

PCC_6605 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/PCC_6605.fna
Z-M001 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/Z-M001.fna
PCC_8305 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/PCC_8305.fna
PCC_7418 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/PCC_7418.fna
PCC_10605 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/PCC_10605.fna
NIES-3709 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/NIES-3709.fna
PCC_6308 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/PCC_6308.fna
NIES-3708 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/NIES-3708.fna
NIES-4102 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/NIES-4102.fna
PCC_7437 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/PCC_7437.fna
NIES-3757 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/fna/only_chromosomes/NIES-3757.fna

Creo un ambiente para correr cactus en la carpeta de cactus Esto que sigue lo hago en la carpeta de cactus

cd cactus
virtualenv -p python3 cactus_env
echo "export PATH=$(pwd)/bin:\$PATH" >> cactus_env/bin/activate
echo "export PYTHONPATH=$(pwd)/lib:\$PYTHONPATH" >> cactus_env/bin/activate
source cactus_env/bin/activate
python3 -m pip install -U setuptools pip
python3 -m pip install -U .
python3 -m pip install -U -r ./toil-requirement.txt

Corro cactus:

cactus ./geminocystis ./examples/evolverGeminocystis.txt ./evolverGeminocystis.hal

como resultado tenemos el archivo evolverPico.hal

Pasamos dicho archivo a formato mfa Usamos a PCC_6308 como referencia:

hal2maf evolverGeminocystis.hal evolverGeminocystis.hal.maf --refGenome PCC_6308 --noAncestors --inMemory

como resultado tenemos el archivo evolverPico.maf

###
##maf version=1 scoring=N/A
# hal ((PCC_6605:0.327138,(Z-M001:0.118673,(PCC_8305:0.0945433,PCC_7418:0.0879512)0.652174:0.0310702)0.965638:0.141524)0.432679:0.022499,((PCC_10605:0.134131,((NIES-3709:0.0631663,PCC_6308:0.0700922)0.739832:0.037213,NIES-3708:0.0878375)0.87798:0.0442458)0.973352:0.136561,(NIES-4102:0.189317,(PCC_7437:0.0304549,NIES-3757:0.0301495)0.986676:0.123867)0.779102:0.0610787)0.4326790:0.022499)Anc00;

a
s   PCC_6308    0   1000    +   4263418 AAACGAACGTACTGCGGCTATTTGGGCTTCCATGGGGCAATCTAGTACCACTACCAGAGAAATTCCACCCCGTGAAATTACTCCAGCACCTGCTACTACTTATACCCCAGAACCAGAAGTTCAACAATTTAATCAACCTATTCGTGGATTGTGGTAATTAATTTTTTGTTGTCATAAACTTCCCTGGGGGTATTGCTCTTTAAAGGTGAGTTTCAGATAAATAAACTTAATATCTTGACAGGCTTTGATACTATAATAATAAAAGTGTCAGAGCTTGTTCTTATTATCGGTAGGTTATTTTAAGTTACCACAAATGAAATCCAATTTTCTTGATTGTTTTTGGGGCAAAAAATCCCTTCTTCCCTTGATTACACTAGGTTTGTGCCTAACTGGTGTAACTGCTTCGATCGCATCCGAATCTACCTTAGCATTAAAAAACGAATTACTGAAAGTAGATTTACCCTCCTTGCCTCCCTTAAAAGAAGCCAATGCTTACCTATATTTTCCCGATGAAGTAGCAACGGTTAAACTTGTATTGAGATTAGGTGCAAAAAAAGTTGAAGTTTATGAACAAGATCAACTAATTGCCAGTTATCCCGTTGCCGTAGGCAAAAAAGGTTGGGAAACTCCCAAAGGGGAATTTGAAATTATCCAAATGATCGAAAATCCCTCTTGGGAAAATCCTTGGACTGGTAAAGTATCCCCCCCGGACCAAATAATCCTTTAGGTGAAAGATGGATCGGATTTTGGAGTAATGGCAAAGATTTTATTGGTTTTCATGGCACTCCCGGAGAACATTTGATCGGTCAAGCCGTATCCCATGGTTGCGTGAGAATGAAAAATAGTGACATCAAAGAACTTTTTAAATTCGTCTCTATGGGTATCCCTGTCACGGTAGTTCAATAATAAGTAGGCGTTGCCCCATAGTCTTTTGATGAGGAAAGGTATTAGGGGTTAGGTTTTAGGTGAAATCATGAAAAAGAAAAAATTTTTAGACTAT

a
s   PCC_6308    1000    1000    +   4263418 GTAATATAAAAAAATAAGCATTTTTAAAAAAGATGTTATCATTTTATGACTAAGGATCAAGAAAGAAGAATACAAGATTTTTTGCCCAAGAAAAGTATGAAAGTGTCATAAGTTCGATCGACGCTATTTTAATAAGTGAAGAAGCTAATAATATTTTAAAAACTTATTTAGCTTTAGCTTATTTTTGTTTAGATAATCAAGACAATTATCAAAGTATTTTACTAGATTTAGTCTTAAAAAGTGACTCAGAAGATTTAATCAATATTGCTCAATGTATCTTTAATATAGCTGATGTTCAAGTAAACAAAGAAAACTTTAATTTAGCCATCGAACTATATCAACAAGGTTTAGAAATTAACCCTCAATACACTCCAGCTTACATTAATTTAGCTCAATTATTCACCGAACAAGGGAATTTTGATTCTGCTATTAGTTTATGGGAAGATTTGATTTTACAATACCCAGATATGATCATTTCTTATGAGCAATTAGGGTTATTATGGCAAAATATTCATGAATATGAAAGTTCCATAAATATTTATCAGCAAGGATTAAATATAGAGTATAATAACTTGAGTATTTTATCAAATTTAGCCTATTGTTATCTTCAAGATAATCAACTTTCTTTAGCCAAAAAATCCTTAGAAAAAATTATTAAACTTGAACCCAATCCGATACAAGCCTATGGAGAATTAGGATATATTTATATTCTCGAAAATAACCTAAATTTAGGAATACAACTATGGCAAAAAATGCTCAAAAAATTTCCTTTTCAAGAATATCTAAATTGGTATAATAATCTTCAAATACCATCAGACAATATCACTTTAAATACTAACTTAATTCGATCGTTACAAAATCATCAAAATCAAGGAGAAATTGCTTTATATATTGCTCATTTGTTATTTAATAAAAAAAATATCTTTTAGCTATTAATTATTATCAAATTGCCTTAGAGAATAATATTGTTGATGAATCTCTTTACTATAACTTAATTCTC

a
s   PCC_6308    2000    81  +   4263418 AGTTTATTCTACACTCAAGAAGCAGTAATATTTCCTGATTTAATTCCTGATCATCTTAACAAATTATATAACATAGATACT

... etc
###

Reformateo las etiquetas de los genomas:

Esta parte lo hago con expresiones regulares en el bloc de notas.

Copio evolverGeminocystis.hal.maf a la carpeta de chromosomes_only y convierto maf a fasta:

mafToFastaStitcher --maf evolverGeminocystis.hal.maf --seqs PCC_6605.fna,Z-M001.fna,PCC_8305.fna,PCC_7418.fna,PCC_10605.fna,NIES-3709.fna,PCC_6308.fna,NIES-3708.fna,NIES-4102.fna,PCC_7437.fna,NIES-3757.fna --breakpointPenalty 5 --interstitialSequence 20 --outMfa evolverGeminocystis.hal.maf.sed.cactus

Trimal

Quito las columnas con muchos gaps:

trimal -in evolverPico.hal.maf.sed.mfa -gt 0.8 -out evolverPico.hal.maf.sed.mfa.trimal

#convierto a philyp

python3 Fasta2Phylip.py only_chromosomes/

Obtengo las coordenadas del palindromo: En esta parte edito la linea28 del código (Spp == ‘PCC_6308’). Spp es la especie que estoy usando como referencia.

python3 AlignmentPalindromeCoords.py only_chromosomes/

El ultimo script tiene como salida el archivo Orthologues_Palindrome_sites.txt el cual contiene 2 columnas, una con el nombre del archivo del ortólogo y otra con el intervalo de sitios en donde se encuentra el palindromo:

###
FILE    PAL START   END
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    1114    1121
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    2332    2339
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    4527    4534
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    6007    6014
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    10468   10475
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    10731   10738
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    11017   11024
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    11826   11833
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    12300   12307
... etc
###

GRAFO

primero cargo las librerias necesarias

library(ggplot2)
library(ggtree)
library(ape)
library(tidyverse)
library(tidytree)
library(phangorn)
library(dplyr)

Preparacion de datos para el grafo

Primero cargo Orthologues_Palindrome_sites.txt del paso anterior.

## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Geminocystis_clade/
Sites <- read.table("fna/Orthologues_Palindrome_sites.txt", sep = "\t", header = TRUE)

Creo la matriz que va a contener unicamente las transiciones

DF <- matrix(0,
             nrow = 0,
             ncol = 2)
colnames(DF) <- c("from","to")

Creo una matriz adicional que tendra las transiciones y una columna extra con la direccion de la transición

{LINKS <- matrix(0,
             nrow = 0,
             ncol = 3)
colnames(LINKS) <- c("from","to","direction")}

Cargo el modelo de evolucion, el metodo de reconstrucción (bayesiano), la filogenia (previamente calculada con orthofinder) la cual enraizo en PCC_8801.

EvolModel = "F81"
Method = "bayes"
Tree = read.tree("SpeciesTree_rooted.txt")#"SpeciesTree_rooted.txt"
#Root = 'PCC_8801'
#Tree <-root(Tree, outgroup = Root, edgelabel = TRUE)
Alignment = "fna/only_chromosomes/evolverGeminocystis.hal.maf.sed.cactus.phy"

Esta filogenia contiene 11 nodos, sin embargo posteriormente para el analisis se quita la raiz y quedan 10.

ggtree(Tree,branch.length="none",) +
  geom_tiplab(color='firebrick', offset = .14)+
  geom_label(aes(label = node),show.legend = FALSE)

En el archivo Orthologues_Palindrome_sites.txt hay 2593 sitios en total. Cargamos el primer sitio el cual corresponde a:

###
FILE    PAL START   END
evolverGeminocystis.hal.maf.sed.cactus.phy  TCGATCGA    1114    1121

###

El numero ORTH correponde al numero de sitio. Por lo tanto cargamos el PATH del alineamiento del primer sitio (Alignment), así como su inicio (SI) y su termino (EI)

ORTH = 1
SI = Sites[ORTH,3]
EI = Sites[ORTH,4]

Alineación multiple

Cargo la alineación multiple

cyanobacterias <- read.phyDat(Alignment, format = "interleaved")

## Cambio el attributo site.patter a FALSE para poder acceder a todos los indices
cyanobacterias <- subset(cyanobacterias, site.pattern = FALSE) 
#TreeR = pratchet(cyanobacterias, trace = 0)|>acctran(cyanobacterias)
parsimony(Tree, cyanobacterias)
## [1] 1027363

Filogenia y parametros del modelo

Estimo los parametros para el modelo. en esta parte uso el modelo de evolucion elegido anteriormente (“F81”). Además en esta parte el arbol se desenraiza y ahora tenemos 10 nodos. Por lo tanto, debo cambiar el arbol al nuevo sin raiz para evitar errores.

fit <- pml(Tree, cyanobacterias)
fit <- optim.pml(fit, model=EvolModel, control = pml.control(trace=0))

Tree2 <- fit[["tree"]]
ggtree(Tree2,branch.length="none",) +
  geom_tiplab(color='firebrick', offset = .14)+
  geom_label(aes(label = node),show.legend = FALSE)

Reconstrucción de Estados Ancestrales

Hago la reconstruccion de estados ancestrales para cada sitio. Para esto hay dos métodos: bayesiano o máxima verosimilitud. Anteriormente habia declarado la variable Method como bayes.

  ## Calculo la probabilidad del árbol filogenético dado el alineamiento y el modelo
  if(Method == "bayes"){
    anc.bayes <- ancestral.pml(fit, type="bayes", return = "prob")
    Reconstruction <- anc.bayes
  }
  if(Method == "ml"){
    anc.ml <- ancestral.pml(fit, "ml")
    Reconstruction <- anc.ml
  }

Armado de palindromos ancestrales

Armo los palindromos ancestrales de acuerdo a los estados ancestrales usando los sitios de inicio y fin de cada sitio palindromico (SI y EI). Posteriormente uso dicho intervalo para buscar lel estado de cada nucleotido del intervalo y guardarlo en el vector (PalindromeNucleotidePositions). Tambien extraigo los nombres de cada nodo de la reconstrucción:

PalInterval = SI:EI ## ESTE ES EL INTERVALO DEL SITIO PALINDROMICO
  
## Obtengo las posiciones para la posicion PalInterval
PalindromeNucleotidePositions <- attributes(cyanobacterias)$index[PalInterval]

## Nombre de los nodos
Nodos <- attributes(Reconstruction)$names

Ahora creo una matriz que contendra 9 columnas, la primera contendrá los Nodos y las 8 restantes contendran los nucleotidos correspondientes a su palindromo y prosigo a rellenar la matriz. Sin embargo, en este paso hay un detalle.

La reconstrucción de estado ancestral para cada sitio arroja una probabilidad para cada uno de los 4 sitios posibles (A, T ,G o C). Si hay equivalencia en la probabilidad en los estados, es decir que el estado ancestral de un nucleotido tenga la misma probabilidad de haber estado en cualesquiera de los 4 posibles estados, entonces dicho estado se omite, de lo contrario se toma aquel estado con la mayor probabilidad. Otro punto importante es que hay sitios en los que una especie puede contener deleciones. Para estos casos, si gfalta algún nucleotido sew rellena con un “-”.

  ## Transformo el arbol en data frame para extraer las etiquetas de los nodos.
  a.1 <- as_tibble(Tree2)

  ## NA.NODE corresponde al nodo a partir del cual ya no estan las puntas del arbol
  NA.NODE = length(cyanobacterias)+1
  
  ## Esto es para agregar la etiqueta de los nodos a los que faltan, es decir a los nodos que no son puntas del arbol.
  for (i in NA.NODE:length(Nodos)){ ## ESTO DEPENDE DE LA ESTRUCTURA DEL ARBOL
    a.1$label[i] <- a.1$node[i] # x$node[i-1]
  }
  
  ## Agrego las etiquetas de nodo al arbol
  b.1 <- tibble::tibble(label = Nodos)
  c.1 <- full_join(a.1, b.1, by = 'label')
  NodeLabels <- c.1$label
  
  TipLabels <- Tree2$tip.label
  NodeLabels2 <- NodeLabels[-c(1:length(TipLabels))]

  ## Creo una matriz de 8 x length(Nodos) (8 nucleotidos x numero de nodos) para guardar los resultados
  ScoresMtx <- matrix(0,
                      nrow = length(Nodos),
                      ncol = 8)
  colnames(ScoresMtx) <- c("1","2","3","4","5","6","7","8")
  rownames(ScoresMtx) <- Nodos
  
  
  ## Relleno la matriz con los scores.
  ## ESTE CICLO ES UNICAMENTE PARA LAS PUNTAS DE LOS NODOS
  j = 1
  i = 1
  k = 0
  for (Position in PalindromeNucleotidePositions){
    for (Nodo in TipLabels){
      PositionScores <- Reconstruction[[Nodo]][Position,] ## Extraigo los scores para el nodo "NODO" en la posicion "POSITION" de la secuencia
      if (length(unique(PositionScores)) == 1){
        k=k+1
        WinnerNucleotide = "-"
        ScoresMtx[i,j] <- WinnerNucleotide
        i=i+1 ## contador para cada Nodo
      }else{
        PositionWinnerScore <- max(PositionScores) ## Extraigo el valor mas grande. Este es el valor de la probabilidad de que el nucleotido ancestral sea ese
        WinnerNucleotide <- which(PositionScores == PositionWinnerScore, arr.ind = T) ## pregunto a que numero (Nucleótido) corresponde ese valor
        ScoresMtx[i,j] <- WinnerNucleotide ## Guardo dicho valor en la matriz de resultados
        i=i+1 ## contador para cada Nodo
      }
    }
    j=j+1 ## contador para cada posición
    i=1 ## Reiniciamoa el contador de Nodo para la siguiente posicion
  }
  
  ## Relleno la matriz con los scores
  ## ESTA PARTE ES PARA TODOS LOS NODOS RESTANTES QUE NO SON PUNTAS
  j = 1
  i = length(TipLabels)+1
  k = 0
  #(length(TipLabels)+1):length(NodeLabels)
  for (Position in PalindromeNucleotidePositions){
    for (Nodo in NodeLabels2){
      PositionScores <- Reconstruction[[Nodo]][Position,] ## Extraigo los scores para el nodo "NODO" en la posicion "POSITION" de la secuencia
      if (length(unique(PositionScores)) == 1){
        k=k+1
        WinnerNucleotide = "X"
        ScoresMtx[i,j] <- WinnerNucleotide
        i=i+1 ## contador para cada Nodo
      }else{
        PositionWinnerScore <- max(PositionScores) ## Extraigo el valor mas grande. Este es el valor de la probabilidad de que el nucleotido ancestral sea ese
        WinnerNucleotide <- which(PositionScores == PositionWinnerScore, arr.ind = T) ## pregunto a que numero (Nucleótido) corresponde ese valor
        ScoresMtx[i,j] <- WinnerNucleotide ## Guardo dicho valor en la matriz de resultados
        i=i+1 ## contador para cada Nodo
      }
    }
    j=j+1 ## contador para cada posición
    i=length(TipLabels)+1 ## Reiniciamoa el contador de Nodo para la siguiente posicion
  }

Ahora decodificamos el palindromo a nucleotidos para esto creo una nueva matriz que contendra dos columnas: nodos y sus correspondientes palíndromos.

## Creo un matriz de 1x14 para guardar la etiqueta (palindromo) de cada nodo
PalsMtx <- matrix(0,
                  nrow = length(Nodos),
                  ncol = 1)
colnames(PalsMtx) <- c("Palindrome")
rownames(PalsMtx) <- Nodos

## Extraigo cada palindromo, cambio el codigo de numero por el codigo de letras y lo guardo en la matriz 1x14
i=0
j=0
for (Nodo in Nodos){
  j=j+1
  Palindrome <- as.character(ScoresMtx[j,]) # Este es el palindromo
  pal =""
  for (nuc in Palindrome){
    if(nuc== 1){nuc = 'A'}
    if(nuc== 2){nuc = 'C'}
    if(nuc== 3){nuc = 'G'}
    if(nuc== 4){nuc = 'T'}
    pal = paste(pal,nuc,sep = "")
  }
  PalsMtx[j,1] <- pal
}

Guardo los nodos y los palindromos de cada nodo en 2 vectores. Los cuales voy a usar para anotar el arbol filogenetico.

l <- as.list(PalsMtx[,1]) 
Nodes = names(l)
NodePals = as.character(PalsMtx[,1])

Anotacion del arbol filogenético

Para poder hacer esto, primero transformo el arbol filogenetico en una tabla y agrego 2 columnas: (label) que contendra el nombre de cada nodo y (Palindromo) que contendra el palindromo de dicho nodo.

## Transformo el arbol en data frame para poder agregar los palindromos de cada nodo
x <- as_tibble(Tree2)

## Extraigo el Boostrap por si se ocupa
NA.NODE = length(cyanobacterias)+1
BootStrap <- x$label[NA.NODE:length(Nodes)]
for (i in NA.NODE:length(Nodes)){ ## ESTO DEPENDE DE LA ESTRUCTURA DEL ARBOL
  x$label[i] <- x$node[i] # x$node[i-1]
}
## Agrego los datos de los palindromos al arbol
d <- tibble::tibble(label = Nodes,
                    Palindromo = NodePals)
y <- full_join(x, d, by = 'label')

Hasta aqui tenemos una filogenoa anotada con su palindromo ancestral en cada nodo:

ggtree(as.treedata(y),branch.length="none") +
    geom_tiplab(color='firebrick', offset = .14)+
    geom_label(aes(label = Palindromo,color=Palindromo),show.legend = FALSE)+
    aes(color=Palindromo)+
    guides(fill = guide_legend(override.aes = list(color = NA)), 
           color = FALSE, 
           shape = 1)

Matriz de transiciones para el grafo

Ahora necesito crear una matriz que tenga las tranmsiciones entre cada nodo. Para esto extraigo las direcciones de los nodos desde el arbol filogenetico previamente anotado es decir las columnas from y to y las guardo en una nueva matriz:

## Esta es la configuracion del arbol
from<-y$label#c("336-3","7","7","8","8","9","9","10","10")
to<-y$parent#c("7","NIES-3974","8","PCC_6303","9","PCC_7716","10","NIES-4071","NIES-4105")

LinksMtx <- matrix(0,
                   nrow = length(from),
                   ncol = 2)
colnames(LinksMtx) <- c("from","to")

for (i in 1:length(from)){
  LinksMtx[i,1] = from[i]
  LinksMtx[i,2] = to[i]
}

Luego extraigo las columnas label y Palindromo y las guardo en otra matriz. Dicha matriz me servira como diccionario para obtener las transiciones en cada nodo. Posteriormente la convierto a data frame.

## AQUI VOY A CREAR LA MATRIZ QUE ME SERVIRA COMO DICCIONARIO PARA RELLENAR LA MATRIZ DE TRANSICIONES
Mat <- y[,4:5]

Mat2 <- as.data.frame(Mat) ## CONVIERTO A DF
Mat2
##        label Palindromo
## 1   PCC_6605   --------
## 2     Z-M001   --------
## 3   PCC_8305   --------
## 4   PCC_7418   --------
## 5  PCC_10605   --------
## 6  NIES-3709   --------
## 7   PCC_6308   TCGATCGA
## 8  NIES-3708   --------
## 9  NIES-4102   --------
## 10  PCC_7437   --------
## 11 NIES-3757   --------
## 12        12   TCGATCGA
## 13        13   TCGATCGA
## 14        14   TCGATCGA
## 15        15   TCGATCGA
## 16        16   TCGATCGA
## 17        17   TCGATCGA
## 18        18   TCGATCGA
## 19        19   TCGATCGA
## 20        20   TCGATCGA

Ahora creo una nueva matriz que va a tener todas las transiciones. Para esto uso las ultimas dos matrices. La matriz LinksMtx que tiene las direcciones entre nodos la paso a palindromos con ayuda de Mat2 la cual me sirve como diccionario y convierto esta ultima matriz (LinksMtx2) a dataframe.

## AQUI CREO LA MATRIZ DE TRANSICIONES.
LinksMtx2 <- matrix(0,
                    nrow = length(from),
                    ncol = 2)
colnames(LinksMtx2) <- c("from","to")

## AQUI RELLENO LA MATRIZ DE TRANSICIONES USANDO COMO DICCIONARIO A LA MATRIZ ANTERIOR (Mat2)
for (i in 1:length(from)) {
  for (j in 1:length(from)){
    if (LinksMtx[i,1]==Mat2[j,1]){
      LinksMtx2[i,1] = Mat2[j,2]
    }
    if (LinksMtx[i,2]==Mat2[j,1]){
      LinksMtx2[i,2] = Mat2[j,2]
    }
  }
}
df <- as.data.frame(LinksMtx2) ## CONVIERTO A DF LA MATRIZ DE TRANSICIONES

Creo una matriz adicional la cual es identica a la anterior salvo que le agrego las direcciones de cada transicion, las cuales podrian ser de utilidad.

directions<-paste(from,to,sep="--") ## ESTE VECTOR CONTIENE LAS DIRECCIONES
Links0 = as.data.frame(LinksMtx2) ## HAGO UNA COPIA DE LA MATRIZ DE TRANSICIONES
Links0$direction<-directions ## AGREGO COMO COLUMNA NUEVA A "directions"

Al final de este paso tendre dos matrices DF y LINKS

DF contiene solo las transiciones

DF <- rbind(DF,df) ## UNO LA MATRIZ ACTUAL DE TRANSICIONES CON LA ANTERIOR PARA CREAR UNA SOLA
DF
##        from       to
## 1  -------- TCGATCGA
## 2  -------- TCGATCGA
## 3  -------- TCGATCGA
## 4  -------- TCGATCGA
## 5  -------- TCGATCGA
## 6  -------- TCGATCGA
## 7  TCGATCGA TCGATCGA
## 8  -------- TCGATCGA
## 9  -------- TCGATCGA
## 10 -------- TCGATCGA
## 11 -------- TCGATCGA
## 12 TCGATCGA TCGATCGA
## 13 TCGATCGA TCGATCGA
## 14 TCGATCGA TCGATCGA
## 15 TCGATCGA TCGATCGA
## 16 TCGATCGA TCGATCGA
## 17 TCGATCGA TCGATCGA
## 18 TCGATCGA TCGATCGA
## 19 TCGATCGA TCGATCGA
## 20 TCGATCGA TCGATCGA

y LINKS las transiciones con sus respectivas direcciones.

LINKS<-rbind(LINKS,Links0) ## UNO LA MATRIZ ACTUAL DE TRANSICIONES(CON DIRECCIONES) CON LA ANTERIOR PARA CREAR UNA SOLA
LINKS
##        from       to     direction
## 1  -------- TCGATCGA  PCC_6605--12
## 2  -------- TCGATCGA    Z-M001--13
## 3  -------- TCGATCGA  PCC_8305--14
## 4  -------- TCGATCGA  PCC_7418--14
## 5  -------- TCGATCGA PCC_10605--16
## 6  -------- TCGATCGA NIES-3709--18
## 7  TCGATCGA TCGATCGA  PCC_6308--18
## 8  -------- TCGATCGA NIES-3708--17
## 9  -------- TCGATCGA NIES-4102--19
## 10 -------- TCGATCGA  PCC_7437--20
## 11 -------- TCGATCGA NIES-3757--20
## 12 TCGATCGA TCGATCGA        12--12
## 13 TCGATCGA TCGATCGA        13--12
## 14 TCGATCGA TCGATCGA        14--13
## 15 TCGATCGA TCGATCGA        15--12
## 16 TCGATCGA TCGATCGA        16--15
## 17 TCGATCGA TCGATCGA        17--16
## 18 TCGATCGA TCGATCGA        18--17
## 19 TCGATCGA TCGATCGA        19--15
## 20 TCGATCGA TCGATCGA        20--19

Finalmente repito esto para los 2592 sitios palindromicos restantes:

#EvolModel = "F81"
#Method = "bayes"
#Tree = read.tree("SpeciesTree_rooted.txt")#"SpeciesTree_rooted.txt"
#Root = '336-3'
#Tree <-root(Tree, outgroup = Root, edgelabel = TRUE)

for (ORTH in 2:length(Sites[,1])){ ## HAY 2593 ortólogos
  #Alignment = paste0("fna/only_chromosomes/",Sites[ORTH,1])
  SI = Sites[ORTH,3]
  EI = Sites[ORTH,4]

  ## Cargo la alineación multiple
  #cyanobacterias <- read.phyDat(Alignment, format = "interleaved")
  
  ## Cambio el attributo site.patter a FALSE para poder acceder a todos los indices
  #cyanobacterias <- subset(cyanobacterias, site.pattern = FALSE) 
  #TreeR = pratchet(cyanobacterias, trace = 0)|>acctran(cyanobacterias)
  #parsimony(Tree, cyanobacterias)
  
  ## Estimo los parametros para el modelo
  #fit <- pml(Tree, cyanobacterias)
  #fit <- optim.pml(fit, model=EvolModel, control = pml.control(trace=0))
  
  #Tree2 <- fit[["tree"]]
  
  ## Calculo la probabilidad del árbol filogenético dado el alineamiento y el modelo
  #if(Method == "bayes"){
  #  anc.bayes <- ancestral.pml(fit, type="bayes", return = "prob")
  #  Reconstruction <- anc.bayes
  #}
  #if(Method == "ml"){
  #  anc.ml <- ancestral.pml(fit, "ml")
  #  Reconstruction <- anc.ml
  #}
  
  #plotAnc(Tree2, anc.bayes, 1)
  PalInterval = SI:EI ## ESTE ES EL INTERVALO DEL SITIO PALINDROMICO
  
  ## Obtengo las posiciones para la posicion PalInterval
  PalindromeNucleotidePositions <- attributes(cyanobacterias)$index[PalInterval]
  
  ## Nombre de los nodos
  Nodos <- attributes(Reconstruction)$names
  
  ########################################################################################################
  ## Transformo el arbol en data frame para poder agregar los palindromos de cada nodo
  a.1 <- as_tibble(Tree2)
  ## Extraigo el Boostrap por si se ocupa
  NA.NODE = length(cyanobacterias)+1
  for (i in NA.NODE:length(Nodos)){ ## ESTO DEPENDE DE LA ESTRUCTURA DEL ARBOL
    a.1$label[i] <- a.1$node[i] # x$node[i-1]
  }
  ## Agrego los datos de los palindromos al arbol
  b.1 <- tibble::tibble(label = Nodos)
  c.1 <- full_join(a.1, b.1, by = 'label')
  NodeLabels <- c.1$label
  
  TipLabels <- Tree2$tip.label
  NodeLabels2 <- NodeLabels[-c(1:length(TipLabels))]

  ## Creo una matriz de 8 x length(Nodos) (8 nucleotidos, numero de nodos) para guardar los resultados
  ScoresMtx <- matrix(0,
                      nrow = length(Nodos),
                      ncol = 8)
  colnames(ScoresMtx) <- c("1","2","3","4","5","6","7","8")
  rownames(ScoresMtx) <- Nodos
  
  
  ## Relleno la matriz con los scores
  j = 1
  i = 1
  k = 0
  for (Position in PalindromeNucleotidePositions){
    for (Nodo in TipLabels){
      PositionScores <- Reconstruction[[Nodo]][Position,] ## Extraigo los scores para el nodo "NODO" en la posicion "POSITION" de la secuencia
      if (length(unique(PositionScores)) == 1){
        k=k+1
        WinnerNucleotide = "-"
        ScoresMtx[i,j] <- WinnerNucleotide
        i=i+1 ## contador para cada Nodo
      }else{
        PositionWinnerScore <- max(PositionScores) ## Extraigo el valor mas grande. Este es el valor de la probabilidad de que el nucleotido ancestral sea ese
        WinnerNucleotide <- which(PositionScores == PositionWinnerScore, arr.ind = T) ## pregunto a que numero (Nucleótido) corresponde ese valor
        ScoresMtx[i,j] <- WinnerNucleotide ## Guardo dicho valor en la matriz de resultados
        i=i+1 ## contador para cada Nodo
      }
    }
    j=j+1 ## contador para cada posición
    i=1 ## Reiniciamoa el contador de Nodo para la siguiente posicion
  }
  
  ## Relleno la matriz con los scores
  j = 1
  i = length(TipLabels)+1
  k = 0
  #(length(TipLabels)+1):length(NodeLabels)
  for (Position in PalindromeNucleotidePositions){
    for (Nodo in NodeLabels2){
      PositionScores <- Reconstruction[[Nodo]][Position,] ## Extraigo los scores para el nodo "NODO" en la posicion "POSITION" de la secuencia
      if (length(unique(PositionScores)) == 1){
        k=k+1
        WinnerNucleotide = "X"
        ScoresMtx[i,j] <- WinnerNucleotide
        i=i+1 ## contador para cada Nodo
      }else{
        PositionWinnerScore <- max(PositionScores) ## Extraigo el valor mas grande. Este es el valor de la probabilidad de que el nucleotido ancestral sea ese
        WinnerNucleotide <- which(PositionScores == PositionWinnerScore, arr.ind = T) ## pregunto a que numero (Nucleótido) corresponde ese valor
        ScoresMtx[i,j] <- WinnerNucleotide ## Guardo dicho valor en la matriz de resultados
        i=i+1 ## contador para cada Nodo
      }
    }
    j=j+1 ## contador para cada posición
    i=length(TipLabels)+1 ## Reiniciamoa el contador de Nodo para la siguiente posicion
  }
  ########################################################################################################
  
  ## EN ESTA PARTE ARMAREMOS EL PALINDROMO PARA CADA NODO
  
  ## Creo un matriz de 1x14 para guardar la etiqueta (palindromo) de cada nodo
  PalsMtx <- matrix(0,
                    nrow = length(Nodos),
                    ncol = 1)
  colnames(PalsMtx) <- c("Palindrome")
  rownames(PalsMtx) <- Nodos
  
  ## Extraigo cada palindromo, cambio el codigo de numero por el codigo de letras y lo guardo en la matriz 1x14
  i=0
  j=0
  for (Nodo in Nodos){
    j=j+1
    Palindrome <- as.character(ScoresMtx[j,]) # Este es el palindromo
    pal =""
    for (nuc in Palindrome){
      if(nuc== 1){nuc = 'A'}
      if(nuc== 2){nuc = 'C'}
      if(nuc== 3){nuc = 'G'}
      if(nuc== 4){nuc = 'T'}
      pal = paste(pal,nuc,sep = "")
    }
    PalsMtx[j,1] <- pal
  }
  
  ## ESTA PATRTE ES PARA GUARDAR LOS PALINDROMOS DE CADA NODO
  l <- as.list(PalsMtx[,1]) 
  Nodes = names(l)
  NodePals = as.character(PalsMtx[,1])
  
  ## Transformo el arbol en data frame para poder agregar los palindromos de cada nodo
  x <- as_tibble(Tree2)
  
  ## Extraigo el Boostrap por si se ocupa
  NA.NODE = length(cyanobacterias)+1
  BootStrap <- x$label[NA.NODE:length(Nodes)]
  for (i in NA.NODE:length(Nodes)){ ## ESTO DEPENDE DE LA ESTRUCTURA DEL ARBOL
    x$label[i] <- x$node[i] # x$node[i-1]
  }
  ## Agrego los datos de los palindromos al arbol
  d <- tibble::tibble(label = Nodes,
                      Palindromo = NodePals)
  y <- full_join(x, d, by = 'label')
  
  #------------------------
  
  ## Esta es la configuracion del arbol
  from<-y$label#c("336-3","7","7","8","8","9","9","10","10")
  to<-y$parent#c("7","NIES-3974","8","PCC_6303","9","PCC_7716","10","NIES-4071","NIES-4105")
  
  LinksMtx <- matrix(0,
                     nrow = length(from),
                     ncol = 2)
  colnames(LinksMtx) <- c("from","to")
  
  for (i in 1:length(from)){
    LinksMtx[i,1] = from[i]
    LinksMtx[i,2] = to[i]
  }
  
  ## AQUI VOY A CREAR LA MATRIZ QUE ME SERVIRA COMO DICCIONARIO PARA RELLENAR LA MATRIZ DE TRANSICIONES
  Mat <- y[,4:5]
  #Mat <- Mat %>% ## AQJUI QUITO EL RENGLON 7 QUE NO CONTIENE NADA. ESTE CORRESPONDE A LA RAIZ DEL ARBOL
  #  filter(!is.na(Palindromo))
  Mat
  
  Mat2 <- as.data.frame(Mat) ## CONVIERTO A DF
  Mat2
  
  ## AQUI CREO LA MATRIZ DE TRANSICIONES.
  LinksMtx2 <- matrix(0,
                      nrow = length(from),
                      ncol = 2)
  colnames(LinksMtx2) <- c("from","to")
  
  ## AQUI RELLENO LA MATRIZ DE TRANSICIONES USANDO COMO DICCIONARIO A LA MATRIZ ANTERIOR (Mat2)
  for (i in 1:length(from)) {
    for (j in 1:length(from)){
      if (LinksMtx[i,1]==Mat2[j,1]){
        LinksMtx2[i,1] = Mat2[j,2]
      }
      if (LinksMtx[i,2]==Mat2[j,1]){
        LinksMtx2[i,2] = Mat2[j,2]
      }
    }
  }
  
  directions<-paste(from,to,sep="--") ## ESTE VECTOR CONTIENE LAS DIRECCIONES
  Links0 = as.data.frame(LinksMtx2) ## HAGO UNA COPIA DE LA MATRIZ DE TRANSICIONES
  Links0$direction<-directions ## AGREGO COMO COLUMNA NUEVA A "directions"
  Links0
  LINKS<-rbind(LINKS,Links0) ## UNO LA MATRIZ ACTUAL DE TRANSICIONES(CON DIRECCIONES) CON LA ANTERIOR PARA CREAR UNA SOLA
  
  df <- as.data.frame(LinksMtx2) ## CONVIERTO A DF LA MATRIZ DE TRANSICIONES
  DF <- rbind(DF,df) ## UNO LA MATRIZ ACTUAL DE TRANSICIONES CON LA ANTERIOR PARA CREAR UNA SOLA
  #print(paste0("Archivo ",ORTH, " de ",length(Sites$FILE),"."))
  write.table(LINKS,file="LINKS.txt", sep="\t", row.names = FALSE)
}

Ahora las matrices DF y LINKS contienen todas las transiciones de todos los sitios. En total son 25930 transiciones ya que son 10 transiciones por arbol y tenemos 2593 arboles.

length(DF$from)
## [1] 58260
length(LINKS$from)
## [1] 58260

Filtrado de datos para el GRAFO

Primero cuento el numero de veces que se repite cada linea de la matriz de transiciones. El conteo para cada transicion sera el peso de cada vertice. Luego cambio el nombre de la columna counts a weight, paso su contenido a numeros enteros y la convierto a tibble. En total tengo 25930 transiciones:

library(data.table)
RowCts <- setDT(DF)[,list(Count=as.numeric(.N)),names(DF)] ## CUENTO EL NUMERO DE OCURRENCIAS DE CADA TRANSICIÓN #7211
colnames(RowCts) <- c("from","to","weight") 
RowCts = transform(RowCts, weight = as.integer(weight)) ## CONVIERTO A ENTEROS LA COLUMNA "weight"
RowCts <-as_tibble(RowCts) ## CONVIERTO A TIBBLE
RowCts #7211
## # A tibble: 3,604 × 3
##    from     to       weight
##    <chr>    <chr>     <int>
##  1 -------- TCGATCGA  22297
##  2 TCGATCGA TCGATCGA  23123
##  3 -------- TTGATCGA    105
##  4 TCGATCGA TTGATCGA     27
##  5 TTGATCGA TTGATCGA    135
##  6 -------- TCGATCAA     97
##  7 TCGATCGA TCGATCAA     31
##  8 TCGATCAA TCGATCAA    136
##  9 -------- GCAATCGA      4
## 10 GCGATCGA GCGATCGA    239
## # ℹ 3,594 more rows

Filtro conteos bajos (<=10) y loops. Al final me quedo con 63 transiciones.

RowCts<-RowCts%>% ##QUITO LOS LOOPS. ES DECIR QUITO AQUELLAS TRANSICIONES QUE VAYAN ASI MISMA
  filter(from!=to)
length(RowCts$from) #5604
## [1] 2854
RowCts<-RowCts%>% ## QUIRO AQUELLAS TRANSICIONES CON UN PESO MENOR A 10
  filter(weight>=7)
length(RowCts$from) #63
## [1] 290

Extraigo todos los nodos posibles y los agrego a una matriz que va a contener los nodos del grafo. Adiconalmente le agrego una columna score que me indicara que tan parecido es el aplindromo a HIP1 (despues la quito porque no creo que sea un buen índice)

GraphNodes <- c(RowCts$from, RowCts$to) ## OBTENGO LOS PALINDROMOS DE CADA NODO
length(GraphNodes) #126
## [1] 580
GNMtx <- matrix(0, ## CREO LA MATRIZ QUE CONTENDRÁ LAS NODOS PARA EL GRAFO
                nrow = length(GraphNodes),
                ncol = 2)
colnames(GNMtx) <- c("pal","score") ## ETIQUETO LAS COLUMNAS DE LA MATRIZ

for (i in 1:length(GraphNodes)){ ## RELLENO LA MATRIZ
  GNMtx[i,1] = GraphNodes[i] ## AGREGO LOS NODOS
  GNMtx[i,2] = RecordLinkage::levenshteinSim(GraphNodes[i], 'GCGATCGC')*1000 ## AGREGO UN SCORE QUE SERÁ EL PORCENTAJE DE IDENTIDAD CON HIP1 MULTIPLICADO POR 1000
}

GNMtx <- as.data.frame(GNMtx) ## CONVIERTO LA MATRIZ DE NODOS A DF
GNMtx <- GNMtx[!duplicated(GNMtx), ] ## ELIMINO LOS NODOS DUPLICADOS

Creo un vector que enumerará los palindromos. Dicha enumeracion correspondera a su ID

palindromes<-GNMtx[,1] #34  ## EXTRAIGO LOS NODOS DE LA MATRIZ DEL GRAFO. ESTO LO HAGO PARA POSTERIAMENTE ENUMERARLAS Y ASOCIAR CADA NUMERO COMO ID

ids<-c() ## CREO UN VECTOR VACIO QUE CONTENDRÁ LOS ID'S
for (k in 1:length(palindromes)){ ## ENUMERO LOS PALINDROMOS. ESTOS NUMEROS SERAN LOS IDS
  ids <-append(ids, k)
}
ids
##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
##  [19]  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
##  [37]  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
##  [55]  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
##  [73]  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
##  [91]  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108
## [109] 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
## [127] 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
## [145] 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
## [163] 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## [181] 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
## [199] 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
## [217] 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
## [235] 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
## [253] 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
length(ids) #34
## [1] 269
length(GNMtx[,1]) #34
## [1] 269

Creo una matriz nueva que contendrá los nodos del grafo pero con su ID para cada palindromo:

Nodes2 <-cbind(GNMtx,ids) ## CREO UNA MATRIZ NUEVA CON LA MATRIZ DE NODOS Y LOS IDS. ESTA LA VOY A USAR COMO UN DICCIONARIO ID-NODO
length(Nodes2$pal)
## [1] 269
Nodes2<-Nodes2[,c(3,1,2)] ## REORDENO LAS COLUMNAS DE LA MATRIZ
colnames(Nodes2)<-c("id","label","score") ## ETIQUETO LAS COLUMNAS

Creo la matriz de vertices del grafo y tambien sustituyo el nombre de los nodos por sus respectivos ID’S. Adicionalmente cambio los scores que tengan 0 por 1’s ya que hay errores. Sin embargo, la columna score no la uso mas adelante.

GEMtx <- RowCts ## ESTA VA A SER LA MATRIZ DE VERTICES
Edges2 <- as.matrix(RowCts)#RowCts ESTA ES OTRA MATRIZ DE VERTICES DE LA CUAL VOY A SUSTITUIR LOS NODOS POR SUS IDS

for (n in 1:length(Nodes2[,2])){ ##AQUI CAMBIARE LOS PALINDROMOS DE LOS NODOS POR IDS
  Edges2[Edges2==Nodes2[n,2]] <- as.numeric(Nodes2[n,1])
  #print(n)
}

for (n in 1:length(Nodes2[,2])){ ## HANDLE 0's
  Nodes2[Nodes2==0] <- as.numeric(1)
  #print(n)
}

VISUALIZACIÓN DEL GRAFO

Cargo las librerias

library(visNetwork)
library(networkD3)
library(tidygraph)
library(ggraph)

Para empezar cargo los nodos y los vertices. En esta parte quito la columna score de los nodos.

nodes = as_tibble(Nodes2)
nodes = select(nodes, -score)

edges = as_tibble(Edges2)
edges = transform(edges, from = as.integer(from))
edges = transform(edges, to = as.integer(to))
edges = transform(edges, weight = as.integer(weight))
edges = as_tibble(edges)

Creo un grafo dirigido agregando un peso de acuerdo a la columna weight.

routes_tidy <- tbl_graph(nodes = nodes,
                         edges = edges,
                         directed = TRUE)

routes_tidy
## # A tbl_graph: 269 nodes and 290 edges
## #
## # A directed simple graph with 1 component
## #
## # A tibble: 269 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGA
## 3     3 TCGATCAA
## 4     4 TCGATCGC
## 5     5 GCGATCGC
## 6     6 GCGATCAA
## # ℹ 263 more rows
## #
## # A tibble: 290 × 3
##    from    to weight
##   <int> <int>  <int>
## 1     1     2  22297
## 2     1     7    105
## 3     2     7     27
## # ℹ 287 more rows
routes_tidy %>% 
  activate(edges) %>% 
  arrange(desc(weight))
## # A tbl_graph: 269 nodes and 290 edges
## #
## # A directed simple graph with 1 component
## #
## # A tibble: 290 × 3
##    from    to weight
##   <int> <int>  <int>
## 1     1     2  22297
## 2     1    38    144
## 3     1    11    137
## 4     1     4    113
## 5     1    10    112
## 6     1     7    105
## # ℹ 284 more rows
## #
## # A tibble: 269 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGA
## 3     3 TCGATCAA
## # ℹ 266 more rows

Visualizo el grafo:

ggraph(routes_tidy, layout = "graphopt") + 
  geom_node_point() +
  geom_edge_link(aes(width = weight), alpha = 0.8) + 
  scale_edge_width(range = c(0.2, 2)) +
  geom_node_text(aes(label = label), repel = TRUE) +
  labs(edge_width = "Times") +
  theme_graph()

GRAFOS INTERACTIVOS

Este es un grafo sin peso en los vertices.

visNetwork(nodes, edges)

Este es el mismo grafo pero con peso en sus vertices

edges <- mutate(edges, width = weight/5 + 1)
visNetwork(nodes, edges) %>% 
  visIgraphLayout(layout = "layout_with_fr") %>% 
  visEdges(arrows = "middle")

Este es otro grafo de fuerzas con peso en sus vertices y te muestra las conceciones de cada nodo de manera mas visual

nodes_d3 <- mutate(nodes, id = id - 1)
edges_d3 <- mutate(edges, from = from - 1, to = to - 1)

forceNetwork(Links = edges_d3, Nodes = nodes_d3,
             Source = "from", Target = "to", 
             NodeID = "label", Group = "id", Value = "weight", 
             opacity = 1, fontSize = 16, zoom = TRUE)

Este ultimo grafo muestra el grafo anterior pero de una forma mas analizable.

sankeyNetwork(Links = edges_d3, Nodes = nodes_d3, Source = "from", Target = "to", 
              NodeID = "label", Value = "weight", fontSize = 16, unit = "TIME(s)")

GRAFOS POR DIRECCION

WGHT = 4

PCC_6605–12

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 20 nodes and 19 edges
## #
## # A rooted tree
## #
## # A tibble: 19 × 4
##    from    to direction    weight
##   <int> <int> <chr>         <int>
## 1     1     2 PCC_6605--12   2246
## 2     1     5 PCC_6605--12     19
## 3     1    11 PCC_6605--12     16
## 4     1     6 PCC_6605--12     13
## 5     1     3 PCC_6605--12     12
## 6     1    10 PCC_6605--12     12
## # ℹ 13 more rows
## #
## # A tibble: 20 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGA
## 3     3 TTGATCGA
## # ℹ 17 more rows

Z-M001–13

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 21 nodes and 20 edges
## #
## # A rooted tree
## #
## # A tibble: 20 × 4
##    from    to direction  weight
##   <int> <int> <chr>       <int>
## 1     1     2 Z-M001--13   2248
## 2     1     5 Z-M001--13     17
## 3     1    11 Z-M001--13     17
## 4     1     6 Z-M001--13     15
## 5     1     3 Z-M001--13     12
## 6     1     4 Z-M001--13     12
## # ℹ 14 more rows
## #
## # A tibble: 21 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGA
## 3     3 TTGATCGA
## # ℹ 18 more rows

}

PCC_8305–14

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 23 nodes and 23 edges
## #
## # A directed acyclic simple graph with 1 component
## #
## # A tibble: 23 × 4
##    from    to direction    weight
##   <int> <int> <chr>         <int>
## 1     1     2 PCC_8305--14   2248
## 2     1    11 PCC_8305--14     18
## 3     1     5 PCC_8305--14     16
## 4     1     6 PCC_8305--14     15
## 5     1    10 PCC_8305--14     14
## 6     1     3 PCC_8305--14     12
## # ℹ 17 more rows
## #
## # A tibble: 23 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGA
## 3     3 TTGATCGA
## # ℹ 20 more rows

PCC_7418–14

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 23 nodes and 22 edges
## #
## # A rooted tree
## #
## # A tibble: 22 × 4
##    from    to direction    weight
##   <int> <int> <chr>         <int>
## 1     1     2 PCC_7418--14   2249
## 2     1    11 PCC_7418--14     18
## 3     1     5 PCC_7418--14     15
## 4     1     6 PCC_7418--14     15
## 5     1    10 PCC_7418--14     14
## 6     1     3 PCC_7418--14     12
## # ℹ 16 more rows
## #
## # A tibble: 23 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGA
## 3     3 TTGATCGA
## # ℹ 20 more rows

PCC_10605–16

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 12 nodes and 12 edges
## #
## # A directed acyclic simple graph with 1 component
## #
## # A tibble: 12 × 4
##    from    to direction     weight
##   <int> <int> <chr>          <int>
## 1     1     4 PCC_10605--16   2236
## 2     1     9 PCC_10605--16     14
## 3     1     8 PCC_10605--16     11
## 4     1     6 PCC_10605--16      9
## 5     1     5 PCC_10605--16      7
## 6     1     7 PCC_10605--16      7
## # ℹ 6 more rows
## #
## # A tibble: 12 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGC
## 3     3 GCGATCGC
## # ℹ 9 more rows

}

NIES-3709–18

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 5 nodes and 5 edges
## #
## # A directed acyclic simple graph with 1 component
## #
## # A tibble: 5 × 4
##    from    to direction     weight
##   <int> <int> <chr>          <int>
## 1     1     4 NIES-3709--18   2122
## 2     1     2 NIES-3709--18      5
## 3     1     5 NIES-3709--18      4
## 4     2     4 NIES-3709--18      4
## 5     3     4 NIES-3709--18      4
## #
## # A tibble: 5 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TTGATCGA
## 3     3 ACGATCGA
## # ℹ 2 more rows

PCC_6308–18

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 12 nodes and 11 edges
## #
## # A rooted tree
## #
## # A tibble: 11 × 4
##    from    to direction    weight
##   <int> <int> <chr>         <int>
## 1     1     9 PCC_6308--18     27
## 2     1     8 PCC_6308--18     24
## 3     1     3 PCC_6308--18     18
## 4     1     6 PCC_6308--18     16
## 5     1     2 PCC_6308--18     14
## 6     1     4 PCC_6308--18     13
## # ℹ 5 more rows
## #
## # A tibble: 12 × 2
##      id label   
##   <int> <chr>   
## 1     1 TCGATCGA
## 2     2 TTGATCGA
## 3     3 TCGATCAA
## # ℹ 9 more rows

NIES-3708–17

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 11 nodes and 12 edges
## #
## # A directed acyclic simple graph with 1 component
## #
## # A tibble: 12 × 4
##    from    to direction     weight
##   <int> <int> <chr>          <int>
## 1     1     6 NIES-3708--17   2205
## 2     1    10 NIES-3708--17     12
## 3     1     8 NIES-3708--17     11
## 4     5     6 NIES-3708--17     10
## 5     2     6 NIES-3708--17      8
## 6     1     5 NIES-3708--17      7
## # ℹ 6 more rows
## #
## # A tibble: 11 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCAA
## 3     3 GCGATCAA
## # ℹ 8 more rows

}

NIES-4102–19

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 17 nodes and 16 edges
## #
## # A rooted tree
## #
## # A tibble: 16 × 4
##    from    to direction     weight
##   <int> <int> <chr>          <int>
## 1     1     2 NIES-4102--19   2250
## 2     1     5 NIES-4102--19     18
## 3     1     9 NIES-4102--19     16
## 4     1     6 NIES-4102--19     15
## 5     1     3 NIES-4102--19     12
## 6     1     8 NIES-4102--19     12
## # ℹ 10 more rows
## #
## # A tibble: 17 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGA
## 3     3 TTGATCGA
## # ℹ 14 more rows

PCC_7437–20

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 17 nodes and 16 edges
## #
## # A rooted tree
## #
## # A tibble: 16 × 4
##    from    to direction    weight
##   <int> <int> <chr>         <int>
## 1     1     2 PCC_7437--20   2246
## 2     1     5 PCC_7437--20     19
## 3     1     9 PCC_7437--20     16
## 4     1     6 PCC_7437--20     14
## 5     1     3 PCC_7437--20     13
## 6     1     8 PCC_7437--20     12
## # ℹ 10 more rows
## #
## # A tibble: 17 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGA
## 3     3 TTGATCGA
## # ℹ 14 more rows

NIES-3757–20

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 17 nodes and 16 edges
## #
## # A rooted tree
## #
## # A tibble: 16 × 4
##    from    to direction     weight
##   <int> <int> <chr>          <int>
## 1     1     2 NIES-3757--20   2247
## 2     1     5 NIES-3757--20     19
## 3     1     9 NIES-3757--20     16
## 4     1     6 NIES-3757--20     14
## 5     1     3 NIES-3757--20     13
## 6     1     8 NIES-3757--20     12
## # ℹ 10 more rows
## #
## # A tibble: 17 × 2
##      id label   
##   <int> <chr>   
## 1     1 --------
## 2     2 TCGATCGA
## 3     3 TTGATCGA
## # ℹ 14 more rows

}

13–12

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 2 nodes and 1 edges
## #
## # A rooted tree
## #
## # A tibble: 1 × 4
##    from    to direction weight
##   <int> <int> <chr>      <int>
## 1     1     2 13--12         2
## #
## # A tibble: 2 × 2
##      id label   
##   <int> <chr>   
## 1     1 TTGACCGC
## 2     2 TTGATCGC

14–13 WEIGHT <= 2

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 99 nodes and 50 edges
## #
## # A rooted forest with 49 trees
## #
## # A tibble: 50 × 4
##    from    to direction weight
##   <int> <int> <chr>      <int>
## 1     1    50 14--13         1
## 2     2    51 14--13         1
## 3     3    52 14--13         1
## 4     4    53 14--13         1
## 5     5    54 14--13         1
## 6     6    55 14--13         1
## # ℹ 44 more rows
## #
## # A tibble: 99 × 2
##      id label   
##   <int> <chr>   
## 1     1 ATTTTCAA
## 2     2 GGCTTCGA
## 3     3 CCGTTCTA
## # ℹ 96 more rows

15–12 WEIGHT <= 1

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 120 nodes and 66 edges
## #
## # A rooted forest with 54 trees
## #
## # A tibble: 66 × 4
##    from    to direction weight
##   <int> <int> <chr>      <int>
## 1     1    61 15--12         1
## 2     2    62 15--12         1
## 3     2    63 15--12         1
## 4     3    64 15--12         1
## 5     2    65 15--12         1
## 6     4    66 15--12         1
## # ℹ 60 more rows
## #
## # A tibble: 120 × 2
##      id label   
##   <int> <chr>   
## 1     1 TCGATCAA
## 2     2 TCGATCGA
## 3     3 TAGATAAA
## # ℹ 117 more rows

}

16–15 weight>=3

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 4 nodes and 2 edges
## #
## # A rooted forest with 2 trees
## #
## # A tibble: 2 × 4
##    from    to direction weight
##   <int> <int> <chr>      <int>
## 1     1     3 16--15         3
## 2     2     4 16--15         3
## #
## # A tibble: 4 × 2
##      id label   
##   <int> <chr>   
## 1     1 TCGATCGC
## 2     2 GCGATCGA
## 3     3 TCGATCGA
## # ℹ 1 more row

17–16

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 3 nodes and 2 edges
## #
## # A rooted tree
## #
## # A tibble: 2 × 4
##    from    to direction weight
##   <int> <int> <chr>      <int>
## 1     1     3 17--16         8
## 2     1     2 17--16         7
## #
## # A tibble: 3 × 2
##      id label   
##   <int> <chr>   
## 1     1 TCGATCGA
## 2     2 TCGATCGC
## 3     3 GCGATCGA

18–17

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 10 nodes and 9 edges
## #
## # A rooted tree
## #
## # A tibble: 9 × 4
##    from    to direction weight
##   <int> <int> <chr>      <int>
## 1     1     5 18--17        10
## 2     1     7 18--17        10
## 3     1     4 18--17         8
## 4     1     6 18--17         8
## 5     1     9 18--17         5
## 6     1     2 18--17         4
## # ℹ 3 more rows
## #
## # A tibble: 10 × 2
##      id label   
##   <int> <chr>   
## 1     1 TCGATCGA
## 2     2 TCTATCGA
## 3     3 TCGCTCGA
## # ℹ 7 more rows

}

19–15 weight>=1

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 119 nodes and 63 edges
## #
## # A rooted forest with 56 trees
## #
## # A tibble: 63 × 4
##    from    to direction weight
##   <int> <int> <chr>      <int>
## 1     1    64 19--15         1
## 2     2    65 19--15         1
## 3     3    66 19--15         1
## 4     4    67 19--15         1
## 5     5    68 19--15         1
## 6     6    69 19--15         1
## # ℹ 57 more rows
## #
## # A tibble: 119 × 2
##      id label   
##   <int> <chr>   
## 1     1 TCGCTCAA
## 2     2 GAAATAAA
## 3     3 TTTTTCGT
## # ℹ 116 more rows

20–19 WEIGHT <= 2

ggtree(Tree2) +
  geom_tiplab(color='firebrick', offset = .10, hjust=2.0)+
  geom_label(aes(label = node),show.legend = TRUE)

## [1] 9226
## # A tbl_graph: 2 nodes and 1 edges
## #
## # A rooted tree
## #
## # A tibble: 1 × 4
##    from    to direction weight
##   <int> <int> <chr>      <int>
## 1     1     2 20--19         2
## #
## # A tibble: 2 × 2
##      id label   
##   <int> <chr>   
## 1     1 TTGACCGA
## 2     2 TTGATCGA