Metadata-Version: 2.4
Name: bwo_io
Version: 0.2.0
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3 :: Only
Classifier: Programming Language :: Rust
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Typing :: Typed
Requires-Dist: numpy
Requires-Dist: gwseq-io>=0.4.1
Requires-Dist: build ; extra == 'dev'
Requires-Dist: twine ; extra == 'dev'
Requires-Dist: maturin>=1.7,<2.0 ; extra == 'dev'
Requires-Dist: bwo-io[test] ; extra == 'dev'
Requires-Dist: pytest ; extra == 'test'
Requires-Dist: mypy ; extra == 'test'
Provides-Extra: dev
Provides-Extra: test
License-File: LICENSE
Summary: Compressed genome coverage tracks (.bwo), with bigWig conversion (Rust)
Author-email: Arthur Gouhier <ajgouhier@gmail.com>
License-Expression: BSD-3-Clause
Requires-Python: >=3.11
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Repository, https://github.com/ajgouhier/bwo_io

# bwo_io

Compressed genome coverage tracks: the `.bwo` format, with bigWig conversion.
A Python package with a Rust core.

- **Lossless, 3–8× smaller than bigWig.** A coverage track is stored as runs of
  whole read counts and a scale factor (value = count × scale), coded with an
  adaptive range coder that uses the read length when the track has one.
- **Lossy with a guarantee, up to 10–90× smaller.** The error of the running
  sum stays within a bound you choose, so the mean of any window is close to
  the original, and totals stay within the bound.
- **Raw counts, any normalization.** Write a pileup's raw counts with a scale
  such as CPM; read CPM by default, raw counts on request, and change the
  scale in place.
