Metadata-Version: 2.4
Name: gwseq_io
Version: 0.2.0
Requires-Dist: numpy
Requires-Dist: build ; extra == 'dev'
Requires-Dist: twine ; extra == 'dev'
Requires-Dist: maturin>=1.7,<2.0 ; extra == 'dev'
Provides-Extra: dev
License-File: LICENSE
Summary: Python library for processing bigWig, bigBed, BAM and HiC files (Rust)
Author-email: Arthur Gouhier <ajgouhier@gmail.com>
License-Expression: BSD-3-Clause
Requires-Python: >=3.9
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Repository, https://github.com/ajgouhier/gwseq_io

# gwseq_io

Python library for processing bigWig, bigBed, BAM and HiC files.
Backed by a Rust core via PyO3.


## Installation

```
pip install gwseq-io
```

Requires numpy, installed automatically as a dependency.

Only a source distribution is published, so pip builds the extension on the
installing machine. That needs a Rust toolchain, 1.85 or newer.


## Usage

### Open bigWig, bigBed, BAM and HiC files for reading

```python
reader = gwseq_io.open(path, mode, parallel, zoom_correction, file_buffer_size, max_file_buffer_count, index_path)

with gwseq_io.open("path/to/file.bigwig") as reader: # .bigbed .bam .hic
    ...
```

Parameters:
- `mode` Opening mode. May be omitted as "r" (read) by default.
- `parallel` Number of parallel file handles and processing threads. Use -1 for recommended (one per core, capped at 12). -1 by default.
- `zoom_correction` Scaling factor for automatic zoom level selection based on bin size. Only for bigWig files. 1/3 by default.
- `file_buffer_size` Size in bytes of each file buffer for caching file reads. Use -1 for recommended (32768 or 1048576 for URLs). -1 by default.
- `max_file_buffer_count` Maximum number of file buffers to keep in cache. Use -1 for recommended (128). -1 by default.
- `index_path` Path of the index. Only for BAM files, where it defaults to the path of the file with ".bai" appended — including for a URL, whose index is fetched from `<url>.bai` over the same connection. An index is optional, but reading entries needs one.

Common attributes and methods:
- `close` Give back the `parallel` file handles and the threads the reader holds for as long as it lives. Calling it twice is harmless, and a reader is a context manager, so leaving the `with` block above closes it. Reading through a closed reader raises, and the headers it read at open stay readable. A reader that is never closed gives everything back when it is collected instead.
- `closed` Whether `close` has run.

Attributes for bigWig and bigBed files:
- `main_header` General file formatting info.
- `zoom_headers` Zooms levels info (reduction level and location).
- `auto_sql` BED entries declaration (only in bigBed).
- `total_summary` Statistical summary of entire file values (coverage, sums and extremes).
- `chr_sizes` Map of chromosome IDs and their sizes.
- `type` Either "bigwig" or "bigbed".

Attributes for BAM files:
- `header` Header lines, each a dict of its "type" (the two letters after the @) and its "fields".
- `chr_sizes` Map of reference IDs and their sizes.
- `is_indexed` Whether the index was found and read. Reading entries needs it.
- `index_error` Why the index is absent, when it is. Empty when it loaded, and empty as well when the file simply has none.

Attributes for HiC files:
- `header` `footer` General file info.
- `chr_sizes` Map of chromosome IDs and their sizes.
- `normalizations` Available normalizations.
- `units` Available units.
- `bin_sizes` Available bin sizes.

### Read bigWig and bigBed values

```python
values = reader.read_values(chr_ids, starts, ends, centers, span, ...)

values = reader.read_values(chr_ids=["chr1", "chr1"], starts=[1000, 1100], ends=[1100, 1200])
values = reader.read_values(chr_ids=["chr1", "chr1"], starts=[1000, 1100], span=100)
values = reader.read_values(chr_ids=["chr1", "chr1"], ends=[1100, 1200], span=100)
values = reader.read_values(chr_ids=["chr1", "chr1"], centers=[1050, 1150], span=100)
values = reader.read_values(chr_ids=["chr1", "chr1"], centers=[1050, 1150], span=100, strands=["+", "-"])
```

