make_primate_mating_system_tree.R

#!/usr/bin/env Rscript

library(ggtree)
library(ape)
library(tidytree)
library(scales)
library(colorspace)
library(ggplot2)

prim = read.tree("morgan.tre")

# Get outgroup node index
out.ind = as_tibble(prim)[as_tibble(prim)$label == "HLnycHum2",]$parent

# Root the tree
prim.rt = root(phy = prim, node = out.ind)

# Hack to fix one node with branch length input wrong?
prim.rt$edge.length[is.nan(prim.rt$edge.length)] = 0.0022076425

# Bring in species info originally from Google Sheet
sp.info = read.csv("mating_system_tree_info.csv")

# Extract the tip labels from the tree
tip.lbls = prim.rt$tip.label

# Create a named vector for easy lookup
name.lookup = setNames(sp.info$species.name, sp.info$id)

# Update the tip labels if the species name is not an empty string
new.lbls = sapply(tip.lbls, function(x) ifelse(name.lookup[x] != "", name.lookup[x], x))

# Assign the new labels back to the tree
prim.rt$tip.label = new.lbls

# Create a data frame that matches the tree tips to the variable
tip.data = data.frame(label = new.lbls, variable = sp.info$mating.sys[match(tip.lbls, sp.info$id)])

# Ensure the order of tip.data matches the tree tip labels order
tip.data = tip.data[match(prim.rt$tip.label, tip.data$label), ]

# Rotate nodes to ensure pairs are correctly ordered
rotate_pairs = function(tree, labels) {
  pairs = grep(" 1$", labels, value = TRUE)
  for (pair in pairs) {
    base_name = gsub(" 1$", "", pair)
    pair_2 = paste0(base_name, " 2")
    if (pair_2 %in% labels) {
      node = MRCA(tree, which(tree$tip.label == pair), which(tree$tip.label == pair_2))
      y1 = ycoord(tree, pair)
      y2 = ycoord(tree, pair_2)
      if (y1 < y2) {  # Rotate only if "2" is above "1"
        tree = rotate(tree, node)
      }
    }
  }
  return(tree)
}

# Get y-coordinates of tips
ycoord = function(tree, label) {
  tree_view = ggtree(tree)
  tip_info = tree_view$data[tree_view$data$isTip, ]
  return(tip_info$y[tip_info$label == label])
}

# Rotate the nodes in the tree
prim.rt = rotate_pairs(prim.rt, prim.rt$tip.label)

# Define custom colors
point.cols = c("Multimale" = "blue", "Unimale" = "red", "Unset" = "grey")
label.cols = c("Multimale" = darken("blue", 0.3), "Unimale" = darken("red", 0.3), "Unset" = darken("grey", 0.3))

# Hide outgroup
prim.rt = drop.tip(prim.rt, "Nycticeius humeralis")

p = ggtree(prim.rt) %<+% tip.data +
  geom_tiplab(aes(label = paste0("italic('", label, "')"),
                  color = variable),
              parse = TRUE, offset = 0.001, show.legend = FALSE) +
  geom_tippoint(aes(fill = variable, color = variable),
                size = 3, stroke = 1.0,
                shape = 21, show.legend = TRUE) +
  guides(fill = guide_legend(override.aes = list(stroke = 0.7))) +
  scale_fill_manual( name = "Mating System",
                     values = point.cols) +
  scale_color_manual(name = "Mating System",
                     values = label.cols,
                     guide = "none") +
  theme(legend.position = "bottom", 
        legend.box = "horizontal",
        legend.title = element_text(size = 10, face = "bold"),
        legend.background = element_rect(color = "black",
                                         fill = NA),
        plot.margin = margin(5, 20, 5, 5)) +
  ggplot2::xlim(0, 0.09)

# Save the plot
ggsave(p, filename="mating_sys_tree.pdf", width = 8, height = 10)