From 74129d5063653f8b2b36f5f82bb1a494e2656c9b Mon Sep 17 00:00:00 2001 From: Lokesh9106 <2400040120@kluniversity.in> Date: Sun, 23 Nov 2025 20:04:34 +0530 Subject: [PATCH 1/3] docs: Add comprehensive Rarefaction section to transformation chapter Added new section 12.3 Rarefaction to address issue #823. Changes include: - Introduction to rarefaction with rarefyAssay() and niter parameter - Subsection on using rarefaction with alpha diversity (addAlpha) - Subsection on using rarefaction with beta diversity (addMDS) - Function comparison explaining differences between: * addAlpha() vs getAlpha() * runMDS() vs addMDS() Includes practical code examples demonstrating iterative rarefaction with niter=100. --- inst/pages/transformation.qmd | 63 +++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/inst/pages/transformation.qmd b/inst/pages/transformation.qmd index c39fc546..b8330e6b 100644 --- a/inst/pages/transformation.qmd +++ b/inst/pages/transformation.qmd @@ -141,6 +141,69 @@ than the minimum abundance value before transformation. Some tools, like values. See [@sec-differential-abundance]. ::: +## Rarefaction {#sec-rarefaction} + +Another approach to control uneven sampling depths is to apply rarefaction with `rarefyAssay()`, which resamples the samples to an equal number of reads. This remains controversial, however, and strategies to mitigate the information loss in rarefaction have been proposed [@Schloss_2024a; @Schloss_2024b]. Moreover, this practice has been discouraged for the analysis of differentially abundant microorganisms [@McMurdie_and_Holmes_2014]. + +Rarefaction can be performed iteratively by using the `niter` parameter in `rarefyAssay()`. This creates multiple rarefied versions of the data, which can help account for the stochasticity introduced by random subsampling. The resulting rarefied assays can then be used for downstream analyses such as alpha and beta diversity calculations. + +### Using rarefaction with alpha diversity + +When calculating alpha diversity indices, you can apply rarefaction iteratively and then compute diversity metrics across the rarefied replicates. The `addAlpha()` function can work with rarefied data: + +```{r} +#| label: rarefaction-alpha +#| eval: false + +# Load example data +library(mia) +data("Tengeler2020") +tse <- Tengeler2020 + +# Get minimum read depth for rarefaction +min_reads <- min(colSums(assay(tse, "counts"))) + +# Perform iterative rarefaction +tse <- rarefyAssay( + tse, + method = "subsample", + sample = min_reads, + niter = 100 +) + +# Calculate alpha diversity on rarefied data +tse <- addAlpha( + tse, + assay_name = "counts_rarefied", + sample = min_reads, + niter = 100 +) +``` + +### Using rarefaction with beta diversity + +Similarly, rarefaction can be applied before calculating beta diversity and performing ordination. The `addMDS()` function can utilize rarefied data for more robust distance calculations: + +```{r} +#| label: rarefaction-beta +#| eval: false + +# Perform MDS ordination on rarefied data +tse <- addMDS( + tse, + assay_name = "counts_rarefied", + method = "bray", + niter = 100 +) +``` + +### Function comparison + +**`addAlpha()` vs `getAlpha()`**: Both functions calculate alpha diversity indices, but `addAlpha()` stores the results directly into the `colData` of the TreeSummarizedExperiment object, while `getAlpha()` returns the diversity values as a separate vector or matrix. Use `addAlpha()` when you want to keep all data together in one object, and `getAlpha()` when you need the diversity values for immediate use in other calculations. + +**`runMDS()` vs `addMDS()`**: The `runMDS()` function calculates multidimensional scaling coordinates and returns them as a separate matrix, whereas `addMDS()` calculates the MDS coordinates and stores them directly into the `reducedDim` slot of the TreeSummarizedExperiment object. Using `addMDS()` is generally preferred as it maintains all results within the same data object, making downstream analyses and visualization more straightforward. + + ## Transformations in practice Below, we apply relative transformation to counts table. From f067bbb742f1a0a34cf5210fc36c73e85be5e756 Mon Sep 17 00:00:00 2001 From: Daena Rys Date: Mon, 4 May 2026 11:40:22 +0300 Subject: [PATCH 2/3] up --- inst/pages/transformation.qmd | 70 ++++++++++++++++++++++++++++++----- 1 file changed, 60 insertions(+), 10 deletions(-) diff --git a/inst/pages/transformation.qmd b/inst/pages/transformation.qmd index b1733274..d89cc102 100644 --- a/inst/pages/transformation.qmd +++ b/inst/pages/transformation.qmd @@ -145,15 +145,15 @@ values. See [@sec-differential-abundance]. Another approach to control uneven sampling depths is to apply rarefaction with `rarefyAssay()`, which resamples the samples to an equal number of reads. This remains controversial, however, and strategies to mitigate the information loss in rarefaction have been proposed [@Schloss_2024a; @Schloss_2024b]. Moreover, this practice has been discouraged for the analysis of differentially abundant microorganisms [@McMurdie_and_Holmes_2014]. -Rarefaction can be performed iteratively by using the `niter` parameter in `rarefyAssay()`. This creates multiple rarefied versions of the data, which can help account for the stochasticity introduced by random subsampling. The resulting rarefied assays can then be used for downstream analyses such as alpha and beta diversity calculations. +Rarefaction can be performed iteratively by using the `niter` parameter in `rarefyAssay()`. This creates multiple rarefied versions of the data, which can help account for the stochasticity introduced by random subsampling. The resulting rarefied assays can then be used for downstream analyses such as alpha and beta diversity calculations. For alpha and beta diversity, the same repeated subsampling can also be done directly within `addAlpha()` and `addMDS()` by setting `niter`. ### Using rarefaction with alpha diversity -When calculating alpha diversity indices, you can apply rarefaction iteratively and then compute diversity metrics across the rarefied replicates. The `addAlpha()` function can work with rarefied data: +When calculating alpha diversity indices, you can first create rarefied assays with `rarefyAssay()` and reuse them later. This is useful when you want to inspect or store the rarefied data itself: ```{r} #| label: rarefaction-alpha -#| eval: false +#| eval: true # Load example data library(mia) @@ -170,39 +170,89 @@ tse <- rarefyAssay( sample = min_reads, niter = 100 ) +``` + +If you only need alpha diversity, `addAlpha()` can do the iterative subsampling directly from the count assay. It stores the values in the object, while `getAlpha()` returns them separately: + +```{r} +#| label: rarefaction-alpha-direct +#| eval: true + +# Reload example data so this chunk is self-contained. +library(mia) +data("Tengeler2020") +tse <- Tengeler2020 -# Calculate alpha diversity on rarefied data +# Calculate alpha diversity with iterative rarefaction tse <- addAlpha( tse, - assay_name = "counts_rarefied", - sample = min_reads, + assay.type = "counts", + index = "shannon", niter = 100 ) ``` ### Using rarefaction with beta diversity -Similarly, rarefaction can be applied before calculating beta diversity and performing ordination. The `addMDS()` function can utilize rarefied data for more robust distance calculations: +Similarly, rarefaction can be applied before calculating beta diversity and performing ordination. If you want to store a rarefied assay first, `rarefyAssay()` still works as above. For ordination, `addMDS()` can also perform the repeated subsampling directly from the count assay: ```{r} #| label: rarefaction-beta -#| eval: false +#| eval: true -# Perform MDS ordination on rarefied data +# Reload example data so this chunk is self-contained. +library(mia) +data("Tengeler2020") +tse <- Tengeler2020 + +# Perform MDS ordination with iterative rarefaction tse <- addMDS( tse, - assay_name = "counts_rarefied", + assay.type = "counts", method = "bray", niter = 100 ) ``` +When you combine iterative rarefaction with a transformation, keep the input as counts and pass the transformation through `transf` instead of pre-transforming the assay. The same pattern can be used for different transformations by swapping the helper function, for example a relative-abundance or clr-style transformation. + +```{r} +#| label: rarefaction-beta-transf +#| eval: true + +# Reload example data so this chunk is self-contained. +library(mia) +library(vegan) +data("Tengeler2020") +tse <- Tengeler2020 + +# Define a custom transformation function. +clr <- function(x) { + vegan::decostand(x, method = "clr", pseudocount = 1) +} + +# Apply the transformation after rarefaction and before the beta diversity calculation. +tse <- addMDS( + tse, + assay.type = "counts", + FUN = getDissimilarity, + method = "euclidean", + niter = 100, + sample = min(colSums(assay(tse, "counts"))), + transf = clr, + replace = TRUE, + name = "MDS_clr_rarefied" +) +``` + ### Function comparison **`addAlpha()` vs `getAlpha()`**: Both functions calculate alpha diversity indices, but `addAlpha()` stores the results directly into the `colData` of the TreeSummarizedExperiment object, while `getAlpha()` returns the diversity values as a separate vector or matrix. Use `addAlpha()` when you want to keep all data together in one object, and `getAlpha()` when you need the diversity values for immediate use in other calculations. **`runMDS()` vs `addMDS()`**: The `runMDS()` function calculates multidimensional scaling coordinates and returns them as a separate matrix, whereas `addMDS()` calculates the MDS coordinates and stores them directly into the `reducedDim` slot of the TreeSummarizedExperiment object. Using `addMDS()` is generally preferred as it maintains all results within the same data object, making downstream analyses and visualization more straightforward. +The same add/get pattern is also used for other ordination workflows where available, for example `runNMDS()` / `addNMDS()` in the community similarity chapter. + ## Transformations in practice From 4c7a37a88c42dcb9ade334e0b6713367ede90a9e Mon Sep 17 00:00:00 2001 From: Tuomas Borman Date: Fri, 18 Sep 2026 13:38:55 +0300 Subject: [PATCH 3/3] Updated to rarefaction --- inst/pages/alpha_diversity.qmd | 44 ++++-- inst/pages/community_similarity.qmd | 37 +---- inst/pages/quality_control.qmd | 4 +- inst/pages/transformation.qmd | 215 +++++++++++++--------------- 4 files changed, 139 insertions(+), 161 deletions(-) diff --git a/inst/pages/alpha_diversity.qmd b/inst/pages/alpha_diversity.qmd index 4ceb23ea..58473ab3 100644 --- a/inst/pages/alpha_diversity.qmd +++ b/inst/pages/alpha_diversity.qmd @@ -123,16 +123,6 @@ modern 16S data, which commonly features denoising and removal of singletons [@Deng2024]. ::: -## Rarefaction {#sec-rarefaction} - -Uneven sequencing read depth can bias diversity metrics. Rarefaction attempts -to normalize read depth by repeatedly subsampling reads and computing -an averaged metric[@Schloss2024rarefaction1;@Schloss2024rarefaction2]. -Rarefaction for alpha diversity metrics can be applied by specifying the optional -arguments `sample` and `niter` in `addAlpha()`. Rarefaction is discussed in -more detail in [@sec-rarefaction_ordination] where it is applied in -the context of beta diversity. - ## Examples ### Calculate alpha diversity measures {#sec-estimate-diversity} @@ -295,6 +285,40 @@ In summary, soil samples exhibit the highest biodiversity, while human-derived samples show less variation in terms of species composition and evolutionary diversity. +### Rarefaction {#sec-rarefaction-alpha} + +To control uneven sampling depths, we can do rarefaction +(see [@sec-rarefaction]). Rarefaction for alpha diversity metrics can be applied +by specifying the optional arguments `sample` and `niter` in `addAlpha()`. +`niter` controls the number of rarefaction iterations (choose a sufficiently +large value for stable estimates). `sample` sets the rarefaction depth, i.e., +the number of reads sampled per iteration from each sample. + +```{r} +#| label: alpha_rarefaction +tse <- addAlpha( + tse, + assay.type = "counts", + index = "observed_richness", + name = "richness_w_rarefaction", + niter = 10 +) +``` + +Next compare the results with and without rarefaction. + +```{r} +#| label: compare:rarefaction +plotBoxplot(tse, col.var = "observed_richness", fill.by = "SampleType") + + plotBoxplot(tse, col.var = "richness_w_rarefaction", fill.by = "SampleType") + + plot_layout(guides = "collect") +``` + +These results illustrate that uneven sequencing depth affects the diversity +estimates: applying rarefaction changes the values compared with the +non-rarefied analysis, indicating that part of the apparent differences was +driven by sampling effort rather than biology. + ### Statistical analysis of alpha diversity measures {#sec-stats-diversity} We can then analyze the statistical significance. We use the non-parametric diff --git a/inst/pages/community_similarity.qmd b/inst/pages/community_similarity.qmd index e958a087..64621039 100644 --- a/inst/pages/community_similarity.qmd +++ b/inst/pages/community_similarity.qmd @@ -426,42 +426,7 @@ plotReducedDim(tse, "unifrac", colour_by = "Group") ### Rarefaction to mitigate impacts of uneven sequencing effort {#sec-rarefaction_ordination} -The sequencing depth of a sample refers to the number of metagenomic reads -obtained from the sequencing process. It is common to find significant variation -in sequencing depth between samples. For instance, the samples of the -*TreeSummarizedExperiment* dataset `GlobalPatterns` show up to a 40-fold -difference in the number of metagenomic reads. Variation in sequencing depth -across samples in a study can bias the calculation of alpha and beta diversity -metrics [@Schloss2023]. - -```{r} -#| label: rarefaction - -# Calculate the list of sequencing depths across samples -sequencing_depths <- colSums(assay(tse)) - -# Calculate variation between highest and lowest sequencing depth -depth_variation <- max(sequencing_depths) / min(sequencing_depths) -depth_variation -``` - -To address uneven sequencing effort, rarefaction aims to normalize metagenomic -reads counts using repeated subsampling. The user first chooses the rarefaction -depth and a number of iterations N. All the samples with metagenomic read counts -below the specified depth are removed and then metagenomic reads are randomly -drawn for the samples left to get subsampled reads at equal read depth. -Then a beta diversity metric is calculated with the subsampled data and the -process is iterated N times. Finally, beta diversity is estimated taking the mean -of all beta diversity values calculated on subsampled data. - -There has been a long-lasting controversy surrounding the use of rarefaction in -microbial ecology. The main concern is that subsampling would omit data -[@mcmurdie2014waste; @Schloss2024rarefaction2]. However, if the subsampling -process is repeated a sufficient number of times, and if the rarefaction depth -is set to the lowest metagenomic reads count found across all samples, no data -will be omitted. Moreover, Patrick Schloss has demonstrated that rarefaction is -"the only method that could control for variation in uneven sequencing effort -when measuring commonly used alpha and beta diversity metrics" [@Schloss2023]. +Background on rarefaction is provided in [@sec-rarefaction]. Let us first convert the count assay to centered log ratio (clr) assay and calculate MDS with Aitchison distance without rarefaction: diff --git a/inst/pages/quality_control.qmd b/inst/pages/quality_control.qmd index 2981b273..0384657a 100644 --- a/inst/pages/quality_control.qmd +++ b/inst/pages/quality_control.qmd @@ -164,7 +164,7 @@ the data as they usually represent sequencing artifacts. See [@sec-subset_prev] and [@sec-agglomerate_prev] for more details on prevalence filtering and agglomeration. -### Library size +### Library size {#sec-library-size} Library size refers to the total number of counts found in a single sample. The returned tables in [@sec-qc-summarize] showed that samples exhibit lots of @@ -205,7 +205,7 @@ sufficient. In case of insufficient sampling depth, one might consider filtering the data based on library size ([@sec-subset-library-size]). To control uneven sampling depths, one should apply data transformation or -apply rarefaction. These both approaches are discussed in [@sec-assay-transform]. +apply rarefaction. These both approaches are discussed in [@sec-rarefaction]. ### Contaminant sequences diff --git a/inst/pages/transformation.qmd b/inst/pages/transformation.qmd index 44032174..c19b9e96 100644 --- a/inst/pages/transformation.qmd +++ b/inst/pages/transformation.qmd @@ -146,119 +146,6 @@ than the minimum abundance value before transformation. Some tools, like values. See [@sec-differential-abundance]. ::: -## Rarefaction {#sec-rarefaction} - -Another approach to control uneven sampling depths is to apply rarefaction with `rarefyAssay()`, which resamples the samples to an equal number of reads. This remains controversial, however, and strategies to mitigate the information loss in rarefaction have been proposed [@Schloss_2024a; @Schloss_2024b]. Moreover, this practice has been discouraged for the analysis of differentially abundant microorganisms [@McMurdie_and_Holmes_2014]. - -Rarefaction can be performed iteratively by using the `niter` parameter in `rarefyAssay()`. This creates multiple rarefied versions of the data, which can help account for the stochasticity introduced by random subsampling. The resulting rarefied assays can then be used for downstream analyses such as alpha and beta diversity calculations. For alpha and beta diversity, the same repeated subsampling can also be done directly within `addAlpha()` and `addMDS()` by setting `niter`. - -### Using rarefaction with alpha diversity - -When calculating alpha diversity indices, you can first create rarefied assays with `rarefyAssay()` and reuse them later. This is useful when you want to inspect or store the rarefied data itself: - -```{r} -#| label: rarefaction-alpha -#| eval: true - -# Load example data -library(mia) -data("Tengeler2020") -tse <- Tengeler2020 - -# Get minimum read depth for rarefaction -min_reads <- min(colSums(assay(tse, "counts"))) - -# Perform iterative rarefaction -tse <- rarefyAssay( - tse, - method = "subsample", - sample = min_reads, - niter = 100 -) -``` - -If you only need alpha diversity, `addAlpha()` can do the iterative subsampling directly from the count assay. It stores the values in the object, while `getAlpha()` returns them separately: - -```{r} -#| label: rarefaction-alpha-direct -#| eval: true - -# Reload example data so this chunk is self-contained. -library(mia) -data("Tengeler2020") -tse <- Tengeler2020 - -# Calculate alpha diversity with iterative rarefaction -tse <- addAlpha( - tse, - assay.type = "counts", - index = "shannon", - niter = 100 -) -``` - -### Using rarefaction with beta diversity - -Similarly, rarefaction can be applied before calculating beta diversity and performing ordination. If you want to store a rarefied assay first, `rarefyAssay()` still works as above. For ordination, `addMDS()` can also perform the repeated subsampling directly from the count assay: - -```{r} -#| label: rarefaction-beta -#| eval: true - -# Reload example data so this chunk is self-contained. -library(mia) -data("Tengeler2020") -tse <- Tengeler2020 - -# Perform MDS ordination with iterative rarefaction -tse <- addMDS( - tse, - assay.type = "counts", - method = "bray", - niter = 100 -) -``` - -When you combine iterative rarefaction with a transformation, keep the input as counts and pass the transformation through `transf` instead of pre-transforming the assay. The same pattern can be used for different transformations by swapping the helper function, for example a relative-abundance or clr-style transformation. - -```{r} -#| label: rarefaction-beta-transf -#| eval: true - -# Reload example data so this chunk is self-contained. -library(mia) -library(vegan) -data("Tengeler2020") -tse <- Tengeler2020 - -# Define a custom transformation function. -clr <- function(x) { - vegan::decostand(x, method = "clr", pseudocount = 1) -} - -# Apply the transformation after rarefaction and before the beta diversity calculation. -tse <- addMDS( - tse, - assay.type = "counts", - FUN = getDissimilarity, - method = "euclidean", - niter = 100, - sample = min(colSums(assay(tse, "counts"))), - transf = clr, - replace = TRUE, - name = "MDS_clr_rarefied" -) -``` - -### Function comparison - -**`addAlpha()` vs `getAlpha()`**: Both functions calculate alpha diversity indices, but `addAlpha()` stores the results directly into the `colData` of the TreeSummarizedExperiment object, while `getAlpha()` returns the diversity values as a separate vector or matrix. Use `addAlpha()` when you want to keep all data together in one object, and `getAlpha()` when you need the diversity values for immediate use in other calculations. - -**`runMDS()` vs `addMDS()`**: The `runMDS()` function calculates multidimensional scaling coordinates and returns them as a separate matrix, whereas `addMDS()` calculates the MDS coordinates and stores them directly into the `reducedDim` slot of the TreeSummarizedExperiment object. Using `addMDS()` is generally preferred as it maintains all results within the same data object, making downstream analyses and visualization more straightforward. - -The same add/get pattern is also used for other ordination workflows where available, for example `runNMDS()` / `addNMDS()` in the community similarity chapter. - - ## Transformations in practice Below, we apply relative transformation to counts table. @@ -367,6 +254,108 @@ this new dataset is stored into `altExp`. altExp(tse, "philr") ``` +## Rarefaction {#sec-rarefaction} + +The number of taxa we observe in a sample depends strongly on how deeply the +sample was sequenced. Sequencing depth (or library size, see +[@sec-library-size]) refers to the number of reads obtained for a sample. +In practice, sequencing depth often varies substantially between samples, even +within the same study. + +A common way to visualize this relationship is with a rarefaction curve, which +shows how the number of observed features increases as a function of the number +of reads sampled. + +```{r} +#| label: rarefaction_curve +library(mia) +library(vegan) + +data(GlobalPatterns) +tse <- GlobalPatterns + +mat <- tse |> assay("counts") |> t() +rarecurve(mat, step = 1000) +``` + +If a sample has not been sequenced deeply enough for its curve to approach +saturation, its observed richness will be artificially low. This is not a +biological signal: it is a measurement artifact caused by limited sampling. + +This becomes problematic when samples are sequenced at different depths. Samples +with larger library sizes tend to show more observed taxa simply because more +reads were collected, not because their microbial communities are truly more +diverse. Consequently, variation in sequencing depth can bias commonly used +alpha and beta diversity metrics [@Schloss2023]. + +To inspect the depth variation in the data, we can compute and visualize total +counts per sample: + +```{r} +#| label: library_size_variation +library(scater) +library(miaViz) + +tse <- addPerCellQCMetrics(tse, assay.type = "counts") + +plotBoxplot(tse, col.var = "total", x = "SampleType") +``` + +We can also quantify the spread in sequencing depth: + +```{r} +#| label: depth_variation + +# Calculate variation between highest and lowest sequencing depth +depth_variation <- max(tse[["total"]]) / min(tse[["total"]]) +depth_variation +``` + +This dataset shows up to a `r depth_variation |> round()`-fold difference in +sequencing depth across samples. + +To mitigate bias from uneven sequencing depths, *rarefaction* has been proposed +[@Schloss_2024a; @Schloss_2024b]. Rarefaction standardizes sequencing depth by +randomly subsampling reads so that all retained samples have the same total +number of reads. + +Conceptually, the workflow is: + +1. Choose a rarefaction depth (often the minimum library size across samples, here `r tse[["total"]] |> min()`). +2. Discard samples with fewer reads than this depth. +3. For each remaining sample, randomly subsample exactly that many reads. + +In `mia`, rarefaction can be applied directly to an assay: + +```{r} +#| label: simple_rarefaction +tse_rarefy <- rarefyAssay( + tse, + assay.type = "counts" +) +tse_rarefy +``` + +There has been a long-lasting controversy surrounding the use of rarefaction in +microbial ecology. The main concern is that subsampling would omit data +[@mcmurdie2014waste; @Schloss2024rarefaction2], which is why rarefaction is +generally discouraged for differential abundance testing performed on a single +rarefied table [@McMurdie_and_Holmes_2014]. + +However, this concern is reduced when rarefaction is treated as a repeated +resampling step rather than a one-off subsample. By rarefying to the minimum +library size and repeating the procedure many times, we compute the diversity +metric (e.g., alpha or beta diversity) for each iteration and average across +iterations, yielding a stable estimate at a common sequencing effort +[@Schloss2024rarefaction1; @Schloss2024rarefaction2]. +Patrick Schloss has demonstrated that this repeated rarefaction procedure is +"the only method that could control for variation in uneven sequencing effort +when measuring commonly used alpha and beta diversity metrics" [@Schloss2023]. + +Section [@sec-rarefaction-alpha] demonstrates rarefaction for calculating alpha +diversity, and Section [@sec-rarefaction_ordination] extends the approach to +beta diversity. + ::: callout-tip ## Summary