- **bigWig in and out** through [gwseq_io](https://pypi.org/project/gwseq_io/),
  whose region syntax and binning the reads share.
- **Smaller bigWigs** for tools that only read bigWig: the same lossy bound,
  written back as a bigWig, 2–4× smaller at 10 bp, 15–40× at 1 bp.

```pycon
>>> import bwo_io
>>> report = bwo_io.from_bigwig("signal.bw", "signal.bwo")            # lossless
>>> report = bwo_io.from_bigwig("signal.bw", "small.bwo", preset="balanced")
>>> r = bwo_io.open("signal.bwo")
>>> r.read_values(["chr1:1,000,001-1,100,000"], bin_size=1000)        # CPM, as in the bigWig
>>> r.quantify(("chr1", 0, 50_000), reduce="sum", scale=1)            # raw counts
>>> bwo_io.to_bigwig("signal.bwo", "back.bw")
>>> report = bwo_io.compress_bigwig("signal.bw", "signal.small.bw")      # a bigWig, balanced
```

Install with `pip install bwo_io` (Python ≥ 3.11; the source distribution
compiles the Rust core, which needs a Rust toolchain ≥ 1.85).

## Sizes

On chr1 + chrX of four test tracks, against the bigWig (see
`local/planning/2/benchmarks.md` in the repository):

| Track | lossless | `fine` | `balanced` | `peaks` | `compact` |
|---|---|---|---|---|---|
| 1 bp, 50 bp reads | 8.0× | 29× | 48× | 68× | 90× |
| 1 bp, 50 bp reads | 7.2× | 25× | 39× | 51× | 75× |
| 10 bp, fixed 160 bp fragments | 4.8× | 6.9× | 9.4× | 14× | 15× |
| 10 bp, paired-end | 3.1× | 4.8× | 6.5× | 8.4× | 10× |

## The value model

A track's values are whole counts times one scale factor. Converting a
bigWig infers the factor: the value of one read, found as the candidate that
puts the most bases on whole multiples. Values that are not whole counts (a
short last bin's mean, a negative value) are kept exactly as exceptions; a
track that is not counts at all is stored as exact values, without the gain:
fewer than 99.9% of bases fit, or the only step that fits is too small for
6 digits to tell neighbouring counts apart (a fold-enrichment or p-value
track with many distinct values). A lossless conversion keeps every count;
each value is within the source's own precision (bigWig values went through
6 significant digits of text, so within about 5 × 10⁻⁶).

Above a few tens of thousands of reads per base, 6 digits may no longer pin
a value to one count; the conversion checks this per value and warns
(`stats.ambiguous_bases` in the report).

**Gaps.** By default, bases without data are stored as zero coverage
(`gaps="zero"`), which suits coverage and lets lossy merging cross them.
`gaps="keep"` keeps them apart from zero; reads then see them as gaps
(`def_value` applies).

## Lossy coding

At every base, the running sum of the stored counts stays within **eps** of
the original's (eps in counts × bases). So the mean of any window of W bases
is within 2·eps / W counts, and chromosome totals within eps. Optional
modifiers:

- `cap`: no base may be off by more than `cap` counts;
- `keep_zeros`: zero stretches stay exactly zero (`True`), or only those of at
  least N bases (`keep_zeros=N`).

| Preset | eps | cap | 100 bp mean within | 1 kb mean within |
|---|---|---|---|---|
| `lossless` (default) | — | — | exact | exact |
| `fine` | 50 | — | 1 count | 0.1 count |
| `balanced` | 100 | — | 2 counts | 0.2 count |
| `peaks` | 250 | 3 | 3 counts (cap) | 0.5 count |
| `compact` | 250 | — | 5 counts | 0.5 count |

On 13 test tracks, `balanced` changed the correlation between tracks by at
most 0.002 and kept replicate fold changes at r ≥ 0.999; `peaks` keeps peaks
and single bases best at the same size. The file records the bound and the
realized error, and `bwo-io verify --against source.bw` checks it at every
base.

## Smaller bigWigs

`compress_bigwig` writes a bigWig from a bigWig with the same value model
and lossy bound as `.bwo` (preset `balanced` by default), for genome
browsers and any tool that reads bigWig. The output holds the source's
units (or raw counts with `scale=1`); `verify=True` reads it back and checks
the bound at every base against the source.

On chr1 + chrX of four test tracks, against the source bigWig (see
`local/planning/2/benchmarks.md` in the repository):

| Track | `lossless` | `fine` | `balanced` | `peaks` | `compact` |
|---|---|---|---|---|---|
| 1 bp, 50 bp reads | 1.3× | 8.5× | 16× | 27× | 38× |
| 1 bp, 50 bp reads | 1.5× | 8.5× | 15× | 23× | 39× |
| 10 bp, fixed 160 bp fragments | 1.3× | 1.5× | 2.3× | 4.6× | 4.7× |
| 10 bp, paired-end | 1.4× | 1.2× | 1.8× | 2.6× | 3.4× |

A smoothed bigWig costs about 5 bytes per run, so it gains less than `.bwo`
(3–5× larger than the `.bwo` file at the same preset). `lossless` merges
equal runs and repacks (zlib 9, 4096 items per block), which already makes
most bigWigs smaller. On a binned track, gwseq_io stores unsmoothed bins
without coordinates, so a light preset can come out larger than `lossless`
(`fine` on the paired-end track above). Keeping gaps (`gaps="keep"`)
stops smoothing at every gap and can cost 2–5× on tracks with many.

`to_bigwig` takes the same presets, to write a small bigWig from a `.bwo`
file; the bound is then measured against the file's values (a lossy file's
bound and this one add up).

## Python API

```python
bwo_io.from_bigwig(input, output, *, preset="lossless", eps=None, cap=None, keep_zeros=False,
                   gaps="zero", read_length="auto", scale_name=None, chromosomes=None,
                   name=None, parallel=None, progress=None) -> dict
bwo_io.to_bigwig(input, output, *, scale=None, preset="lossless", eps=None, cap=None, keep_zeros=False,
                 chromosomes=None, compression_level=9, items_per_slot=4096, parallel=None,
                 progress=None) -> dict
bwo_io.compress_bigwig(input, output, *, preset="balanced", eps=None, cap=None, keep_zeros=False,
                       gaps="zero", scale=None, chromosomes=None, compression_level=9,
                       items_per_slot=4096, zoom_share=0.2, verify=False, parallel=None,
                       progress=None) -> dict
bwo_io.rescale(path, scale, *, name=None)
bwo_io.open(path, "r") -> BwoReader
bwo_io.open(path, "w", **writer_options) -> BwoWriter
```

