""":class:`PositionScore` -- one value per genomic position.
The kind whose records carry no reference or alternative allele, so a
position is the whole key. Adds the position-run reads and the binned and
aggregated region queries built on them.
"""
from __future__ import annotations
import copy
from collections.abc import Generator, Iterator, Sequence
from typing import (
Any,
ClassVar,
)
from gain.genomic_resources.repository import (
GenomicResource,
)
from gain.genomic_resources.score_def import (
GenomicScoreDef,
ScoreValue,
)
from gain.utils.regions import (
calc_bin_begin,
calc_bin_end,
calc_bin_index,
)
from ..aggregators import (
AGGREGATOR_SCHEMA,
Aggregator,
PositionScoreAggregationQuery,
)
from .aggregation import (
AggregateValue,
build_region_aggregator,
distinct_score_ids,
resolve_aggregator_name,
score_def_for,
)
from .base import GenomicScore
from .records import clip_span
[docs]
class PositionScore(GenomicScore):
"""Position-based genomic score resource.
A position score gives one or more scores for a genomic position or a
run of positions. Unlike an allele score, it ignores the reference and
alternative alleles. The data lives in a tabix-indexed file, a bigWig
file or an in-memory table.
Typical resources are:
- conservation scores, such as ``phastCons`` and ``phyloP``
- mappability scores
- GC content
- recombination rates
- any metric that depends only on the genomic position
A score is built closed. Call ``open()`` before the first read.
Examples
--------
.. code-block:: python
import gain.grr as grr
repo = grr.build_genomic_resource_repository()
score = grr.build_position_score_from_resource_id(
"hg38/scores/phastCons100way", repo)
with score.open() as score:
values = score.get_scores_at_position("chr1", 12345)
for begin, end, row in score.fetch_region_segments_scores(
"chr1", 10000, 20000):
print(f"{begin}-{end}: {row}")
Attributes:
resource (GenomicResource): The underlying genomic resource.
resource_id (str): The identifier of the resource.
config (dict): The validated configuration of the score.
score_definitions (dict[str, GenomicScoreDef]): The definition of
each score, by score id.
HISTOGRAM_COUNT_UNIT (str): The unit of the histogram counts. For a
position score it is ``"base pairs"``.
DEFAULT_AGGREGATORS (dict[str, str | None]): The aggregator name
that a region read uses for each value type, when the
configuration of the score names none. A value of ``None``
means that the type takes no default.
"""
# Aggregating the values of a region is the resource's job since
# gain#1131: ``get_scores_in_region_agg`` reduces a region to one value
# per query, and ``gain.annotation.position_score_annotator`` asks for
# that rather than folding records of its own. What the kind
# contributes to the reduction is ``record_weight``: how many queried
# bases a record covers, and so how many times its value counts.
HISTOGRAM_COUNT_UNIT: ClassVar[str] = "base pairs"
# A region of positions reduces by ``mean``: each position's value counts
# once per base pair it covers (see :meth:`record_weight`).
DEFAULT_AGGREGATORS: ClassVar[dict[str, str | None]] = {
"float": "mean",
"int": "mean",
"str": "list",
"bool": None,
}
def __init__(self, resource: GenomicResource):
if resource.get_type() != "position_score":
raise ValueError(
"The resource provided to PositionScore should be of "
f"'position_score' type, not a '{resource.get_type()}'")
super().__init__(resource)
# The only kind whose weight is not 1. That there is exactly one value
# per position -- what a position score PROMISES -- is not stated here
# but in the rules this kind is registered under in
# :mod:`gain.genomic_resources.statistics.record_validation`, the only
# places that enforce it. This is a MEASURE.
#
# Elementwise, as the base requires: handed the position columns of a
# whole batch, the same expression answers that batch's weights.
@classmethod
def record_weight(cls, left: int, right: int) -> int:
"""Return how many times a record counts: once per base pair.
This is a hook for the aggregation engine, not a user call.
Parameters
----------
left : int
The first position of the record.
right : int
The last position of the record.
Returns
-------
int
The number of base pairs the record covers.
:meta private:
"""
return right - left + 1
[docs]
@staticmethod
def get_schema() -> dict[str, Any]:
"""Return the configuration schema of a position score.
The schema is the one of :class:`GenomicScore`, with an
``aggregator`` key added to the definition of each score.
Returns
-------
dict[str, Any]
A new Cerberus schema dictionary, which the caller may modify.
"""
schema = copy.deepcopy(GenomicScore.get_schema())
scores_schema = schema["scores"]["schema"]["schema"]
scores_schema["aggregator"] = AGGREGATOR_SCHEMA
return schema
# -- The logical read plane (#727) -------------------------------------
#
# On this plane a position score is a function from a genomic position
# to a record of named score values, defined on the per-position
# expansion: a gap in coverage is not an absence, it is a run of
# positions whose value is ``None``, and it counts. It is computed by
# walking segments -- only the two region reads materialise a value per
# position; the aggregating and binning reads fold one
# ``(values, run_length)`` pair per segment, so cost stays proportional
# to record count.
def _plane_read_defs(
self, chrom: str, start: int, end: int,
scores: Sequence[str] | None,
) -> list[GenomicScoreDef]:
"""Refuse what this region read cannot serve, and resolve its scores.
The preamble every read that takes a MANDATORY region shares: the
span guard, then the request refusals, then the definitions in the
order the caller asked for them. Here rather than repeated per
read so that the two cannot drift in what they refuse -- the deeper
fold, of :meth:`_guard_region_span` into :meth:`_region_read_defs`,
is declined in that method's docstring because it would also refuse
``fetch_*`` requests on all three kinds, which this does not touch.
The DEFINITIONS are what travels on, not the ids resolved out of
them: a read that has resolved its request carries the answer down
to :meth:`~.base.GenomicScore._fetch_segments_for_defs` rather than
handing back ids for the segment read to resolve a second time
(gain#1282). A caller wanting the ids reads them off the defs.
The aggregating and binning reads do NOT share it: their scores
come from ``distinct_score_ids`` over the queries, not from
``scores``, and a contig this score never mentions is uncovered
there rather than refused (gain#1211).
"""
self._guard_region_span(start, end)
return self._region_read_defs(chrom, scores)
def _fetch_segments_for_defs(
self, chrom: str, start: int, end: int,
score_defs: list[GenomicScoreDef],
) -> Generator[
tuple[int, int, list[ScoreValue]], None, None]:
""":meth:`~.base.GenomicScore.fetch_region_segments_scores`,
pre-resolved.
The same read, from the same records, through the same per-kind
transform (:meth:`~.base.GenomicScore._score_segments`) -- entered
one step lower, with the definitions this plane already holds
instead of ids for that method to resolve a second time.
A read that reaches here has run :meth:`_plane_read_defs`, or the
aggregating pair's own resolver, and would otherwise pay the whole
refusal preamble twice: the open check, the contig screen and the id
resolution. The second answer can only agree with the first, so it
was waste, and gain#1282 removed it. The contig half was the
expensive one when that landed -- a walk of the ordered contig list,
O(contigs) on a tabix table (gain#1173); it is a set lookup since
gain#1304, and the id resolution is now the larger of the two. The
repeat is skipped either way, because a second answer that can only
agree is waste at any price.
No refusal is weakened, because this is not the only one that runs:
every caller resolves eagerly first and reaches here only after.
What is skipped is a REPEAT, never the first answer.
On :class:`PositionScore` rather than on the base, though it uses
only base methods: a door that skips the eager refusals may be
opened only by a read that has already made them, and this plane's
reads are the ones that have. Same placement argument as
:meth:`~.base.GenomicScore._read_defs_for_any_contig`, which is on
the base and has to say in prose who may call it; this does not
need to. The fragment plane wants nothing here -- it never
resolves before reading, so it has no definitions to carry and
pays the preamble exactly once through the public entry.
Like the method it shadows, a plain function returning a generator:
what the caller resolved was refused at ITS call, and nothing moved
to the first ``next()``.
No ``score_filter``, unlike the method it shadows: none of this
plane's reads select records. One that grew a filter would take
the parameter here and forward it to
:meth:`~.base.GenomicScore.fetch_records`, where it belongs
(gain#1272).
"""
return self._score_segments(
self.fetch_records(chrom, start, end), score_defs)
def _position_runs(
self, chrom: str, start: int, end: int,
score_defs: list[GenomicScoreDef],
) -> Generator[tuple[list[ScoreValue] | None, int], None, None]:
"""Yield ``(values, run_length)`` runs tiling ``[start, end]``.
The per-position expansion, run-length encoded: every position of
the region belongs to exactly one run, uncovered positions to a
``None`` run. Where two records cover one position the FIRST
answers -- a later record contributes only the positions the
earlier one did not -- so the total run length always equals the
region width, and accumulated weight can never exceed it.
This region read fetched from the table; :meth:`_runs_from_segments`
is the same encoding over a segment stream a caller supplies.
Takes DEFINITIONS, not ids: every caller has already resolved the
request -- through :meth:`_plane_read_defs`, or through the
aggregating reads' own resolver -- and
:meth:`_fetch_segments_for_defs` is the segment read that does not
resolve it a second time (gain#1282).
"""
return self._runs_from_segments(
self._fetch_segments_for_defs(chrom, start, end, score_defs),
start, end)
@staticmethod
def _runs_from_segments(
segments: Iterator[tuple[int, int, list[ScoreValue]]],
start: int, end: int,
) -> Generator[tuple[list[ScoreValue] | None, int], None, None]:
"""Run-length encode ``segments`` over ``[start, end]``.
:meth:`_position_runs` as a function OF a segment stream, so what
the region is read FROM is the caller's to choose. An empty stream
is a whole region of uncovered positions -- one ``(None, width)``
run out of the tail below -- which is what lets
:meth:`_aggregating_runs` answer a contig the score never mentions
without a second spelling of that shape (gain#1211).
"""
cursor = start
for left, right, values in segments:
span = clip_span(left, right, start, end)
if span is None:
continue
left, right = span
if right < cursor:
continue
# cursor >= start, so this also subsumes the clip's left clamp.
left = max(left, cursor)
if left > cursor:
yield None, left - cursor
yield values, right - left + 1
cursor = right + 1
if cursor <= end:
yield None, end - cursor + 1
@staticmethod
def _expand_position_runs(
runs: Iterator[tuple[list[ScoreValue] | None, int]],
n_scores: int,
) -> Generator[tuple[ScoreValue | None, ...], None, None]:
"""Expand run-length encoded runs to one tuple per position."""
for values, length in runs:
row: tuple[ScoreValue | None, ...] = (
tuple(values) if values is not None else (None,) * n_scores)
for _ in range(length):
yield row
[docs]
def get_score_at_position(
self, chrom: str, pos: int,
score: str | None = None,
) -> ScoreValue | None:
"""Return the value of one score at one position.
This is the singular form of :meth:`get_scores_at_position`.
Parameters
----------
chrom : str
The chromosome name.
pos : int
The 1-based position.
score : str | None
The id of the score. ``None`` is accepted only when the
resource defines exactly one score.
Returns
-------
ScoreValue | None
The value of the score. ``None`` means that no record covers
the position, or that the record has no value for the score.
Raises
------
ValueError
If the score is closed, the chromosome is unknown, the score id
is unknown, or ``score`` is ``None`` and the resource defines
more than one score.
"""
return self.get_scores_at_position(
chrom, pos, [self._resolve_single_score(score)])[0]
[docs]
def get_scores_at_position(
self, chrom: str, pos: int,
scores: Sequence[str] | None = None,
) -> tuple[ScoreValue | None, ...]:
"""Return the values of the scores at one position.
Parameters
----------
chrom : str
The chromosome name.
pos : int
The 1-based position.
scores : Sequence[str] | None
The ids of the scores. ``None`` asks for every score that the
resource defines, in definition order.
Returns
-------
tuple[ScoreValue | None, ...]
One value for each requested score, in the order of ``scores``.
A ``None`` value means that no record covers the position, or
that the record has no value for that score. The tuple is
empty when ``scores`` is an empty sequence.
Raises
------
ValueError
If the score is closed, the chromosome is unknown, or a score
id is unknown.
"""
# Read straight off _position_runs rather than through the region
# read's per-position expansion: a one-position region is exactly
# ONE run -- no run is ever yielded empty, and the run lengths sum
# to the region width -- so there is nothing to expand and nothing
# to index a position out of.
#
# The run is DRAINED rather than abandoned; that is simpler to
# write, not load-bearing, since gain#1120 moved the tabix buffer
# prune into a ``finally`` (``fragment.py`` cross-references this
# paragraph, and ``test_a_walk_of_point_reads_leaves_the_tabix_
# buffer_pruned`` pins the drained half).
#
# Going straight to the runs rather than delegating to
# get_scores_in_region is worth ~4% on this read, which the
# position annotator pays once per substitution.
score_defs = self._plane_read_defs(chrom, pos, pos, scores)
# The single-element unpack both DRAINS the generator and asserts
# the one-run invariant the docstring claims; a second run would
# raise here rather than be silently dropped.
(values, _), = self._position_runs(chrom, pos, pos, score_defs)
return (
tuple(values) if values is not None
else (None,) * len(score_defs))
[docs]
def get_score_in_region(
self, chrom: str, start: int, end: int,
score: str | None = None,
) -> Generator[ScoreValue | None, None, None]:
"""Yield the value of one score at each position of a region.
This is the singular form of :meth:`get_scores_in_region`.
Parameters
----------
chrom : str
The chromosome name.
start : int
The first 1-based position of the region.
end : int
The last 1-based position of the region.
score : str | None
The id of the score. ``None`` is accepted only when the
resource defines exactly one score.
Yields
------
ScoreValue | None
One value for each position, in position order. A ``None``
value means that no record covers the position, or that the
record has no value for the score.
Raises
------
ValueError
If the score is closed, the region is not valid, the chromosome
is unknown, or the score id is unknown. The function raises
before the first value, not during the iteration.
"""
rows = self.get_scores_in_region(
chrom, start, end, [self._resolve_single_score(score)])
return (row[0] for row in rows)
[docs]
def resolve_aggregation_queries(
self, queries: Sequence[PositionScoreAggregationQuery],
) -> list[tuple[str, str, ScoreValue]]:
"""Resolve each query to a score id, an aggregator name and a value.
Use it to find out whether a query is answerable without reading.
The position score annotator calls it once when the pipeline loads.
Parameters
----------
queries : Sequence[PositionScoreAggregationQuery]
The queries. Each names a score, an aggregator, and a
replacement value for the uncovered positions.
Returns
-------
list[tuple[str, str, ScoreValue]]
One tuple for each query, in the order of ``queries``: the id
of the score, the name of the aggregator, and the
``none_value_replacement`` of the query. The list is empty
when ``queries`` is empty.
Raises
------
ValueError
If a score id is unknown, a ``none_value_replacement`` has a
type that the score cannot take, or a query names no aggregator
and the value type of the score has no default one.
"""
# A query asks the same two questions every aggregation request
# does -- which score, and what reduces it -- so they are asked
# where they are answered for every surface, in :mod:`.aggregation`
# (score_def_for, resolve_aggregator_name).
#
# What a query asks BESIDES is the third: a none_value_replacement
# must be of a type the score can mean, following
# ``validate_aggregator``'s precedent. It is judged BETWEEN the
# other two -- after the score is known, since its value type is
# what judges the replacement, and before an aggregator is looked
# for, so that a query wrong in both ways is answered about the
# value it named rather than the one it left out. That order is a
# decision and not an accident of composition; it is pinned by
# ``test_a_query_invalid_several_ways_reports_the_first_ground``.
#
# Public, and stopping at the NAME, because asking whether a query
# is answerable is a question a caller may have without wanting to
# read (gain#1131). ``PositionScoreAnnotator`` asks it once when
# the pipeline loads, so an attribute naming no aggregator for a
# ``bool`` score is refused there rather than on the first region
# that reaches it. Building the accumulators is the READ's
# business -- ``_resolve_aggregation_queries`` adds them, per call,
# which is what keeps each read's accumulators its own and an
# annotator stateless.
#
# It lives on this kind rather than on GenomicScore because of its
# middle step: ``none_value_replacement`` is a field only a
# ``PositionScoreAggregationQuery`` carries, a position score being
# the only kind with an uncovered position to speak for. Should a
# fragment or an allele score come to want a resolver of its own
# (gain#1124, gain#1132), the score_def_for /
# resolve_aggregator_name pair is the part that generalises -- both
# already live in :mod:`.aggregation`, kind-neutral, for that
# reason -- and the replacement validation is the part that does
# not.
resolved = []
for query in queries:
score_def = score_def_for(
query.score,
score_definitions=self.score_definitions,
resource_id=self.resource_id)
self._validate_none_value_replacement(
query.score, score_def.value_type,
query.none_value_replacement)
aggregator = resolve_aggregator_name(
query.aggregator, score_def,
resource_id=self.resource_id)
resolved.append((
query.score, aggregator, query.none_value_replacement))
return resolved
def _resolve_aggregation_queries(
self, queries: Sequence[PositionScoreAggregationQuery],
) -> list[tuple[str, Aggregator, ScoreValue]]:
"""The read's form of :meth:`resolve_aggregation_queries`.
The same triples with the aggregator NAME replaced by a freshly
built accumulator -- one per call, never shared between reads.
"""
return [
(
score_id,
build_region_aggregator(
score_id, aggregator,
resource_id=self.resource_id),
none_value_replacement,
)
for score_id, aggregator, none_value_replacement
in self.resolve_aggregation_queries(queries)
]
# Which python types a none_value_replacement may have per score value
# type. A ``bool`` is deliberately not a valid int or float
# replacement, exactly as a bool-typed score is not a numeric one.
_NONE_VALUE_REPLACEMENT_TYPES: ClassVar[dict[str, tuple[type, ...]]] = {
"float": (int, float),
"int": (int,),
"str": (str,),
"bool": (bool,),
}
def _validate_none_value_replacement(
self, score_id: str, value_type: str | None,
none_value_replacement: ScoreValue,
) -> None:
"""Refuse a none_value_replacement of a type the score cannot mean."""
if none_value_replacement is None or value_type is None:
return
allowed = self._NONE_VALUE_REPLACEMENT_TYPES.get(value_type, ())
if isinstance(none_value_replacement, bool) and bool not in allowed:
pass
elif isinstance(none_value_replacement, allowed):
return
raise ValueError(
f"none_value_replacement {none_value_replacement!r} for score "
f"{score_id!r} of resource {self.resource_id!r} does not "
f"match its value type {value_type!r}")
[docs]
def get_score_in_region_agg(
self, chrom: str, start: int, end: int,
score: str | None = None,
aggregator: str | None = None,
none_value_replacement: ScoreValue | None = None,
) -> AggregateValue:
# pylint: disable=too-many-positional-arguments
"""Reduce a region to one value for one score.
This is the singular form of :meth:`get_scores_in_region_agg`.
Parameters
----------
chrom : str
The chromosome name.
start : int
The first 1-based position of the region.
end : int
The last 1-based position of the region.
score : str | None
The id of the score. ``None`` is accepted only when the
resource defines exactly one score.
aggregator : str | None
The name of the aggregator. ``None`` selects the aggregator
that the configuration of the score names, or the default one
for its value type.
none_value_replacement : ScoreValue | None
The value that stands for an uncovered position. ``None``
means that an uncovered position is skipped by the aggregator.
Returns
-------
AggregateValue
The aggregated value. It is ``None`` when the aggregator has
no value to report, for example when no position of the region
has a value. The aggregator sets the type of the value:
``list`` gives a list and ``value_count`` gives a dict.
Raises
------
ValueError
If the score is closed, the region is not valid, the score id
is unknown, no aggregator can be chosen, or
``none_value_replacement`` does not match the value type.
"""
return self.get_scores_in_region_agg(
chrom, start, end, [
PositionScoreAggregationQuery(
self._resolve_single_score(score),
aggregator, none_value_replacement),
])[0]
[docs]
def get_scores_in_region_agg(
self, chrom: str, start: int, end: int,
queries: Sequence[PositionScoreAggregationQuery],
) -> tuple[AggregateValue, ...]:
"""Reduce a region to one value for each query.
Each position of the region counts once. A record counts once for
each position of the region that it covers.
Parameters
----------
chrom : str
The chromosome name.
start : int
The first 1-based position of the region.
end : int
The last 1-based position of the region.
queries : Sequence[PositionScoreAggregationQuery]
The queries. Each names a score, an aggregator and a
replacement value for the uncovered positions.
Returns
-------
tuple[AggregateValue, ...]
One value for each query, in the order of ``queries``. A value
is ``None`` when the aggregator has no value to report. The
aggregator sets the type of its value: ``list`` gives a list
and ``value_count`` gives a dict. The tuple is empty when
``queries`` is empty.
Raises
------
ValueError
If the score is closed, the region is not valid, or a query is
not valid (see :meth:`resolve_aggregation_queries`).
Notes
-----
A chromosome that the resource never mentions is not an error. It
answers as a region that no record covers.
"""
# Defined on the per-position expansion -- an uncovered position is
# a ``None``, and with a ``none_value_replacement`` set it counts --
# but computed by walking segments, so cost stays proportional to
# record count. Where two records cover one position the first
# answers, so accumulated weight never exceeds the region width.
#
# A contig this score never mentions is uncovered rather than
# refused (gain#1211): it answers as a window on a contig the score
# has but does not cover. That is this read and the binned one
# only; a read that materialises positions still refuses an unknown
# contig.
self._guard_region_span(start, end)
targets, score_defs = self._resolve_aggregation_query_targets(queries)
for values, length in self._aggregating_runs(
chrom, start, end, score_defs):
for column, aggregator, none_value_replacement in targets:
value = values[column] if values is not None else None
if value is None:
value = none_value_replacement
aggregator.add(value, length)
return tuple(
aggregator.get_final() for _, aggregator, _ in targets)
def _resolve_aggregation_query_targets(
self, queries: Sequence[PositionScoreAggregationQuery],
) -> tuple[
list[tuple[int, Aggregator, ScoreValue]],
list[GenomicScoreDef]]:
"""Resolve queries to per-run fold targets, and the fetch columns.
One fetch serves every query: each DISTINCT score is fetched once,
and each query folds the column its score landed in -- so one score
may be requested twice with different aggregators, as the allele and
fragment aggregating reads allow too. Which scores those are, and in
what order, is :func:`~.aggregation.distinct_score_ids`, the derivation
the aggregating reads hand to
:func:`~.aggregation.fold_region_segments` as its column list: the
list both names what is fetched and indexes what comes back, so a
second spelling of it that ordered the scores differently would
have every aggregator quietly reading its neighbour's column.
The resolved DEFINITIONS are what comes back, and they are what the
read hands down: this is the aggregating pair's resolution of its
request, so re-deriving ids from them for a lower layer to resolve
again is the waste gain#1282 removed. The column index is taken
off the same list, in the same order, so the two still cannot
disagree about which column a query folds.
No contig reaches here: what the aggregating reads resolve is which
scores to fetch and how to fold them, which is a property of the
queries alone. Whether the contig exists is
:meth:`_aggregating_runs`' question, and it answers it by choosing a
run source rather than by refusing (gain#1211).
"""
resolved = self._resolve_aggregation_queries(queries)
score_defs = self._read_defs_for_any_contig(
distinct_score_ids(sid for sid, _, _ in resolved))
column_of = {
score_def.score_id: i for i, score_def in enumerate(score_defs)}
targets = [
(column_of[score_id], aggregator, none_value_replacement)
for score_id, aggregator, none_value_replacement in resolved
]
return targets, score_defs
[docs]
def get_score_in_bins(
self, chrom: str, start: int, end: int, bin_size: int,
score: str | None = None,
aggregator: str | None = None,
none_value_replacement: ScoreValue | None = None,
) -> Generator[tuple[int, int, ScoreValue], None, None]:
# pylint: disable=too-many-positional-arguments
"""Yield the aggregated value of one score for each bin of a region.
This is the singular form of :meth:`get_scores_in_bins`.
Parameters
----------
chrom : str
The chromosome name.
start : int
The first 1-based position of the region.
end : int
The last 1-based position of the region.
bin_size : int
The width of a bin in positions. It must be at least 1.
score : str | None
The id of the score. ``None`` is accepted only when the
resource defines exactly one score.
aggregator : str | None
The name of the aggregator. ``None`` selects the aggregator
that the configuration of the score names, or the default one
for its value type.
none_value_replacement : ScoreValue | None
The value that stands for an uncovered position. ``None``
means that an uncovered position is skipped by the aggregator.
Yields
------
tuple[int, int, ScoreValue]
The first position, the last position and the aggregated value
of a bin, in position order. The value is ``None`` when the
aggregator has no value to report for the bin.
Raises
------
ValueError
If the score is closed, the region is not valid, ``bin_size`` is
less than 1, the score id is unknown, no aggregator can be
chosen, or ``none_value_replacement`` does not match the value
type.
"""
bins = self.get_scores_in_bins(
chrom, start, end, bin_size, [
PositionScoreAggregationQuery(
self._resolve_single_score(score),
aggregator, none_value_replacement),
])
return (
(bin_start, bin_end, values[0])
for bin_start, bin_end, values in bins)
[docs]
def get_scores_in_bins(
self, chrom: str, start: int, end: int, bin_size: int,
queries: Sequence[PositionScoreAggregationQuery],
) -> Generator[tuple[int, int, tuple[ScoreValue, ...]], None, None]:
"""Yield the aggregated values of each bin of a region.
The bins follow the global grid that starts at position 1, so the
results of adjacent regions tile and are comparable across calls.
The first and the last bin are clipped to the region. The function
emits every bin of the region, also a bin that no record touches.
Parameters
----------
chrom : str
The chromosome name.
start : int
The first 1-based position of the region.
end : int
The last 1-based position of the region.
bin_size : int
The width of a bin in positions. It must be at least 1.
queries : Sequence[PositionScoreAggregationQuery]
The queries. Each names a score, an aggregator and a
replacement value for the uncovered positions.
Yields
------
tuple[int, int, tuple[ScoreValue, ...]]
The first position, the last position and one aggregated value
for each query of a bin, in position order. A value is
``None`` when the aggregator has no value to report.
Raises
------
ValueError
If the score is closed, the region is not valid, ``bin_size`` is
less than 1, or a query is not valid (see
:meth:`resolve_aggregation_queries`).
Notes
-----
A chromosome that the resource never mentions is not an error.
The function emits every bin of it, as for a region that no record
covers.
"""
# Edge bins are clipped to the query, so the yielded bounds name
# exactly what was aggregated. A segment straddling a bin boundary
# contributes its weight to each bin it touches, split at the
# boundary. The bin grid is the global one (``calc_bin_index`` /
# ``calc_bin_begin`` / ``calc_bin_end``).
#
# A contig this score never mentions is uncovered rather than
# refused (gain#1211): every bin of it is emitted, with the bounds
# a covered contig would yield. That is this read and
# get_scores_in_region_agg only; a read that materialises positions
# still refuses an unknown contig.
self._guard_region_span(start, end)
if bin_size < 1:
raise ValueError(
f"genomic score <{self.resource_id}> asked for bins of "
f"size {bin_size}; a bin holds at least one position")
targets, score_defs = self._resolve_aggregation_query_targets(queries)
return self._binned_runs(
self._aggregating_runs(chrom, start, end, score_defs),
start, end, bin_size, targets)
def _aggregating_runs(
self, chrom: str, start: int, end: int,
score_defs: list[GenomicScoreDef],
) -> Iterator[tuple[list[ScoreValue] | None, int]]:
"""The run source the two aggregating reads fold.
A contig the score holds records for is read from the table; one it
never mentions is read from an EMPTY segment stream (gain#1211). A
genome-wide fold over a track that skips a chromosome is the normal
case, not an error, and it needs no shape of its own: a region no
segment covers is already one ``(None, width)`` run out of
:meth:`_runs_from_segments`, so the absent contig answers as the
wholly-uncovered window it is, by the same code, with the width
computed once. Both branches are that one encoder over a different
stream -- which is why they cannot drift, and why a degenerate span
would fail identically in both rather than only in one.
The BRANCH lives here, above :meth:`_position_runs`, and not inside
it. Not because the per-position read would otherwise see it -- it
would not; that read refuses through :meth:`_region_read_defs`
before it ever builds these runs. The reason is the caller after
next: :meth:`_position_runs` is the plain region read, and a future
consumer composing it has no way to say "and I did mean to allow an
absent contig". Left down there the exemption would be the default
and silent; up here each read opts in by choosing this source, which
is the same reason :meth:`_read_defs_for_any_contig` is a second
door rather than a parameter on the first.
Choosing this source also means the backend is never asked for the
absent contig, so the refusals below -- the per-kind record
transform's, and each table's own -- stay exactly as they are and
are simply not reached. The contig is screened once per read here,
which is the question :meth:`_region_read_defs` used to ask on this
path; it is the same one screen, moved, not a new one.
That was true of the ABSENT contig from the start, and of the
covered one only since gain#1282: this screen chose the source, and
then the segment read under :meth:`_position_runs` resolved the
request over again and asked about the contig a second time. It
takes resolved definitions now, so the one screen claimed above is
the one that happens, on either branch.
Called a *scan* until gain#1304, which is what it was: a walk of the
ordered contig list, costing more the further down the list the
contig sat and most of all for one that was absent -- which on this
branch is the case the read exists to serve. It is a set lookup
now, so what is worth saying about the count is only that the read
asks once.
"""
if not self.has_chromosome(chrom):
return self._runs_from_segments(iter(()), start, end)
return self._position_runs(chrom, start, end, score_defs)
@staticmethod
def _binned_runs(
runs: Iterator[tuple[list[ScoreValue] | None, int]],
start: int, end: int, bin_size: int,
targets: list[tuple[int, Aggregator, ScoreValue]],
) -> Generator[tuple[int, int, tuple[ScoreValue, ...]], None, None]:
"""Fold position runs into grid bins, splitting at boundaries."""
bin_idx = calc_bin_index(bin_size, start)
pos = start
for values, length in runs:
remaining = length
while remaining > 0:
bin_stop = calc_bin_end(bin_size, bin_idx)
take = min(remaining, bin_stop - pos + 1)
for column, aggregator, none_value_replacement in targets:
value = values[column] if values is not None else None
if value is None:
value = none_value_replacement
aggregator.add(value, take)
pos += take
remaining -= take
if pos > bin_stop:
yield (
max(calc_bin_begin(bin_size, bin_idx), start),
bin_stop,
tuple(agg.get_final() for _, agg, _ in targets))
for _, aggregator, _ in targets:
aggregator.clear()
bin_idx += 1
if calc_bin_begin(bin_size, bin_idx) <= end:
yield (
max(calc_bin_begin(bin_size, bin_idx), start),
end,
tuple(agg.get_final() for _, agg, _ in targets))
[docs]
def get_scores_in_region(
self, chrom: str, start: int, end: int,
scores: Sequence[str] | None = None,
) -> Generator[tuple[ScoreValue | None, ...], None, None]:
"""Yield the values of the scores at each position of a region.
Parameters
----------
chrom : str
The chromosome name.
start : int
The first 1-based position of the region.
end : int
The last 1-based position of the region.
scores : Sequence[str] | None
The ids of the scores. ``None`` asks for every score that the
resource defines, in definition order.
Yields
------
tuple[ScoreValue | None, ...]
Exactly ``end - start + 1`` tuples, in position order. A tuple
has one value for each requested score. A ``None`` value means
that no record covers the position, or that the record has no
value for that score.
Raises
------
ValueError
If the score is closed, the region is not valid, the chromosome
is unknown, or a score id is unknown. The function raises
before the first tuple, not during the iteration.
"""
score_defs = self._plane_read_defs(chrom, start, end, scores)
return self._expand_position_runs(
self._position_runs(chrom, start, end, score_defs),
len(score_defs))