# Bio Applied Functional Annotation

> Profile MetaCyc pathways/genes from metagenomes with HUMAnN3, test differential abundance via MaAsLin2, detect AMR genes with AMRFinderPlus. Use for metagenomic functional annotation, pathway abundance, or AMR profiling.

- Skill: `pavel-kravchenko/bio-applied-functional-annotation` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-functional-annotation`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-functional-annotation/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-functional-annotation

---


# Functional Annotation of Metagenomes

## When to Use

- Quantifying microbial gene families and MetaCyc pathways from shotgun metagenomic reads (HUMAnN3 output: `_genefamilies.tsv`, `_pathabundance.tsv`, `_pathcoverage.tsv`).
- Comparing pathway/gene-family abundance between conditions (e.g. Healthy vs IBD vs CRC) with covariate-adjusted GLMs (MaAsLin2).
- Screening assembled contigs or reads for antimicrobial resistance (AMR) genes, virulence factors, or CAZymes.
- Normalizing HUMAnN abundance tables (RPK to CPM) or joining/splitting per-sample tables for cohort-level analysis.

## Version Compatibility

- HUMAnN ≥ 3.6 (MetaPhlAn ≥ 4.0, Bowtie2 ≥ 2.4, DIAMOND ≥ 2.0) — ChocoPhlAn + UniRef90 databases.
- MaAsLin2 ≥ 1.16 (Bioconductor/R ≥ 4.2).
- AMRFinderPlus ≥ 3.12 with CARD/NCBI reference database ≥ 2024.
- Python ≥ 3.10 with pandas ≥ 2.0, numpy ≥ 1.26, statsmodels ≥ 0.14, seaborn ≥ 0.13.

## Prerequisites

- Reads already quality-filtered and host/contaminant-decontaminated (see `bio-read-qc-quality-filtering`, `bio-metagenomics-kraken-classification` for decontamination).
- `humann`, `metaphlan`, `bowtie2`, `diamond` installed and on PATH; ChocoPhlAn nucleotide DB and UniRef90 DIAMOND DB downloaded via `humann_databases`.
- R with `Maaslin2` (Bioconductor) for differential abundance; `amrfinder` (+ `amrfinder_update`) for AMR screening.

## HUMAnN3 Workflow

**Goal:** turn paired-end decontaminated FASTQ into per-sample gene-family and pathway abundance tables.
**Approach:** merge paired reads into one file (HUMAnN treats input as a single unpaired stream), run the three-stage pipeline (MetaPhlAn4 taxonomic profiling → Bowtie2 pangenome alignment → DIAMOND translated search of unmapped reads against UniRef90), then renormalize RPK to CPM for cross-sample comparability.

```bash
cat decontam_1.fastq.gz decontam_2.fastq.gz > merged.fastq.gz

humann \
    --input merged.fastq.gz \
    --output humann3_out/ \
    --threads 16 \
    --metaphlan-options '--bowtie2db metaphlan4_db' \
    --nucleotide-database chocophlan_db/ \
    --protein-database uniref90_diamond/

# Normalize RPK -> CPM (copies per million) for cross-sample comparison
humann_renorm_table \
    --input humann3_out/sample_pathabundance.tsv \
    --output humann3_out/sample_pathabundance_cpm.tsv \
    --units cpm

# Join all samples into one cohort table, then split stratified/unstratified rows
humann_join_tables --input humann3_out/ \
    --output all_samples_pathabundance.tsv --file_name pathabundance_cpm
humann_split_stratified_table \
    --input all_samples_pathabundance.tsv --output humann3_stratified/
```

`_genefamilies.tsv` and `_pathabundance.tsv` are in RPK (reads per kilobase, corrects for gene length); after CPM renormalization values are proportional to the fraction of reads mapping to that feature. Each pathway/gene family is reported **unstratified** (cohort total) and **stratified** by contributing species (`PathwayID|Species`). `_pathcoverage.tsv` gives the fraction of pathway reactions detected (capped at 1) — high coverage with low abundance means the pathway's genes are present but lowly expressed/covered by few reads.

## Visualizing and Loading Pathway Abundance

**Goal:** load a joined HUMAnN3 pathway table and visualize cohort-level patterns.
**Approach:** log2-transform CPM values (features span orders of magnitude) and heatmap by sample/condition; verify group differences with a per-pathway bar plot.

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns


def load_pathway_table(path: str) -> pd.DataFrame:
    """Load a HUMAnN3 joined/unstratified pathway abundance TSV (features x samples)."""
    df = pd.read_csv(path, sep='\t', header=0, index_col=0, comment='#')
    return df


def plot_pathway_heatmap(path_df: pd.DataFrame, ax=None):
    """Log2(CPM+1) heatmap of pathway abundance across samples."""
    log2_df = np.log2(path_df + 1)
    short_labels = [p.split(':')[1].strip()[:35] if ':' in p else p[:35] for p in path_df.index]
    ax = sns.heatmap(log2_df, ax=ax, cmap='YlOrRd', yticklabels=short_labels,
                      cbar_kws={'label': 'log2(CPM+1)'}, linewidths=0.3)
    ax.set_title('HUMAnN3 pathway abundances')
    return ax


# Demo with simulated data shaped like real HUMAnN3 output
rng = np.random.default_rng(42)
pathways = [
    'GLYCOLYSIS-E-D: Glycolysis (Embden-Meyerhof-Parnas)',
    'PYRUVATE-FERMENTATION-PWY: Pyruvate fermentation to acetate and lactate',
    'FASYN-ELONG-PWY: Fatty acid elongation, saturated',
    'BUTYRATE-BIOSYNTHESIS-I: Butyrate biosynthesis I',
    'PROPIONATE-BIOSYNTHESIS-I: Propionate biosynthesis I',
    'FOLSYN-PWY: Folate biosynthesis',
]
n_samples = 12
conditions = ['Healthy'] * 4 + ['IBD'] * 4 + ['CRC'] * 4
base = rng.lognormal(mean=6, sigma=2, size=(len(pathways), n_samples))

ibd_idx = [i for i, c in enumerate(conditions) if c != 'Healthy']
but_idx = next(i for i, p in enumerate(pathways) if 'Butyrate' in p)
base[but_idx, ibd_idx] *= 0.3   # butyrate biosynthesis reduced in IBD/CRC
base[0, ibd_idx] *= 1.5          # glycolysis up

path_df = pd.DataFrame(base, index=pathways,
                        columns=[f'{c[:3]}_{i % 4 + 1}' for i, c in enumerate(conditions)])

fig, ax = plt.subplots(figsize=(10, 5))
plot_pathway_heatmap(path_df, ax=ax)
plt.tight_layout()
plt.show()
```