Parameters:
- `chr_ids` `starts` `ends` `centers` Chromosome IDs, starts, ends and centers of the locations. Both `starts` `ends`, or one of `starts` `ends` `centers` with `span`, may be specified.
- `span` Reading window in bp relative to `starts`, `ends` or `centers`. Only one of the three may be given with it. Not by default.
- `strands` Strand of each location, as "+" or "-" ("." and "" count as "+"). The values of a "-" location are reversed, so that every location reads from its own start. All "+" by default.
- `bin_size` Reading bin size in bp. May vary in output if locations have variable spans or `bin_count` is specified. 1 by default.
- `bin_count` Output bin count. Inferred as max location span / bin size by default.
- `bin_mode` Method to aggregate bin values, all three per base of the bin rather than per record of the file: "mean" is the base-weighted mean, "sum" the value summed over each base it covers, and "count" the bases of the bin carrying data — the bin's width where the file covers it fully. "mean" by default. For a bigBed the value of a bin is the depth of coverage its entries make over that bin.
- `full_bin` Extend locations ends to overlapping bins if true. Not by default.
- `def_value` Default value to use when no data overlap a bin. 0 by default.
- `zoom` BigWig zoom level to use. Use full data if -1, or auto-detect if -2 by taking the coarsest level whose bin size is under `bin_size` times `zoom_correction` (may be the full data). Full data by default.
- `progress` Function called during extraction with the extracted and the total coverage in bp. Use the default callback if true. None by default.

Returns a numpy float32 array of shape (locations, bin count).

### Quantify bigWig and bigBed values

```python
values = reader.quantify(chr_ids, starts, ends, centers, span, ...)
```

Parameters:
- `chr_ids` `starts` `ends` `centers` `span` `bin_size` `full_bin` `def_value` `zoom` `progress` Identical to `read_values` method.
- `reduce` Method to aggregate values over span. Either "mean", "sd", "sem", "sum", "count", "min", "max", "l1norm" or "l2norm". "mean" by default.

Notes:
- Read through a zoom level (`zoom` other than -1), `l1norm` is a **lower bound** rather than the exact `Σ|x|`. A zoom record stores a min, a max, a sum and a sum of squares, and `Σ|x|` cannot be recovered from those in general — so it is exact where the record's sign is not in doubt (`min >= 0`, or `max <= 0`, which is every non-negative track) and `|Σx|` where the record straddles zero. Every other reduction is exact at every level.

Returns a numpy float32 array of shape (locations).

### Profile bigWig and bigBed values

```python
values = reader.profile(chr_ids, starts, ends, centers, span, ...)
```

Parameters:
- `chr_ids` `starts` `ends` `centers` `span` `strands` `bin_size` `bin_count` `bin_mode` `full_bin` `def_value` `zoom` `progress` Identical to `read_values` method. A "-" location takes part in the profile reversed, as it would come out of `read_values`.
- `reduce` Method to aggregate values over locations. Either "mean", "sd", "sem", "sum", "count", "min", "max", "l1norm" or "l2norm". "mean" by default.

Returns a numpy float32 array of shape (bin count).

### Iterate over all bigWig and bigBed values

```python
iterator = reader.iter_all_values(...)

iterator = reader.iter_all_values(bin_size=10)
for values in iterator:
    ...
for (chr_id, start, end), values in zip(iterator.locs, iterator):
    ...
```

Parameters:
- `chr_ids` Only walk these chromosomes. All by default.
- `bin_mode` `full_bin` `def_value` `zoom` `progress` Identical to `read_values` method. `full_bin` decides whether the partial bin a chromosome ends on is walked at all.
- `span` Window in bp for each step. 1,000,000 by default.

Returns an iterator over successive windows, each one a numpy float32 array of shape (bins), in chromosome then coordinate order. `len(iterator)` gives the number of windows, and `iterator.locs` the region of each, so the nth array covers `locs[n]`:

Notes:
- A window never spans two chromosomes and no bin straddles a window boundary, so concatenating the windows of a chromosome gives exactly what `read_values` gives for the whole of it at the same bin size. For a bigBed the values are the pileup of its entries.
- An iterator may be walked more than once. `__iter__` hands back a fresh cursor over the same plan, so a second `for` loop reads the file again and `zip(iterator.locs, iterator)` works every time. The plan is shared rather than copied, so an extra pass costs a little over a hundred bytes and the reads it makes. The same holds for every `iter_*` method of every reader.

### Read bigBed entries

```python
entries = reader.read_entries(chr_ids, starts, ends, centers, span, ...)
```

Parameters:
- `chr_ids` `starts` `ends` `centers` `span` `progress` Identical to `read_values` method.
- `col_count` Only read this number of columns (eg, 3 for chr, start and end). Must be 0 (all) or at least 3. The columns left out are never parsed, so a narrower read is a cheaper one. All by default.

Returns a list (locations) of list of entries (dict with at least "chr", "start" and "end" keys).

### Read all bigBed entries

```python
entries = reader.read_all_entries(...)
```

Parameters:
- `chr_ids` Only extract data from these chromosomes. All by default.
- `col_count` Identical to `read_entries` method.

Returns a list of entries (as in `read_entries`).

### Iterate over all bigBed entries

