omicverse-microbiome-16s-amplicon-dada2

name

omicverse-microbiome-16s-amplicon-dada2

description

16S rRNA amplicon analysis from raw FASTQs to a samples × ASVs AnnData with 7-rank SINTAX taxonomy, plus the canonical alpha / beta / ordination / DA stack. Use when running the vsearch-or-DADA2 amplicon pipeline (`ov.alignment.amplicon_16s_pipeline`), when ingesting an existing OTU/ASV count table with `build_amplicon_anndata`, or when computing Shannon / Bray-Curtis / PCoA / Wilcoxon-DA on the resulting AnnData.

OmicVerse Microbiome — 16S amplicon end-to-end

Goal

Take paired-end Illumina FASTQs from a 16S rRNA gene survey through to an analysis-ready AnnData (samples × ASVs, with 7-rank SINTAX taxonomy in var and the ASV centroid sequence stored alongside), plus the de-facto downstream stack: alpha-diversity (Shannon / observed OTUs / Simpson) → beta-diversity (Bray-Curtis distance matrix) → PCoA / NMDS ordination → differential abundance (Wilcoxon by default; pyDESeq2 / ANCOM-BC available).

Cover two backends behind the same wrapper:

• vsearch / UNOISE3 — the default; fast (single-threaded C, parallelised via jobs); slightly more reads retained per sample on noisy datasets.
• DADA2 — Bayesian denoising; longer to fit but produces tighter ASV-level resolution. Same wrapper, just backend='dada2'.

Plus the stepwise API for users who want to inspect or substitute individual stages (merge → filter → dereplicate → UNOISE3 → uchime3 → SINTAX → OTU table → AnnData).

Quick Workflow

1. Reference DB: fetch a SINTAX-formatted 16S reference once with ov.alignment.fetch_rdp(db_dir=...) — returns the path to the .fa.gz you'll feed all subsequent calls.
2. Sample metadata: build a pd.DataFrame keyed by sample-id with at least the phenotype/factor column (here called 'group'). The wrapper writes this into adata.obs.
3. One-shot pipeline: call ov.alignment.amplicon_16s_pipeline(fastq_dir=..., workdir=..., db_fasta=..., backend='vsearch'|'dada2', sample_metadata=meta, ...). Auto-discovers R1/R2 from FASTQ names, runs the full chain, writes intermediates under workdir, returns an AnnData.
4. Validate ingest: confirm adata.shape == (n_samples, n_ASVs), 'phylum' in adata.var.columns (and the rest of the 7-rank), int(adata.X.sum()) is reasonable (>50 % of input reads typically retained).
5. Alpha diversity: ov.micro.Alpha(adata, rarefy_depth=min_depth).run(metrics=['shannon', 'observed_otus', 'simpson']) — rarefies first, writes per-sample metrics into adata.obs.
6. Beta + ordination: ov.micro.Beta(adata, rarefy_depth=min_depth).run(metric='braycurtis') then ov.micro.Ordinate(adata, dist_key='braycurtis').pcoa(n=3). Distance matrix lives at adata.obsp['braycurtis']; PCoA coords at adata.obsm['braycurtis_pcoa']; variance-explained at adata.uns['micro']['braycurtis_pcoa_var'].
7. Differential abundance: ov.micro.DA(adata).wilcoxon(group_key='group', group_a, group_b, rank='phylum'|'genus'|...). Returns DataFrame with effect size + BH-FDR.

Interface Summary

Reference DB:

• ov.alignment.fetch_rdp(db_dir=None, overwrite=False) → str — alias for fetch_sintax_ref('rdp_16s_v18'). ~6.8 MB. Returns path to .fa.gz.
• ov.alignment.fetch_sintax_ref(db_name, db_dir=None) → str — switch DB ('silva138', 'gtdb', etc.).

Pipeline wrapper (one-shot path):

