library(ggtree)
library(ape)
library(treeio)
library(ggplot2)
library(aplot)
library(grafify)

setwd('~/Desktop/diatom_genomes/diatom-genomes/cyclostephanoid_genome_report/')

##### plot the phylogenetic tree #####

tree <- read.tree('busco.prot.concat.partition.rooted.tree')

d <- data.frame(label=c('AJA232-27_Discostella_pseudostelligera',
                        'AJA228-03_Cyclostephanos_tholiformis',
                        'AJA276-08_Praestephanos_triporus',
                        'CCMP1335_Cyclotella_nana',
                        'CCMP332_Cyclotella_cryptica',
                        'CNS00166_Skeletonema_tropicum',
                        'CNS00100_Skeletonema_marinoi',
                        'CNS00243_Skeletonema_costatum'),
                label2=c('Discostella pseudostelligera',
                         'Cyclostephanos tholiformis',
                         'Praestephanos triporus',
                         'Cyclotella nana',
                         'Cyclotella cryptica',
                         'Skeletonema tropicum',
                         'Skeletonema marinoi',
                         'Skeletonema costatum'))

tree2 <- rename_taxa(tree, d, label, label2)

d2 <- data.frame(label=d$label2,
                 habitat=c(rep('freshwater',3),rep('euryhaline',2),'marine','euryhaline','marine'))
d2
  

t1 <- ggtree(tree2, right=T, ladderize=T) %<+% d2 +
  geom_tiplab(size=4, aes(color=habitat), fontface=2) +
  #scale_color_manual(values=c('#56B4E9','black')) +
  #guides(color='none') +
  theme(legend.position='bottom',
        legend.title=element_blank(),
        legend.text=element_text(size=12)) +
  geom_treescale(x=0, y=0) +
  scale_x_continuous(limits=c(0,0.65)) +
  scale_color_grafify()
t1
t2 <- ggtree::rotate(t1, 10)
t2


##### plot genome improvement #####

table <- data.frame(species=c('C. tholiformis', 'C. tholiformis', 
                              'D. pseudostelligera', 'D. pseudostelligera', 
                              'P. triporus', 'P. triporus'),
                    status=c('v1','v2','v1','v2','v1','v2'),
                    num_contigs=c(19362, 914, 5284, 316, 15639, 1854),
                    N50_length=c(6.9, 259.7, 14.8, 194.4, 7.3, 50.8))

table$status <- factor(table$status, levels=c('v1','v2'))

p1 <- ggplot(table, aes(y=N50_length, x=status)) +
  geom_bar(stat='identity', fill='grey50') +
  facet_wrap(~species) +
  theme_bw() +
  theme(panel.grid=element_blank(),
        strip.background=element_rect(fill='white'),
        strip.text=element_text(color='black', size=10, face='bold')) +
  labs(x='', y='N50 length (kb)') +
  geom_text(aes(label=N50_length), vjust=-0.5) +
  scale_y_continuous(limits=c(0,300))
p1

p2 <- ggplot(table, aes(y=num_contigs, x=status)) +
  geom_bar(stat='identity', fill='grey50') +
  facet_wrap(~species) +
  theme_bw() +
  theme(panel.grid=element_blank(),
        strip.text=element_blank(),
        strip.background=element_blank()) +
  labs(x='genome version', y='number of contigs') +
  geom_text(aes(label=num_contigs), vjust=-0.5) +
  scale_y_continuous(limits=c(0,25000))
p2

p3 <- plot_list(p1, p2, nrow=2, labels=c('(b)','(c)'))
p3

p4 <- plot_list(t2, p3, ncol=2, widths=c(0.5,0.5), labels=c('(a)',''))
p4

ggsave('Figure_1.pdf', p4, height=5, width=10)
ggsave('Figure_1.png', p4, height=5, width=10)