mpath-pseudotime 0.2.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.)

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. (Optional) BrdU filtering. 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. 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).

Why no manual input.bed anymore? Earlier versions required hand-merging the methylation calls and the WGBS ratios into a single bed. mpath metrics now does that merge internally (joining on chrom + position), so you pass the modkit calls and the WGBS bed directly. (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 (escape hatch)

If you prefer to do the CpG↔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.

License

MIT — see LICENSE.