if(!"devtools" %in% installed.packages()){
  install.packages("devtools")
}
devtools::install_github("gmteunisse/fantaxtic")

require("fantaxtic")
require("phyloseq")
require("tidyverse")
require("magrittr")
require("ggnested")
require("knitr")
require("gridExtra")

data(GlobalPatterns)
top_asv <- top_taxa(GlobalPatterns, n_taxa = 10)
plot_nested_bar(ps_obj = top_asv$ps_obj,
                top_level = "Phylum",
                nested_level = "Species")

top_nested <- nested_top_taxa(GlobalPatterns,
                              top_tax_level = "Phylum",
                              nested_tax_level = "Species",
                              n_top_taxa = 3, 
                              n_nested_taxa = 3)
plot_nested_bar(ps_obj = top_nested$ps_obj,
                top_level = "Phylum",
                nested_level = "Species")

top_species <- top_taxa(GlobalPatterns,
                        n_taxa = 10, 
                        tax_level = "Species")
top_species$top_taxa %>%
  mutate(abundance = round(abundance, 3)) %>%
  kable(format = "markdown")

top_grouped <- top_taxa(GlobalPatterns,
                        n_taxa = 1,
                        grouping = "SampleType")
top_grouped$top_taxa %>%
  mutate(abundance = round(abundance, 3)) %>%
  kable(format = "markdown")

top_max <- top_taxa(GlobalPatterns,
                        n_taxa = 10,
                        FUN = max)
top_max$top_taxa %>%
  mutate(abundance = round(abundance, 3)) %>%
  kable(format = "markdown")

top_nested <- nested_top_taxa(GlobalPatterns,
                              top_tax_level = "Phylum",
                              nested_tax_level = "Species",
                              n_top_taxa = 3, 
                              n_nested_taxa = 3, 
                              nested_merged_label = "NA and other <tax>")
top_nested$top_taxa %>%
  mutate(top_abundance = round(top_abundance, 3),
         nested_abundance = round(nested_abundance, 3)) %>%
  kable(format = "markdown")

plot_nested_bar(top_nested$ps_obj,
                top_level = "Phylum",
                nested_level = "Species",
                nested_merged_label = "NA and other <tax>",
                legend_title = "Phylum and species") +
  facet_wrap(~SampleType,
             scales = "free_x") +
  labs(title = "Relative abundances of the top 3 species for each of the top 3 phyla") +
  theme(plot.title = element_text(hjust = 0.5, 
                                  size = 8, 
                                  face = "bold"),
        legend.key.size = unit(10, 
                               "points"))

plot_nested_bar(top_nested$ps_obj,
                top_level = "Phylum",
                nested_level = "Species",
                nested_merged_label = "NA and other",
                palette = c(Bacteroidetes = "red", 
                            Proteobacteria = "blue"),
                merged_clr = "black")

# Order samples by the total abundance of Proteobacteria
sample_order <- psmelt(top_nested$ps_obj) %>%
  data.frame() %>%
  
  # Calculate relative abundances
  group_by(Sample) %>%
  mutate(Abundance = Abundance / sum(Abundance)) %>%
  
  # Sort by taxon of interest
  filter(Phylum == "Proteobacteria") %>%
  group_by(Sample) %>%
  summarise(Abundance = sum(Abundance)) %>%
  arrange(Abundance) %>% 
  
  # Extract the sample order
  pull(Sample) %>%
  as.character()

# Plot
plot_nested_bar(top_nested$ps_obj,
                "Phylum",
                "Species",
                sample_order = sample_order,
                nested_merged_label = "NA and other")

# Get the top taxa
top_level <- "Phylum"
nested_level <- "Species"
sample_order <- NULL
top_asv <- top_taxa(GlobalPatterns, n_taxa = 10)

# Create names for NA taxa
ps_tmp <- top_asv$ps_obj %>%
  name_na_taxa()

# Add labels to taxa with the same names
ps_tmp <- ps_tmp %>%
  label_duplicate_taxa(tax_level = nested_level)

# Generate a palette basedon the phyloseq object
pal <- taxon_colours(ps_tmp,
                     tax_level = top_level)

# Convert physeq to df
psdf <- psmelt(ps_tmp)

# Move the merged labels to the appropriate positions in the plot:
# Top merged labels need to be at the top of the plot,
# nested merged labels at the bottom of each group
psdf <- move_label(psdf = psdf,
                   col_name = top_level,
                   label = "Other",
                   pos = 0)
psdf <- move_nested_labels(psdf,
                           top_level = top_level,
                           nested_level = nested_level,
                           top_merged_label = "Other",
                           nested_label = "Other",
                           pos = Inf)

# Reorder samples
if(!is.null(sample_order)){
  if(all(sample_order %in% unique(psdf$Sample))){
    psdf <- psdf %>%
      mutate(Sample = factor(Sample, levels = sample_order))
  } else {
    stop("Error: not all(sample_order %in% sample_names(ps_obj)).")
  }

}

# Generate a base plot
p <- ggnested(psdf,
              aes_string(main_group = top_level,
                         sub_group = nested_level,
                         x = "Sample",
                         y = "Abundance"),
              main_palette = pal) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_nested(theme_light) +
  theme(axis.text.x = element_text(hjust = 1, vjust = 0.5, angle = 90))

# Add relative abundances
p <- p + geom_col(position = position_fill())
p

# Fill in names for NA taxa, including their rank
ps_tmp <- name_na_taxa(top_asv$ps_obj)
tax_table(ps_tmp) %>%
  kable(format = "markdown")

# Leave the rank out and alter the label
ps_tmp <- name_na_taxa(top_asv$ps_obj, 
                       include_rank = F, 
                       na_label = "NA <tax>")
