gwseq_io

Python library for processing bigWig, bigBed, BAM, CRAM 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, CRAM and HiC files for reading

reader = gwseq_io.open(path, ...)

with gwseq_io.open("path/to/file.bigwig") as reader: # .bigbed .bam .cram .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 and CRAM files, where it defaults to the path of the file with ".bai" or ".crai" 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. A local CRAM with no .crai beside it is indexed by walking its container headers.
• reference Path or URL of the reference FASTA. may be plain or bgzip-compressed (less efficient). Needs a ".fai" (and ".gzi" if compressed) beside it, written by samtools faidx mm10.fa or bgzip -k mm10.fa && samtools faidx mm10.fa.gz. Only for CRAM files. A CRAM without reference opens and reads everything but the alignments sequences.

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 and CRAM files:

• type Either "bam" or "cram".
• 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 CRAM files only:

• version The format version as a "major.minor" string, eg "3.1".
• reference Path of the reference FASTA in use, or None when none resolved.
• reference_error Why there is no reference, or empty if resolved.

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

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. The default callback writes the time of day, what it has read of what it has to read, the percentage and an ETA to stderr — redrawn in place five times a second when stderr is a terminal, and written a line at a time at 1, 5, 10, 25, 50, 75 and 100% when it is redirected to a file or a pipe.

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

Quantify bigWig and bigBed values

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:

• reduce can't be set to "l1norm" if zoom is not -1.

Returns a numpy float32 array of shape (locations).

Profile bigWig and bigBed values

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

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.
• window Window size 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

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.
• as_cols Return one column per field instead of one entry per row. False by default.

Returns a list (locations) of list of BedEntry or dict of columns if as_cols is True.

The columns are typed from the file's own autoSql: start, end and every column the file declares int, uint or float come back as a numpy array, the rest as a list of str. A column whose values do not fit what it was declared to be stays text — UCSC declares itemRgb a uint and then writes 255,0,0 into it — so the keys and their lengths never depend on the data, only the dtype does, and only for a file that contradicts itself.

A BedEntry carries chr, start and end as attributes, and every column the file has by subscript: entry["name"]. keys(), len(), in and iteration work as they do on a dict, so dict(entry) builds one, and entry.to_dict() does the same in one step.

Read all bigBed entries

entries = reader.read_all_entries(...)

Parameters:

• chr_ids Only extract data from these chromosomes. All by default.
• col_count as_cols Identical to read_entries method.

Returns a list (locations) of list of BedEntry or dict of columns if as_cols is True.

Iterate over all bigBed entries

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

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

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 and CRAM entries

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.
• tags Decode the optional fields of every alignment. True by default.
• sequence Decode the alignments sequences. Only for CRAM files. True by default.
• as_cols Return one column per field instead of one SamEntry per row. False by default.
• tag_cols Optional fields to return as columns of their own, by two-letter name, e.g. ["NM", "MD"]. Needs as_cols and tags. Empty by default.

Returns a list (locations) of list of SamEntry, in the order the locations were given, or of dict of columns if as_cols is True. Each location gets its own full list, so two overlapping locations both report the alignments they share.

With as_cols, the keys are to_dict's — chr, start, end, read_name, flag, mapping_quality, cigar, sequence, qualities, next_chr, next_start, template_length, bai_bin, reference_length, query_length and the ten is_* flags — with the numeric ones as numpy arrays, plus one column per name in tag_cols. There is no tags column: a dict per alignment is the cost this avoids, and which tags an alignment carries is data, so a column per tag found would make the shape depend on it. An alignment without a named tag gets None. pd.DataFrame(columns) takes the result directly, and does so about six times faster than a list of to_dict()s.

Read all BAM and CRAM entries

entries = reader.read_all_entries(...)

Parameters:

• chr_ids Only extract data from these references. All by default.
• filter tags sequence as_cols tag_cols progress Identical to read_entries method.

