Set up the workspace

This analysis uses many open source packages that are available on CRAN and Bioconductor.

To run this analysis yourself, download the supplemental material here along with the R markdown file. Save these files in the same directory on your computer.

# load libraries
library(phyloseq)
library(ggplot2) 
library(tidyverse)
library(RColorBrewer)
library(reshape2)
library(dplyr)
library(knitr)
options(digits=2)
library(kableExtra)
library(gridExtra)
library(ggpubr)

# save session info (packages and versions loaded)
session <- sessionInfo()

Create functions to use in this analysis

We created a function called expCompBarPlot to display the expected mock microbial ASVs in color and unexpected contaminants sequences in grayscale.

# Create function to plot bar plots with contaminants in grey scale and expected mock microbial sequences in color

# Create function to plot bar plots with contaminants in grey scale and expected mock microbial sequences in color

 expCompBarPlot <- function(physeq, exp_taxa, title){
   ## physeq - phyloseq object that will be plotted
   ## exp_taxa - taxa that are expected to be in the mock community 
   ## title - title for plot
#set up data_table
  data_table <- as.data.frame(t(physeq@otu_table))

  data_table$reference = FALSE
  data_table$reference[rownames(data_table) %in% exp_taxa] = TRUE
  sample_names <- sample_names(physeq)
  data_table$id <- paste0('ASV_', 1:nrow(data_table))
  dilution_labels <- sample_data(physeq)$Dilutions

  set.seed(444)

# define the colors to use for reference and non-reference OTUs/ASVs
  ref_colors <- brewer.pal(sum(data_table$reference), "Paired")
  other_colors <- sample(grey.colors(5, start = 0.5, end = 0.9), sum(!data_table$reference), replace = TRUE)

# add a color variable to the data table
  data_table$color <- rep(NA, nrow(data_table))
  data_table$color[data_table$reference] <- ref_colors
  data_table$color[!data_table$reference] <- other_colors

# reshape the data table into a ggplot-friendly format, by gathering samples into a single column called "count"

  color_gg <- data_table %>% select(id, sample_names, color) %>% gather("sample", "count", sample_names)
  legend_color <- c(bright = ref_colors[2], dull = other_colors[2])
  data_gg <- data_table %>% gather("sample", "count", sample_names)

  data_gg <- inner_join(data_gg,color_gg)

# create the composition bar plot
  comp_bar <- ggplot(data_gg, aes(x = sample, y = count)) +
    geom_col(aes(fill = color, group = reference, alpha = ifelse(reference, "bright", "dull")), width = 0.7, position = position_fill()) +
    scale_fill_identity(guide = FALSE) +
    scale_alpha_manual(name = "Sequence type",
                     labels = c("expected", "other"),
                     values = c(bright = 1, dull = 1),
                     guide = guide_legend(override.aes = list(fill = c(ref_colors[4], "#AEAEAE")),
                                          keywidth = NULL, keyheight = NULL)) +
    labs(title = title, x = "Sample", y = "Relative Abundance") +
    theme(legend.position = "right", legend.title = element_text(size = 10),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        plot.title = element_text(size = 12))
  return(comp_bar)
 }

# Define functions to evaluate filter performance
eval_filter <- function(physeq, mock_taxa, recovered_otus, removed_otus, filter_method){
  #  % of mock community ASVs correctly classified as mock community ASVs 
  true_neg <- rowSums(recovered_otus[,colnames(recovered_otus) %in% mock_taxa])
  # % of mock community incorrectly classified as mock community ASVs
  false_neg <- rowSums(recovered_otus[,!colnames(recovered_otus) %in% mock_taxa])
  #  identify non-mock community ASVs correctly classified as not belonging to mock   community
  true_pos <- rowSums(removed_otus[,!colnames(removed_otus) %in% mock_taxa])
  #  identify mock community ASVs incorrectly classified as not belonging to mock   community
  false_pos <- rowSums(removed_otus[,colnames(removed_otus) %in% mock_taxa])

  profile <- rbind(false_neg, false_pos,true_neg,true_pos)

  long_profile <- melt(data = profile, 
                 id.vars = rownames(), 
                 variable.name = colnames(), 
                 value.name = "Abundance"
                )
  names(long_profile)[names(long_profile)=="Var1"] <- "SequenceClass"
  customPalette <- c('#969696','#bdbdbd', '#1B9E77', '#D95F02')

  # Figures
  classificationPlot <- ggplot(long_profile, aes(x = Var2, y = Abundance)) + 
    geom_col(aes(fill = SequenceClass), width = 0.7, position = position_fill())  +
    scale_fill_manual(values=customPalette) + theme(text = element_text(size=12)) +
    labs(x = "Sample", y = 'Proportion of Reads') +
      ggtitle(paste0('Sequence classification for \n  ', filter_method)) + 
      theme(legend.position = "right", legend.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 12),
        plot.title = element_text(size = 12))

  compositionPlot <- expCompBarPlot(physeq,mock_taxa,paste0('Recovered profile after \n ',filter_method))
  # plot figures
  grid.arrange(compositionPlot[[1]], compositionPlot[[2]],classificationPlot, ncol = 3)
}

