3.r

3.r

# Script para trabajar con Alineamientos
# [Goku Nivel 7]

Sys.setenv(http_proxy="http://danilo:[email protected]:3128/")
Sys.setenv(https_proxy="http://danilo:[email protected]:3128/")

# Paquetes a instalar!
source("https://bioconductor.org/biocLite.R")
biocLite("Biostrings")
biocLite("msa") #update all = a
#lazy loading NO KIDDING!

# Librerias
library(Biostrings)
library(seqinr)
library(phangorn)
library(msa)
library(ape)

# Una de las funciones más básicas en el análisis de secuencias es el alineamiento de secuencias.

# Una alineación de secuencias toma dos o más secuencias que se sospecha son similares y las "alinea", 
# una encima de la otra, haciendo coincidir residuos similares (bases) e insertando espacios (guiones) donde sea necesario.
# Las alineaciones de secuencia también pueden incluir "desajustes" - residuos o bases que no coinciden pero son ortólogas 
# lo que significa que se sabe que se encuentran en la misma posición en ambas secuencias a pesar de no ser idénticas.

# Las alineaciones de secuencia generalmente son el primer paso en cualquier análisis comparativo.

# Secuencia de "consenso": la secuencia de acuerdo, o la secuencia en común entre todas las especies representadas.
# Las alineaciones no tienen que ser entre secuencias de diferentes especies: algunos genes dentro de la misma 
# especie también compartirán similitudes.

# La alineación de secuencia más simple implica solo dos secuencias: esto se llama "pairwise alignment"

# Las alineaciones por pares se pueden usar para analizar dos secuencias directamente.
# Las alineaciones por pares también pueden ser un paso intermedio para crear "alineaciones de secuencias múltiples" más complejas o MSA.

# Vamos a comenzar alineando dos secuencias en R y luego realizando algunos análisis básicos sobre los resultados.

# Utilizaremos el paquete Biostrings de Bioconductor para este ejercicio,
# entonces, debe asegurarse de que su script comience con la cadena de carga estándar de Bioconductor,
# así como también con el comando para instalar el paquete Biostrings.

# Ahora, carguemos algunas secuencias para alinear.
#Clickto download the file prok.fasta and place it in your working directory for RStudio:
#https://prod-edxapp.edx-cdn.org/assets/courseware/v1/c309caedb1530152a48468b74c2f03b3/asset-v1:USMx+BIF003x+1T2018+type@asset+block/prok.fasta

# Leer el archivo prok.fasta en una nueva variable llamada prokaryotes
# Ahora agreguemos la siguiente línea en la parte inferior de su script:
prokaryotes <- read.fasta(file = "prok.fasta", seqtype = "DNA")

# Deberias ver que la variable "prokaryotes" aparece en el enviroment.  

# Leer las primeras dos secuencias en prokaryotes en seq1 y seq2 como "character strings"
# Necesitamos dividir el archivo en secuencias individuales primero para nuestro alineamiento de pares de secuencias.
# Para hacer eso, creemos dos variables nuevas, seq1 y seq2 agregando las siguientes líneas:
seq1 <- as.character(prokaryotes[[1]])
seq1 = paste(seq1, collapse = "")

seq2 <- as.character(prokaryotes[[2]])
seq2 = paste(seq2, collapse = "")

# El primer comando en cada par simplemente analiza la primera secuencia de procariotas en su propia variable.
# El segundo comando convierte eso en una cadena de texto (carácter) simple.
## Esto es necesario para que Biostrings alinee las secuencias.

# Alinear seq1 y seq2 utilizando la configuracion por defecto de Biostrings
pairwiseAlignment(pattern = seq2, subject = seq1)
# Tenga en cuenta que le está diciendo que hizo un alineamiento global, más sobre este último.
# No importa en qué orden pones las secuencias para este ejercicio, 
# pairwiseAlignment alineará de seq1 a seq2.
# Modifiquemos el código y guardemos esa alineación en una variable pairalign cambiando la última línea a:
pairalign <- pairwiseAlignment(pattern = seq2, subject = seq1)
# Vamos a ejecutar el código anterior y luego escribir en la consola
summary(pairalign)
#para ver el resumen de las estadisticas de este alneamiento

# Convertir el alineamiento a FASTA y guardar
# Ahora que tenemos una alineación, podemos exportarla de R como archivo FASTA,
# preservando los guiones insertados (indels)

