Data Exploration and Parameter Selection for Core Identification
Rarefaction metrics and visualization
We are not going to talk about the importance of rarefaction as it is not the goal of this vignette, but if you are interested you should read Schloss (2024a), Schloss (2024b), Hong et al. (2022), and Weiss et al. (2017) as a start.
In BRCore, to identify the ideal rarefaction cutoff we
can plot a simple series of visuals to guide us in the decision. We can,
more or less, estimate how well a DNA sample will sequence, but there
are always some samples that fail to sequence (or sequence badly) and
some others that will sequence “too well” for reasons that are outside
our control or by chance.
In the end we will need to decide how many samples we are willing to discard and how much diversity we want to retain in our data. To do that we can use some help by plotting some diagnostics plots as shown below.
bcse_metrics <- add_rarefaction_metrics(data = bcse)
bcse_metrics
#> phyloseq-class experiment-level object
#> otu_table() OTU Table: [ 2861 taxa and 50 samples ]
#> sample_data() Sample Data: [ 50 samples by 7 sample variables ]
#> tax_table() Taxonomy Table: [ 2861 taxa by 10 taxonomic ranks ]
#> refseq() DNAStringSet: [ 2861 reference sequences ]rarefaction_plot <- plot_rarefaction_metrics(bcse_metrics)
#> ℹ Processing 50 samples
#>
ℹ Generating rarefaction diagnostic plots
✔ Rarefaction diagnostic plots generated successfully
#> ℹ Generating rarefaction diagnostic plots
✔ Generating rarefaction diagnostic plots [1.2s]
print(rarefaction_plot)
Figure 1: Rarefaction metrics. a-b, histograms of samples. c, Good’s coverage per total number of sequences in a sample. d-e, log10 of sequences in a sample. f, samples ranked by sequence reads in a sample.
Note: The Good’s coverage” is a statistics used in ecology that estimates the proportion of species in a community that are represented in a sample, based on the number of species encountered and the total number of individuals sampled. Essentially, it quantifies how well a sample represents the overall diversity of a population or ecosystem. Good’s coverage is calculated as \((1-(n1/N)\), where \(n1\) is the number of unique ASV/OTUs (i.e. species) found only once, and \(N\) is the total number of sequence reads (i.e. individuals); a high percentage (e.g., >95%) means most reads are from common taxa, suggesting sufficient sampling, while low coverage indicates many rare, potentially missed, taxa.
Rarefy the data and explore variance propagation
Originally developed by Sanders (1968) rarefaction has been adopted widely in ecology and microbiome research. It consists of randomly resampling the same number of species (ASV/OTUs) multiple times across all samples to estimate the expected number of species that would be observed if all samples were identical. It prevents samples with more individuals (e.g., more sequencing reads) from appearing more diverse simply because they were sampled more in depth.
Because rarefaction is a random process, it adds noise to the data. Each resampling event is slightly different from others due to chance. This can be problematic, especially when a very low rarefaction cutoff is selected, when, for example, we want to compare 2 treatments. With lower rarefaction cutoffs, the greater the chance you have that rarefaction iterations will differ from one another, as also shown by Cameron et al. (2021).
We developed functions that can perform and store all the rarefaction iterations to evaluate (and visualize) the noise generated in the process, as shown below:
bcse_rarefied_list <-
multi_rarefy(
physeq_obj = bcse,
depth_level = 1000,
num_iter = 3,
.as = "list",
set_seed = 7642
)
#>
#> ── Rarefaction iterations starting... ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#>
#> ── Input Validation ──
#>
#> ✔ Input phyloseq object is valid!
#> ℹ Seed: 7642
#> ℹ Input (matrix/df dim): 47 samples x 2861 taxa
#> ℹ Rarefaction depth: 1000
#> ℹ Iterations: 3
#> ℹ taxa_are_rows: TRUE
#> ℹ OTU matrix/df rownames head: bcse50, bcse69, bcse73, bcse191, bcse82, bcse102
#> ℹ OTU matrix/df colnames head: OTU_427, OTU_11, OTU_253, OTU_148, OTU_3, OTU_78
#> ℹ Row sums summary: Min=1193, Max=107643, Median=25209
#>
#> ── Rarefaction Results ──
#>
#> ── Sample Removal
#> ! 3 samples removed (depth < 1000)
#> ! Samples removed: "bcse108, bcse105, bcse110"
#>
#> ── Taxa Removal
#> ✔ No taxa removed.
#> ! Taxa are not removed across iterations to maintain consistent dimensions.
#> Downstream analyses should handle zero-abundance taxa appropriately.
#>
#> ── Data Sparsity
#> ℹ Returning list of data frames for each iteration.
#> • Rarefied matrix (across 3 iterations):
#> • Min: 130570 zeros (97.1% sparsity) out of 134467 entries
#> • Max: 130663 zeros (97.17% sparsity) out of 134467 entries
#> • Avg: 130608 zeros (97.13% sparsity) out of 134467 entries
#>
#> ── Final Data Dimensions
#> ✔ Output: 3 iterations with 50 unique samples
#> • Samples per iteration:
#> • Min: 47
#> • Max: 47
#> • Non-zero taxa per iteration:
#> • Min: 1046
#> • Max: 1098
#> • Avg: 1073Let’s explore the results.
class(bcse_rarefied_list)
#> [1] "list"
str(bcse_rarefied_list[[1]], list.len = 10) # Dimensions of iteration #1
#> 'data.frame': 47 obs. of 2861 variables:
#> $ OTU_427 : int 0 0 0 0 0 0 0 0 0 0 ...
#> $ OTU_11 : int 0 0 0 0 0 0 0 0 0 0 ...
#> $ OTU_253 : int 0 0 0 0 0 0 0 0 0 0 ...
#> $ OTU_148 : int 0 1 0 11 2 0 0 0 0 78 ...
#> $ OTU_3 : int 261 86 651 272 249 114 192 73 96 320 ...
#> $ OTU_78 : int 0 0 0 0 0 0 0 0 0 0 ...
#> $ OTU_58 : int 0 0 0 0 0 0 0 0 0 0 ...
#> $ OTU_152 : int 0 0 0 0 0 0 0 0 0 0 ...
#> $ OTU_16 : int 2 0 0 1 1 0 0 1 0 0 ...
#> $ OTU_35 : int 0 0 0 0 0 0 0 0 0 0 ...
#> [list output truncated]We can plot the variance (the noise we talked above) across the 3
iterations of rarefaction. We can use the
plot_variance_propagation() function to visualize this
variance in Alpha diversity across samples and iterations and compare
Raw (non-rarefied) vs. Rarefied data.
rarefaction_variance_plot <- plot_variance_propagation(
physeq_obj = bcse,
rarefied = bcse_rarefied_list,
q = 0,
group_var = "Crop",
group_color = "Plot"
)
#> ✔ Input phyloseq object is valid!
#>
#> ── Rarefaction Variance Propagation Visualization ────────────────────────────────────────────────────────────────────────────────────────────────────────
#> ℹ Hill number order selected, q= 0
#> ℹ Number of rarefaction iterations, n_iter= 3
#> ℹ Comparison plot generated!
print(rarefaction_variance_plot)
Figure 2: Alpha diversity variance propagation (variance between rarefaction iteration) across 3 iterations of rarefaction. Points represent samples. The x-axis represent the grouping variable (Crop). The y-axis represent the alpha diversity metric (q = 0, i.e. richness) calculated on the samples. The color of the points represent the ‘Plot’ variable, which is a nested variable within ‘Crop’.
Review the rarefaction and variance plots to determine the
ideal sampling depth. You are looking for a “sweet spot” that
maintains high diversity and samples while keeping variance across
iterations low. Pay extra attention when rarefaction cutoff is quite
low, as show by Cameron et al. (2021). Note this
depth, as it will serve as the cutoff for the
identify_core() step in your pipeline.
Recreate a phyloseq object with a rarefied
otu_table()
This is an optional step.
The identify_core() function accepts a phyloseq object,
whether rarefied or non-rarefied. When supplying a phyloseq object with
a non-rarefied otu_table, you must provide a list of
rarefied OTU tables that result from multi_rarefy(). In our
case, we generated bcse_rarefied_list.
The results from add_rarefaction_metrics(),
plot_rarefaction_metrics(), and multi_rarefy()
are meant to explore the data and inform the chosen parameters for
identify_core(). Function identify_core() will
calculate the distance matrix for all iterations supplied and rank OTUs
based on the average distance across iterations.
If you have a phyloseq object with a rarefied
otu_table() or if you want to use a specific rarefaction
iteration generated by multi_rarefy(), you also have these
options. The choice is up to you and depends on your preferences and
the specific goals of your analysis. The
update_otu_table() function allows for replacing the
otu_table() with the rarefied table you prefer.
bcse_rare_single <- update_otu_table(
physeq_obj = bcse,
rarefied_otus = bcse_rarefied_list,
iteration = 1 # Speficify which iteration to use for the updated OTU table
)
#> ℹ Extracting iteration 1 from list of 3 iterations.
#> ! Phyloseq object and rarefied otu_table sample names are NOT identical. Check samples removed by rarefaction below.
#> ! 3 samples removed due to rarefaction: "bcse108, bcse105, bcse110"
#> ℹ Building phyloseq object with 47 samples and 2861 taxa
#> ✔ Update complete!Let’s proceed to the core analysis.
