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.

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.

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.

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):

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.

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.