Metadata-Version: 2.4
Name: countmut
Version: 0.1.4
Summary: Unified ultra-fast strand-aware mutation counter (C backend + samtools-style -e/-p filters)
Author-email: Ye Chang <yech1990@gmail.com>
License-Expression: MIT
Keywords: bioinformatics,bam,pileup,mutation,bisulfite,sequencing,genomics
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
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: Programming Language :: Python :: 3.14
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Topic :: Software Development :: Libraries :: Python Modules
Classifier: Topic :: Scientific/Engineering :: Information Analysis
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: pysam>=0.21.0
Requires-Dist: rich>=13.0.0
Requires-Dist: click>=8.0.0
Requires-Dist: rich-click>=1.6.0
Requires-Dist: numpy>=1.21.0
Provides-Extra: dev
Requires-Dist: pytest>=7.0.0; extra == "dev"
Requires-Dist: pytest-cov>=4.0.0; extra == "dev"
Requires-Dist: ruff>=0.14.0; extra == "dev"
Requires-Dist: black>=22.0.0; extra == "dev"
Provides-Extra: test
Requires-Dist: pytest>=7.0.0; extra == "test"
Requires-Dist: pytest-cov>=4.0.0; extra == "test"
Dynamic: license-file

# CountMut

> **Unified ultra-fast strand-aware mutation counter** — C backend, Python wrapper.

CountMut counts base/substitution ratios from BAM files with **fast C core** and
a thin Python wrapper. It fuses the two classic ways of walking a BAM:

- **pileup-based** (`bam_mplp_auto` / pysam pileup) — fast, sees indels/ref-skips, general.
- **read-walk** (countmut's "no pileup") — walk reads directly, only touch the target sites.

Both produce **identical** output, and the tool can process whole genomes in
parallel (threads).

## Why it's fast & correct

- The hot loop (BAM read, pileup, per-(site,strand) base counting, mate-overlap
  dedup, quality/conversion classification) is in **C** (`backend/countmut_core`,
  built on the self-contained htslib subset from lh3/minipileup).
- **Strand-aware** (countmut biological-strand rule for paired-end reads).
- **Paired-end overlap dedup**: at an overlapping position a fragment is counted
  once, choosing the best mate by `(mapq, read1, base-qual)` — the thing
  minipileup gets wrong.
- **Parallel**: divides the genome into bins and processes them across threads
  (`--threads`).
- **Memory-clean** (verified under AddressSanitizer).

## Install

```bash
pip install -e .
# or, to prebuild the C core:
make backend
```

## Quick start

```bash
# strand-aware A->G mutation count (bisulfite / m6A style)
countmut -i in.bam -r ref.fa -o mut.tsv --ref-base A --mut-base G

# per-site base counts (perbase/mpileup style)
countmut -i in.bam -r ref.fa --mode base -o depth.tsv

# alleles -> VCF (minipileup style)
countmut -i in.bam -r ref.fa --mode allele --vcf -o allele.vcf
```

## Modes

| Mode | Output |
|------|--------|
| `mutation` | `chrom pos strand motif u0 u1 u2 m0 m1 m2 [o0 o1 o2]` (strand-aware substitution table) |
| `base` | `chrom pos [strand] ref depth a c g t n [ins del ref_skip fail]` |
| `allele` | `chrom pos ref depth ref_count alt alt_count`, or VCF with `--vcf` |

## Filtering with expressions (`-e` / `-p`)

Filtering is done with **samtools-style filter expressions** — there are no
separate `--min-mapq`/`--min-baseq`/`--trim-*` flags; write them as expressions
instead.

* `-e, --expression <STR>` — per-base **read** filter (samtools SAM fields).
* `-p, --pile-expression <STR>` — per-**site** filter (pileup fields).

Grammar is the samtools `filter=STRING` expression language (C-style precedence,
`&&`/`||`/`!`, bit fields, tags, regex). See `docs/filter_grammar.md`.

```bash
# keep high-quality, non-5prime, properly paired reads
countmut -i x.bam -r ref.fa -e "mapq >= 20 && bq >= 20 && dist5 >= 2 && flag & PROPER_PAIR"

# restrict to one RG group
countmut -i x.bam -r ref.fa -e "tag('RG') == 'sampleA'"

# report only A-reference sites with depth >= 5 and > 2 G alleles
countmut -i x.bam -r ref.fa -p "ref == 'A' && depth >= 5 && g > 2"
```

Read variables: `mapq`, `flag` (+ `flag.dup`, `flag.unmap`, ...), `qname`, `pos`,
`endpos`, `pnext`, `rname`, `mrname`, `tlen`, `qlen`, `rlen`, `ncigar`, `seq`,
`qual`, `sclen`, `hclen`, `bq`, `dist5`/`dist3`, `strand`, `[NM]`/`[RG]` tags,
`avg(qual)`, `exists([NM])`, `sqrt(mapq)`, ...

Site variables: `depth`, `pos`, `ref`, `a c g t n`, `ins`, `del`, `ref_skip`, `fail`.

> When `-e`/`-p` is given, counting runs on the Python engine (the C core cannot
> evaluate strings). Without expressions, the fast C backend is used.

## Engine selection

`--engine {auto|read-walk|pileup}` (default `auto`):

- `auto` → read-walk for `mutation` (targeted sites), pileup for `base`/`allele`.
- `read-walk` / `pileup` → force a strategy.

## Options

```
-i/--input, -r/--reference, -o/--output
--mode {mutation,base,allele}   --engine {auto,read-walk,pileup}
--region, --threads/-t
--ref-base, --mut-base, --pad, --save-rest
--split-strand, --count-indels, --min-depth, --min-allele-support, --vcf
-e/--expression, -p/--pile-expression
```

## Design

```
countmut/
  cli.py               rich CLI (routes to the C core directly)
  backend.py           builds/loads the C binary and calls it
  model.py             FilterConfig / MutationConfig / StrandConfig / EngineConfig
  core.py              the original pure-Python countmut (legacy reference)
backend/
  countmut_core.c      computation core: read-walk AND pileup engines,
                       mate-overlap dedup, -e/-p Lua filters (embedded lua5.4)
  countmut_expr.c      Lua filter-expression evaluator (-e / -p)
  countmut_core_main.c CLI wrapper
  Makefile             builds the `countmut_core` binary (links lua5.4)
```

Both BAM-walk strategies (`--engine read-walk` and `--engine pileup`) are
implemented in C; Python is a thin wrapper and does no counting.  `-e`/`-p`
filters are embedded-Lua expressions (see `docs/filter_grammar.md`), evaluated
in C identically by both engines.

Both walks fill the same per-site, per-strand base counts and emit byte-identical
output -- the two strategies are interchangeable, with identical results (and
identical `-e`/`-p` filtering) whether you pick `read-walk` or `pileup`.

## References this tool learns from

- [minipileup](https://github.com/lh3/minipileup) — pileup walk, filters, allele counting
- [perbase](https://github.com/sstadick/perbase) / [pbr](https://github.com/brentp/pbr) — mate-aware overlap dedup, base counts
- [countmut](https://github.com/y9c/countmut) — biological strand, bisulfite NS/Zf/Yf tiers
- [mpileup](https://github.com/y9c/mpileup) / [cpup](https://github.com/y9c/cpup) — base-count output
- samtools `--input-fmt-option filter=STRING` — the expression grammar

## License

MIT
