Exploring sequence and phylogenetic diversity for a taxonomic group of interest.
In metagenomeFeatures: Exploration of marker-gene sequence taxonomic annotations

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE
)
suppressPackageStartupMessages(library(metagenomeFeatures)) 
suppressPackageStartupMessages(library(ggtree))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(forcats))

Vignette overview

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. As such, 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. A MgDb-class object with the greengenes database version 13.8 85\% similarity OTUs (gg85) is included in the r Biocpkg("metagenomeFeatures") package.

Other databases are available as Bioconductor AnnotationData.

Greengenes Release 13.5: r Biocpkg("greengenesMgDb13.5")
Ribosomal Database Project Release 11.5: r Biocpkg("ribosomaldatabaseproject11.5MgDb")
Silva 128.1: r Biocpkg("silva128.1MgDb")

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 r Biocpkg("metagenomeFeatures") package, the vignette uses r CRANpkg("tidyverse") packages r CRANpkg("dplyr") and r CRANpkg("tidyr") for data manipulation, r CRANpkg("ggplot2") and for data visualization. The r Biocpkg("ggtree") Bioconductor package is used to plot the phylogenetic tree.

library(metagenomeFeatures)
library(dplyr)
library(forcats)
library(ggplot2)
library(ggtree)

Loading and subsetting the database

The gg85 database is loaded using get_gg13.8_85MgDb().
Next the 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. The 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()

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

Taxonomic Analysis

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

## 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 r sum(gamma_df$Count) OTUs assigned to r nlevels(gamma_df$Genus) genera in the Class Gammaproteobacteria. The number of OTUs assigned to specific genera, range from r sort(gamma_df$Count, decreasing = TRUE)[2] to r min(gamma_df$Count) (Fig. \@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 a genus is represented in the database. For example no sequences present in the database are assigned to the genus Shigella and only r escherichia_count is assigned to Escherichia. OTUs only assigned to the family, g__, is the most abundant group, r gamma_df$Count[gamma_df$Genus == "g__"].

 ggplot(gamma_trim_df) + 
    geom_bar(aes(x = Genus, y = Count), stat = "identity") + 
    labs(y = "Number of OTUs") + 
    coord_flip() + 
    theme_bw()

Phylogenetic Diversity

We will use the r Biocpkg("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. \@ref(fig:annoTree)).

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

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. When working with large subsets of the database or the full database, the Bioconductor package r Biocpkg("DECIPHER") provides functions for analyzing sequences in an SQLite database. The MgDb-class uses r Biocpkg("DECIPHER") to format the sequences as an SQLite database. The MgDb-class sequence slot is a connection to the SQLite database, therefore the r Biocpkg("DECIPHER") package can be used to analyze the dataset reference sequence data as well.