
##################################################
##################################################
##                                              ##
##          Plot trees of fish species          ##
##                                              ##
##################################################
##################################################

# R version used:
# > R.version
# "R version 4.4.2 (2024-10-31)"

# Clear previous data
rm(list = ls())

# Load packages
require(ape)
require(dplyr)
require(ggnewscale)
require(ggplot2)
require(ggpmisc)
require(ggpubr)
require(ggpp)
require(ggrepel)
require(ggtree)
require(grDevices)
require(grid)
require(gridSVG)
require(grImport2)
require(openxlsx)
require(phylobase)
require(phylosignal)
require(phytools)
require(rphylopic)
require(rsvg)
require(tibble)
require(tidytree)

# Loading data file
my_fish <- openxlsx::read.xlsx("Supplementary_file_S4.xlsx")

# Tree obtained from:
# Additional file 2 (complete tree in newick format), in
# Betancur-R et al. (2017)
# https://doi.org/10.1186/s12862-017-0958-3

# Loading tree
fishtree_orig <- ape::read.tree("Fish_tree.tre")

# Ladderize order of tips in tree
fishtree <- ape::ladderize(fishtree_orig, right = FALSE)

# Making sure the tree is ultrametric
ape::is.ultrametric(fishtree)


##################################################
# Modifying data file
##################################################

# Selecting only unique species names, no duplicates
my_fish_unique <- my_fish %>%
  dplyr::distinct(.keep_all = TRUE)

# Remove species not in tree from data file
store_values_dump1 <- NULL
for (i in 1:nrow(my_fish_unique)) {
  temp <- sum(grepl(my_fish_unique[i,3],
                    fishtree$tip.label))
  store_values_dump1 <- c(store_values_dump1, temp)
}

my_fish_unique_species_in_tree <- my_fish_unique %>%
  dplyr::filter(store_values_dump1 > 0)

# Store branch tip number of species in both tree and data file
tree_tip_numbers <- NULL
for (i in 1:nrow(my_fish_unique_species_in_tree)) {
  temp <- grep(my_fish_unique_species_in_tree[i,3],
               fishtree$tip.label)
  tree_tip_numbers <- c(tree_tip_numbers, temp)
}

my_fish_unique_species_in_tree_2 <- cbind(my_fish_unique_species_in_tree,
                                          tree_tip_numbers)

# Store exact character strings of species names from tree
my_fish_unique_species_in_tree_2$tree_tips_to_include <- fishtree$tip.label[my_fish_unique_species_in_tree_2$tree_tip_numbers]

# Store names of species not in tree
store_values_dump2 <- NULL
for (i in 1:nrow(my_fish_unique)) {
  temp <- sum(grepl(my_fish_unique[i,3],
                    fishtree$tip.label))
  store_values_dump2 <- c(store_values_dump2, temp)
}

my_fish_unique_species_NOT_in_tree <- my_fish_unique %>%
  dplyr::filter(store_values_dump2 < 1)


##################################################
# Adding missing species to tree
##################################################

# Species missing from tree:
# Catostomus_commersonii
# Catostomus_discobolus
# Catostomus_latipinnis
# Xyrauchen_texanus
# Clarias_gariepinus
# Clupea_pallasii
# Cottus_bairdii
# Liza_klunzingeri
# Oncorhynchus_clarkii
# Oncorhynchus_mykiss
# Salvelinus_namaycush
# Epinephelus_marginatus
# Acanthopagrus_australis

# Species missing from tree without representation from the same genus in tree:
# Clarias_gariepinus (Family: Clariidae)
# Clupea_pallasii (Family: Clupeidae)
# Xyrauchen_texanus (Family: Catostomidae)

node <- which(fishtree$tip.label == "Catostomidae_Catostomus_plebius_I2585-18")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Catostomus_commersonii",
                               where = node,
                               position = (8.3052488/2)) # 8.3052488 is the edge length from the
# tip to the internal node, where Catostomidae branches from other cypriniform fish in the original tree

node <- which(fishtree$tip.label == "Catostomidae_Catostomus_plebius_I2585-18")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Catostomus_discobolus",
                               where = node,
                               position = (8.3052488/2)) # 8.3052488 is the edge length from the
# tip to the internal node, where Catostomidae branches from other cypriniform fish in the original tree

