Hi all,
I know this has been talked about this forum many, many times, so I don't want to beat a dead horse, but I do want to organize my thoughts into one place. There's way too much information for beginners like myself to look at. My intentions for this post are to 1) put myself out there to be corrected if I'm misunderstanding something, 2) provide a compiled list of literature/plugins/forum threads for others to reference in the future, and 3) poll the forum to see if there's a community favorite.
To be clear, I'm not trying to advocate for any particular method... I don't even know which method to use for my own analysis
. I'm also not a biostatician, so I'm struggling to make a decision. It seems like every source I read advocates for something different. It's worth noting that I work with 16S data, not RNA-seq. I also haven't really used R yet. (Still learning how it works.)
Pipeline Context: You imported your data to QIIME2. You denoised with dada2 or deblur. You performed your taxonomic classification and aligned your sequences. You're looking at your feature-table summarize results and questioning everything before you move onto analysis. Sample A has 1k reads, sample B has 50k reads, and sample C has 100k reads. How do we control for the frequency of reads in each sample so that none of them have a stronger impact on our analysis than the others?
Quick Explanations (as I understand them - please correct me if I'm wrong and add additions if there's anything I left out):
-
Rarefying (aka single subsampling as @zippyzo and @colinbrislawn call it):
-
How it works: You choose a minimum number of reads that a sample must have to be included in your analysis, then qiime2 randomly picks that number of reads out of each sample to use for diversity comparisons. qiime2 filters out all the excess reads.
-
Why: Ensures that all samples have the same number of reads so that alpha and beta analyses can be performed. This helps weigh high-frequency and low-frequency sample data equally in your analysis. You can also discard samples with too few reads from your analysis without making a new metadata file or going back through earlier steps.
-
Command that uses it: qiime diversity core-metrics --p-sampling-depth
- Note: --p-sampling depth can be chosen based on lots of different information, such as alpha rarefaction curves, the minimum read count in your samples, a very general predetermined number (i.e. 2,000), the sample type, the questions you're trying to answer, the type of downstream analysis, and your confidence in the quality of the reads, and the number of samples in each treatment group... in other words, this is highly variable between studies. There is not a standardized approach that I've encountered.
-
Example Tutorial: Moving Pictures Tutorial
Pros Cons The number of reads is even across samples. You potentially exclude a ton of reads from your analysis, so the significance of your results can change depending on the sample depth. (See McMurdie & Holmes 2014) Tons of downstream analysis techniques use rarified feature tables (ie alpha and beta diversity). Sampling depth isn't standardized across studies, so results could be hard to interpret. Very commonly used across literature, so you're less likely to be questioned for using it if a person with a non-metagenomics background reviews your work Adds artificial zeroes* to your data Increases the false discovery rate in datasets with a larger number of samples, 10x+ uneven library sizes, or compositional effects (see Weiss 2017). -
*Artificial zeroes: Artificial zeroes are made in many ways. For example, when a taxon is present in your sample but its abundance is below the limit of detection by your sequencing method, your data might reflect that there is absolutely none of that taxon in the sample. This is an artificial zero. When rarefaction excludes reads beyond the sampling depth, and a rare taxon is filtered out, this also creates an artificial zero. It is really hard to determine which zeroes are real and which ones are artificial, so it is generally not recommended to add more zeroes with your sequence depth control method. (@jwdebelius discusses that here.)
Rarefaction (aka multiple subsampling/bootstrapping):
- How it works: Works generally the same as rarefying in the sense that you choose a sampling depth, but the rarefying is done more than once. qiime2 randomly picks out the number of reads from each of your samples
nnumber of times to makennumber of rarefied tables withnnumber of diversity analyses. You can merge your rarefied tables together and find the average diversity metrics in downstream steps. (Graphic summary from @gregcaporaso here). - Why: Like rarefying, rarefaction helps to make sure that all your samples are the same size for alpha and beta diversity analyses. However, rarefaction attempts to capture a more complete picture of what your samples look like as a whole by decreasing the chance that an important feature will be captured with one random subsampling and not another.
- Plugin that uses it: q2-boots
- Example Tutorial: Gut-to-Soil Axis Tutorial
Pros Cons Similar to rarefying, the number of reads are made equal across samples The use of rarefaction still excludes a number of reads from analysis, so the statistical impact of your results can still be affected by sampling depth. (See McMurdie & Holmes 2014). Tons of downstream analysis techniques use rarefied feature tables (ie alpha and beta diversity) Sample depth isn't standardized across studies, so results could be hard to interpret. Rarefaction uses repeated subsampling to reduce the chances of missing important taxa or over-emphasizing unimportant taxa. Some sources say that the diversity metrics resulting from repeated subsampling are only negligibly different from single subsampling, so it may not be necessary. Another con to add later --- adds artificial zeroes to your data.
-
-
Center Log Ratio (CLR)
- How it works: Divides each feature's read count by the mean of all counts in the sample, then takes the logarithm. This transforms the data so that each feature's abundance is relative to the mean rather than an absolute quantity.
- Why: Because of the nature of sequencing technology available today, sequencing data typically reflects the relative quantities of organisms in a given sample rather than the exact quantities of organisms that were present. This is at odds with how reads are perceived by computer software (which treats the frequency of a feature as an absolute number). This miscommunication often complicates the analysis process and perpetuates a lack of clarity in literature, as discussed in this forum thread, Bindels (2025), Weiss (2017), and many other places. To combat this miscommunication, the CLR technique aims to transform the data into relative values. This helps to make sure that downstream analyses are not overly influenced by compositionality (ie the influence of one taxon's abundance in a sample on the relative abundance of another; see Gloor 2017 for more details).
- Plugin that uses it: gemelli (formerly DEICODE)
- Note: ANCOM uses an additive version of the log-ratio, so the technique is similar to CLR in the sense that the numbers are transformed into relative quantities. See Quinn 2019 for a better explanation of this concept.
- Followup note: ALDEx2 also uses the log ratio approach, but rather than using the mean of all features as a reference value, it uses the mean of a pre-specified feature (Lin & Peddada 2020)
- Example Tutorial: RPCA-moving-pictures.ipyb
- Pros:
- Treats relative abundance data as relative abundance data.
- Since a sampling depth is not typically used, it does not reduce the library size or throw out reads
- Cons:
- Not recommended in combination with rarifaction because that adds more artificial zeroes to the data. Zeroes don't work with log transformations.+
- Makes some features' relative abundance values negative, which means you can't use the transformed feature-table output for alpha/beta diversity analyses.
- +This could potentially be mitigated by adding a pseudocount to your data before performing the transformation, but then you run the risk of turning a real zero into a false non-zero abundance. As mentioned earlier, it's hard to know which zeroes are real and which ones are caused by shallow sequencing. Some plugins have pre-built-in zero handling (as @jwdebelius mentions in connection to DEICODE here), so you also don't want to double the pseudocount already added.
-
Conditional Filtering:
- How it works: You determine the minimum threshold of abundance/prevalence that a feature must have to remain in your analysis. For example, in your farmer's market, workers must report at least 1% of apples in at least 10% of your vegetable carts for you to feel confident that the apples are actually there. Any feature that doesn't meet this threshold is filtered out before you do your analysis. (Math explained here by @jwdebelius, or you can see the effects of 10% prevalence filtering on the identification of significant taxa in a graph from Nearing 2022)
- Why: This helps to cut out features that aren't being detected consistently across samples. See @jwdebelius's comment here. In a way, this helps cut out features you think might be misidentified/contaminants/unimportant for your analysis. It could also change the sensitivity to differentially abundant taxa in your downstream tests (metagenomeSeq, MaAsLin, corncob, ANCOM, ALDEx, etc.), although that could be either a good or bad thing depending on your perspective.
- Command that uses it/example: feature-table filter-features-conditionally
- Pros:
- Doesn't introduce more zeroes to the data.
- It cuts out the features that we are least confident in rather than cutting down the number of frequencies to a specific sampling depth. This helps to preserve statistical power.
- Cons:
- Not discussed in many forums here.
- Not standardized practice across studies (at least as far as I know).
- It's hard to know if the things you're cutting out are actually contaminants/misidentified species or if they're just really rare taxa. Rare taxa definitely get cut out with this mechanism, so this could tank your diversity metrics.
-
Combining these methods within the same pipeline (ex - Rarefaction + CLR):
- Generally not recommended because many methods add artificial zeroes, which aren't easy to distinguish from real zeroes. See notes from @jwdebelius's replies on another thread.
-
Using RNA-Seq sequence depth control methods (like DESeq2 or edgeR) for metagenomics data
- Not generally recommended because they inflate the false discovery rates, although ALDEx2 might be an exception (Quinn 2019, Lin & Peddada 2020, and Nearing 2022).
- See the linked sources for further explanations. They talk about how RNA-seq methods usually assume that the absolute abundances of transcripts are pretty equal across samples, and only a few transcripts will be different. In contrast, rRNA datasets usually have heterogenous abundances across samples, so the typical control methods used for RNA-seq apply false assumptions to the data. ALDEx2 controls for this false-positive rate by acknowledging the uncertainty of low-frequency samples, but it's still a little hard to use because you have to pre-specify a reference taxon.
- Not generally recommended because they inflate the false discovery rates, although ALDEx2 might be an exception (Quinn 2019, Lin & Peddada 2020, and Nearing 2022).
Here's a battle of the minds presented in literature:
Generally Anti-Rarefying:
- McMurdie & Holmes (2014) - Notes that rarefying data causes many reads to be excluded from analysis, so the statistical significance of findings are affected and there is a high rate of false discovery. Although highly critiqued by Schloss and others, they bring up good points to consider regarding data exclusion.
Generally Pro-Rarefaction:
- Schloss (2024) (and his preceding paper, Schloss 2023) - Notes the potential errors in the McMurdie and Holmes study design, then performs a ton of analyses to show the effects of rarefaction on various diversity metrics and dataset types. Rarefaction impacts some metrics more than others.
Other Literature (formal & informal):
- Gloor (2017) - Microbiome datasets are compositional: And this is not optional
- Lifechanging paper for considering this question. Worth the time to read.
- Helped me understand why q2-composition is probably named the way it is.
- Weiss (2017) - Normalization and microbial differential abundance strategies depend upon data characteristics
- Also discusses the "scaling" method, parametric tests, and normalization techniques, which I did not touch on here.
- Willis & Bunge (2015) - Estimating diversity via frequency ratios
- Discusses the breakaway package in R (which has a qiime2 plugin). I did not touch on it here because I'm running out of steam.
- Quinn (2019)
- Discusses the importance of choosing a method based on the data type (rRNA vs cDNA vs RNA-seq, etc.), then discusses their RNA-seq pipeline.
- Lin & Peddada (2020) - Analysis of microbial compositions: a review of normalization and differential abundance analysis
- This review does a nice job of explaining the math underlying differential abundance analysis techniques. Also explains scaling, differential ranking, and balance-based techniques, which I didn't cover here.
- Nearing (2022) - Microbiome differential abundance methods produce different results across 38 datasets
- Provides analyses of many datasets with many different methods. Nice way to see the results of different techniques compared directly.
- Tedersoo (2022) - Best practices in metabarcoding of fungi: from experimental design to results
- While this focuses on fungal research, the study design components from extraction to analysis are clearly explained, and many of the same challenges encountered in 16S metagenomics are discussed.
- Bindels (2025) - A blueprint for contemporary studies of microbiomes
- Clarifies a lot of the language & practices used in metagenomics literature & endorses the STREAMS guidelines for future works.
- Nice Reddit threads:
- Rarefaction vs other normalization (Lots of community discussion & more links)
- Which transformation method to use with microbiome counts data for addressing differences in sequencing depth and compositionality? (OP discusses a lot of other biostatistical methods I didn't touch on here because I'm not a trained biostatician.)
Finally, just so I can get a vibe check from the community, which method of sampling depth control do you prefer for 16s analysis?
- Rarefying
- Rarefaction
- Control Log Ratio (CLR)
- Conditional Filtering
- Frequency Ratios
- Mixture of these methods (same pipeline)
- Other method not described here (I'll describe in reply)
- Idk, I'm new here
- No favorite / Depends on the study