Bio Microbiome Functional Prediction

Version Compatibility

Reference examples tested with: PICRUSt2 2.5+, pandas 2.2+.

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

• CLI: <tool> --version then <tool> --help to confirm flags
• Python: pip show <package> then help(module.function) to check signatures

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

The bundled REFERENCE is the version that matters. PICRUSt2 ships its own reference tree, alignment, and per-genome trait tables (16S/KO/EC/Pfam/COG/TIGRFAM) for ~20,000 genomes - there is no separate multi-GB download as in HUMAnN. That fixed reference is the entire ceiling on accuracy, and NSTI is computed against THAT reference. The reference is versioned with the PICRUSt2 release; record the release so a prediction is reproducible.

Functional Prediction with PICRUSt2

"Predict the functional pathways from my 16S data" -> Place ASVs on a reference genome tree and report the gene content of their nearest sequenced relatives - because PICRUSt2 never sees a functional gene from the sample, it infers POTENTIAL from who-is-there.

• CLI: picrust2_pipeline.py -s asv_seqs.fna -i asv_table.biom -o picrust2_out -p 8 --max_nsti 2

Scope: PREDICTED gene-content potential from amplicon ASVs. MEASURED shotgun function (HUMAnN) -> metagenomics/functional-profiling. ASV table + rep-seqs come from amplicon-processing + taxonomy-assignment. The QIIME2 q2-picrust2 path -> qiime2-workflow. DA of predicted tables -> differential-abundance (compositional, same theory as metagenomics/abundance-estimation). Reading KO/pathway lists -> pathway-analysis/go-enrichment.

The Single Most Important Modern Insight -- Predicted Function Is Taxonomy Re-Encoded, Not a Measurement

PICRUSt2 never sequences a functional gene from the sample. It places each ASV on a tree of ~20,000 reference genomes, asks "what genes do this ASV's nearest sequenced relatives carry?", and reports that guess as the community's function. Every output number is the gene content of OTHER organisms, weighted by how much of the 16S resembles them. Three corollaries each common misuse violates:

1. Potential, never activity. The honest claim is "increased butyrate-production capacity," never "increased butyrate production / upregulated / more active." The inference chain is 16S abundance -> who-is-there -> relatives' genomes -> predicted gene presence -> predicted copy number -> potential. Transcription and flux are further measurement layers away.
2. Accuracy IS reference coverage. Predicted-vs-shotgun Spearman tops out around 0.8 in the densely-referenced human gut (Douglas 2020); in soil, marine, sediment, and novel environments the references are sparse, NSTI rises, and the prediction collapses toward a restatement of taxonomy.
3. Circularity. Predicted function is a DETERMINISTIC function of the ASV table - same input, same KO table, every time. So a "functional difference" between groups is the TAXONOMIC difference projected through a fixed lookup, not independent corroboration. Predicted-function analysis is hypothesis-generating; a functional conclusion needs shotgun or metatranscriptomics.

Organize the analysis around defending these three (report NSTI, forbid activity verbs, do not double-count taxonomy as a second finding), not around listing flags.

The Predicted < Measured Ladder

State this explicitly when reporting. Each rung is a real measurement layer above the one below:

predicted potential (PICRUSt2, this skill) < measured gene carriage (HUMAnN / shotgun, metagenomics/functional-profiling) < measured transcription (metatranscriptomics) < measured flux (fluxomics/metabolomics).

PICRUSt2 is the bottom rung. A PICRUSt2 path_abun_unstrat table looks identical to a HUMAnN pathabundance table - same MetaCyc IDs, same shape, same MinPath logic - but HUMAnN counts reads that actually aligned to genes in the sample, while PICRUSt2 reports relatives' genomes. Never merge or compare the two as interchangeable.

How PICRUSt2 Builds the Table (the five-stage interpolation)

Every stage is an inference layered on the previous, so error compounds (Douglas 2020 Nat Biotechnol 38:685; original Langille 2013 Nat Biotechnol 31:814). PICRUSt2's headline advance over PICRUSt1 is that it places arbitrary de-novo ASVs (not closed-reference Greengenes OTUs):

1. Placement. Align ASVs to the reference alignment (hmmalign), place into the reference tree with EPA-ng (Barbera 2019 Syst Biol 68:365), resolve placements with gappa (Czech 2020 Bioinformatics 36:3263). Default --placement_tool epa-ng (alternative sepp).
2. Hidden-state prediction (HSP). For each placed ASV - whose genome was never sequenced - interpolate the copy number of every gene family from neighboring reference genomes using castor (Louca & Doebeli 2018 Bioinformatics 34:1053). Default --hsp_method mp (maximum parsimony; pic is faster but mp is the recommended default).
3. 16S copy-number normalization. Divide each ASV's abundance by its PREDICTED 16S copy number (a genome with seven 16S copies is over-counted 7x in an amplicon survey). The correction is itself an HSP output, not a measurement.
4. Metagenome inference. Per-sample gene abundance = sum over ASVs of (ASV abundance / 16S copies) x (predicted gene copy number). Emits KO/EC/Pfam tables.
5. Pathway inference. Map genes to MetaCyc reactions and call pathways with MinPath (Ye & Doak 2009 PLoS Comput Biol 5:e1000465). Pathway COVERAGE (--coverage, presence confidence) and pathway ABUNDANCE are different questions - do not conflate.