node <- which(fishtree$tip.label == "Catostomidae_Catostomus_plebius_I2585-18")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Catostomus_latipinnis",
                               where = node,
                               position = (8.3052488/2)) # 8.3052488 is the edge length from the
# tip to the internal node, where Catostomidae branches from other cypriniform fish in the original tree

node <- which(fishtree$tip.label == "Cottidae_Cottus_carolinae_E00281")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Cottus_bairdii",
                               where = node,
                               position = 0.001)

node <- which(fishtree$tip.label == "Mugilidae_Liza_richardsonii_E00808")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Liza_klunzingeri",
                               where = node,
                               position = 0.001)

node <- which(fishtree$tip.label == "Salmonidae_Oncorhynchus_nerka_mykiss_E00437")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Oncorhynchus_clarkii",
                               where = node,
                               position = (19.947525/2)) # 19.947525 is the edge length from the
# tip to the internal node, where Oncorhynchus branches from other salmonid fish in the original tree

node <- which(fishtree$tip.label == "Salmonidae_Oncorhynchus_nerka_mykiss_E00437")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Oncorhynchus_mykiss",
                               where = node,
                               position = (19.947525/2)) # 19.947525 is the edge length from the
# tip to the internal node, where Oncorhynchus branches from other salmonid fish in the original tree

node <- which(fishtree$tip.label == "Salmonidae_Salvelinus_alpinus_G01473")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Salvelinus_namaycush",
                               where = node,
                               position = 0.001)

node <- which(fishtree$tip.label == "Serranidae_Epinephelus_maculatus_E00549")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Epinephelus_marginatus",
                               where = node,
                               position = 0.001)

node <- which(fishtree$tip.label == "Sparidae_Acanthopagrus_catenula_E00953")
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Acanthopagrus_australis",
                               where = node,
                               position = 0.001)

node <- which(fishtree$tip.label == "Clariidae_Channallabes_apus_P02_D10_RG08") # Same taxonomic Family rather than Genus
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Clarias_gariepinus",
                               where = node,
                               position = 0.001)

node <- which(fishtree$tip.label == "Clupeidae_Jenkinsia_lamprotaenia_G01339") # Same taxonomic Family rather than Genus
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Clupea_pallasii",
                               where = node,
                               position = 0.001)

node <- which(fishtree$tip.label == "Catostomidae_Ictiobus_niger_I2706-187") # Choice of location becomes somewhat
# arbitrary for this species, given that
# multiple equally unrelated species of Catostomidae are represented in the tree
fishtree <- phytools::bind.tip(fishtree,
                               tip.label = "Xyrauchen_texanus",
                               where = node,
                               position = 0.001)


##################################################
##################################################
##                                              ##
##              Plot of small tree              ##
##                                              ##
##################################################
##################################################

##################################################
# Modifying tree
##################################################

fishtree_small <- fishtree

# Renaming tip labels
old_names <- c(my_fish_unique_species_in_tree_2$tree_tips_to_include,
               my_fish_unique_species_NOT_in_tree$species_name)
new_names <- c(my_fish_unique_species_in_tree_2$species_name,
               my_fish_unique_species_NOT_in_tree$species_name)

fishtree_small$tip.label <- new_names[match(fishtree_small$tip.label, old_names)]

# Pruning tree, dropping the species that are not in the data file
species_to_keep_1 <- new_names
to.drop_1 <- setdiff(fishtree_small$tip.label, species_to_keep_1)
fishtree_small_cut <- ape::drop.tip(fishtree_small, to.drop_1)

# Renaming internal nodes in tree to have them be different from the numbering of the tips
fishtree_small_cut$node.label <- c(101:(length(fishtree_small_cut$node.label) + 100))

# Ladderizing tree for visualization in plot
fishtree_small_cut_ladderized <- ape::ladderize(fishtree_small_cut, right = FALSE)


##################################################
# Obtaining species icons for plot
##################################################

# Function to override errors due to missing data when retrieving icons
image_function <- function (X) {
  return(tryCatch(rphylopic::get_uuid(X, n = 1), error = function(e) "NA"))
}

# Obtain image IDs (uuid's) and reshape data for plot
image_IDs <- lapply(new_names, image_function) # Executing function to collect image IDs
image_IDs_2 <- unlist(image_IDs)
image_IDs_final <- as.data.frame(image_IDs_2)

