# Bio Applied Qiime2 16s

> Run QIIME2 16S amplicon workflows: import FASTQ, DADA2 denoise to ASVs, SILVA taxonomy, alpha/beta diversity, ANCOM-BC. Use when analyzing 16S/amplicon microbiome data or .qza/.qzv pipelines.

- Skill: `pavel-kravchenko/bio-applied-qiime2-16s` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-qiime2-16s`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-qiime2-16s/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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:

```text
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.

```bash
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
```

```python
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.

```bash
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).

```bash
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.

```bash
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.

```bash
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.

```python
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