Profiling.md

Profiling

Classify reads and estimate species abundances, producing a per sample abundance table (relative abundance, host removed). Four profilers are supported and swap through one flag.

Profiler	Notes
metaphlan4	marker gene, SGB database. Primary for the current analysis.
kraken2bracken	k mer classification with Bracken reestimation
ganon2	HIBF index, abundances corrected for genome size
MetaPhlAn3	not run here; the iHMP distribution ships precomputed profiles that we parse

Database locations come from mbcausal.yaml. See Configuration.

Profile a cohort

mbcausal profile -i <fastq_dir> -o $REPORTS/ihmp \
  --profiler metaphlan4 --platform illumina --max-workers 8

--platform sets the fastq suffix and pairing, plus per profiler presets. Use ont for long reads, which are single end.

Scale with --max-workers, not --threads. Roughly 5 to 6 threads per sample is the useful maximum for both kraken2 and metaphlan. Note that kraken2 shares one memory mapped database across workers, about 92 GB in total, whereas metaphlan and bowtie2 load their own index per process, about 20 GB times the number of workers.

Output layout

<reports_dir>/
  <label>-reports/<sample>/...          one directory per sample
  <label>-logs/<sample>.{log,err,time}  .time records wall clock and peak RAM

Existing outputs are not overwritten. Pass --overwrite to redo a sample.

Reports to a cohort table

mbcausal build-cohort -i $REPORTS/ihmp -o $RESULTS \
  --cohort ihmp --depth-from bracken --min-depth 150000000

Writes <cohort>-<tool>/ holding abundance.csv, sample_table.csv with a qc_pass column, and qc.csv. One keep list is computed per cohort from a single depth source and applied to every profiler, so profilers are compared on identical samples.

The sidecar abundance.meta.json records the profiler, the database, the parameters and the tool version.

Notes

• parse drops Homo sapiens and renormalizes to microbial relative abundance, so rows sum to 1. The unclassified fraction is discarded.
• Downstream code should load through load_cohort, which applies qc_pass by default so no model can train on QC failing samples.