Returns a list of SamEntry (as in read_entries), or a dict of columns if as_cols is True. Unplaced alignments are left out, the index reaching an alignment only through the reference it sits on.

Iterate over BAM and CRAM entries

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 tags sequence as_cols tag_cols progress Identical to read_entries method.
• sort_locations Read the locations in reference and position order, and report them in that order. False by default.

Returns an iterator over locations, each one a list of SamEntry. 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 and CRAM entries

iterator = reader.iter_all_entries(...)

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

Parameters:

• chr_ids filter tags sequence as_cols tag_cols progress Identical to read_all_entries method.
• window Window size in bp for each step. 1,000,000 by default.

Returns an iterator over successive windows, each one a list of SamEntry, 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 and CRAM entries

A SamEntry is one alignment, from a BAM or a CRAM file. CRAM records are rebuilt into exactly what a BAM would have held.

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.
• sequence (str) SEQ, unpacked from its 4-bit encoding, or "*" when the record carries none. A CRAM stores a mapped read as its differences from the reference, so with no reference or with sequence=False every base that came from the reference reads as "N" and the ones the file carries itself do not: a soft clip, an inserted base, a substitution. "MD" and "NM" go with them, being rebuilt from the reference too unless the file stored them.
• qualities (str) QUAL as phred+33, or "*" when the record carries no qualities at all.
• next_chr (str) RNEXT, or "*" 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 tags reading parameter 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.

Methods:

• tag(name) One optional field by its two-letter name, or None. Reads the record's bytes for that one field rather than building the whole tags dict.
• to_dict() 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.

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.
• Where a CRAM reference to decode the aligments sequences is looked for, in order: 1. the reference argument; 2. the slice's own embedded reference, when the file carries one; 3. the "UR" field of the header's "@SQ lines; 4. "REF_CACHE" and "REF_PATH", by the "M5" checksum the header states.

Read HiC values

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

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

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 are autoSql's: "string", "int", "uint" and "float" are the ones to reach for, and "lstring", "char", "short", "ushort", "byte", "ubyte", "bigint", "double", an array of any of them ("int[blockCount]"), an enum(...) or a set(...) are accepted so that a file read here can be written back with the columns it came with. {"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

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

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

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

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

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.

Pileup BAM or CRAM to bigWig

gwseq_io.pileup_reads(input_path, output_path, ...)

gwseq_io.pileup_reads("in.bam", "out.cpm.bigwig", normalize="cpm") # ChIP
gwseq_io.pileup_reads("in.bam", "out.cpm.bigwig", shift_reads=-36, transform_reads="force_extend", extend_size=80, normalize="cpm") # ATAC
gwseq_io.pileup_reads("in.bam", "out.cpm.bigwig", bin_size=1, transform_reads="splice", normalize="cpm") # RNA-seq

Parameters:

• input_path Path to input BAM or CRAM file. Must be sorted and indexed.
• output_path Path to output bigWig file.
• bin_size Output bin size in bp. 10 by default.
• filter_reads Filters out PCR duplicates, and reads whose MAPQ is below the given value (0 to drop duplicates only). No filtering by default.
• max_fragment Maximum fragment length in bp for paired reads. Longer pairs are discarded. 2000 by default.
• shift_reads Shift reads (in bp) relative to their 5' extremity. Cannot be used with transform_reads set to "extend". 0 by default.
• transform_reads Extend reads from their 5' (shifted) extremity to mate (if paired) or extend_size bp (if unpaired) if "extend", force extension to extend_size bp (even if paired) if "force_extend", do not extend reads and splice them (ie, only exons) if "splice". "extend" by default.
• extend_size Extend reads (in bp) from their 5' (shifted) extremity or to read length if set to 0, only applied if transform_reads is "extend" or "force_extend". 150 by default.
• def_value Default value to use for bins without coverage, or None to write nothing for them at all. 0 by default.
• normalize Normalization mode, either "count" (raw interval counts, no scaling), "cpm" (counts per million counted intervals), "bpm" (bin as its share of the total signal, per million) or "enrichment" (fold over the geometric mean of the bins that have coverage). "count" by default (no normalization).
• scale A further multiplier, applied after the normalization. 1.0 by default.
• window_size Bases a window covers, which is how much of the input is read and binned at a time. Use -1 for recommended, which derives it from the bin size. -1 by default.
• hold_bytes Bytes of binned runs a normalizing pileup may keep in memory rather than read the input a second time. -1 for the default (1 GiB, about 160 million reads); 0 to always read twice. -1 by default.
• items_per_slot Identical to open in write mode.
• compression_level zlib level for the data blocks, 0 to 9. 3 by default, where the rest of the library writes 6: on a pileup's own sections level 3 is both smaller and faster than 6.
• progress True for the default progress report described under read_values, or a callable taking (done, total). None by default.
• parallel Number of threads. Use -1 for recommended (one per core, capped at 12). -1 by default.