# Podemos hacer eso con otro comando de Biostrings llamado writeXStringSet

# Primero necesitamos convertir la alineación a cadenas, usando un comando diferente de Biostrings 
# llamado BStringSet - agregue esta línea a su código:
pairalignString = BStringSet(c(toString(subject(pairalign)), toString(pattern(pairalign))))
# Observe cómo tuvimos que invocar tanto el "patrón" como el tema por separado
# Ahora agreguemos esta línea:
writeXStringSet(pairalignString, "aligned.txt", format = "FASTA")
# Ejecutamos estas líneas y abrimos el archivo de salida desde su pestaña ARCHIVOS para verlo
# Notarás un archivo FASTA bien formateado

# Crear un DotPlot utilizando seqinr
# Uno de los visuales de comparación de secuencias más antiguos es un gráfico de puntos.
# En su forma más simple, se genera un gráfico de puntos colocando un punto en la posición (i, j) si 
# el número de carácter i en la primera secuencia es igual que el número de carácter j en la segunda secuencia.
# Los diagramas se pueden usar para ver visualmente cómo son secuencias similares, y para ver si hay
# duplicación, inversiones, repeticiones, etc. entre los dos.

# seqinr incluye una rutina simple para crear gráficos de puntos a partir de dos secuencias llamadas dotPlot.
# Usemos secuencias de proteínas para esta comparación, ya que los gráficos de puntos con ADN se vuelven "desordenados" rápidamente.
# Haga clic para descargar el archivo cox1multi.fasta y guárdelo en su directorio de trabajo:
#https://prod-edxapp.edx-cdn.org/assets/courseware/v1/6ebb427f5f53343cbb7af97945366d02/asset-v1:USMx+BIF003x+1T2018+type@asset+block/cox1multi.fasta

# Primero cargar las nuevas secuencias
# Agreguemos estas lineas a tu script:
coxgenes <- read.fasta(file = "cox1multi.fasta", seqtype = "AA")
cox1 <- as.character(coxgenes[[1]])
cox2 <- as.character(coxgenes[[2]])

# Recuerde que esto carga los datos usando el comando read.fasta de seqinr y luego analiza las 
# dos primeras secuencias de este FASTA multisecuencia como datos de caracteres en cox1 y cox2 - 
# asegúrese de entender esto antes de continuar.

# Un simple DotPlot
# Ahora generemos el gráfico de puntos más simple posible: agregue esta línea y ejecútela:
dotPlot(cox1, cox2, main = "Human vs Mouse Cox1 Dotplot")
# Esto genera una gráfica de puntos simple y ser paciente ya que las gráficas de puntos tardan unos minutos en construirse.

# Haz clic para obtener un buen repaso sobre cómo interpretar las gráficas de puntos de tu clase de bioinformatica de ADN.
#https://ugene.net/wiki/pages/viewpage.action?pageId=4227426



# Un dotplot mejorcito con ventanas
# Hagamos una gráfica de puntos con los siguientes parámetros:
dotPlot(cox1, cox2, wsize = 3, wstep = 3, nmatch = 3, main = "Human vs Mouse Cox1 Dotplot\nwsize = 3, wstep = 3, nmatch = 3")
# Recuerde ir y venir entre las gráficas con las flechas 
# izquierda / derecha en la parte superior izquierda del panel PLOTS y ver cómo se diferencian
# El comando \ n en la línea anterior le dice a la rutina de trazado que mueva el texto que lo sigue a una nueva 
# línea, mientras que \ es un "carácter de escape" estándar utilizado para insertar controles en los comandos

# Un dotplot con ventana es "más fácil" de interpretar

# Probemos un último ejemplo: en este solo usaremos los primeros 100 residuos de las secuencias de aminoácidos

# Un dotplot de solo los primeros 100 residuos de las dos secuencias
# Recordemos que en R podemos especificar un subconjunto de un vector usando corchetes, 
# mientras que cox1 le dice a R que use todo el vector, cox1 [1: 100] le dice a R que use los residuos 1 a 100.
# Así que nuestra nueva línea de trazos debería verse así:
dotPlot(cox1[1:100], cox2[1:100], wsize = 3, wstep = 3, nmatch = 3, main = "Human vs Mouse Cox1 Dotplot\nwsize = 3, wstep = 3, nmatch = 3")


