Exploring sequence and phylogenetic diversity for a taxonomic group of interest.

  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. A MgDb object with this database (gg85) is included in the metagenomeFeatures package. 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.

  1. Loading and subsetting the database.
  2. Taxonomic analysis.
  3. Phylogenetic diversity.
  4. Sequence diversity.

Loading and subsetting the database

First we need to load the metagenomeFeatues package. The tidyverse packages, dplyr, tidyr, and ggplot are used display the results.


The gg85 database is loaded using get_gg13.8_85MgDb().

gg85 <- get_gg13.8_85MgDb()

Other databases are avilable as Bioconductor AnnotationData.

Next the 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")

Taxonomic Analysis

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() + 

Phylogenetic Diversity

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).


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")

Sequence Diversity

The 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.


System Information

s_info <- devtools::session_info()

Package Versions


Try the metagenomeFeatures package in your browser

Any scripts or data that you put into this service are public.

metagenomeFeatures documentation built on Nov. 1, 2018, 4:31 a.m.