image_IDs_final$species_number <- c(1:nrow(image_IDs_final)) # Adding unique identifier for each species
phylopic_all <- data.frame(name = new_names,
                           uid = image_IDs_final$image_IDs_2) # Using data frame structure needed for 'geom_tiplab()' function

phylopic_info <- phylopic_all %>% dplyr::filter(uid != "NA") # Filter out NAs, i.e. species without icons


##################################################
# Final processing of names and icons
##################################################

# Removing underscores from species names in data file
phylopic_info2 <- phylopic_info %>%
  dplyr::mutate(name = gsub("_", " ", name))

# Removing underscores from species names in tree
fishtree_small_cut_ladderized2 <- fishtree_small_cut_ladderized
placeholder_for_names <- as.data.frame(fishtree_small_cut_ladderized$tip.label)
colnames(placeholder_for_names) <- c("species_names")
removed_underscores <- placeholder_for_names %>%
  dplyr::mutate(species_names = gsub("_", " ", species_names))
fishtree_small_cut_ladderized2$tip.label <- removed_underscores$species_names

# To classify polytomies as true polytomies
fishtree_small_cut_ladderized3 <- ape::di2multi(fishtree_small_cut_ladderized2,
                                                tol = 0.002)

# Splitting up images to be flipped or rotated
phylopic_info3 <- phylopic_info2 %>%
  dplyr::mutate(
    unedited = case_when(
      name != "Cyclopterus lumpus" &
        name != "Sander vitreus" &
        name != "Pseudopleuronectes americanus" &
        name != "Clarias gariepinus" &
        name != "Acipenser fulvescens" &
        name != "Oncorhynchus mykiss" ~ uid),
    flipped = case_when(
      name == "Cyclopterus lumpus" |
        name == "Sander vitreus" |
        name == "Pseudopleuronectes americanus" |
        name == "Clarias gariepinus" |
        name == "Acipenser fulvescens" ~ uid),
    rotated = case_when(
      name == "Oncorhynchus mykiss" ~ uid))

# Loading local .SVG images of fish
file1 <- "Salmo_salar_Blue.svg"
svg1 <- grImport2::readPicture(rawToChar(rsvg::rsvg_svg(file1)))
file2 <- "Leuciscus_leuciscus_Blue.svg"
svg2 <- grImport2::readPicture(rawToChar(rsvg::rsvg_svg(file2)))
file3 <- "Catostomus_latipinnis_Blue.svg"
svg3 <- grImport2::readPicture(rawToChar(rsvg::rsvg_svg(file3)))
file4 <- "Epinephelus_marginatus_Blue.svg"
svg4 <- grImport2::readPicture(rawToChar(rsvg::rsvg_svg(file4)))



##################################################
# Main plot
##################################################

small_tree_plot <- ggtree::ggtree(fishtree_small_cut_ladderized3) %<+%
  phylopic_info3 +
  ggtree::geom_tiplab(offset = 60, size = 2) +
  ggtree::geom_rootedge(rootedge = 10) +
  xlim(-10, 535) +
  rphylopic::geom_phylopic(aes(uuid = unedited, x = 380),
                           alpha = 0.7,
                           color = 'steelblue',
                           width = 50) +
  rphylopic::geom_phylopic(aes(uuid = flipped, x = 380),
                           alpha = 0.7,
                           color = 'steelblue',
                           horizontal = TRUE, # to horizontally flip / mirror the image
                           width = 50) +
  rphylopic::geom_phylopic(aes(uuid = rotated, x = 380),
                           alpha = 0.7,
                           color = 'steelblue',
                           angle = 7, # to rotate image
                           width = 50) +
  ggpp::annotate(geom = "grob_npc", npcx = 0.744, npcy = 0.403,
                 vp.width = 0.095, vp.height = 0.095,
                 label = grImport2::pictureGrob(svg1,
                                                ext = "gridSVG",
                                                gpFUN = function(gp) gpar(col = NA,
                                                                          fill = "steelblue",
                                                                          alpha = 0.84))) +
  ggpp::annotate(geom = "grob_npc", npcx = 0.744, npcy = 0.173,
                 vp.width = 0.095, vp.height = 0.095,
                 label = grImport2::pictureGrob(svg2,
                                                ext = "gridSVG",
                                                gpFUN = function(gp) gpar(col = NA,
                                                                          fill = "steelblue",
                                                                          alpha = 0.84))) +
  ggpp::annotate(geom = "grob_npc", npcx = 0.744, npcy = 0.271,
                 vp.width = 0.095, vp.height = 0.095,
                 label = grImport2::pictureGrob(svg3,
                                                ext = "gridSVG",
                                                gpFUN = function(gp) gpar(col = NA,
                                                                          fill = "steelblue",
                                                                          alpha = 0.84))) +
  ggpp::annotate(geom = "grob_npc", npcx = 0.744, npcy = 0.896,
                 vp.width = 0.095, vp.height = 0.095,
                 label = grImport2::pictureGrob(svg4,
                                                ext = "gridSVG",
                                                gpFUN = function(gp) gpar(col = NA,
                                                                          fill = "steelblue",
                                                                          alpha = 0.84)))


