#requirements: tree files (Newick format) for five major orders ('COITree_5major.tre') and all 'other' orders ('COITree_Others.tre')
#requirements: file TreeData.txt, which contains all associated metadata in tab-delimited format

#set working directory and install required packages
setwd("/PATH/TO/WORKING/DIRECTORY")
require(ggplot2)
require(scales)
require(grid)
require(gridExtra)
require(phytools)
require(ggtreeExtra)
require(ggtree)
require(RColorBrewer)

####################################
####### FIGURE 8 - MAIN ORDERS #####
####################################

#load tree file and extra data file
tree <- read.tree("COITree_5major.tre")
data.in <- read.csv("TreeData.txt", sep = "\t", header = TRUE)

#reduce data to only species contained within tree (if necessary)
data <- data.in[data.in$Tips %in% tree[4][[1]],]

#set order of Orders for plot legend
data$Order <- factor(data$Order, levels = c("Lepidoptera","Hymenoptera","Coleoptera","Diptera","Hemiptera","Blattodea","Ephemeroptera","Megaloptera","Neuroptera","Odonata","Orthoptera","Phasmatodea","Plecoptera","Siphonaptera","Strepsiptera","Thysanoptera","Trichoptera"))

#add custom colours
colours <- c("#cc0011","#ffaa00","#66cc33","#443333","#0099cc")

#Create tree object and attach associated data
###NUMT Count###
g.tree <- ggtree(tree, aes(color = Order), branch.length="none", layout = "fan", open.angle = 7.5, size = 0.1) %<+% data
rotate_tree(g.tree +
  scale_color_manual(breaks = levels(data$Order), values = colours) +
  geom_fruit(geom = geom_bar,
             mapping = aes(x = NUMT_Count_Total, fill = Order),
             orientation = "y",
             stat = "identity",
             pwidth = 1,
             grid.params = list(vline = F, alpha = 0.5),
             axis.params = list(axis = "x", 
                                nbreak = 6,
                                text.size = 7.5, text.angle = 0,
                                hjust = -0.1)) +
  #geom_text(aes(label=node), hjust=-.3, size = 1) +
  geom_cladelabel(node = 1108, label = "5.7 ± 9.6", color = "#cc0011", fill = "white", fontsize = 6, offset = 7, barsize = NA, geom = 'label', align = T) + #lep
  geom_cladelabel(node = 906, label = "17.7 ± 32.1", color = "#ffaa00", fill = "white", fontsize = 6, offset = 10, barsize = NA, geom = 'label', align = T) + #hym
  geom_cladelabel(node = 1220, label = "9.7 ± 12.7", color = "#66cc33", fill = "white", fontsize = 6, offset = 10, barsize = NA, geom = 'label', align = T) + #col
  geom_cladelabel(node = 823, label = "8 ± 15.4", color = "#443333", fill = "white", fontsize = 6, offset = 13, barsize = NA, geom = 'label', align = T) + #dip
  geom_cladelabel(node = 674, label = "23 ± 45.1", color = "#0099cc", fill = "white", fontsize = 6, offset = 10, barsize = NA, geom = 'label', align = T) + #hem
  scale_fill_manual(values = colours) +
  guides(color = "none", fill = guide_legend(override.aes = list(size = 3))) +
  labs(title = "NUMT Count") +
  theme(legend.key=element_rect(fill=NA),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 15, face = "bold"),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = -10),
        plot.margin = unit(c(0,0,0,0),"cm")),
  90)

ggsave("Figure 8_NUMTCountMain.pdf", width = 18, height = 16, units = "in", dpi = 300)

