Reconstrucción Ancestral de Secuencias Completas - ORTÓLOGOS - TGGCGCCA - PCC_7716 - ROOTED

Este pipeline es para hacer una reconstrucción ancestral de secuencias y crear un grafo con toda la información de las transiciones.

Obtención de datos

En este pipe voy a utilizar un clado de 6 especies de calothrix:

File	Accesion
Calothrix_sp_336-3_336-3.gbff	GCF_000734895.2
Calothrix_sp_NIES-3974_NIES-3974.gbff	GCF_002368395.1
Calothrix_sp_NIES-4071_NIES-4071.gbff	GCF_002368455.1
Calothrix_sp_NIES-4105_NIES-4105.gbff	GCF_002368415.1
Calothrix_sp_PCC_6303_PCC_6303.gbff	GCF_000317435.1
Calothrix_sp_PCC_7716_PCC_7716.gbff	GCF_019977735.1

Ortólogos

Primero obtengo los ortólogos usando get_homologues:

## /home/lalibelulalo/PIPELINES_2023/ASR_ORTHOLOGUES/Callothrix_clade/
get_homologues.pl -d gbff -t 0 -M -n PPN

Obtengo el pangenoma:

## /home/lalibelulalo/PIPELINES_2023/ASR_ORTHOLOGUES/Callothrix_clade/
compare_clusters.pl -o pangenoma -m -d gbff_homologues/CalothrixspNIES-3974NIES-3974_f0_0taxa_algOMCL_e0_

Obtengo una lista de ortólogos (orthologues.list) la cual contenga solo aquellos en los que dicho ortólogo esta en las 6 especies del clado.

## /home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS
awk -F "," '{if($15==1 && $16==1 && $17==1 && $18==1 && $19==1 && $20==1) {print $1}}' ../ASR_ORTHOLOGUES/Callothrix_clade/pangenoma/pangenome_matrix_t0.tr.csv | sed 's/.faa/.fna/g' >orthologues.list

Conteo de palindromos

Hago un conteo de palindromos para saber cuales ortólogos contienen los palindromos de interés. Para eso uso orthologues.list y una lista con los palindromos de interes: pals.list.

Para esto uso parte del código del conteo de palindromos.

Pals.list contiene los palindromos a contar. Estos van en un a sola linea y separados por comas. El numero final corresponde al orden del modelo de markov (0-3) sin embargo es irrelevante para este paso ya que solo queremos saber el numero observado de palindromos.

## /home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS
python3 CountPalsInOrthologues.py ../ASR_ORTHOLOGUES/Callothrix_clade/gbff_homologues/CalothrixspNIES-3974NIES-3974_f0_0taxa_algOMCL_e0_ pals.list orthologues.list 3

Filtrado del conteo

Para empezar escojo la especie con mayor conteo de palindromos segun la filogenia anotada (PCC_7716). Sin embargo tambien me aseguro de que dicha especie sea la que tenga el número de sitios:

## /home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS/PALINDROMES/TGGCGCCA
Counts <- read.table("../Markov_count_TGGCGCCA_2023-5-7_18hrs43mins_Octanuc_.txt", sep = "\t", header = TRUE)

Spp <- Counts$spp
Spp <- unique(Spp)

i=0
for (n in 1:length(Spp)){ 
  for (m in 1:length(Counts$spp)) {
    if (Spp[n] == Counts[m,2]){
      i=i+Counts[m,5]
    }
  }
  print(paste0(Spp[n]," ",i))
  i=0
}

## [1] "Calothrix_sp_336/3 3"
## [1] "Calothrix_sp_NIES-3974 0"
## [1] "Calothrix_sp_NIES-4071 670"
## [1] "Calothrix_sp_NIES-4105 671"
## [1] "Calothrix_sp_PCC_6303 1"
## [1] "Calothrix_sp_PCC_7716 795"

Aqui podemos ver que efectivamente Calothrix_sp_PCC_7716 es el que tiene mayores sitios (571). Solo quiero aquellos donde Calothrix_sp_PCC_7716 tenga un conteo mayor o igual a 1 ya que en ese sitio es donde voy a reconstruir la secuencia ancestral

## /home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS/PALINDROMES/TGGCGCCA/
awk '{if ($5>=1 && $2=="Calothrix_sp_PCC_7716"){print $1}}' ../Markov_count_TGGCGCCA_2023-5-7_18hrs43mins_Octanuc_.txt | sed 's/.fna//g' >Ortologos_TGGCGCCA_Calothrix_sp_PCC_7716.txt

Ortologos_TGGCGCCA_Calothrix_sp_PCC_7716.txt contiene 667 sitios los cuales estan presentes en las 6 especies y tienen por lo menos una ocurrencia del palindromo TGGCGCCA

###
5355_caspase_family_prote..
5357_hypothetical_protein
5357_hypothetical_protein
5358_OmpA_family_protein
5360_DUF11_domain-contain..
5362_trpA
5367_larC
5373_DUF751_family_protei..
5377_trpE
5382_DegT-DnrJ-EryC1-StrS..
###

Preparación de ortólogos para la alineción

Primero creo una carpeta que contendrá unicamente los ortologos que vamos a usar. Es decir aquellos que estan presentes en todas las especies a analizar y que contienen por lo menos un sitio con palíndromo.

## /home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS/PALINDROMES/TGGCGCCA/
mkdir Orthologues_TGGCGCCA_PCC_7716  ## Creo Carpeta
cd ../../../ASR_ORTHOLOGUES/Callothrix_clade/gbff_homologues/CalothrixspNIES-3974NIES-3974_f0_0taxa_algOMCL_e0_/ ## Entro a la carpeta de homólogos
for word in $(cat ../../../../CODON_MUTATIONS/PALINDROMES/TGGCGCCA/Ortologos_TGGCGCCA_Calothrix_sp_PCC_7716.txt); do cp $word.fna ../../../../CODON_MUTATIONS/PALINDROMES/TGGCGCCA/Orthologues_TGGCGCCA_PCC_7716; done ## copio los ortólogos de la lista a la carpeta

for word in $(cat ../../../../CODON_MUTATIONS/PALINDROMES/TGGCGCCA/Ortologos_TGGCGCCA_Calothrix_sp_PCC_7716.txt); do cp $word.faa ../../../../CODON_MUTATIONS/PALINDROMES/TGGCGCCA/Orthologues_TGGCGCCA_PCC_7716; done

Filtrado de paralogos

voy a crear dos carpetas una que contenga archivos con paralogos y otra que solo tenga ortólogos

