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 4047b03a..a3dfedbb 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 b69abe10..48757e04 100644 --- a/inst/pages/quality_control.qmd +++ b/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 diff --git a/inst/pages/transformation.qmd b/inst/pages/transformation.qmd index f4070d47..35e68c71 100644 --- a/inst/pages/transformation.qmd +++ b/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