Metadata-Version: 2.4
Name: mpath-pseudotime
Version: 0.2.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.)

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

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

## License

MIT — see [LICENSE](LICENSE).
