library(tidyverse)   # para manipular tablas
library(phytools)    # para leer y graficar arboles
library(formattable) # para formatear tablas
# install.packages("formattable", dependencies = TRUE)

# Primero cargar datos
# session > set working directory > to source file location

# 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)
# 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 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

#################################################
# ¿cuántos genes ARO hay en los Ortólogos? #(#7)
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 == "N2") %>%   
  count(Species_Tree_Node)  #                    # (#8a)

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

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

###########################
# PLOT de GENE TREES # (#9)
  
# 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()