Bio Workflows Microbiome Pipeline

Version Compatibility

Reference examples tested with: DADA2 1.30+, cutadapt 4.6+, phyloseq 1.46+, vegan 2.6+, ALDEx2 1.34+, ANCOMBC 2.4+, QIIME2 2024.2+, PICRUSt2 2.5+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

The error model is a PER-RUN artifact and the reference database is a versioned dependency, not just a tool: run learnErrors once per sequencing run, never pooled; record the SILVA/GTDB/UNITE release and the classifier training region, the tree-build method (SEPP/Greengenes2 vs de novo), the rarefaction depth, and the PICRUSt2 reference release alongside every result.

Microbiome Pipeline

"Run my 16S amplicon study end to end from FASTQ" -> Stage primer removal, per-run denoising, region-matched taxonomy, a placed tree, declared-depth diversity, and a consensus differential-abundance result - deferring every per-step decision to the owning microbiome skill, because each stage's parameters silently set what the next stage can find.

Complete workflow from demultiplexed amplicon FASTQ to a confidence-graded differential-abundance result. This is an ORCHESTRATION skill: it sequences the stages and routes the science to the six microbiome category skills; it does not re-teach the per-step decisions.

Scope: amplicon (16S/ITS) reads -> ASV table -> taxonomy -> diversity + consensus DA + optional predicted function. Shotgun/WGS reads -> workflows/metagenomics-pipeline. Each method choice -> its owning microbiome skill (links below). QIIME2 artifact/provenance route -> microbiome/qiime2-workflow.

The Single Most Important Modern Insight -- The Pipeline Is a Chain of Modeling Choices, Not a Conveyor Belt

The feature table, the diversity number, and the differential-taxa list are not observations of the community - each is the residue of a knob turned at an EARLIER stage, before any result existed. Turn the knobs differently and the answer changes. The orchestration discipline is therefore to declare and defend each stage's choice, in order, because every stage silently constrains the next:

1. What is stripped and where reads are truncated decides what is detectable. Primers left on corrupt the error model and inflate chimeras; truncating below the merge-overlap budget erases taxa by arithmetic. These are set before any ASV exists (microbiome/amplicon-processing).
2. The error model is per-run, the tree is a model, and the depth is a sample-deletion knob. Pooling runs into one learnErrors denoises wrong; a de novo tree from short reads injects topology noise into UniFrac; a sampling depth above some samples' totals silently drops the lowest-biomass samples (microbiome/diversity-analysis).
3. Rarefy for diversity, never for DA; and which taxa are "significant" depends more on the tool than the biology. Keep the raw counts, rarefy only into the diversity branch, and report a CONSENSUS of >=2 compositionally-aware tools rather than one tool's hit list (Nearing 2022 Nat Commun 13:342; microbiome/differential-abundance).
4. PICRUSt2 predicts potential, not activity, and is circular with taxonomy. Predicted function is the ASV table re-encoded through a fixed lookup - hypothesis-generating, gated on NSTI, never "more active" (microbiome/functional-prediction).

Pipeline Stages

Each stage hands its decisions to the owning skill; this table is the routing map, not a parameter sheet.

QIIME2 artifact/provenance alternative for the whole chain: microbiome/qiime2-workflow (qiime cutadapt -> dada2 denoise-paired -> feature-classifier classify-sklearn -> fragment-insertion sepp -> diversity core-metrics-phylogenetic -> composition ancombc -> q2-picrust2). Same decisions, .qza/.qzv provenance ecosystem.

Workflow Overview

Demultiplexed amplicon FASTQ (per sample, per run; primers/region known)
    |
    v
[0. Read QC]                         FastQC/MultiQC -> read-qc/quality-reports
    |
    v
[1. Remove primers]  cutadapt -g FWD -G REV --discard-untrimmed   (BEFORE truncation)
    |
    v
[2. Per-RUN DADA2]   filterAndTrim -> learnErrors (one run) -> dada -> mergePairs
    |                (repeat per sequencing run; never pool FASTQs into one learnErrors)
    v
[3. Combine + chimeras]  mergeSequenceTables(run1, run2, ...) -> ONE removeBimeraDenovo
    |
    v
[4. Taxonomy]        region-matched SILVA/GTDB/UNITE classifier -> genus (16S)
    |
    v
[5. phyloseq + tree] otu+tax+sample_data + SEPP/Greengenes2 placed tree (NOT de novo)
    |
    +--> [6. Diversity]  rarefy_even_depth(declared depth) -> alpha; UniFrac/Bray + adonis2 + betadisper
    |        (report the depth AND the dropped-sample list)
    |
    +--> [7. Differential abundance]  UNRAREFIED counts -> ALDEx2 AND ANCOM-BC2/LinDA -> CONSENSUS
    |
    +--> [8. Functional prediction]  PICRUSt2 (optional) -> KO/MetaCyc POTENTIAL, report NSTI
    |
    v
ASV table + taxonomy + diversity + consensus DA hits (+ predicted potential)