##################################################
##################################################
##                                              ##
##          Plot of full (massive) tree         ##
##                                              ##
##################################################
##################################################

# Species that are not fish but are in the tree:
# Western clawed frog: "Pipidae_Xenopus_tropicalis_G01511"
# Gray short-tailed opossum: "Didelphidae_Monodelphis_domestica_G01378"
# House mouse: "Muridae_Mus_musculus_G01382"
# Humans: "Hominidae_Homo_sapiens_G01328"

# Pruning tree, dropping four tetrapod (non-fish) species
species_to_keep <- fishtree$tip.label[!(fishtree$tip.label == "Pipidae_Xenopus_tropicalis_G01511" |
                                          fishtree$tip.label == "Didelphidae_Monodelphis_domestica_G01378" |
                                          fishtree$tip.label == "Muridae_Mus_musculus_G01382" |
                                          fishtree$tip.label == "Hominidae_Homo_sapiens_G01328")]
to.drop <- setdiff(fishtree$tip.label, species_to_keep)
fishtree_cut <- ape::drop.tip(fishtree, to.drop)

# To classify polytomies as true polytomies
fishtree_cut2 <- ape::di2multi(fishtree_cut, tol = 0.002)

# Renaming tip labels
fishtree2 <- fishtree_cut2
fishtree2$node.label <- rep("", length(fishtree2$node.label))

# Creating data frame where information on presence/absence is included
species_in_tree <- as.data.frame(cbind(Species = fishtree2$tip.label,
                                       In_study = "No", In_study_as_number = 0))
species_in_tree$In_study_as_number <- as.numeric(species_in_tree$In_study_as_number)
species_in_data_set <- c(my_fish_unique_species_in_tree_2$tree_tips_to_include,
                         my_fish_unique_species_NOT_in_tree$species_name)
species_in_data_set2 <- as.data.frame(cbind(Species = species_in_data_set,
                                            In_study = "Yes", In_study_as_number = 1))
species_in_data_set2$In_study_as_number <- as.numeric(species_in_data_set2$In_study_as_number)
excluding_data_set_species <- dplyr::anti_join(species_in_tree,
                                               species_in_data_set2,
                                               by = "Species")
dat_with_scores <- rbind(excluding_data_set_species,
                         species_in_data_set2)

dat_with_scores2 <- dat_with_scores %>%
  dplyr::mutate(salmonids = case_when(
    Species == "Salmonidae_Coregonus_clupeaformis_G01263" |
      Species == "Salmonidae_Thymallus_brevirostris_G01502" |
      Species == "Salmonidae_Oncorhynchus_nerka_mykiss_E00437" |
      Species == "Salmonidae_Salvelinus_alpinus_G01473" |
      Species == "Salmonidae_Salmo_salar_G01522" |
      Species == "Salmonidae_Parahucho_perryi_G01418" |
      Species == "Salvelinus_namaycush" |
      Species == "Oncorhynchus_mykiss" |
      Species == "Oncorhynchus_clarkii" ~ 1,
    TRUE ~ 0))
dat_with_scores2$salmonids <- as.factor(dat_with_scores2$salmonids)

x <- tidytree::full_join(tidytree::as_tibble(fishtree2),
                         dat_with_scores2,
                         by = c("label" = "Species"))
tree2 <- tidytree::as.treedata(x)

# Obtain node index of MRCA node for salmonids and suckers.
# Selecting two ancestor species of the most basal node of e.g. salmonids in the tree,
# where that node is their most recent common ancestor.
back_to_phylo <- ape::as.phylo(tree2)
index_for_MRCA_node <- ape::getMRCA(back_to_phylo, tip = c("Salmonidae_Coregonus_clupeaformis_G01263",
                                                           "Salmonidae_Salmo_salar_G01522"))
