Cacher is the code snippet organizer for pro developers

library(tidyverse)
library(ggtree)
library(ape)
library(readxl)
library(viridis)

# Get tree data
tree_data <- read_xlsx(path = "clade_table_with_region.xlsx") %>%
  rename(label = ID) %>%
  mutate(group = case_when(age >= 50000 ~ "50+",
                           (age < 50000 & age >= 26500) ~ "Pre-LGM",
                           (age < 26500 & age >= 19000) ~ "LGM",
                           (age < 19000 & age >= 11600) ~ "Late-Glacial",
                           (age < 11600 & age >= 8200) ~ "Early-Holocene",
                           (age < 8200 & age >= 6000) ~ "Mid-Holocene",
                           (age < 6000 & age > 0) ~ "Late-Holocene",
                           TRUE ~ "Unknown")
  )

# Read tree and massage
tree <- read.tree(file = "final_tree_with_supports.nwk")
tree$tip.label <- gsub("\\{\\d+\\}", "", tree$tip.label)

# Massage tree data
data <- tree_data %>% 
  select(label, type, region, group) %>%
  mutate(present = 1) %>%
  complete(label, type, group, fill = list(present = 0)) %>%
  mutate(region = ifelse(region == "NA", NA, region))


# Create function to get data for heatmaps
get_data_gheatmap <- function(X){
  data %>%
    filter(group == X) %>% 
    mutate(region = ifelse(present == 0 , NA, region)) %>% 
    select(label, type, region) %>%
    distinct() %>% 
    pivot_wider(names_from = type, values_from = region) %>%
    as.data.frame() %>% 
    column_to_rownames("label")
} 

# Define groups
groups <- c("Unknown", "Late-Holocene", "Mid-Holocene", "Early-Holocene", "Late-Glacial", "LGM", "Pre-LGM", "50+")

# Get data for heatmaps
data_gheatmap <- map(groups, get_data_gheatmap)
names(data_gheatmap) <- groups

# Plot base tree
p <- ggtree(tree, layout = "fan", open.angle = 50) %<+% (data %>% filter(present == 1)) +
   geom_tiplab(size = 0,
                align = TRUE,
                linesize = 0.2,
                linetype = "dotted",
                color = "black") +
  theme(legend.position = "top",
        legend.key = element_blank()) # no keys

# Get some colors and names for color vector
colors <- c("#31688EFF", "#FDE725FF", "#35B779FF", "#440154FF")

names(colors) <- data %>% 
  select(region) %>%
  filter(!is.na(region)) %>% 
  distinct() %>% 
  deframe() 

# Plot heatmaps in the tree
offset <-  0.0012
p1 <- NULL

# Initial heatmap (Unknown)
p1 <- gheatmap(p, data_gheatmap$Unknown, width=.1, offset= 0, color = "black") + 
  theme(legend.position = "none") +
  scale_fill_manual(values = colors, na.value = grey.colors(8)[1], alpha = 0.6)

for (i in 2:length(data_gheatmap)) {
  p1 <- gheatmap(p1, data_gheatmap[[i]], width=.1, offset = offset, color = "black") + 
    theme(legend.position = "top") +
    
    scale_fill_manual(values = colors)
  offset <- offset + 0.001
}

# Rotate tree
p1 <- rotate_tree(p1, angle = 89)

# Save tree
ggsave(p1, filename = "tree.pdf", width = 8, height = 8)