Metadata-Version: 2.4
Name: mpath-pseudotime
Version: 0.3.0
Summary: MPATH: methylation pseudotime analysis for Oxford Nanopore long reads
Project-URL: Homepage, https://github.com/downinglab/mpath
Project-URL: Repository, https://github.com/downinglab/mpath
Project-URL: Issues, https://github.com/downinglab/mpath/issues
Author: Nandor Laszik
License-Expression: MIT
License-File: LICENSE
Keywords: bioinformatics,epigenetics,methylation,nanopore,pseudotime
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.9
Requires-Dist: joblib
Requires-Dist: matplotlib
Requires-Dist: numpy
Requires-Dist: pandas
Requires-Dist: polars
Requires-Dist: scikit-learn
Requires-Dist: scipy
Requires-Dist: tqdm
Provides-Extra: dev
Requires-Dist: build; extra == 'dev'
Requires-Dist: pytest; extra == 'dev'
Requires-Dist: ruff; extra == 'dev'
Description-Content-Type: text/markdown

# 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

```bash
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/`](matlab/) for reproducibility of
the published figures; its full run instructions and input-table guidelines are in
[`matlab/README.md`](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:
   ```bash
   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:
   ```bash
   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):
   ```bash
   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](https://github.com/MBoemo/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`:

```bash
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`

```bash
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):

```bash
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:

```bash
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:

```bash
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`](examples/pipeline.md) for a full
worked example, and [`examples/nextflow/`](examples/nextflow/) for a runnable
reference pipeline (modkit → mpath, version-pinned) that starts from an aligned
BAM produced by a standard dorado workflow.

## Development

```bash
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](https://github.com/googleapis/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](https://doi.org/10.1101/2025.03.03.641287).

## License

MIT — see [LICENSE](LICENSE).
