library(tidyverse)
library(phytools)
library(formattable)
# install.packages("formattable", dependencies = TRUE)

# Primero cargar datos
# session > set working directory > to source file location
library(rstudioapi)
current_path = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(current_path ))
getwd()

# CARGAR DATOS
Lstats <- read_tsv("Stat_Overall.tsv")          # estadísticas Globales
LstatsPerSpp <- read_tsv("StatsPerSpecies.tsv") # estadísticas x spp
LAnot <- readxl::read_xlsx("eggNOGanotacion.xlsx")  # Anotación todos los ortólogos  
Dupl <- read_tsv("Duplications_per_Orthogroup.tsv") #  Duplicación en cada Ortólogo
DetalleDupl <- read_tsv("DuplicationsDetalle.tsv") # 
DuplxSpp <- read_tsv("DuplwSpp.tsv")  # Duplicaciones x spp
Overlap <- read_tsv("Orthogroups_SpeciesOverlaps.tsv")
ltree  <- read.tree("SpeciesTree_rooted.txt") # árbol inferido con Ortólogos
labeltree  <- read.tree("SpeciesTree_rooted_node_labels.txt") # árbol inferido con Ortólogos

# ver el árbol (#1)
plotTree(ltree) # plot tree, función de phytools
plot(labeltree, show.node.label = TRUE) # plot + node labels

# PRIMERO EXPLORAMOS EL 'SUMMARY' #(#2)
Lstats  # nos da el summary de la búsqueda de ortólogos [OrtoFinder]

# ¿Cual es el número de genes en  Ortogrupos?
# ¿Número de ortogrupos de copia simple?

# PLOTS      
# "¿Que spp tiene mas duplicaciones?" #(#3)
DuplxSpp %>% arrange(desc(Duplications)) # es parte del clado lento o rápido?

# suma de todas las duplicaciones dentro de cada spp
DuplxSpp %>% summarise(sum(Duplications))
DuplTree   <- DuplxSpp %>% column_to_rownames("SpeciesNode") # genero objeto para plot 
plotTree.barplot(ltree,DuplTree)  # plot Nº duplicaciones, se ve un patrón?

# Overlap de ortologos..  #(#4)
# ¿que spp comparten más ortólogos? ¿coincide esto con el abrol de spp?
# misma tabla esta en Carpeta
Overlap <- read_tsv("Orthogroups_SpeciesOverlaps.tsv")
formattable(Overlap, list(area(col = C_jadinii:W_anomalus) ~ 
                            color_tile( "white", "red"))) # colores opcionales

# PLOTS Estadisticos x Spp  #(#5)
LstatsPerSpp # todas estadisticas x spp , ver x ej numero de genes

# genero vectores (otro objeto) para hacer plot despues
GenesTree  <- LstatsPerSpp %>% select(Species, Number_of_genes) %>% 
  column_to_rownames("Species") # GENES
OGTree <- LstatsPerSpp %>% select(Species, Number_of_genes_in_orthogroups) %>% 
  column_to_rownames("Species") # Ortogrupos
OGSppTree <- LstatsPerSpp %>% select(Species, Number_of_orthogroups_containing_species) %>% 
  column_to_rownames("Species") # Ortogrupos en estas spp

############################# (#5)
plotTree.barplot(ltree, GenesTree, args.barplot = list(main="N° de genes", 
                                                       col = "peachpuff1"))   # Nº de genes

plotTree.barplot(ltree, OGTree,  args.barplot = list(main="N° de genes en ortogrupos", 
                                                     col = "violet"))      # número genes en los ortogrupos

#Cargo anotaciones y lista de genes duplicados
LAnot # Anotaciones de todo (menos los single copy, estan en .xlsx también)
Dupl  # es solo una lista de duplicaciones de cada ortogrupo

# Anotacion inicial de TODOS los ortólogos, en todo el dataset...
# ejemplo con una Alcohol deshidrogenase mencionada en el artículo  #(#6)
LAnot %>% filter(grepl("SPE", Preferred_name)) # filtro, se pueden probar otros 

# junto Annotaciones +  lista de duplicaciones
JAnot <- LAnot %>% select(Orthologue, Preferred_name, Description, COG_category) %>% 
  left_join(Dupl, by = join_by(Orthologue == Orthogroup)) %>%  # juntamos dos tablas (o tibbles) 'Joint Anotado'
  rename(DA = `Duplications(all)`) %>% rename(name = Preferred_name) %>% rename(COG = COG_category) # renombramos para hacer mas leible todo
Dupl # número de ortologos 949

# ¿cuantos genes con duplicaciones tienen anotación?  #(#7)
JAnot %>% 
  arrange(desc(DA)) %>% # ordeno DA x descendiente 
  filter(0 != DA) # filtro valores NO zero

# ¿En que categorías COG estan presentes estos genes duplicados? #(#8)
#¿Que función cumplen?  #(#8) 
JAnot %>% arrange(desc(DA)) %>% filter(0 != DA) %>% # filtro de nuevo
  count(COG) %>% arrange(desc(n)) # cuento y ordeno x esa columna 

# ¿cuántos genes ARO hay en los Ortólogos? #(#9)
JAnot %>% filter(grepl("ARO", name)) # filtro

#¿Cuantas duplicaciones en todo el género Hanseniaspora?" ¿Que nodo corresponde a este ancestro?
DetalleDupl %>% filter(Species_Tree_Node == "N6") %>%   
  count(Species_Tree_Node)  #                    # (#11a)

#¿Cuantas duplicaciones en las otras levaduras?" # (#11b)
DetalleDupl %>% filter(Species_Tree_Node == "N1") %>%   
  count(Species_Tree_Node)  # es igual a tabla DuplTree 

# LAnot %>% filter(grepl("0000174", Orthologue))

DetalleAnot <- LAnot %>% select(Orthologue, Preferred_name, Description, COG_category, PFAMs) %>% 
  right_join(DetalleDupl, by = join_by(Orthologue == Orthogroup))

# PLOT de GENE TREES # (#16)
  
# que genes ARO se encuentran duplicados y en que genomas?
DetalleAnot %>% filter(grepl("ARO", Preferred_name))

read.tree("Gene_Trees/OG0000184_tree.txt") %>% plotTree()
read.tree("Gene_Trees/OG0000340_tree.txt") %>% plotTree()

# Cuantos ortologos con PDC1 se encuentran duplicados y en que genomas?
DetalleAnot %>% filter(grepl("PDC", Preferred_name))
read.tree("Gene_Trees/OG0000024_tree.txt") %>% plotTree()