Functional Profiling: HUMAnN and PICRUSt2

Taxonomic profiling answers who is in a community; functional profiling asks what those organisms can do by quantifying gene families and metabolic pathways. On shotgun metagenomes you measure this directly with HUMAnN, which maps reads to protein families. On 16S amplicon data you have no genes to map, so PICRUSt2 predicts likely gene content from the marker gene instead.

Running HUMAnN

HUMAnN uses a tiered search: it first profiles species with MetaPhlAn, aligns reads with Bowtie2 against the pangenomes of detected species (ChocoPhlAn), then runs a translated DIAMOND search of the leftovers against UniRef90. It reports UniRef90 gene family abundances in reads-per-kilobase (length-normalized but not depth-normalized) and MetaCyc pathway abundances reconstructed from the underlying reactions.

# Input reads should already be quality-trimmed and host-depleted.
# HUMAnN ignores read pairing, so concatenate mates into one file first:
#   cat sample_R1.fastq.gz sample_R2.fastq.gz > sample.fastq.gz
# Databases: humann_databases --download chocophlan full $DBDIR
#            humann_databases --download uniref uniref90_diamond $DBDIR
humann --input sample.fastq.gz \
       --output humann_out \
       --output-basename sample \
       --threads 8

Output files will be written to: /home/user/humann_out
Running metaphlan ........
Total species selected from prescreen: 12
Running diamond ........
Computing gene families ...
Computing pathways abundance and coverage ...
Output files created:
/home/user/humann_out/sample_genefamilies.tsv
/home/user/humann_out/sample_pathabundance.tsv
/home/user/humann_out/sample_pathcoverage.tsv

Regroup and renormalize

# RPK depends on sequencing depth, so renormalize before comparing samples.
# Regroup UniRef90 gene families into KEGG Orthology (needs utility_mapping DB):
#   humann_databases --download utility_mapping full $DBDIR
humann_regroup_table \
  --input humann_out/sample_genefamilies.tsv \
  --groups uniref90_ko \
  --output sample_ko.tsv

# Renormalize RPK to copies per million (use --units relab for relative abundance):
humann_renorm_table \
  --input sample_ko.tsv \
  --units cpm \
  --output sample_ko_cpm.tsv

# Merge many single-sample tables into one matrix for downstream stats:
humann_join_tables --input humann_out --output all_genefamilies.tsv --file_name genefamilies

PICRUSt2 from 16S

PICRUSt2 places your ASV sequences into a reference tree of genomes with known gene content, then uses hidden-state prediction to estimate the KEGG Orthologs, EC numbers, and MetaCyc pathways each ASV is expected to carry, weighted by the feature table counts. It also reports a Nearest Sequenced Taxon Index (NSTI) per ASV: the phylogenetic distance to its closest reference genome, which tells you how trustworthy each prediction is.

# Inputs: representative ASV sequences (FASTA) + feature table (BIOM or TSV).
# ASVs with NSTI above --max_nsti (default 2.0) are dropped as too distant.
picrust2_pipeline.py \
  -s asv_seqs.fna \
  -i asv_table.biom \
  -o picrust2_out \
  -p 8

# Attach human-readable labels to the predicted KO abundances:
add_descriptions.py \
  -i picrust2_out/KO_metagenome_out/pred_metagenome_unstrat.tsv.gz \
  -m KO \
  -o picrust2_out/KO_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz