mbcausal

Disease agnostic pipeline for causal inference from longitudinal gut microbiome profiles: profile shotgun or long read metagenomes, build QC'd per cohort abundance tables, CLR transform them, and run interpretable representation learning plus lagged causal discovery over clinical outcomes. Developed on the iHMP/IBDMDB inflammatory bowel disease cohort.

Full documentation is in the wiki. This file is the quick reference: install it, run the four stages, find the output. The wiki goes further on each stage, on the simulation study, and on reproducing the paper.

Installation	env and editable install
Configuration	paths and run settings
Profiling	reads to abundance tables
Preprocessing	cohort tables, QC and CLR
Discovery	factors, screening and lagged effects
Simulation	cohorts where the answer is known
Reproduction	the commands behind the paper

Setup

make env                # create the conda env (conda env create -f environment.yml)
conda activate mbcausal
make install            # editable install (pip install -e .)

Pipeline

The examples below write to $REPORTS and $RESULTS. Set them to wherever your profiler output and cohort tables should live, matching reports_dir and preprocess_dir in your mbcausal.yaml:

export REPORTS=/path/to/reports RESULTS=/path/to/preprocess

Stages chain through standard artifact files. Every subcommand takes -i/-o directories and writes a run log to <out>/logs/<cmd>_<host>_<timestamp>.log.

profile ──▶ build-cohort ──▶ preprocess ──▶ discover
(fastq)     (per cohort       (CLR feature   (lagged factor
            tables + QC)      matrix)        causal discovery)

1. profile: fastq to per sample profiler reports

# short reads (Illumina) with Kraken2+Bracken
mbcausal profile -i <fastq_dir> -o $REPORTS/ihmp \
  --profiler kraken2bracken --platform illumina --max-workers 8

# long reads (ONT) with ganon2 (or metaphlan4)
mbcausal profile -i <fastq_dir> -o $REPORTS/<cohort> \
  --profiler ganon2 --platform ont --max-workers 4

Profilers: kraken2bracken, metaphlan4, ganon2. --platform ont sets long read fastq suffixes and enables MetaPhlAn4 --long_reads. Scale throughput with --max-workers rather than --threads when the host is shared.

2. build-cohort: reports to model input tables, with QC

Merges each tool's reports into <cohort>-<tool>/{abundance.csv, sample_table.csv, qc.csv}. sample_table.csv carries a shared qc_pass flag: it drops technical replicates and samples sequenced too shallowly, and drops the same samples across all of a cohort's tools, so comparisons between profilers are like for like.

# iHMP: qc_pass depth comes from bracken's read counts
mbcausal build-cohort -i $REPORTS/ihmp -o $RESULTS \
  --cohort ihmp --depth-from bracken --min-depth 150000000

# a cohort whose depth comes from the metadata `bases` column instead
mbcausal build-cohort -i $REPORTS/<cohort> -o $RESULTS \
  --cohort <name> --min-depth 100000000

(build-abundance -i <reports> --profiler X is the step below this one, merging a single tool into an abundance.csv with no metadata and no QC.)

How qc_pass is decided, and why preprocess needs its own output subdirectory: Preprocessing.

3. preprocess: cohort tables to a CLR feature matrix

Applies the prevalence filter and CLR to a built cohort, writing the feature_matrix.csv and the aligned sample_table.csv that discover reads. Filters on QC by default; --no-qc keeps every sample.

mbcausal preprocess -i $RESULTS/ihmp-bracken \
  -o $RESULTS/ihmp-bracken/clr --outcome fecalcal

4. discover: causal discovery and robustness

Fits interpretable NMF factors and estimates each factor's lagged within patient effect on a time varying outcome, as outcome(t+lag) ~ factor(t) | baseline outcome + confounders with standard errors clustered by patient. It then runs the robustness suite: seed stability, specification curve, permutation null, Mundlak decomposition.

mbcausal discover -i $RESULTS/ihmp-bracken/clr \
  --outcome fecalcal --unit patient_id --time week \
  --confounders antibiotics --confounder-lags 1 --lag 1

Writes the suite's tables to <input>/discovery/, overridable with -o. Common knobs: --n-factors, --subgroup diagnosis=CD,UC, --max-gap, --n-seeds, --n-perm, --steps. See mbcausal discover -h for the full set.

Simulation

mbcausal sim runs the simulation study: cohorts where the planted answer is known, put through the identical workflow.

mbcausal sim list                          # the available experiments
mbcausal sim generate --root <dir>         # write simulated cohorts, four stages each
mbcausal sim run rtm_demo --n-rep 300      # one experiment; flags pass through

Generators live in mbcausal.sim.generate and do not import the rest of the package, so a defect in one cannot cancel itself out in the other.

The defaults are calibrated to iHMP and nothing else. The generator is not specific to a cohort, but every parameter it ships with was fitted to the iHMP MetaPhlAn4 tables and to calprotectin as the outcome. Using it to argue about another cohort means refitting those values first.

What gets planted and how a run is scored: Simulation.

Using the tables in model code

Downstream code loads a built cohort through a single entry point, which applies qc_pass by default, so no model code can train on samples that failed QC by accident:

from mbcausal.data import load_cohort

cohort = load_cohort("preprocess/ihmp-bracken")   # QC clean by default
# cohort.X = abundance (samples x species), cohort.meta = metadata
# load_cohort(dir, qc=False) returns the complete, unfiltered record

Configuration

Settings (paths, threads and so on) resolve in this order: per run argument, then env var (MBCAUSAL_*), then mbcausal.yaml, then the default.

• Copy mbcausal.example.yaml to mbcausal.yaml and edit it. The real file is gitignored.
• Relative paths anchor to the project root; absolute paths are used as they are.
• The defaults are relative placeholders, so a fresh clone writes under the project. Point reports_dir, preprocess_dir and the database paths at real locations before running.

Outputs

• reports_dir: raw profiler output, as <tool>-reports/<sample>/ and <tool>-logs/<sample>.{log,err,time}, where .time holds wall clock and peak RAM.
• preprocess_dir: one <cohort>-<tool>/ per cohort and tool, holding abundance.csv with its abundance.meta.json provenance, sample_table.csv with qc_pass, and qc.csv. preprocess adds feature_matrix.csv; discover writes a discovery/ table set.
• Every subcommand writes logs/<cmd>_<host>_<timestamp>.log.
• Existing outputs are not overwritten unless --overwrite.

Develop

make test        # pytest
make lint        # ruff
make typecheck   # mypy
make check       # all three

More

Everything above in more depth, plus the exact commands behind the paper's numbers, is in the wiki.

Status

The upstream stages (profiling, cohort tables, CLR) and the causal layer (NMF factors and the lagged_factor discovery and robustness model) are implemented and tested. GraphSAGE is implemented as an auxiliary encoder and is never used as a causal exposure, since its dimensions carry no named taxa. Autoencoder, additional causal methods, prediction and evaluation are stubs.