Stage 1 to 3: Primers, Per-Run Denoising, Chimeras

Goal: Turn demultiplexed reads into one chimera-free ASV table, with primers off first and the error model fit per run.

Approach: Strip primers with cutadapt before any quality step; run the DADA2 block (filter -> learnErrors -> dada -> mergePairs) ONCE PER sequencing run; merge the run-level tables by exact sequence string; remove chimeras once on the combined table. Decisions (truncLen budget, maxEE, pooling mode, ITS handling) live in microbiome/amplicon-processing.

# Primers come OFF before truncation - leftover primer bases corrupt the error model and look chimeric.
# -g = forward primer (5' on R1), -G = reverse primer (5' on R2); --discard-untrimmed drops primerless pairs.
cutadapt -g GTGYCAGCMGCCGCGGTAA -G GGACTACNVGGGTWTCTAAT --discard-untrimmed \
    -o trimmed_R1.fastq.gz -p trimmed_R2.fastq.gz sample_R1.fastq.gz sample_R2.fastq.gz

library(dada2)
# Run this block ONCE PER sequencing run. truncLen is a merge-overlap budget, not just a quality cut:
# truncLen_F + truncLen_R >= amplicon_length + ~12 (mergePairs minOverlap). See amplicon-processing.
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen = c(240, 160),
                     maxEE = c(2, 2), truncQ = 2, maxN = 0, rm.phix = TRUE,
                     compress = TRUE, multithread = TRUE)
errF <- learnErrors(filtFs, multithread = TRUE)   # fit THIS run only - never pool runs
errR <- learnErrors(filtRs, multithread = TRUE)
mergers <- mergePairs(dada(filtFs, err = errF, multithread = TRUE), filtFs,
                      dada(filtRs, err = errR, multithread = TRUE), filtRs)
seqtab_run <- makeSequenceTable(mergers)
# Combine per-run tables (exact-sequence string is the join key), THEN one chimera removal:
seqtab <- mergeSequenceTables(seqtab_run1, seqtab_run2)   # single run: skip; pass seqtab_run
seqtab_nochim <- removeBimeraDenovo(seqtab, method = 'consensus', multithread = TRUE)

A merge cliff (near-zero merged column) or a large READ fraction lost as chimeric is a stage-1/2 knob problem, not bad data - see microbiome/amplicon-processing. For low-biomass studies, sequence negative/positive controls and run decontam on the ASV table before downstream analysis (microbiome/amplicon-processing; metagenomics/contamination-controls).

Stage 4 to 6: Taxonomy, Tree, Diversity

Goal: Label ASVs against a region-matched reference, assemble a phyloseq object with a PLACED tree, and summarize diversity at a depth that retains samples.

Approach: Assign taxonomy with a classifier trained on the amplicon region (report genus for 16S, not species); filter host mitochondria/chloroplast features (universal 16S primers amplify them); build the phyloseq object and attach a SEPP/Greengenes2 placed tree rather than a de novo tree from short reads; rarefy ONLY into the diversity branch at a declared depth, report the dropped samples, and pair adonis2 with betadisper. Decisions live in microbiome/taxonomy-assignment and microbiome/diversity-analysis.

library(phyloseq); library(vegan)
# minBoot 50 = DADA2/RDP default for reads <=250 nt; ranks below it return NA, not a guess.
taxa <- assignTaxonomy(seqtab_nochim, 'silva_nr99_v138.1_train_set.fa.gz', minBoot = 50, multithread = TRUE)
ps <- phyloseq(otu_table(seqtab_nochim, taxa_are_rows = FALSE), tax_table(taxa),
               sample_data(metadata), phy_tree(placed_tree))   # placed_tree from SEPP/GG2, NOT de novo
ps <- subset_taxa(ps, is.na(Order) | Order != 'Chloroplast')      # drop host organelle 16S before diversity/DA
ps <- subset_taxa(ps, is.na(Family) | Family != 'Mitochondria')

# Diversity branch ONLY: rarefy to a DECLARED depth (not min(sample_sums)) and report who was dropped.
depth <- 10000   # choose on the alpha-rarefaction plateau; below this samples are dropped - report them
ps_rare <- rarefy_even_depth(ps, sample.size = depth, rngseed = 42, replace = FALSE)
adonis2(UniFrac(ps_rare, weighted = TRUE) ~ Group, data = data.frame(sample_data(ps_rare)), permutations = 999)
permutest(betadisper(UniFrac(ps_rare, weighted = TRUE), sample_data(ps_rare)$Group))   # location vs dispersion

SEPP placement is qiime fragment-insertion sepp; a de novo align-to-tree-mafft-fasttree is acceptable only when no reference package fits the marker, and unweighted UniFrac on it must be treated as suspect (microbiome/diversity-analysis).

Stage 7: Consensus Differential Abundance

Goal: Identify differentially abundant taxa as a confidence-graded consensus, not one tool's volcano plot.