```python
iterator = reader.iter_all_entries(...)

iterator = reader.iter_all_entries()
for entries in iterator:
    ...
for (chr_id, start, end), entries in zip(iterator.locs, iterator):
    ...
```

Parameters:
- `chr_ids` `col_count` `progress` Identical to `read_all_entries` method.
- `span` Identical to `iter_all_values` method.

Returns an iterator over successive windows, each one a list of entries (as in `read_entries`), in chromosome then coordinate order. `len(iterator)` gives the number of windows, and `iterator.locs` the region of each.

Notes:
- A window never spans two chromosomes, and an entry reaching over a window boundary is reported by the window it starts in. Concatenating the windows gives exactly what `read_all_entries` returns, in the same order.

### Convert bigWig to bedGraph or WIG

```python
reader.to_bedgraph(output_path, ...)
reader.to_wig(output_path, ...)
```

Parameters:
- `output_path` Path to output file.
- `chr_ids` Only extract data from these chromosomes. All by default.
- `bin_size` `bin_mode` `full_bin` `def_value` `zoom` `progress` Identical to `read_values` method.
- `merge_bins` Whether adjacent bins with the same value are merged into one interval. Only for bedgraph output. True by default.

Notes:
- The values written are the ones `read_values` gives for the same chromosomes at the same settings, so a default export writes one interval per base and `bin_size=10000` writes one per 10,000 bases. `full_bin` decides whether the shorter last bin a chromosome ends on is written at all, as it does in `iter_all_values`.
- A bin no data reaches holds `def_value`, so the gaps of a bigWig come out as intervals of 0 by default. `def_value=float("nan")` is how the covered part alone is asked for: a bin holding NaN is left out of the file, which is also what keeps a NaN out of text no reader of either format would accept.
- A fixedStep WIG section carries one value per line and has no way to say "and again", which is why `merge_bins` is bedGraph's alone. `to_wig` opens a new section wherever the run of bins breaks: a change of chromosome, a bin left out, and the shorter last bin of a `full_bin` export.

### Convert bigBed to BED

```python
reader.to_bed(output_path, ...)
```

Parameters:
- `output_path` `chr_ids` `progress` Identical to `to_bedgraph` and `to_wig` methods.
- `col_count` Only write this number of columns (eg, 3 for chr, start and end). All by default.

### Read BAM entries

```python
entries = reader.read_entries(chr_ids, starts, ends, centers, span, ...)
```

Parameters:
- `chr_ids` `starts` `ends` `centers` `span` `progress` Identical to bigWig `read_values` method.
- `filter` Drop unmapped alignments, improperly paired reads, secondary and supplementary records, and anything marked as failing quality control or as a duplicate. True by default.
- `parse_tags` Keep the optional fields of every alignment. They are only decoded when read, so this costs one small copy per alignment. True by default.

Returns a list (locations) of list of `BamEntry`, in the order the locations were given. Each location gets its own full list, so two overlapping locations both report the alignments they share.

Notes:
- Needs an index, as `is_indexed` reports. Reading without one raises, naming either the index that was not found or the error it gave.

### Read all BAM entries

```python
entries = reader.read_all_entries(...)
```

Parameters:
- `chr_ids` Only extract data from these references. All by default.
- `filter` `parse_tags` `progress` Identical to `read_entries` method.

Returns a list of `BamEntry` (as in `read_entries`). Unplaced alignments are left out, the index reaching an alignment only through the reference it sits on.

### Iterate over BAM entries

```python
iterator = reader.iter_entries(chr_ids, starts, ends, centers, span, ...)

iterator = reader.iter_entries(chr_ids=["chr1", "chr1"], starts=[1000, 1100], ends=[1100, 1200])
for entries in iterator:
    ...
for loc_index, entries in zip(iterator.order, iterator):
    ...
```

Parameters:
- `chr_ids` `starts` `ends` `centers` `span` `filter` `parse_tags` `progress` Identical to `read_entries` method.
- `sort_locations` Read the locations in reference and position order, and report them in that order. Can be more efficient. False by default.

Returns an iterator over locations, each one a list of `BamEntry`. `len(iterator)` gives the number of locations, and `iterator.order` the request index of each, so the nth list belongs to location `order[n]`:

### Iterate over all BAM entries

```python
iterator = reader.iter_all_entries(...)

for entries in reader.iter_all_entries():
    ...
```

Parameters:
- `chr_ids` `filter` `parse_tags` `progress` Identical to `read_all_entries` method.
- `span` Window in bp for each step. 1,000,000 by default.

Returns an iterator over successive windows, each one a list of `BamEntry`, in reference then coordinate order. `len(iterator)` gives the number of windows.

Notes:
- A window never spans two references, and an alignment reaching over a window boundary is reported by the window it starts in. Concatenating the windows gives exactly what `read_all_entries` returns, in the same order.