filter_results <- function(physeq_original, physeq_filter, physeq_removed, mock_taxa,filter_method){
  ## physeq_original <- the original phyloseq object containing actual counts
  ## physeq_filter <- the filtered phyloseq object containing actual counts
  ## physeq_removed <- a phyloseq object of removed actual counts
  ## mock_taxa <- list of mock taxa ASV sequences
  
  # extract out original normalized otu table
  physeq_norm <- transform_sample_counts(physeq_original,function(x) 100* x/sum(x))
  
  original_otus <- as.matrix(as.data.frame(physeq_original@otu_table))
  # identify the orginal proportion of contaminants
  contaminants_original <- rowSums(original_otus[,!colnames(original_otus) %in% mock_taxa])
  # identify the orginal proportion of mock community ASVs
  mock_original <- rowSums(original_otus[,colnames(original_otus) %in% mock_taxa])
  
  # Normalize the filtered physloseq object to relative abundance (each sample sums to 100)
  physeq_filter_norm <- transform_sample_counts(physeq_filter,function(x) 100* x/sum(x))
  
  # subset out the otu table of recovered otus (otus that are kept)
  recovered_otus <- as.matrix(as.data.frame(physeq_filter@otu_table))
  
  # create a subset of removed otus for evaluation
  removed_otus <- as.matrix(as.data.frame(physeq_removed@otu_table))
  
  #  % of mock community ASVs correctly classified as mock community ASVs 
  true_neg <- rowSums(recovered_otus[,colnames(recovered_otus) %in% mock_taxa])
  # % of mock community incorrectly classified as non-mock community ASVs
  false_neg <- rowSums(recovered_otus[,!colnames(recovered_otus) %in% mock_taxa])
  #  identify non-mock community ASVs correctly classified as not belonging to mock   community
  true_pos <- rowSums(removed_otus[,!colnames(removed_otus) %in% mock_taxa])
  #  identify mock community ASVs incorrectly classified as not belonging to mock   community
  false_pos <- rowSums(removed_otus[,colnames(removed_otus) %in% mock_taxa])
  
  sensitivity <- true_pos/(true_pos + false_neg)
  specificty <- true_neg/(true_neg + false_pos)
  accuracy <- (true_pos + true_neg) / (false_pos + true_pos + false_neg + true_neg)
  prevalence <- (true_pos + false_neg) / (false_pos + true_pos + false_neg + true_neg)

  ## proportion of contaminants removed (of all total contaminant ASVs)
  contaminants_removed = (rowSums(removed_otus[,!colnames(removed_otus) %in% mock_taxa])/ contaminants_original) * 100
    ## proportion of mock removed (of all total mock ASVs)
  mock_ASVs_removed = (rowSums(removed_otus[,colnames(removed_otus) %in% mock_taxa])/ mock_original) * 100
  
  ## total amount of conatminants remaining in ne
  #contaminants_remaining =  rowSums(recovered_otus[,!colnames(recovered_otus) %in% mock_taxa])
  # calculate alpha diverity and summary of abundances
  diversity <- estimate_richness(physeq_filter, measures = c('Observed','Shannon','InvSimpson'))
  rel_abundance <- as.data.frame(physeq_filter_norm@otu_table)
  mock_abundance <- rel_abundance[, colnames(rel_abundance) %in% mock_taxa]
  total_mock_abundance <- rowSums(mock_abundance)
  con_abundance <- rowSums(rel_abundance[,!colnames(rel_abundance) %in% mock_taxa])
  # return results
  results <- cbind(contaminants_removed, mock_ASVs_removed,con_abundance, total_mock_abundance, diversity, mock_abundance,sensitivity , specificty, accuracy, prevalence, true_pos, true_neg, false_pos, false_neg)
  # add filter_method to results table
    results <- results %>%
    mutate(method =  filter_method) %>%
    mutate(sample_names = rownames(results))

  return(results)
  }

