Evaluating sequencing-depth uncertainty with Sibyl

Overview

Sibyl is an R package for evaluating for evaluating how sequencing-depth, and thus uncertainty on community composition, affects beta-diversity ordinations. It does this through a repeated-ordination framework in which alternative representations of the same microbial count data are repeatedly generated, ordinated, aligned, and summarized. The package currently supports two uncertainty-generation workflows: repeated rarefaction and Dirichlet Monte Carlo sampling through ALDEx2.

In the rarefaction workflow, these plausible representations are generated by repeated subsampling at user-defined sequencing-depth thresholds. In the probability-based workflow, they are generated by Dirichlet Monte Carlo sampling through ALDEx2 and analysed using Aitchison distances.

The first part of this vignette focuses on the rarefaction workflow, because it also supports threshold selection. The final section demonstrates how the same downstream framework can be applied to Dirichlet Monte Carlo compositions.

Regardless of which uncertainty generator is used, Sibyl follows the same downstream framework:

Generate alternative community representations → compute ordinations → align replicate configurations by Procrustes transformation → summarize sample positions using consensus coordinates → visualize uncertainty with ordination clouds.

When performing Principal Coordinates Analysis (PCoA), or other types of ordinations, it is necessary to choose a rarefaction threshold which does not affect the structure of the data during exploratory analysis. This is one application the framework supports, not its only purpose.

The two main user-controlled ideas in Sibyl are the uncertainty generator and the number of replicate representations. The uncertainty generator defines how alternative versions of the community table are produced. In the rarefaction workflow, these are generated by repeated subsampling at one or more sequencing-depth thresholds. In the probability-based workflow, they are generated by Dirichlet Monte Carlo sampling. The number of replicates or draws controls how many alternative representations are generated, and therefore how finely uncertainty in ordination space is estimated.

In the rarefaction workflow, users additionally define the sequencing-depth threshold, or range of thresholds, to evaluate.

Keep in mind that increasing the number of replicates requires more computational power. The main limiting factor is the memory needed to store a large number of ordination results.

Example data

The data contained in the package and used for illustration comes from Illumina MiSeq paired-end sequencing of regions V3-V4 of the 16S genes. The samples used in the study are adult mosquitoes from Burkina Faso, a more thorough description is in the published paper from Buck et al. (2016). The other data provided in the package is a collection of mosquito larvae from Ethiopia from Assentato et al. (2025), which provides a wider array of sequence depths, from low to high.

Rarefaction workflow

Threshold selection is one application of the repeated-ordination framework: by repeatedly subsampling a dataset at a given sequencing depth, Sibyl estimates how much a sample’s ordination position is affected by that choice of depth.

Sample completeness and candidate threshold range

After loading the package the first step is to perform the calculation of accumulation (rarefaction) curves. This allows a rough understanding of how the sequencing depth relates to the completeness of each sample.

The function performs the calculation of rarefaction curves with the selected step size (a coarser step size requires less computational power).

The Abundance Coverage Estimator value (ACE) is calculated for each sample as a measure of sample richness, and a generic accumulation curve is fitted to the data, with its asymptote being the ACE value. This functionality is particularly useful to give a visual impression of where the curve “break-point” is, or where the curve is approaching the plateau, fixed as 75% of the ACE value.

Since in reality these samples never reach their ACE value, which is always way over their actual sequencing depth (due to it being an estimation), the “plateau” 75% threshold is to be considered just a guidance estimate. It is not an empiric measure, therefore should not be treated as such.

accumulation_result <- accumulation_test(adults, step = 5)

The resulting list contains:

accumulation_plot: Aggregated plot for all samples with fitted accumulation curves.
threshold_density: 75% thresholds density plot.
individual_plots: Individual rarefaction plots with fitted accumulation curves.

Let’s take a look at the first 4 individual accumulation plots:

ggpubr_installed <- requireNamespace("ggpubr", quietly = TRUE)
if (!ggpubr_installed) {
  message("Install 'ggpubr' to reproduce the vignette examples: install.packages('ggpubr')")
}

ggpubr::ggarrange(plotlist = accumulation_result$individual_plots[1:4],
                  nrow = 2,
                  ncol = 2)

The density plot is a visual summary of the range at which most of the 75% thresholds calculated on the dataset fell.

accumulation_result$threshold_density

Together, the results produced from this test can be used to choose an appropriate threshold range to perform the rest of the analysis.

Ideally, while prioritizing sample completeness, we want to make sure to select a threshold which allows us to keep most of our samples. In this particular case, to not render the dataset unbalanced.

Repeated threshold testing with `test_threshold()`

We proceed by running now the main testing function in the Sibyl package.