## Alineamiento Local vs Global 

# De sus clases de ADN y proteínas, debe recordar la diferencia entre las alineaciones locales y globales:
# En una alineación local, está haciendo coincidir su consulta con una subcadena (fragmento) de la secuencia de asunto.

# En una alineación global, realiza una alineación de extremo a extremo entre los dos.

# Puede terminar con una gran cantidad de gaps en la alineación global si los tamaños de la consulta y el tema son diferentes).
# Las alineaciones locales también pueden tener gaps.

# De hecho, hay una tercera opción: la alineación de superposición.
# Las alineaciones de superposición se usan cuando intenta ensamblar secuencias superpuestas, por ejemplo, 
# de múltiples ejecuciones de secuencia en un ensamblaje de genoma.

# Las alineaciones de superposición se enfocan en hacer las mejores alineaciones entre el final de una secuencia y el final de otra.

# Puede usar type = "<type>" para indicar pairwiseAlignment qué tipo de alineación está usando, p. type = "local".

# Hay otras opciones para la alineación por parejas que le permiten controlar la matriz de sustitución a
# usar (submat = "<tipo>"), así como las penalizaciones de apertura de espacio y extensión (gapOpening = X, gapExtension = Y).

# Por ahora puede experimentar con esto por su cuenta: usaremos varios ajustes más adelante en el curso.

## Ahora un poco de MAS (Multiple Sequence Alignments)

# Una Alineación de Secuencia Múltiple, o MSA, es exactamente lo que parece: una alineación
# de muchas (más de dos) secuencias.

# Si bien podríamos usar pairwiseAlignment de Biostrings para alinear múltiples secuencias de forma individual 
# con nuestro tema, esto no es eficiente ni tan preciso como el uso de un algoritmo "verdadero" de MSA como CLUSTALW.
# Recuerde instalar msa del paquete Bioconductor: esta es la implementación del Bioconductor de CLUSTALW.

# Leer los datos utilizando readXXStringset de Biostrings
# Usemos los Biostrings readAAStringSet y readDNAStringSet para leer nuestros archivos fasta en forma de juegos de caracteres
# Esto esencialmente hace lo que hicimos antes usando los comandos as.character y pegar, pero lo hace todo por nosotros

# Sin embargo, siempre hay muchas maneras de hacer lo mismo en R: vale la pena experimentar
# Añada las siguientes líneas a su secuencia de comandos y ejecútelas:
coxAA <- readAAStringSet("cox1multi.fasta")
prokDNA <- readDNAStringSet("prok.fasta")
# Estas dos líneas se leen exactamente en los mismos archivos fasta que utilizamos anteriormente, pero crean StringSets

#StringSets son colecciones de secuencias procesadas y pueden ser de varios tipos:
  #DNAStringSet
  #RNAStringSet
  #AAStringSet
  #BStringSet
# Los primeros tres se escriben para permitir solo un tipo específico de secuencia, p. ADN, ARN o AA
# A BStringSet puede contener cualquier tipo de secuencia.

# Utilizar un StringSet tipeado permite a otros análisis "saber" qué tipo de secuencia esperar, 
# lo que significa que no tenemos que declarar el tipo en muchas llamadas.

# Alinear las secuencias de AA en coxAA
# Ahora usemos los StringSets que creamos como entrada a la rutina msa:
coxAligned <- msa(coxAA)
# # Esto alinea las secuencias de proteína cox1 utilizando CLUSTALW y los parámetros de alineación predeterminados.

# Del mismo modo, podemos alinear las secuencias de ADN agregando y ejecutando esta línea:
prokAligned <- msa(prokDNA)

# Los conjuntos de datos coxAligned y prokAligned deberían estar ahora en su panel ENVIRONMENT, 
# véalos escribiendo su nombre en la consola (uno a la vez) y luego presione RETURN.
coxAligned
prokAligned

# La pantalla en la consola solo muestra parte de la alineación; podría ser 
# bueno poder ver las alineaciones en forma no truncada.

# Podemos aprovechar el propio comando de impresión del paquete msa para mostrar la alineación de una manera más amigable.

# Ingrese esta línea en la consola y presione RETURN:
print(prokAligned, show="complete")