cd ../../../../CODON_MUTATIONS/PALINDROMES/TGGCGCCA/
mkdir Only_ORTHOLOGUES
mkdir PARALOGUES

Para filtrar aquellos con paralogos uso el siguiente script de python:

from Bio import SeqIO
import re
import os
import sys
import shutil

SEQUENCES_PATH = sys.argv[1] #'Orthologues_GCGATCGC_336_3/' ## Path de la carpeta con los ortólogos.
#output_file = 'Orthologues_Palindrome_sites.txt' ## Nombre del archivo de salida
#output = open (output_file, 'w') ## Abrimos el archivo de salida

if SEQUENCES_PATH.endswith("/"): ## Esta parte la pongo por si corro este codigo en la terminal
    SequencesPath = re.sub('/', '', SEQUENCES_PATH) ## De este modo evito errores en el argumento path 
    SequencesDir = str(SEQUENCES_PATH)
else:
    SequencesPath = SEQUENCES_PATH
    SequencesDir = str("".join ([SEQUENCES_PATH,'/']))
Orthologues = [x for x in os.listdir(SequencesDir) if x.endswith(".fna")] ## creo un arreglo con todos los ortólogis de la carpeta

for Orthologue in Orthologues:
    i = 0
    FNA = str("".join ([SequencesDir,Orthologue])) ## creo el path completo con el path de la carpeta de ortólogos y el path del ortólogo
    FAA = re.sub('.fna', '.faa', FNA)
    OrthologueA = re.sub('.fna', '.faa', Orthologue)
    for record in SeqIO.parse(open(FNA),'fasta'):
        i += 1
    #print ('{}\t{}'.format(Orthologue,i))
    if i == 6:
        shutil.copyfile(FNA,str("".join (['Only_ORTHOLOGUES/',Orthologue]))) ## si hasy unicamente 6 entradas entonces no hay paralogos
        shutil.copyfile(FAA,str("".join (['Only_ORTHOLOGUES/',OrthologueA])))
    else:
        shutil.copyfile(FNA,str("".join (['PARALOGUES/',Orthologue]))) ## Si hay mas de 6 entradas hay parálogos
        shutil.copyfile(FAA,str("".join (['PARALOGUES/',OrthologueA])))

python3 ../../FiltradoParalogos.py Orthologues_TGGCGCCA_PCC_7716/

Esto me dio como resultado 1044 archivos sin paralogos y 194 archivos con paralogos los cuals omitiré por ahora.

Ahora edito el archivo fasta para que el encabezado incluya unicamente la especie

## /home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS/PALINDROMES/TGGCGCCA/Only_ORTHOLOGUES/
cd Only_ORTHOLOGUES/
for f in *.fna; do awk -F "|" '{if (NR%2){print ">"$2}else{print $1}}' $f | sed 's/ /_/g' | sed 's/\.//g' | sed 's/\[//g' | sed 's/\]//g' | sed 's/\//-/g' | sed 's/Calothrix_sp_//g' >$f.awk1;done
for f in *.faa; do awk -F "|" '{if (NR%2){print ">"$2}else{print $1}}' $f | sed 's/ /_/g' | sed 's/\.//g' | sed 's/\[//g' | sed 's/\]//g' | sed 's/\//-/g' | sed 's/Calothrix_sp_//g' >$f.awk1;done

Alineación Múltiple

Hago la alineación multiple con MAFFT

## /home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS/PALINDROMES/TGGCGCCA/Only_ORTHOLOGUES/
for f in *faa.awk1; do mafft $f >$f.mafft;done
ls *.fna | sed 's/.fna//g' >../only.orthologues.txt

for f in $(cat ../only.orthologues.txt); do perl ../../../pal2nal.pl $f.faa.awk1.mafft $f.fna.awk1 -output fasta >$f.codon.alg.fasta;done

Cambio de fasta a PHYLIP

## /home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS/PALINDROMES/TGGCGCCA/
cd ..
python3 ../../Fasta2Phylip.py Only_ORTHOLOGUES/

Obtengo las coordenadas del palindromo en la alineación usando el siguiente codigo de python:

from Bio import SeqIO
import re
import os
import sys

SEQUENCES_PATH = sys.argv[1] #'RECONSTRUCCIONES/' #'Only_ORTHOLOGUES/' #
pattern = sys.argv[2] #'GCGATCGC' #
SPP = sys.argv[3]
OUT = sys.argv[4]
#SEQUENCES_PATH = 'gbff/Only_ORTHOLOGUES/' ## Path de la carpeta con los ortólogos.

output_file = '' ## Nombre del archivo de salida
output_file = '.'.join(['Orthologues_Palindrome_sites.AllFrames',OUT,SPP,pattern,'txt'])
output = open (output_file, 'w') ## Abrimos el archivo de salida
output.write('FILE\tSpp\tSTART\tEND\tReadingFrame\tOrthLength\tPAL\tAA\n')

output_fileM1 = '.'.join(['Orthologues_Palindrome_sites.M1',OUT,SPP,pattern,'txt'])
outputM1 = open (output_fileM1, 'w') ## Abrimos el archivo de salida
outputM1.write('FILE\tSpp\tSTART\tEND\tReadingFrame\tOrthLength\tPAL\tAA\n')

output_fileM2 = '.'.join(['Orthologues_Palindrome_sites.M2',OUT,SPP,pattern,'txt'])
outputM2 = open (output_fileM2, 'w') ## Abrimos el archivo de salida
outputM2.write('FILE\tSpp\tSTART\tEND\tReadingFrame\tOrthLength\tPAL\tAA\n')

output_fileM3 = '.'.join(['Orthologues_Palindrome_sites.M3',OUT,SPP,pattern,'txt'])
outputM3 = open (output_fileM3, 'w') ## Abrimos el archivo de salida
outputM3.write('FILE\tSpp\tSTART\tEND\tReadingFrame\tOrthLength\tPAL\tAA\n')

if SEQUENCES_PATH.endswith("/"): ## Esta parte la pongo por si corro este codigo en la terminal
    SequencesPath = re.sub('/', '', SEQUENCES_PATH) ## De este modo evito errores en el argumento path 
    SequencesDir = str(SEQUENCES_PATH)
else:
    SequencesPath = SEQUENCES_PATH
    SequencesDir = str("".join ([SEQUENCES_PATH,'/']))