### BAM entries

A `BamEntry` is one alignment.

Attributes:
- `chr` (str) Reference the alignment sits on.
- `start` `end` (int) 0-based half-open span, the end derived from the cigar. Equal for an alignment covering no reference.
- `read_name` (str) QNAME.
- `flag` (int) FLAG, as the raw bitfield.
- `mapping_quality` (int) MAPQ.
- `cigar` (str) CIGAR, eg "10S80M10S". A cigar of more than 65 535 operations does not fit the record's own field, so it is stored in a `CG` optional field and read from there; the placeholder the record carries in its place is never returned, and the `CG` field is left in `tags`.
- `sequence` (str) SEQ, unpacked from its 4-bit encoding, or "*" when the record carries none.
- `qualities` (str) QUAL as phred+33, or "*" when the record carries no qualities at all. A record with only some of them missing spells those "*" in place.
- `next_chr` (str) RNEXT, "*" when the mate sits on no reference.
- `next_start` (int) PNEXT.
- `template_length` (int) TLEN.
- `bai_bin` (int) Index bin the record declares itself in.
- `reference_length` `query_length` (int) Bases of the reference the alignment covers, and of the read its cigar consumes.
- `tags` (dict) Optional fields by two-letter tag, in the order the record stores them. Typed as the file types them: character, integer, float, string, or list of integers or floats. Empty when `parse_tags` is off.
- `is_paired` `is_proper_pair` `is_mapped` `is_next_mapped` `is_reverse` `is_next_reverse` `is_first_in_pair` `is_last_in_pair` `is_secondary_or_supplementary` `is_failed_qc_or_duplicate` (bool) The `flag` bits, decoded. Finer ones are yours to mask off `flag`.

Notes:
- `cigar`, `sequence`, `qualities` and `tags` are decoded the first time they are read and kept afterwards, so an alignment read for its coordinates never pays for the rest of it. The optional fields are only walked when `tags` is read, so a record with malformed ones reads fine and raises there.
- `to_dict()` returns every field as a plain dict under the same names, for pandas, for serialising, or for sending to another process, an alignment itself not being picklable.

### Read HiC values

```python
values = reader.read_values(chr_ids, starts, ends, ...)
```

Parameters:
- `chr_ids` `starts` `ends` Chromosome IDs, starts and ends of the two locations.
- `bin_size` Input bin size or -1 to use the smallest. Must be available in the file. Smallest by default.
- `bin_count` Approximate output bin count. Takes precedence over `bin_size` if specified by selecting the closest bin size resulting in `bin_count`. Not specified by default.
- `exact_bin_count` Resize output to match `bin_count` (if specified). Not by default.
- `full_bin` Extend locations ends to overlapping bins if true. Not by default.
- `def_value` Default value to use when no data overlap a bin. 0 by default. A bin holding a contact the file cannot value — one the chosen `normalization` has no factor for, or an expected value of zero — comes back NaN instead, that being a different answer from not having been observed at all.
- `triangle` Skip symmetrical data if true. Not by default. On one chromosome a hic file stores one side of the diagonal only, and `triangle` reads that side alone rather than mirroring it, so a window lying entirely on the *other* side comes back empty — `starts=[15_000_000, 10_000_000]` gives `def_value` throughout where the mirrored request gives the data. Leave it off unless you know which side your window is on.
- `min_distance` `max_distance` Min and max distance in bp from diagonal for contacts to be reported. All by default.
- `normalization` Either "none" or any normalization available in the file, such as "kr", "vc" or "vc_sqrt". "none" by default.
- `mode` Either "observed", "oe" (observed/expected) or "expected". "observed" by default.
- `unit` Either "bp" or "frag". "bp" by default.
- `save_to` Save output to this .npz path (under "values" key) and return nothing. Not by default.

Returns a numpy float32 array of shape (loc 1 bins, loc 2 bins).

### Read HiC sparse values

```python
values = reader.read_sparse_values(chr_ids, starts, ends, ...)
```

Parameters:
- `chr_ids` `starts` `ends` `bin_size` `bin_count` `exact_bin_count` `full_bin` `triangle` `min_distance` `max_distance` `normalization` `mode` `unit` `save_to` Identical to `read_values` method. There is no `def_value`: a cell this does not list is one no contact reached, which is what a sparse matrix says by leaving it out.

Returns a COO sparse matrix as a dict with keys:
- `values` Values as a numpy float32 array.
- `row` Values rows indices as a numpy uint32 array.
- `col` Values columns indices as a numpy uint32 array.
- `shape` Shape of the dense array as a tuple.

Convert in python using `scipy.sparse.csr_array((x["values"], (x["row"], x["col"])), shape=x["shape"])`.