Approach: On the UNrarefied counts (rarefying discards information DA needs), filter rare features, run ALDEx2 plus a second compositionally-aware tool (ANCOM-BC2 or LinDA), gate ALDEx2 on q AND effect size, and report the intersection as high-confidence. Tool mechanics live in microbiome/differential-abundance.

library(ALDEx2); library(ANCOMBC)
ps_filt <- filter_taxa(ps, function(x) sum(x > 0) >= 0.10 * nsamples(ps), TRUE)   # >=10% prevalence; declared knob
counts <- as.matrix(otu_table(ps_filt)); if (!taxa_are_rows(ps_filt)) counts <- t(counts)   # taxa in ROWS, integer counts
groups <- as.character(sample_data(ps_filt)$Group)

ax <- aldex(counts, groups, mc.samples = 128, test = 't', effect = TRUE, denom = 'all')
sig_aldex <- rownames(ax)[ax$we.eBH < 0.05 & abs(ax$effect) > 1]   # q AND |effect|>1 (between-group diff exceeds within-condition dispersion); NOT p alone

# Multi-run study: carry run as a batch covariate -> fix_formula = 'run_id + Group'. (ALDEx2 cannot
# take a covariate via aldex(); use aldex.clr() + aldex.glm() on a model matrix for the run-adjusted arm.)
ab <- ancombc2(data = ps_filt, fix_formula = 'Group', p_adj_method = 'BH',   # default is 'holm' - set BH deliberately
               prv_cut = 0.10, group = 'Group', struc_zero = TRUE, pseudo_sens = TRUE)$res
dcol <- grep('^diff_Group', names(ab), value = TRUE)[1]   # coefficient = variable+factor level, verbatim case (e.g. 'Grouptreated')
sig_ancombc <- ab$taxon[ab[[dcol]] & ab[[sub('^diff_', 'passed_ss_', dcol)]]]    # significant AND pseudo-count-robust

confident <- intersect(sig_aldex, sig_ancombc)   # high-confidence; union = exploratory; name both tools

Stage 8: Functional Prediction (optional)

Goal: Summarize predicted community functional POTENTIAL, framed as hypothesis-generating and gated on NSTI.

Approach: Run PICRUSt2 on the rep-seqs and ASV table, report the NSTI distribution and the read fraction dropped at --max_nsti 2, and restrict every claim to "potential" - never "activity". Decisions live in microbiome/functional-prediction; for MEASURED function use shotgun (metagenomics/functional-profiling).

# Predicts KO/MetaCyc POTENTIAL from who-is-there - never measured genes, never activity.
picrust2_pipeline.py -s asv_seqs.fna -i asv_table.biom -o picrust2_out -p 8 --max_nsti 2 --hsp_method mp
# Report mean/median NSTI and the ASV+read fraction dropped at NSTI>2 (marker_predicted_and_nsti.tsv.gz).

Per-Stage Failure Modes

Primers left on before truncation (stage 1)

Trigger: running filterAndTrim/learnErrors on reads still carrying primers. Mechanism: degenerate primer bases read as sequencing error, corrupting the error model and masquerading as chimeras. Symptom: wrong error fit, a large READ fraction lost as chimeric, inflated ASV count. Fix: cutadapt --discard-untrimmed first; order is primers -> filter -> learnErrors (microbiome/amplicon-processing).

Pooling runs into one error model (stage 2)

Trigger: concatenating multiple runs' FASTQs before learnErrors. Mechanism: one error model is fit to a mixture of run-specific error structures. Symptom: ASVs that vanish or appear when runs are split. Fix: per-run inference, then mergeSequenceTables, then one chimera removal; carry run as a batch covariate into DA.

Merge cliff from over-truncation (stage 2)

Trigger: truncLen_F + truncLen_R below amplicon length + ~12. Mechanism: denoised pairs no longer overlap enough to merge. Symptom: near-zero merged column, misread as "low diversity". Fix: compute the overlap budget first; keep length and loosen maxEE on the reverse for long amplicons.

De novo tree on short reads (stage 5)

Trigger: UniFrac/Faith PD on a MAFFT+FastTree tree from ~250 bp reads. Mechanism: short reads give an unstable topology and arbitrary midpoint root. Symptom: unweighted-UniFrac separation that vanishes under SEPP. Fix: use SEPP fragment-insertion or Greengenes2 placement (microbiome/diversity-analysis).

Sampling-depth sample massacre (stage 6)

Trigger: a rarefaction depth above some samples' totals (or min(sample_sums)). Mechanism: samples below the depth are silently dropped; the lost ones skew low-biomass. Symptom: fewer points in the PCoA than samples in the metadata. Fix: pick the depth from the alpha-rarefaction plateau, report the dropped-sample list, confirm at a nearby depth.

Rarefied table reused for DA (stage 6 -> 7)

Trigger: feeding ps_rare into the DA tools. Mechanism: rarefaction discards count information the compositional model needs. Symptom: underpowered or distorted DA. Fix: keep raw counts; rarefy only into the diversity branch; run DA on the unrarefied ps.