# Puede hacer lo mismo para mostrar la alineación coxAligned si lo desea.
print(coxAligned, show="complete")

# Dependiendo de qué tan grande sea su pantalla, verá una gran parte de la alineación, 
# muy bien formateada y también puede desplazarse hacia arriba y hacia abajo.

# Observe que hay una nueva línea para el CONSENSO (con) de todas las secuencias.

# Naturalmente, hay muchas opciones disponibles para crear alineaciones de secuencias múltiples.
# Por ejemplo, podemos elegir qué algoritmo alineará las secuencias con:
msa(prokDNA, "ClustalW")
msa(prokDNA, "ClustalOmega")
msa(prokDNA, "Muscle")

# Cada uno de estos le dará resultados ligeramente diferentes, ya que cada algoritmo es ligeramente 
# diferente en la forma en que alinea las secuencias.

# Cada uno también comparte un grupo de parámetros de comando comunes, así como sus propios parámetros de comando únicos.

# Para ser breces, solo vamos a ver los parámetros comunes.

# Haga clic en el enlace de MSA para obtener más información; puede consultar la documentación de MSA:
#https://www.bioconductor.org/packages/devel/bioc/vignettes/msa/inst/doc/msa.pdf

# La documentación también analiza la rutina msaPrettyPrint, que le permite generar impresiones codificadas
# por colores y sombreadas; esto requiere que tenga instalado un binario LaTex funcional,que difiere del
# sistema operativo al sistema operativo, por lo que no vamos a verlo aquí. Siéntase libre de explorar por su cuenta.

# Los parámetros comunes que pueden usar los 3 algoritmos son:
 #cluster
 #gapOpening
 #gapExtension
 #Maxiters
 #SubstitutionMatrix
 #Order type
 #Verbose
# Cada uno se puede agregar de forma separada por comas al comando msa, p. msa (coxAA, cluster = "upgma").
# Estos son opcionales y pueden considerarse características "avanzadas".

#A menudo puede ser necesario exportar un MSA:
#P.ej. si desea verlo en una aplicación de escritorio como Seaview. Abrir Seqview en biolinux!
# O si necesita usarlo en un paquete de análisis que no sea R, p. un paquete de filogenia.

# Biostrings tiene dos formas de exportar alineaciones:
 #FASTA
#PHYLIP (PHYLIP es un paquete filogenético que exploraremos más adelante).

#Nota IMPORTANTE:
# Para exportar nuestras alineaciones como FASTA, primero debemos convertirlas en un StringSet con el que Biostrings pueda trabajar.
#Stringsets son colecciones de biostrings individuales.

# Exportar alineamientos como FASTA
# Comience escribiendo esta línea en su secuencia de comandos, esto guardará la alineación de 
# la secuencia procariota como conjunto de cadenas:
prokAlignStr = as(prokAligned, "DNAStringSet")
#Note que estamos especificando el tipo de stringset.
  # Esto tiene que coincidir con el tipo de alineación de entrada.
  # Nota que estamos usando una "mejor práctica" de usar nombres de variables que intentan indicar qué contiene la variable, por ejemplo prokAlignStr es nuestra secuencia procariota, alineada, guardada como conjunto de cuerdas ...

# Vamos a agregar esta línea a continuación:
writeXStringSet(prokAlignStr, file = "prokAligned.fasta")
# Esto simplemente escribe el stringset que creó anteriormente como un archivo FASTA en su directorio de trabajo
# Vamos a ejecutar estas dos líneas, luego vamos a tu directorio de trabajo y abrimos el archivo para verlo

# Puede abrir prokAligned.fasta de dos maneras:
  # Con un editor de texto como BBEdit o Sublime Text
  # Con un visor FASTA como Seaview

# Nota IMPORTANTE:
# Debes abrir el archivo original (prok.fasta) y compararlo con tu versión alineada (prokAligned.fasta) y
# ver las diferencias. Ahora revisemos las vistas de alineación usando BBEdit y Seaview.

# De forma similar, puede exportar la alineación de aminoácidos contenida en coxAligned, 
# asegurándose de especificar el tipo correcto de stringset:
coxAlignStr = as(coxAligned, "AAStringSet")
writeXStringSet(coxAlignStr, file = "coxAligned.fasta")