### Open bigWig and bigBed files for writing

```python
writer = gwseq_io.open(path, mode, type, chr_sizes, genome, fields, items_per_slot, compression_level)

with gwseq_io.open("path/to/file.bigwig", "w", genome="mm10") as writer:
    ...
```

Parameters:
- `mode` Opening mode. Must be set to "w" (write).
- `type` Type of file to write. Either "bigwig" or "bigbed". "bigwig" by default.
- `chr_sizes` Map of chromosome IDs and their sizes. Every written coordinate is checked against it, and a written ID is resolved against its keys the way a read one is, so "1" and "chr1" reach the same entry. Inferred from what is written if omitted, a chromosome then ending where its last value or entry does. None by default.
- `genome` Genome ID to get the chromosome IDs and their sizes, as `get_chr_sizes` returns them. May not be set with `chr_sizes`. None by default.
- `fields` Entries keys and their types. Only for bigBed files. The first three must be the coordinates `chr`, `start`, `end` (or the aliases `chr_id`, `chrom`, `chromStart` and `chromEnd`). Anything else is refused when the file is written. Types may be "string", "int", "uint" or "float". {"chr": "string", "start": "uint", "end": "uint", "name": "string"} by default.
- `items_per_slot` Values or entries one block holds, and records one zoom block holds. Use -1 for recommended (1024 for bigWig, 512 for bigBed, as the UCSC writers use). -1 by default.
- `compression_level` zlib level for the data and zoom blocks, or 0 to leave them uncompressed. 6 by default.
- `parallel` Number of threads compressing blocks. Use -1 for recommended (one per core, capped at 12). -1 by default.

Attributes:
- `path` `type` `closed` Path being written, "bigwig" or "bigbed", and whether `close` has run.
- `chr_sizes` Chromosome sizes as they will be written, in the order the chromosomes were written.
- `section_count` `section_counts` Sections, or blocks of entries, written so far — in total and, for a bigWig, by encoding ("bedgraph", "varstep", "fixedstep"). Each section takes whichever of the three costs the fewest bytes. Counted when a block is placed in the file, not when it is handed over, so with `parallel > 1` a section still being compressed is not in the tally yet; the counts are complete once `close()` has run.
- `entry_count` `fields` Entries written so far, and the columns they are written with (bigBed only).
- `skipped_count` Values dropped for not being finite.
- `clipped_count` Values cut back to the end of a declared chromosome they hung over.

Notes:
- `close`, which is what leaving the `with` block runs, finishes the file and stops the compression threads. Nothing on disk is a bigWig or bigBed until it returns, and nothing of the writer's is still running once it has. A writer dropped without it is closed by the garbage collector instead, which is later and not up to you.
- Only chromosomes that were actually written go into the file, so `chr_sizes` and `genome` are a bounds check and a spelling of the names rather than a list of what the file will contain.
- A bigWig value that starts inside a declared chromosome and ends past it is written up to that end rather than refused, since a chromosome is rarely a whole number of bins long. One that starts at or past the end raises, as does any bigBed entry running past it.

### Write bigWig values

```python
writer.write_value(chr_id, start, end, value)
writer.write_values(chr_id, start, span, values)

writer.write_value("chr1", start=1000, end=1010, value=0.1)
writer.write_values("chr1", start=1000, span=10, values=[0.1, 0.3, 0.2, 0.1])
```

Parameters (write_value):
- `chr_id` `start` `end` Chromosome ID, start and end of value.
- `value` Location value.

Parameters (write_values):
- `chr_id` `start` Chromosome ID and start of first value.
- `span` Window in bp of each successive locations relative to their starts, so `values[n]` covers `[start + n * span, start + (n + 1) * span)`.
- `values` Locations values, as a list or a numpy array.

Notes:
- Values must be pooled by chromosome, added in order and without overlap.
- For better performance, hand successive locations of one span over in one `write_values` call.
- Each section is written in the narrowest of the three encodings that holds it: fixedStep while every value shares one span and one step, four bytes a value; variableStep once the starts turn irregular, eight; bedGraph once the spans differ too, twelve.
- A NaN or infinite value is not written. It leaves a gap, which is what a bigWig means by a base carrying no data, and a reader fills it with the `def_value` it was asked for. `skipped_count` counts them.
- A chromosome that received only non-finite values still enters the file's chromosome list, with a size of 1 where its size was not declared — the writer sizes an undeclared chromosome from what reached it, and nothing did. Declare `chr_sizes` if a chromosome has to keep its real length whatever lands on it.
- A value hanging over the end of a declared chromosome is written up to that end, which is what makes a span that does not divide a chromosome ordinary rather than an error. `clipped_count` counts them. A value starting at or past the end raises, and in a `write_values` run only the last value can hang over, so a run reaching whole values past the end raises too.

