Binning position scores with `binning_tool`

Many analyses that compare genome-wide tracks against each other — correlating conservation with chromatin accessibility, clustering a set of ATAC tracks, feeding tracks to a model — begin the same way: cut the genome into fixed-size bins and reduce every track to one number per bin. binning_tool does that first step in one command. It takes a run definition in YAML naming the bins and the position-score resources to bin, and writes one HDF5 file holding a bins × tracks matrix with the bin coordinates and the provenance of every track stored beside it.

binning_tool run.yaml

Only position_score resources can be binned in this version. Each track is one score of one resource, reduced over each bin by an aggregator such as mean or max.

The run definition

A run definition has three parts: the reference genome, the bins block, and a list of binners entries.

input_reference_genome: hg38/genomes/GRCh38-hg38

bins:
  bin_size: 10240
  regions: [chr1, chr2, "chrX:1-50000000"]

binners:
- position_score_binner:
    resource_query: "hg38/scores/*"
    search_term: conservation
- position_score_binner:
    resource_query: hg38/scores/phyloP100way
    aggregator: max
    none_value_replacement: 0.0

Every key is checked. A mistyped key anywhere in the file — aggregtor for aggregator — is an error that names the entry, never a silently applied default.

Reference genome

input_reference_genome (required) names the reference genome resource whose chromosomes define the grid. It is the run definition’s business alone: there is no command-line flag for it, and a genome offered by the genomic context is not consulted. The chromosome lengths decide how many bins there are and what each one spans, so one run definition always describes one matrix — whoever runs it, and wherever.

The GRR itself is never named in the run definition; it comes from -g, --grr-directory, or the default GRR definition, as for every GAIn tool (see Genomic resources and repositories). One run definition therefore runs unchanged on a laptop against a cached GRR and on a cluster node against a node-local one.

Bins

bin_size (required) is the width of every bin in base pairs.

regions (optional) lists what to bin, in GAIn’s region notation: chr21 is the whole chromosome and chr21:20000000-22000000 is a window with inclusive bounds. Omit regions to bin every chromosome of the reference genome in genome order. The rows of the output follow the listed order, so listing the chromosomes you want is also how to leave out alternate and unplaced contigs. A window that ends before it starts, a region beyond its chromosome’s length, a chromosome the genome does not have, and two regions that overlap are all errors reported before any work runs, naming the offending entries.

Binner entries

binners (required) is a list of entries. Each entry is a one-key mapping whose key is the binner kind; the only kind in this version is position_score_binner, and its value takes these keys:

resource_query (required): A repository search that selects the resources to bin, resolved through the same search that grr_browse uses: an exact resource id, or a glob over ids such as hg38/scores/*, optionally followed by a bracketed label filter. Only position_score resources are taken from the matches, so a glob that also matches allele scores or other resource types is safe. Matches are ordered by resource id. A query that matches no position score is an error naming the entry, so a typo cannot silently produce a file with fewer tracks.
search_term (optional): A full-text filter conjoined with resource_query, in the syntax of grr_browse -s (see Genomic resources and repositories). It needs a repository that carries a full-text index; the public IossifovLab GRR does, and a directory GRR has one once grr_manage repo-index has been run on it. A search_term on a repository without an index is an error naming the entry.
aggregator (optional): How the positions inside a bin are reduced to one number. The default is the aggregator the score itself declares as its position_aggregator in its resource configuration, so a conservation score and a coverage track are each reduced the way their authors intended. Setting aggregator on an entry overrides that for every resource the entry matches. Only aggregators that produce a number are accepted, and only numeric (int or float) scores can be binned; a string-typed score, or an aggregator such as join or list that builds a string or a list, is refused at parse time. The numeric aggregators are max, min, mean, median and count. To bin one resource under two aggregators, list it in two entries.
none_value_replacement (optional): A value fed to the aggregator for every position no record covers. By default there is no replacement: a bin that no record covers is NaN in the output, so missing data cannot masquerade as a measurement in a correlation. Set a replacement, typically 0.0, when “no signal” is genuinely zero for the track, as it is for a coverage track.

A resource that defines more than one score is refused, with its scores listed: a track is exactly one score, and this version offers no key to pick among several.

Track names

Every entry expands into one track per matched resource, and a track is named by its resource id. Tracks follow entry order, and within an entry the order of resource ids. When the same resource occurs more than once in the expanded list — typically the same resource matched by a glob in one entry and named again with another aggregator in a second — every track of that resource gets :<aggregator> appended, whichever entry came first, so hg38/scores/phyloP100way:mean and hg38/scores/phyloP100way:max sit side by side. Two entries that would still produce one and the same track (same resource, same aggregator) are refused, naming both entries: nothing in the file could tell the two columns apart. Names in the output are therefore always unique. There is no key to rename a track; renaming is a one-line change on the tracks table after the fact.

Coordinates and the grid

Bin coordinates are 1-based and inclusive, the convention of every GAIn reader and of the region notation, so a bin’s chrom:start-end can be pasted into any other GAIn tool unchanged.

Bins follow a global grid anchored at position 1 of every chromosome: with a bin size of 10240 the bins are 1-10240, 10241-20480 and so on, regardless of where a region starts. A window that does not start on a grid boundary therefore has its first and last bins clipped to the window — the bin bounds always name exactly what was aggregated — and bins from different runs, regions, or bin sizes that divide each other line up and can be joined by coordinate.

Running the tool

binning_tool RUN_DEFINITION [-o OUTPUT] [-w WORK_DIR] [--keep-work-dir] [--dry-run]
             [--task-budget BP] [-g GRR] [-j N] [--force] ...

The only positional argument is the run definition. -o names the HDF5 file to write; by default it is the run definition’s path with an .h5 suffix, beside it, so binning_tool run.yaml writes run.h5.

What the run describes — the genome, the grid and the tracks — is the run definition’s business; the command line only says where to find the GRR, where to put the output, and how to run the work.

Dry run

--dry-run reads the run definition, resolves every query against the GRR, checks every rule described above, prints the list of tracks the run would produce — name, resource id, score id and aggregator — together with the number of regions, bins and tasks, and exits without writing anything. Use it to see what a query matched, how large the matrix will be and how the work will be cut before committing cluster time.

Parallelism

The work is split into one task per track and bundle of regions, followed by one task that assembles the HDF5 file. Consecutive regions of the run definition are packed, in order, into bundles of at most --task-budget bases (default 50,000,000); a region is never split, so a chromosome longer than the budget is a task of its own, while the hundreds of alternate and unplaced contigs of a human genome — under two percent of its bases — pack into a handful of tasks instead of one each. The budget only decides how many tasks there are: the file is the same whatever its value.

Both ends of the dial are reachable. A budget of 1 is one task per region — a region is never split, so the smallest budget that cuts at all cuts everywhere. A budget of 0 or less is no budget at all: the whole run goes into one task per track. A task opens its track’s resource once, whatever the budget, so the whole-run task opens it once for the run.

Which end helps depends on how many tracks a run has, because a track’s work is never split across tasks at budget 0: the achievable parallelism is then exactly the number of tracks. A run of hundreds of tracks on a few workers loses nothing and sheds tasks; a run of two or three tracks with -j 8 leaves most of the workers idle and should keep the default. Note too that a task is the unit a rerun repeats: at budget 0 a task that fails part-way recomputes its whole track next time, where at the default a rerun keeps every bundle that finished.

The tasks run through the same task graph as the annotation tools, so the same flags apply: -j N sets the number of workers, -N names a configured dask cluster, and --task-log-dir keeps a log per task. The genomic-context flags (-g, --grr-directory, -R) and the verbosity flags (--verbose, --logfile) are likewise the shared ones. Run binning_tool --help for the full list.

Work directory and reruns

Each task writes one column chunk per region of its bundle into a work directory; the final task assembles the file from those chunks, region by region, so at no point is the whole matrix in memory — a genome-wide run at a small bin size is possible on an ordinary node. The work directory is -w; by default it is a sibling of the output named after it (run_work next to run.h5), and the task-status directory lives inside it. A work directory the tool created is removed after a successful run; --keep-work-dir keeps it.

A run that was interrupted resumes from the tasks it had finished when it is started again with the same work directory and budget: a task whose chunks are all present is skipped, a task missing any of its chunks is computed again in full, and the file is assembled if it is missing. A run whose chunks and output are all present does nothing. The chunks are keyed by everything that decides their values — resource, score, aggregator, replacement, bin size and region — so two run definitions sharing a work directory share exactly the chunks they compute identically. A rerun does not notice that a resource has changed underneath it; pass --force to recompute every chunk, for example after a resource was updated in the GRR.

The output file

The output is one HDF5 file in a plain layout that any HDF5 reader understands without a library beyond h5py. It holds three datasets:

/values: float64, shape (n_bins, n_tracks). NaN where a bin has no data. Stored in row blocks and gzip-compressed, so reading every track for one chromosome is one contiguous read and the NaN- and zero-heavy tracks compress well.
/bins: A compound (structured) dataset of shape (n_bins,) with fields chrom (fixed-length bytes, sized to the longest chromosome name in the run), start and end (int64, 1-based inclusive).
/tracks: A compound dataset of shape (n_tracks,) with fields name, resource_id, score_id and aggregator (variable-length UTF-8 strings) and none_value_replacement (float64, NaN when the entry set none).

Row i of /bins describes row i of /values, and row j of /tracks describes column j. The root attributes record what the run was, so that two files can be checked for comparability before they are joined:

Attribute	Content
`input_reference_genome`	The resource id of the reference genome the grid was built on.
`bin_size`	The bin width in base pairs.
`regions`	The binned regions, in row order, each as `chrom:start-end`.
`coordinates`	`"1-based-inclusive"`.
`gain_version`	The GAIn version that wrote the file.
`created`	The UTC time the file was written, in ISO 8601.

Every value in /values is a float64, including the result of an integer-valued aggregator such as count or the max of an integer score: the matrix has one dtype, and HDF5 has no null for integers. A bin no record covers is NaN unless the entry set none_value_replacement, in which case the replacement was fed to the aggregator for every uncovered position and the bin holds the aggregate of that.

Two notes on scores stored as bigWig. First, a bigWig stores float32, so binning the bigWig form and the text (tabix) form of the same score does not give bit-identical matrices; the two agree to roughly 1e-8 through a mean over 10 kb bins. Second, for a binning run the storage format of a score dominates its cost: binning the same per-base conservation score over chromosome 21 took about three times longer from its tabix form than from its bigWig form, while the scope of the tabix index (whole genome or one chromosome) made no difference. When a score is published in both forms, prefer the bigWig for binning.

Reading the file back

With h5py alone:

import h5py

with h5py.File("binning_run.h5", "r") as h5:
    values = h5["values"][()]       # numpy float64 array, bins x tracks
    bins = h5["bins"][()]           # structured array: chrom, start, end
    tracks = h5["tracks"][()]       # structured array: name, resource_id, ...
    attrs = dict(h5.attrs)

print(values.shape, attrs["bin_size"], attrs["coordinates"])
print(tracks["name"][:3])
print(bins[:3])

With pandas, decoding the fixed-length chromosome names once:

import h5py
import pandas as pd

with h5py.File("binning_run.h5", "r") as h5:
    tracks = pd.DataFrame(h5["tracks"][()])
    for column in ("name", "resource_id", "score_id", "aggregator"):
        tracks[column] = tracks[column].str.decode("utf-8")
    bins = pd.DataFrame(h5["bins"][()])
    bins["chrom"] = bins["chrom"].str.decode("utf-8")
    values = pd.DataFrame(
        h5["values"][()], columns=tracks["name"],
        index=pd.MultiIndex.from_frame(bins))

print(values[["hg38/scores/phyloP100way:mean",
              "hg38/scores/phyloP100way:max"]].head())

values is then a DataFrame indexed by (chrom, start, end) with one column per track, and tracks explains each column.

A worked example

The run definition below bins the conservation scores of the public IossifovLab GRR over two megabases of chromosome 21 (binning_run.yaml):

input_reference_genome: hg38/genomes/GRCh38-hg38

bins:
  bin_size: 10240
  regions: ["chr21:20000000-22000000"]

binners:
- position_score_binner:
    resource_query: "hg38/scores/*"
    search_term: conservation
- position_score_binner:
    resource_query: hg38/scores/phyloP100way
    aggregator: max

The first entry is a glob over every resource under hg38/scores/ narrowed by the full-text term conservation; of the fourteen resources there, eight are position scores — the phastCons and phyloP scores computed over 7, 20, 30 and 100 species — and those are what the entry expands to. The second entry names hg38/scores/phyloP100way again under max, so that resource is binned twice and both of its tracks carry their aggregator in their name.

Check what the run would produce first:

binning_tool binning_run.yaml --dry-run

tracks:
  hg38/scores/phastCons100way       hg38/scores/phastCons100way     phastCons100way mean
  hg38/scores/phastCons20way        hg38/scores/phastCons20way      phastCons20way  mean
  hg38/scores/phastCons30way        hg38/scores/phastCons30way      phastCons30way  mean
  hg38/scores/phastCons7way hg38/scores/phastCons7way       phastCons7way   mean
  hg38/scores/phyloP100way:mean     hg38/scores/phyloP100way        phyloP100way    mean
  hg38/scores/phyloP20way   hg38/scores/phyloP20way phyloP20way     mean
  hg38/scores/phyloP30way   hg38/scores/phyloP30way phyloP30way     mean
  hg38/scores/phyloP7way    hg38/scores/phyloP7way  phyloP7way      mean
  hg38/scores/phyloP100way:max      hg38/scores/phyloP100way        phyloP100way    max
regions: 1
bins: 196

Then run it:

binning_tool binning_run.yaml -j 4

This writes binning_run.h5 beside the run definition: a 196 × 9 matrix. The window chr21:20000000-22000000 does not start on a grid boundary, so the first bin is clipped to chr21:20000000-20008960 and the last to chr21:21995521-22000000, while every bin in between is a full 10240 bp.

The eight scores are per-base tracks stored as bigWig files of 6 to 10 GB each, read over HTTP by range, so this two-megabase run takes well under a minute. A whole-genome run over the same eight tracks would fold about three billion positions per track — on the order of forty minutes of CPU per track — which is why the example restricts regions: try a window first, then widen it.

Note

Reading a bigWig directly from a remote GRR needs a pyBigWig built with remote (libcurl) support. The conda package of GAIn ships one; the pyBigWig wheel on PyPI does not. If a run fails with this pyBigWig build has no remote-file support, install GAIn from conda, or point the tool at a local copy of the repository with --grr-directory. A GRR definition with a cache_dir (see Getting started on CLI) also works, but it downloads each bigWig in full before opening it — about 60 GB for the eight scores of this example — whatever the size of the window.

Reading the result back with the pandas snippet above gives:

name                     hg38/scores/phyloP100way:mean  hg38/scores/phyloP100way:max
chrom start    end
chr21 20000000 20008960                      -0.014157                         3.356
      20008961 20019200                      -0.006858                         3.539
      20019201 20029440                       0.001797                         3.694
      20029441 20039680                      -0.017875                         2.491
      20039681 20049920                       0.006743                         2.455

Binning position scores with binning_tool