Returns a map of what was read and written: bin_size, window_size, normalize, the scale every value was multiplied by, how many passes over the input it took (1 unless the runs did not fit in hold_bytes), read_count, interval_count, bin_sum, covered_bins, item_count, section_count, section_counts, skipped_count, clipped_count and chr_sizes. read_count is what survived the filter and interval_count what was piled up, so the gap between them is what transform_reads dropped — improper pairs, unmapped mates, fragments longer than max_fragment.

Command line

python -m gwseq_io <command> <arguments...>

python -m gwseq_io pileup_reads in.bam out.cpm.bigwig -n cpm # ChIP
python -m gwseq_io pileup_reads in.bam out.cpm.bigwig -s -36 -t force_extend -e 80 -n cpm # ATAC
python -m gwseq_io pileup_reads in.bam out.cpm.bigwig -b 1 -t splice -n cpm # RNA-seq

One command so far, pileup_reads, which is the function above with its arguments spelled as flags: -b --bin-size, -f --filter-reads, -s --shift-reads, -t --transform-reads, -e --extend-size, -d --def-value, -n --normalize, -p --parallel, and --max-fragment, --scale, --window-size, --hold-bytes, --items-per-slot and --compression-level in full only. --no-def-value is def_value=None, which leaves an uncovered bin out of the file altogether. Negative values need no escaping: -s -36 is a shift of thirty-six bases against the read.

Get genome chromosome sizes

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)

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)
│       ├── cram/        # CRAM reader (containers, codecs, slices, records, references)
    │       ├── pileup/     # BAM and CRAM into a bigWig: read transforms, bins, runs
    │       ├── 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 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

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

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

# 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

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

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:

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

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:

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

Tests

cargo test --workspace                 # 657 tests + 5 doctests, ~15 s
cargo test --release --workspace       # the same, in the profile that ships
cargo test --release --workspace -- --ignored   # the decompression bomb, ~2 s
python crates/gwseq-io-py/python/api_smoke.py   # the Python surface

--ignored is one test — a gibibyte of zeros deflated, to prove the real inflate limit is the one wired in — and it wants --release: the same check in a debug build is minutes of arithmetic rather than the two seconds it takes here.

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. pileup.rs piles that BAM up and checks every bin against a brute-force pileup written in the test, which shares no code with the library.
• Properties (crates/gwseq-io/tests/suite/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:

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.

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

The pileup against the tool it ports

local/pileup_reads/ is the Python implementation pileup_reads was written from — a front end over samtools view and wigToBigWig. It is the specification, and local/probe/pileupcheck.py runs both over the same BAM and compares every bin of every chromosome:

python local/probe/pileupcheck.py                 # the README's three recipes,
                                                  # plus count and a filter
python local/probe/pileupcheck.py --chr chr19     # one chromosome, quicker

It needs samtools and a BAM in local/test_data/, and nothing else: the Python is asked for a bedGraph, so neither UCSC binary is involved.