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_browseuses: an exact resource id, or a glob over ids such ashg38/scores/*, optionally followed by a bracketed label filter. Only resources of the entry’s kind are taken from the matches —position_scoreresources for aposition_score_binner,fragment_scoreresources for afragment_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 ofgrr_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 oncegrr_manage repo-indexhas been run on it. Asearch_termon 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
aggregatorthe score itself declares in its resource configuration, so a conservation score and a coverage track are each reduced the way their authors intended. Settingaggregatoron an entry overrides that for every resource the entry matches. Only aggregators that produce a number are accepted, and only numeric (intorfloat) scores can be binned; a string-typed score, or an aggregator such asjoinorlistthat builds a string or a list, is refused at parse time. The numeric aggregators aremax,min,mean,median,sum,productandcount. 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
NaNin the output, so missing data cannot masquerade as a measurement in a correlation. The exception isaggregator: count: an uncovered bin is0.0, because a count of nothing is zero. Set a replacement, typically0.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, defaulttrue)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. Withfalse, 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 constantNAME.{cell_score_id, cell_meta_column, group_meta_column}— one track per group of a metadata table, whichmetanames.cell_score_idis the fragment score carrying a fragment’s cell barcode,cell_meta_columnthe table column holding the barcode andgroup_meta_columnthe column holding the cell’s group. Each key is optional and defaults to the convention’s —cell,barcodeandclass— sogroup: {}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, withpool: 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 thestrscoreS; see Grouping by a raw value.metadoes 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 (intorfloat) fragment scoreS, which every matched resource must have;{value: V}is the constant numberVfor 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_lengthorcoverage_profile. See Modes.aggregator— the reduction,sumwhen 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 ornull. Onlyfragment_lengthandcoverage_profiletake it, and each has its own default. See Uncovered bases.
With the default
value,suminfragment_startmode 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— thedata_frameresourceID;resource_label: LABEL— thedata_frameresource that the fragment resource’s labelLABELnames;file_name: PATH— a local file; see Local metadata files.file_format(csv,tsvorexcel) andfile_separatorgo with it and only with it. Withoutfile_formatthe format follows the file’s suffix:.tsvis read astsv,.xlsand.xlsxasexcel, and anything else ascsv. Only.xlsxis readable out of the box: a legacy.xlsfile also needs thexlrdpackage, 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 columnCequal to the literalV, or{column: C, label: L}, the columnCequal to the fragment resource’s labelL— so one filter selects each sample’s own rows. A row must satisfy every conjunct; nofilterselects every row. A label whose value is a list cannot be a filter operand.metagoes only with the metadata form ofgroup(or with nogroupon 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_lengthA 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_profileThe 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 |
|
|
|
|
|
|
|---|---|---|---|---|---|---|
|
|
|
— |
1.0 |
1.0 |
0.0 |
|
|
|
|
5.0 |
12.0 |
0.0 |
|
|
|
|
5.0 |
12.0 |
0.0 |
|
|
|
— |
1.0 |
1.0 |
NaN |
|
|
|
|
1.0 |
1.0 |
NaN |
|
|
|
|
0.5 |
0.857 (12/14) |
0.0 |
|
|
|
|
0.5 |
1.2 |
0.0 |
|
|
|
|
1.0 |
1.5 |
NaN |
How to read the table:
fragment_startputs each fragment in the bin of its start: one in1-10and one in11-20.fragment_lengthadds 5 bp of the first fragment to1-10, and 4 bp of the first and 8 bp of the second to11-20. The mean of the value 1 over these base pairs is 1.0. Withuncovered_value: 0, the five uncovered bases of1-10and the two of11-20also add 0, so the means are 5/10 and 12/14.coverage_profilemakes the profile 1 on6-10, 2 on11-14, 1 on15-18and 0 elsewhere. With the defaultuncovered_value: 0, the mean of11-20is the mean depth over the full bin, 12/10. Withuncovered_value: null, it is the mean depth over the eight covered bases, 12/8.The sums of
fragment_lengthandcoverage_profileagree, 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_lengthcounts it two times, andcoverage_profilecounts 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 |
|---|---|---|
|
not applicable |
refused |
|
|
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, ameanis 2 and asumover 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, useaggregator: sumwithvalue: {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 |
|
fragment resource label |
this sample’s key in that table |
|
metadata table column |
the sample key |
|
metadata table column |
the cell barcode |
|
metadata table column |
the cell’s group |
|
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 examplesc/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 |
With |
|---|---|---|
fragment, |
|
|
fragment, |
|
|
position score |
|
|
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:
/valuesRead it with
gain.binning.h5_reader.read_values().float64, shape(n_bins, n_tracks).NaNwhere a bin has no data. Stored in row blocks and gzip-compressed, so reading every track for one chromosome is one contiguous read and theNaN- and zero-heavy tracks compress well./binsRead it with
gain.binning.h5_reader.read_bins(). A compound (structured) dataset of shape(n_bins,)with fieldschrom(fixed-length bytes, sized to the longest chromosome name in the run),startandend(int64, 1-based inclusive)./tracksRead 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,selectorandbinner_defare variable-length UTF-8 strings.binner_indexandtrack_indexareint64.binner_typeThe kind of the entry that produced the track, without the
_binnersuffix:position_scoreorfragment_score.binner_indexThe position of the entry in the
binnerslist of the run definition, counted from 0.track_indexThe 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_idsThe 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.
selectorWhat picks the track within its resources: the score id of a position-score track, the group of a fragment track.
binner_defThe 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 holdsaggregator,none_value_replacementandname. A fragment definition holdsaggregator,mode(fragment_start,fragment_lengthorcoverage_profile),uncovered_value,score_id(nullfor a track that counts fragments),resource_query,poolandname.nameis the optionalnameof 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 |
|---|---|
|
|
|
The resource id of the reference genome the grid was built on.
A file without |
|
The bin width in base pairs. |
|
The binned regions, in row order, each as |
|
|
|
The GAIn version that wrote the file. |
|
The UTC time the file was written, in ISO 8601. |
|
Only when the run read a local metadata file: the absolute path of each local file the run read. |
|
Only with |
|
Only with |
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_toolfile as tab-separated text