### Write bigBed entries

```python
writer.write_entry(chr_id, start, end, ...)

writer.write_entry("chr1", start=1000, end=1010, fields={"name": "read#1"})
```

Parameters:
- `chr_id` `start` `end` Chromosome ID, start and end of entry.
- `fields` Map of additional fields as specified in file `fields`. The first three declared fields are the coordinates and are written from `start` and `end`, so naming one here is an error. A declared field left out is written empty for a string and 0 for a number.

Notes:
- Entries must be pooled by chromosome and added in order of their start. Unlike bigWig values, they may overlap and nest freely.
- `end` may equal `start`, which BED allows and is how an insertion is written. Such an entry covers no base, so it adds nothing to the coverage the summary and the zoom levels describe, and it is read back by the location containing the base it names. `end` before `start` raises.
- The summary statistics a bigBed carries, and its zoom levels, describe the depth of coverage its entries make, as the format asks: a base under three entries counts once towards the bases covered and three towards the sum.

### Convert bedGraph or WIG to bigWig

```python
gwseq_io.convert_to_bigwig(input_path, output_path, ...)

gwseq_io.convert_to_bigwig("track.bedgraph.gz", "track.bigwig", genome="mm10")
```

Parameters:
- `input_path` Path to input bedGraph or WIG file. May be gzipped.
- `output_path` Path to output file.
- `bin_size` Force a specified bin size in the output. A bin holds the base-weighted mean of what falls in it, so a 500 bp interval counts for five hundred times what a 1 bp one does, and a bin nothing covers is left as a gap. Takes bins as is by default.
- `chr_sizes` `genome` `items_per_slot` `compression_level` `parallel` Identical to `open` in write mode.
- `progress` Function called during conversion. Takes the bytes read and the total size of the input as parameters. Use default callback function if true. None by default.

Returns a map of `format` ("bedgraph" or "wig", as it was sniffed), `line_count`, `item_count`, `skipped_count`, `clipped_count` (values cut back to the end of their chromosome) and `chr_sizes` as written.

Notes:
- Which of the two formats the input is comes from its content, not from its name: the first line that is neither blank, a comment, nor a `track` or `browser` declaration decides. A `fixedStep` or `variableStep` line makes it a WIG, four columns of chromosome, start, end and value a bedGraph, and anything else is refused. The format is settled once, so a file holding both is refused as well.
- WIG coordinates are 1-based and bedGraph ones 0-based half-open. `step` and `span` both default to 1; a declaration with no `chrom`, or a `fixedStep` with no `start`, is an error rather than a guess.
- Values must be pooled by chromosome, in order and without overlap, as `write_value` asks — what `sort -k1,1 -k2,2n` gives. Input that is not raises, naming the line. Nothing is sorted or spooled, so a conversion of any size holds a megabyte of input and one open section.
- Nothing on disk is a bigWig until the call returns, exactly as for a writer.

### Convert BED to bigBed

```python
gwseq_io.convert_to_bigbed(input_path, output_path, ...)

gwseq_io.convert_to_bigbed("peaks.bed.gz", "peaks.bigbed", genome="mm10")
```

Parameters:
- `input_path` Path to input BED file. May be gzipped.
- `output_path` Path to output file.
- `chr_sizes` `genome` `items_per_slot` `compression_level` `progress` Identical to `convert_to_bigwig`.
- `fields` Entries keys and their types, as in `open` in write mode. Taken from the standard BED columns — `chrom`, `chromStart`, `chromEnd`, `name`, `score`, `strand`, `thickStart`, `thickEnd`, `itemRgb`, `blockCount`, `blockSizes`, `blockStarts`, then `field13` and up — for however many columns the first record has, by default.

Returns a map as `convert_to_bigwig` does, with `format` always "bed".

Notes:
- Lines are split on tabs, a BED being tab-delimited and its `name` column being allowed to hold spaces. Every record must carry the same number of columns as the first one, a bigBed storing one shape of record.
- A BED carries no column names of its own, so a file whose columns are named otherwise needs `fields` to keep them.
- Entries must be pooled by chromosome and in order of their start, but may overlap and nest freely, as `write_entry` allows.

### Convert SAM to BAM

Not implemented. The name exists and raises `Unsupported` (a
`NotImplementedError`).

```python
gwseq_io.convert_to_bam(input_path, output_path)
```

Parameters:
- `input_path` Path to input SAM file. May be gzipped.
- `output_path` Path to output file.

### Get genome chromosome sizes

```python
gwseq_io.get_chr_sizes(genome, ...)
```

Parameters:
- `genome` Genome name (eg, "mm10").
- `full` Include unplaced chromosomes if true. Not by default.