ov.alignment.amplicon_16s_pipeline(
    fastq_dir: Optional[str] = None,
    samples: Optional[Sequence[(name, fq1, fq2|None)]] = None,
    workdir: Optional[str] = None,
    db_fasta: Optional[str] = None, *,
    primer_fwd: Optional[str] = None, primer_rev: Optional[str] = None,
    backend: str = 'vsearch',          # or 'dada2'
    threads: int = 4,
    jobs: Optional[int] = None,
    merge_max_diffs: int = 10, merge_min_overlap: int = 16,
    filter_max_ee: float = 1.0, filter_min_len: int = 0, filter_max_len: int = 0,
    derep_min_uniq: int = 2,
    unoise_alpha: float = 2.0, unoise_minsize: int = 2,
    chimera_removal: bool = True,
    otutab_identity: float = 0.97,
    sintax_cutoff: float = 0.8, sintax_strand: str = 'both',
    sample_metadata: [pd.DataFrame] = ,
    overwrite:  = ,
) → AnnData

Stepwise vsearch helpers (under ov.alignment.vsearch.*):

• merge_pairs(samples, output_dir, max_diffs=10, min_overlap=16, threads=8, jobs=4).
• filter_quality(merge_res, output_dir, max_ee=1.0, threads=8, jobs=4).
• dereplicate(filt_res, output_dir, min_uniq=2, threads=8).
• unoise3(uniques_fasta, output_dir, alpha=2.0, minsize=2, threads=8).
• uchime3_denovo(asv_fasta, output_dir).
• sintax(asv_fasta, db_fasta, output_dir, cutoff=0.8, strand='both', threads=8).
• usearch_global(combined_fastq, asv_fasta, output_dir, identity=0.97, threads=8).

Composer:

• ov.alignment.build_amplicon_anndata(otutab_tsv, asv_fasta, sintax_tsv=None, sample_metadata=None, sample_order=None) → AnnData — convert pre-computed OTU table + ASV centroids + SINTAX TSV into the canonical AnnData.

Diversity / ordination / DA (ov.micro):

• ov.micro.Alpha(adata, rarefy_depth=None, seed=0).run(metrics=('shannon', 'observed_otus'), write_to_obs=True, tree_key='tree') → pd.DataFrame. Convenience methods: .shannon(), .observed() return individual pd.Series. Faith PD via 'faith_pd' requires a phylogenetic tree at adata.uns[tree_key] (see the phylogeny skill).
• ov.micro.Beta(adata, rarefy_depth=None, seed=0).run(metric='braycurtis'|'jaccard'|'aitchison'|'unifrac_unweighted'|'unifrac_weighted', rarefy=None, tree_key='tree', write_to_obsp=True) → pd.DataFrame. UniFrac requires unifrac package + a tree.
• ov.micro.Ordinate(adata, dist_key='braycurtis').pcoa(n=3, write_to_obsm=True) → pd.DataFrame (eigen-decomposition of the distance matrix; proportion_explained() gives variance fractions). .nmds(n=2, random_state=0, write_to_obsm=True) for non-metric MDS.
• ov.micro.DA(adata).wilcoxon(group_key, group_a=None, group_b=None, rank=None, relative=True, min_prevalence=0.1) → pd.DataFrame with log2fc, pvalue, padj. Per-feature Mann-Whitney U on relative abundances.
• ov.micro.DA(adata).deseq2(group_key, group_a=None, group_b=None, rank=None, min_prevalence=0.1, alpha=0.05) — pyDESeq2 NB-GLM on raw counts.
• ov.micro.DA(adata).ancombc(group_key, rank=None, min_prevalence=0.1, pseudocount=1.0) — ANCOM-BC via skbio≥0.7.1.

Preprocessing helpers:

• ov.micro.rarefy(adata, depth=None, seed=0, drop_shallow=True, save_original=True, copy=False) — subsample counts; original counts cached at adata.layers['raw_counts'] when save_original=True.
• ov.micro.filter_by_prevalence(adata, min_prevalence=0.1, min_count=1, copy=False) — drop rare features.
• ov.micro.collapse_taxa(adata, rank='genus', unassigned_label='Unassigned') — sum ASV counts within taxonomic rank.
• ov.micro.clr(adata, layer_out='clr', copy=False) / ilr(...) — compositional transforms (post pseudo-count).

Boundary

Inside scope:

