Binning infrastructure

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

Two kinds of resource can be binned. A position_score resource gives one track per resource: one score, reduced over each bin by an aggregator such as mean or max. A fragment_score resource — the fragments of a single-cell ATAC sample, or any other collection of intervals — gives one track per group of fragments: every fragment, or the fragments of each cell type, or the fragments of each cell, counted or summed over the bin, or reduced as a depth of coverage.

The run definition

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

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
- fragment_score_binner:
    resource_query: "sc/atac_fragments/*"

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

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, one of two: position_score_binner bins position_score resources, and fragment_score_binner bins fragment_score resources. Both kinds select their resources the same way:

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 resources of the entry’s kind are taken from the matches — position_score resources for a position_score_binner, fragment_score resources for a fragment_score_binner — so a glob that also matches other resource types is safe. Matches are ordered by resource id. A query that matches no resource of the kind 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.

The other keys depend on the kind.

Position-score entries

A position_score_binner entry makes one track per matched resource and takes, besides resource_query and search_term, these keys:

aggregator (optional)

How the positions inside a bin are reduced to one number. The default is the aggregator the score itself declares 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, sum, product 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. The exception is aggregator: count: an uncovered bin is 0.0, because a count of nothing is zero. Set a replacement, typically 0.0, when “no signal” is genuinely zero for the track, as it is for a coverage track.

name (optional)

The base of the entry’s track names, in place of the resource id: a track is named <name>/<resource id>. Use it to bin one resource in two entries that would otherwise produce the same track name; see Track names.

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.

Fragment-score entries

A fragment_score_binner entry reduces, per bin, the values of the fragments of the matched resources. The entry’s mode decides which fragments reach a bin and with what weight (see Modes). The default mode, fragment_start, counts the fragments that start in each bin.

For a resource that follows the single-cell convention described below, the query is all an entry needs:

- fragment_score_binner:
    resource_query: "sc/atac_fragments/*"

Here sc/atac_fragments/donor1 and sc/atac_fragments/donor2 are two samples of a single-cell ATAC study, each labelled with its cell metadata table, sc/cell_meta_data, and with its sample. The entry pools both samples and counts the fragments that start in each bin per cell type, one track per cell type. To sum the read pairs behind the fragments instead, add value: {score_id: count}. The same entry with every default spelled out, and with the two samples kept apart, is:

- fragment_score_binner:
    resource_query: "sc/atac_fragments/*"
    pool: false
    group:
      cell_score_id: cell
      cell_meta_column: barcode
      group_meta_column: class
    value:
      value: 1
    aggregate:
      mode: fragment_start
      aggregator: sum
    meta:
      resource_label: cell_meta_resource_id
      filter:
      - column: sample_id
        label: sample_id

Besides resource_query and search_term, an entry takes these keys:

pool (optional, default true)

With true, every matched resource is merged by position into one set of tracks: a pooled track counts the fragments of all the samples together. With false, each resource is binned on its own and gets its own tracks.

group (optional)

How fragments are split into tracks, in exactly one of three forms; giving keys of two forms at once is an error.

  • {group: NAME} — one track, holding every fragment, whose group is the constant NAME.

  • {cell_score_id, cell_meta_column, group_meta_column} — one track per group of a metadata table, which meta names. cell_score_id is the fragment score carrying a fragment’s cell barcode, cell_meta_column the table column holding the barcode and group_meta_column the column holding the cell’s group. Each key is optional and defaults to the convention’s — cell, barcode and class — so group: {} groups by the convention throughout. The tracks are the distinct non-empty groups of the rows each job selects, in sorted order: the rows of every resource of a pooled entry together, or, with pool: false, of each resource on its own, so two samples may get different sets of tracks.

  • {group_score_id: S} — one track per distinct value of the str score S; see Grouping by a raw value. meta does not go with this form.

Without group, the resource’s labels decide (see Defaults).

value (optional, default {value: 1})

What one fragment adds, in exactly one of two forms. {score_id: S} is the value of the numeric (int or float) fragment score S, which every matched resource must have; {value: V} is the constant number V for every fragment.

aggregate (optional, default {mode: fragment_start, aggregator: sum})