Returns a map of chromosome IDs and their sizes, sorted by chromosome ID as a string, so `chr10` comes before `chr2`. Genomes that are not bundled, and any call with `full`, are fetched from `api.genome.ucsc.edu` and cached for the process lifetime.

### Exceptions

Everything this library raises descends from `gwseq_io.Error`, and each leaf
also inherits the built-in that means the same thing, so ordinary Python
handling keeps working without knowing the hierarchy exists.

```
Error                       everything gwseq_io raises
├── InvalidRequest          the call asked for something impossible  (ValueError)
│   ├── UnknownChromosome   ... a name the file does not carry
│   └── ReaderClosed        ... through a reader that has been closed
├── InvalidFile             the file is not one this reads
│   └── CorruptFile         ... it is, and its bytes contradict each other
├── SourceError             the bytes did not arrive                 (OSError)
│   └── HttpError           ... from a URL, and the status says why
└── Unsupported             a real feature of the format, not implemented
                                                        (NotImplementedError)
```

```python
try:
    values = reader.read_values(chr_ids, starts, ends, bin_size=bin_size)
except gwseq_io.UnknownChromosome as e:
    ...     # the message lists what the file does carry
except gwseq_io.InvalidRequest:
    ...     # a bin size of zero, an end before its start, a bad reduction
except gwseq_io.SourceError:
    ...     # the file went away, or the network did
```

Notes:
- `UnknownChromosome` says what the file does hold, and — when the name is
  longer than the file's own name field can store — says that too, and names
  the chromosome the first characters spell. A file written from names too long
  for its field stores them truncated, and that is what the caller is looking
  at.
- `CorruptFile` carries the offset the contradiction was found at.
- A conversion complaining about its *input* raises `InvalidFile`, naming the
  line: a column that will not parse, a record out of order, a coordinate past
  the end of its chromosome. `InvalidRequest` is what the *call* asked for — a
  negative `bin_size`, a genome that does not exist.
- `InvalidRequest` is a `ValueError` and `SourceError` an `OSError`, so
  `except (ValueError, OSError)` still covers the common cases.


## Dev notes

### Project layout

```
gwseq_io/
├── Cargo.toml              # Rust workspace (three crates, shared versions)
├── pyproject.toml          # PEP 517 build config (maturin)
├── dist.py                 # build, check and publish the source distribution
├── ARCHITECTURE.md         # the design, and why each piece is shaped as it is
├── docs/                   # the file format specifications this implements
├── archives/               # previous releases, zipped
├── fuzz/                   # cargo-fuzz targets, for longer runs than the tests take
├── tools/                  # the generator for the bundled genome table
└── crates/
    ├── gwseq-io/           # the library — no Python in it
    │   └── src/
    │       ├── bbi/        # bigWig / bigBed reader, writer and text converters
    │       ├── hic/        # HiC reader
    │       ├── bam/        # BAM reader (header, BGZF, BAI index, records)
    │       ├── genomic/    # loci, bins and chromosome maps
    │       ├── genomes/    # built-in genome chromosome sizes
    │       └── source/     # byte sources: local, HTTP, cache, gzip, and sinks
    ├── gwseq-io-py/        # the PyO3 extension and the gwseq_io package
    │   ├── src/            # bindings, one module per format
    │   └── python/         # Python package entry-point
    └── gwseq-io-cli/       # `gwseq`, a front end that is not a binding
```

The split is the point: `gwseq-io` has no PyO3 in it, so the library is usable
from Rust and the binding layer is thin enough to read. `gwseq-io-cli` exists to
keep it honest — anything the CLI cannot reach is a feature that only exists as
a Python argument.

[ARCHITECTURE.md](ARCHITECTURE.md) is the design in detail — the source layer,
the concurrency model, the extraction kernels — and records the traps this code
was built against, most of which are only visible once you have hit them.

### Build from source

