Skip to content

Reproduction.md

Jump to bottom
Wenyu (Eddy) Huang edited this page Aug 5, 2026 · 1 revision

Reproduction

The commands behind the paper's numbers. Each was recovered from the logs/ directory of the run that produced the result, so it is what ran rather than what was meant to run.

Every run log opens with the commit it ran at. Check there rather than trusting this page if the two ever disagree.

Set your paths first

The commands below use these. Point them wherever your data lives.

export RAW=/path/to/HMP2/raw                 # downloaded fastq
export REPORTS=/path/to/reports              # profiler output
export PREPROCESS=/path/to/preprocess        # cohort tables and CLR matrices
export SIMDATA=/path/to/data/sim             # generated simulation cohorts
export RUNS=/path/to/discover                # discovery output

Real data

0. Get the data

Raw shotgun metagenomes for the iHMP IBD cohort come from IBDMDB (https://ibdmdb.org/results), served over Globus. The full cohort is 1338 samples and roughly 1.7 TB.

Two things the pipeline expects afterwards:

  • fastq named <sample>_R1.fastq.gz and <sample>_R2.fastq.gz in one directory, which is what $RAW points at. The distribution ships one .tar per sample, so unwrap them.
  • the metadata CSV somewhere ihmp_metadata can point at. IHMPLoader reads it, and build-cohort joins it onto the abundances.

1. Profile

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

MetaPhlAn 4.2.4, database mpa_vJun23_CHOCOPhlAnSGB_202403, recorded in abundance.meta.json beside the table. 1331 samples, 1578 species.

The run log for this stage did not survive, so this command is reconstructed from the documented form and matches the recorded provenance. Every other command on this page is quoted from a log.

2. Build the cohort

mbcausal build-cohort -i $REPORTS/ihmp -o $PREPROCESS \
  --cohort ihmp --tools metaphlan4 --depth-from bracken

--depth-from bracken matters: the qc_pass keep list comes from one depth source for the whole cohort, so every profiler is filtered to identical samples. See Preprocessing.

The MetaPhlAn3 cohort was built the same way with --tools metaphlan3, reading the profiles the iHMP distribution ships rather than reports of our own.

3. Preprocess

mbcausal preprocess -i $PREPROCESS/ihmp-metaphlan4 \
  -o $PREPROCESS/ihmp-metaphlan4/clr-fecalcal-p00 --outcome fecalcal

Note the subdirectory. preprocess also writes sample_table.csv, so writing into the cohort directory would replace the complete record with the QC filtered one. p00 in the name records the prevalence filter, which is 0.0 here, the default.

4. Discover, the reported run

mbcausal discover -i $PREPROCESS/ihmp-metaphlan4/clr-fecalcal-p00 \
  -o $RUNS/nmf_7 \
  --outcome fecalcal --unit patient_id --time week \
  --subgroup diagnosis=CD,UC \
  --confounders antibiotics age \
  --n-factors 10 --target-factor nmf_7 \
  --score-scale rank --n-seeds 10 --n-perm 200

Three flags carry weight:

  • --score-scale rank applies the rank inverse normal transform to the factor scores. Raw NMF scores are skewed enough that cluster robust standard errors break down, rejecting at 33% against a nominal 5% under a known null. Every reported number depends on this.
  • --subgroup diagnosis=CD,UC restricts to participants with IBD.
  • --target-factor nmf_7 names the factor for the single factor panels. It does not affect the family the FDR correction runs over.

The resolved settings are written to spec.json in the output directory. Read that rather than reconstructing from the command line.

The reported figures were produced by two runs rather than one. The NMF loadings and the pair gap panel come from a later run made without --deny-genera, which used --steps gaps and so wrote the factorization but no effect estimates. The effect estimates come from an earlier full run that did pass --deny-genera Chlamydia. That flag changed nothing: no factor is led by that genus, none was flagged, dropped_factors is empty, and the nmf_loadings.csv and nmf_scores.csv of the two runs are identical to machine precision. The command above therefore reproduces both.

5. The MetaPhlAn3 comparison

mbcausal discover -i $PREPROCESS/ihmp-metaphlan3/clr-fecalcal-p00 \
  -o $RUNS/mpa3 \
  --outcome fecalcal --unit patient_id --time week \
  --subgroup diagnosis=CD,UC \
  --confounders antibiotics age \
  --n-factors 10 --target-factor nmf_6 \
  --score-scale rank --n-seeds 10 --n-perm 200

Identical settings, different profiler. The target factor differs because factor numbering is not portable between profilers.

Analysis specification

Every setting the reported run resolved to, from its spec.json. A dagger marks a value that came from a default rather than the command line.

Setting Value
Profiler MetaPhlAn4 4.2.4, mpa_vJun23_CHOCOPhlAnSGB_202403
Factors K / init / seed 10 / NNDSVDa / 0
Max NMF iterations 1000 †
Score scale rank inverse normal (van der Waerden)
Outcome (transform) fecalcal (log1p) †
Unit / time patient_id / week
Lag / max pair gap 1 visit / none †
Adjustment set antibiotics, age, at the same visit
Confounder lags and windows none †
Subgroup diagnosis in {CD, UC}
Exposures / errors z scored / CR1 clustered by unit †
FDR level q < 0.1, Benjamini Hochberg over the eligible family †
Dominance threshold / deny list 0.9 / empty †
Target factor nmf_7
ICC ceiling 0.75 †
Prevalence filter 0.0, nothing dropped †

--score-scale rank is not a default. It was raw before 2026-07-30 and every number produced under raw is unusable, so a spec.json showing raw marks an old run.

What the run should report back

If you rerun step 4 and these disagree, something upstream has changed:
factors built 10
passing the dominance screen 10, none dropped
eligible for the within patient design 7
outcome ICC 0.668
median factor ICC 0.68

n_significant in spec.json refers to the pooled forward scan, not the within patient arm the paper reports. Read mundlak_all.csv for the latter.

Simulation

1. Generate the cohorts

mbcausal sim generate --root $SIMDATA

Defaults give 10 effect seeds, 10 no_effect seeds and 5 disease seeds. The paper uses the first twenty. Generation is deterministic given the parameters and seed, both recorded per cohort in 2_cohort/abundance.meta.json, which is authoritative because the generator commit was dirty when the cohorts were written.

2. Preprocess and discover, per cohort

mbcausal preprocess -i $SIMDATA/effect/seed05/2_cohort \
  -o $RUNS/sim-effect-seed05 --outcome fecalcal

mbcausal discover -i $RUNS/sim-effect-seed05 \
  -o $RUNS/sim-effect-seed05 \
  --outcome fecalcal --unit patient_id --time week \
  --confounders antibiotics age \
  --n-factors 10 \
  --score-scale rank --n-seeds 10 --n-perm 200

Repeat for each of effect/seed00 to seed09 and for no_effect likewise. The settings match the real data run except for --subgroup, which does not apply, and --target-factor, which would presuppose the answer.

3. Score against the planted truth

Scoring compares each run against its truth/truth.json: how concentrated the best factor is on the planted species, how many factors were called significant, how many of those were the planted set, and how many were not. The forward screen and the within patient arm are scored separately, because they select different factors on most cohorts.

Note

The scoring script is not in this repo yet. It sits on the feat/sim-minimal branch as src/mbcausal/sim/evaluate/score.py and should be merged before this page is relied on.

Figures

scripts/plots.py builds the figure set from the CSVs a run writes, and scripts/plots_workbench.ipynb is a cell per figure over the same module.

export MBCAUSAL_RUN=$RUNS/nmf_7
export MBCAUSAL_SAMPLES=$PREPROCESS/ihmp-metaphlan4/clr-fecalcal-p00/sample_table.csv
export MBCAUSAL_SIM_RUNS=$RUNS
export MBCAUSAL_FIGS=/path/to/figures
python scripts/plots.py $MBCAUSAL_RUN --out $MBCAUSAL_FIGS

The paper's simulation forest is the mundlak_forest figure of sim-effect-seed05.

See also

Clone this wiki locally