# Exportar a PHYLIP
# Abrir con Seaview
# Además del formato FASTA, puede exportar el formato PHYLIP utilizando Biostrings.
# Muchas aplicaciones pueden leer el formato # PHYLIP, difiere de FASTA en dos formas:
  # 1.Tiene un encabezado formal que define el número de secuencias, el número de posiciones en la alineación.
  # 2. Muestra las secuencias en un formato de bloque intercalado, donde cada bloque contiene los mismos residuos posicionados para CADA secuencia.

# Para exportar la alineación de cox en formato PHYLIP, agregue esta línea a su secuencia de comandos:
write.phylip(coxAligned, "coxAligned.phylip")
# Tenga en cuenta que NO tuvo que convertir la alineación a stringset.
# El comando write.phylip puede trabajar directamente con la alineación; el único otro parámetro 
# que debe enviar es el nombre del archivo de salida.

# Nota IMPORTANTE:
# Puedes abrir este archivo con BBEdit, Sublime Text o Seaview.

##¿Para qué sirven las MSA?

# Ahora que hemos creado Alineamientos de secuencia múltiple, ¿qué podemos hacer con ellos?

# Recordar que un MSA es una colección de secuencias similares donde los residuos ortólogos están en una columna.
  # Las secuencias pueden ser el mismo gen de diferentes especies
    #O
  # Las secuencias pueden ser genes diferentes (pero relacionados) de la misma especie o diferente

# En términos más simples, los MSA se usan para observar las similitudes y diferencias entre los grupos de secuencias.
# Ejemplo: Reconstrucción filogenética.


## Reconstruccion Filogenetica

# Recuerde de sus cursos anteriores que la filogenia de un grupo de secuencias es simplemente una
# representación de la relación de las secuencias, en otras palabras, qué tan similares son entre sí.

# Esto puede usarse para determinar la relación real en términos de historia evolutiva 
#(ascendencia de especie o ascendencia genética en el caso de genes duplicados).

# Los tres enfoques más simples para la reconstrucción filogenética en R son la distancia, 
# la parsimonia máxima y la máxima verosimilitud (maximun likehood).

# Comencemos con lo más simple: distancia
# En la parsimonia a distancia se calcula una matriz de distancia para cada posible combinación de pares de
# secuencias: esto se basa en la cantidad de diferencias que tienen, y a menudo se "corrige" para múltiples visitas.

#Nota IMPORTANTE:
#Para obtener más información sobre Reconstrucción filogenética, haga clic en el enlace Filogenia.
#https://en.wikipedia.org/wiki/Distance_matrices_in_phylogeny

# El paquete MSA nos permitirá convertir nuestras alineaciones a los diversos formatos que necesitaremos.

# Debe tener en cuenta que hay muchas formas de construir un árbol de distancia: el paquete PHYLIP, 
# para el que creamos una exportación de nuestra alineación anteriormente, es el más utilizado.

# Puedes usar secuencias de ADN o aminoácidos para todos los enfoques que vamos a ver aquí, pero debes tener 
# en cuenta que la mayoría de la reconstrucción filogenética usará genes elegidos específicamente porque 
# son buenos para la filogenia debido a su reloj tasa de mutación y restricciones evolutivas variables 
# (por ejemplo, 16S para procariotas o 18S para eucariotas).

# Para este ejemplo, usemos nuestro alineamiento de ADN procariota ya que las secuencias muestran más variación
# que las secuencias de aminoácidos que tenemos.

# Convertir alineamiento prokariotic a formato seqinr
# Comencemos por convertir nuestra alineación procariota al formato seqinr:
prokAligned2 <- msaConvert(prokAligned, type = "seqinr::alignment")

# Generar una matrix de distancia using seqinr
# Ahora podemos generar una matriz de distancia usando seqinr con el paquete ape
prokdist <- dist.alignment(prokAligned2, "identity")
# Esta línea usa dist.alignment para calcular la matriz de distancia de nuestro alineamiento prokAligned2 - 
# y las distancias se calculan en función de si las posiciones son idénticas (identidad).
# Si desea ver la matriz de distancia, escriba
prokdist

# Ahora podemos usar esta matriz para construir un árbol filogenético de distancia básica.

# Utilizaremos el método de unión de vecinos, que es un enfoque aglomerativo que utiliza las distancias por
# pares para agrupar las secuencias más cercanas primero, y luego agrega iterativamente secuencias a este par.