| Dependency | Version | Notes |
|---|---|---|
| Python | ≥ 3.9 | the interpreter alone — PyO3 declares the C API in Rust, so no `Python.h` is involved |
| Rust | ≥ 1.85 | via [rustup](https://rustup.rs); `cargo` comes with it |
| maturin | ≥ 1.7 | installed automatically as a build requirement |
| numpy | any | a runtime dependency, installed by pip |

maturin is resolved by pip from `pyproject.toml`, so it never needs to be
installed by hand. All that has to come from the OS is a Rust toolchain and a
linker.

#### Prerequisites

```bash
xcode-select --install
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh

# Debian / Ubuntu
sudo apt install python3-venv
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh

# Fedora / RHEL
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh

# Arch
sudo pacman -S --needed python rustup && rustup default stable
```

```powershell
# Windows
winget install Python.Python.3.12
winget install Rustlang.Rustup
```

Nothing in those lists is a compiler: rustup brings its own linker driver on
every platform, and there are no C sources to build. No `-dev` / `-devel`
Python package either — the extension is built against the stable ABI and
includes no Python header. On
Windows the import library the official installer ships is what the linker
wants, and it is there by default.

#### Build

```bash
# 1. Create and activate a virtual environment (recommended)
python -m venv .venv
source .venv/bin/activate      # Windows: .venv\Scripts\activate

# 2. Build and install the package in editable mode, with the dev extras
pip install -e ".[dev]"
```

Iterating on the Rust is faster through maturin directly, which rebuilds only
what changed:

```bash
pip install maturin
maturin develop --release
```

`maturin` installs into `CONDA_PREFIX` when one is set, whatever `VIRTUAL_ENV`
says — `unset CONDA_PREFIX` first if a conda environment is active, or the
build lands somewhere you are not importing from.

To build the source distribution instead of installing in place:

```bash
python dist.py                 # build it
python dist.py --check         # and compile it in a clean venv
```

`dist.py --upload` publishes to PyPI, `--test-upload` to TestPyPI; both compile
the tarball in a throwaway environment first, since building an sdist never runs
cargo and that check is the only thing standing between a compile error and the
index. Nothing but the source is published: no wheel is uploaded, so every
install compiles from it.

#### Working on the Rust alone

```bash
cargo build -p gwseq-io -p gwseq-io-cli   # the library and the CLI
cargo test --workspace
cargo test --release --workspace          # see ARCHITECTURE.md §12 on why both
cargo clippy --workspace --all-targets --all-features
cargo fmt --all --check
```

`cargo build --workspace` is not in that list on purpose: it tries to *link* the
extension, and a `cdylib` full of undefined Python symbols only links on macOS
with `-undefined dynamic_lookup`, which maturin passes and plain cargo does not.
Build the extension with maturin, and `cargo check -p gwseq-io-py` when you only
want the type errors.

If the checkout is on a synced drive (OneDrive, Dropbox, iCloud), keep the build
out of it — `target/` is thousands of files that will be indexed and uploaded
whatever the ignore file says:

```bash
export CARGO_TARGET_DIR="$HOME/.cache/gwseq_io/target"
```

### Tests

```bash
cargo test --workspace                 # 366 tests + 4 doctests, ~25 s
cargo test --release --workspace       # the same, in the profile that ships
cargo test --workspace -- --ignored    # the decompression bomb, ~16 s
python crates/gwseq-io-py/python/api_smoke.py   # the Python surface
```

**Everything above is self-contained** — no fixtures, no network, no second
checkout. That is deliberate: a test that only runs on one machine is a test
that stops running. The HTTP source is covered by a server the suite starts
itself on a loopback port the OS picks, which is the one way to exercise the
case that matters there: a server that ignores `Range`.

Three layers, in the order they catch things:

- **Unit tests, in-crate.** Every format module has `#[cfg(test)]` tests over
  byte literals — a hand-built R-tree node, a truncated header, a BGZF block.
  Where a corrupt-input bug should surface.
- **Round trips**, over files the tests build. `roundtrip.rs` writes with the
  writer and reads with the reader: a file sniffs as what was written, a walk
  concatenates to a whole-chromosome read, a closed reader keeps its headers
  and refuses to read, and answers do not change with the thread count.
  `bam_roundtrip.rs` and `hic_roundtrip.rs` build a real BAM (with its BAI) and
  a real hic byte by byte, since neither format has a writer to round-trip
  through, and read them back through the public API.
- **Properties** (`crates/gwseq-io/tests/properties.rs`), over inputs proptest
  chooses. The unit tests check the cases someone thought of; these check the
  ones nobody did, and shrink a failure to the smallest input that still shows
  it. `PROPTEST_CASES=5000` for a longer run; the shrunk seeds that once failed
  are committed in `tests/properties.proptest-regressions` and re-run first.

Every parser is fuzzed on each `cargo test` — `crate::fuzz` runs all twenty of
them over random bytes, corrupted headers and every truncation of a valid file,
through both a bare source and the block cache every reader sits behind, in a
second or two. `fuzz/` holds cargo-fuzz targets for a longer look:

```bash
cargo test -p gwseq-io --lib fuzz
cargo +nightly fuzz run open_any       # needs cargo-fuzz
```

#### Benchmarks

The read benchmarks want a real file and look for one in `local/test_data/`,
which is not committed; each skips itself when its fixture is absent. The write
benchmarks generate their input and run anywhere.

```bash
cargo bench -p gwseq-io --bench extract
cargo bench -p gwseq-io --bench writer
```