index_for_MRCA_node2 <- ape::getMRCA(back_to_phylo, tip = c("Catostomidae_Ictiobus_niger_I2706-187",
                                                            "Catostomidae_Thoburnia_rhothoeca_I2703-183"))

# Color branches based on whether they are downstream of the node
tree2_grouped <- tidytree::groupClade(tree2, c(index_for_MRCA_node, index_for_MRCA_node2))

# Renaming species labels for plotting
tree2_renamed <- tree2_grouped
tree2_renamed@phylo$tip.label <- tree2_renamed@phylo$tip.label[match(tree2_renamed@phylo$tip.label,
                                                                     species_in_data_set2$Species)]
tree2_renamed@phylo$tip.label[!is.na(tidytree::node.label(tree2_renamed@phylo, node = "external"))] <- c("   ")


##################################################
# Plot
##################################################

large_tree_plot1 <- ggtree::ggtree(tree2_renamed,
                                   size = 0.5,
                                   layout = "circular",
                                   aes(color = group)) +
  scale_color_manual(values = c("0" = "black", "1" = "orange", "2" = "red"),
                     na.value = "black", guide = "none") +
  ggtree::geom_hilight(node = index_for_MRCA_node, fill = NA, alpha = 0.9)

# Plot tree
large_tree_plot2 <- large_tree_plot1 +
  ggrepel::geom_label_repel(aes(label = label),
                            max.overlaps = Inf,
                            segment.color = "steelblue",
                            label.r = 0.5,
                            label.size = 1,
                            min.segment.length = 0.0000001,
                            force = 10,
                            col = alpha(c("steelblue"), 1),
                            fill = alpha(c("steelblue"), 1),
                            xlim = c(-Inf, Inf),
                            ylim = c(-Inf, Inf)) +
  ggnewscale::new_scale_color() +
  ggnewscale::new_scale_fill() +
  scale_color_manual(values = c("Yes" = "black"),
                     na.value = NA, guide = "none") +
  scale_fill_manual(values = c("Yes" = "steelblue"),
                    na.value = NA) +
  theme(legend.position = "none",
        plot.margin = unit(c(10, 10, 10, 10), "mm")) +
  ggpubr::theme_transparent() +
  xlim(0, 463)


##################################################
# Combine and save plots
##################################################

# Open separate XQuartz plotting window
grDevices::x11(width = 7, height = 7) # You may need to install XQuartz for this
large_tree_plot2

# Save large tree - using set.seed() for replication of random placement of 'ggrepel' labels
set.seed(1);ggplot2::ggsave("large_tree.svg",
                            bg = "transparent")

dev.off()

# Import the saved large tree
file <- "large_tree.svg"
svg <- grImport2::readPicture(rawToChar(rsvg::rsvg_svg(file)))
vector_image_to_add <- grImport2::pictureGrob(svg, ext = "gridSVG")

# Adjust relative location of large tree in plotting window
image_parameters <- grid::viewport(x = unit(0.20, "npc"),
                                   y = unit(0.80, "npc"),
                                   width = unit(0.42, "npc"),
                                   height = unit(0.42, "npc"))
vector_image_to_add_2 <- grid::grobTree(vector_image_to_add, vp = image_parameters)

# Open separate XQuartz plotting window
grDevices::x11(width = 7, height = 7)
small_tree_plot
grid::grid.draw(vector_image_to_add_2)

# Save image of large and small tree combined
grDevices::dev.print(pdf,
                     file = "Combined_plot.pdf")
dev.off()


##################################################
##################################################
##                                              ##
##        Test of phylogenetic signal           ##
##                                              ##
##################################################
##################################################

# Use species names as row names
dat_with_scores_species_row_names <- dat_with_scores %>%
  tibble::remove_rownames() %>%
  tibble::column_to_rownames(var = "Species")

dat_for_analysis <- dat_with_scores_species_row_names %>%
  dplyr::select(!In_study)

# Creating phylosig object
phylosig_object <- phylobase::phylo4d(fishtree2, dat_for_analysis)

# Running model on phylogenetic signal
phylosig_mod <- phylosignal::phyloSignal(p4d = phylosig_object,
                                         method = "all")
phylosig_mod