• Full vsearch / UNOISE3 pipeline.
• Full DADA2 pipeline (same wrapper, backend='dada2').
• Stepwise vsearch composability + build_amplicon_anndata from external outputs.
• 7-rank SINTAX taxonomy assignment with the shipped RDP DB (or any SINTAX-format reference).
• Alpha (Shannon / Simpson / Chao1 / observed OTUs) and beta (Bray-Curtis / Jaccard / Aitchison) diversity.
• PCoA / NMDS ordination on a stored distance matrix.
• Wilcoxon DA at any taxonomic rank.

Outside scope — separate skill:

• DA method comparison across Wilcoxon / pyDESeq2 / ANCOM-BC — see omicverse-microbiome-da-comparison.
• Cross-cohort meta-analysis (combine_studies, meta_da) — see omicverse-microbiome-meta-analysis.
• Phylogenetic tree construction + UniFrac + Faith PD — see omicverse-microbiome-phylogeny.
• Paired microbiome + metabolomics integration (MMvec, paired Spearman / CCA) — see omicverse-micro-metabol-paired.
• Shotgun metagenomics — ov.alignment doesn't ship a metagenomics pipeline; this skill is amplicon-only.

Branch Selection

backend='vsearch' vs 'dada2'

• vsearch (default): faster (~minutes for the mothur SOP demo), ships a single binary, parallelises across samples (jobs). Use for production runs and CI.
• DADA2: Bayesian denoising; tighter ASV resolution at the cost of longer runtime. Use when noise floor matters (e.g. low-biomass samples) or when reviewers expect DADA2 specifically.
• Both write the same canonical AnnData layout — switching backends doesn't change downstream code.

Primer trimming (primer_fwd / primer_rev)

• Set both to the actual primer sequences for cutadapt to run first.
• Pass None only when you've confirmed primers are already trimmed (e.g. the mothur SOP test dataset).
• Mismatch: cutadapt errors out with a useful message — don't silence it.

filter_max_ee (expected error rate)

• 1.0 — strict (vsearch default; loses ~30 % of low-quality reads).
• 2.0 — DADA2-tutorial default (more permissive); needed on Q20+ MiSeq runs to retain enough reads.
• Empirical: try 1.0 first, drop to 2.0 if fewer than 50 % of input reads survive the filter.

unoise_minsize and derep_min_uniq

• unoise_minsize=2 (defaults to 2) — UNOISE3 retains an ASV only if it appears at least this many times across the whole derep table. Lower (1) for low-biomass studies; higher (4–8) to suppress sequencing noise on deep cohorts.
• derep_min_uniq is the per-sample equivalent before pooling.

sintax_cutoff (taxonomic confidence)

• 0.8 (default, RDP recommendation) — bootstrap support threshold; ranks below this become empty strings.
• 0.5 — looser; useful on novel environmental samples where genus-level resolution isn't available in any DB.
• Lower than 0.5 — the published "do-not-go-below" floor; results are not reproducible across DB versions.

otutab_identity

• 0.97 — classical OTU clustering threshold (≈ species-level for 16S V3-V4); use when you specifically want OTU-style buckets, not ASVs.
• For pure ASV workflows, this stage is just used to build the per-sample count matrix from the original reads — keep at 0.97 for compatibility with downstream tools.

Alpha metrics

• shannon — diversity (counts both richness and evenness); always report this.
• observed_otus — pure richness, depth-sensitive (always rarefy first).
• simpson — diversity weighted toward dominant taxa.
• chao1 — depth-corrected richness estimator; reports unobserved-taxa correction.
• faith_pd — phylogenetic diversity; needs a tree (phylogeny skill).

Beta metrics

• braycurtis — abundance-weighted, the de-facto 16S default.
• jaccard — presence/absence; complement to Bray-Curtis when interpretation differs.
• aitchison — CLR-Euclidean; correct for compositional data, but harder to interpret intuitively.
• unifrac_* — phylogenetically aware; needs a tree (phylogeny skill).

Ordination — PCoA vs NMDS

• PCoA: linear, fast, gives variance-explained — use when distances are ≈ Euclidean (Aitchison) or for the canonical "PC1/PC2 with %variance" plot.
• NMDS: rank-based, no eigenvalues; less distorted on Bray-Curtis / Jaccard. Report the final stress (in adata.uns['micro']).

DA method

