mpath-pseudotime 0.3.0

MPATH: methylation pseudotime analysis for Oxford Nanopore long reads

MPATH

Methylation pseudotime analysis for deciphering epigenetic cross-talk across sub-cell-cycle timescales.

MPATH calculates pseudotime for long-read methylation data from Oxford Nanopore sequencing. The full flow of data is:

Nanopore BAM (MM/ML tags)
  -> modkit extract calls         (per-CpG methylation calls)
  -> mpath metrics                (read-level metrics; merges in the WGBS reference)
  -> mpath pca fit / apply        (PCA separation of nascent vs mature reads)
  -> pseudotime score

Installation

pip install mpath-pseudotime

This installs the mpath command-line tool and the mpath Python package. No MATLAB required — the PCA has been ported to Python. (The original MATLAB PCA script is retained, unchanged, under matlab/ for reproducibility of the published figures; its full run instructions and input-table guidelines are in matlab/README.md.)

Requires Python ≥ 3.9. Dependencies (numpy, pandas, polars, scipy, joblib, tqdm, matplotlib, scikit-learn) are installed automatically.

Preliminary setup

Before MPATH, process the raw Nanopore data into per-CpG methylation calls. Exact syntax depends on your dorado / modkit / samtools versions.

1. Basecall + align with a model that emits 5mCG tags, then sort:

   dorado basecaller hac input.pod5 --modified-bases 5mCG_5hmCG > calls.bam
   dorado aligner hg19.fa calls.bam > reads_aligned.bam
   samtools sort reads_aligned.bam > reads_aligned_sorted.bam
   

2. Remove chimeric reads and index:

   samtools view -b -F0x900 reads_aligned_sorted.bam > reads_nonchimeric.bam
   samtools index reads_nonchimeric.bam
   

3. Extract read-level methylation calls with modkit (developed against modkit v0.4.x). The --read-calls output is the file MPATH consumes; it can be large, so gzip it (MPATH reads .gz directly):

   modkit extract reads_nonchimeric.bam extract_full.tsv \
     --read-calls calls_cpg.tsv \
     --ref hg38.fa \
     --include-bed CG_motifs.bed   # optional: restrict to CpG-motif positions
   pigz calls_cpg.tsv
   

   (On newer modkit the equivalent is modkit extract calls reads.bam calls_cpg.tsv.) The output has a header row; MPATH resolves the columns it needs (read_id, chrom, ref_position, call_code, ref_strand, and fail if present) by name, so added/reordered columns across modkit versions are fine.

4. BrdU filtering When working with a new cell line, we recommend performing a BrdU pulse-chase experiment to create a pseudotime PCA projection. For BrdU-labeled data, filter to BrdU-positive reads first. BrdU can be called with DNAscent; average the per-read BrdU probabilities to score each read.

5. Obtain a WGBS reference bed — the expected methylation ratio for each CpG in your cell type. You supply your own. For the MPATH study, we used WGBS, but any bulk methylation measurement in stable cells can be used. The expected form is a simple tab-separated bed; the common 4-column chrom, start, end, ratio works out of the box:

   chrom   start   end   ratio
   chr1    1061    1062  0.823
   

   The ratio column is selected with -wgbs_column (0-based; default 3), and the start must use the same genome build as modkit's ref_position. You do not need to know the file's exact coordinate convention (0- vs 1-based, which base of the CpG dyad/strand it anchors on) — MPATH auto-probes it (see below). The one thing it can't infer is the ratio scale (0.5 could be a ratio or a rounded percent): ratios must be 0–1, or pass --wgbs-scale 0-100 for percentages (MPATH warns if your values look like percentages).

If you'd rather do the intersection yourself — e.g. with bedtools in a pipeline — feed a pre-merged bed via -path_input_bed; see Pre-merged input below.)

Coordinate auto-probe + intersection QC

WGBS beds vary in convention and you usually can't tell from the file. So instead of asking you to know, mpath metrics reads a sample of CpGs and measures the match rate for each candidate coordinate offset (-1/0/+1) and strand-collapse mode, then picks the best and prints what it found:

WGBS alignment probe (offset, collapse_strands -> match):
  offset=+0  collapse_strands=True   ->  100.0%
  offset=+0  collapse_strands=False  ->  53.3%
  offset=-1  collapse_strands=False  ->  46.7%
  ...
  chosen: offset=+0, collapse_strands=True (100.0%)
WGBS intersection: matched 43630/43630 CpGs (100.0%); 0/609 reads had no WGBS overlap.

The correct convention reveals itself as a sharp jump in the matched fraction (this is self-verifying, unlike the ratio scale). If even the best alignment matches poorly, MPATH warns (it never aborts) and names the likely cause — chromosome naming (chr1 vs 1), genome build, or the wrong ratio column/scale. Force a specific convention with --wgbs-offset / --wgbs-collapse if you ever need to override the probe.

Pre-merged input

If you prefer to do the WGBS intersection yourself (e.g. bedtools intersect inside a pipeline), skip the merge entirely and pass a pre-merged 7-column bed — chrom, start, stop, strand, read_id, methylation(0/1), wgbs_ratio, grouped by read_id:

mpath metrics -path_input_bed merged.bed -path_output_csv metrics.csv -p 8

In this mode MPATH does no merge, probe, or QC — it just computes metrics on what you provide. -path_input_bed is mutually exclusive with -path_calls/-path_wgbs.

Metrics: mpath metrics

mpath metrics \
  -path_calls calls_cpg.tsv.gz \
  -path_wgbs wgbs.bed \
  -wgbs_column 3 \
  -path_output_csv metrics.csv \
  -p 8 -min_cpgs 3 -bin_limits 0,100,1000,5000,10000 --use_full_matrix

For each read, MPATH compares all CpG pairs to produce a simple matching coefficient (SMC), a uniformity score, and a Pearson correlation — overall, per genomic-distance bin, and for nearest-neighbour CpGs.

Arguments

name	description	type	required	default
-path_calls	modkit read-calls TSV (.tsv or .tsv.gz)	str	yes*	—
-path_wgbs	WGBS reference bed (.bed or .bed.gz)	str	yes*	—
-path_input_bed	pre-merged 7-col bed (alternative to calls+wgbs)	str	yes*	—
-path_output_csv	output metrics CSV	str	yes	—
-wgbs_column	0-based column of the WGBS ratio	int	no	3
-min_cpgs	minimum CpGs on a read to compute metrics	int	no	3
-bin_limits	comma-separated distance-bin limits	str	no	0,100,1000,5000,10000
-batch_size	approx CpGs per batch (controls RAM)	int	no	1e8
-p	number of parallel processes	int	no	1
--use_full_matrix	use the full CpG-pair matrix (vs one triangle)	flag	no	off
--include-hydroxy	count 5hmC (h) calls as methylated	flag	no	off
--keep-unmatched-wgbs	keep CpGs absent from the WGBS bed (NaN ratio)	flag	no	off
--wgbs-scale	scale of the WGBS ratio: 0-1 or 0-100	str	no	0-1
--wgbs-offset	coordinate offset for the join: auto/-1/0/1	str	no	auto
--wgbs-collapse	map -strand CpGs to + dyad anchor: auto/on/off	str	no	auto

* Provide either -path_calls + -path_wgbs (on-the-fly merge) or -path_input_bed (pre-merged); the two modes are mutually exclusive.

Notes. -bin_limits 0,100,1000,5000,10000 defines bins [0,100], [100,1000], [1000,5000], [5000,10000]. Larger -batch_size and more -p processes use more RAM (multiplicatively); tune to your machine. By default only CpGs present in the WGBS reference are used (inner join); pass --keep-unmatched-wgbs to keep the rest. --wgbs-offset/--wgbs-collapse are auto-probed (above) unless you force them.

Output

A wide CSV: one row per read, with the per-read columns read_id, read_wgbs_distance (mean squared difference of read vs WGBS ratio), read_meth_ratio, followed by <metric>_<bin> columns for each metric (num_pairs, smc, uniformity, pearson_r, pearson_p) and each bin (bin_all, bin_0 … bin_N, bin_closest).

Not all metrics need to feed the PCA. Different combinations may give better nascent/mature separation — choose them with mpath pca fit --columns.

PCA: mpath pca

Fit a PCA on labelled nascent + mature metric tables (run mpath metrics separately on each):

mpath pca fit \
  --nascent nascent_metrics.csv \
  --mature mature_metrics.csv \
  --out-dir pca_out/

This writes, into pca_out/:

• nascent_scores.csv, mature_scores.csv — input tables with PCA1…PCAk appended,
• pca_model.npz — the fitted transform (loadings, mean, feature columns),
• coefficients.png, explained.png, scatter.png, scatter_histogram.png — diagnostics.

By default the PCA uses every numeric metric column except identifiers and the num_pairs_* counts; restrict it with --columns smc_bin_all,uniformity_bin_all,.... Like MATLAB's pca, data is mean-centred and not scaled (pass --standardize to z-score). PC signs may be flipped relative to MATLAB; the separation is unchanged.

Downstream analysis of unlabelled data

Compute metrics for an unlabelled dataset, then project it into the PCA space that was fit on the labelled data:

mpath pca apply \
  --model pca_out/pca_model.npz \
  --input unlabelled_metrics.csv \
  --out unlabelled_scores.csv

Docker

A version-pinned image is published to GHCR for containerized pipelines:

docker run --rm -v "$PWD":/data ghcr.io/downinglab/mpath:latest \
  metrics -path_calls /data/methylation_calls.tsv.gz -path_wgbs /data/wgbs.bed \
          -path_output_csv /data/metrics.csv -p 4

The image contains MPATH only; run the upstream tools (dorado, modkit, samtools) in their own steps. See examples/pipeline.md for a full worked example, and examples/nextflow/ for a runnable reference pipeline (modkit → mpath, version-pinned) that starts from an aligned BAM produced by a standard dorado workflow.

Development

git clone https://github.com/downinglab/mpath.git
cd mpath
pip install -e ".[dev]"
pytest         # test suite
ruff check     # lint

Releases are automated with release-please: merging the rolling release PR bumps the version, updates CHANGELOG.md, tags the release, and publishes to PyPI + GHCR. The multi-Python test matrix runs on that release PR.

Citation

If you use MPATH for your research, please cite the publication (preprint currently):

A. Trinh, N. Akhtar, K. Bonsu, N. Laszik, A. Mendelevich, T. Wen, J. L. P. Morival, K. E. Diune, M. Frazeur, J. E. Vega, A. A. Gimelbrant, E. L. Read, T. L. Downing, Methylation pseudotime analysis for label-free profiling of the temporal chromatin landscape with long-read sequencing. [Preprint] (2025). https://doi.org/10.1101/2025.03.03.641287.

License

MIT — see LICENSE.