How a bin reduces the values of its fragments, in three keys:

  • mode — which fragments reach a bin, and with what weight: fragment_start (the default), fragment_length or coverage_profile. See Modes.

  • aggregator — the reduction, sum when omitted. The aggregators that each mode accepts are in Aggregators in each mode. One aggregator applies to every resource of the entry.

  • uncovered_value — what each base that no fragment covers adds to its bin, a number or null. Only fragment_length and coverage_profile take it, and each has its own default. See Uncovered bases.

With the default value, sum in fragment_start mode counts fragments. value: {score_id: count} sums the read pairs of a single-cell resource.

name (optional)

The base of the entry’s track names, in place of the resource query of a pooled entry or the resource id of an unpooled one. Use it to bin one query in two entries, for example under two modes. See Track names.

meta (optional)

Where the table mapping barcodes to groups comes from: exactly one of

  • resource_id: ID — the data_frame resource ID;

  • resource_label: LABEL — the data_frame resource that the fragment resource’s label LABEL names;

  • file_name: PATH — a local file; see Local metadata files. file_format (csv, tsv or excel) and file_separator go with it and only with it. Without file_format the format follows the file’s suffix: .tsv is read as tsv, .xls and .xlsx as excel, and anything else as csv. Only .xlsx is readable out of the box: a legacy .xls file also needs the xlrd package, which GAIn does not install, and without it the run is refused when it resolves.

and filter, a list of conjuncts selecting the rows of one sample. Each is {column: C, value: V}, the column C equal to the literal V, or {column: C, label: L}, the column C equal to the fragment resource’s label L — so one filter selects each sample’s own rows. A row must satisfy every conjunct; no filter selects every row. A label whose value is a list cannot be a filter operand. meta goes only with the metadata form of group (or with no group on a resource the convention groups). Without it, the convention’s table is read.

There is no none_value_replacement key on this kind. In the fragment_length and coverage_profile modes, uncovered_value says what a base without fragments adds. In fragment_start mode, the aggregator decides what an empty bin holds (see Empty bins).

Modes

aggregate.mode decides which fragments reach a bin, and with what weight. Mode, value and group are independent: each mode operates with each value form and each grouping, pooled or unpooled.

fragment_start (the default)

A fragment adds its value, with a weight of 1, to the bin that holds its start position. It adds nothing to other bins, however far it extends. Adjacent bins and adjacent regions see each fragment exactly once. Thus the column of a fragment count sums to the number of fragments, and pooling several resources is exact.

fragment_length

A fragment adds its value to each bin that it overlaps. The weight is the number of base pairs of the overlap. A fragment that crosses a region boundary is clipped to each region, so each base pair adds to one region only.

coverage_profile

The tool first makes the coverage profile of each track: at each position, the sum of the values of the track’s fragments that overlap that position. Each position of a bin then adds its profile value, with a weight of 1 bp. A position where the profile is 0 is uncovered. With pooled resources, the profile of a group adds the fragments of every resource.

Two fragments in three bins

The run definition below bins two fragments, chr1:6-14 and chr1:11-18, on a 30-base chr1 with bins of 10. Each fragment adds the default value 1. The entries bin one resource under the three modes, so each entry gets a name to keep its tracks apart:

reference_genome: genome

bins:
  bin_size: 10

binners:
- fragment_score_binner:
    resource_query: frags/worked
    name: start
- fragment_score_binner:
    resource_query: frags/worked
    name: length
    aggregate: {mode: fragment_length}
- fragment_score_binner:
    resource_query: frags/worked
    name: profile
    aggregate: {mode: coverage_profile}
- fragment_score_binner:
    resource_query: frags/worked
    name: start
    aggregate: {aggregator: mean}
- fragment_score_binner:
    resource_query: frags/worked
    name: length
    aggregate: {mode: fragment_length, aggregator: mean}
- fragment_score_binner:
    resource_query: frags/worked
    name: length_bg
    aggregate: {mode: fragment_length, aggregator: mean, uncovered_value: 0}
- fragment_score_binner:
    resource_query: frags/worked
    name: profile
    aggregate: {mode: coverage_profile, aggregator: mean}
- fragment_score_binner:
    resource_query: frags/worked
    name: covered
    aggregate: {mode: coverage_profile, aggregator: mean, uncovered_value: null}

The run writes these tracks:

Track

mode

aggregator

uncovered_value

1-10

11-20

21-30

start:all:sum

fragment_start

sum

—

1.0

1.0

0.0

length:all:sum

fragment_length

sum

null

5.0

12.0

0.0

profile:all:sum

coverage_profile

sum

