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. .. code-block:: bash 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. .. _binning-run-definition: The run definition ------------------ A run definition has three parts: the reference genome, the ``bins`` block, and a list of ``binners`` entries. .. code-block:: yaml 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 :doc:`/grr`). 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 :doc:`/grr`). 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 ``/``. 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. .. _binning-fragment-score-entries: 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: .. code-block:: yaml - 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: .. code-block:: yaml - 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: .. code-block:: yaml 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: .. list-table:: :header-rows: 1 :widths: 22 22 10 16 10 10 10 * - 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. .. list-table:: :header-rows: 1 :widths: 25 40 35 * - 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: .. code-block:: text 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}``: .. code-block:: text 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: .. code-block:: text 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: .. list-table:: :header-rows: 1 :widths: 30 40 30 * - 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 ``: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``): .. code-block:: text 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. .. code-block:: yaml - 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, ``:``. 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. .. _binning-local-metadata-files: 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: .. code-block:: yaml - 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 ``:``; in the plain case, with no grouping, ``:all``. * A fragment track of a pooled entry is named ``:``, 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 ``::``. An entry's ``name`` key replaces the base of these names. Without ``name``, the names are as above. .. list-table:: :header-rows: 1 :widths: 30 35 35 * - Entry - Without ``name`` - With ``name: N`` * - fragment, ``pool: true`` - ``:`` - ``N:`` * - fragment, ``pool: false`` - ``:`` - ``N/:`` * - position score - ```` - ``N/`` 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 ``:`` 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 ``:`` suffix cannot tell modes apart. The refusal tells you to add ``name`` to one of the entries: .. code-block:: text 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 :func:`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 :func:`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 :func:`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: .. list-table:: :header-rows: 1 :widths: 30 70 * - 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 --------------------- :func:`gain.binning.h5_reader.read_tracks` and :func:`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``. :func:`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. :func:`gain.binning.h5_reader.track_labels` returns one display label for each track: .. code-block:: python 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 ``#`` appended, where the position is the column of the track. To print the tables without Python, use :doc:`cli/binning_tool_view`. With ``h5py`` alone: .. code-block:: python 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 (:download:`binning_run.yaml `): .. literalinclude:: files/binning_run.yaml :language: yaml 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: .. code-block:: bash binning_tool binning_run.yaml --dry-run .. code-block:: text 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 :ref:`binning-tool-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: .. code-block:: bash 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 :doc:`/gain_getting_started_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: .. code-block:: text 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: .. code-block:: text 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: * :doc:`binning_tool `: bin the scores that a run definition names into one HDF5 file * :doc:`binning_tool_view `: print the tables of a ``binning_tool`` file as tab-separated text