Bio Applied Qiime2 16s

QIIME2 16S Amplicon Workflow

When to Use

• Analyzing 16S rRNA amplicon sequencing data (paired-end FASTQ) from gut, soil, or environmental samples
• Building an import → denoise → classify → diversity → differential-abundance pipeline with QIIME2 .qza/.qzv artifacts
• Choosing DADA2 truncation lengths from quality profiles, or picking a SILVA classifier for a primer pair
• Computing alpha diversity (Shannon, Faith's PD), beta diversity (Bray-Curtis, UniFrac), and PCoA ordination
• Testing which taxa differ between groups (e.g., healthy vs. disease) while accounting for compositional bias

Version Compatibility

QIIME2 2024.5–2024.10 (q2cli), DADA2 plugin (q2-dada2), SILVA 138 reference, ANCOM-BC via q2-composition. Python analysis snippets use numpy ≥1.26, pandas ≥2.0, scipy ≥1.11, scikit-learn ≥1.3, statsmodels ≥0.14.

Prerequisites

• Install QIIME2 in its own conda env: conda env create -n qiime2 --file qiime2-amplicon-2024.10-py310-linux-conda.yml
• Paired-end demultiplexed FASTQ files + a sample manifest CSV; a metadata TSV with sample IDs and group labels
• Familiarity with amplicon sequencing basics (primers, ASV vs. OTU) — see bio-applied-microbial-diversity

Import Reads

Manifest CSV maps sample IDs to absolute FASTQ paths:

sample-id,absolute-filepath,direction
Sample1,/data/Sample1_R1.fastq.gz,forward
Sample1,/data/Sample1_R2.fastq.gz,reverse
Sample2,/data/Sample2_R1.fastq.gz,forward
Sample2,/data/Sample2_R2.fastq.gz,reverse

Goal: get raw FASTQ into a provenance-tracked QIIME2 artifact and pick DADA2 truncation lengths. Approach: import with the manifest, summarize quality, then read off the position where median quality drops below Q20.

qiime tools import \
    --type 'SampleData[PairedEndSequencesWithQuality]' \
    --input-path manifest.csv \
    --output-path reads.qza \
    --input-format PairedEndFastqManifestPhred33

# Check quality to determine DADA2 truncation lengths
qiime demux summarize --i-data reads.qza --o-visualization reads_summary.qzv
# Open reads_summary.qzv at view.qiime2.org — per-position quality boxplots

import numpy as np


def suggest_truncation(mean_qual, positions, q_threshold=25):
    """Return the first read position where mean quality drops below q_threshold.

    mean_qual: 1D array of mean Phred quality per position (as reported by
    `qiime demux summarize`). positions: matching 1-based position array.
    """
    below = np.where(mean_qual < q_threshold)[0]
    return int(positions[below[0]]) if below.size else int(positions[-1])


positions = np.arange(1, 251)
qual_mean_fwd = 37 - 0.03 * positions - 3 * np.exp(-(250 - positions) / 30)
qual_mean_rev = 36 - 0.06 * positions - 5 * np.exp(-(250 - positions) / 20)
trunc_f = suggest_truncation(qual_mean_fwd, positions)
trunc_r = suggest_truncation(qual_mean_rev, positions)
print(f'--p-trunc-len-f {trunc_f} --p-trunc-len-r {trunc_r}')

Denoising with DADA2

DADA2 produces ASVs (exact error-corrected sequences, single-nucleotide resolution) rather than OTUs (97% similarity clusters). ASVs are reproducible across studies — the same ASV sequence means the same organism in any dataset.

Truncation rule: forward + reverse must overlap ≥20 bp after truncation for merging. For V3-V4 (~460 bp amplicon), 220 + 200 = 420 bp combined gives sufficient overlap. After denoising, verify ≥75% of input reads pass filtering and ≥60% merge successfully.

qiime dada2 denoise-paired \
    --i-demultiplexed-seqs reads.qza \
    --p-trim-left-f 13 --p-trim-left-r 13 \
    --p-trunc-len-f 220 --p-trunc-len-r 200 \
    --p-n-threads 0 \
    --o-table feature_table.qza \
    --o-representative-sequences rep_seqs.qza \
    --o-denoising-stats stats.qza

Taxonomic Classification with SILVA

SILVA 138 is the standard 16S reference. Use the pre-trained Naive Bayes classifier matched to your primer pair and read length (e.g., silva-138-99-nb-classifier.qza for full-length, or a region-specific classifier trained on your primer pair).

qiime feature-classifier classify-sklearn \
    --i-classifier silva138_nb_515_806_classifier.qza \
    --i-reads rep_seqs.qza \
    --o-classification taxonomy.qza \
    --p-n-jobs 1

qiime taxa barplot \
    --i-table feature_table.qza --i-taxonomy taxonomy.qza \
    --m-metadata-file metadata.tsv --o-visualization taxa_barplot.qzv

# Remove host/contaminant sequences before diversity analysis
qiime taxa filter-table \
    --i-table feature_table.qza --i-taxonomy taxonomy.qza \
    --p-exclude Mitochondria,Chloroplast,Eukaryota \
    --o-filtered-table feature_table_filtered.qza

Confidence threshold: default --p-confidence 0.7 means assignments below 70% bootstrap confidence are reported as unclassified. Species-level accuracy is poor; genus-level is reliable.

Diversity Analysis

Phylogenetic metrics (Faith's PD, UniFrac) need a tree built from the ASV sequences first.

qiime phylogeny align-to-tree-mafft-fasttree \
    --i-sequences rep_seqs.qza \
    --o-alignment aligned_rep_seqs.qza \
    --o-masked-alignment masked_aligned.qza \
    --o-tree unrooted_tree.qza \
    --o-rooted-tree rooted_tree.qza \
    --p-n-threads 8

qiime diversity core-metrics-phylogenetic \
    --i-phylogeny rooted_tree.qza \
    --i-table feature_table_filtered.qza \
    --p-sampling-depth 5000 \
    --m-metadata-file metadata.tsv \
    --output-dir diversity_metrics/

Produces alpha diversity (faith_pd_vector.qza, shannon_vector.qza, observed_features_vector.qza, evenness_vector.qza) and beta diversity distance matrices (bray_curtis, unweighted_unifrac, weighted_unifrac) plus PCoA ordinations, ready for qiime diversity beta-group-significance (PERMANOVA).

Differential Abundance with ANCOM-BC

Relative-abundance data are compositional (they sum to a constant total depth), so a real increase in one taxon forces apparent decreases in others. ANCOM-BC (Analysis of Compositions of Microbiomes with Bias Correction) estimates and corrects each sample's sampling-fraction bias before testing.

qiime composition ancombc \
    --i-table feature_table_filtered.qza \
    --m-metadata-file metadata.tsv \
    --p-formula 'condition' \
    --o-differentials differentials.qza

qiime composition da-barplot \
    --i-data differentials.qza \
    --p-significance-threshold 0.05 \
    --o-visualization da_barplot.qzv

Output columns: lfc (log-fold change), se, W (test statistic), p_val, q_val (BH-adjusted). Filter on q_val < 0.05 for significant taxa.

def call_significant_asvs(ancom_df, q_threshold=0.05):
    """Filter an ANCOM-BC differentials table to significant, sorted hits.

    ancom_df: DataFrame with columns 'lfc', 'q_val' (as exported from
    `differentials.qza`). Returns hits sorted by |lfc| descending.
    """
    sig = ancom_df[ancom_df['q_val'] < q_threshold].copy()
    return sig.reindex(sig['lfc'].abs().sort_values(ascending=False).index)

Pitfalls

• Rarefaction depth: --p-sampling-depth drops every sample below that read count entirely — check the DADA2 feature table summary first to avoid losing samples
• Compositional data: never run a t-test/ANOVA directly on relative abundances; use ANCOM-BC or CLR-transform first
• Classifier mismatch: a SILVA classifier trained on the wrong primer region/read length gives systematically biased taxonomy calls
• Multiple testing: always apply FDR correction (Benjamini-Hochberg, q_val) across the hundreds of ASVs tested

See Also

• bio-microbiome-qiime2-workflow — end-to-end QIIME2 pipeline wrapper
• bio-microbiome-diversity-analysis — deeper alpha/beta diversity statistics
• bio-microbiome-differential-abundance — ANCOM-BC, ALDEx2, and other DA methods compared
• bio-applied-taxonomic-profiling — shotgun metagenomics alternative to 16S