def translate(seq):
    table = {
        'ATA':'I', 'ATC':'I', 'ATT':'I', 'ATG':'M',
        'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACT':'T',
        'AAC':'N', 'AAT':'N', 'AAA':'K', 'AAG':'K',
        'AGC':'S', 'AGT':'S', 'AGA':'R', 'AGG':'R',
        'CTA':'L', 'CTC':'L', 'CTG':'L', 'CTT':'L',
        'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCT':'P',
        'CAC':'H', 'CAT':'H', 'CAA':'Q', 'CAG':'Q',
        'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGT':'R',
        'GTA':'V', 'GTC':'V', 'GTG':'V', 'GTT':'V',
        'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCT':'A',
        'GAC':'D', 'GAT':'D', 'GAA':'E', 'GAG':'E',
        'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGT':'G',
        'TCA':'S', 'TCC':'S', 'TCG':'S', 'TCT':'S',
        'TTC':'F', 'TTT':'F', 'TTA':'L', 'TTG':'L',
        'TAC':'Y', 'TAT':'Y', 'TAA':'_', 'TAG':'_',
        'TGC':'C', 'TGT':'C', 'TGA':'_', 'TGG':'W',
        '---':'-',
    }
    protein =""
    if len(seq)%3 == 0:
        for i in range(0, len(seq), 3):
            codon = seq[i:i + 3]
            protein+= table[codon]
    return protein

    
Orthologues = [x for x in os.listdir(SequencesDir) if x.endswith(".phy")] ## creo un arreglo con todos los ortólogis de la carpeta
#pattern = '[cC][-]*[gG][-]*[gG][-]*[cC][-]*[gG][-]*[cC][-]*[cC][-]*[gG]'
pattern = re.sub('G', '[gG][-]*', pattern)
pattern = re.sub('C', '[cC][-]*', pattern)
pattern = re.sub('T', '[tT][-]*', pattern)
pattern = re.sub('A', '[aA][-]*', pattern)
pattern = pattern[:-4]