0

5.0

12.0

0.0

start:all:mean

fragment_start

mean

—

1.0

1.0

NaN

length:all:mean

fragment_length

mean

null

1.0

1.0

NaN

length_bg:all

fragment_length

mean

0

0.5

0.857 (12/14)

0.0

profile:all:mean

coverage_profile

mean

0

0.5

1.2

0.0

covered:all

coverage_profile

mean

null

1.0

1.5

NaN

How to read the table:

  • fragment_start puts each fragment in the bin of its start: one in 1-10 and one in 11-20.

  • fragment_length adds 5 bp of the first fragment to 1-10, and 4 bp of the first and 8 bp of the second to 11-20. The mean of the value 1 over these base pairs is 1.0. With uncovered_value: 0, the five uncovered bases of 1-10 and the two of 11-20 also add 0, so the means are 5/10 and 12/14.

  • coverage_profile makes the profile 1 on 6-10, 2 on 11-14, 1 on 15-18 and 0 elsewhere. With the default uncovered_value: 0, the mean of 11-20 is the mean depth over the full bin, 12/10. With uncovered_value: null, it is the mean depth over the eight covered bases, 12/8.

  • The sums of fragment_length and coverage_profile agree, because both add each covered base pair once for each fragment over it. The two modes differ in the weight of a base that two fragments cover: fragment_length counts it two times, and coverage_profile counts it one time, with the value 2.

Uncovered bases

aggregate.uncovered_value is the value of each base of a region that no fragment of the track covers. The tool adds it to the bin with a weight of 1 bp for each uncovered base. In coverage_profile mode, “uncovered” means a profile value of 0.

Mode

Default

Other values

fragment_start

not applicable

refused

fragment_length

null: the uncovered bases add nothing

any number

coverage_profile

0: an uncovered base has depth 0

null or any number

When the value is a number:

  • The tool fills only the bases of the region. A bin is already clipped to its region (see Coordinates and the grid), so the bin holds only bases of the region.

  • A chromosome that a resource does not have is uncovered throughout. Each of its bins holds the aggregate of the uncovered value over its bases: with uncovered_value: 2, a mean is 2 and a sum over a 10 bp bin is 20.

  • No bin is empty, so the empty-bin rule of Empty bins does not apply.

With null, the empty-bin rule applies. The value must be finite. NaN is refused, because the output file records null as NaN.

fragment_start refuses the key, because a bin in this mode reduces only the fragments that start in it:

binners[0].aggregate: uncovered_value does not apply in fragment_start mode, where a bin reduces only the fragments starting in it; give mode: fragment_length or coverage_profile, or drop uncovered_value

Aggregators in each mode

fragment_start accepts sum, count, mean, max, min, median and product. fragment_length and coverage_profile accept sum, mean, max, min and median. These two modes refuse two aggregators when the run definition is read:

  • count, because with a weight it counts base pairs, not fragments. For the number of covered base pairs, use aggregator: sum with value: {value: 1}:

    binners[0].aggregate: aggregator 'count' would count base pairs, not fragments, in fragment_length mode; for the covered base pairs give aggregator: sum with value: {value: 1}
    
  • product, because with a weight it raises a value to the power of a length:

    binners[0].aggregate: aggregator 'product' would weigh a value by its overlap as a power in coverage_profile mode; use one of max, mean, median, min, sum
    

The single-cell convention

A single-cell fragment resource that names its metadata the conventional way needs no group and no meta. The convention is:

Where

What

Name

fragment resource label

the id of the data_frame resource holding the cell metadata

cell_meta_resource_id

fragment resource label

this sample’s key in that table

sample_id

metadata table column

the sample key

sample_id

metadata table column

the cell barcode

barcode

metadata table column

the cell’s group

class

and the barcode is the fragment score cell. The default meta is therefore the table that cell_meta_resource_id names, restricted to the rows whose sample_id column equals the resource’s sample_id label — spelled out, {resource_label: cell_meta_resource_id, filter: [{column: sample_id, label: sample_id}]}.

All of it is checked before any work runs, with an error naming the entry and the resource: the label exists and names a data_frame resource the repository has; the table has the columns the entry uses; every filter label exists on the resource; no barcode occurs on two of one resource’s selected rows; and the selected rows name at least one group. A table that predates the convention is used through the explicit group and meta keys.

Defaults