`parallel` sets the number of threads, which never changes the output. By
default, `from_bigwig` and `compress_bigwig` scan 8 chromosomes at once
while collecting values and read 4 ahead of the encoder, whose runs are held
until it takes them, and `from_bigwig` encodes on all cores. A given
`parallel` replaces all three; a lower one saves memory on dense 1 bp
tracks. In `compress_bigwig` and `to_bigwig` it also sets gwseq_io's bigWig
writer threads.

`BwoReader` has `chr_sizes`, `metadata`, `scale`, `scale_name`, `lossy`,
`error_bound`, `zoom_bins`, and:

- `read_values(regions, *, span, anchor, bin_size, bin_count, snap_to_grid, full_bin, bin_mode="mean",
  def_value=0.0, zoom="full", scale=None)`: an array of shape (regions, bins);
- `quantify(regions, *, span, anchor, reduce="mean", def_value=0.0, zoom="full", scale=None)`: one
  value per region;
- `read_runs(chr, start=0, end=None, *, counts=False, scale=None)`: runs as starts, ends and values
  (float32), or integer counts with `counts=True` (−1 marks a kept gap);
- `iter_runs(*, chromosomes=None, counts=False, scale=None)`: every run, chunk by chunk;
- `layout()`: per chromosome, grid, tile span and bytes.

Regions, binning parameters and reductions (mean, sd, sem, sum, count, min,
max, l1norm, l2norm) are gwseq_io's own, as are the rules for gaps and
`def_value`. `zoom="auto"` reads the zoom levels (bins from 4 kb, ×4) when a
bin is at least 3 times a level's bin. `scale=` applies another factor for
one call.

`BwoWriter(path, *, chr_sizes, scale=None, scale_name=None, preset="lossless", eps=None, cap=None,
keep_zeros=False, gaps="zero", read_length="auto", metadata=None, parallel=None)` writes:

- `write_counts(chr, starts, ends, counts)`: raw counts, e.g. a pileup;
- `write_runs(chr, starts, ends, values)`;
- `write_values(chr, start, values, span=1)`: one value per `span` bases.

Each chromosome is written once, in order within itself. Nothing is at
`path` until `close()` returns the report; a `with` block abandons the file
on an exception.

```python
with bwo_io.open("pileup.bwo", "w", chr_sizes=sizes, scale=1e6 / mapped_reads, scale_name="CPM") as w:
    for chrom, starts, ends, counts in pileup:
        w.write_counts(chrom, starts, ends, counts)
```

Every error is a `bwo_io.Error`: `CorruptFile`, `InvalidRequest`,
`UnknownChromosome`, `UnsupportedVersion`; missing files raise
`FileNotFoundError`.

## Command line

```text
bwo-io from-bigwig IN.bw OUT.bwo [--preset P] [--eps E] [--cap C] [--keep-zeros [N]]
                                  [--gaps zero|keep] [--read-length R] [--scale-name NAME]
                                  [--parallel N]
bwo-io to-bigwig IN.bwo OUT.bw [--scale S] [--preset P] [--eps E] [--cap C] [--keep-zeros [N]]
                               [--parallel N]
bwo-io compress IN.bw OUT.bw [--preset P] [--eps E] [--cap C] [--keep-zeros [N]] [--gaps zero|keep]
                              [--scale S] [--compression-level L] [--verify] [--json]
                              [--parallel N]
bwo-io rescale FILE.bwo SCALE [--name NAME]
bwo-io info FILE.bwo [--chromosomes] [--json]
bwo-io verify FILE.bwo [--against SRC.bw]
```

Also `python -m bwo_io`.

## Format and design

`docs/format.md` specifies the file byte by byte; `docs/architecture.md`
explains the code and the measurements behind the design.