# Generar un Arbol de Distancia Neighbor-Joining
prokTree <- nj(prokdist)
#Esto crea el árbol y lo almacena en prokTree.

#Para ver el árbol puede escribir:
plot(prokTree)

# Nota IMPORTANTE:
# Haga clic para obtener más información sobre Neighbor Joining.
#https://en.wikipedia.org/wiki/Neighbor_joining

# Los árboles de distancia son "simplistas", no reflejan necesariamente una ascendencia evolutiva real, ya que
# se centran exclusivamente en la similitud absoluta, pero dan una representación muy rápida y rápida de la probable filogenia.

# Un enfoque más centrado en la evolución es la máxima parsimonia, que en realidad analiza los caracteres
# compartidos entre secuencias y maximiza la similitud,al tiempo que minimiza el número de 
# pasos necesarios para construir una historia filogenética entre los "cades" o especies.

# La máxima parsimonia es mucho más intensiva desde el punto de vista computacional y 
# solo puede realizarse por un número limitado (por ejemplo, 30 o menos) secuencias,
# mientras que la parsimonia a distancia es muy rápida y puede observar secuencias ilimitadas.

# Arbol de Máxima Parsimonia
# Afortunadamente, nuestro conjunto de datos es lo suficientemente pequeño como para que podamos hacer 
# un análisis de parsimonia de manera fácil y rápida.
# Recuerde instalar el paquete phangorn: install.packages ("phangorn")
# Naturalmente, también debería agregar la biblioteca (phangorn) a la sección de bibliotecas de su secuencia de comandos y ejecutar eso.

# Ahora podemos generar los datos filogenéticos, o phyDat, para nuestro árbol de parsimonia.
# Normalmente usarías la llamada PhyDat de phangorn para hacer eso desde una matriz de secuencias, pero
# también podemos usar msa para crear los datos correctos inmediatamente usando esta línea:
prokAligned3 <- msaConvert(prokAligned, type = "phangorn::phyDat")

# Ahora podemos alimentar esos datos en pratchet, que es la rutina en phangorn que genera árboles de parsimonia:
ParsTree <- pratchet(prokAligned3)
# pratchett realiza automáticamente reordenamientos de ramas para encontrar #trees "mejores" (más parsimoniosos) - 
# para este simple conjunto de datos de solo 6 secuencias, la respuesta vuelve bastante rápido. 

# Cuantas más secuencias haya en su conjunto de datos, más tardará.
# Ahora tiene una variable llamada ParsTree: puede trazar esto escribiendo
plot(ParsTree)
# Su árbol de parsimonia y árboles de distancia deberían ser idénticos en topología.
# Este no siempre será el caso; Cuantas más especies examines, más complicadas se vuelven las topologías.


## Arbol de Maxima Verosimilitud
# Crear un test estadistico inicial de la matrix de distancia con los datos
# Una variación en el enfoque de distancia es el método de máxima verosimilitud.

# Máxima verosimilitud le permite usar una matriz de distancia (como anteriormente) pero también incluye 
# la alineación de secuencia completa y calcula la probabilidad estadística de un árbol dados los datos completos.

# Como tal máxima verosimilitud es mucho más computacionalmente intensa, pero también puede proporcionar reconstrucciones más precisas que los enfoques de distancia "convencionales".

# El primer paso en la máxima verosimilitud es calcular la probabilidad utilizando la función pml.
fit <- pml(prokTree, prokAligned3)
# Recuerda que prokTree es un árbol de distancia que une vecinos, y prokAligned3 es la alineación completa.

# Ahora que hemos evaluado la probabilidad estadística, o ajuste, del árbol usando pml, podemos optimizar
#este árbol usando optim.pml.

# optim.pml requiere que proporcionemos el "modelo" de evolución que creemos que se ajustará mejor a nuestros datos.

# El modelo que utilizaremos para este ejercicio es el modelo de Jukes-Cantor, que asume que todas las bases
# se encuentran en números iguales, y que todas las mutaciones son igualmente probables (por ejemplo, A a T o C o G, etc.).

# Otro modelo "popular" es K80, o el modelo de dos parámetros de Kimura, que establece que las mutaciones 
# de transición (A <-> G) tendrán una velocidad diferente de las transversales (C <-> T).

