16S rRNA Sequencing Data Analysis: A Complete Guide

16S rRNA sequencing data analysis converts raw amplicon reads into a profile of which bacteria are in a sample and how abundant they are. A modern workflow removes primers, filters and denoises reads into amplicon sequence variants (ASVs), assigns each ASV a taxonomy against a reference database such as SILVA, and then compares samples with diversity and differential-abundance statistics. This guide walks through each step, the decisions that matter, and the mistakes that most often ruin a study.

What is 16S rRNA sequencing?

The 16S ribosomal RNA gene is about 1,500 bp long and is present in every bacterium and archaeon. It contains conserved stretches, where universal PCR primers can bind, separated by nine hypervariable regions (V1–V9) whose sequence differs between taxa. A 16S study amplifies one or two of those regions, such as V4 or V3-V4, from all the DNA in a sample and sequences the amplicons, usually on an Illumina MiSeq.

The result is a FASTQ file per sample, or a pair of files for paired-end runs, containing tens of thousands of reads. On its own a FASTQ file says nothing about biology. The analysis turns it into a table of organisms and abundances you can interpret.

What are the steps in a 16S analysis?

Step 1: Check your input files

Most failed 16S analyses fail here, before any biology. Confirm that:

• Each sample has its own FASTQ file or R1/R2 pair. Demultiplexing by barcode is normally done by the sequencing facility.
• Paired files really are pairs, with matching read counts, and R1 is the forward read. Swapped R1/R2 files are a common, silent error.
• The reads start with the primers you expect. If you analyze V3-V4 reads with V4 primers, nearly every read is discarded.
• Each sample has enough reads. Aim for at least 10,000 reads per sample after filtering.

Step 2: Remove primers

Primer sequences are synthetic. They contain ambiguity codes (for example N, W or Y) and do not reflect the organism's real sequence, so they must be removed before denoising. Cutadapt (Martin, 2011) finds the primer at the start of each read, removes it, and can discard reads where the primer is missing, which are usually off-target products.

Step 3: Filter and truncate reads

Illumina read quality falls toward the end of each read, especially in R2. Truncating reads where quality drops improves denoising, but there is a hard limit: the forward and reverse reads must still overlap enough to merge into one sequence. DADA2 needs at least 12 bp of overlap by default; 20 bp is a safer margin.

For example, the V3-V4 insert between the 341F and 805R primers is about 427 bp. With a 20 bp overlap, the truncated R1 and R2 must together reach 447 bp. That works with 2×300 bp reads and cannot work with 2×150 bp reads, however they are truncated. See choosing a 16S region for each region's requirements.

Step 4: Denoise reads into ASVs

Denoising separates real biological sequences from sequencing errors. DADA2 (Callahan et al., 2016) learns the error profile of each run, then tests whether each less-abundant sequence is more likely to be an error derived from a more abundant one or a real variant. It then merges read pairs and removes chimeras, the artificial hybrid sequences created during PCR.

The output is an ASV table: rows are exact sequences, columns are samples, and values are read counts. ASVs replaced the older practice of clustering reads into 97% OTUs; ASVs vs OTUs explains why.

Always check how many reads survive each stage. Losing most reads at the merge step usually means the reads are too short for the region, or were truncated too aggressively.

Step 5: Assign taxonomy

Each ASV is classified against a curated reference database. SILVA (Quast et al., 2013) is the most widely used for 16S; release 138.2 is current. The standard method is a naive Bayesian classifier (Wang et al., 2007), which reports a confidence for each rank and stops at the deepest rank it is confident about.

Expect reliable assignments to genus. Species-level assignment is possible only when an ASV exactly matches a single species in the database, which is uncommon for short regions such as V4. Be cautious with species names from 16S data, and note that SILVA follows current nomenclature: for example Lactobacillus was split into several genera in 2020, so L. fermentum appears as Limosilactobacillus fermentum.

Step 6: Build a phylogeny

Aligning ASVs (for example with MAFFT) and inferring a tree (for example with FastTree) lets you use phylogeny-aware metrics, such as Faith's phylogenetic diversity and UniFrac distances, which count closely related ASVs as more similar than distant ones.