Since most of our 75% thresholds, according to the accumulation analysis, falls in the range between 1000-2000 (and since we are interested in seeing how low we can push this specific dataset), we select a range going from 100 to 2000.

Optionally, the repeat number can be a vector of different values. In that case the testing over different thresholds takes place at different repeat amounts, and it is possible to directly explore how the amount of subsampling events affects rarefaction and sample position in the ordination.

The choice of a range to test your data on should always be made considering the compromise between sample completeness (how well described is the sample microbial community according to the previous analysis), and how many samples are excluded or included when changing rarefaction threshold.

threshold_test_result <- test_threshold(adults,
                                        repeats = 100,
                                        t_min = 100,
                                        t_max = 2000,
                                        t_step = 20,
                                        group = "location",
                                        cores = 4)

This function produces a series of results which can be visualized immediately or through auxiliary functions.

The Calinski-Harabasz (CH) index

That is probably the most important result. It shows how the Calinski-Harabasz index (CH) changes when increasing rarefaction threshold. The index is a ratio of the between-cluster separation (BCSS) to the within-cluster dispersion (WCSS), normalized by their number of degrees of freedom. A higher index means that the between cluster separation is getting bigger and clusters are less spread out. The index is calculated based on the consensus coordinates obtained from the repeated subsampling. The clusters are defined as the “group” variable used in the function “threshold_test_result”.

We are looking for a point in which the CH index is close to a plateau, which indicates that a higher sequencing effort (a higher rarefaction threshold) does not affect the ordination.

threshold_test_result$index_plot

In this other example, we can see what happens when there are samples in the dataset which do not have enough reads to meet the rarefaction threshold.

We get a sudden drop in the value of the CH index which indicates that samples have been removed, and there has been either a decrease of BCSS or an increase of WCSS (or both). However, this is a good indication of when samples are getting removed, and if it is important to keep all samples in the analysis, we might want to operate a choice that maximize the CH index.

This has been performed on a different dataset, accessible (when package is not loaded) with:

data(larvae)

Ordination clouds across thresholds

The ordination_plots element contains all the individual ordinations for the selected repetition number across the whole threshold range to test.

Let’s take a look at some of them (lowest threshold, median and highest):

# Isolate the repeat number we need: 
ordination_plots <- threshold_test_result$ordination_plots$`repeat_number 100`

It is possible to see how the ordination changes according to the rarefaction threshold, and how the highest value is not necessarily making a big difference when compared to the median. The biggest impact lies in the initial increase of the threshold, but after a critical point the sample clouds are stable and thus not really affected by the stochastic subsampling procedure.

Thus, this also means that going under a specific rarefaction value might impact our ordination because of an artificial increase of point dispersal, which in turn introduces a bias creating an artificial overlap of two groups which are normally separated.

# Print first, middle and last plot:
n <- length(ordination_plots)
# Get indices
indices <- unique(c(1, ceiling(n / 2), n))
selected_plots <- lapply(indices, function(i) ordination_plots[[i]])

ggpubr_installed <- requireNamespace("ggpubr", quietly = TRUE)
if (!ggpubr_installed) {
  message("Install 'ggpubr' to reproduce the vignette examples: install.packages('ggpubr')")
}
ggpubr::ggarrange(plotlist = selected_plots,
                  ncol = 3,
                  common.legend = TRUE)

Average pairwise distance (APD)

The avg_distances element is a list of dataframes which contains the values of average pairwise distances (APD) for each cloud of points (subsampling ordination results) generated from a single sample.

A decreasing trend means that the average pairwise distance of the cloud gets smaller.

Measuring APD is especially valuable to give further insight into when we reach a point of diminishing returns on the ordination structure. When increasing the threshold does not increase APD, it means that a sample position is less affected by repeating the subsampling.

APD_plots <- threshold_test_result$avg_distances$repeat_number_100
avg_pairwise_dist_plot(APD_plots)

Interpreting the rarefaction workflow

Overall in this worked example, data suggests that a threshold as low as 750 does not substantially affect the ordination structure (according to CH index), and does not affect the single sample position when repeating subsampling (looking at the APD plots). Despite that, when considering the accumulation curve we can also consider that 750 reads is most likely not enough to describe accurately our samples.

The rarefaction workflow allowed us to explore the effect of different rarefaction thresholds, and to confirm that lowering a threshold to include more samples does not necessarily affect our exploratory analysis, which is still valid for relatively low threshold numbers. This is one application of Sibyl’s repeated-ordination framework, focused specifically on threshold selection. The same framework can be used more generally to evaluate ordination uncertainty, as shown next.

Probability-based workflow