• wilcoxon — non-parametric, fast, no distributional assumptions. Default for moderate cohorts (n>=10/group).
• deseq2 — NB-GLM on raw counts; better-calibrated than Wilcoxon, more conservative with small counts. Requires pip install pydeseq2.
• ancombc — bias-corrected ANCOM; closest to the compositional ground-truth but slowest. Requires skbio>=0.7.1. See the DA-comparison skill for trade-offs.

Input Contract

• Pipeline ingest: paired-end Illumina FASTQs with predictable R1/R2 naming (e.g. <sample>_R1.fastq.gz); one sample per pair. The wrapper auto-discovers from fastq_dir; pass samples=[(name, fq1, fq2)] for explicit control.
• workdir: writable directory, >1 GB free space recommended for moderate cohorts. The wrapper does not fall back to $HOME or /tmp — pass an explicit path.
• db_fasta: SINTAX-formatted 16S FASTA (fetch_rdp returns a valid one). Skipping (db_fasta=None) means no taxonomy in var.
• sample_metadata: pd.DataFrame indexed by sample-id (matching FASTQ-derived names exactly). The phenotype column is conventionally called 'group' to match ov.micro.DA's defaults.
• After ingest, adata.X is a sparse (or dense) integer count matrix — never log-transform before alpha / beta / DA; those functions handle scaling internally.

Minimal Execution Patterns

# vsearch one-shot pipeline
import omicverse as ov
import pandas as pd

ov.plot_set()

DB_FASTA = ov.alignment.fetch_rdp(db_dir='./db/rdp')

meta = pd.DataFrame({
    'group': ['Early', 'Early', 'Late', 'Late', 'Mock'],
    'day':   [0, 1, 7, 8, 0],
}, index=['F3D0', 'F3D1', 'F3D7', 'F3D8', 'Mock'])

adata = ov.alignment.amplicon_16s_pipeline(
    fastq_dir='./raw/MiSeq_SOP',
    workdir='./run_vsearch',
    db_fasta=DB_FASTA,
    threads=8, jobs=4,
    primer_fwd=None, primer_rev=None,    # SOP is pre-trimmed
    filter_max_ee=1.0,
    unoise_minsize=2,
    sintax_cutoff=0.8,
    sample_metadata=meta,
)
print(adata)                       # samples × ASVs
print(int(adata.X.sum()), 'reads mapped')

# DADA2 backend — same wrapper
adata_dada2 = ov.alignment.amplicon_16s_pipeline(
    samples=[('F3D0', 'raw/F3D0_R1.fastq.gz', 'raw/F3D0_R2.fastq.gz'), ...],
    workdir='./run_dada2',
    db_fasta=DB_FASTA,
    backend='dada2',
    primer_fwd=None, primer_rev=None,
    filter_max_ee=2.0,        # DADA2 typically wants more permissive EE
    sintax_cutoff=0.8,
    threads=4,
    sample_metadata=meta,
)

# Stepwise (vsearch) — when you need to inspect / substitute a stage
import numpy as np

merge_res = ov.alignment.vsearch.merge_pairs(samples, output_dir='./step/merged',
                                             max_diffs=10, min_overlap=16,
                                             threads=8, jobs=4)
filt_res  = ov.alignment.vsearch.filter_quality(merge_res, output_dir='./step/filtered',
                                                max_ee=1.0, threads=8, jobs=4)
derep     = ov.alignment.vsearch.dereplicate(filt_res, output_dir='./step/derep',
                                             min_uniq=2, threads=8)
unoise    = ov.alignment.vsearch.unoise3(derep['uniques'], output_dir='./step/asv',
                                         alpha=2.0, minsize=2, threads=8)
nochim    = ov.alignment.vsearch.uchime3_denovo(unoise['asv'], output_dir='./step/asv')
tax       = ov.alignment.vsearch.sintax(nochim['asv'], db_fasta=DB_FASTA,
                                        output_dir='./step/taxonomy',
                                        cutoff=0.8, strand='both', threads=8)
otutab    = ov.alignment.vsearch.usearch_global(derep['combined'], nochim['asv'],
                                                output_dir='./step/otutab',
                                                identity=0.97, threads=8)

adata_step = ov.alignment.build_amplicon_anndata(
    otutab_tsv=otutab['otutab'],
    asv_fasta=nochim['asv'],
    sintax_tsv=tax['tsv'],
    sample_metadata=meta,
    sample_order=[s[]  s  samples],
)

