Step 15 — PICRUSt2 (Functional Prediction)

Script: scripts/mbx_picrust_run.sh

Companion files in this folder: - 15_picrust2.html — same content with copy buttons. - 15_picrust2.pptx — slide deck for the talk.

Why this step matters

Every step up to now has answered "what's there" questions. Step 15 answers a "what is it doing" question — at least, predicted from the phylogenetic placement of each ASV.

16S sequencing tells us identity but nothing direct about function. PICRUSt2 (Douglas et al. 2020) fills that gap by exploiting the fact that closely-related organisms tend to share gene families. The algorithm places every ASV on a reference phylogenetic tree, looks up the known gene content of its nearest reference genomes, and infers the expected gene-family abundance per sample. The output is four functional abundance tables:

• KO (KEGG Orthology) — individual gene families.
• EC (Enzyme Commission) — enzymatic activities.
• COG (Clusters of Orthologous Groups) — broad functional categories.
• MetaCyc pathways — high-level metabolic pathways.

Step 15's job is to run PICRUSt2 reliably (it lives in its own conda env — see "fallbacks"), flag samples whose ASVs are too distantly related to the reference to trust the prediction, run differential abundance across treatments on each database, and produce a self-contained HTML report so the user can read the functional result alongside the taxonomic ones.

What the script does in one sentence

It runs PICRUSt2 in a pre-baked conda env to predict KO/EC/COG/MetaCyc abundances per sample, computes per-sample NSTI (the phylogenetic distance from each ASV to its nearest reference genome, the official trustworthiness metric), runs Kruskal-Wallis + Dunn differential abundance per categorical variable, optionally correlates genera × pathways with the Step 14 results, and writes a self-contained HTML report.

The algorithm, step by step

1. Gate on Step 11

First, same READY_FOR_DIVERSITY gate. PICRUSt2 uses raw counts (non-rarefied) like ANCOMBC2 — it needs absolute count information for its NSTI weighting.

2. Use the pre-baked mbx_picrust2 conda env

Then the script confirms a conda env named mbx_picrust2 exists. If the user is running inside the Docker image, the env is baked in at build time (see Phase 3.3 in CHANGELOG.md) — picrust2=2.6.3 from bioconda. If running on a fresh host, the script creates the env on first use (~2 GB download). All picrust2_*.py commands are invoked via conda run -n mbx_picrust2, so the script never modifies the active QIIME2 env.

3. Export the rep-seqs and feature table to PICRUSt2-friendly formats

Next the script converts the Step-4 representative_sequences.qza to a plain FASTA and the Step-7 filtered feature table to a BIOM/TSV pair. Those are PICRUSt2's required inputs.

4. Run place_seqs.py

Now PICRUSt2 places every ASV on its reference 20,000-genome tree using epa-ng (placement) + gappa (re-formatting):

place_seqs.py \
  -s rep_seqs.fasta \
  -o placed_seqs.tre \
  -p $THREADS

This is the most computationally expensive sub-step (a few minutes per 1,000 ASVs).

5. Compute NSTI per ASV

Then PICRUSt2 computes the Nearest Sequenced Taxon Index for each ASV — the phylogenetic distance from the ASV to its closest reference genome on the tree. A low NSTI (≤ 0.15) means the prediction is well-grounded in a real reference; a high NSTI (> 1.0) means the ASV is so distantly related to anything sequenced that the predicted gene content is essentially a guess.

The script computes the weighted mean NSTI per sample (read-count- weighted, not just ASV-mean) and **flags any sample with weighted NSTI

1.0** in the final report. Those samples are excluded from the per- variable DA analysis by default.

6. Predict gene-family abundance — three databases

For each of KO, EC, COG:

metagenome_pipeline.py \
  -i feature_table.biom \
  -m marker_nsti_predicted.tsv.gz \
  -f <KO|EC|COG>_predicted.tsv.gz \
  --strat_out          \
  -o <KO|EC|COG>_metagenome/

Output is per-sample predicted abundance of every KO/EC/COG term.

7. Predict MetaCyc pathway abundance

Then the script runs pathway_pipeline.py to roll the EC predictions into MetaCyc pathway abundances. MetaCyc is the standard human-readable pathway database — much friendlier than raw KO codes for biology-focused readers.

8. Per-variable Kruskal-Wallis + Dunn on each database

For every categorical metadata variable (auto-detected, same as Steps 9/10) the script runs:

• Kruskal-Wallis on each KO / EC / COG / MetaCyc-pathway abundance.
• Benjamini-Hochberg across the gene families.
• Dunn pairwise on KW hits.
• Compact letter displays on significant pathways.

Outputs: DA_KO_KW_<Variable>.xlsx, DA_EC_KW_<Variable>.xlsx, DA_COG_KW_<Variable>.xlsx, DA_metacyc_pathways_KW_<Variable>.xlsx, each with the per-term p, p_adj, group means, and Dunn pairs for significant terms.

9. Cross-reference with Step 14 (when available)

If Step 14 ran, the script computes Spearman correlations between every (significant differentially-abundant genus from Step 14) and every (significant differentially-abundant MetaCyc pathway). The result is correlation_DAgenera_x_DApathways_<Variable>.{tsv,png,pdf} — a heatmap that lets the reader say "the increase in pathway X is driven by the increase in genus Y".

10. Plots

• Stacked-bar top-20 MetaCyc pathways per variable.
• Heatmap of differentially-abundant pathways with log-fold-change + significance stars.
• Functional PCoA based on Bray-Curtis of MetaCyc abundance — the functional analogue of Step 13's taxonomic PCoA.
• Genus × pathway correlation heatmap (when Step 14 is available).

PNG + SVG always; PDF on --publication-figures.

11. Self-contained HTML report

Finally the script writes picrust2_report.html, embedding every plot as inline <img src="data:image/..."> so the user can email the single file. Per-sample NSTI flagging is shown at the top so the reader can see immediately which samples had unreliable predictions.

Default parameters and why they are what they are

When and why we fall back to defaults

What the output file looks like

Summary_picrust2_NSTI.xlsx:

DA_metacyc_pathways_KW_<Variable>.xlsx (per-pathway KW results):

Plus the cross-referenced correlation_DAgenera_x_DApathways_<Variable>.tsv when Step 14 is present, and the self-contained HTML report.

Takeaway

Step 15 turns the taxonomic 16S story into a functional story — phylogenetic placement + nearest-relative gene lookup + per-sample NSTI trust scoring. The script keeps PICRUSt2 in its own conda env so it doesn't poison the QIIME2 env, and bakes that env into the Docker image so the pipeline is fully offline-capable.

Sources

• The script: mbXPro/scripts/mbx_picrust_run.sh
• PICRUSt2: Douglas et al. (2020), PICRUSt2 for prediction of metagenome functions, Nature Biotechnology 38:685–688.
• NSTI: Langille et al. (2013), Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences, Nature Biotechnology 31:814.
• MetaCyc: Caspi et al. (2020), The MetaCyc database of metabolic pathways and enzymes, NAR 48:D445.
• KEGG: Kanehisa et al. (2017), KEGG: new perspectives on genomes, pathways, diseases and drugs, NAR 45:D353.