Without group, grouping is on if and only if the resource carries the cell_meta_resource_id label: the label is the declaration, and nothing is guessed from the names of a resource’s scores. A labelled resource is grouped as with group: {}; any other resource is the one group all, so its track is named <resource id>:all. A resource that has cell and count scores but no label is binned as the single all track, and --dry-run warns about it.

Without value and aggregate, an entry counts the fragments that start in each bin, as value: {value: 1} and aggregate: {mode: fragment_start, aggregator: sum}, whatever scores the resource has. The read pairs behind each fragment are summed only on request, with value: {score_id: count}.

A pooled entry is one job reading one table, so its resources must agree: on whether they carry cell_meta_resource_id (when group is omitted), and on the one metadata table their meta names. Resources that disagree are an error listing both sides; give the setting explicitly, or pool: false.

Empty bins

A bin that no fragment of a group reaches holds the aggregator’s empty value: 0 for count and sum, which have an answer for no fragments, and NaN for mean, max, min, median and product, which have none. A chromosome a resource has no fragments on is treated exactly like empty bins, by the same rule, so one sample lacking a contig neither fails the run nor turns a pooled column into NaN.

This rule applies in fragment_start mode, and in the two other modes with uncovered_value: null. With a numeric uncovered_value, no bin is empty (see Uncovered bases).

Dropped fragments

With a grouping, a fragment reaches no track when its barcode is not among its resource’s selected rows, when its row’s group is empty, or — when grouping by a raw value — when its value is missing or is not one the score’s histogram lists. Such fragments are dropped, counted per resource and region, and the count is logged at INFO level (-v):

INFO:gain.binning.fragment_binner:sc/atac_fragments/donor1 chr1:1-1000: 2 fragments dropped, their barcode absent from the cell metadata rows or of no group

Grouping by a raw value

group: {group_score_id: S} makes one track per distinct value of the str score S. Grouped by the cell score, that is one track per cell: a per-cell matrix.

- fragment_score_binner:
    resource_query: "sc/atac_fragments/*"
    group: {group_score_id: cell}

The values are not discovered by reading the data: they are read, when the entry is resolved, from the score’s full categorical histogram in the resource’s statistics. The statistics must therefore be built (grr_manage resource-stats) and, in a repository that keeps them in DVC, pulled; a resource without the full histogram, or with a truncated one, is an error naming it. A score whose values pass the default categorical histogram’s limit of 100 values keeps the full histogram only when its configuration declares histogram: {type: categorical}. The grouping score may not be the value score.

Unpooled, a track’s group is the bare value. Pooled, the same barcode recurs in every 10x sample while naming a different cell, so the group is prefixed with the resource’s sample_id label, <sample_id>:<value>. Every pooled resource must therefore carry a sample_id label of one value, without a :, and no two pooled resources may share one; any of these is an error naming the resources, which suggests pool: false.

Mind the size before running: a per-cell matrix has one column per cell, which is thousands of columns per sample, and the full histogram of a cell score is itself tens of megabytes per resource. --dry-run reports the group count read from the histograms, so the width of the matrix is known before it is built.

Local metadata files

meta: {file_name: PATH} reads the metadata table from a local file instead of a data_frame resource, with the same column conventions. A relative PATH is read from the directory of the run definition, not from the current directory:

- fragment_score_binner:
    resource_query: "sc/atac_fragments/*"
    meta:
      file_name: cells.csv
      filter:
      - {column: sample_id, label: sample_id}

A run that reads a local file can be repeated only where that file is, so --dry-run warns about it, and the output file records which file it was: its absolute path, size and modification time, in the metadata_files, metadata_file_sizes and metadata_file_mtimes root attributes (see The output file).

Track names

