16S rRNA reference databases are most often used to annotate OTUs or sequence variants from marker-gene survey datasets.
The reference databases are also used to study the diversity of specific taxonomic groups.
MgDb-class objects for 16S rRNA sequence reference databases can be used to explore the sequence and phylogenetic diversity of taxonomic groups of interest.
In this vignette, we will explore the 16S rRNA diversity of the Gammaproteobaceria class using the Greengenes 16S rRNA database version 13.8 clustered at 85% similarity.
MgDb-class object with the greengenes database version 13.8 85% similarity OTUs (
gg85) is included in the metagenomeFeatures package.
Other databases are available as Bioconductor AnnotationData.
We will first load the
MgDb-class object, then select the taxa of interest, and finally explore the Gammaproteobaceria class phylogenetic diversity, and taxonomic composition.
In addition to using the metagenomeFeatures package, the vignette uses tidyverse packages dplyr and tidyr for data manipulation, ggplot2 and for data visualization.
The ggtree Bioconductor package is used to plot the phylogenetic tree.
library(metagenomeFeatures) library(dplyr) library(forcats) library(ggplot2) library(ggtree)
gg85 database is loaded using
mgDb_select() function is used to subset the database, the
key arguments is used to define the taxonmic group(s) of interest and
keytype is used to define the taxonomic level for the group(s) of interest.
With the subsetted database you can analyze the taxonomy, sequences, and phylogeny for the taxonomic group of interest.
mgDb_select() function returns a subset of the database taxonomic, sequence, or phylogenetic data, as well as a named list with any two or all three data types.
gg85 <- get_gg13.8_85MgDb() #> Found more than one class "phylo" in cache; using the first, from namespace 'metagenomeFeatures' #> Also defined by 'tidytree' #> Found more than one class "phylo" in cache; using the first, from namespace 'metagenomeFeatures' #> Also defined by 'tidytree' #> Found more than one class "phylo" in cache; using the first, from namespace 'metagenomeFeatures' #> Also defined by 'tidytree' #> Found more than one class "phylo" in cache; using the first, from namespace 'metagenomeFeatures' #> Also defined by 'tidytree' gamma_16S <- mgDb_select(gg85, type = "all", keys = "Gammaproteobacteria", keytype = "Class")
We want to know how many genera are in the Class Gammaproteobacteria and how many 85% similarity OTUs are in each genera in the Greengenes 13.8 85% database. We first create a data frame with the desired information then plot the results.
## Per genus count data gamma_df <- gamma_16S$taxa %>% group_by(Genus) %>% summarise(Count = n()) %>% ungroup() %>% mutate(Genus = fct_reorder(Genus, Count)) #> `summarise()` ungrouping output (override with `.groups` argument) ## Count info for text 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)
For the Greengenes database there are a total of 247 OTUs assigned to 56 genera in the Class Gammaproteobacteria.
The number of OTUs assigned to specific genera, range from 9 to 1 (Fig. 1).
As this database is preclustered to 85% similarity the number of OTUs per genus is more of an indicator of genera diversity than how well a genus is represented in the database.
For example no sequences present in the database are assigned to the genus Shigella and only 1 is assigned to Escherichia.
OTUs only assigned to the family,
g__, is the most abundant group, 154.
ggplot(gamma_trim_df) + geom_bar(aes(x = Genus, y = Count), stat = "identity") + labs(y = "Number of OTUs") + coord_flip() + theme_bw()
We will use the ggtree package to visualize the phylogenetic tree and annotate the tips with Genera information. Of the genera with more than 3 representative OTUs Stenotrophomonas is the only genus with unclassified and other OTUs within the clade (Fig. 2).
genus_lab <- paste0("g__", c("","Stenotrophomonas", "Pseudomonas")) genus_anno_df <- gamma_16S$taxa %>% group_by(Genus) %>% mutate(Count = n()) %>% ungroup() %>% mutate(Genus_lab = if_else(Count > 3, Genus, "")) ggtree(gamma_16S$tree) %<+% genus_anno_df + ## Add taxa name for unlabeled Stenotrophomonas branch geom_tippoint(aes(color = Genus_lab), size = 3) + scale_color_manual(values = c("#FF000000","darkorange", "blue", "green", "tan", "black")) + theme(legend.position = "bottom") + labs(color = "Genus")