Tool Taxonomy

FAPROTAX is conceptually different: a curated taxon-to-function lookup answering "is this community nitrifying?", saying nothing about uncharacterized taxa (they map to nothing). Use FAPROTAX for biogeochemical-cycle questions, PICRUSt2 for a predicted KO/pathway profile - different questions.

Decision Tree by Scenario

Run the Pipeline

# Full pipeline: place -> HSP -> 16S-normalize -> metagenome -> MetaCyc pathways
picrust2_pipeline.py \
    -s asv_seqs.fna \              # representative ASV sequences (FASTA)
    -i asv_table.biom \            # ASV abundance table (BIOM or TSV; samples as columns)
    -o picrust2_out \
    -p 8 \
    --hsp_method mp \              # maximum parsimony (recommended default); pic is faster but not recommended
    --max_nsti 2 \                 # ASVs ABOVE 2 are DROPPED before inference; report how many (see below)
    --verbose
# Key outputs (gzipped):
#   KO_metagenome_out/pred_metagenome_unstrat.tsv.gz   KEGG ortholog abundances
#   EC_metagenome_out/pred_metagenome_unstrat.tsv.gz   EC-number abundances
#   pathways_out/path_abun_unstrat.tsv.gz              MetaCyc pathway abundances
#   marker_predicted_and_nsti.tsv.gz                   per-ASV 16S copies + metadata_NSTI (the quality file)

--stratified additionally emits per-ASV contribution tables (large, much slower); --per_sequence_contrib is only meaningful with it. --coverage adds pathway coverage (a different question from abundance). The unrolled per-step scripts are place_seqs.py -> hsp.py -i {16S,KO,EC} -m mp [-n] -> metagenome_pipeline.py --max_nsti 2 -> pathway_pipeline.py; --max_nsti filtering and 16S normalization happen in metagenome_pipeline.py. add_descriptions.py -m METACYC attaches human-readable names.

Report NSTI (mandatory)

Goal: Quantify how much of the prediction is extrapolation and how much of the sampled community the NSTI gate discarded, so the result is interpretable.

Approach: Read marker_predicted_and_nsti.tsv.gz, summarize the metadata_NSTI distribution, and report the number of ASVs AND the fraction of READS dropped at --max_nsti 2 (a study that loses 40% of reads predicted function for a different community than it sampled).

import pandas as pd

nsti = pd.read_csv('picrust2_out/marker_predicted_and_nsti.tsv.gz', sep='\t')   # cols: sequence, metadata_NSTI
asv_counts = pd.read_csv('asv_table.tsv', sep='\t', index_col=0)                # ASVs x samples
nsti = nsti.set_index('sequence')
reads_per_asv = asv_counts.sum(axis=1)

max_nsti = 2.0   # PICRUSt2 default; ASVs above this are dropped before metagenome inference
dropped = nsti.index[nsti['metadata_NSTI'] > max_nsti]
reads_dropped_frac = reads_per_asv.reindex(dropped).sum() / reads_per_asv.sum()
print(f'mean NSTI {nsti.metadata_NSTI.mean():.3f}  median {nsti.metadata_NSTI.median():.3f}')
print(f'ASVs dropped at NSTI>{max_nsti}: {len(dropped)}/{len(nsti)}  reads dropped: {reads_dropped_frac:.1%}')

Per-Method Failure Modes

Claiming activity or expression

Trigger: reporting "increased butyrate production" / "upregulated" / "more metabolically active." Mechanism: PICRUSt2 measured no genes and no transcripts - only inferred gene presence from relatives' genomes. Symptom: a results sentence with an activity verb on a predicted pathway. Fix: restrict every claim to "potential" / "predicted capacity"; for activity, cite metatranscriptomics, not this skill.

NSTI ignored or under-reported

Trigger: accepting the default --max_nsti 2 filter without reporting the distribution or the dropped fraction. Mechanism: the filter silently deletes the most novel/under-referenced ASVs - exactly the organisms an environmental study cares about. Symptom: a predicted-function result with no NSTI numbers in the methods. Fix: report mean/median NSTI, the distribution, and the ASV AND read fraction dropped (the helper above); treat high mean NSTI as a red flag that the result is mostly extrapolation.

Wrong environment (reference coverage too sparse)

Trigger: running PICRUSt2 on soil/marine/sediment/plant/novel hosts and reporting fine-grained KO differences. Mechanism: the ~20k-genome reference tree is gut/host-biased; sparse references mean high NSTI and predictions interpolated from distant relatives. Symptom: high mean NSTI yet confident KO/pathway tables. Fix: report the environment and NSTI; prefer FAPROTAX for the broad biogeochemical question, or do real shotgun. "Relatively better than other predictors" (per the paper) is not "trustworthy in absolute terms."