k=0
RF1=0
RF2=0
RF3=0
for orthologue in Orthologues:
    file = re.sub('\.fna\.awk1\.mafft\.phy', '', orthologue)
    FNA = str("".join ([SequencesDir,orthologue]))
    spps = [str(record.description) for record in SeqIO.parse(open(FNA),'phylip')]
    sequencesDict = {str(record.description):str(record.seq) for record in SeqIO.parse(open(FNA),'phylip')}
    Sites = [match.span()[0] for match in re.finditer(pattern, sequencesDict[SPP])]
    EndSites = [match.span()[1] for match in re.finditer(pattern, sequencesDict[SPP])]
    j=0
    for site in Sites:
        k += 1
        for spp in spps:    
            end = EndSites[j]
            kmer = sequencesDict[spp][site:end]
            OrthLength = len(sequencesDict[spp])
            start = site
            mod1 = (start+1+2) % 3
            mod2 = (start+1+1) % 3
            mod3 = (start+1+0) % 3
            if mod1 == 0:
                RF = 1
                RF1+=1
            elif mod2 == 0:
                RF = 2
                RF2+=1
            elif mod3 == 0:
                RF = 3
                RF3+=1
                
            if len(kmer)==8 and RF ==1:
                AAstart = int(((start+3)/3)-1)
                AAend = AAstart+3
                AAseq = sequencesDict[spp]
                AA = translate(AAseq)[AAstart:AAend]
                codon1 = sequencesDict[spp][site:site+3]
                codon2 = sequencesDict[spp][site+3:site+3+3]
                codon3 = sequencesDict[spp][site+3+3:site+3+3+3]
                word = " ".join ([codon1,codon2,codon3])
                print ('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(file,spp,start+1,end,RF,OrthLength,word,AA))
                output.write('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(file,spp,start+1,end,RF,OrthLength,word,AA))
                outputM1.write('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(file,spp,start+1,end,RF,OrthLength,word,AA))
            elif len(kmer)==8 and RF ==2:
                AAstart = int(((start+2)/3)-1)
                AAend = AAstart+3
                AAseq = sequencesDict[spp]
                AA = translate(AAseq)[AAstart:AAend]
                codon1 = sequencesDict[spp][site-1:site-1+3]
                codon2 = sequencesDict[spp][site-1+3:site-1+3+3]
                codon3 = sequencesDict[spp][site-1+3+3:site-1+3+3+3]
                word = " ".join ([codon1,codon2,codon3])
                print ('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(file,spp,start+1,end,RF,OrthLength,word,AA))
                output.write('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(file,spp,start+1,end,RF,OrthLength,word,AA))
                outputM2.write('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(file,spp,start+1,end,RF,OrthLength,word,AA))
            elif len(kmer)==8 and RF ==3:
                AAstart = int(((start+1)/3)-1)
                AAend = AAstart+4
                AAseq = sequencesDict[spp]
                AA = translate(AAseq)[AAstart:AAend]
                codon1 = sequencesDict[spp][site-2:site-2+3]
                codon2 = sequencesDict[spp][site-2+3:site-2+3+3]
                codon3 = sequencesDict[spp][site-2+3+3:site-2+3+3+3]
                codon4 = sequencesDict[spp][site-2+3+3+3:site-2+3+3+3+3]
                word = " ".join ([codon1,codon2,codon3,codon4])
                print ('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(file,spp,start+1,end,RF,OrthLength,word,AA))
                output.write('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(file,spp,start+1,end,RF,OrthLength,word,AA))
                outputM3.write('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n'.format(file,spp,start+1,end,RF,OrthLength,word,AA))
        j += 1
        print ('_________________________________________________________________________\n')
        #k = k + len(Sites)
output.close()
outputM1.close()
outputM2.close()
outputM3.close()
print ("TOTAL: {} sitios".format(k))

print ("Hay {} sitios en el marco de lectura 1.".format(RF1/len(spps)))
print ("Hay {} sitios en el marco de lectura 2.".format(RF2/len(spps)))
print ("Hay {} sitios en el marco de lectura 3.".format(RF3/len(spps)))

python3 ../../AlignmentPalindromeCoordsAA.py Only_ORTHOLOGUES/ TGGCGCCA PCC_7716 FIRST

El ultimo script tiene como salida 4 archivos: - Orthologues_Palindrome_sites.M1.FIRST.PCC_7716.TGGCGCCA.txt

• Orthologues_Palindrome_sites.M2.FIRST.PCC_7716.TGGCGCCA.txt

• Orthologues_Palindrome_sites.M3.FIRST.PCC_7716.TGGCGCCA.txt

• Orthologues_Palindrome_sites.AllFrames.FIRST.PCC_7716.TGGCGCCA.txt

Los primero 3 archivos corresponden a los 3 marcos de lectura respectivamente y el ultimo corresponde a los 3 marcos en uno solo archivo.

###
FILE    Spp START   END ReadingFrame    OrthLength  PAL AA
6222_NAD-P--dependent_oxi...codon.alg.fasta.phy 336-3   344 351 2   870 GAA GCA CCT EAP
6222_NAD-P--dependent_oxi...codon.alg.fasta.phy NIES-3974   344 351 2   870 GAA GCA CCT EAP
6222_NAD-P--dependent_oxi...codon.alg.fasta.phy NIES-4071   344 351 2   870 GAG GCG CCT EAP
6222_NAD-P--dependent_oxi...codon.alg.fasta.phy NIES-4105   344 351 2   870 GAG GCG CCT EAP
6222_NAD-P--dependent_oxi...codon.alg.fasta.phy PCC_6303    344 351 2   870 GAA GCC CCC EAP
6222_NAD-P--dependent_oxi...codon.alg.fasta.phy PCC_7716    344 351 2   870 GAG GCG CCT EAP
9935_cobyrinate_a-c-diami...codon.alg.fasta.phy 336-3   1256    1263    2   1461    CAT TCT CCG HSP
9935_cobyrinate_a-c-diami...codon.alg.fasta.phy NIES-3974   1256    1263    2   1461    TCC ACA CCT STP
9935_cobyrinate_a-c-diami...codon.alg.fasta.phy NIES-4071   1256    1263    2   1461    AAG GCG CCT KAP
... etc
###

RECONSTRUCCION DE SECUENCIA COMPLETA

En esta parte reconstruyo las secuencias ancestrales para cada nodo en la filogenia. Despues creo un archivo fasta con todas las secuencias de cada nodo alineadas.

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

#setwd("/home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS/PALINDROMES/TGGCGCCA/")
Sites <- read.table("Orthologues_Palindrome_sites.AllFrames.FIRST.PCC_7716.TGGCGCCA.txt", sep = "\t", header = TRUE)
PATH = 'Only_ORTHOLOGUES/'
Sites<-Sites%>% 
  filter(Spp=='PCC_7716')

EvolModel = "F81"
Method = "bayes"
Tree = read.tree("../../SpeciesTree_rooted.txt")

v = 0
#ORTH = 1
for (ORTH in 1:length(Sites[,1])){
  ## CARGO EL ALINEAMIENTO
  Alignment =paste0(PATH,Sites[ORTH,1])
  cyanobacterias <- read.phyDat(Alignment, format = "interleaved")
  cyanobacterias <- subset(cyanobacterias, site.pattern = FALSE) 
  parsimony(Tree, cyanobacterias)
  
  ## CARGO EL MODELO
  fit <- pml(Tree, cyanobacterias,model=EvolModel)
  #fit <- optim.pml(fit, model=EvolModel, control = pml.control(trace=0),)
  
  ## CARGO La NUEVA FILOGENIA SIN RAIZ
  #Tree2 <- fit[["tree"]]
  
  ## HAGO LA RECONSTRUCCIÓN
  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
  }
  
  ## LEO EL SITIO DE INICIO Y FINAL DE LA SECUENCIA
  ## ES DECIR EL INTERVALO COMPLETO DEL ORTOLOGO
  PalInterval0 = 1:Sites[ORTH,6]
  
  ## Obtengo las posiciones para la posicion PalInterval
  PalindromeNucleotidePositions0 <- attributes(cyanobacterias)$index[PalInterval0]
  
  ## Nombre de los nodos
  Nodos0 <- attributes(Reconstruction)$names
  
  ################## EN ESTA PARTE ARMAREMOS LOS PALINDROMOS DE LAS PUNTAS ###################
  #a.1 <- as_tibble(Tree2)
  a.1 <- as_tibble(Tree)
  ## Extraigo el Boostrap por si se ocupa
  NA.NODE = length(cyanobacterias)+1
  for (i in NA.NODE:length(Nodos0)){ ## ESTO DEPENDE DE LA ESTRUCTURA DEL ARBOL
    a.1$label[i] <- a.1$node[i] 
  }
  ## Agrego los datos de los palindromos al arbol
  b.1 <- tibble::tibble(label = Nodos0)
  c.1 <- full_join(a.1, b.1, by = 'label')
  NodeLabels <- c.1$label
  
  #TipLabels <- Tree2$tip.label
  TipLabels <- Tree$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
  ScoresMtx0 <- matrix(0,
                       nrow = length(Nodos0),
                       ncol = length(PalInterval0))
  colnames(ScoresMtx0) <- PalInterval0
  rownames(ScoresMtx0) <- Nodos0
  ## Relleno la matriz con los scores
  j0 = 1
  i0 = 1
  k0 = 0
  for (Position in PalindromeNucleotidePositions0){
    for (Nodo0 in TipLabels){
      PositionScores0 <- Reconstruction[[Nodo0]][Position,] ## Extraigo los scores para el nodo "NODO" en la posicion "POSITION" de la secuencia
      if (length(unique(PositionScores0)) == 1){
        k0=k0+1
        WinnerNucleotide0 = "-"
        ScoresMtx0[i0,j0] <- WinnerNucleotide0
        i0=i0+1 ## contador para cada Nodo
      }else{
        PositionWinnerScore0 <- max(PositionScores0) ## Extraigo el valor mas grande. Este es el valor de la probabilidad de que el nucleotido ancestral sea ese
        WinnerNucleotide0 <- which(PositionScores0 == PositionWinnerScore0, arr.ind = T) ## pregunto a que numero (Nucleótido) corresponde ese valor
        ScoresMtx0[i0,j0] <- WinnerNucleotide0 ## Guardo dicho valor en la matriz de resultados
        i0=i0+1 ## contador para cada Nodo
      }
    }
    j0=j0+1 ## contador para cada posición
    i0=1 ## Reiniciamoa el contador de Nodo para la siguiente posicion
  }
  ################## EN ESTA PARTE ARMAREMOS LOS PALINDROMOS DE LOS NODOS ###################
  ## Relleno la matriz con los scores
  j0 = 1
  i0 = length(TipLabels)+1
  k0 = 0
  #(length(TipLabels)+1):length(NodeLabels)
  for (Position0 in PalindromeNucleotidePositions0){
    for (Nodo0 in NodeLabels2){
      PositionScores0 <- Reconstruction[[Nodo0]][Position0,] ## Extraigo los scores para el nodo "NODO" en la posicion "POSITION" de la secuencia
      if (length(unique(PositionScores0)) == 1){
        k0=k0+1
        WinnerNucleotide0 = "X"
        ScoresMtx[i0,j0] <- WinnerNucleotide
        i0=i0+1 ## contador para cada Nodo
      }else{
        PositionWinnerScore0 <- max(PositionScores0) ## Extraigo el valor mas grande. Este es el valor de la probabilidad de que el nucleotido ancestral sea ese
        WinnerNucleotide0 <- which(PositionScores0 == PositionWinnerScore0, arr.ind = T) ## pregunto a que numero (Nucleótido) corresponde ese valor
        ScoresMtx0[i0,j0] <- WinnerNucleotide0 ## Guardo dicho valor en la matriz de resultados
        i0=i0+1 ## contador para cada Nodo
      }
    }
    j0=j0+1 ## contador para cada posición
    i0=length(TipLabels)+1 ## Reiniciamoa el contador de Nodo para la siguiente posicion
  }
  ################## EN ESTA PARTE ARMAREMOS LOS PALINDROMOS DE LOS NODOS ###################
  
  ## EN ESTA PARTE ARMAREMOS EL PALINDROMO PARA CADA NODO
  
  ## Creo un matriz de 1x14 para guardar la etiqueta (palindromo) de cada nodo
  PalsMtx0 <- matrix(0,
                     nrow = length(Nodos0),
                     ncol = 1)
  colnames(PalsMtx0) <- c("Sequence")
  rownames(PalsMtx0) <- Nodos0
  
  ## 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 Nodos0){
    j=j+1
    Palindrome <- as.character(ScoresMtx0[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 = "")
    }
    PalsMtx0[j,1] <- pal
  }
  ## CREO EL ARCHIVO FASTA CON LA RECONSTRUCCION CONMPLETA DE SECUENCIAS
  FASTA <-c()
  for (i in 1:length(PalsMtx0[,1])){
    header <- paste0('>',rownames(PalsMtx0)[i])
    seq <- PalsMtx0[i,1]
    FASTA <- rbind(FASTA,header)
    FASTA <- rbind(FASTA,seq)
  }
  
  ## ESCRIBO EL ARCHIVO FASTA
  write.table(FASTA,file=paste0(PATH,Sites[ORTH,1],'.RECONSTRUCTION.fasta'),sep="\t",row.names = FALSE,col.names = FALSE,quote = FALSE)
  #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #print(paste0("Sitio ",ORTH, " de ",length(Sites$FILE),"."))
}

Ahora creo una nueva carpeta llamada RECONSTRUCCIONES y muevo las reconstrucciones ahi. Posteriormente cambio el formato fasta a phylip:

## /home/lalibelulalo/PIPELINES_2023/CODON_MUTATIONS/PALINDROMES/TCGCGCGA/
mkdir RECONSTRUCCIONES
mv Only_ORTHOLOGUES/*.RECONSTRUCTION.fasta RECONSTRUCCIONES
python3 ../../Fasta2Phylip.py RECONSTRUCCIONES/

Ahora vuelvo a correr el script de conteos:

python3 ../../AlignmentPalindromeCoordsAA.py RECONSTRUCCIONES/ TGGCGCCA PCC_7716 SECOND

ESTE ES EL RESULTADO:

File= "Orthologues_Palindrome_sites.AllFrames.SECOND.PCC_7716.TGGCGCCA.txt"
table = read.table(File,header = TRUE,row.names = NULL,sep = "\t")
#table2 = table[ , c(1,2,5,7,8)]  
rmarkdown::paged_table(table,options = list(rows.print = 20, cols.print = 5))

FRECUENCIA DE LOS PEPTIDOS

Primero leo el archivo de resultados del segundo conteo y extraigo solo las columnas de ESPECIE (Spp), MARCO DE LECTURA (ReadingFrame) y PEPTIDO (AA):

File= "Orthologues_Palindrome_sites.AllFrames.SECOND.PCC_7716.TGGCGCCA.txt"
table = read.table(File,header = TRUE,sep = "\t",row.names = NULL)
tableMA = table[ , c(2,5,8)]
tableMA = as.data.frame(tableMA)

MARCO DE LECTURA 1

Filtro el marco de lectura 1 y cuento los peptidos por especie:

tableM1<-tableMA%>%
  filter(ReadingFrame==1)
AACountsM1 <- data.table::setDT(tableM1)[,list(Counts=as.numeric(.N)),names(tableM1)]
AACountsM1 <- AACountsM1%>% 
  filter(Counts>=2)

grafico:

ggplot(AACountsM1, aes(x = reorder(Spp, -Counts), y = Counts, fill = AA)) + 
  geom_bar(position="dodge",stat = "identity") +
  ggtitle("Marco de lectura 1") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Peptido")+
  geom_text(aes(label = Counts, group = AA), position = position_dodge(0.9), vjust = -0.3, size = 2.0)

Figura 1. Grafico de densidad que muestra la frecuencia de los peptidos en el Marco de Lectura 1.

MARCO DE LECTURA 2

tableM2<-tableMA%>%
  filter(ReadingFrame==2)
AACountsM2 <- data.table::setDT(tableM2)[,list(Counts=as.numeric(.N)),names(tableM2)]
AACountsM2 <- AACountsM2%>% 
  filter(Counts>=5)
ggplot(AACountsM2, aes(x = reorder(Spp, -Counts), y = Counts, fill = AA)) + 
  geom_bar(position="dodge",stat = "identity") +
  ggtitle("Marco de lectura 2") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Peptido")+
  geom_text(aes(label = Counts, group = AA), position = position_dodge(0.9), vjust = -0.3, size = 2.5)

Figura 2. Grafico de densidad que muestra la frecuencia de los peptidos en el Marco de Lectura 2.

MARCO DE LECTURA 3

tableM3<-tableMA%>%
  filter(ReadingFrame==3)
AACountsM3 <- data.table::setDT(tableM3)[,list(Counts=as.numeric(.N)),names(tableM3)]
AACountsM3 <- AACountsM3%>% 
  filter(Counts>=8)
ggplot(AACountsM3, aes(x = reorder(Spp, -Counts), y = Counts, fill = AA)) + 
  geom_bar(position="dodge",stat = "identity") +
  ggtitle("Marco de lectura 3") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Peptido")+
  geom_text(aes(label = Counts, group = AA), position = position_dodge(0.9), vjust = -0.3, size = 3.0)

Figura 3. Grafico de densidad que muestra la frecuencia de los peptidos en el Marco de Lectura 3.

MUTACIONES POR CODON Y MARCO DE LECTURA

Primero extraigo el nombre de los nodos parentales del de la filogenia:

x <- as_tibble(Tree)
UNROOTED = 7:10
ROOTED = 7:11
for (i in ROOTED){
  x$label[i] <- x$node[i]
}
FromTo = as.data.frame(x[ , c(1,4)])
colnames(FromTo) <- c("from","to")

Marco de Lectura 1

Primero filtro el Marco de lectura 1

tableRF1<-table%>%
  filter(ReadingFrame==1)

Agrego 3 columnas que contendrán: el nodo parental, secuencia parental y AA parental.

tableRF1$parent = NA
tableRF1$parentPAL = NA
tableRF1$parentAA = NA

Quito espacio de cada codon

for (i in 1:length(tableRF1[,1])){
  #print(paste0('Linea ',i,' de ',length(tableRF1[,1]),'.'))
  tableRF1[i,7] <- gsub(' ','',tableRF1[i,7])
}

Agrego ancestros (relleno las nuevas columnas)

for (i in 1:length(tableRF1[,1])){
  #print(paste0('Linea ',i,' de ',length(tableRF1[,1]),'.'))
  for (j in 1:length(FromTo[,2])){
    if (tableRF1[i,2]=='336-3'|tableRF1[i,2]=='NIES-3974'|tableRF1[i,2]=='PCC_6303'|tableRF1[i,2]=='PCC_7716'|tableRF1[i,2]=='NIES-4071'){
      tableRF1[i,9] <- tableRF1[i+6,2]
      tableRF1[i,10] <- tableRF1[i+6,7]
      tableRF1[i,11] <- tableRF1[i+6,8]
    }
    else if (tableRF1[i,2]=='NIES-4105'){
      tableRF1[i,9] <- tableRF1[i+5,2]
      tableRF1[i,10] <- tableRF1[i+5,7]
      tableRF1[i,11] <- tableRF1[i+5,8]
    }
    else if (tableRF1[i,2]=='8'|tableRF1[i,2]=='9'|tableRF1[i,2]=='10'|tableRF1[i,2]=='11'){
      tableRF1[i,9] <- tableRF1[i-1,2]
      tableRF1[i,10] <- tableRF1[i-1,7]
      tableRF1[i,11] <- tableRF1[i-1,8]
    }
  }
}

Quito el nodo 7 ya que es el nodo mas antiguo y no tiene datos “ancestrales”.

tableRF1<-tableRF1%>%
  filter(Spp!='7')

Creo el archivo con la nueva informacion (RF1.RECONSTRUCTION.rooted.parents.txt):

write.table(tableRF1,file='RF1.RECONSTRUCTION.rooted.parents.txt',sep="\t",row.names = FALSE,col.names = FALSE,quote = FALSE)

Ahora voy agregar mas información a la tabla con el siguiente script. Este script agrega 9 columnas mas para los marcos de lectura 1 y 2. Para el marco de lectura 3 agrega 11. Las columnas agregadas son: - C1Mutations (Mutaciones DEL SITIO en el CODON 1)

• C2Mutations (Mutaciones DEL SITIO en el CODON 2)

• C3Mutations (Mutaciones DEL SITIO en el CODON 3)

• TotalMutations (Mutaciones totales DEL SITIO)

• BL1 (Score (BLOSUM62) de la sustitucion del AA del Codon 1)

• BL2 (Score (BLOSUM62) de la sustitucion del AA del Codon 2)

• BL3 (Score (BLOSUM62) de la sustitucion del AA del Codon 3)

• AASubs (Código de 3 (Marcos de lectura 1 y 2) o 4 (Marco de lectura 3) letras ):

  • S (Mismo AA)
  • s (Diferente AA, score BLOSUM62>=0)
  • N (Diferente AA)
  • D (Deleción)

• SType (Dependiendo de la combinacion de letras (S,s,N,D)):

  • Conservative (la secuencia de AA cambió pero tiene similitud de acuerdo al score de BLOSUM62)
  • ConservativeNoSiteMut (la secuencia de AA cambió pero tiene similitud de acuerdo al score de BLOSUM62. Sin embargo, el sitio no sufrió mutaciones)
  • Deletion (La secuencia de AA tiene sufrio 1 o mas deleciones)
  • NoMutation (La secuencia de AA no sufrio mutaciones)
  • NoSynonym (La secuencia de AA cambió)
  • NoSynonymNoSiteMut (La secuencia de AA cambió. Sin embargo, el sitio no sufrió mutaciones.)
  • Synonym (El sitio sufrió mutaciones. Sin embargo, la secuencia de AA no cambió.)

Corro el script

python3 ../CodonMutations.py RF1.RECONSTRUCTION.rooted.parents.txt 1 codon_mutations_RF1.rooted.txt TGGCGCCA

Como resultado del script anterior obtengo el archivo codon_mutations_RF1.rooted.txt Cargo dicho archivo en R:

FileRF1 = "codon_mutations_RF1.rooted.txt"
tableRF1 = read.table(FileRF1,header = TRUE,sep = "\t",row.names = NULL)

Extraigo unicamente la columna de especie (Spp) y SType:

tableM1.SType = tableRF1[ , c(2,20)]

Cuento los tipos de sustitución para cada Nodo:

tableM1.SType <- data.table::setDT(tableM1.SType)[,list(Counts=as.numeric(.N)),names(tableM1.SType)]

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

Grafico:

ggplot(tableM1.SType, aes(x = reorder(Spp, -Counts), y = Counts, fill = SType)) + 
  geom_bar(position="dodge",stat = "identity") +
  ggtitle("Marco de lectura 1") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Type")+
  geom_text(aes(label = Counts, group = SType), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 4. Grafico de densidad que muestra la frecuencia del tipo de sustituciónes en el Marco de Lectura 1.

Marco de Lectura 2

python3 ../CodonMutations.py RF2.RECONSTRUCTION.rooted.parents.txt 2 codon_mutations_RF2.rooted.txt TGGCGCCA

FileRF2 = "codon_mutations_RF2.rooted.txt"
tableRF2 = read.table(FileRF2,header = TRUE,sep = "\t",row.names = NULL)

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

tableM2.SType = tableRF2[ , c(2,20)]
tableM2.SType <- data.table::setDT(tableM2.SType)[,list(Counts=as.numeric(.N)),names(tableM2.SType)]
ggplot(tableM2.SType, aes(x = reorder(Spp, -Counts), y = Counts, fill = SType)) + 
  geom_bar(position="dodge",stat = "identity") +
  ggtitle("Marco de lectura 2") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Type")+
  geom_text(aes(label = Counts, group = SType), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 5. Grafico de densidad que muestra la frecuencia del tipo de sustituciónes en el Marco de Lectura 2.

Marco de lectura 3

python3 ../CodonMutations.py RF3.RECONSTRUCTION.rooted.parents.txt 3 codon_mutations_RF3.rooted.txt TGGCGCCA

FileRF3 = "codon_mutations_RF3.rooted.txt"
tableRF3 = read.table(FileRF3,header = TRUE,sep = "\t",row.names = NULL)

tableM3.SType = tableRF3[ , c(2,22)]
tableM3.SType <- data.table::setDT(tableM3.SType)[,list(Counts=as.numeric(.N)),names(tableM3.SType)]

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

ggplot(tableM3.SType, aes(x = reorder(Spp, -Counts), y = Counts, fill = SType)) + 
  geom_bar(position="dodge",stat = "identity") +
  ggtitle("Marco de lectura 3") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Type")+
  geom_text(aes(label = Counts, group = SType), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 6. Grafico de densidad que muestra la frecuencia del tipo de sustituciones en el Marco de Lectura 3.

Solo sitios TGGCGCCA

Agrego dos columnas a la tabla:

• AncestorSite (sequencia parental)

• ActualSite (sequencia actual)

• AncestorType (SITE o NoSITE)

• HIPMutPos1-8 (indica con un 1 o 0 si hay una mutacion en las posiciones 1-8)

• HIPMutTypePos1-8 (indica si la mutacion en las posiciones 1-8 son transiciones o transversiones)

Tipos de sustitución para sitios con secuencias TGGCGCCA en sus ancestros

tableM1.SType = tableRF1[ , c(2,20)]
tableM1.SType <- data.table::setDT(tableM1.SType)[,list(Counts=as.numeric(.N)),names(tableM1.SType)]
ggplot(tableM1.SType, aes(x = reorder(Spp, -Counts), y = Counts, fill = SType)) + 
  geom_bar(position="dodge",stat = "identity") +
  ggtitle("Marco de lectura 1 - Nodos con TGGCGCCA en el nodo parental") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Type")+
  geom_text(aes(label = Counts, group = SType), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 7. En esta figura se muestra el tipo de sustitucion que tuvieron aquellos sitios que en su nodo parental eran secuencias TGGCGCCA en el Marco de Lectura 1. Por ejemplo, para el nodo 8. La mayoria de los sitios TGGCGCCA del nodo parental (Nodo 7) se conservaron en el Nodo 8, mientras que una parte de esos sitios (169) tuvieron sustituciones sinonimas y 210 tuvieron sustituciones No sinonimas.

Cargo la tabla

FileRF1 = "codon_mutations_RF1.rooted.txt"
tableRF1 = read.table(FileRF1,header = TRUE,sep = "\t",row.names = NULL)

Filtro dejando solo aquellos casos en los que TGGCGCCA fue el palíndromo parental

tableRF1<-tableRF1%>%
  filter(AncestorType=='SITE')

####Quito los nodos correspondientes a la especie PCC_7716 ya que todos sus sitios son AGGCGCCT

#tableRF1<-tableRF1%>%
#  filter(Spp!='PCC_7716')

Quito los casos que tengan deleciones

tableRF1<-tableRF1%>%
  filter(SType!='Deletion')

Quito los casos que no hayan tenido cambios

#tableRF1<-tableRF1%>%
#  filter(SType!='NoMutation')

Extraigo los nombres de los nodos

Spps <-unique(tableRF1[ , c(2)])

Marco de Lectura 1

TGGCGCCA per site mutations

Creo otra tabla a partir de las columnas con las mutaciones por sitio (columnas 24 a 31)

HIPmut <- tableRF1[ , c(2,25,26,27,28,29,30,31,32)]

Creo una tabla con número de mutaciones en cada nucleotido para cada nodo.

## Creo un nuevo dataframe vacío
df = c()

## Para cada nodo hago la suma de las mutaciones  para cada posicion del octanucleótido y lo agrego al nuevo dataframe
for (i in 1:length(Spps)){
  spp <- Spps[i] # Especie en cuestión
  HIPmut.spp <- HIPmut%>%
    filter(Spp==spp) # Filtro la especie en cuestion de la nueva tabla
  for (j in 1:8){
    # Este nuevo vector "comprime" las mutaciones por sitio de cada nodo en un solo renglón haciendo una suma de todos los renglones
    Node.spp <- summarise(HIPmut.spp,Nodo=unique(spp),NucPos=paste0('Nuc',j),Count=sum(HIPmut.spp[j+1]))
    # Agrego este nuevo renglon al nuevo dataframe
    df <- rbind(df, Node.spp)
  }
}

ggplot(df, aes(x = reorder(Nodo, -Count), y = Count, fill = NucPos)) + 
  geom_bar(position="dodge",stat = "identity") +
  scale_fill_brewer(palette="Set2")+
  ggtitle("Marco de lectura 1 - Nodos con TGGCGCCA en el nodo parental - Mutaciones por nucleotido") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Position")+
  geom_text(aes(label = Count, group = NucPos), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 8. En esta figura se muestra la cantidad de mutaciones por cada nucleótido de TGGCGCCA en el Marco de lectura 1

TCGCGCGA per mutation Type

TR <- tableRF1[ , c(2,33,34,35,36,37,38,39,40)]

df2 = c()
for (i in 1:length(Spps)){
  spp <- Spps[i]
  TR.spp <- TR%>%
  filter(Spp==spp)
  for (j in 1:length(TR.spp[,1])){
    for (k in 1:8){
      vector <- c(TR.spp[j,1],TR.spp[j,k+1])
      df2 <- rbind(df2, vector)
    }
  }
}
df2<- as.data.frame(df2)
df2.SMType <- data.table::setDT(df2)[,list(Counts=as.numeric(.N)),names(df2)]
colnames(df2.SMType) <- c("Nodo","SMType","Counts")

ggplot(df2.SMType, aes(x = reorder(Nodo, -Counts), y = Counts, fill = SMType)) + 
  geom_bar(position="dodge",stat = "identity") +
  scale_fill_brewer(palette="Set3")+
  ggtitle("Marco de lectura 1 - Nodos con TGGCGCCA en el nodo parental - Tipo de sustituciones por nucleotido") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Type")+
  geom_text(aes(label = Counts, group = SMType), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 9. En esta figura se muestra el número de tipos de sustitucion (transición o transversión) para cada nodo en el Marco de lectura 1

Marco de Lectura 2

## Cargo la tabla
FileRF2 = "codon_mutations_RF2.rooted.txt"
tableRF2 = read.table(FileRF2,header = TRUE,sep = "\t",row.names = NULL)


tableRF2<-tableRF2%>%
  filter(AncestorType=='SITE')
#tableRF2<-tableRF2%>%
#  filter(Spp!='PCC_7716')
tableRF2<-tableRF2%>%
  filter(SType!='Deletion')
tableRF2<-tableRF2%>%
  filter(SType!='NoMutation')

Spps <-unique(tableRF2[ , c(2)])

## HIP per site mutations
HIPmutB <- tableRF2[ , c(2,25,26,27,28,29,30,31,32)]
dfB = c()

for (i in 1:length(Spps)){
  spp <- Spps[i] # Especie en cuestión
  HIPmutB.spp <- HIPmutB%>%
    filter(Spp==spp) # Filtro la especie en cuestion de la nueva tabla
  for (j in 1:8){
    # Este nuevo vector "comprime" las mutaciones por sitio de cada nodo en un solo renglón haciendo una suma de todos los renglones
    Node.spp <- summarise(HIPmutB.spp,Nodo=unique(spp),NucPos=paste0('Nuc',j),Count=sum(HIPmutB.spp[j+1]))
    # Agrego este nuevo renglon al nuevo dataframe
    dfB <- rbind(dfB, Node.spp)
  }
}

## Gráfico
ggplot(dfB, aes(x = reorder(Nodo, -Count), y = Count, fill = NucPos)) + 
  geom_bar(position="dodge",stat = "identity") +
  scale_fill_brewer(palette="Set2")+
  ggtitle("Marco de lectura 2 - Nodos con TGGCGCCA en el nodo parental - Mutaciones por nucleotido") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Position")+
  geom_text(aes(label = Count, group = NucPos), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 10. En esta figura se muestra la cantidad de mutaciones por cada nucleótido de TGGCGCCA en el Marco de lectura 2

## HIP per mutation Type
TRB <- tableRF2[ , c(2,33,34,35,36,37,38,39,40)]
df2B = c()
for (i in 1:length(Spps)){
  spp <- Spps[i]
  TRB.spp <- TRB%>%
    filter(Spp==spp)
  for (j in 1:length(TRB.spp[,1])){
    for (k in 1:8){
      vector <- c(TRB.spp[j,1],TRB.spp[j,k+1])
      df2B <- rbind(df2B, vector)
    }
  }
}
df2B<- as.data.frame(df2B)
df2B.SMType <- data.table::setDT(df2B)[,list(Counts=as.numeric(.N)),names(df2B)]
colnames(df2B.SMType) <- c("Nodo","SMType","Counts")

ggplot(df2B.SMType, aes(x = reorder(Nodo, -Counts), y = Counts, fill = SMType)) + 
  geom_bar(position="dodge",stat = "identity") +
  scale_fill_brewer(palette="Set3")+
  ggtitle("Marco de lectura 2 - Nodos con TGGCGCCA en el nodo parental - Tipo de transiciones por nucleotido") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Type")+
  geom_text(aes(label = Counts, group = SMType), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 11. En esta figura se muestra el número de tipos de sustitucion (transición o transversión) para cada nodo en el Marco de lectura 2

Marco de Lectura 3

## Cargo la tabla
FileRF3 = "codon_mutations_RF3.rooted.txt"
tableRF3 = read.table(FileRF3,header = TRUE,sep = "\t",row.names = NULL)


tableRF3<-tableRF3%>%
  filter(AncestorType=='SITE')
#tableRF3<-tableRF3%>%
#  filter(Spp!='PCC_7716')
tableRF3<-tableRF3%>%
  filter(SType!='Deletion')
#tableRF3<-tableRF3%>%
#  filter(SType!='NoMutation')

Spps <-unique(tableRF3[ , c(2)])

## HIP per site mutations
HIPmutC <- tableRF3[ , c(2,27,28,29,30,31,32,33,34)]
dfC = c()

for (i in 1:length(Spps)){
  spp <- Spps[i] # Especie en cuestión
  HIPmutC.spp <- HIPmutC%>%
    filter(Spp==spp) # Filtro la especie en cuestion de la nueva tabla
  for (j in 1:8){
    # Este nuevo vector "comprime" las mutaciones por sitio de cada nodo en un solo renglón haciendo una suma de todos los renglones
    Node.spp <- summarise(HIPmutC.spp,Nodo=unique(spp),NucPos=paste0('Nuc',j),Count=sum(HIPmutC.spp[j+1]))
    # Agrego este nuevo renglon al nuevo dataframe
    dfC <- rbind(dfC, Node.spp)
  }
}

## Gráfico
ggplot(dfC, aes(x = reorder(Nodo, -Count), y = Count, fill = NucPos)) + 
  geom_bar(position="dodge",stat = "identity") +
  scale_fill_brewer(palette="Set2")+
  ggtitle("Marco de lectura 3 - Nodos con TGGCGCCA en el nodo parental - Mutaciones por nucleotido") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Position")+
  geom_text(aes(label = Count, group = NucPos), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 12. En esta figura se muestra la cantidad de mutaciones por cada nucleótido de TGGCGCCA en el Marco de lectura 3

## HIP per mutation Type
TRC <- tableRF3[ , c(2,35,36,37,38,39,40,41,42)]
df2C = c()
for (i in 1:length(Spps)){
  spp <- Spps[i]
  TRC.spp <- TRC%>%
    filter(Spp==spp)
  for (j in 1:length(TRC.spp[,1])){
    for (k in 1:8){
      vector <- c(TRC.spp[j,1],TRC.spp[j,k+1])
      df2C <- rbind(df2C, vector)
    }
  }
}
df2C<- as.data.frame(df2C)
df2C.SMType <- data.table::setDT(df2C)[,list(Counts=as.numeric(.N)),names(df2C)]
colnames(df2C.SMType) <- c("Nodo","SMType","Counts")

ggplot(df2C.SMType, aes(x = reorder(Nodo, -Counts), y = Counts, fill = SMType)) + 
  geom_bar(position="dodge",stat = "identity") +
  scale_fill_brewer(palette="Set3")+
  ggtitle("Marco de lectura 3 - Nodos con TGGCGCCA en el nodo parental - Tipo de transiciones por nucleotido") +
  xlab("Nodo") + ylab("Conteo") + labs(fill = "Substitution Type")+
  geom_text(aes(label = Counts, group = SMType), position = position_dodge(0.9), vjust = -0.3, size = 3.5)

Figura 13. En esta figura se muestra el número de tipos de sustitucion (transición o transversión) para cada nodo enel Marco de lectura 3.