knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
MgDb objects for 16S rRNA sequence reference databases can be used to explore sequnce and phylogenetic diversity of taxonomic groups of interest.
In this vignette we will explore the 16S rRNA diversity for the Enterobacteriace family.
The Greengenes 16S rRNA database version 13.8 clusted at the 85\% threshold will be used in this vignette.
MgDb object with this database (
gg85) is included in the
We will first load the MgDb class object, then select the taxa of interest, and finally provide a explor the phylogenetic and sequence diversity and taxonomic composition of the Enterobacteriace family.
The database exploration will consist of four parts.
First we need to load the
The tidyverse packages, dplyr, tidyr, and ggplot are used display the results.
gg85 database is loaded using
gg85 <- get_gg13.8_85MgDb()
Other databases are avilable as Bioconductor AnnotationData.
mgDb_select() function is first 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.
The select function returns single object with a subset of the database taxonomic, sequence, or phylogenetic data, as well as a named list with any two or all three data types.
gamma_16S <- mgDb_select(gg85, type = "all", keys = "Gammaproteobacteria", keytype = "Class")
We want to know how many genera in the Order Enterobacteriaceae there are sequences for in the Greengenes 85% OTU database as well as how many sequences are assigned to each genera. We will first create a dataframe with the desired information then a plot to summarize the results.
## 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)
For the Greengenes database there are a total of
r total_otus OTUs assigned to
r nlevels(gamma_df$Genus) genera in the Class Gammaproteobacteria.
The number of OTUs assigned to specific genera, range from 76 to 1 \@ref(fig:generaCount).
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 the genera is represented in the database.
For example no sequences present in the database are assigned to the genus Shigella and only
r escherichia_count are assigned to Escherichia.
OTUs only assigned to the family,
g__, is the most abundant group,
r no_genus_assignment, many of which are likely from these two genera.
Next we will use the phylogenetic information to evaluate this assumption.
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.
A number of OTUs unassigned at the genus level are in the same clade as OTUs assigned to the Genera Escherichia and a related Genera Salmonella \@ref(fig:annoTree).
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")
mgDb_select function returns the subsetted sequence data as a biostring object.
The sequence data can be used for additional analysis such as, multiple sequencing alignment, primer region extraction, or PCR primer design.
s_info <- devtools::session_info() print(s_info$platform)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.