# Nota IMPORTANTE:
# Haga clic para obtener más información sobre los modelos de sustitución evolutiva.
#https://en.wikipedia.org/wiki/Models_of_DNA_evolution


# Tratar de mejorar el arbol usando modelo Jukes-Cantor
# Para ejecutar optim.pml y optimizar el árbol que creó en forma, agregue esta línea a su secuencia de comandos y ejecútelo:
fitJC <- optim.pml(fit, model = "JC", rearrangement = "stochastic")
# En este comando, la reorganización indica qué tipo de cambio de rama intentará el algoritmo 
# cuando se intenta mejorar los tres: las opciones son estocásticas, NNI (vecino más cercano) y trinquete.

# Como con muchas cosas, debes experimentar: prueba diferentes modelos y tipos de reorganización y
# observa cómo difieren los árboles resultantes.

# Una buena regla general en filogenia es hacer muchos árboles diferentes a partir de los mismos datos, usando
# diferentes enfoques, y compararlos: cuando todos los árboles están de acuerdo, puede confiar en los resultados.

# Si miras la consola, verás que el ajuste en realidad no mejora.

# Una vez más, este es un conjunto de datos "simples" con solo 6 especies; no hay mucho margen de mejora.
##### Intenta escribir un comando que use K80 como modelo.

# Todo lo que queda ahora es ver el árbol de máxima verosimilitud como se muestra a continuación:
#Escribe
plot(fitJC)

# Si crea el árbol K80, trace eso y compare los resultados.

# Verá que los árboles se ven diferentes, pero solo porque la topología está invertida, 
# las relaciones reales son exactamente las mismas, solo las distancias entre las especies 
# representadas por la longitud de las ramas son diferentes.

# Lo último que debemos discutir con respecto a nuestras reconstrucciones filogenéticas es Bootstrapping.

# Bootstrapping es un método de submuestreo de datos que crea de manera efectiva nuevos conjuntos de datos (generalmente 100)
# del mismo tamaño que el original al muestrear aleatoriamente columnas de los datos existentes.

# Una vez que se construyen los nuevos conjuntos de datos, todos están sujetos al análisis y los resultados se fusionaron.

# Bootstrapping se puede utilizar para ver con qué fuerza el conjunto de datos en conjunto respalda las conclusiones, 
# esencialmente una medida estadística.

# Click para mas informacion en bootstrapping.
#https://en.wikipedia.org/wiki/Bootstrapping_(statistics)


#Agregar Bootstrap al arbol ML optimizado
# Bootstrapping se puede hacer con cualquiera de los algoritmos que hemos examinado,
# pero dado que esencialmente significa analizar 100 conjuntos de datos del mismo tamaño del original, 
# utilizaremos nuestro árbol existente de máxima verosimilitud Jukes-Cantor en lugar del árbol de parsimonia.

# Vamos a agregar esta línea a su script y ejecutarlo:
bootstrapped <- bootstrap.pml(fitJC, bs = 100, optNni = TRUE, multicore = TRUE, control = pml.control(trace=0))
# bootstrap.pml realiza un bootstrap paramétrico; en este caso, esto es más rápido que generar los 
# datasets por separado y analizarlos, dado que ya tenemos un árbol de ML ajustado, pero si desea generar
# los sets, puede usar bootstrap.phyDat y luego ejecutar el análisis original en esos conjuntos de datos.

# Tenga en cuenta que estamos utilizando nuestro árbol de ML ajustado, con 100 repeticiones de arranque 
# y el intercambio de vecino más cercano para el intercambio de ramas.


# Plot con bootstrap
# La variable bootstrapped ahora se encuentra en el panel ENVIRONMENT, todo lo que tenemos que hacer es visualizarlo.

# Podemos hacer eso usando este comando - agrégalo a tu script y presiona RUN:
plotBS(midpoint(fitJC$tree), bootstrapped, p = 50, type = "p")
# Esto usa la rutina de trazado de bootstrap phangorn y toma todos los árboles jC ajustados generados en el
# conjunto de datos bootstrapped y muestra resultados donde el soporte (número de árboles que acuerdan entre
# los 100) es mayor que 50; el tipo controla la visualización del árbol (aquí usamos p para el filograma).

# Hay muchas más opciones disponibles para todos los comandos que hemos visto 
# siempre puede usar el panel AYUDA en RStudio para ver las opciones y explorar más por su cuenta.