bracken_plot.Rmd

---
title: "bracken_plot"
output: html_document
date: "2023-07-28"
---

# load required libraries
```{r}
# R version 4.2.3
library(stats) # v4.2.3
library(tidyverse) # v2.0.0
library(svglite) # v2.1.1
```

# create plotting function
```{r}
plot_bracken <- function(file, top = NA, pal, otherColor){
    
    # file is the bracken data
    # top is an integer, trims the data to the number of top taxa (by median) in a given level for better plotting
    ## top = NA is the default and will plot all taxa in a given level
    # pal is a palette of hex colors
    # otherColor is the color used to represent taxa below median abundance threshold
    
    levdat <-
      dplyr::select(file,
                    -contains(c("_num","taxonomy_"))) |>
      tidyr::pivot_longer(cols = !name,
                          names_to = "sample",
                          values_to = "fraction") |>
      dplyr::mutate(sample = str_remove_all(sample, "_frac"))
    
    level <-
      dplyr::select(file,
                    contains("taxonomy_lvl")) |>
      unique()
    
    newpal <-
      paste0("#",
             stringr::str_split(pal, ",") |> unlist()
      )
    
    if (level == "K") {
        label <- "Kingdom"
    } else if (level == "P") {
        label <- "Phylum"
    } else if (level == "C") {
        label <- "Class"
    } else if (level == "O") {
        label <- "Order"
    } else if (level == "F") {
        label <- "Family"
    } else if (level == "G") {
        label <- "Genus"
    } else if (level == "S") {
        label <- "Species"
    } else if (level == "S1") {
        label <- "Strain"
    } else {
        message("Bracken taxonomy_lvl needs to be K, P, C, O, F, G, S, or S1")
    }
    
    if (!is.na(top) && top < length(unique(levdat$name))) {
        
        top_select <-
          stats::aggregate(x = levdat$fraction,
                           by = list(levdat$name),
                           FUN = median) |>
          purrr::set_names(c("group","value")) |>
          dplyr::arrange(dplyr::desc(value)) |>
          dplyr::slice(1:top) |>
          dplyr::pull(var = group)
        
        levdat <- dplyr::filter(levdat, levdat$name %in% top_select)
        
        fill_unk <-
          stats::aggregate(x = levdat$fraction,
                           by = list(levdat$sample),
                           FUN = sum) |> 
          purrr::set_names(c("sample","sum")) |> 
          dplyr::mutate_if(is.numeric, round, 3)
        
        fill_unk$residual <- 1 - fill_unk$sum
        
        for (i in 1:nrow(fill_unk)) {
            
            levdat <- rbind(levdat, c("other", fill_unk[i,1], fill_unk[i,3]))
            
        }
        
        levdat <- transform(levdat, fraction = as.numeric(fraction))
        levdat$name <- factor(levdat$name, levels = unique(levdat$name))
        
        fill_vals <- c(rep_len(newpal, nrow(unique(levdat[1]))-1), otherColor)
        
    } else {
        
        fill_vals <- c(rep_len(newpal, nrow(unique(levdat[1]))))
        
    }
    
    print(nrow(levdat))
    
    ggplot2::ggplot(data = levdat) +
      ggplot2::geom_bar(mapping = ggplot2::aes(x = sample,
                                                 y = fraction,
                                                 fill = name),
                        position = "fill",
                        stat = "identity",
                        width = 0.75) + 
      ggplot2::scale_fill_manual(paste(label),
                                 values = fill_vals) +
      ggplot2::theme_classic() +
      ggplot2::scale_x_discrete("") +
      ggplot2::scale_y_continuous("Relative Abundance") +
      ggplot2::theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 1, size = 12),
              axis.text.y = element_text(size = 10),
              axis.title.x = element_text(),
              axis.title.y = element_text(size = 13),
              legend.position = "right")
    
    
    
}
```

# create plot
```{r}
# specify coltypes to prevent F character (family) from being read as logical FALSE
notLogical <- cols(taxonomy_lvl = col_character())
p1 <- plot_bracken(file = readr::read_tsv(file = "example_data_bracken/example_bracken_merged_F.txt",
                                          col_names = TRUE,
                                          col_types = notLogical),
                   top = 10,
                   pal = "dcdcdc,7f0000,006400,708090,808000,483d8b,000080,9acd32,8b008b,66cdaa",
                   otherColor = "gray")

p1 
```

# save as svg, pdf
```{r}
ggsave(filename = "plots/bracken_plot_C.svg",
       plot = p1, units = "in",
       width = 10,
       height = 6)

ggsave(filename = "plots/bracken_plot_C.pdf",
       plot = p1, units = "in",
       width = 10,
       height = 6)
```