Data Preprocessing - mucosal-immunology-lab/microbiome-analysis GitHub Wiki
- Data preprocessing
Please refer to the relevant pipeline for processing of raw sequencing reads:
- 16S rRNA amplicon sequencing: DADA2 pipeline
- Shotgun metagenomics: Sunbeam pipeline
After completing the steps in the DADA2 pipeline, you should have three input files for generation of a phyloseq object (a container object to hold your taxa counts, information, and sample data):
-
seqtab_nochim.rds: the counts data for each ASV after the removal of chimeras. -
taxonomy_species.rds: the taxonomic lineage information for each ASV, as assigned using the SILVA database. -
tree.rds: a de novo phylogenetic tree generated in the DADA2 pipeline that encompasses all of your ASVs.
Firstly, load in your input files. You may wish to change the sample names (i.e. column names) of your seqtab_nochim.rds file, but be sure that the sample names for your metadata at the next step match (or are changed to match).
The highest_ClassID() function is provided as a script here
# Load bacterial count table seqtab_nochim.rds and corresponding taxonomic classification taxonomy_species.rds
bact_seqtab <- t(data.frame(readRDS(here::here('data', 'dada2', 'seqtab_nochim.rds'))))
bact_taxonomy <- data.frame(readRDS(here::here('data', 'dada2', 'taxonomy_species.rds')))
bact_tree <- readRDS(here::here('data', 'dada2', 'tree.rds'))
# Add a taxonomy column 'ID' to the taxonomy for visualisation purposes (highest level of classification and a unique ASV number)
bact_taxonomy$ID <- highest_ClassID(bact_taxonomy)Your metadata should be prepared with columns in the same order as your OTU table, with identical column names.
To ensure everything is in the correct order, we can also match the sequencing data with the metadata and subset the dataset to keep only bact_seqtab samples with corresponding metadata.
# Keep only metadata rows that match samples in the bact_seqtab
metadata <- metadata[which(rownames(metadata) %in% colnames(bact_seqtab)),]
# Now, keep only bact_seqtab samples with corresponding metadata
bact_seqtab <- bact_seqtab[, which(colnames(bact_seqtab) %in% rownames(metadata))]
# Take a look at the dimensions of the remaining seqtab
dim(bact_seqtab)phyloseq is an R package to import, store, analyse, and graphically display complex phylogenetic sequencing data that has already been processed.
An example for its creation is given below.
# Get data ready for phyloseq creation
otu_mat <- as.matrix(bact_seqtab[, sort(colnames(bact_seqtab))])
tax_mat <- as.matrix(bact_taxonomy)
samples_df <- metadata[colnames(bact_seqtab),]
# Use data to create the individual phyloseq object elements
OTU <- otu_table(otu_mat, taxa_are_rows = TRUE)
TAX <- tax_table(tax_mat)
samples <- sample_data(samples_df)
sample_names(samples) <- colnames(otu_mat)
# Import data into a phyloseq object (omit tree if using shotgun)
bact_data_raw <- phyloseq(OTU, TAX, samples, phy_tree(bact_tree$tree))
# Remove temp files
rm(list = c('otu_mat', 'tax_mat', 'samples_df', 'OTU', 'TAX'))Given that much of the quality control and filtering of our shotgun metagenomic data has already taken place using the Sunbeam pipeline, and we also have taxonomic assignments using the combined power of the Kraken2 and Bracken tools, we do not require as much local pre-processing of data in R before downstream analysis.
The first step is to load in your all_samples_kraken2.csv or all_samples_bracken.csv file, and extract the OTU ID and consensus lineage information. This will help set you up for preparation of a phyloseq container object to hold your data.
We provide a custom function here (kraken2_preprocess.R for producing a taxonomy table and OTU table with unique identifiers from your output file. You may wish to correct the column names after this function however.
The kraken2_preprocess() function will return a list with two elements: firstly an object called kraken2_tax_table which holds the taxonomy table, and secondly one called kraken2_otu_table which holds the OTU read count information (these can be individually extracted using the typical $ notation).
# Read in and pre-process raw kraken2 data
bact_kraken2_input <- kraken2_preprocess(filepath = here::here('input', 'shotgun_data', 'all_samples_kraken2.csv'))
# Extract individual elements
kraken2_otu_table <- bact_kraken2_input$kraken2_otu_table
kraken2_tax_table <- bact_kraken2_input$kraken2_tax_tableYour metadata should be prepared with row names in the same order as the sample names in your OTU table, with identical names.
This will ensure your metadata will be matched to the count data.
phyloseq is an R package to import, store, analyse, and graphically display complex phylogenetic sequencing data that has already been processed.
An example for its creation is given below.
# Get data ready for phyloseq creation
otu_mat <- as.matrix(kraken2_otu_table)
tax_mat <- as.matrix(kraken2_tax_table)
samples_df <- sample_metadata
# Use data to create the individual phyloseq object elements
OTU <- otu_table(otu_mat, taxa_are_rows = TRUE)
TAX <- tax_table(tax_mat)
samples <- sample_data(samples_df)
sample_names(samples) <- colnames(otu_mat)
# Import data into a phyloseq object
bact_data_raw <- phyloseq(OTU, TAX, samples)
# Retain only bacterial taxa (you could also choose another kingdom here)
bact_data_raw <- subset_taxa(bact_data_raw, kingdom == 'Bacteria')
# Remove temp files
rm(list = c('otu_mat', 'tax_mat', 'samples_df', 'OTU', 'TAX'))We will use the R decontam package to remove ASVs/taxa that are more prevalent in our negative controls than in real samples.
The isContaminant function also offers other methods for detection and removal of contaminants, and these can be viewed in the package documentation. This method uses the prevalence (presence/absence across samples) of each sequence feature in true positive samples compared to those in negative controls to identify contaminants.
We will first remove super low read count samples, and then summarise our negative control samples in a logical variable called is_neg, where TRUE indicates negative controls and FALSE indicates real samples.
We can then create a data.frame object to inspect which features have been identified as contaminants.
# Remove very low count samples
minreadsThreshold <- 250
nsamples(bact_data_raw)
bact_data_raw <- prune_samples(colSums(otu_table(bact_data_raw)) > minreadsThreshold, bact_data_raw)
nsamples(bact_data_raw)
# Re-identify the negative controls
sample_data(bact_data_raw)$is_neg <- sample_data(bact_data_raw)$group == 'Negative_control'
# Determine contaminants based on prevalence in negative controls
contamdf_prev <- isContaminant(bact_data_raw, method = 'prevalence', neg = 'is_neg', threshold = 0.5)
table(contamdf_prev$contaminant)
# See which ASVs were contaminants
(contaminants <- data.frame(tax_table(bact_data_raw)[contamdf_prev$contaminant, 8], row.names = NULL))Now that we have this data, we can visually inspect the identified contaminants using a presence/absence scatter plot, with the prevalence in negative controls along the x-axis, and prevalence in true samples along the y-axis.
Once we have made this plot, we can remove both the contaminant features and the negative control samples from our phyloseq object. Then we will be ready for filtering and normalisation. In this way, the contaminant read counts won't be counted towards the reads we are actually interested in.
# Make phyloseq object of presence-absence in negative controls and true samples
ps_pa <- transform_sample_counts(bact_data_raw, function(abund) 1 * (abund > 0))
ps_pa_pos <- phyloseq::subset_samples(ps_pa, is_neg == FALSE)
ps_pa_neg <- phyloseq::subset_samples(ps_pa, is_neg == TRUE)
# Make data.frame of prevalence in positive and negative samples
df_pa <- data.frame('pa_pos' = taxa_sums(ps_pa_pos),
'pa_neg' = taxa_sums(ps_pa_neg),
'contaminant' = contamdf_prev$contaminant)
rownames(df_pa) <- tax_table(ps_pa)[,8]
df_pa$contaminant <- ifelse(df_pa$pa_pos == 0 & df_pa$pa_neg > 0, TRUE, df_pa$contaminant)
# Visualise the contaminants
(decontam_plot <- ggplot(data = df_pa, aes(x = pa_neg, y = pa_pos)) +
geom_point(aes(colour = contaminant, fill = contaminant),
shape = 21, alpha = 0.5,
position = position_jitter(width = 0.05, height = 0.05)) +
scale_fill_npg(name = 'Contaminant') +
scale_colour_npg(name = 'Contaminant') +
labs(title = 'Decontam results',
x = 'Prevalence (Negative Controls)',
y = 'Prevalence (True Samples)')
)
ggsave(here::here('figures', 'contaminant_scatter.pdf'), decontam_plot,
width = 16, height = 12, units = 'cm')
# Remove contaminants
bact_data_clean <- prune_taxa(!df_pa$contaminant, bact_data_raw)
# Keep only samples
bact_data_clean <- phyloseq::subset_samples(bact_data_clean, BT_group != '5_Negative_control')
saveRDS(bact_data_clean, here::here('output', 'bact_data_clean_decontam.rds'))An example of the contaminant scatter plot is shown here:
At this stage, we should now have a single, decontaminated, unified phyloseq object that contains all of our count data, taxonomic information, and sample metadata:
-
bact_data_raw: the rawphyloseqobject before filtering and normalisation (can be derived from either 16S rRNA amplicon or shotgun metagenomic sequencing data).
We can now perform some filtration steps to remove samples with low read counts, and also taxa with very few reads or those that are only present in a small number of samples.
# Remove samples with less than [minreadsThreshold] reads
minreadsThreshold <- 10000 # originally set to 10000
bact_data_samples <- prune_samples(colSums(otu_table(bact_data_raw)) > minreadsThreshold, bact_data_raw)
# Filter out taxa based on a minimum number of reads [detectionThreshold] and prevalence [prevalenceThreshold]
detectionThreshold <- 0 # originally set to 0 (i.e. 0 reads minimum)
prevalenceThreshold <- 0.1 # originally set to 0.1 (i.e. found in 10% of samples)
bact_data_filtered <- core(bact_data_samples,
detection = detectionThreshold,
prevalence = prevalenceThreshold,
include.lowest = FALSE)
# Remove samples with less than [minreadsThreshold] reads, after the prevalence filtering.
minreadsThreshold <- 5000 # originally set to 5000
bact_data_filtered <- prune_samples(colSums(otu_table(bact_data_filtered)) > minreadsThreshold, bact_data_filtered)
# Generate a summary of filtration
data_summary <- data.frame('Summary' = c('OTUs/genera', 'Samples'),
'Bacteria_before_filtering' = dim(otu_table(bact_data_raw)),
'Bacteria_after_filtering' = dim(otu_table(bact_data_filtered)))
# Print the table as a kable object
kable(data_summary) %>%
kable_styling(bootstrap_options = 'striped', full_width = FALSE)
# Save the filtered phyloseq to an .rds file
saveRDS(bact_data_filtered, here::here('output', 'bact_data_filtered.rds'))We can add diversity index information at this time, so it will be incorporated into any normalised datasets. The Shannon, Chao1, and Simpson indices are all provided as part of the phyloseq package. For the Hill diversity index, you will need to install the hillR package.
The rotate_df() function can be found in the code section here
# Estimate Chao1 index (richness)
sample_data(bact_data_filtered)$Chao1 <- estimate_richness(bact_data_filtered, split = TRUE, measures = c('Chao1'))$Chao1
# Estimate Shannon index and retrieve the number of ASVs (observed) for bacteria
sample_data(bact_data_filtered)$Shannon <- estimate_richness(bact_data_filtered, split = TRUE, measures = c('Shannon'))$Shannon
# Estimate Simpson index
sample_data(bact_data_filtered)$Simpson <- estimate_richness(bact_data_filtered, split = TRUE, measures = c('Simpson'))$Simpson
# Estimate Hill diversity with a phylogenetic tree
library(hillR)
sample_data(bact_data_filtered)$Hill_phylo <- hill_phylo(comm = data.frame(otu_table(bact_data_filtered)) %>% rotate_df(), tree = bact_data_filtered@phy_tree, q = 0)
# Or Hill diversity without a phylogenetic tree
sample_data(bact_data_filtered)$Hill_taxa <- hill_taxa(comm = data.frame(otu_table(bact_data_filtered)) %>% rotate_df(), q = 0)MMUPHin is an R package implementing meta-analysis methods for microbial community profiles. It has interfaces for: a) covariate-controlled batch and study effect adjustment, b) meta-analytic differential abundance testing, and meta-analytic discovery of c) discrete (cluster-based) or d) continuous unsupervised population structure.
Install the tool using Bioconductor:
BiocManager::isntall('MMUPHin')We will first extract the OTU counts table and metadata as data.frame objects from our phyloseq object.
# Recover phyloseq object
bact_data_filtered <- readRDS('output', 'bact_data_filtered.rds')
# Prepare data for batch correction
bact_abd <- data.frame(bact_data_filtered@otu_table)
bact_meta <- data.frame(bact_data_filtered@sam_data)We can now correct the data according to batch identity, as well as study condition or treatment group for example (i.e. some other covariate). We can then replace the OTU table component of our phyloseq object.
# Batch correct the OTU data, taking group into account as a covariate
fit_adjust_batch <- adjust_batch(feature_abd = bact_abd,
batch = 'batch',
covariates = 'group',
data = bact_meta,
control = list(verbose = TRUE))
# Replace OTU data with batch-corrected data
bact_data_bc <- bact_data_filtered
bact_data_bc@otu_table <- otu_table(fit_adjust_batch$feature_abd_adj, taxa_are_rows = TRUE)
# Estimate Shannon index and retrieve the number of ASVs (observed) for bacteria
sample_data(bact_data_bc)$Shannon <- estimate_richness(bact_data_bc, split = TRUE, measures = c('Shannon'))$Shannon
# Estimate Chao1 index (richness)
sample_data(bact_data_bc)$Chao1 <- estimate_richness(bact_data_bc, split = TRUE, measures = c('Chao1'))$Chao1
# Estimate Simpson index
sample_data(bact_data_bc)$Simpson <- estimate_richness(bact_data_bc, split = TRUE, measures = c('Simpson'))$Simpson
# Estimate Hill diversity with a phylogenetic tree
library(hillR)
sample_data(bact_data_bc)$Hill_Shannon_phylo <- hill_phylo(comm = data.frame(otu_table(bact_data_bc)) %>% rotate_df(), tree = bact_data_bc@phy_tree, q = 1)
# Or Hill diversity without a phylogenetic tree
sample_data(bact_data_bc)$Hill_Shannon_taxa <- hill_taxa(comm = data.frame(otu_table(bact_data_bc)) %>% rotate_df(), q = 1)
# Save phyloseq
saveRDS(bact_data_bc, here('output', 'bact_data_bc.rds'))We should also use PERMANOVA analysis to determine whether batch correction was successful. We will use the adonis2 function from the vegan package.
# Obtain the counts data and rotate it for downstream SAMPLE-wise PERMANOVA
bact_counts <- t(data.frame(otu_table(bact_data_filtered)))
bact_counts_bc <- t(data.frame(otu_table(bact_data_bc)))
# Replace the zero values with a very small number to avoid NaN values in the dissimilarity matrix
bact_counts[bact_counts == 0] <- 1e-12
bact_counts_bc[bact_counts_bc == 0] <- 1e-12
# Generate the Bray-Curtis dissimilarity matrix (or the appropriate distance metric)
bact_bray_dist <- vegdist(bact_counts, method = 'bray')
bact_bray_dist_bc <- vegdist(bact_counts_bc, method = 'bray')
# Run PERMANOVA
bact_bray_permanova <- adonis2(bact_bray_dist ~ batch + group, data = data.frame(sample_data(bact_data_filtered)), permutations = 1e6)
perm_label <- paste0('SS: ', round(bact_bray_permanova$SumOfSqs[1], 4),
'\nR2: ', round(bact_bray_permanova$R2[1], 4),
'\np: ', round(bact_bray_permanova$`Pr(>F)`[1], 4))
bact_bc_bray_permanova <- adonis2(bact_bray_dist_bc ~ batch + group, data = data.frame(sample_data(bact_data_bc)), permutations = 1e6)
perm_bc_label <- paste0('SS: ', round(bact_bc_bray_permanova$SumOfSqs[1], 4),
'\nR2: ', round(bact_bc_bray_permanova$R2[1], 4),
'\np: ', round(bact_bc_bray_permanova$`Pr(>F)`[1], 4))
# Bacterial ordination
(bact_ord_plot <- plot_ordination(bact_data_filtered, ordinate(bact_data_filtered, 'PCoA', 'bray', weighted = TRUE)) +
stat_ellipse(aes(fill = batch), geom = 'polygon', type = 't', level = 0.9, alpha = 0.2) +
scale_shape_identity() +
geom_point(aes(fill = batch), shape = 21, size = 3) +
scale_fill_jama(name ='Batch') +
annotate('text', x = 0.4, y = 0.3, label = perm_label) +
theme_pubr(legend = 'right') +
labs(title = 'Bacteria PCoA - Bray-Curtis distance',
subtitle = 'Pre batch-correction with MMUPHin'))
(bact_bc_ord_plot <- plot_ordination(bact_data_bc, ordinate(bact_data_bc, 'PCoA', 'bray', weighted = TRUE)) +
stat_ellipse(aes(fill = batch), geom = 'polygon', type = 't', level = 0.9, alpha = 0.2) +
scale_shape_identity() +
geom_point(aes(fill = batch), shape = 21, size = 3) +
scale_fill_jama(name ='Batch') +
annotate('text', x = 0.3, y = -0.4, label = perm_bc_label) +
theme_pubr(legend = 'right') +
labs(title = 'Bacteria PCoA - Bray-Curtis distance',
subtitle = 'Post batch-correction with MMUPHin'))
p_arr <- ggarrange(bact_ord_plot, bact_bc_ord_plot, nrow = 1, common.legend = TRUE, legend = 'right')
ggsave(here('figures', 'ordination', 'bact_PCoA_unifrac_batchCorrected.pdf'), p_arr,
width = 28, height = 14, units = 'cm')# Generate the unifrac dissimilarity matrix (or the appropriate distance matrix)
bact_unifrac_dist <- phyloseq::distance(bact_data_filtered, method = 'wunifrac', type = 'samples')
bact_bc_unifrac_dist <- phyloseq::distance(bact_data_bc, method = 'wunifrac', type = 'samples')
# Run PERMANOVA
bact_unifrac_permanova <- adonis2(bact_unifrac_dist ~ location, data = data.frame(sample_data(bact_data_filtered)), permutations = 1e6)
perm_label <- paste0('SS: ', round(bact_unifrac_permanova$SumOfSqs[1], 4),
'\nR2: ', round(bact_unifrac_permanova$R2[1], 4),
'\np: ', round(bact_unifrac_permanova$`Pr(>F)`[1], 4))
bact_bc_unifrac_permanova <- adonis2(bact_bc_unifrac_dist ~ location, data = data.frame(sample_data(bact_data_bc)), permutations = 1e6)
perm_bc_label <- paste0('SS: ', round(bact_bc_unifrac_permanova$SumOfSqs[1], 4),
'\nR2: ', round(bact_bc_unifrac_permanova$R2[1], 4),
'\np: ', round(bact_bc_unifrac_permanova$`Pr(>F)`[1], 4))
# Bacterial ordination
(bact_ord_plot <- plot_ordination(bact_data_filtered, ordinate(bact_data_filtered, 'PCoA', 'unifrac', weighted = TRUE)) +
stat_ellipse(aes(fill = batch), geom = 'polygon', type = 't', level = 0.9, alpha = 0.2) +
scale_shape_identity() +
geom_point(aes(fill = batch), shape = 21, size = 3) +
scale_fill_jama(name ='Batch') +
annotate('text', x = 0.4, y = 0.3, label = perm_label) +
theme_pubr(legend = 'right') +
labs(title = 'Bacteria PCoA - Weighted Unifrac',
subtitle = 'Pre batch-correction with MMUPHin'))
(bact_bc_ord_plot <- plot_ordination(bact_data_bc, ordinate(bact_data_bc, 'PCoA', 'unifrac', weighted = TRUE)) +
stat_ellipse(aes(fill = batch), geom = 'polygon', type = 't', level = 0.9, alpha = 0.2) +
scale_shape_identity() +
geom_point(aes(fill = batch), shape = 21, size = 3) +
scale_fill_jama(name ='Batch') +
annotate('text', x = 0.3, y = -0.4, label = perm_bc_label) +
theme_pubr(legend = 'right') +
labs(title = 'Bacteria PCoA - Weighted Unifrac',
subtitle = 'Post batch-correction with MMUPHin'))
p_arr <- ggarrange(bact_ord_plot, bact_bc_ord_plot, nrow = 1, common.legend = TRUE, legend = 'right')
ggsave(here('figures', 'ordination', 'bact_PCoA_unifrac_batchCorrected.pdf'), p_arr,
width = 28, height = 14, units = 'cm')An example of the batch correction with 16S sequencing from two different sequencing locations is shown below, with weighted UniFrac distance.
Prior to conducting multivariate analysis, it is important to consider the structure of the data. First, high-throughput sequencing data has high variability in library size and therefore differences between samples do not demonstrate true biological diversity. Second, shotgun counts are sparse, with most components containing a zero count. To accommodate for the high variability and sparseness, we require normalisation techniques to improve downstream statistical analysis.
We will use 3 different normalisation methods:
Centered log ration transformation (CLR): This normalisation is robust to compositionality. This is the preferred transformation when using some multivariate approaches such as sPLS-DA.
Cumulative Sum Scaling transformation (CSS): This normalisation is robust to compositionality and has been specifically developed for microbiome data.
Log Cumulative Sum Scaling transformation (logCSS): The log + 1 transformation of CSS-normalised data.
# Perform centred log ratio transformation
bact_data_CLR <- microbiome::transform(bact_data_filtered, transform = 'clr')
saveRDS(bact_data_CLR, here('output', 'bact_data_CLR.rds'))
# Perform cumulative sum scaling transformation
bact_data_CSS <- bact_data_filtered
otu_table(bact_data_CSS) <-
otu_table(MRcounts(cumNorm(phyloseq_to_metagenomeSeq(bact_data_filtered), p = 0.5)), taxa_are_rows = TRUE)
saveRDS(bact_data_CSS, here('output', 'bact_data_CSS.rds'))
# Perform cumulative sum scaling tranformation followed by log transformation
bact_data_logCSS <- microbiome::transform(bact_data_CSS, transform = 'log')
saveRDS(bact_data_logCSS, here('output', 'bact_data_logCSS.rds'))At this point, we can also write each of the OTU tables to .csv files using the otu_to_csv() function provided here.
# Write OTU tables to .csv file
otu_to_csv(bact_data_filtered, here::here('output', 'otu_table_bact_data_filtered.csv'))
otu_to_csv(bact_data_CLR, here::here('output', 'otu_table_bact_data_CLR.csv'))
otu_to_csv(bact_data_CSS, here::here('output', 'otu_table_bact_data_CSS.csv'))
otu_to_csv(bact_data_logCSS, here::here('output', 'otu_table_bact_data_logCSS.csv'))We may also want to generate agglomerated forms of the phyloseq at different taxonomic levels. This will allow us to interrogate differences at the genus or family levels, for example.
Double check whether your taxonomy table uses the capitalised form of the taxonomic level, i.e. 'genus' or 'Genus'.
# Prepare genus-agglomerated phyloseq objects
bact_genus_filtered <- tax_glom(bact_data_filtered, taxrank = 'Genus')
bact_genus_logCSS <- tax_glom(bact_data_logCSS, taxrank = 'Genus')
# Save objects to RDS files
saveRDS(bact_genus_filtered, here::here('output', 'bact_genus_filtered.rds'))
saveRDS(bact_genus_logCSS, here::here('output', 'bact_genus_logCSS.rds'))
# Write OTU tables to .csv file
otu_to_csv(bact_genus_filtered, here::here('output', 'otu_tables', 'otu_table_bact_genus_filtered.csv'))
otu_to_csv(bact_genus_logCSS, here::here('output', 'otu_tables', 'otu_table_bact_genus_logCSS.csv'))If you used this repository in a publication, please mention its url.
- Copyright (c) 2023 Mucosal Immunology Lab, Monash University, Melbourne, Australia.
- Authors: M. Macowan