The Dirichlet workflow uses the same downstream repeated-ordination framework as the rarefaction workflow, but allows users to examine compositional uncertainty without discarding reads through rarefaction.

Instead of repeatedly subsampling reads, dirichlet_ordination() performs repeated Monte Carlo sampling from the Dirichlet posterior distribution estimated by ALDEx2::aldex.clr(). Each Monte Carlo draw represents one plausible compositional state of the observed count data, given its Dirichlet posterior. These replicate compositions are ordinated using Aitchison distances, aligned, and summarized with consensus coordinates and uncertainty clouds as per the rarefaction workflow.

Ordination clouds should be interpreted as probability-based compositional uncertainty, arising from how confidently the observed counts reflect the underlying community proportions.

rarefaction_result <- repeated_rarefaction(adults,
                                           repeats = 50,
                                           threshold = 150,
                                           group = "location",
                                           colorb = "location",
                                           cloud = TRUE,
                                           ellipse = TRUE)


dirichlet_result <- dirichlet_ordination(adults,
                                         draws = 50,
                                         group = "location",
                                         colorb = "location",
                                         cloud = TRUE,
                                         ellipse = TRUE)

ggpubr_installed <- requireNamespace("ggpubr", quietly = TRUE)
if (!ggpubr_installed) {
  message("Install 'ggpubr' to reproduce the vignette examples: install.packages('ggpubr')")
}

ggpubr::ggarrange(rarefaction_result$plot + ggtitle("Rarefaction"),
                  dirichlet_result$plot + ggtitle("Dirichlet"),
                  ncol = 2, common.legend = TRUE, legend = "right")

Both panels show the same samples, with consensus coordinates highlighted in red. Output can be compared conceptually with the rarefaction-based ordination clouds, but since the two workflows use different distance measures (Bray-Curtis for rarefaction, Aitchison for the Dirichlet approach) and different sources of stochasticity, the absolute orientation and scale of the two ordinations should not be overinterpreted.

What is comparable is the relative spread of each sample’s cloud of repeats/draws: samples with lower sequencing depth tend to show a wider cloud under repeated rarefaction, since low depth means each subsample discards proportionally more of the observed community. The Dirichlet-based approach instead captures uncertainty in the estimated community composition given the observed counts, rather than uncertainty introduced by repeatedly discarding reads through rarefaction.

Using both approaches side by side gives a more complete picture of how much of the variability observed in an ordination is attributable to the chosen subsampling scheme, versus how much reflects genuine uncertainty in the compositional estimate itself.

Practical notes

Number of repeats/draws: more replicates give a finer estimate of ordination uncertainty, at the cost of more computation. A small number (below 5) is not enough to draw confidence ellipses, and Sibyl will warn and skip them automatically.
Memory and runtime: each replicate requires its own ordination and Procrustes alignment. The main limiting factor for larger datasets or higher repeat/draw counts is usually available memory to store all the replicate results, not raw computation time.
Sample retention: in the rarefaction workflow, samples whose read count falls below the chosen threshold are removed before ordination, and Sibyl issues a warning listing them. The Dirichlet workflow does not discard samples this way, since it does not depend on a sequencing-depth threshold.
Comparing across distance metrics: rarefaction-based ordinations (Bray-Curtis) and Dirichlet-based ordinations (Aitchison) are not directly comparable in absolute orientation or scale. Compare the relative spread of sample clouds between the two workflows, not their exact coordinates.

Data references

Buck, M., Nilsson, L., Brunius, C. et al. Bacterial associations reveal spatial population dynamics in Anopheles gambiae mosquitoes. Sci Rep 6, 22806 (2016).
Assentato, L., Nilsson, L., Brunius, C. et al. The type of environment has a greater impact on the larval microbiota of Anopheles arabiensis than on the microbiota of their breeding water, FEMS Microbiology Ecology, Volume 101, Issue 1, (2025).

Lorenzo Assentato

2026-07-20

Overview

Example data

Rarefaction workflow

Sample completeness and candidate threshold range

Repeated threshold testing with `test_threshold()`

The Calinski-Harabasz (CH) index

Ordination clouds across thresholds

Average pairwise distance (APD)

Interpreting the rarefaction workflow

Probability-based workflow

Practical notes

Data references

Evaluating sequencing-depth uncertainty with Sibyl

Lorenzo Assentato

2026-07-20

Overview

Example data

Rarefaction workflow

Sample completeness and candidate threshold range

Repeated threshold testing with test_threshold()

The Calinski-Harabasz (CH) index

Ordination clouds across thresholds

Average pairwise distance (APD)

Interpreting the rarefaction workflow

Probability-based workflow

Practical notes

Data references

Repeated threshold testing with `test_threshold()`