tax_table(ps_tmp) %>%
  kable(format = "markdown")

# Function to plot the colours in a palette
plot_colours <- function(pal){
  pal %>%
    data.frame(name = names(.),
               colour = .) %>%
    ggplot(aes(x = 1,
               y = name,
               fill = colour,
               label = paste(name, "-", colour))) +
    geom_tile() +
    geom_text() +
    scale_fill_identity() +
    scale_x_discrete(expand = c(0,0)) +
    scale_y_discrete(expand = c(0,0)) +
    theme(axis.text = element_blank(),
          axis.title = element_blank(),
          axis.ticks = element_blank(),
          plot.title = element_text(hjust = 0.5,
                                    size = 12))
}

# Get top taxa and generate a palette
pal <- taxon_colours(top_asv$ps_obj, tax_level = "Phylum")
p1 <- plot_colours(pal) +
  ggtitle("Default")

# Generate a palette with a different base_clr
pal2 <- taxon_colours(top_asv$ps_obj, tax_level = "Phylum", base_clr = "blue")
p2 <- plot_colours(pal2) +
  ggtitle("Base colour blue")

# Provide a custom incomplete palette
pal3 <- taxon_colours(top_asv$ps_obj, 
                     tax_level = "Phylum", 
                     palette = c(Cyanobacteria = "blue", 
                                 Bacteroidetes = "pink"))
p3 <- plot_colours(pal3) +
  ggtitle("Incomplete custom palette")
grid.arrange(p1, p2, p3, nrow = 1)

# Label the lowest non-NA level
ps_tmp <- label_duplicate_taxa(top_asv$ps_obj, 
                               tax_level = "Genus")
tax_table(ps_tmp)  %>%
  kable(format = "markdown")

# Use ASVs as ids rather than counts
ps_tmp <- label_duplicate_taxa(top_asv$ps_obj,
                               tax_level = "Genus",
                               asv_as_id = T, 
                               duplicate_label = "<tax> ASV <id>")
tax_table(ps_tmp) %>%
  kable(format = "markdown")

# Turn physeq object into a dataframe
ps_tmp <- name_na_taxa(top_asv$ps_obj)
ps_tmp <- label_duplicate_taxa(ps_tmp, tax_level = "Species")
psdf <- psmelt(ps_tmp)
levels(as.factor(psdf$Phylum))
#> [1] "Actinobacteria" "Bacteroidetes"  "Cyanobacteria"  "Firmicutes"    
#> [5] "Other"          "Proteobacteria"

# Move the other label to the start, and Bacteroidetes to the end
psdf <- move_label(psdf, col_name = "Phylum", label =  "Other", pos = 0)
psdf <- move_label(psdf, col_name = "Phylum", label =  "Bacteroidetes", pos = Inf)
levels(psdf$Phylum)
#> [1] "Other"          "Actinobacteria" "Cyanobacteria"  "Firmicutes"    
#> [5] "Proteobacteria" "Bacteroidetes"

levels(as.factor(psdf$Species))
#>  [1] "Haemophilusparainfluenzae"       "Other"                          
#>  [3] "Streptococcussanguinis"          "Unknown ACK-M1 (Family)"        
#>  [5] "Unknown Bacteroides (Genus) 1"   "Unknown Bacteroides (Genus) 2"  
#>  [7] "Unknown Dolichospermum (Genus)"  "Unknown Neisseria (Genus)"      
#>  [9] "Unknown Providencia (Genus)"     "Unknown Stramenopiles (Order) 1"
#> [11] "Unknown Stramenopiles (Order) 2"
psdf <- move_nested_labels(psdf, 
                           top_level = "Phylum", 
                           nested_level = "Species",
                           top_merged_label = "Other",
                           nested_label = "Unknown Bacteroides 1",
                           pos = Inf)
levels(psdf$Species)
#>  [1] "Unknown ACK-M1 (Family)"         "Unknown Bacteroides (Genus) 1"  
#>  [3] "Unknown Bacteroides (Genus) 2"   "Unknown Dolichospermum (Genus)" 
#>  [5] "Unknown Stramenopiles (Order) 1" "Unknown Stramenopiles (Order) 2"
#>  [7] "Streptococcussanguinis"          "Haemophilusparainfluenzae"      
#>  [9] "Unknown Neisseria (Genus)"       "Unknown Providencia (Genus)"    
#> [11] "Other"

# Generate a base plot
p <- ggnested(psdf,
              aes_string(main_group = top_level,
                         sub_group = nested_level,
                         x = "Sample",
                         y = "Abundance"),
              main_palette = pal) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_nested(theme_light) +
  theme(axis.text.x = element_text(hjust = 1, vjust = 0.5, angle = 90))

# Add relative abundances
p <- p + geom_col(position = position_fill(), width = 0.5)
p

# Create a boxplot instead of a barplot
psdf_rel <- psdf %>%
  group_by(Sample) %>%
  mutate(Abundance = Abundance / sum(Abundance)) %>%
  ungroup() %>%
  filter(Phylum != "Other")
p <- ggnested(psdf_rel,
              aes_string(main_group = top_level,
                         sub_group = nested_level,
                         x = "Phylum",
                         y = "Abundance",
                         grouping = "Species"),
              main_palette = pal,
              main_keys = T) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_nested(theme_light) +
  theme(axis.text.x = element_text(hjust = 1, vjust = 0.5, angle = 90))
p + geom_boxplot(alpha = 0.5)

ps_tmp <- get_top_taxa(GlobalPatterns, n = 10)
ps_tmp <- name_taxa(ps_tmp, label = "Unkown", species = T, other_label = "Other")
fantaxtic_bar(ps_tmp, color_by = "Phylum", label_by = "Species", other_label = "Other")