## ----setup, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----message = FALSE-----------------------------------------------------
library(metagenomeFeatures) 
library(tidyverse)

## ------------------------------------------------------------------------
gg85 <- get_gg13.8_85MgDb()

## ------------------------------------------------------------------------
gamma_16S <- mgDb_select(gg85, 
                          type = "all", 
                          keys = "Gammaproteobacteria", 
                          keytype = "Class")

## ------------------------------------------------------------------------
## Per genus count data
gamma_df <- gamma_16S$taxa %>% 
    group_by(Genus) %>% 
    summarise(Count = n()) %>% 
    ungroup() %>% 
    mutate(Genus = fct_reorder(Genus, Count)) 

## Count info for text 
total_otus <- sum(gamma_df$Count)
no_genus_assignment <- gamma_df$Count[gamma_df$Genus == "g__"]
escherichia_count <- gamma_df$Count[gamma_df$Genus == "g__Escherichia"]

## excluding unassigned genera and genera with only one assigned sequence
gamma_trim_df <- gamma_df %>% filter(Genus != "g__", Count > 1)

## ----generCount, fig.cap = "Number of seqeunces assigned to genera in the Class Gammaproteobacteria. Only Genera with more than one assigned sequence are shown.", fig.width = 6, fig.height = 6----
 ggplot(gamma_trim_df) + 
    geom_bar(aes(x = Genus, y = Count), stat = "identity") + 
    labs(y = "Number of OTUs") + 
    coord_flip() + 
    theme_bw() 

## ----annoTree, fig.cap = "Phylogenetic tree for Gammaproteobacteria class OTU representative sequences. OTUs not classified to the genus level are annotated with orange dots, and _Escherichia_ and _Salmonella_ OTUs are annotated with blue and green dots respectively.", fig.height = 8, fig.width = 6, message = FALSE----
library(ggtree)

genus_anno_df <- gamma_16S$taxa %>% 
    group_by(Genus) %>% 
    mutate(Count = n()) %>% 
    ungroup() %>% 
    mutate(Genus_lab = if_else(Genus %in% paste0("g__", c("","Escherichia", "Salmonella")), Genus, ""))

ggtree(gamma_16S$tree) %<+% 
    genus_anno_df + 
    geom_tippoint(aes(color = Genus_lab), size = 3) + 
    scale_color_manual(values = c("#FF000000","darkorange","blue", "darkgreen")) + 
    theme(legend.position = "bottom") + 
    labs(color = "Genus")

## ------------------------------------------------------------------------
gamma_16S$seq

## ------------------------------------------------------------------------
s_info <- devtools::session_info()
print(s_info$platform)

## ------------------------------------------------------------------------
s_info$packages