summary_table <- function(results, filter_method){
  
  caption_text = paste0(filter_method,' summary')
  drop_col <- c('true_neg','true_pos','false_neg','false_pos','method')
  results <- results %>% select(-one_of(drop_col))
  kable(t(results), digits = 2, caption = caption_text) %>%
    kable_styling(bootstrap_options = c("striped", "hover"), full_width = F) %>%
    group_rows("Percent Removed", 1,2) %>%
    group_rows("Percent Remaining after Contaminant Removal", 3,4) %>%
    group_rows("Diversity Estimates after Contaminant Removal", 5,7) %>% 
    group_rows("Mock Abundances after Contaminant Removal", 8,ncol(results) - 2)

}
# define function for evaluateing abundance filtering
abundance_filter_eval <- function(physeq, abundance_filter){
  # normalize physeq to 100
  physeq_norm <- transform_sample_counts(physeq,function(x) 100* x/sum(x))
  # subset out the otu table, name rec_otu for recovered otus (otus that are kept)
  # name relative abundance otu table rec_otu_rel
  rec_otu_rel <- as.matrix(as.data.frame(physeq_norm@otu_table))
  rec_otu <- as.matrix(as.data.frame(physeq@otu_table))
  
  # remove low abundance sequences per sample
  rec_otu_rel[rec_otu_rel <= abundance_filter] <- 0
  rec_otu[rec_otu_rel == 0] <- 0
  ps_filt <- physeq
  otu_table(ps_filt) <- otu_table(rec_otu, taxa_are_rows = FALSE)

  # create a subset of removed otus for evaluation
  rem_otu_rel <- as.matrix(as.data.frame(physeq_norm@otu_table))
  rem_otu <- as.matrix(as.data.frame(physeq@otu_table))
  rem_otu[rem_otu_rel > abundance_filter] <- 0
  ps_rem <- physeq
  otu_table(ps_rem) <- otu_table(rem_otu, taxa_are_rows = FALSE)

  # plot results
  eval_filter(ps_filt, mock_taxa, rec_otu,rem_otu, paste0('Abundance less than ', abundance_filter, '% removed'))
  
  # return results
  results <- filter_results(physeq, ps_filt, ps_rem, mock_taxa,paste0('Abundance filter, ',abundance_filter) )
  return(results)
  }

Load the dataset

This analysis assumes that you have the RData file and R Markdown file saved in the same directory.

## Load the dataset
load("mockDilutionsPrep.RData")

This work outlines each step required for evaluating the impact of contaminants on 16S rRNA gene experiments with samples that have varying starting material. The conclusions from this analysis is that contaminant ASVs increase with decreasing starting material and lead to distorted microbial community profiles and overestimated alpha diversity measurements.

While contaminants are prevalent in low microbial biomass samples, and potentially unavoidable, all hope is not lost. There are a few computational approaches to identify and remove contaminant sequences from 16S rRNA gene sequencing data. To see how we evaluated these different methods for identifying and removing contaminants from this data set, see: * Removing contaminant ASVs with decontam * Removing contaminant ASVs with SourceTracker * Evaluating SourceTracker results

This RMarkdown file provides the code for evaluating how unrecommended methods - removing all sequences that appear in a blank control and an abundance filter would perform for removing contaminants.

Create a summary of results before contaminant removal

# calculate alpha diverity and summary of abundances
  diversity <- estimate_richness(mock_ps, measures = c('Observed','Shannon','InvSimpson'))

## Warning in estimate_richness(mock_ps, measures = c("Observed", "Shannon", : The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.

  rel_abundance <- as.data.frame(mock_ps_norm@otu_table)
  mock_abundance <- rel_abundance[, colnames(rel_abundance) %in% mock_taxa]
  total_mock_abundance <- rowSums(mock_abundance)
  con_abundance <- rowSums(rel_abundance[,!colnames(rel_abundance) %in% mock_taxa])
  
  ## Add true positive/true negative ect for 'uncorrected plot'
  mock_otus <- as.matrix(as.data.frame(mock_ps@otu_table))
  #  % of mock community ASVs correctly classified as mock community ASVs 
  true_neg <- rowSums(mock_otus[,colnames(mock_otus) %in% mock_taxa])
  # % of mock community incorrectly classified as mock community ASVs
  false_neg <- rowSums(mock_otus[,!colnames(mock_otus) %in% mock_taxa])
  
  # true_pos and false_pos do not exist  

  
  # return results
  results_original <- cbind(con_abundance, total_mock_abundance, diversity, mock_abundance, true_neg, false_neg)
  # add filter_method to results table
    results_original <- results_original %>%
    mutate(method =  'Original Data') %>%
    mutate(sample_names = rownames(results_original))
    
summary_table(results_original,'Original results - no contaminant removal')

## Warning: Unknown columns: `true_pos`, `false_pos`