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, 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

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

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

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

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.

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

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.

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

Module functions

@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. 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.
• 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.

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

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).
• 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

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

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

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

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

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

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.

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

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

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

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): 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").

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

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

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

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.

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

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):
  • 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.

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

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.

>>> 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)]

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.

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

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.

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

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.

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

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.

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

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.

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

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

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

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

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

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

$ 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:

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