Reconstrucción Ancestral de Secuencias - GENOMAS - CGGCGCCG - PCC_7716
Eduardo PM
2023-03-13
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.
## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/fna
for f in *.fna; do seqkit split --by-id $f;doneLuego 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 :
Por ejemplo:
En este arbol la primera y la ultima etiqueta de nodo son identicas, esto ocasionó problemas a la hora de ejecutar cactus:
###
(336-3:0.101119,(NIES-3974:0.218161,(PCC_6303:0.172103,(PCC_7716:0.058496,(NIES-4071:8.62091e-07,NIES-4105:8.62091e-07)0.986266:0.0633482)0.960653:0.109454)0.507795:0.0389476)1:0.101119);
Para que el programa funcionara unicamente cambié el 8.62091e-07 por 8.620912e-07:
(336-3:0.101119,(NIES-3974:0.218161,(PCC_6303:0.172103,(PCC_7716:0.058496,(NIES-4071:8.62091e-07,NIES-4105:8.620912e-07)0.986266:0.0633482)0.960653:0.109454)0.507795:0.0389476)1:0.101119);
(336-3:0.101119,(NIES-3974:0.218161,(PCC_6303:0.172103,(PCC_7716:0.058496,(NIES-4071:8.62091e-07,NIES-4105:8.620912e-07)0.986266:0.0633482)0.960653:0.109454)0.507795:0.0389476)1:0.1011191);
336-3 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/fna/only_chromosomes/336-3.fna
NIES-4071 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/fna/only_chromosomes/NIES-4071.fna
PCC_6303 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/fna/only_chromosomes/PCC_6303.fna
NIES-3974 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/fna/only_chromosomes/NIES-3974.fna
NIES-4105 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/fna/only_chromosomes/NIES-4105.fna
PCC_7716 /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/fna/only_chromosomes/PCC_7716.fnaCreo un ambiente para correr cactus en la carpeta de cactus Esto que sigue lo hago en la carpeta de cactus
## /home/lalibelulalo/software/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.txtCorro cactus:
## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade
cd ../../PIPELINES_2023/ASR_GENOMES/Callothrix_clade/
cactus ./calothrix ./evolverCalothrix.txt ./evolverCalothrix.halcomo resultado tenemos el archivo evolverPico.hal
Convierto a MAF
Para este paso necesitamos una referencia. Esta referencia será el genoma con mayor conteo de palindromos. Para esto corremos el siguiente script:
## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade
python3 CountPalsInOrthologues.py fna/only_chromosomes/ pals.list genomes.list 3## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade
File = "../../Markov_count_CGGCGCCG_2023-5-9_14hrs14mins.txt"
table = read.table(File,header = TRUE,sep = "\t")
table = table[order(table$obs,decreasing = TRUE),]
#knitr::kable(table, align = "lccrr")
rmarkdown::paged_table(table,options = list(rows.print = 10, cols.print = 4))Usamos a PCC_7716 como referencia:
Pasamos dicho archivo a formato MAF
## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade
hal2maf evolverCalothrix.hal evolverCalothrix.PCC_7716.hal.maf --refGenome PCC_7716 --noAncestors --inMemorycomo resultado tenemos el archivo evolverCallothrix.PCC_7716.hal.maf
###
##maf version=1 scoring=N/A
# hal (336-3:0.101119,(NIES-3974:0.218161,(PCC_6303:0.172103,(PCC_7716:0.058496,(NIES-4071:8.62091e-07,NIES-4105:8.62091e-07)0.986266:0.0633482)0.960653:0.109454)0.507795:0.0389476)1:0.101119)Anc0;
a
s 336-3.336-3 0 1000 + 6283267 GCACTGCCCTACTGGAAGCACTTCGTGGCACTTACTGCCAACATGTGGTTCAAGTCACAGAAGAACCCTTCATTTCTTTACAAACTATTCCCCTTCACCGGGGACGGAAACTTGTCAATTATTGCGCCTCCAAAACCACAGATGGAAACCACTTTACAAACTATTCCCCTTCACCGGGGACGAAAACGCTTGCATTAGCCCACAGGTCAGTGGACTGACTAATTATTTCCACAACTATGGGCAAACTAGCTTACAGCCGCTACCATCTGAGTTGACAGGGAATCAAAAATCATAGTATTATAATCTGGACAAGGAAGAACTATCTTTCAATAATAAACCTGTGAATATTTCTCCATCTCCAAAGCAGCCGCGTGTCTCATGGCAGGCGGCTTTCTCCATGTTGATGTACATTTTCATGTATCTACCAATATTGGTGCTTGCCTTCTACAGCTTCAACAATTCAGCTTACAGTGCTGGTTGGCAAGGATTCACCCTCGATTGGTATCGCCAGCTATTTGCGGATGAACGCATCTTATTAGCGACCAAAAACAGTCTCCTAGTAGCCAGTTGTGCGGTGGGGATTGCGGCAGTATTTGGAACCCTGATGGCAGTGGGATTAGCTCGATATCAGTTTCCCGGCAAAGGTTTATATAAAGGAGTCAGTTACCTACCCTTAATTATTCCCGACATTGCGATCGCCGTTGCCACCCTAGTTTTTCTGGCTGCTTTTGCCATACCTCTAAGTCTATGGACAATCATCGCCGCCCATGTGGTTTTCTGTCTTGCCTACATTGCCCTAGTTGTCTCCTCCCGACTAACAAACCTGGATTCCCATTTAGAAGAAGCCGCCCTAGATTTAGGTGCTACACCTACCCAAGCATTTATTAAAGTTTTAATACCCCAGTTAATGCCGGGTATAATTGCAGGCTGTCTGTTAGCTTTTATCCTTAGCTTAGATGACTTTCTGATTGCCAGTTTTACCTCCGGTAGCGGCTCTAGT
a
s 336-3.336-3 1000 1000 + 6283267 ACCCTGCCCATGGAAATTTTTAGTCGCCTGCGAACAGGTGTTAAACCAGACATAAATGCCCTCAGTGTTATCCTAATTCTCGTATCAGCAATTGTAGCTTTTATTGCCGAATCAATTAGAATTGCTGGAGACAAAAAGCACAGATAAAAAAGTATTTATACCACAATCTCTATCATCGTAATTCCTCTGAGCCACCAGGCAAACATTCTGCAATGCCATAGACAATAAATACTTATTTCATATAAGACTTAATACTTATAATAAAAATACCCCTCTCCTGAGTCAGGCACAGTGCCCCTGGGCAGTTGGGGGAAATTAATAATTACAACTGACGAGATAAATCACAATCCTATTGACAAAATCACTCCCCATCAACTATTATTCTTTCTGTAATGGAAATCCATGCGCTCATATATATTCATATTTAGAAAATTACCAGTACAGATAAGCAAATAAATAGCATTTCCCATCTGCATTTTTCTCAACACAGAACCTCAGTAGATGAAGTAAATCTATAGATTTGTACAATTCTGTTTTCAAATAGATAATATTGAGCCTCACTTAGGCAAGATGCGAAGCTCACTATGCAAATCTGTCAAAATCCCCATTGCTCCAATCCCTTCAATGCTGACGGCAACAGATTTTGCATCAGTTGCGGGCAAAGCAACTTTGGTGAATTACTGCGGAACCGTTATCGTGTTCTGAAATTATTAGGAGAAGGTGGTTTTAGTAGAACCTACGCTACAGAAGATGCTGATAGACTCAATGCACCCTGTGTAATTAAACAATTTTTCCCCCAGGTGCATGGAACATCAGAACGTACTAAAGCCGCAGAATTATTTAAAGAAGAAGCCAAGCGTCTATATGAATTGGGAGAAAATCATATTCAAATCCCTCGACTATTGGCATTTTTTGAACAAGGTTCCAGTTTATATTTAGTTCAAGAATTTATTCAAGGAAAAACCCTCTTAGAAGAGGTACAAAAGCAGCATTTTACAGA
a
s 336-3.336-3 2000 1000 + 6283267 AGATAAAGTTAGGGAAATTTTAACTGAGTTATTACCTGTTCTCGATTTTATTCACACCAATAACGTGATTCATCGAGATATTAAACCAGAAAATATTATTCGCCGAGATACAGACGGTAAGCTCGTATTAATCGATTTTGGTGGAGCCAAACAAGTTACCCAAACAAGTTTAGCCAGACAAGCGACAGTTCTCTATACAATTGGTTATGCTCCCAGTGAGCAAATGGCAGGATTTGCTTGCCATGCCAGTGATATCTATGCTTTAGGTGTGACTTGTGTACGCTTACTTACCAGATGTTTGCCACTCCATGATTCATCAGGGCAAGTTTACGATCGCCTTTACGATCCCATGAGTGCTAAATGGCTATGGCGACAACGCATCCAGGAACAAAACATCACAATCAGTGAAGAATTAGGGCATATCCTGGATAAGATGCTGAAAAATTTAGCCTCAGAAAGATATCAATCTGCTGCGGTAATTATGAGAGAGCTACAAAATAGCAAAACCCTATCCCCTGGAGCAAATTCTTCAGAAACACAACCGATAATTGCAATTCCTCGAAATCCATCGAATTTTACCAACTTAGACACCTTTGAATTTGACGTAGTGACAATTGATAGTCATAGTCAGATAGTCAATCGGAACCGTACTGTTGCCCAATATTATGCAGAAGAATTGGGACATGACGTAATTATGGATATGGTAGCAGTACCTGGTGGCTCTTTTATGATGGGTTCCAAGGAAGATGAAGGAGATGCAGACGAACGTCCCCAACACTTAGTAGCAATTGAGCCTTTCTTCATTGGTAGATATCCTGTTACCCAAGCCCAATGGCGAGCAGTGGCATCTTTGCCCAAAGTCAATCAACCCCTCAACCCCTATCCGTCGAAATTCAAAGGTAGCAATCGCCCAGTAGAAAATGTATCTTGGTACGAAGCAGTGGAATTTTGCGCTCGACTAGAGCAAAAAACCGGGCGCAACTATCGCTTACCCAGTGAG
... etc
###Reformateo las etiquetas de los genomas:
## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade
sed 's/PCC_7716.PCC_7716/PCC_7716/g' evolverCalothrix.PCC_7716.hal.maf | sed 's/336-3.336-3/336-3/g' | sed 's/NIES-3974.NIES-3974/NIES-3974/g' | sed 's/NIES-4071.NIES-4071/NIES-4071/g' | sed 's/NIES-4105.NIES-4105/NIES-4105/g' | sed 's/PCC_6303.PCC_6303/PCC_6303/g' >evolverCalothrix.PCC_7716.hal.maf.sedCopio evolverCallothrix.hal.maf.sed a la carpeta de chromosomes_only y convierto maf a fasta:
## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/fna/only_chromosomes/
cd fna/only_chromosomes/
mafToFastaStitcher --maf ../../evolverCalothrix.PCC_7716.hal.maf.sed --seqs 336-3.fna,NIES-3974.fna,NIES-4071.fna,NIES-4105.fna,PCC_6303.fna,PCC_7716.fna --breakpointPenalty 5 --interstitialSequence 20 --outMfa evolverCalothrix.PCC_7716.hal.maf.sed.fastaConvierto a philyp
## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/PALINDROMES/CGGCGCCG/
cd ../../PALINDROMES/CGGCGCCG/
python3 ../../Fasta2Phylip.py ../../fna/only_chromosomes/evolverCalothrix.PCC_7716.hal.maf.sed.fastaObtengo las coordenadas del palindromo
## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade/PALINDROMES/CGGCGCCG/
python3 ../../AlignmentPalindromeCoords.py evolverCalothrix.PCC_7716.hal.maf.sed.fasta.phy CGGCGCCG PCC_7716El 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:
## /home/lalibelulalo/PIPELINES_2023/ASR_GENOMES/Callothrix_clade
File = "Orthologues_Palindrome_sites.txt"
table = read.table(File,header = TRUE,sep = "\t")
rmarkdown::paged_table(table,options = list(rows.print = 10, cols.print = 4))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/Callothrix_clade/PALINDROMES/CGGCGCCG/
Sites <- read.table("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")}Creo una matriz adicional que tendra las deleciones:
{DELECIONES <- matrix(0,
nrow = 0,
ncol = 4)
colnames(DELECIONES) <- c("ortólogo","sitio","deleciones","species")}Cargo el modelo de evolucion, el metodo de reconstrucción (bayesiano), la filogenia (previamente calculada con orthofinder).
EvolModel = "F81"
Method = "bayes"
Tree = read.tree("../../SpeciesTree_rooted.txt")
Alignment = "evolverCalothrix.PCC_7716.hal.maf.sed.fasta.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 1094 sitios en total. Cargamos el primer sitio el cual corresponde a:
table[1,]## FILE PAL START END
## 1 evolverCalothrix.PCC_7716.hal.maf.sed.fasta.phy CGGCGCCG 1269 1276El 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)
parsimony(Tree, cyanobacterias)## [1] 3241393Filogenia 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)$namesAhora 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]
}
## 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
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 336-3 --------
## 2 NIES-3974 --------
## 3 PCC_6303 --------
## 4 PCC_7716 CGGCGCCG
## 5 NIES-4071 CGGCGCCG
## 6 NIES-4105 CGGCGCCG
## 7 7 CGGCGCCG
## 8 8 CGGCGCCG
## 9 9 CGGCGCCG
## 10 10 CGGCGCCGHago un conteo de las deleciones y las guardo en un vector que contendrá el ortologo, sitio y numero de deleciones:
d = 0
SPP = c()
for (i in 1:length(Mat2[,1])) {
if(Mat2[i,2] == '--------'){
d = d+1
spp = Mat2[i,1]
SPP = c(SPP,spp)
}
}
Spp <- paste(SPP, collapse = ",")
deleciones <-c()
if(d != 0){
SITE = paste0(Sites[ORTH,3],":",Sites[ORTH,4])
deleciones <- c(Sites[ORTH,1],SITE,d,Spp)
}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 TRANSICIONESCreo 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 -------- CGGCGCCG
## 2 -------- CGGCGCCG
## 3 -------- CGGCGCCG
## 4 CGGCGCCG CGGCGCCG
## 5 CGGCGCCG CGGCGCCG
## 6 CGGCGCCG CGGCGCCG
## 7 CGGCGCCG CGGCGCCG
## 8 CGGCGCCG CGGCGCCG
## 9 CGGCGCCG CGGCGCCG
## 10 CGGCGCCG CGGCGCCGy 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 -------- CGGCGCCG 336-3--7
## 2 -------- CGGCGCCG NIES-3974--7
## 3 -------- CGGCGCCG PCC_6303--8
## 4 CGGCGCCG CGGCGCCG PCC_7716--9
## 5 CGGCGCCG CGGCGCCG NIES-4071--10
## 6 CGGCGCCG CGGCGCCG NIES-4105--10
## 7 CGGCGCCG CGGCGCCG 7--7
## 8 CGGCGCCG CGGCGCCG 8--7
## 9 CGGCGCCG CGGCGCCG 9--8
## 10 CGGCGCCG CGGCGCCG 10--9Finalmente repito esto para los 1093 sitios palindromicos restantes.
TABLA DE DELECIONES
File= "DELECIONES.txt"
table = read.table(File,header = TRUE,sep = "\t")
rmarkdown::paged_table(table,options = list(rows.print = 10, cols.print = 4))Ahora las matrices DF y LINKS contienen todas las transiciones de todos los sitios. En total son 10940 transiciones ya que son 10 transiciones por arbol y tenemos 1094 arboles.
length(DF$from)## [1] 10940length(LINKS$from)## [1] 10940length(DELECIONES[,1])## [1] 981Filtrado 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.
library(data.table)
RowCts <- setDT(DF)[,list(Count=as.numeric(.N)),names(DF)] ## CUENTO EL NUMERO DE OCURRENCIAS DE CADA TRANSICIÓN #6867
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 #6867## # A tibble: 1,531 × 3
## from to weight
## <chr> <chr> <int>
## 1 -------- CGGCGCCG 2880
## 2 CGGCGCCG CGGCGCCG 3961
## 3 -------- CGGCGACG 5
## 4 CTGCGACG CGGCGACG 1
## 5 CGGCGCCG CGGCGTCG 2
## 6 CGACGTCG CGACGTCG 2
## 7 CGGCGACG CGGCGACG 16
## 8 CGGCGTCG CGGCGACG 2
## 9 CGACGTCG CGGCGTCG 1
## 10 -------- GAGTCAGG 3
## # ℹ 1,521 more rowsFiltro conteos bajos (<=7) y loops. Al final me quedo con 14 transiciones.
RowCts<-RowCts%>% ##QUITO LOS LOOPS. ES DECIR QUITO AQUELLAS TRANSICIONES QUE VAYAN ASI MISMA
filter(from!=to)
length(RowCts$from) #1144## [1] 1144RowCts<-RowCts%>% ## QUITO AQUELLAS TRANSICIONES CON UN PESO MENOR A 7
filter(weight>=7)
length(RowCts$from) #14## [1] 14Extraigo 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) #28## [1] 28GNMtx <- 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 DUPLICADOSCreo 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 10length(ids) #15## [1] 10length(GNMtx[,1]) #15## [1] 10Creo 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] 10Nodes2<-Nodes2[,c(3,1,2)] ## REORDENO LAS COLUMNAS DE LA MATRIZ
colnames(Nodes2)<-c("id","label","score") ## ETIQUETO LAS COLUMNASCreo 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: 10 nodes and 14 edges
## #
## # A directed acyclic simple graph with 1 component
## #
## # A tibble: 10 × 2
## id label
## <int> <chr>
## 1 1 --------
## 2 2 CGGCGCCG
## 3 3 CGGCGCCA
## 4 4 TGGCGCCG
## 5 5 AGGCGCCG
## 6 6 GGGCGCCG
## # ℹ 4 more rows
## #
## # A tibble: 14 × 3
## from to weight
## <int> <int> <int>
## 1 1 2 2880
## 2 2 3 32
## 3 2 4 32
## # ℹ 11 more rowsroutes_tidy %>%
activate(edges) %>%
arrange(desc(weight))## # A tbl_graph: 10 nodes and 14 edges
## #
## # A directed acyclic simple graph with 1 component
## #
## # A tibble: 14 × 3
## from to weight
## <int> <int> <int>
## 1 1 2 2880
## 2 1 4 69
## 3 1 3 63
## 4 1 10 41
## 5 2 3 32
## 6 2 4 32
## # ℹ 8 more rows
## #
## # A tibble: 10 × 2
## id label
## <int> <chr>
## 1 1 --------
## 2 2 CGGCGCCG
## 3 3 CGGCGCCA
## # ℹ 7 more rowsVisualizo 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
Las siguientes figuras corresponden a cada una de las direcciones.
WGHT = 4336-3–7
ggtree(Tree2) +
geom_tiplab(color='firebrick', offset = .10, hjust=3.0)+
geom_label(aes(label = node),show.legend = TRUE)
## [1] 3023## # A tbl_graph: 8 nodes and 7 edges
## #
## # A rooted tree
## #
## # A tibble: 7 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 1 2 336-3--7 649
## 2 1 3 336-3--7 23
## 3 1 4 336-3--7 21
## 4 1 8 336-3--7 13
## 5 1 5 336-3--7 6
## 6 1 7 336-3--7 6
## # ℹ 1 more row
## #
## # A tibble: 8 × 2
## id label
## <int> <chr>
## 1 1 --------
## 2 2 CGGCGCCG
## 3 3 TGGCGCCG
## # ℹ 5 more rowsNIES-3974–7
ggtree(Tree2) +
geom_tiplab(color='firebrick', offset = .10, hjust=3.0)+
geom_label(aes(label = node),show.legend = TRUE)
## [1] 3023## # A tbl_graph: 8 nodes and 7 edges
## #
## # A rooted tree
## #
## # A tibble: 8 × 2
## id label
## <int> <chr>
## 1 1 --------
## 2 2 CGGCGCCG
## 3 3 TGGCGCCG
## 4 4 CGGCGCCA
## 5 5 AGGCGCCG
## 6 6 CGGCGCCC
## # ℹ 2 more rows
## #
## # A tibble: 7 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 1 2 NIES-3974--7 649
## 2 1 3 NIES-3974--7 24
## 3 1 4 NIES-3974--7 21
## # ℹ 4 more rows## # A tbl_graph: 8 nodes and 7 edges
## #
## # A rooted tree
## #
## # A tibble: 7 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 1 2 NIES-3974--7 649
## 2 1 3 NIES-3974--7 24
## 3 1 4 NIES-3974--7 21
## 4 1 8 NIES-3974--7 13
## 5 1 5 NIES-3974--7 6
## 6 1 7 NIES-3974--7 6
## # ℹ 1 more row
## #
## # A tibble: 8 × 2
## id label
## <int> <chr>
## 1 1 --------
## 2 2 CGGCGCCG
## 3 3 TGGCGCCG
## # ℹ 5 more rowsNIES-4071–10
ggtree(Tree2) +
geom_tiplab(color='firebrick', offset = .10, hjust=3.0)+
geom_label(aes(label = node),show.legend = TRUE)
## [1] 3023## # A tbl_graph: 10 nodes and 10 edges
## #
## # A directed multigraph with 9 components
## #
## # A tibble: 10 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 2 1 NIES-4071--10 473
## 2 1 1 NIES-4071--10 258
## 3 3 3 NIES-4071--10 34
## 4 4 4 NIES-4071--10 27
## 5 9 9 NIES-4071--10 15
## 6 8 8 NIES-4071--10 8
## # ℹ 4 more rows
## #
## # A tibble: 10 × 2
## id label
## <int> <chr>
## 1 1 CGGCGCCG
## 2 2 --------
## 3 3 CGGCGCCA
## # ℹ 7 more rowsNIES-4105–10
ggtree(Tree2) +
geom_tiplab(color='firebrick', offset = .10, hjust=3.0)+
geom_label(aes(label = node),show.legend = TRUE)
## [1] 3023## # A tbl_graph: 10 nodes and 10 edges
## #
## # A directed multigraph with 9 components
## #
## # A tibble: 10 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 2 1 NIES-4105--10 473
## 2 1 1 NIES-4105--10 258
## 3 3 3 NIES-4105--10 34
## 4 4 4 NIES-4105--10 27
## 5 9 9 NIES-4105--10 15
## 6 8 8 NIES-4105--10 8
## # ℹ 4 more rows
## #
## # A tibble: 10 × 2
## id label
## <int> <chr>
## 1 1 CGGCGCCG
## 2 2 --------
## 3 3 CGGCGCCA
## # ℹ 7 more rowsPCC_6303–8
ggtree(Tree2) +
geom_tiplab(color='firebrick', offset = .10, hjust=3.0)+
geom_label(aes(label = node),show.legend = TRUE)
## [1] 3023## # A tbl_graph: 9 nodes and 7 edges
## #
## # A rooted forest with 2 trees
## #
## # A tibble: 7 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 1 3 PCC_6303--8 636
## 2 1 4 PCC_6303--8 22
## 3 1 5 PCC_6303--8 21
## 4 1 9 PCC_6303--8 13
## 5 1 8 PCC_6303--8 6
## 6 1 6 PCC_6303--8 5
## # ℹ 1 more row
## #
## # A tibble: 9 × 2
## id label
## <int> <chr>
## 1 1 --------
## 2 2 CGACGACG
## 3 3 CGGCGCCG
## # ℹ 6 more rows}
PCC_7716–9
ggtree(Tree2) +
geom_tiplab(color='firebrick', offset = .10, hjust=3.0)+
geom_label(aes(label = node),show.legend = TRUE)
## [1] 3023## # A tbl_graph: 8 nodes and 8 edges
## #
## # A directed multigraph with 1 component
## #
## # A tibble: 8 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 1 1 PCC_7716--9 746
## 2 1 2 PCC_7716--9 30
## 3 1 3 PCC_7716--9 30
## 4 1 7 PCC_7716--9 17
## 5 1 6 PCC_7716--9 9
## 6 1 4 PCC_7716--9 6
## # ℹ 2 more rows
## #
## # A tibble: 8 × 2
## id label
## <int> <chr>
## 1 1 CGGCGCCG
## 2 2 CGGCGCCA
## 3 3 TGGCGCCG
## # ℹ 5 more rows10–9
ggtree(Tree2) +
geom_tiplab(color='firebrick', offset = .10, hjust=3.0)+
geom_label(aes(label = node),show.legend = TRUE)
## [1] 3023## # A tbl_graph: 8 nodes and 9 edges
## #
## # A directed multigraph with 7 components
## #
## # A tibble: 9 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 1 1 10--9 735
## 2 2 2 10--9 29
## 3 3 3 10--9 28
## 4 7 7 10--9 15
## 5 6 6 10--9 9
## 6 4 4 10--9 6
## # ℹ 3 more rows
## #
## # A tibble: 8 × 2
## id label
## <int> <chr>
## 1 1 CGGCGCCG
## 2 2 CGGCGCCA
## 3 3 TGGCGCCG
## # ℹ 5 more rows9–8
ggtree(Tree2) +
geom_tiplab(color='firebrick', offset = .10, hjust=3.0)+
geom_label(aes(label = node),show.legend = TRUE)
## [1] 3023## # A tbl_graph: 8 nodes and 9 edges
## #
## # A directed multigraph with 6 components
## #
## # A tibble: 9 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 1 1 9--8 664
## 2 2 2 9--8 23
## 3 3 3 9--8 23
## 4 7 7 9--8 13
## 5 6 6 9--8 8
## 6 1 8 9--8 6
## # ℹ 3 more rows
## #
## # A tibble: 8 × 2
## id label
## <int> <chr>
## 1 1 CGGCGCCG
## 2 2 TGGCGCCG
## 3 3 CGGCGCCA
## # ℹ 5 more rows8–7
ggtree(Tree2) +
geom_tiplab(color='firebrick', offset = .10, hjust=3.0)+
geom_label(aes(label = node),show.legend = TRUE)
## [1] 3023## # A tbl_graph: 11 nodes and 11 edges
## #
## # A directed multigraph with 11 components
## #
## # A tibble: 11 × 2
## id label
## <int> <chr>
## 1 1 CGGCGCCG
## 2 2 CGGCGACG
## 3 3 TGGCGCCG
## 4 4 CGGCGCCA
## 5 5 AGGCGCCG
## 6 6 CGCCGCCG
## # ℹ 5 more rows
## #
## # A tibble: 11 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 1 1 8--7 650
## 2 2 2 8--7 6
## 3 3 3 8--7 23
## # ℹ 8 more rows## # A tbl_graph: 11 nodes and 11 edges
## #
## # A directed multigraph with 11 components
## #
## # A tibble: 11 × 4
## from to direction weight
## <int> <int> <chr> <int>
## 1 1 1 8--7 650
## 2 3 3 8--7 23
## 3 4 4 8--7 21
## 4 9 9 8--7 13
## 5 8 8 8--7 10
## 6 2 2 8--7 6
## # ℹ 5 more rows
## #
## # A tibble: 11 × 2
## id label
## <int> <chr>
## 1 1 CGGCGCCG
## 2 2 CGGCGACG
## 3 3 TGGCGCCG
## # ℹ 8 more rows