docs: Add comprehensive Rarefaction section to transformation chapter

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

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:

inst/pages/quality_control.qmd

@@ -165,7 +165,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
@@ -206,7 +206,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
 

inst/pages/transformation.qmd

@@ -254,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