Step 7: Measure diversity

Alpha diversity describes a single sample:

Beta diversity compares samples. Bray-Curtis distance uses abundances; Jaccard uses presence and absence. Principal coordinates analysis (PCoA) plots the distances, and PERMANOVA (Anderson, 2001) tests whether groups, such as treatment and control, differ in overall composition.

Rarefaction curves show whether sequencing depth was enough to capture each sample's diversity. Whether to rarefy before computing diversity is debated (McMurdie & Holmes, 2014); report what you did.

Step 8: Test for differential abundance

Sequencing gives relative abundances: each sample's counts sum to a total set by the sequencer, not by the biology. This compositional nature means that a simple t-test on relative abundances can report false differences (Gloor et al., 2017). Methods built for compositional data, such as ANCOM-BC2 (Lin & Peddada, 2024), correct for this and should be used to claim that a taxon differs between groups.

Step 9: Predict function

16S data does not measure genes, but PICRUSt2 (Douglas et al., 2020) can predict likely metabolic pathways from the reference genomes of related organisms. Treat the results as hypotheses. Predictions are least reliable for environments, such as soil, whose organisms are poorly represented by sequenced genomes.

How many reads per sample do you need?

There is no single number, but most studies aim for at least 10,000 high-quality reads per sample after filtering. Low-diversity samples such as some infant gut samples need fewer; high-diversity soils need more. Check the rarefaction curves: if they are still climbing steeply, the sample was under-sequenced. Good's coverage is not a useful check on DADA2 output, because denoising removes the singletons it is calculated from.

Why do 16S analyses fail?

• Wrong primers selected: almost all reads are discarded at primer removal.
• Reads too short for the region: read pairs cannot merge, so most are lost.
• Swapped R1 and R2 files: primers are not found where expected.
• Too few reads: diversity estimates become unreliable.
• No negative controls: in low-biomass samples, reagent contaminants can dominate the results (Salter et al., 2014). Sequence extraction and PCR blanks alongside your samples.
• Batch effects: samples from different runs or extraction kits differ for technical reasons. Randomize samples across runs.

Which tools can you use for 16S analysis?

Running this workflow on BioAnalysis.ca

BioAnalysis.ca runs every step above automatically. You upload FASTQ files and choose a primer preset, and the pipeline checks primers, read orientation and merge length before it runs. It then returns an interactive report, figures, the ASV and taxonomy tables and a methods paragraph you can adapt from the how to cite page. Data is processed and stored in Canada.

References

1. Anderson MJ. A new method for non-parametric multivariate analysis of variance. Austral Ecology. 2001;26:32–46.
2. Bolyen E, Rideout JR, Dillon MR, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nature Biotechnology. 2019;37:852–857.
3. Callahan BJ, McMurdie PJ, Rosen MJ, et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nature Methods. 2016;13:581–583.
4. Douglas GM, Maffei VJ, Zaneveld JR, et al. PICRUSt2 for prediction of metagenome functions. Nature Biotechnology. 2020;38:685–688.
5. Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. Microbiome datasets are compositional: and this is not optional. Frontiers in Microbiology. 2017;8:2224.
6. Lin H, Peddada SD. Multigroup analysis of compositions of microbiomes with covariate adjustments and repeated measures. Nature Methods. 2024;21:83–91.
7. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal. 2011;17:10–12.
8. McMurdie PJ, Holmes S. Waste not, want not: why rarefying microbiome data is inadmissible. PLoS Computational Biology. 2014;10:e1003531.
9. Quast C, Pruesse E, Yilmaz P, et al. The SILVA ribosomal RNA gene database project. Nucleic Acids Research. 2013;41:D590–D596.
10. Salter SJ, Cox MJ, Turek EM, et al. Reagent and laboratory contamination can critically impact sequence-based microbiome analyses. BMC Biology. 2014;12:87.
11. Schloss PD, Westcott SL, Ryabin T, et al. Introducing mothur. Applied and Environmental Microbiology. 2009;75:7537–7541.
12. Wang Q, Garrity GM, Tiedje JM, Cole JR. Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Applied and Environmental Microbiology. 2007;73:5261–5267.