## Differential Pathway Analysis with MaAsLin2

MaAsLin2 fits a GLM (with optional random effects) per feature against metadata variables, handling continuous/categorical covariates, repeated measures, and compositionality via TSS/CLR normalization.

```r
library(Maaslin2)

pathway_table <- read.table("humann3_stratified/all_pathabundance_unstratified.tsv",
                             sep = "\t", header = TRUE, row.names = 1, comment.char = "#")
metadata <- read.table("metadata.tsv", sep = "\t", header = TRUE, row.names = 1)

fit <- Maaslin2(
  input_data = t(pathway_table), input_metadata = metadata,
  output = "maaslin2_out/",
  fixed_effects = c("condition"), random_effects = c("batch"),
  transform = "LOG", normalization = "TSS",
  min_prevalence = 0.1, min_abundance = 0.0001
)
# Output: maaslin2_out/significant_results.tsv with columns feature, metadata, coef, pval, qval
```

For a quick Python-side sanity check of an already-fit result table (log2FC + BH-corrected q-values):

```python
from scipy import stats
from statsmodels.stats.multitest import multipletests


def call_significant_pathways(log2fc: np.ndarray, se: np.ndarray, df: int = 6,
                               alpha: float = 0.05) -> pd.DataFrame:
    """BH-adjust per-pathway t-test p-values from effect size + standard error."""
    t_stats = log2fc / se
    pvals = stats.t.sf(np.abs(t_stats), df=df) * 2  # two-sided
    _, qvals, _, _ = multipletests(pvals, method='BH')
    return pd.DataFrame({'log2FC': log2fc, 'pval': pvals, 'qval': qvals})


rng = np.random.default_rng(55)
n_pw = 8
log2fc = np.concatenate([rng.normal(-1.8, 0.5, 2), rng.normal(1.5, 0.6, 2),
                          rng.normal(0, 0.4, n_pw - 4)])
se = np.abs(rng.normal(0.3, 0.1, n_pw)) + 0.1
diff_df = call_significant_pathways(log2fc, se)
print(diff_df[diff_df['qval'] < 0.05].to_string())
```

## AMR Gene Detection

Gene family abundances (or assembled contigs) can also be screened for antimicrobial resistance genes (CARD/ResFinder/ARG-ANNOT), virulence factors (VFDB), or CAZymes.

```bash
amrfinder --nucleotide megahit_assembly/final.contigs.fa \
    --output amr_report.txt --threads 8 --plus
```

CARD resistance categories: **Intrinsic** (naturally present, e.g. AmpC), **Acquired** (HGT-transferred, e.g. blaCTX-M, mcr-1), **Mutational** (point mutations, e.g. gyrA quinolone resistance). Report AMR abundance as RPKM or presence/absence per sample, and check prevalence (fraction of samples with RPKM above a detection threshold) before calling a gene "present" in a cohort.

## Pitfalls

- **Merging paired reads**: HUMAnN takes one unpaired input stream — concatenate R1+R2 rather than running twice, or you double-count/miss cross-pairing information.
- **RPK vs CPM**: never compare raw RPK across samples with different sequencing depth — always renormalize with `humann_renorm_table --units cpm` first.
- **Stratified vs unstratified rows**: aggregate stats (MaAsLin2, heatmaps) should use the unstratified table; mixing both inflates feature counts.
- **Coverage != abundance**: `_pathcoverage.tsv` coverage is capped at 1 and reflects reaction completeness, not read depth — don't treat it as an abundance proxy.
- **Coordinate systems**: BED is 0-based half-open; VCF/GFF are 1-based inclusive — mixing them when annotating AMR hits back to contigs causes off-by-one errors.
- **Multiple testing**: thousands of pathways/gene families require FDR correction (Benjamini-Hochberg), not raw p-values.

## See Also

- `bio-metagenomics-functional-profiling` — alternative functional profiling tools and workflows.
- `bio-metagenomics-amr-detection` — dedicated AMR gene calling and interpretation.
- `bio-metagenomics-abundance-estimation` — taxonomic/functional abundance normalization.
- `bio-differential-expression-de-results` — general patterns for interpreting differential-abundance output tables.