Every track has a name, and names in the output are always unique.

  • A position-score track is named by its resource id.

  • A fragment track of an unpooled entry is named <resource id>:<group>; in the plain case, with no grouping, <resource id>:all.

  • A fragment track of a pooled entry is named <resource_query>:<group>, the query string as written in the entry, which is the one name every pooled entry has — for example sc/atac_fragments/*:ODC. Grouped by a raw value, the group carries its sample, so a pooled track is <resource_query>:<sample_id>:<value>.

An entry’s name key replaces the base of these names. Without name, the names are as above.

Entry

Without name

With name: N

fragment, pool: true

<resource_query>:<group>

N:<group>

fragment, pool: false

<resource id>:<group>

N/<resource id>:<group>

position score

<resource id>

N/<resource id>

An unpooled entry keeps the resource id in each name, so each resource still has its own tracks. For example, an entry with resource_query: "sc/atac_fragments/*" and name: depth gives the pooled tracks depth:ODC and depth:PVALB, and the same entry with pool: false gives depth/sc/atac_fragments/donor1:ODC and its siblings.

Tracks follow entry order, and within an entry the order of resource ids and then of groups. When the same name and group occur 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 set gets :<aggregator> appended, whichever entry came first, so hg38/scores/phyloP100way:mean and hg38/scores/phyloP100way:max sit side by side, as do sc/atac_fragments/*:ODC:sum and sc/atac_fragments/*:ODC:mean. Two entries that would still produce one and the same track (same name, group and aggregator) are refused, naming both entries: nothing in the file could tell the two columns apart. This happens when two entries over one query differ only in their mode, because the :<aggregator> suffix cannot tell modes apart. The refusal tells you to add name to one of the entries:

binners[4] and binners[5] both produce the track 'sc/atac_fragments/*:ODC:mean'; a resource may be binned once per aggregator and name -- add name to one of them to tell the tracks apart

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.

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

Read it with gain.binning.h5_reader.read_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

Read it with gain.binning.h5_reader.read_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

Read it with gain.binning.h5_reader.read_tracks(). A compound dataset of shape (n_tracks,), one row per track, with the same fields for every binner kind. binner_type, resource_ids, selector and binner_def are variable-length UTF-8 strings. binner_index and track_index are int64.

binner_type

The kind of the entry that produced the track, without the _binner suffix: position_score or fragment_score.

binner_index

The position of the entry in the binners list of the run definition, counted from 0.

track_index

The position of the track among all the tracks of that entry, in output order, counted from 0. The pair (binner_index, track_index) is unique in a file.

resource_ids

The resources the track was computed from, comma-separated in resource-id order: one id for a position-score track and for an unpooled fragment track, every matched id for a pooled one.

selector

What picks the track within its resources: the score id of a position-score track, the group of a fragment track.

binner_def

The definition of the track that the other fields do not carry, as canonical JSON (sorted keys, no extra whitespace). A missing value is JSON null. A position-score definition holds aggregator, none_value_replacement and name. A fragment definition holds aggregator, mode (fragment_start, fragment_length or coverage_profile), uncovered_value, score_id (null for a track that counts fragments), resource_query, pool and name. name is the optional name of the entry.

The file has no track name. The name is in memory while the tool runs, where it keeps the tracks apart and names the chunk files.

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

tracks_layout

2, the layout of /tracks described above. A file without this attribute has an older layout.

reference_genome

The resource id of the reference genome the grid was built on. A file without tracks_layout keeps the older name input_reference_genome for this attribute.

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.

metadata_files

Only when the run read a local metadata file: the absolute path of each local file the run read.

metadata_file_sizes

Only with metadata_files: each file’s size in bytes, in the same order.

metadata_file_mtimes

Only with metadata_files: each file’s modification time, UTC, in ISO 8601, in the same order.

Note

The /tracks layout changed after the 2026.9 releases. A file written by an earlier release has the fields name, resource_ids, group, score_id, aggregator and more, or in the oldest files a resource_id field. A file written now has the root attribute tracks_layout with the value 2 and the fields described above. Tell the layouts apart with h5.attrs.get("tracks_layout"): absent means an older layout. Do not sniff a field, and do not use the gain_version root attribute: a development build can report an old version and already write the new layout. In the same change the root attribute input_reference_genome became reference_genome; read the name that matches tracks_layout. The tool does not convert files written by earlier releases.

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. For a position-score track, 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. A count track is the exception: an uncovered bin is 0.0, because a count of nothing is zero. For a fragment track, an empty bin holds the aggregator’s empty value (see Empty bins).

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

gain.binning.h5_reader.read_tracks() and gain.binning.h5_reader.read_bins() read /tracks and /bins as pandas DataFrames. Each takes a path or an open h5py.File. The columns come from the file, the text fields are str, and the index is the position in /values. read_tracks raises ValueError for a file with a missing or unsupported tracks_layout. Both functions raise ValueError for a table that is missing, not a dataset, or not compound. A path that is missing or not an HDF5 file raises the OSError of h5py. gain.binning.h5_reader.read_values() reads /values as a numpy float64 matrix of shape (n_bins, n_tracks), not as a DataFrame. It reads the whole matrix into memory, and it raises ValueError for a /values that is missing, not a dataset, or not two-dimensional. gain.binning.h5_reader.track_labels() returns one display label for each track:

from gain.binning.h5_reader import (
    read_bins, read_tracks, read_values, track_labels)

values = read_values("binning_run.h5")  # numpy matrix, bins x tracks
bins = read_bins("binning_run.h5")
labels = track_labels(read_tracks("binning_run.h5"))

column = labels.index("hg38/scores/phyloP100way:mean")
on_chr1 = (bins["chrom"] == "chr1").to_numpy()
print(values.shape)
print(values[on_chr1, column][:5])

values is a numpy matrix, not a DataFrame. Row i belongs to row i of bins and column j to labels[j], which is the label of row j of the tracks. The label of a position-score track is its resource id and its aggregator, for example hg38/scores/phyloP100way:mean. The label of a fragment track is its resource, the group and the aggregator, for example frags/s1:B:sum. A name of the entry prefixes the resource id of a position-score track and of an unpooled fragment track, and replaces the resource query in the label of a pooled fragment track. Tracks with equal labels get #<position> appended, where the position is the column of the track.

To print the tables without Python, use Printing a binning file with binning_tool_view.

With h5py alone:

import json

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: binner_type, ...
    attrs = dict(h5.attrs)

assert attrs["tracks_layout"] == 2  # absent: a file of an older layout
print(values.shape, attrs["bin_size"], attrs["coordinates"])
print(tracks["binner_type"][:3], tracks["selector"][:3])
print(json.loads(tracks["binner_def"][-1]))
print(bins[:3])

A worked example

The run definition below bins the conservation scores of the public IossifovLab GRR over two megabases of chromosome 21, together with the structural variants of gnomAD v4.1 starting in each bin (binning_run.yaml):

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
- fragment_score_binner:
    resource_query: hg38/cnv_collections/gnomAD.v4.1_Genome_SV
    pool: false
    group: {group_score_id: deletion_duplication}

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.

The third entry bins a fragment_score resource: the gnomAD v4.1 genome structural-variant collection, whose records are intervals. It is grouped by the raw value of its deletion_duplication score, whose statistics list two values, deletion and duplication, so it makes two tracks. Without value and aggregate, each track counts the variants that start in each bin. A single resource carries no sample_id label, so the entry sets pool: false and the groups are the bare values.

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
  hg38/cnv_collections/gnomAD.v4.1_Genome_SV:deletion       hg38/cnv_collections/gnomAD.v4.1_Genome_SV              sum
  hg38/cnv_collections/gnomAD.v4.1_Genome_SV:duplication    hg38/cnv_collections/gnomAD.v4.1_Genome_SV              sum
regions: 1
bins: 196
tasks: 10
binners[2]: fragment_score_binner
    resources: hg38/cnv_collections/gnomAD.v4.1_Genome_SV
    tables: none
    groups: 2
    tracks: 2
    mode: fragment_start
    uncovered_value: null

A job is what one task binds to, and a task is one job over one bundle of regions (see Parallelism). Nine position-score tracks are nine jobs, and the fragment entry’s two tracks are one job, so the single region makes ten tasks.

Then run it:

binning_tool binning_run.yaml -j 4

This writes binning_run.h5 beside the run definition: a 196 × 11 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, where each column is labelled resource_ids:selector:aggregator, gives:

                         hg38/scores/phyloP100way:phyloP100way:mean  hg38/scores/phyloP100way: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

and the two fragment tracks, selected the same way, count the deletions and duplications that start in each bin — 826 and 179 over the whole window:

                         hg38/cnv_collections/gnomAD.v4.1_Genome_SV:deletion:sum  hg38/cnv_collections/gnomAD.v4.1_Genome_SV:duplication:sum
chrom start    end
chr21 20000000 20008960                                                                     2.0                                                         0.0
      20008961 20019200                                                                     2.0                                                         3.0
      20019201 20029440                                                                     4.0                                                         0.0
      20029441 20039680                                                                     1.0                                                         0.0
      20039681 20049920                                                                     4.0                                                         0.0

Command line tools

GAIn provides two command-line tools for binning:

  • binning_tool: bin the scores that a run definition names into one HDF5 file

  • binning_tool_view: print the tables of a binning_tool file as tab-separated text