Metadata-Version: 2.4
Name: gwseq_io
Version: 0.3.0
Classifier: Development Status :: 4 - Beta
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: build ; extra == 'dev'
Requires-Dist: twine ; extra == 'dev'
Requires-Dist: maturin>=1.7,<2.0 ; extra == 'dev'
Requires-Dist: gwseq-io[test] ; extra == 'dev'
Requires-Dist: pytest ; extra == 'test'
Requires-Dist: mypy ; extra == 'test'
Requires-Dist: pandas ; extra == 'test'
Requires-Dist: scipy ; extra == 'test'
Provides-Extra: dev
Provides-Extra: test
License-File: LICENSE
Summary: Python library for processing bigWig, bigBed, BAM, CRAM and HiC files (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/gwseq_io

# gwseq_io

Read and write the files genome-scale sequencing produces — bigWig, bigBed,
BAM, CRAM and HiC — from Python or the command line. The parsing, indexing and
binning are done in Rust, and results are handed to numpy without a copy.

- **Read** bigWig and bigBed values binned over any regions, bigBed entries,
  BAM and CRAM alignments (with or without an index), and HiC contact matrices,
  from a local path or an `http(s)://` URL.
- **Write** bigWig, bigBed and BAM files, and convert bedGraph, WIG, BED, SAM,
  BAM and CRAM into them.
- **Pile up** a BAM or CRAM into a coverage bigWig (ChIP, ATAC, RNA-seq).
- **`gwseq-io`**, a command line over all of it.

Every example in this file is run by the test suite
(`crates/gwseq-io-py/tests/test_readme.py`), and every signature is checked
against the package's type stub.


## Install

- `pip install gwseq-io` installs the package and the `gwseq-io` command. numpy
  comes with it.
- Only a source distribution is published, so pip compiles the extension on the
  installing machine. That needs a Rust toolchain, **1.85 or newer**, from
  [rustup](https://rustup.rs), and Python 3.11 or newer. No C source is
  compiled and no Python headers are needed: the extension is pure Rust, built
  against the stable ABI, and needs only the linker Rust itself uses.
- `cargo install --path crates/gwseq-io-cli`, from a checkout, builds the same
  `gwseq-io` command without Python.


## Quick start

```pycon
>>> import gwseq_io
>>> with gwseq_io.open("signal.bw", "w", chr_sizes={"chr1": 10_000}) as writer:
...     writer.write_values("chr1", 0, 100, [0.5, 1.0, 1.5, 2.0])
>>> reader = gwseq_io.open("signal.bw")
>>> reader
<gwseq_io.BigWigReader chromosomes=1>
>>> reader.read_values("chr1:1-400", bin_size=100)
array([[0.5, 1. , 1.5, 2. ]], dtype=float32)
>>> reader.quantify(["chr1:1-200", "chr1:201-400"])
array([0.75, 1.75], dtype=float32)
```

The same from the command line:

```console
$ gwseq-io read signal.bw -r chr1:1-400 --bin-size 100
chr1    0    400    0.5    1    1.5    2
$ gwseq-io quantify signal.bw -r chr1:1-200 -r chr1:201-400
chr1    0    200    0.75
chr1    200  400    1.75
```


## Concepts

### Paths and URLs

Every path argument is a `str` or an `os.PathLike`. A reader also takes an
`http://` or `https://` URL, which it reads with range requests: only the
blocks a call needs are fetched, through a cache of blocks each reader keeps
(`block_size`, `cache_bytes`). A BAM's or CRAM's index is looked for at the
URL with `.bai`, `.csi` or `.crai` appended.

Reads in flight are bounded by `io_parallel`. For a local file it is the most
reads at once, 8 by default. For a URL it is where a limit **shared by every
reader of that host** starts, 16 by default: the limit halves when the server
answers 429 or 503 or resets a connection, and grows back while it does not,
up to 128 (or up to the `io_parallel` given). A server that throttles needs no
setting.

### Regions and coordinates

Every read by region takes a `RegionsLike`:

- a string, **1-based and inclusive** as samtools, UCSC and IGV write one: `"chr1:1,001-2,000"` is the thousand bases from 1000 to 1999 counted from 0. `"chr1"` is the whole chromosome, `"chr1:1,001"` runs to its end, and `,` and `_` group digits;
- a `(chr, start, end)` or `(chr, start, end, strand)` tuple, **0-based and half-open**, as every number in this package is;
- a sequence of strings and tuples, in any mix;
- a mapping or a DataFrame of columns: `chr`, `start` and `end`, or one of `start`, `end` and `center` with `span`, and optionally `strand`. Numbers are **0-based and half-open**. Aliases like chromosome, chrom, chr_id, chromStart and plural forms are accepted.

`span` re-cuts each region to that width, from the point `anchor` names:
`"start"`, `"end"` or `"center"` (genomic, not strand-aware). `anchor` may be
left out when the regions carry only one of start, end and center.

```pycon
>>> reader = gwseq_io.open("signal.bw")
>>> reader.quantify("chr1:101-200")
array([1.], dtype=float32)
>>> reader.quantify(("chr1", 100, 200))
array([1.], dtype=float32)
>>> reader.quantify({"chr": ["chr1"], "start": [100], "end": [200]})
array([1.], dtype=float32)
>>> reader.read_values({"chr": ["chr1", "chr1"], "center": [100, 300]}, span=100, bin_size=50)
array([[0.5, 1. ],
       [1.5, 2. ]], dtype=float32)
```

A walk over a file reports its regions as `Region`, a `NamedTuple` of `chr`,
`start`, `end` and `strand` (`"+"`, `"-"`, or `"."` for none).

The types the signatures below use:

```pyi
import os
from collections.abc import Callable, Mapping, Sequence
from typing import Any, Literal

import numpy as np
import numpy.typing as npt

StrPath = str | os.PathLike[str]
Progress = bool | Callable[[int, int], object] | None
Flags = int | str
Zoom = Literal["full", "auto"] | int
Anchor = Literal["start", "end", "center"]
Reduction = Literal["mean", "sd", "sem", "sum", "count", "min", "max", "l1norm", "l2norm"]
RegionTuple = tuple[str, int, int] | tuple[str, int, int, str]
RegionsLike = str | RegionTuple | Sequence[str | RegionTuple] | Mapping[str, Any] | Any  # Any: a DataFrame
ColumnDict = dict[str, Any]
```

### Threads

Two process-wide pools do the work, shared by every reader and writer:

- **the CPU pool** decodes, bins and compresses. It has one thread per core, up
  to 64, unless `GWSEQ_IO_THREADS` or `set_threads(n)` says otherwise; either
  must come before the pool is first used.
- **the IO pool** waits on files and servers. It starts threads on demand, up to
  128 (`GWSEQ_IO_IO_THREADS` or `set_io_threads(n)`, changeable at any time),
  and they exit after 30 s idle.

`parallel` caps how much of the CPU pool one call uses: a reader takes it once,
at open, and a writer or a conversion as an argument. `None` is the whole pool.
A child made by `os.fork` starts with fresh pools.

`set_max_output_bytes(n)` bounds what one call may return, 4 GiB by default: a
read whose result would be larger raises `InvalidRequest` before it reads
anything.

### Ctrl-C and progress

Every read, conversion and pileup, and the `close()` of a bigWig or bigBed
writer, runs on the pool while the calling thread waits, so **Ctrl-C** stops
it: the work is cancelled, the call raises `KeyboardInterrupt`, and a file
being written is discarded, leaving what was at its path before. The reader
stays usable.

`progress` is `None` (nothing), `True` (a bar on standard error, redrawn in
place on a terminal and a line at a time otherwise), or a callable taking
`(done, total)`. The callable runs on the calling thread, at most every
~100 ms, and once more with `done == total` when the call finishes; an
exception it raises cancels the call and is re-raised. The units are the
call's own, most often bytes of the input; `done / total` is the share done.

```pycon
>>> seen = []
>>> summary = gwseq_io.pileup_reads("reads.bam", "coverage.bw", progress=lambda done, total: seen.append((done, total)))
>>> seen[-1][0] == seen[-1][1]
True
```

### Errors

Everything the package raises is a `gwseq_io.Error`, and each class also
derives from the built-in exception that means the same thing, so
`except ValueError` and `except OSError` keep working:

```text
Error
├── InvalidRequest (ValueError)   the call asked for something impossible
│   ├── UnknownChromosome         .chromosome .available
│   ├── Closed                    a reader or writer used after close()
│   └── WriterFailed              a writer used after a write failed
├── InvalidFile                   not a file this reads            .filename
│   └── CorruptFile               its bytes contradict each other  .offset
├── SourceError (OSError)         the bytes did not arrive, or did not land
│   ├── FileNotFound (FileNotFoundError)
│   ├── PermissionDenied (PermissionError)
│   └── HttpError                 .url .status
└── Unsupported (NotImplementedError)   a feature of the format this does not implement
```

A CRAM slice that does not match the reference it is decoded against raises
`CorruptFile`. A conversion that finds its *input* malformed raises
`InvalidFile` naming the line.

```pycon
>>> reader = gwseq_io.open("signal.bw")
>>> reader.quantify("chrZ:1-100")
Traceback (most recent call last):
    ...
gwseq_io.UnknownChromosome: Chromosome chrZ not found (available: chr1)
>>> try:
...     gwseq_io.open("missing.bw")
... except FileNotFoundError as error:
...     print(error)
[Errno 2] No such file or directory: 'missing.bw'
```


## Reference

The examples run in a directory holding `signal.bw` (the quick start's),
`reads.bam` with its `.bai`, `reads.cram` with `ref.fa`, the bigBeds
`features.bb` and `genes.bb`, and `contacts.hic`: small files of the test suite,
listed under [Development](#development).

### Module functions

```pyi
@overload
def open(
    path: StrPath,
    mode: Literal["r"] = "r",
    *,
    parallel: int | None = ...,
    io_parallel: int | None = ...,
    zoom_correction: float = ...,
    block_size: int | None = ...,
    cache_bytes: int | None = ...,
    index: StrPath | None = ...,
    reference: StrPath | None = ...,
) -> BigWigReader | BigBedReader | BamReader | CramReader | HiCReader: ...
@overload
def open(
    path: StrPath,
    mode: Literal["w"],
    *,
    type: Literal["bigwig", "bigbed", "bam"] | None = None,
    **options: Any,
) -> BigWigWriter | BigBedWriter | BamWriter: ...
def get_chr_sizes(genome: str, *, full: bool = False) -> dict[str, int]: ...
def resolve_regions(
    regions: RegionsLike,
    *,
    chr_sizes: Mapping[str, int] | str | None = None,
    span: int | None = None,
    anchor: Anchor | None = None,
    bin_size: int | None = None,
    bin_count: int | None = None,
    snap_to_grid: bool | Literal["auto"] = "auto",
    full_bin: bool = False,
) -> dict[str, npt.NDArray[Any]]: ...
def set_threads(n: int) -> None: ...
def set_io_threads(n: int) -> None: ...
def set_max_output_bytes(n: int) -> None: ...
def compress_sparse_by_color(
    values: npt.ArrayLike, row: npt.ArrayLike, col: npt.ArrayLike, color_count: int
) -> list[npt.NDArray[np.uint32]]: ...
```

- `open(path)` reads the format from the file's first bytes, not its name, and
  returns the reader of that format; `options` are that class's keywords (the
  `...` above), and one it does not take is refused. `open(path, "w")` returns
  the writer `type` names, or the one the extension names: `.bw`/`.bigwig`,
  `.bb`/`.bigbed`, `.bam`.
- `resolve_regions` says what a binned read would compute, without a file:
  see [Resolving regions](#resolving-regions). `chr_sizes` is a `{name:
  size}` mapping or a bundled assembly's name (no network); with `None`,
  nothing is checked or clipped.
- `get_chr_sizes` gives an assembly's chromosome sizes: bundled for mm10, mm39,
  hg19, hg38, dm6, ce11 and sacCer3, asked of the UCSC API for any other and
  for `full=True`, which keeps the unplaced and alternate contigs.
- `set_threads`, `set_io_threads` and `set_max_output_bytes`: see
  [Threads](#threads).
- `compress_sparse_by_color` groups the cells of a sparse matrix by colour, as
  runs of adjacent cells: one array of `(row, first col, last col)` triples per
  colour, `color_count` colours of equal width spread over the values.

```pycon
>>> type(gwseq_io.open("reads.bam")).__name__
'BamReader'
>>> gwseq_io.open("reads.cram", reference="ref.fa").type
'cram'
>>> gwseq_io.open("reads.bam", zoom_correction=0.5)
Traceback (most recent call last):
    ...
gwseq_io.InvalidRequest: zoom_correction is not an option of a bam reader, which is what reads.bam needs (its options: index, parallel, io_parallel, block_size, cache_bytes)
>>> list(gwseq_io.get_chr_sizes("sacCer3").items())[:2]
[('chrI', 230218), ('chrII', 813184)]
```

### Readers

Every reader is one class per format, registered with the abstract classes
`Reader`, `BbiReader` (bigWig, bigBed) and `AlignmentReader` (BAM, CRAM), so
`isinstance(reader, gwseq_io.BbiReader)` works. A reader is frozen, may be used
from several threads at once, and pickles (it is opened again on its path). It
holds a file handle and a block cache until `close()` or the end of a `with`
block; after that a read raises `Closed`, and what was read at open stays
readable.

```pyi
class Reader:
    path: str
    type: str  # "bigwig", "bigbed", "bam", "cram" or "hic"
    closed: bool
    chr_sizes: dict[str, int]  # in the file's order
    def close(self) -> None: ...
```

Every reader's constructor takes:

- `parallel`: the most of the CPU pool one call uses. `None` is all of it.
- `io_parallel`: reads in flight (see [Paths and URLs](#paths-and-urls)).
- `block_size`, `cache_bytes`: the cache in front of the file. Blocks are
  32 KiB for a local file and 1 MiB for a URL, 128 of them by default.

### bigWig and bigBed

```pyi
class BigWigReader:
    def __init__(
        self,
        path: StrPath,
        *,
        parallel: int | None = None,
        io_parallel: int | None = None,
        zoom_correction: float = 0.3333333333333333,
        block_size: int | None = None,
        cache_bytes: int | None = None,
    ) -> None: ...

class BbiReader:
    main_header: BbiHeader
    zoom_headers: list[ZoomHeader]  # finest first
    total_summary: TotalSummary
    def resolve_regions(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        bin_size: int | None = None,
        bin_count: int | None = None,
        snap_to_grid: bool | Literal["auto"] = "auto",
        full_bin: bool = False,
    ) -> dict[str, npt.NDArray[Any]]: ...
    def read_values(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        bin_size: int | None = None,
        bin_count: int | None = None,
        snap_to_grid: bool | Literal["auto"] = "auto",
        full_bin: bool = False,
        bin_mode: Reduction = "mean",
        def_value: float = 0.0,
        zoom: Zoom = "full",
        progress: Progress = None,
    ) -> npt.NDArray[np.float32]: ...
    def quantify(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        reduce: Reduction = "mean",
        def_value: float = 0.0,
        zoom: Zoom = "full",
        progress: Progress = None,
    ) -> npt.NDArray[np.float32]: ...
    def profile(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        bin_size: int | None = None,
        bin_count: int | None = None,
        snap_to_grid: bool | Literal["auto"] = "auto",
        full_bin: bool = False,
        bin_mode: Reduction = "mean",
        def_value: float = 0.0,
        zoom: Zoom = "full",
        reduce: Reduction = "mean",
        progress: Progress = None,
    ) -> npt.NDArray[np.float32]: ...
    def iter_all_values(
        self,
        *,
        chromosomes: Sequence[str] | None = None,
        window: int = 1000000,
        bin_size: int | None = None,
        bin_mode: Reduction = "mean",
        full_bin: bool = False,
        def_value: float = 0.0,
        zoom: Zoom = "full",
        progress: Progress = None,
    ) -> RegionIterator[npt.NDArray[np.float32]]: ...
    def to_bedgraph(
        self,
        output: StrPath,
        *,
        chromosomes: Sequence[str] | None = None,
        bin_size: int | None = None,
        bin_mode: Reduction = "mean",
        full_bin: bool = False,
        def_value: float = 0.0,
        zoom: Zoom = "full",
        merge_bins: bool = True,
        progress: Progress = None,
    ) -> None: ...
    def to_wig(
        self,
        output: StrPath,
        *,
        chromosomes: Sequence[str] | None = None,
        bin_size: int | None = None,
        bin_mode: Reduction = "mean",
        full_bin: bool = False,
        def_value: float = 0.0,
        zoom: Zoom = "full",
        progress: Progress = None,
    ) -> None: ...
```

- `read_values` gives an array of shape `(regions, bins)`; `quantify` one value
  per region, reduced over exactly the bases it names; `profile` one value per
  bin, the regions binned as `read_values` bins them and each column reduced by
  `reduce` over the regions, NaN left out.
- `iter_all_values` walks whole chromosomes a `window` at a time; each step is a
  1-D array of the window's bins. The windows start at multiples of `bin_size`,
  so laid end to end they are **bit-identical** to one read of the whole
  chromosome.
- `to_bedgraph` and `to_wig` write whole chromosomes as text, the values
  `read_values` gives at the same settings. A bin holding NaN is left out of the
  file, so `def_value=float("nan")` writes the covered part alone. `merge_bins`
  writes a run of equal bins as one bedGraph line.
- On a **bigBed**, the values are the depth of coverage its entries make.
- `zoom_correction` scales the choice `zoom="auto"` makes (below).

```pycon
>>> reader = gwseq_io.open("signal.bw")
>>> reader.read_values(["chr1:1-400", ("chr1", 100, 500, "-")], bin_size=200)
array([[0.75, 1.75],
       [1.  , 1.25]], dtype=float32)
>>> reader.profile([("chr1", 0, 200), ("chr1", 200, 400)], bin_size=100)
array([1. , 1.5], dtype=float32)
>>> walk = reader.iter_all_values(bin_size=1000, window=5000)
>>> len(walk), walk.regions[0]
(2, Region(chr='chr1', start=0, end=5000, strand='.'))
>>> for region, values in walk.items():
...     print(region.start, values)
0 [0.5 0.  0.  0.  0. ]
5000 [0. 0. 0. 0. 0.]
>>> reader.to_bedgraph("signal.bedgraph")
>>> print(open("signal.bedgraph").read(), end="")
chr1	0	100	0.5
chr1	100	200	1
chr1	200	300	1.5
chr1	300	400	2
chr1	400	10000	0
```

#### Binning

A **gap base** is a base in `[0, chr_size)` that no record covers. Bases
outside the chromosome are neither data nor gap: a bin crossing the chromosome
end is simply shorter.

| Parameter | Applies to | Default | Meaning |
|---|---|---|---|
| `bin_size` | read_values, profile, walks, exports | 1 when `bin_count` is unset | Bin width, laid from each region's own start (its 5′ end for "−"). A row has `floor(span/bin_size)` bins, or `ceil` with `full_bin`. Every region must give the same count, or the call raises naming two spans and suggesting `bin_count`. A size of 0 is refused. Not a quantify parameter. |
| `bin_count` | read_values, profile | unset | Each region is cut into exactly `n` bins of `span/n`, which may be fractional; weights are exact. The way to read regions of unequal span into one matrix. Refused with `bin_size`, `full_bin` or `snap_to_grid=True`. |
| `snap_to_grid` | read_values, profile | "auto" | Moves each row onto the genome-wide grid of `bin_size` (its multiples from base 0): the region's 5′ end goes to the nearest multiple, the lower one on a tie, and the row of `n` bins moves with it, so `n` still comes from the span. `True` always, `False` never, `"auto"` when a region has at least 10 bins, where the move (at most half a bin) is at most 5% of it. Every region of a call has the same `n`, so its rows all snap or none does. Ignored with `bin_count`, refused there as `True`. |
| `full_bin` | read_values, profile, walks, exports | False | Keeps a last partial bin at full width, clipped at the chromosome end. In walks and exports it keeps a chromosome's short last bin. |
| `zoom` | every value read, bigBed included | "full" | "auto" picks the coarsest level with `reduction ≤ w × zoom_correction`, where `w` is the narrowest output bin (for quantify, the region's span). An int picks a level. |
| `def_value` | every value read | 0 | The value of every gap base. **Finite:** gap bases join every reduction (UCSC's "mean0"). **NaN:** gap bases are left out, and a bin or region with no data is NaN. |
| `bin_mode` | read_values, profile, walks, exports | "mean" | Any reduction (mean, sd, sem, sum, count, min, max, l1norm, l2norm), over the bases of one bin: the same function quantify applies over a region. |
| `count` | everywhere | — | Covered bases, the sum of record weights. It never includes def-filled bases. |
| `reduce` (quantify) | quantify | "mean" | Over the region's bases, read exactly `[start, end)` and never snapped. |
| `reduce` (profile) | profile | "mean" | Folds each region's finished bin value per column, skipping NaN. |
| l1norm with zoom | quantify, `bin_mode` | — | An explicit level is refused; "auto" resolves to full data. |
| strand "−" | read_values, profile | — | Bins are laid from the end, so bin 0 is the 5′ end, and any remainder falls at the 3′ side. A snapped "−" row moves its end, the 5′ end. Ignored by quantify and entries. |
| NaN in a file | everywhere | — | A gap in every reduction. `sd` and `sem` use the population (n) variance. |

```pycon
>>> reader = gwseq_io.open("signal.bw")
>>> reader.read_values([("chr1", 50, 400)], bin_size=100)
array([[0.75, 1.25, 1.75]], dtype=float32)
>>> reader.read_values([("chr1", 50, 400)], bin_size=100, full_bin=True)
array([[0.75, 1.25, 1.75, 1.  ]], dtype=float32)
>>> reader.read_values([("chr1", 50, 400)], bin_size=100, snap_to_grid=True)
array([[0.5, 1. , 1.5]], dtype=float32)
>>> reader.read_values([("chr1", 30, 1030)], bin_size=100)
array([[0.5, 1. , 1.5, 2. , 0. , 0. , 0. , 0. , 0. , 0. ]], dtype=float32)
>>> reader.read_values([("chr1", 30, 1030)], bin_size=100, snap_to_grid=False)
array([[0.65, 1.15, 1.65, 1.4 , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ]],
      dtype=float32)
>>> reader.read_values(["chr1:1-400", "chr1:1-200"], bin_count=2)
array([[0.75, 1.75],
       [0.5 , 1.  ]], dtype=float32)
>>> reader.read_values("chr1:1-500", bin_size=250, def_value=float("nan"))
array([[0.9      , 1.8333334]], dtype=float32)
>>> reader.quantify("chr1:1-500", reduce="count")
array([400.], dtype=float32)
```

#### Resolving regions

`resolve_regions` gives what a binned read would compute, without reading
anything: the regions after `span` and `anchor`, whether `snap_to_grid` moved
them, and every bin. It takes the arguments of `read_values` and refuses what
it refuses. Every reader has it, reading against its own chromosomes, and so
does the module, given `chr_sizes`.

The result is a dict of numpy arrays, one row per region in the order given:

| Key | Shape | What |
|---|---|---|
| `chr`, `strand` | (regions,) | the chromosome as the file names it, and the strand given (`"."` for none) |
| `start`, `end` | (regions,) | the region after `span` and `anchor` |
| `snapped`, `shift` | (regions,) | whether the row was snapped, and how far its 5′ end moved |
| `bin_start`, `bin_end` | (regions, bins) | each cell's bases, in `read_values`' column order, clipped to the chromosome; whole with `bin_size`, float with `bin_count` |
| `bin_bases` | (regions, bins) | the bases each cell is reduced over |

`read_values(...)[i, j]` is the reduction of `bin_start[i, j]` to
`bin_end[i, j]`.

```pycon
>>> reader = gwseq_io.open("signal.bw")
>>> resolved = reader.resolve_regions([("chr1", 30, 1030)], bin_size=100)
>>> resolved["snapped"], resolved["shift"]
(array([ True]), array([-30]))
>>> resolved["bin_start"][0, :3], resolved["bin_end"][0, :3]
(array([  0, 100, 200]), array([100, 200, 300]))
>>> gwseq_io.resolve_regions(("chr1", 9950, 10_150), chr_sizes={"chr1": 10_000}, bin_size=100)["bin_bases"]
array([[50.,  0.]])
```

#### bigBed entries

```pyi
class BigBedReader:
    def __init__(
        self,
        path: StrPath,
        *,
        parallel: int | None = None,
        io_parallel: int | None = None,
        zoom_correction: float = 0.3333333333333333,
        block_size: int | None = None,
        cache_bytes: int | None = None,
    ) -> None: ...
    auto_sql: dict[str, str]
    extra_indexes: list[str]
    def read_entries(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        col_count: int | None = None,
        as_cols: bool = False,
        progress: Progress = None,
    ) -> list[list[BedEntry]] | list[ColumnDict]: ...
    def read_all_entries(
        self,
        *,
        chromosomes: Sequence[str] | None = None,
        col_count: int | None = None,
        as_cols: bool = False,
        progress: Progress = None,
    ) -> list[BedEntry] | ColumnDict: ...
    def iter_all_entries(
        self,
        *,
        chromosomes: Sequence[str] | None = None,
        col_count: int | None = None,
        as_cols: bool = False,
        window: int = 1000000,
        progress: Progress = None,
    ) -> RegionIterator[list[BedEntry]] | RegionIterator[ColumnDict]: ...
    def search(
        self,
        values: str | Sequence[str],
        *,
        field: str = "name",
        col_count: int | None = None,
        as_cols: bool = False,
    ) -> list[BedEntry] | ColumnDict: ...
    def to_bed(
        self,
        output: StrPath,
        *,
        chromosomes: Sequence[str] | None = None,
        col_count: int | None = None,
        progress: Progress = None,
    ) -> None: ...
```

- `auto_sql` is each column and the type the file declares for it, in the file's
  order. `extra_indexes` are the columns the file carries an index on: the
  fields `search` looks entries up by. `search` gives the entries whose `field`
  is one of `values`, in the file's order.
- `read_entries` gives one list of entries per region; an entry belongs to every
  region it overlaps. `read_all_entries` and `iter_all_entries` give every entry
  of the named chromosomes (all of them by default) once, in file order; a
  walk's entry belongs to the window it starts in.
- `col_count` is how many columns to read, the three coordinates included;
  `None` reads them all, and the columns left out are never parsed.
- `as_cols=True` gives one dict of columns instead of a list of entries: the
  numeric columns as numpy arrays, the rest as lists of `str`, typed from the
  autoSql. `pandas.DataFrame(columns)` takes it as it is.

```pycon
>>> reader = gwseq_io.open("features.bb")
>>> reader.extra_indexes
['name']
>>> entries = reader.read_entries("chr1:1-3000")[0]
>>> entries
[BedEntry(chr1:846-4649, 5 columns), BedEntry(chr1:1336-3623, 5 columns), BedEntry(chr1:2330-2492, 5 columns)]
>>> entries[1]["name"], entries[1]["signal"]
('LOC472203_1', 205.43)
>>> reader.read_entries("chr1:1-3000", as_cols=True)[0]["start"]
array([ 846, 1336, 2330])
>>> reader.search("LOC472203_1")
[BedEntry(chr1:1336-3623, 5 columns)]
>>> reader.quantify("chr1:1-3000", reduce="max")
array([3.], dtype=float32)
```

A `BedEntry` holds the coordinates as attributes and every column read by key.
Its values are typed from the autoSql: `int` for an integer column, `float` for
`float` and `double`, `str` for everything else and for a cell that does not
read as the number its column declares. It is frozen and pickles; `keys()`,
`len()`, `in` and iteration work as on a dict, so `dict(entry)` builds one.

```pyi
class BedEntry:
    chr: str
    start: int
    end: int
    def keys(self) -> list[str]: ...
    def to_dict(self) -> dict[str, Any]: ...
```

```pycon
>>> gwseq_io.open("genes.bb").read_entries("chr1:1-400", col_count=6)[0][0].to_dict()
{'chr': 'chr1', 'start': 323, 'end': 3048, 'name': 'tx0', 'score': 633, 'strand': '.'}
```

### BAM and CRAM

```pyi
class BamReader:
    def __init__(
        self,
        path: StrPath,
        *,
        index: StrPath | None = None,
        parallel: int | None = None,
        io_parallel: int | None = None,
        block_size: int | None = None,
        cache_bytes: int | None = None,
    ) -> None: ...
    def read_entries(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
        require_flags: Flags = 0,
        min_mapq: int = 0,
        drop_improper_pairs: bool = True,
        tags: bool = True,
        as_cols: bool = False,
        tag_cols: Sequence[str] = (),
        progress: Progress = None,
    ) -> list[list[SamEntry]] | list[ColumnDict]: ...
    def read_all_entries(
        self,
        *,
        chromosomes: Sequence[str] | None = None,
        unplaced: bool = False,
        exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
        require_flags: Flags = 0,
        min_mapq: int = 0,
        drop_improper_pairs: bool = True,
        tags: bool = True,
        as_cols: bool = False,
        tag_cols: Sequence[str] = (),
        progress: Progress = None,
    ) -> list[SamEntry] | ColumnDict: ...
    def iter_entries(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        sort_regions: bool = False,
        exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
        require_flags: Flags = 0,
        min_mapq: int = 0,
        drop_improper_pairs: bool = True,
        tags: bool = True,
        as_cols: bool = False,
        tag_cols: Sequence[str] = (),
        progress: Progress = None,
    ) -> RegionIterator[list[SamEntry]] | RegionIterator[ColumnDict]: ...
    def iter_all_entries(
        self,
        *,
        chromosomes: Sequence[str] | None = None,
        unplaced: bool = False,
        window: int = 1000000,
        exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
        require_flags: Flags = 0,
        min_mapq: int = 0,
        drop_improper_pairs: bool = True,
        tags: bool = True,
        as_cols: bool = False,
        tag_cols: Sequence[str] = (),
        progress: Progress = None,
    ) -> RegionIterator[list[SamEntry]] | RegionIterator[ColumnDict]: ...

class AlignmentReader:
    header: list[dict[str, Any]]
    is_indexed: bool
    index_path: str | None
    index_kind: Literal["bai", "csi", "crai"] | None
    index_error: str | None
    warnings: list[str]
    def resolve_regions(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        bin_size: int | None = None,
        bin_count: int | None = None,
        snap_to_grid: bool | Literal["auto"] = "auto",
        full_bin: bool = False,
    ) -> dict[str, npt.NDArray[Any]]: ...
```

- `index` is the index's path when it is not beside the file as `<path>.bai`,
  `<path>.csi` or `<path>.crai`. `index_path` and `index_kind` say which one
  is in use; `index_error` says why an index that was there could not be read,
  and is `None` when one loaded or there is none. A local CRAM with no `.crai`
  is indexed from its container headers.
- `read_entries` and `iter_entries` need an index. `read_all_entries` and
  `iter_all_entries` walk whole references and need none (a walk needs an
  index or a file sorted by coordinate). `unplaced=True` adds the alignments
  with no reference after the others, which the default filter drops as
  unmapped. `sort_regions` reads the regions in the file's order, and the
  iterator's `order` says which region of the request each step is.
- `resolve_regions` lays a binned read's bins over the file's references, as
  a bigWig reader's does (see [Resolving regions](#resolving-regions)): the
  bins `read_values` would read from a pileup of this file.
- `header` is the header's lines, each `{"type": "SQ", "fields": {...}}`.
  `warnings` is what was found wrong at open without stopping it: no
  end-of-file marker, an index that could not be read where another could.
- `tags=False` leaves the optional fields unread. `tag_cols` gives the named
  optional fields a column each, with `as_cols`; an alignment without one gets
  `None`.

**The read filter.** Every read of alignments, the pileup and the conversion
to BAM take the same four keywords, and by default drop unmapped, secondary and
supplementary alignments and improper pairs, and keep duplicates and QC
failures:

- `exclude_flags`: drop an alignment carrying any of these flags;
- `require_flags`: drop one not carrying all of these;
- `min_mapq`: drop one whose MAPQ is lower;
- `drop_improper_pairs`: drop a paired read (0x1) that is not a proper pair
  (0x2 unset), whose template length is 0, or whose mate is on another
  reference.

Flags are an `int` or samtools flag names joined by commas (`"UNMAP,DUP"`).

```pycon
>>> reader = gwseq_io.open("reads.bam")
>>> reader.index_kind, reader.chr_sizes
('bai', {'chr1': 200000, 'chr2': 100000, 'chr3': 5000})
>>> reader.read_entries("chr1:1,001-2,300")[0]
[<gwseq_io.SamEntry chr1:1000-1050 plain_fwd>, <gwseq_io.SamEntry chr1:1100-1150 plain_rev>, <gwseq_io.SamEntry chr1:2000-2050 pair_fr>, <gwseq_io.SamEntry chr1:2200-2250 pair_fr>]
>>> [entry.read_name for entry in reader.read_entries("chr1:1,001-2,300", min_mapq=40)[0]]
['plain_fwd', 'pair_fr', 'pair_fr']
>>> reader.read_entries("chr1:6,001-9,100")[0]
[]
>>> len(reader.read_entries("chr1:6,001-9,100", drop_improper_pairs=False)[0])
8
>>> columns = reader.read_entries("chr1:1,001-2,300", as_cols=True, tag_cols=["NM"])[0]
>>> columns["start"], columns["NM"]
(array([1000, 1100, 2000, 2200]), [0, 1, None, None])
>>> walk = reader.iter_entries(["chr2:1-1000", "chr1:1-3000"], sort_regions=True)
>>> walk.order, [len(entries) for entries in walk]
(array([1, 0]), [5, 0])
>>> len(reader.read_all_entries(exclude_flags=0, unplaced=True))
46
```

```pyi
class CramReader:
    def __init__(
        self,
        path: StrPath,
        *,
        index: StrPath | None = None,
        reference: StrPath | None = None,
        parallel: int | None = None,
        io_parallel: int | None = None,
        block_size: int | None = None,
        cache_bytes: int | None = None,
    ) -> None: ...
    version: str
    references: dict[str, dict[str, str | None]]
    def read_entries(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
        require_flags: Flags = 0,
        min_mapq: int = 0,
        drop_improper_pairs: bool = True,
        tags: bool = True,
        sequence: bool = True,
        as_cols: bool = False,
        tag_cols: Sequence[str] = (),
        progress: Progress = None,
    ) -> list[list[SamEntry]] | list[ColumnDict]: ...
    def read_all_entries(
        self,
        *,
        chromosomes: Sequence[str] | None = None,
        unplaced: bool = False,
        exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
        require_flags: Flags = 0,
        min_mapq: int = 0,
        drop_improper_pairs: bool = True,
        tags: bool = True,
        sequence: bool = True,
        as_cols: bool = False,
        tag_cols: Sequence[str] = (),
        progress: Progress = None,
    ) -> list[SamEntry] | ColumnDict: ...
    def iter_entries(
        self,
        regions: RegionsLike,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        sort_regions: bool = False,
        exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
        require_flags: Flags = 0,
        min_mapq: int = 0,
        drop_improper_pairs: bool = True,
        tags: bool = True,
        sequence: bool = True,
        as_cols: bool = False,
        tag_cols: Sequence[str] = (),
        progress: Progress = None,
    ) -> RegionIterator[list[SamEntry]] | RegionIterator[ColumnDict]: ...
    def iter_all_entries(
        self,
        *,
        chromosomes: Sequence[str] | None = None,
        unplaced: bool = False,
        window: int = 1000000,
        exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
        require_flags: Flags = 0,
        min_mapq: int = 0,
        drop_improper_pairs: bool = True,
        tags: bool = True,
        sequence: bool = True,
        as_cols: bool = False,
        tag_cols: Sequence[str] = (),
        progress: Progress = None,
    ) -> RegionIterator[list[SamEntry]] | RegionIterator[ColumnDict]: ...
```

A `CramReader` reads CRAM 2.1, 3.0 and 3.1 and has `BamReader`'s methods, each
with `sequence`. A CRAM stores a mapped read as its differences from the
reference, so the bases come from a FASTA (with its `.fai`, or bgzipped with
its `.gzi` too), looked for per sequence in this order: the bases the file
embeds; `reference`; the FASTA the `@SQ` line's `UR` names; `REF_CACHE` and
`REF_PATH` by the `M5` checksum (a URL in `REF_PATH` only when set
explicitly). `references` says, per sequence, which one was found:
`{"source", "path", "error"}`, where `source` is `"embedded"`, `"given"`,
`"header"`, `"cache"`, `"not_required"` or `"missing"`. A sequence with no
reference, or a read with `sequence=False`, reads every field but the bases the
reference supplies, which read `N` (and `MD`/`NM` unless the file stored them).
The slice's MD5 is checked against the reference, and a mismatch raises
`CorruptFile`. `version` is `"major.minor"`.

```pycon
>>> reader = gwseq_io.open("reads.cram", reference="ref.fa")
>>> reader.version, reader.references["chr1"]["source"]
('3.1', 'given')
>>> entry = reader.read_entries("chr1:1-120")[0][1]
>>> entry.cigar, entry.sequence[:12], entry.tag("MD")
('250M', 'ATACGATGATAC', '78A125C45')
>>> gwseq_io.open("reads.cram").references["chr1"]["source"]
'missing'
>>> gwseq_io.open("reads.cram").read_entries("chr1:1-120")[0][1].sequence[:12]
'NNNNNNNNNNNN'
```

A `SamEntry` is one alignment, frozen; it pickles its record's bytes. A CRAM
record is rebuilt into what a BAM holds.

```pyi
class SamEntry:
    chr: str
    start: int  # 0-based
    end: int  # from the CIGAR; equal to start for an alignment covering no reference
    read_name: str
    flag: int
    mapping_quality: int
    cigar: str
    sequence: str  # "*" when the record holds none
    qualities: str  # phred+33, "*" when the record holds none
    next_chr: str  # "*" when the mate has no reference
    next_start: int
    template_length: int
    bai_bin: int
    reference_length: int
    query_length: int
    tags: dict[str, Any]
    is_paired: bool
    is_proper_pair: bool
    is_mapped: bool
    is_next_mapped: bool
    is_reverse: bool
    is_next_reverse: bool
    is_first_in_pair: bool
    is_last_in_pair: bool
    is_secondary_or_supplementary: bool
    is_failed_qc_or_duplicate: bool
    def tag(self, name: str) -> Any: ...
    def to_dict(self) -> dict[str, Any]: ...
```

`tags` is the optional fields by two-letter name, in the record's order, typed
as the file types them; it is the one attribute that can raise, on a record
whose optional fields are malformed. `tag(name)` reads one field, or `None`.
`to_dict()` has the same keys for every entry, `"tags"` included. With
`as_cols`, the columns are `to_dict()`'s keys but `tags`.

```pycon
>>> entry = gwseq_io.open("reads.bam").read_entries("chr1:2,001-2,050")[0][0]
>>> entry.read_name, entry.flag, entry.is_proper_pair, entry.next_start, entry.template_length
('pair_fr', 99, True, 2200, 250)
```

### HiC

```pyi
class HiCReader:
    def __init__(
        self,
        path: StrPath,
        *,
        parallel: int | None = None,
        io_parallel: int | None = None,
        block_size: int | None = None,
        cache_bytes: int | None = None,
    ) -> None: ...
    header: HiCHeader
    footer: HiCFooter
    normalizations: list[str]
    units: list[str]
    bin_sizes: list[int]
    def resolve_regions(
        self,
        region1: str | tuple[str, int, int],
        region2: str | tuple[str, int, int] | None = None,
        *,
        bin_size: int | None = None,
        bin_count: int | tuple[int, int] | None = None,
        snap_to_grid: bool | Literal["auto"] = "auto",
        full_bin: bool = False,
    ) -> dict[str, Any]: ...
    def read_values(
        self,
        region1: str | tuple[str, int, int],
        region2: str | tuple[str, int, int] | None = None,
        *,
        bin_size: int | None = None,
        bin_count: int | tuple[int, int] | None = None,
        snap_to_grid: bool | Literal["auto"] = "auto",
        full_bin: bool = False,
        def_value: float = 0.0,
        triangle: bool = False,
        min_distance: int | None = None,
        max_distance: int | None = None,
        normalization: str = "none",
        matrix_type: Literal["observed", "oe", "expected"] = "observed",
        unit: Literal["bp", "frag"] = "bp",
        save_to: StrPath | None = None,
        progress: Progress = None,
    ) -> npt.NDArray[np.float32] | None: ...
    def read_sparse_values(
        self,
        region1: str | tuple[str, int, int],
        region2: str | tuple[str, int, int] | None = None,
        *,
        bin_size: int | None = None,
        bin_count: int | tuple[int, int] | None = None,
        snap_to_grid: bool | Literal["auto"] = "auto",
        full_bin: bool = False,
        triangle: bool = False,
        min_distance: int | None = None,
        max_distance: int | None = None,
        normalization: str = "none",
        matrix_type: Literal["observed", "oe", "expected"] = "observed",
        unit: Literal["bp", "frag"] = "bp",
        save_to: StrPath | None = None,
        progress: Progress = None,
    ) -> dict[str, Any] | None: ...
```

A Juicer `.hic` file, versions 6 to 9.

- `read_values` gives the contact matrix of `region1` (the rows) against
  `region2` (the columns; `region1` when left out) as a float32 array.
  `read_sparse_values` gives it as a dict of `"data"` (float32), `"row"` and
  `"col"` (uint32) and `"shape"`: what `scipy.sparse.coo_array((m["data"],
  (m["row"], m["col"])), shape=m["shape"])` takes.
- Each region is binned as a bigWig region is (see [Binning](#binning)):
  - `bin_size` is one of `bin_sizes`, the finest when left out. A region has
    `floor(span / bin_size)` bins from its start, and `full_bin` keeps a last
    partial one, valued over its part inside the chromosome.
  - `snap_to_grid` moves a region's start to the nearest multiple of
    `bin_size`: `True` always, `False` never, `"auto"` (the default) when the
    region has at least 10 bins. Rows and columns decide on their own.
  - `bin_count` cuts the regions into exactly that many bins instead, along
    the longer region for one number (the other gets its share), along each
    for a `(rows, cols)` pair, read from the coarsest resolution giving at
    least that many. It is refused beside `bin_size`, `full_bin` and
    `snap_to_grid=True`.
  - A bin that is one of the file's is its cell as stored. Any other (a region
    off the grid, or a bin count's bins) is the mean of the file's cells it
    covers, weighted by how much of each, NaN left out.
  - From 0.2.x: a region off the grid no longer gives every file bin it
    touches, and `exact_bin_count` is gone, a `bin_count` read being exact.
- `resolve_regions` gives the bins of a read without reading: `"rows"` and
  `"cols"`, each what a bigWig reader's `resolve_regions` gives for one
  region, and `"bin_size"`, the resolution the cells come from.
- `def_value` is the value of a cell no contact reached; a cell the chosen
  normalization cannot value is NaN.
- `triangle` reads the one side of the diagonal the file stores, without
  mirroring it; `min_distance` and `max_distance` keep the contacts that far
  from the diagonal, in bp.
- `normalization` is `"none"` or one of `normalizations`; `matrix_type` is
  `"observed"`, `"oe"` (observed over expected) or `"expected"`. `unit="frag"`
  raises `Unsupported`; the fragment resolutions and restriction sites are in
  `header`.
- With `save_to`, the matrix is written there as a `.npz` and `None` is
  returned: under the key `"values"`, or for a sparse read as
  `scipy.sparse.save_npz` writes one.

```pycon
>>> reader = gwseq_io.open("contacts.hic")
>>> reader.bin_sizes, reader.normalizations
([100000, 25000, 10000, 1000], ['none', 'scale', 'vc', 'vc_sqrt'])
>>> reader.read_values("chr1:1-300,000", bin_size=100_000)
array([[367.,  96.,  50.],
       [ 96., 414., 113.],
       [ 50., 113., 428.]], dtype=float32)
>>> reader.read_values("chr1:1-500,000", "chr2:1-300,000", bin_size=100_000, normalization="vc").shape
(5, 3)
>>> reader.read_values("chr1", bin_count=4).shape
(4, 4)
>>> reader.read_values("chr1:50,001-250,000", bin_size=100_000)
array([[243.25, 168.25],
       [168.25, 267.  ]], dtype=float32)
>>> sparse = reader.read_sparse_values("chr1:1-200,000", bin_size=100_000)
>>> sparse["data"], sparse["row"], sparse["col"], sparse["shape"]
(array([367.,  96.,  96., 414.], dtype=float32), array([0, 0, 1, 1], dtype=uint32), array([0, 1, 0, 1], dtype=uint32), (2, 2))
```

#### Quantifying contacts

```pyi
class HiCReader:
    def quantify(
        self,
        regions1: RegionsLike,
        regions2: RegionsLike | None = None,
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        paired: bool = False,
        bin_size: int | None = None,
        reduce: Reduction = "mean",
        def_value: float = 0.0,
        min_distance: int | None = None,
        max_distance: int | None = None,
        normalization: str = "none",
        matrix_type: Literal["observed", "oe", "expected"] = "observed",
        unit: Literal["bp", "frag"] = "bp",
        progress: Progress = None,
    ) -> npt.NDArray[np.float32]: ...
    def quantify_partners(
        self,
        regions: RegionsLike,
        partners: Sequence[RegionsLike],
        *,
        span: int | None = None,
        anchor: Anchor | None = None,
        bin_size: int | None = None,
        reduce: Reduction = "mean",
        def_value: float = 0.0,
        min_distance: int | None = None,
        max_distance: int | None = None,
        normalization: str = "none",
        matrix_type: Literal["observed", "oe", "expected"] = "observed",
        unit: Literal["bp", "frag"] = "bp",
        progress: Progress = None,
    ) -> list[npt.NDArray[np.float32]]: ...
```

One value per pair of regions: `reduce` over the file's cells under the
rectangle `region1 × region2`, at `bin_size` (the finest when left out).

- The regions are any `RegionsLike`, `span` and `anchor` applied to every
  list; the strand is ignored.
- `quantify` gives every region of `regions1` against every region of
  `regions2` as an array of shape `(len(regions1), len(regions2))`; with
  `regions2` left out, the regions against each other, a symmetric matrix.
  With `paired=True` it zips the two lists, which need the same length, and
  gives one value per pair, as for the rows of a BEDPE.
- `quantify_partners` gives each region of `regions` against its own
  partners: `partners[i]` is a `RegionsLike` of any length, and the result is
  a list of one float32 array per region, `[i][k]` being `regions[i]` against
  `partners[i][k]`. Each value is what `quantify` gives that pair.
- A rectangle is read exactly, never snapped. A file bin a region covers in
  part weighs the share of it the region covers, and a cell the product of
  its two sides' shares, so `mean` is the value `read_values` gives the same
  rectangle as one bin.
- The reductions are those of a bigWig `quantify`, over the cells:

  | `reduce` | Value |
  |---|---|
  | `mean` | The weighted mean of the cells. |
  | `sum` | Σ share × value: on the grid, the plain sum of the cells, so observed counts sum to the contacts between the two regions. |
  | `count` | The cells a contact reached, weighted alike. It never counts a cell filled with `def_value`. |
  | `sd`, `sem` | The weighted population standard deviation, and it over √`count`, `def_value` cells included. |
  | `min`, `max` | Over the cells the rectangle reaches. |
  | `l1norm`, `l2norm` | Σ share × \|value\|, and √(Σ share × value²). |

- `def_value` is the value of a cell no contact reached, and joins every
  reduction but `count`; NaN leaves such cells out. A cell the normalization
  cannot value is left out of every reduction, and a pair with nothing left,
  or with a region that has no base on its chromosome, is NaN.
- On one chromosome the stored triangle is mirrored, the diagonal counted
  once. `min_distance` and `max_distance` keep the cells of such a pair that
  far from the diagonal, in bp; a pair on two chromosomes is quantified whole.
- `normalization`, `matrix_type` and `unit` are those of `read_values`.
- Each block of the file is read once per call, whatever the number of pairs
  it feeds, and the result does not depend on the number of threads.

```pycon
>>> reader = gwseq_io.open("contacts.hic")
>>> reader.quantify(["chr1:1-100,000", "chr1:200,001-300,000"], bin_size=100_000, reduce="sum")
array([[367.,  50.],
       [ 50., 428.]], dtype=float32)
>>> reader.quantify(["chr1:1-100,000"], ["chr1:200,001-300,000"], paired=True, bin_size=100_000)
array([50.], dtype=float32)
>>> reader.quantify("chr1:50,001-250,000", bin_size=100_000)
array([[211.6875]], dtype=float32)
>>> reader.quantify_partners(
...     ["chr1:1-100,000", "chr1:100,001-200,000"],
...     [["chr1:100,001-200,000", "chr1:200,001-300,000"], ["chr1:200,001-300,000"]],
...     bin_size=100_000,
...     reduce="sum",
... )
[array([96., 50.], dtype=float32), array([113.], dtype=float32)]
```

```pyi
class HiCHeader:
    version: int
    genome_id: str
    footer_position: int
    nvi_position: int
    nvi_length: int
    bp_resolutions: list[int]
    frag_resolutions: list[int]
    sites: dict[str, list[int]]
    attributes: dict[str, str]

class HiCFooter:
    byte_count_v5: int
    normalizations: list[str]
    units: list[str]
    master_index: dict[str, HiCIndexItem]

class HiCIndexItem:
    position: int
    size: int
```

### Writers

Every writer is frozen, does not pickle, and is registered with `Writer`.
Nothing at its path is the new file until `close()` returns: the file is written
beside it and moved into place, so a writer that fails, is abandoned, or is
interrupted leaves what the path held before. A `with` block closes the writer
when it ends normally and **abandons** it when it ends with an exception; to
keep the file anyway, call `close()` in your own `except`. A writer that is
never closed is abandoned when it is collected.

```pyi
class Writer:
    path: str
    type: str  # "bigwig", "bigbed" or "bam"
    closed: bool  # finished, or abandoned
    state: Literal["open", "failed", "closed"]
    chr_sizes: dict[str, int]
    def abandon(self) -> None: ...
```

`state` is `"failed"` once a write or the close has failed and the file has
been discarded; a failed writer raises `WriterFailed`, a closed one `Closed`.

#### BigWigWriter

```pyi
class BigWigWriter:
    def __init__(
        self,
        path: StrPath,
        *,
        chr_sizes: Mapping[str, int] | None = None,
        genome: str | None = None,
        items_per_slot: int | None = None,
        compression_level: int = 6,
        layout: Literal["optimal", "as_given"] = "optimal",
        zoom_share: float = 0.2,
        parallel: int | None = None,
    ) -> None: ...
    def write_value(self, chr: str, start: int, end: int, value: float) -> None: ...
    def write_values(self, chr: str, start: int, span: int, values: npt.ArrayLike) -> None: ...
    def close(self, *, progress: Progress = None) -> None: ...
    section_count: int
    section_counts: dict[str, int]
    skipped_count: int
    clipped_count: int
    item_count: int
```

- `chr_sizes`, or an assembly's through `genome`, bound every coordinate, and a
  name written is matched against them as a read one is. Without them a
  chromosome ends where its last value does. Only chromosomes written go into
  the file.
- Values must be pooled by chromosome and written in order without overlap.
  `write_values` writes a run of values `span` bases each, the first at
  `start`, which is stored in a third of what the same values cost one at a
  time. A NaN or infinite value is not written (`skipped_count`): it leaves a
  gap. A value hanging over the end of a declared chromosome is cut back to it
  (`clipped_count`); one starting past it raises.
- `items_per_slot` is the values a block holds (1024, as UCSC's tools);
  `compression_level` 0 to 9; `layout="optimal"` re-cuts the values into the
  sections that store them in the fewest bytes, `"as_given"` writes one item per
  value handed over. `zoom_share` is the most of the data section the finest
  zoom level may take: smaller is a smaller file whose finest summary is
  coarser.
- `close(progress=)` builds the zoom levels and the index; its progress counts
  the bytes the zoom levels read back.

```pycon
>>> writer = gwseq_io.BigWigWriter("draft.bw", chr_sizes={"chr1": 1000})
>>> writer.write_value("chr1", 0, 10, 1.0)
>>> writer.write_values("chr1", 10, 5, [1.0, 2.0, float("nan")])
>>> writer.item_count, writer.skipped_count, writer.state
(3, 1, 'open')
>>> writer.abandon()
>>> writer.state, writer.closed
('closed', True)
```

#### BigBedWriter

```pyi
class BigBedWriter:
    def __init__(
        self,
        path: StrPath,
        *,
        fields: Mapping[str, str] | None = None,
        chr_sizes: Mapping[str, int] | None = None,
        genome: str | None = None,
        items_per_slot: int | None = None,
        compression_level: int = 6,
        zoom_share: float = 0.2,
        extra_index: Sequence[str] = (),
        parallel: int | None = None,
    ) -> None: ...
    def write_entry(
        self, chr: str, start: int, end: int, fields: Mapping[str, Any] | None = None
    ) -> None: ...
    def close(self, *, progress: Progress = None) -> None: ...
    entry_count: int
    fields: dict[str, str]
```

- `fields` names the columns and their autoSql types, the three coordinates
  first; by default `chr`, `start`, `end` and `name`. `"string"`, `"int"`,
  `"uint"` and `"float"` are the usual ones; the rest of autoSql (`lstring`,
  `char[n]`, `double`, arrays, `enum`, `set`) is accepted so a file read here
  can be written back as it was.
- `write_entry`'s `fields` are the columns past the coordinates; one left out
  is written empty, or as zero when numeric, and each value is checked against
  its column's type. Entries must be pooled by chromosome and in order of their
  start, and may overlap and nest.
- `extra_index` names columns to index by their text, for `search`.
- The summary and zoom levels describe the depth of coverage the entries make.

```pycon
>>> fields = {"chr": "string", "start": "uint", "end": "uint", "name": "string", "score": "float"}
>>> with gwseq_io.open("scored.bb", "w", fields=fields, chr_sizes={"chr1": 1000}) as writer:
...     writer.write_entry("chr1", 10, 20, {"name": "a", "score": 1.5})
...     writer.write_entry("chr1", 15, 30, {"name": "b"})
>>> gwseq_io.open("scored.bb").read_all_entries(as_cols=True)
{'chr': ['chr1', 'chr1'], 'start': array([10, 15]), 'end': array([20, 30]), 'name': ['a', 'b'], 'score': array([1.5, 0. ])}
```

#### BamWriter

```pyi
class BamWriter:
    def __init__(
        self,
        path: StrPath,
        *,
        header: str | Sequence[Mapping[str, Any]] | BamReader | CramReader,
        index: bool = True,
        compression_level: int = 6,
        parallel: int | None = None,
    ) -> None: ...
    def write(self, entry: SamEntry) -> None: ...
    def close(self) -> None: ...
    record_count: int
    index_path: str | None
```

- `header` is SAM header text, a reader's `header` (the list of dicts), or a
  `BamReader` or `CramReader` itself, whose header is copied verbatim.
- `write` takes a `SamEntry` as a reader gave it, with every field its record
  holds: an entry read with `tags=False` is written without its optional
  fields, and one read from a CRAM with `sequence=False` without its bases. An
  entry of a file whose references are numbered otherwise is placed on the
  header's reference of the same name; a name the header lacks raises
  `InvalidRequest`, with nothing written.
- With `index=True` (the default) the entries must come sorted by coordinate:
  one out of order raises `InvalidRequest`, naming both reads, and the writer
  stays open. `close()` writes `<path>.bai`, or `<path>.csi` when a reference is
  longer than 512 Mb (`index_path`).
- **Existing indexes are removed.** `close()` deletes an index already beside
  the path (`<path>.bai`, `<path>.csi`, and `<stem>.bai` for a path ending
  `.bam`), since it indexes the file being replaced, whether or not a new one
  is written. A failed or abandoned writer deletes nothing.

```pycon
>>> reader = gwseq_io.open("reads.bam")
>>> with gwseq_io.open("subset.bam", "w", header=reader) as writer:
...     for entry in reader.read_entries("chr1:1,001-2,300")[0]:
...         writer.write(entry)
>>> writer.record_count, writer.index_path
(4, 'subset.bam.bai')
>>> [entry.read_name for entry in gwseq_io.open("subset.bam").read_all_entries()]
['plain_fwd', 'plain_rev', 'pair_fr', 'pair_fr']
```

### Converters

```pyi
def convert_to_bigwig(
    input: StrPath,
    output: StrPath,
    *,
    bin_size: int | None = None,
    chr_sizes: Mapping[str, int] | None = None,
    genome: str | None = None,
    items_per_slot: int | None = None,
    compression_level: int = 6,
    zoom_share: float = 0.2,
    parallel: int | None = None,
    progress: Progress = None,
) -> dict[str, Any]: ...
def convert_to_bigbed(
    input: StrPath,
    output: StrPath,
    *,
    fields: Mapping[str, str] | None = None,
    chr_sizes: Mapping[str, int] | None = None,
    genome: str | None = None,
    items_per_slot: int | None = None,
    compression_level: int = 6,
    zoom_share: float = 0.2,
    extra_index: Sequence[str] = (),
    parallel: int | None = None,
    progress: Progress = None,
) -> dict[str, Any]: ...
def convert_to_bam(
    input: StrPath,
    output: StrPath,
    *,
    reference: StrPath | None = None,
    regions: RegionsLike | None = None,
    exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
    require_flags: Flags = 0,
    min_mapq: int = 0,
    drop_improper_pairs: bool = True,
    index: bool = True,
    compression_level: int = 6,
    parallel: int | None = None,
    progress: Progress = None,
) -> dict[str, Any]: ...
```

Each converter streams its input, which may be gzipped where it is text, and
tells its format from its content, not its name. The output holds the finished
file or what it held before, never part of one. Each returns a dict of what was
read and written. The options they share with a writer are that writer's.

- **`convert_to_bigwig`**: bedGraph or WIG (the first line that is not blank, a
  comment, `track` or `browser` decides). Values must be pooled by chromosome
  and in order without overlap, as `sort -k1,1 -k2,2n` gives. `bin_size`
  averages the values into bins of that width, each weighted by the bases it
  covers.
- **`convert_to_bigbed`**: BED, split on tabs, every record with the first's
  number of columns. `fields` names the columns; by default the standard BED
  names and types, and naming only some standard columns changes their types
  and keeps the rest (`fields={"score": "float"}`). Standard columns are checked
  as UCSC checks them. `extra_index` is `bedToBigBed -extraIndex`.
- **`convert_to_bam`**: SAM (gzipped or not), BAM or CRAM, with `reference` for
  a CRAM as `CramReader` takes it; a CRAM holding a sequence no reference is
  found for is refused before anything is written. `regions` keeps the records
  overlapping any of them, each once, in the input's order; a BAM or CRAM input
  needs its index for it. The read filter applies, with the defaults above:
  pass `exclude_flags=0, drop_improper_pairs=False` to keep every record.
  `index` writes `<output>.bai` (or `.csi`) as `BamWriter` does, and the
  input must then be sorted; as with `BamWriter`, an index already beside
  `output` is removed. Progress is counted in bases of the input's references,
  or bytes of a SAM.

```pycon
>>> with open("track.bedgraph", "w") as text:
...     _ = text.write("chr1\t0\t100\t0.5\nchr1\t100\t300\t1.5\nchr2\t0\t50\t2\n")
>>> gwseq_io.convert_to_bigwig("track.bedgraph", "track.bw", chr_sizes={"chr1": 1000, "chr2": 500})
{'format': 'bedgraph', 'line_count': 3, 'item_count': 3, 'skipped_count': 0, 'clipped_count': 0, 'chr_sizes': {'chr1': 1000, 'chr2': 500}}
>>> with open("peaks.bed", "w") as text:
...     _ = text.write("chr1\t100\t200\tpeak1\t500\t+\nchr1\t150\t400\tpeak2\t900\t-\n")
>>> gwseq_io.convert_to_bigbed("peaks.bed", "peaks.bb", genome="mm10", extra_index=["name"])["item_count"]
2
>>> gwseq_io.open("peaks.bb").search("peak2")
[BedEntry(chr1:150-400, 3 columns)]
>>> gwseq_io.convert_to_bam("reads.cram", "from_cram.bam", reference="ref.fa", regions="chr1:1-2,000")
{'format': 'cram', 'record_count': 90, 'index': 'from_cram.bam.bai', 'index_kind': 'bai', 'chr_sizes': {'chr1': 34000, 'chr2': 8000, 'chr3': 3000}}
>>> gwseq_io.convert_to_bam("reads.bam", "every_record.bam", exclude_flags=0, drop_improper_pairs=False)["record_count"]
60
```

### Pileup

```pyi
def pileup_reads(
    input: StrPath,
    output: StrPath,
    *,
    bin_size: int = 10,
    exclude_flags: Flags = "UNMAP,SECONDARY,SUPPLEMENTARY",
    require_flags: Flags = 0,
    min_mapq: int = 0,
    drop_improper_pairs: bool = True,
    strand: Literal["+", "-"] | None = None,
    max_fragment: int = 2000,
    shift_reads: int = 0,
    transform_reads: Literal["extend", "force_extend", "splice"] = "extend",
    extend_size: int = 150,
    count_fragments_once: bool = False,
    def_value: float | None = 0.0,
    normalize: Literal["count", "cpm", "bpm", "enrichment"] = "count",
    scale: float = 1.0,
    index: StrPath | None = None,
    reference: StrPath | None = None,
    hold_bytes: int | None = None,
    items_per_slot: int | None = None,
    compression_level: int = 6,
    zoom_share: float = 0.2,
    parallel: int | None = None,
    progress: Progress = None,
) -> dict[str, Any]: ...
```

Piles a BAM or CRAM up into a bigWig of `bin_size` bins, each the number of
intervals that touch it. The input must be sorted by coordinate (an unsorted
one is refused, naming `samtools sort`); an index is used when there is one and
is not needed. `index` names it when it is not beside the input, and
`reference` is a CRAM's FASTA.

A read becomes an interval in three steps:

1. **The filter** (the [read filter](#bam-and-cram), same defaults): unmapped,
   secondary and supplementary alignments and improper pairs are dropped;
   duplicates and QC failures are kept.
2. **The strand**, with `strand="+"` or `"-"`: a read is kept when its
   fragment's strand is that one. A fragment's strand is the read's own, flipped
   for the second read of a pair (0x80 set, 0x40 not), so a pair counts on the
   strand of its first read and a single read on its own. **dUTP
   (fr-firststrand) libraries sequence the complement, so there the labels are
   swapped**: `"+"` is the track of the genes on the minus strand.
3. **The transform**, `transform_reads`:
   - `"extend"`: a read of a proper pair with its mate on the same reference,
     a template length and exactly one of 0x40/0x80 is piled up as its
     **fragment**, from the leftmost mate's start over `|TLEN|` bases; a
     fragment longer than `max_fragment`, or starting more than `max_fragment`
     before its read, is dropped (its read still counted).
     Every other read is extended from its 5′ end to `extend_size` bases, or its
     own aligned length at 0. With `count_fragments_once`, a fragment is piled
     up once, by the mate with the lower position (READ1 on a tie), rather than
     once per mate.
   - `"force_extend"`: every read is moved `shift_reads` bases along its own
     direction and extended from its 5′ end to `extend_size`, pairs or not.
     `shift_reads` is refused with `"extend"`.
   - `"splice"`: the segments of each read between the `N` operations of its
     CIGAR (its introns), each moved `shift_reads` bases along the read.

`normalize` is `"count"` (none), `"cpm"` (per million intervals), `"bpm"` (each
bin as its share of the total, per million) or `"enrichment"` (over the
geometric mean of the covered bins); `scale` multiplies after it. Every
normalization but `"count"` is a property of the whole file: the binned result
is held in memory, or the input read twice when it would take more than
`hold_bytes` (1 GiB when `None`; `passes` in the result says which).
`def_value` is written for a bin with no coverage, whatever the normalization;
`None` leaves those bins out of the file. Progress is counted in bytes of the
input. The result's `read_count` is what the filter and the strand kept, and
`interval_count` what was piled up.

```pycon
>>> summary = gwseq_io.pileup_reads("reads.bam", "coverage.bw")
>>> summary["read_count"], summary["interval_count"], summary["passes"]
(43, 40, 1)
>>> gwseq_io.open("coverage.bw").read_values("chr1:2,001-2,300", bin_size=50)
array([[2., 2., 2., 2., 2., 0.]], dtype=float32)
>>> gwseq_io.pileup_reads("reads.bam", "coverage.cpm.bw", normalize="cpm")["scale"]
25000.0
>>> gwseq_io.pileup_reads("reads.bam", "coverage.once.bw", count_fragments_once=True)["interval_count"]
33
>>> gwseq_io.pileup_reads("reads.cram", "cram.bw", reference="ref.fa", def_value=None)["read_count"]
682
```

Recipes: ChIP `normalize="cpm"`; ATAC `shift_reads=-36,
transform_reads="force_extend", extend_size=80, normalize="cpm"`; RNA-seq
`bin_size=1, transform_reads="splice", normalize="cpm"`.

### Walks, regions and headers

`RegionIterator` is what every `iter_*` method returns: one step per region or
window. It is re-iterable (each `for` loop reads again), and an interrupted pass
does not spoil the next. `len()` is the number of steps, `regions[n]` the region
of step `n`, and `items()` pairs each step with its region. `order` is, for
`iter_entries`, the index in the request of each step's region, and `None` for a
walk over whole chromosomes. A walk with `unplaced=True` ends in batches with no
region: `len()` raises `TypeError` and `items()` pairs them with `None`.
`Regions` is a sequence of `Region` with `to_dict()`, the columns
`pandas.DataFrame` takes.

```pyi
class RegionIterator(Generic[T]):
    regions: Regions
    order: npt.NDArray[np.int64] | None
    def items(self) -> Iterator[tuple[Region | None, T]]: ...

class Regions:
    def index(self, value: object) -> int: ...
    def count(self, value: object) -> int: ...
    def to_dict(self) -> dict[str, Any]: ...

class BbiHeader:
    magic: int
    version: int
    zoom_levels: int
    chr_tree_offset: int
    full_data_offset: int
    full_index_offset: int
    field_count: int
    defined_field_count: int
    auto_sql_offset: int
    total_summary_offset: int
    uncompress_buffer_size: int

class ZoomHeader:
    reduction_level: int
    data_offset: int
    index_offset: int

class TotalSummary:
    bases_covered: int
    min_value: float
    max_value: float
    sum_data: float
    sum_squared: float
```

```pycon
>>> walk = gwseq_io.open("features.bb").iter_all_entries(chromosomes=["chrX"], window=30_000)
>>> len(walk), [len(entries) for entries in walk]
(3, [16, 15, 9])
>>> walk.regions.to_dict()["start"]
array([    0, 30000, 60000])
>>> summary = gwseq_io.open("signal.bw").total_summary
>>> summary.bases_covered, summary.min_value, summary.max_value, summary.sum_data
(400, 0.5, 2.0, 500.0)
```


## The command line

`gwseq-io` is installed with the package, and is the same program as
`python -m gwseq_io`. Data goes to standard output as tab-separated text;
summaries, progress and errors go to standard error. An error is one line,
`gwseq-io: <message>`.

| Command | What |
|---|---|
| `info` | A file's format, headers, chromosomes and summary |
| `read` | bigWig or bigBed values over regions: one row per region, one column per bin; every chromosome, one at a time, with no region |
| `quantify` | One value per region |
| `profile` | One value per bin, reduced across regions |
| `entries` | bigBed entries, or BAM/CRAM alignments (`name chr start end flag mapq cigar`); the whole file, with no index needed, when no region is given |
| `search` | The bigBed entries whose indexed column holds one of the values |
| `hic` | A HiC contact matrix over one or two regions |
| `export` | A bigWig or bigBed to bedGraph, WIG or BED |
| `convert` | bedGraph, WIG or BED to bigWig or bigBed; SAM, BAM or CRAM to BAM, by the output's extension |
| `pileup` | A BAM or CRAM to a bigWig (`pileup-reads` and `pileup_reads` are the same command) |
| `chr-sizes` | An assembly's chromosome sizes |

- Global options: `-p/--parallel N`, the most threads one command uses;
  `-q/--quiet`, nothing but errors; `--progress`, a progress bar on standard
  error, drawn by default when it is a terminal. `gwseq-io <command> --help` lists
  a command's options.
- Regions are `-r chr:start-end`, 1-based and inclusive, or a chromosome's
  name, repeatable. `--zoom` is `full`, `auto` or a level's index. Flags are
  samtools names or numbers.
- The flags of `gwseq-io pileup` are `pileup_reads`' keywords, hyphenated
  (`--exclude-flags`, `--count-fragments-once`, `-d none` for
  `def_value=None`), plus `--range-bytes`, the compressed bytes read as one
  range. `-b`, `-s`, `-t`, `-e`, `-d`, `-n` and `-T` are the short forms of
  `--bin-size`, `--shift-reads`, `--transform-reads`, `--extend-size`,
  `--def-value`, `--normalize` and `--reference`.
- The filter of `entries`, `convert` and `pileup` is the read filter, with its
  defaults; `--keep-improper-pairs` is `--drop-improper-pairs false`, and
  `convert --no-filter` keeps every record.
- Exit codes: 0 success, 1 an error, 2 a usage error, 130 interrupted (Ctrl-C,
  which leaves no partial output). A closed pipe on standard output is a
  success.

```console
$ gwseq-io info signal.bw
path           signal.bw
type           bigwig
version        4
chromosomes    1
genome size    10000 bp
data           45 bytes
zoom levels
summary        400 bases covered, min 0.5, max 2, mean 1.25

chr1    10000
$ gwseq-io entries reads.bam -r chr1:1,001-2,300
plain_fwd    chr1    1000    1050    0      60    50M
plain_rev    chr1    1100    1150    16     30    50M
pair_fr      chr1    2000    2050    99     60    50M
pair_fr      chr1    2200    2250    147    60    50M
$ gwseq-io search features.bb LOC472203_1
chr1    1336    3623    LOC472203_1    47    -    205.43    n/a
$ gwseq-io hic contacts.hic chr1:1-300,000 --bin-size 100000
367    96     50
96     414    113
50     113    428
$ gwseq-io pileup reads.bam coverage.cpm.bw -n cpm
bin size       10
reads          43
intervals      40
normalization  cpm x25000 (1 pass)
covered bins   367
items written  59
sections       9
chromosomes    3
$ gwseq-io convert reads.cram from_cram.bam -T ref.fa -r chr1:1-2,000
format           cram
records written  90
references       3
index            from_cram.bam.bai
$ gwseq-io export signal.bw signal.bedgraph --bin-size 200
$ gwseq-io chr-sizes sacCer3
chrI       230218
chrII      813184
...
$ gwseq-io read missing.bw
gwseq-io: error reading missing.bw: No such file or directory (os error 2)
```

Recipes:

```console
$ gwseq-io pileup reads.bam chip.cpm.bw -n cpm -q
$ gwseq-io pileup reads.bam atac.cpm.bw -s -36 -t force_extend -e 80 -n cpm -q
$ gwseq-io pileup reads.bam rna.cpm.bw -b 1 -t splice -n cpm -q
```

`cargo install --path crates/gwseq-io-cli` builds `gwseq-io` alone, as a Rust
binary.


## Development

The workspace holds three crates: `gwseq-io`, the library, with no Python in
it (its own README covers the Rust API); `gwseq-io-cli`, the command line, a
library and a binary; and `gwseq-io-py`, the extension, which links both and
holds the Python package (`crates/gwseq-io-py/python/gwseq_io`) with its
hand-written stub `_gwseq_io.pyi`. [docs/architecture.md](docs/architecture.md) is the
design. `docs/` holds the format specifications implemented.

**Building.** `pip install -e ".[dev]"` in a venv builds and installs the
package with what building, checking and publishing need; `maturin develop
--release` rebuilds only what changed. `cargo build --workspace` does not link
the extension (maturin does): use `cargo check -p gwseq-io-py` for its type
errors. The Rust floor is `rust-version` in `Cargo.toml`, 1.85.

**The check list.** No CI runs anywhere; tests run on this machine through one
script, standard library only:

- `python3 tools/check.py --quick`: `cargo fmt --check`, clippy (all features,
  and the core with none), the debug tests, `cargo check -p gwseq-io-py`.
- `python3 tools/check.py`: all of that, then the release tests, the ignored
  tests, the codec oracle, `cargo doc` with warnings as errors, the 1.85
  toolchain and the Windows and macOS targets when installed, and **the Python
  step**: a release wheel built with `maturin build` into the cargo target
  directory, installed into a throwaway venv there with the `test` extra
  (pytest, mypy, pandas, scipy), and the pytest suite run against it.
  `--python PATH` also runs the suite under another interpreter; the floor
  check is `--python python3.11`, which gets the oldest numpy with wheels for
  it (1.23.2).
- `python3 tools/dist.py --check`: `tools/check.py`, then the sdist built,
  compiled in a clean venv and tested there with the pytest suite.
  `--check-versions` checks that the workspace version, the `gwseq-io` pin and
  `pyproject.toml` agree. `tools/dist.py` also uploads (`--upload`,
  `--test-upload`, `--crate-upload`).

**Tests.**

- `cargo test --workspace`: unit tests in every module, the integration suite
  (`crates/gwseq-io/tests/suite/`: round trips, requests against independent
  oracles, a fuzz pass over every parser, properties), and the command line's
  tests. No network, samtools, Java or C compiler is needed.
- `crates/gwseq-io-py/tests/`: the pytest suite, run against an installed
  build, never the source tree. `conftest.py` copies the committed fixtures
  into each test's own directory. `test_readme.py` runs this file: each
  ```` ```pycon ```` block through doctest, each ```` ```console ```` command
  through the installed `gwseq-io`, each ```` ```pyi ```` block against the stub;
  any other fence but ```` ```text ```` fails. The stub is held to the
  extension by `mypy.stubtest` and `tests/typing/calls.py` by `mypy --strict`.
- Tests that need the large files of `local/test_data/` (not committed) skip
  without them; so do the benchmarks (`cargo bench -p gwseq-io`).

**Fixtures** are committed under `crates/gwseq-io/tests/data/`, excluded from
both packages, and written by deterministic generators run by hand, each with a
`--check` that regenerates and compares:

- `tools/fixtures/make_alignment_fixtures.sh`: the BAMs, indexes and samtools
  views, from hand-written SAM files (samtools 1.23);
- `tools/fixtures/cram_fixtures.py`: CRAM 2.1, 3.0 and 3.1 files in every
  samtools profile, and what samtools reads back from them;
- `tools/fixtures/bbi/make_extra_index.sh`: bigBeds written by UCSC's
  `bedToBigBed` (pinned, sha256-checked download);
- `tools/fixtures/hic/make.sh`: `.hic` files written by juicer_tools 1.22 and
  2.20, and what straw reads out of them.

The examples of this README use `geometry.bam`, `v31_normal_default.cram` with
`ref.fa`, `extra_index.bb`, `bed12_ucsc.bb` and `v9.hic` from there, copied
under the names above, and write `signal.bw` as the quick start does.

**Oracles.** `oracle/` is a crate outside the workspace, the one place that
compiles C: it vendors htscodecs and checks the CRAM codecs against it
(`cd oracle && cargo test`, which `tools/check.py` runs when a C compiler is
installed). `tools/oracle/pileup_oracle.py` computes a pileup from
`samtools view` text, sharing no code with the library.