# Diversity + ordination + DA on the AnnData
min_depth = int(np.asarray(adata.X.sum(axis=1)).min())

# alpha — writes shannon/observed_otus/simpson into adata.obs
ov.micro.Alpha(adata, rarefy_depth=min_depth).run(
    metrics=['shannon', 'observed_otus', 'simpson'],
)

# beta — distance matrix in adata.obsp['braycurtis']
ov.micro.Beta(adata, rarefy_depth=min_depth).run(metric='braycurtis')

# PCoA — coords in adata.obsm['braycurtis_pcoa'], var-fraction in uns
ord_ = ov.micro.Ordinate(adata, dist_key='braycurtis').pcoa(n=3)
pct = ord_.proportion_explained() * 100.0

# DA — Wilcoxon at phylum level
da = ov.micro.DA(adata).wilcoxon(
    group_key='group', group_a='Early', group_b='Late',
    rank='phylum', min_prevalence=0.1,
)
print(da.head(10))

Validation

• After ingest: adata.shape matches the expected (n_samples, n_ASVs). The mothur SOP demo lands at ~22 samples × ~250 ASVs with vsearch defaults.
• Read retention: int(adata.X.sum()) / int(<input read total>) > 0.5 is healthy. Below 0.3 indicates the filter / chimera / derep stages dropped too much; revisit filter_max_ee and unoise_minsize.
• Taxonomy coverage: (adata.var['phylum'] != '').sum() / adata.n_vars should exceed 0.8 — lower means the DB is wrong species range or sintax_cutoff is too strict.
• Alpha rarefaction: pass rarefy_depth=int(adata.X.sum(axis=1).min()) to rarefy down to the shallowest sample. If that's <1000, the cohort is too shallow for trustworthy alpha — flag this rather than silently downsampling.
• Beta rarefaction: same — uneven depths inflate Bray-Curtis distances against the deeper samples.
• After ordination: proportion_explained() is decreasing; PC1 + PC2 typically capture 30–60 % on 16S. If PC1 alone captures >70 %, suspect a single dominant outlier sample.
• After DA: (da['padj'] < 0.05).sum() should be sane for the cohort — for the mothur SOP 22-sample demo, expect ~3–8 phylum-level hits at padj<0.05.
• The DADA2 backend produces different ASV identities than vsearch; downstream interpretations (e.g. specific ASV names) don't transfer between backends. Diversity / ordination / DA patterns should agree.

Resource Map

• See reference.md for compact one-shot + stepwise + diversity snippets.
• See references/source-grounding.md for verified pipeline + ov.micro API signatures and the docstring backfill log (Alpha.shannon / Alpha.observed / Beta.braycurtis / Ordinate.nmds / Ordinate.proportion_explained / DA.deseq2).
• For DA method comparison (Wilcoxon vs pyDESeq2 vs ANCOM-BC), see omicverse-microbiome-da-comparison.
• For cross-cohort meta-analysis on the resulting AnnDatas, see omicverse-microbiome-meta-analysis.
• For phylogenetic tree + UniFrac + Faith PD, see omicverse-microbiome-phylogeny.
• For paired microbiome + metabolomics analysis, see omicverse-micro-metabol-paired.

Examples

• "Run the vsearch one-shot pipeline on the mothur MiSeq SOP FASTQs with filter_max_ee=1.0 and the RDP v18 SINTAX DB; return an AnnData."
• "Re-run the same FASTQs with backend='dada2' and filter_max_ee=2.0; confirm diversity patterns match the vsearch result."
• "Compute Shannon + observed-OTUs alpha at the shallowest depth, then Bray-Curtis beta + 3-D PCoA, and plot the first two PCs colored by obs['group']."
• "Run Wilcoxon DA between Early and Late at phylum level with min_prevalence=0.1 and report the top 10 hits."

References

• Tutorial notebooks:
  • t_16s_amplicon.ipynb — vsearch end-to-end + alpha/beta/PCoA/DA on the mothur SOP.
  • t_16s_dada2.ipynb — DADA2 backend on the same FASTQs.
• Live API verified — see references/source-grounding.md.