###Genome Size###
g.tree <- ggtree(tree, aes(color = Order), branch.length="none", layout = "fan", open.angle = 7.5, size = 0.1) %<+% data
rotate_tree(g.tree +
              scale_color_manual(breaks = levels(data$Order), values = colours) +
              geom_fruit(geom = geom_bar,
                         mapping = aes(x = Genome_Asssembly_Size_Mb, fill = Order),
                         orientation = "y",
                         stat = "identity",
                         pwidth = 1,
                         grid.params = list(vline = F, alpha = 0.5),
                         axis.params = list(axis = "x", 
                                            nbreak = 6,
                                            text.size = 7.5, text.angle = 0,
                                            hjust = -0.1)) +
              #geom_text(aes(label=node), hjust=-.3, size = 1) +
              geom_cladelabel(node = 1108, label = "529 ± 197.6", color = "#cc0011", fill = "white", fontsize = 6, offset = 7, barsize = NA, geom = 'label', align = T) + #lep
              geom_cladelabel(node = 906, label = "330.1 ± 323", color = "#ffaa00", fill = "white", fontsize = 6, offset = 7, barsize = NA, geom = 'label', align = T) + #hym
              geom_cladelabel(node = 1265, label = "561.4 ± 433.4", color = "#66cc33", fill = "white", fontsize = 6, offset = 16, barsize = NA, geom = 'label', align = T) + #col
              geom_cladelabel(node = 823, label = "284.4 ± 233.5", color = "#443333", fill = "white", fontsize = 6, offset = 18, barsize = NA, geom = 'label', align = T) + #dip
              geom_cladelabel(node = 691, label = "659.6 ± 502.4", color = "#0099cc", fill = "white", fontsize = 6, offset = 25, barsize = NA, geom = 'label', align = T) + #hem
              scale_fill_manual(values = colours) +
              guides(color = "none", fill = guide_legend(override.aes = list(size = 3))) +
              labs(title = "Genome Size (Mb)") +
              theme(legend.key=element_rect(fill=NA),
                    legend.text = element_text(size = 14),
                    legend.title = element_text(size = 15, face = "bold"),
                    plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = -10),
                    plot.margin = unit(c(0,0,0,0),"cm")),
            90)

ggsave("Figure 8_GenSizeMain.pdf", width = 18, height = 16, units = "in", dpi = 300)


###################################################
####### FIGURE 8 - OTHER ORDERS ###################
###################################################

#load tree file and extra data file
tree.Others <- read.tree("COITree_Others.tre")

#reduce data to only species contained within tree
data <- data.in[data.in$Tips %in% tree.Others[4][[1]],]

#set order of Orders for plot legend
data$Order <- factor(data$Order, levels = c("Orthoptera","Plecoptera","Blattodea","Odonata","Ephemeroptera","Phasmatodea","Trichoptera","Thysanoptera","Neuroptera","Megaloptera","Strepsiptera","Siphonaptera","Coleoptera","Diptera","Hemiptera","Hymenoptera","Lepidoptera"))

#set custom colours
colours <- c("#cc0011","#ffaa00","#66cc33","#443333","#0099cc","purple","#cc0011","#ffaa00","#66cc33","#443333","#0099cc","purple")

#Create tree object and attach associated data
###NUMT Count###
g.tree <- ggtree(tree.Others, aes(color = Order), branch.length="none", layout = "fan", open.angle = 5, size = 0.1) %<+% data
rotate_tree(g.tree +
              scale_color_manual(breaks = levels(data$Order), values = colours) +
              geom_fruit(geom = geom_bar,
                         mapping = aes(x = NUMT_Count_Total, fill = Order),
                         orientation = "y",
                         stat = "identity",
                         pwidth = 1,
                         grid.params = list(vline = F, alpha = 0.5),
                         axis.params = list(axis = "x", 
                                            nbreak = 6,
                                            text.size = 7.5, text.angle = 0,
                                            hjust = -0.1)) +
              scale_fill_manual(values = colours) +
              guides(color = "none", fill = guide_legend(override.aes = list(size = 3))) +
              labs(title = "NUMT Count") +
              theme(legend.key=element_rect(fill=NA),
                    legend.text = element_text(size = 14),
                    legend.title = element_text(size = 15, face = "bold"),
                    plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = -10),
                    plot.margin = unit(c(0,0,0,0),"cm")),
            90)

ggsave("Figure 8_NUMTCountOthers.pdf", width = 18, height = 16, units = "in", dpi = 300)

###Genome Size###
g.tree <- ggtree(tree.Others, aes(color = Order), branch.length="none", layout = "fan", open.angle = 5, size = 0.1) %<+% data
rotate_tree(g.tree +
              scale_color_manual(breaks = levels(data$Order), values = colours) +
              geom_fruit(geom = geom_bar,
                         mapping = aes(x = Genome_Asssembly_Size_Mb, fill = Order),
                         orientation = "y",
                         stat = "identity",
                         pwidth = 1,
                         grid.params = list(vline = F, alpha = 0.5),
                         axis.params = list(axis = "x", 
                                            nbreak = 6,
                                            text.size = 7.5, text.angle = 0,
                                            hjust = -0.1)) +
              scale_fill_manual(values = colours) +
              guides(color = "none", fill = guide_legend(override.aes = list(size = 3))) +
              labs(title = "Genome Size") +
              theme(legend.key=element_rect(fill=NA),
                    legend.text = element_text(size = 14),
                    legend.title = element_text(size = 15, face = "bold"),
                    plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = -10),
                    plot.margin = unit(c(0,0,0,0),"cm")),
            90)

ggsave("Figure 8_GenSizeOthers.pdf", width = 18, height = 16, units = "in", dpi = 300)

