"""The :class:`GenomicScore` base class.
Everything the three score kinds share: config parsing and score-def
construction, the open/close lifecycle over the position table, and the
record and array reads. The kinds themselves live in :mod:`.position`,
:mod:`.allele` and :mod:`.fragment`, and the factories that dispatch
between them in :mod:`~gain.genomic_resources.genomic_scores.builders`.
Decomposing this class -- so that a kind's author reads the handful of hooks
their kind overrides rather than the whole base -- is gain#1027. Its first
extraction (gain#1044) moved the scoredef lifecycle to
:mod:`~gain.genomic_resources.score_def` and took this module under the
1500-line cap, so the file-scoped ``too-many-lines`` pragma gain#1007 added
here when it restored that cap is gone; its second (gain#1074) moved the
region-aggregation machinery to :mod:`.aggregation`, leaving
:meth:`GenomicScore.aggregate_region` here as the orchestrator that hands
it the per-kind weight rule; its third (gain#1114) moved the
value-extraction seam -- picking the per-record read, and addressing each
score def to a payload column -- to :mod:`.value_extraction`, leaving
:meth:`GenomicScore.open` calling both in the order that seam requires and
the two public per-record getters here. The remaining seams are #1027's
other children.
A score's defs are **finished in place at open**:
:func:`~.value_extraction.resolve_score_indices` writes ``score_index`` onto
the definitions this class already holds rather than handing back new ones
(it says there who reads them). The earlier half of the lifecycle differs:
``finish_scoredefs`` runs inside :meth:`GenomicScore._build_scoredefs`,
before there is a ``score_definitions`` to write onto, and so returns the
mapping.
"""
from __future__ import annotations
from abc import abstractmethod
from collections.abc import Generator, Iterator, Sequence
from types import TracebackType
from typing import (
TYPE_CHECKING,
Any,
ClassVar,
Self,
cast,
)
import numpy as np
from gain import logging
from gain.genomic_resources.bigwig_scores import (
build_bigwig_scoredefs,
validate_bigwig_scoredefs,
)
from gain.genomic_resources.genomic_position_table import (
BigWigTable,
ChromLengthSource,
VCFGenomicPositionTable,
build_genomic_position_table,
)
from gain.genomic_resources.genomic_position_table.record import (
ALT,
CHROM,
POS_BEGIN,
POS_END,
REF,
Record,
)
from gain.genomic_resources.repository import (
GenomicResource,
)
from gain.genomic_resources.resource_errors import inverted_span_error
from gain.genomic_resources.score_def import (
BULK_PARSEABLE_VALUE_TYPES,
GenomicScoreDef,
ScoreValue,
ValueExtractor,
build_genomic_score_schema,
finish_scoredefs,
parse_scoredef_config,
validate_scoredefs,
)
from gain.genomic_resources.score_filter import (
ScoreFilter,
compile_score_filter,
select_records,
)
from gain.genomic_resources.score_resource import (
ScoreResource,
refuse_unfoldable_histograms,
)
from gain.genomic_resources.vcf_scores import (
parse_vcf_scoredefs,
)
from .aggregation import (
build_region_aggregators,
fold_region_segments,
request_score_ids,
resolve_aggregator_requests,
)
from .records import (
RecordArrays,
)
from .value_extraction import (
resolve_score_indices,
select_value_extractor,
)
if TYPE_CHECKING:
# Only ever needed to type the VCF INFO proxies in annotations. pysam
# is a hard runtime dep and is already imported by the VCF table anyway, but
# keeping it behind TYPE_CHECKING makes it unambiguous that the annotations
# cost nothing at runtime.
pass
logger = logging.getLogger(__name__)
# Default rows-per-batch hint for GenomicScore.fetch_region_value_arrays. Big
# enough that the per-batch numpy overhead disappears against the per-row work
# it replaces, small enough that one batch's arrays stay comfortably in cache.
DEFAULT_VALUE_ARRAYS_BATCH_SIZE = 100_000
[docs]
class GenomicScore(ScoreResource[GenomicScoreDef]):
"""Base class for genomic score resources.
GenomicScore provides a unified interface for accessing and managing
genomic annotation scores stored in various formats. It serves as the
foundation for specialized score types including PositionScore (position-
based scores) and AlleleScore (variant-specific scores).
This abstract base class handles:
- Resource configuration validation and normalization
- Score definition management and parsing
- File format abstraction through GenomicPositionTable
- Histogram and statistics management
- Default annotation attribute configuration
- Context manager protocol for resource lifecycle
Score resources can be stored in multiple formats:
- Tabix-indexed files (TSV, BED)
- VCF files (particularly for allele scores)
- BigWig files (for position scores)
- In-memory tables (for testing)
Configuration Structure:
A genomic score resource requires a YAML configuration file
(genomic_resource.yaml) specifying:
- **type**: Resource type (position_score, allele_score)
- **table**: Table configuration with filename, format, and column
mappings for chrom, pos_begin, pos_end (and ref/alt for allele scores)
- **scores**: List of score definitions with id, type, name/index,
description, and optional aggregators
- **default_annotation**: Optional list specifying which scores to
include in default annotations with optional name mappings
- **histograms**: Optional histogram configurations for statistics
Score Definition:
Each score in the resource is defined with:
- **id**: Unique identifier for the score
- **type**: Data type (int, float, str, bool)
- **name/index**: Column name or index in the data file
- **desc**: Human-readable description
- **na_values**: Values to treat as missing/NA (optional)
- **hist_conf**: Histogram configuration for statistics (optional)
- **aggregator**: Default aggregator (optional). How several
values for one annotatable are reduced to one; the default
depends on the resource type and the score's value type.
Usage Pattern:
Genomic scores follow a resource lifecycle pattern:
1. Build/retrieve the resource from a repository
2. Create a score object from the resource
3. Open the score to initialize data access
4. Query scores using fetch methods
5. Close the score to release resources
Example using context manager:
>>> from gain.genomic_resources.genomic_scores import (
... build_score_from_resource_id
... )
>>> score = build_score_from_resource_id("phastCons100way")
>>> with score.open():
... # Score is open and ready to use
... chromosomes = score.get_all_chromosomes()
... scores = score.get_all_scores()
... # Query data...
>>> # Score is automatically closed
Statistics and Histograms:
GenomicScore supports automatic statistics generation including:
- Value distribution histograms
- Min/max ranges for numeric scores
- Category frequencies for categorical scores
- Custom histogram configurations per score
Attributes:
resource (GenomicResource): The underlying genomic resource object
resource_id (str): Unique identifier for the resource
config (dict): Validated and normalized configuration dictionary
table (GenomicPositionTable): Data access abstraction layer
score_definitions (dict[str, GenomicScoreDef]): Mapping of score IDs to
their internal definitions including parsers and metadata
table_loaded (bool): Flag indicating if the table is currently open
Key Methods:
open(): Initialize the score resource for data access
close(): Release resources and close the data table
get_all_scores(): Get list of all available score IDs
get_all_chromosomes(): Get list of all available chromosomes
has_chromosome(): Whether one chromosome is available
get_score_definition(): Get metadata for a specific score
get_default_annotation_attributes(): Get default annotation config
get_histogram(): Load histogram for a score (if available)
get_score_range(): Get value range for a numeric scores
Per-kind Methods:
A kind whose records read as something other than the span they
cover states that ONCE, by overriding:
- _score_segments(): what a region's raw records mean for this kind.
``region_values_from_records`` is the request resolution followed
by it, ``fetch_region_segments_scores`` is THAT applied to
``fetch_records``, and the statistics scan is it applied to
``validate_records(score, fetch_records(...))`` -- so a kind states
its reading once and every consumer gets it (ADR 0008). Override
this and not ``region_values_from_records``: the resolving entry is
shared by every kind, and a read holding an already-resolved
request composes this body without going through it (gain#1282).
- record_weight(): how many times one record's value counts when a
region is aggregated. Every reader goes through it -- the
annotators' ``aggregate_region``, the per-record scan, and the
bulk scan via ``record_weights`` (gain#1095).
- _aggregation_segments(): whether those records are cut down to
the queried window before they are weighed. Unlike the others
this HAS a default -- not clipping -- because it is a
consequence of the weight rule rather than a rule of its own: a
kind that counts a record once counts it wherever it falls.
All but the last have no default. A kind that inherited one would
be weighed by a rule nobody chose for it, which is the failure ADR
0008 exists to undo.
A kind's two validation rules are not on this list, and not on this
class: they are registered per kind in
:mod:`gain.genomic_resources.statistics.record_validation`, whose
default refuses a kind nobody wrote one for (ADR 0027).
See Also:
- PositionScore: For position-based genomic scores
- AlleleScore: For variant-specific genomic scores
- GenomicResource: Base resource abstraction
- GenomicPositionTable: Table format abstraction
"""
# How a value is read off a record. Installed by :meth:`open`, from the
# table's ``yields_records`` claim, and declared here with NO default on
# purpose: a record's payload means two different things -- a raw row or a
# VCF (variant, allele index) pair -- so no single extractor reads both,
# and a default would have to be wrong for one of them. Unset until open()
# routes, an unopened score raises AttributeError rather than silently
# reading a VCF record as a row; open() installs it *before* publishing
# table_loaded, so no caller can observe the gap (see open()).
_extract_value: ValueExtractor
# How a score of this resource type is reduced when a caller reads several
# values for one annotatable, keyed by the score's value type. Declared
# per CLASS because the reduction is a property of the resource type -- a
# position score is aggregated over a region of positions, an allele score
# over the alleles at one -- and a ``GenomicScoreDef`` cannot know which
# kind it belongs to.
#
# Handed to ``score_def.finish_scoredefs``, which applies it to every
# score whose config does not state an ``aggregator:`` -- today every
# deployed score: 0 of 16502 resource configs set one.
DEFAULT_AGGREGATORS: ClassVar[dict[str, str | None]] = {}
def __init__(self, resource: GenomicResource):
self.resource = resource
self.resource_id = resource.resource_id
assert self.resource.config is not None
self.config: dict = self.resource.config
self.config = self.validate_and_normalize_schema(
self.config, resource,
)
self.config["id"] = resource.resource_id
self.table_loaded = False
self.table = build_genomic_position_table(
self.resource, self.config["table"],
)
self.score_definitions = self._build_scoredefs()
[docs]
@staticmethod
def get_schema() -> dict[str, Any]:
"""The config this kind accepts; each kind extends it."""
return build_genomic_score_schema()
def _build_scoredefs(self) -> dict[str, GenomicScoreDef]:
"""Route this resource's definitions through their construction path.
The one piece of the scoredef lifecycle that did NOT move to
:mod:`~gain.genomic_resources.score_def` in gain#1044: it dispatches
on the table's type into ``parse_vcf_scoredefs`` and
``build_bigwig_scoredefs``, and both of those modules import
``score_def``, so hosting this there would close an import cycle.
Everything it calls is a function now, and the class's only
contribution is ``DEFAULT_AGGREGATORS``, passed explicitly -- once,
at the single exit the three routes converge on, which is where
``finish_scoredefs`` documents that it has to be applied.
"""
config_scoredefs = None
if "scores" in self.config:
config_scoredefs = parse_scoredef_config(self.config)
scoredefs: dict[str, GenomicScoreDef]
if isinstance(self.table, VCFGenomicPositionTable):
merge = bool(self.config.get("merge_vcf_scores", False))
scoredefs = parse_vcf_scoredefs(
cast(dict[str, Any], self.table.header),
config_scoredefs,
resource_id=self.resource_id,
merge=merge)
elif config_scoredefs is None:
raise ValueError("No scores configured and not using a VCF")
elif isinstance(self.table, BigWigTable):
scoredefs = build_bigwig_scoredefs(self.config, config_scoredefs)
else:
scoredefs = config_scoredefs
# Refused AFTER finish_scoredefs, where an unstated ``type:``
# becomes ``float`` -- earlier would judge a score by a type it
# does not end up with -- and at the convergence of all three
# construction routes, so a ``scores:`` block, a VCF header and a
# bigWig are held to the one rule (gain#1336).
return refuse_unfoldable_histograms(
finish_scoredefs(scoredefs, self.DEFAULT_AGGREGATORS),
self.resource_id)
[docs]
def get_config(self) -> dict[str, Any]:
"""The configuration, validated and normalized at construction."""
return self.config
[docs]
def get_default_annotation_attributes(self) -> list[Any]:
"""Collect default annotation attributes."""
default_annotation = self.get_config().get("default_annotation")
if default_annotation is None:
return [
{"source": attr, "name": attr}
for attr in self.score_definitions
]
if not isinstance(default_annotation, list):
raise TypeError(
"The default_annotation in the "
f"{self.resource_id} resource is not a list.")
return default_annotation
[docs]
def get_default_annotation_attribute(self, score_id: str) -> str | None:
"""Return default annotation attribute for a score.
Returns None if the score is not included in the default annotation.
Returns the name of the attribute if present or the score if not.
"""
attributes = self.get_default_annotation_attributes()
result = []
for attr in attributes:
if attr["source"] != score_id:
continue
dst = score_id
if "name" in attr:
dst = attr["name"]
result.append(dst)
if result:
return ",".join(result)
return None
[docs]
def close(self) -> None:
"""Close the underlying table and mark the score not open.
:meth:`open` may be called again afterwards.
"""
self.table.close()
self.table_loaded = False
[docs]
def is_open(self) -> bool:
"""Whether :meth:`open` has run and :meth:`close` has not since."""
return self.table_loaded
[docs]
def open(self) -> Self:
"""Open genomic score resource and returns it.
**Validate and route BEFORE opening, and so before publishing.** Every
input to both steps is known at construction -- the table's class, its
``yields_records`` ClassVar, and the score definitions -- so neither
needs the open handle, and two things fall out of that order:
* a refusal costs no handle. Routing after ``table.open()`` would
leave a caller that is not using the ``with`` form holding an opened
pysam handle it cannot reach: ``table_loaded`` would still be
False, so ``close()`` would not have been reached. Raising first
means there is nothing to leak. The bigWig config validation sits
here for exactly that reason.
* ``table_loaded = True`` is what makes this score look open to
everyone else: from that write on, another caller's open() takes the
is_open() early return above and reads ``_extract_value`` straight
away. Routed last, that caller could catch the score
published-but-unrouted, and since the routing has no default at
all, that caller reads an AttributeError. Scores are
shared across threads (the process-wide in-memory fragment-score
cache; gain-web-api's thread pool), so the window is reachable;
this ordering keeps the ROUTING out of it. Pinned by
test_the_score_is_routed_before_it_reports_itself_open.
It does not make open() as a whole safe to race, and does not claim to:
``resolve_score_indices`` still runs after the score has published
itself open, so a caller that catches that window reads a score def
with no ``score_index`` yet. That window is older than this ordering
and untouched by it -- open() is not synchronised, and making it so is
a separate change.
"""
if self.is_open():
logger.info(
"opening already opened genomic score: %s",
self.resource.resource_id)
return self
is_vcf = isinstance(self.table, VCFGenomicPositionTable)
is_bigwig = isinstance(self.table, BigWigTable)
if is_bigwig:
validate_bigwig_scoredefs(
self.resource_id, self.score_definitions)
self._extract_value = select_value_extractor(
score_definitions=self.score_definitions,
table=self.table,
is_vcf=is_vcf,
is_bigwig=is_bigwig)
self.table.open()
self.table_loaded = True
# A bigWig's score config has already been validated, and by a stricter
# rule: ``validate_bigwig_scoredefs`` permits no column addressing at
# all (bar the deprecated ``index: 3``), where this method *demands*
# one whenever the table reports a header. A bigWig table has no
# header to speak of -- but one whose config carries a stray
# ``header:``/``header_mode:`` pair reports one anyway, and those keys
# are ignored for bigWig (see ``genomic_position_table.utils``), so
# they must not decide how the scores are checked either.
#
# A VCF's is skipped too: its scores have no column address to
# check, and what an entry may say about one was already refused
# where the definitions were built (``_refuse_overridden_address``
# says why). The tabular route always has a ``scores:`` block --
# ``_build_scoredefs`` refuses a non-VCF resource without one.
if not (is_bigwig or is_vcf):
validate_scoredefs(self.config, self.table, self.resource)
resolve_score_indices(
self.score_definitions,
is_vcf=is_vcf,
is_bigwig=is_bigwig,
table=self.table,
resource_id=self.resource_id)
return self
def __enter__(self) -> Self:
return self
def __exit__(
self,
exc_type: type[BaseException] | None,
exc_value: BaseException | None,
exc_tb: TracebackType | None,
) -> None:
if exc_type is not None:
logger.error(
"exception while working with genomic score: %s, %s, %s",
exc_type, exc_value, exc_tb)
self.close()
def _get_header(self) -> tuple[Any, ...] | None:
assert self.table is not None
return self.table.header
[docs]
def compile_filter(self, expression: str) -> ScoreFilter:
"""Compile a boolean expression into a filter over this score.
The expression names this resource's own scores and relates them
with ``>``, ``>=``, ``<``, ``<=``, ``==``, ``!=`` and ``in``,
combined with ``not``, ``and`` and ``or``; the result is
passed back to any of the record reads as ``score_filter``.
Raises :class:`~gain.genomic_resources.score_filter.ScoreFilterError`
on an expression that does not parse or that names a score this
resource does not define. See
:func:`~gain.genomic_resources.score_filter.compile_score_filter` for
what compiling settles,
``docs/adr/0017-score-filtering-is-a-score-capability.md`` for why
the capability sits on the score, and
``docs/adr/0018-score-filter-grammar-extension.md`` for the
language's precedence and what a name may contain.
"""
return compile_score_filter(self, expression)
[docs]
def fetch_records(
self,
chrom: str,
pos_begin: int | None,
pos_end: int | None,
*,
score_filter: ScoreFilter | None = None,
) -> Generator[Record, None, None]:
"""Yield the records of a region, optionally filtered.
A caller reads a record's positional fields from its slots
(``record[CHROM]``, ``record[POS_BEGIN]``, ...) and a score value
through :meth:`get_score_value_from_record` or
:meth:`get_score_values_from_record` on this score.
``score_filter`` is a predicate from :meth:`compile_filter`, applied
to each record; only records it accepts are yielded. ``None`` --
the default -- yields what the table yields, and is the whole of
what this method did before filtering became a score capability.
``chrom`` is required, here and throughout the region-read family.
A caller that wants every record of a table asks the table:
``score.table.get_all_records()``.
Nothing here is checked before the first ``next()``, the filter's
ownership included: every backend's
``get_records_in_region`` is itself a generator function, so an
unknown contig has always been reported from the first record read
rather than from the call, and there is no eagerness left to
preserve by structuring this any other way. That is a property of
*this* read rather than a rule for the family: a read that
materialises has no generator body to defer a refusal into, and
:meth:`AlleleScore.fetch_allele_records()
<.allele.AlleleScore.fetch_allele_records>` accordingly refuses from
the call.
"""
records = self.table.get_records_in_region(chrom, pos_begin, pos_end)
yield from select_records(self, records, score_filter)
[docs]
def get_score_value_from_record(
self, record: Record, score_id: str,
) -> ScoreValue:
"""Read one configured score off a record of this score's table."""
return self._extract_value(record, self.score_definitions[score_id])
def _resolve_score_defs(
self, scores: Sequence[str] | None,
) -> list[GenomicScoreDef]:
"""Resolve requested score ids to definitions, refusing unknown ones.
``None`` asks for every score this resource defines. A score id the
resource does not define is a caller error, and it is refused here --
before any data is read -- so the refusal does not depend on whether
the queried region happens to hold a record. A typo answering
differently on a populated contig than on an empty one is the failure
this exists to prevent.
"""
if scores is None:
scores = self.get_all_scores()
unknown = [
score_id for score_id in scores
if score_id not in self.score_definitions
]
if unknown:
raise ValueError(
f"genomic score <{self.resource_id}> does not define "
f"{sorted(unknown)}; it has "
f"{sorted(self.score_definitions)}")
return [self.score_definitions[score_id] for score_id in scores]
def _resolve_single_score(self, score: str | None) -> str:
"""Resolve a singular method's ``score`` argument to one score id.
``None`` means "all the scores this resource has", which a singular
method can honour only when there is exactly one.
Here rather than on the one kind that calls it today: the rule is
about what a SINGULAR read can resolve, which is a property of the
score rather than of the positions it is read at, and the fragment
and allele planes reach for it as they grow singular reads of their
own. Until they do, ``PositionScore`` is its only caller.
"""
if score is not None:
return score
all_scores = self.get_all_scores()
if len(all_scores) != 1:
raise ValueError(
f"genomic score <{self.resource_id}> defines "
f"{sorted(all_scores)}; a singular read can resolve "
f"score=None only when there is exactly one")
return all_scores[0]
[docs]
def get_score_values_from_record(
self, record: Record, score_defs: list[GenomicScoreDef],
) -> list[ScoreValue]:
"""Read several scores off one record, for ALREADY-resolved defs.
The bulk counterpart of :meth:`get_score_value_from_record`: a caller
resolves score names to definitions once per fetch and passes them per
record, so the name->definition lookup stays out of the per-record
loop.
"""
extract = self._extract_value
return [extract(record, score_def) for score_def in score_defs]
[docs]
def supports_region_value_arrays(self, scores: list[str]) -> bool:
"""Whether :meth:`fetch_region_value_arrays` will serve these scores.
Answers the two things a caller can be wrong about: the backend
serves the bulk column-array read, AND every named score is one this
facade can parse. A predicate that answered only the first would say
True for a call that then refuses -- not a capability query but a trap.
It is not a promise the call cannot fail for some OTHER reason. A
score whose configured column index does not exist in its backend's
payload still raises (deliberately -- see ``BigWigTable``), and so
does a closed score or an unknown contig. This answers "is this score
the kind this method serves", not "is every argument valid".
The value-type half is not a consumer's condition leaking in: the
facade parses, so it serves the value types
:meth:`GenomicScoreDef.parse_array()
<gain.genomic_resources.score_def.GenomicScoreDef.parse_array>`
defines a column parse for
(:data:`~gain.genomic_resources.score_def.BULK_PARSEABLE_VALUE_TYPES`)
and no others. What a *consumer* additionally needs stays with the
consumer: the statistics scan also requires a bounded region, and
each of its two entry points asks for its own histogram or min/max
pairing (see
``genomic_scores_impl.scan.can_bulk_histogram`` and
``can_bulk_min_max``). The scan does not re-test the resource
KIND; ADR 0001 records why. What it does
NOT require is a particular record shape: the accumulator reads the
kind's own ``record_weight`` and the scan's door reads the rule
registered for the kind in
:mod:`gain.genomic_resources.statistics.record_validation`, so a
position, allele and fragment score are all served.
Answerable on an UNOPENED score: the table and the score definitions
are both built in ``__init__``, so nothing here touches the file.
"""
if not self.table.supports_value_arrays:
return False
for score_id in scores:
score_def = self.score_definitions.get(score_id)
if score_def is None \
or score_def.value_type not in BULK_PARSEABLE_VALUE_TYPES:
return False
return True
def _require_open(self) -> None:
"""Refuse any read of a score that is not open.
The one home of this message. Reads that also need a known contig
take :meth:`_require_open_and_known_chrom`; the two aggregating
reads of :class:`~.position.PositionScore`, which accept a contig
the score never mentions, take this one alone (gain#1211).
"""
if not self.is_open():
raise ValueError(f"genomic score <{self.resource_id}> is not open")
def _require_open_and_known_chrom(self, chrom: str) -> None:
"""Refuse a region read this score cannot answer at all.
The two conditions every bulk column read shares, stated once for
the readers that widen it. Several OLDER reads in this module spell
the contig half out inline; they are left as they are rather than
swept into this change, and a few of them word it differently on
purpose (an allele read names the resource in it).
"""
self._require_open()
if not self.has_chromosome(chrom):
raise ValueError(
f"{chrom} is not among the available chromosomes.")
def _guard_region_span(self, start: int, end: int) -> None:
"""Refuse a span no genomic region can mean.
A 1-based lower bound and an end that does not precede its start.
There is deliberately no UPPER bound check: the exact chromosome
length is not knowable for most scores, and a position past the end
of the data is simply uncovered (see #727).
Here rather than on the one kind that spelled it first: the rule is
about what a REGION can be, which is a property of neither the
positions a score is read at nor the records it holds. Two kinds
call it -- the position and fragment logical planes, whose region
bounds are mandatory ints. ``AlleleScore`` does not, and that is
not an omission: it has no logical plane, and its record reads take
``int | None`` bounds where ``None`` legitimately means unbounded.
It also stops a backend's own leniency from reaching a reader --
the in-memory table tests its bounds for truthiness, so a ``0``
would otherwise arrive as "unbounded" and answer a caller error
with a whole contig.
Called by each read rather than folded into
:meth:`_region_read_defs`, which is the seam the shared
:meth:`region_values_from_records` already runs for every kind and
which already receives the two positions it ignores. Folding it in
is the deeper placement and is deliberately not taken here: it would
refuse ``fetch_*`` requests that are accepted today, on all three
kinds at once, which is a behaviour change no reader of this slice
asked for. Until that is decided, a read that takes a mandatory
region calls this first.
"""
if start < 1:
raise ValueError(
f"genomic score <{self.resource_id}> asked for a region "
f"with start {start}; positions are 1-based")
if end < start:
raise ValueError(
f"genomic score <{self.resource_id}> asked for a region "
f"whose end {end} precedes its start {start}")
def _value_arrays_refusal_reason(self) -> str:
"""Why :meth:`supports_region_value_arrays` said no, for a raiser.
The two halves of that predicate, worded for a caller who ignored it.
Stated here rather than at each raise site so the reason cannot drift
from the predicate it explains, nor between the readers that widen it
(:meth:`AlleleScore.fetch_region_allele_arrays`).
"""
if not self.table.supports_value_arrays:
return (
f"its {type(self.table).__name__} backend leaves "
f"supports_value_arrays False")
return (
"not every requested score has a value type the column "
f"parse serves {sorted(BULK_PARSEABLE_VALUE_TYPES)}")
[docs]
def fetch_region_value_arrays(
self,
chrom: str,
pos_begin: int | None,
pos_end: int | None,
scores: list[str],
*,
batch_size: int = DEFAULT_VALUE_ARRAYS_BATCH_SIZE,
) -> Generator[
RecordArrays, None, None]:
"""Fetch a region as column arrays, without building a record per row.
The bulk counterpart of :meth:`fetch_records`, for a caller that scans a
whole region and wants columns rather than rows -- statistics, above
all. Each batch is ``(pos_begin, pos_end, {score_id: values})``: the
one-based position arrays, plus one array of **parsed** values per
requested score.
**Values are parsed, by the same contract the per-record read uses.**
Each column goes through :meth:`GenomicScoreDef.parse_array()
<gain.genomic_resources.score_def.GenomicScoreDef.parse_array>`, whose
agreement with the per-value :meth:`GenomicScoreDef.parse_value()
<gain.genomic_resources.score_def.GenomicScoreDef.parse_value>` is
pinned by test_parse_array_agrees_with_parse_value_fuzz. So NA
sentinels and unparseable cells arrive as that score's non-value,
whatever the backend stores underneath: a ``float`` or ``int`` score
yields ``float64`` with ``nan`` for no value, a ``str`` score an
``object`` array with ``None``. **The array's dtype follows the
score's declared type, not the backend's** -- a caller reading several
scores in one batch can be handed both shapes.
That parse is why a value type the definition cannot parse as a column
is refused, and why :meth:`supports_region_value_arrays` asks about
the scores and not only about the backend.
**It does NOT clip.** A record overlapping the region's start is
yielded whole, exactly as :meth:`fetch_records` yields it; trimming to
``[pos_begin, pos_end]`` is the caller's, because what a partial
overlap means depends on what the caller is computing.
``batch_size`` is a HINT. A backend whose read granularity is fixed by
its own windowing -- ``BigWigTable``, whose batches are sized by its
adaptive fetch window -- ignores it.
Each score id gets an array of its own -- the parse builds one per
id, so two ids sharing a payload column do not alias.
The guards below run when this method is CALLED, not on the first
``next()`` -- which is why the streaming half lives in
``_value_array_batches`` rather than a ``yield`` here.
"""
if not self.supports_region_value_arrays(scores):
# Refuse here rather than let the call reach the table. A VCF
# table INHERITS the tabix implementation, so an unguarded call
# does not fail cleanly -- it trips that method's
# ``assert isinstance(self.pysam_file, pysam.TabixFile)`` and
# yields a message-less AssertionError (nothing at all under
# ``python -O``). Probing this capability by catching is therefore
# not viable; ask supports_region_value_arrays() first.
reason = self._value_arrays_refusal_reason()
raise TypeError(
f"genomic score <{self.resource_id}> does not serve "
f"fetch_region_value_arrays for {sorted(scores)}: {reason}. "
f"Ask supports_region_value_arrays(scores) before calling.")
self._require_open_and_known_chrom(chrom)
return self._value_array_batches(
self._score_column_indexes(scores),
(chrom, pos_begin, pos_end), batch_size)
def _value_array_batches(
self,
columns: dict[str, int],
region: tuple[str, int | None, int | None],
batch_size: int,
) -> Generator[
RecordArrays, None, None]:
"""Stream the batches for an already-validated request.
Split out so :meth:`fetch_region_value_arrays` is a plain function and
its guards fire when it is CALLED. Were it a generator itself, every
one of those checks would be deferred to the first ``next()``, so a
caller that built the generator and passed it elsewhere would be handed
the refusal at some arbitrary later point, far from the mistake.
"""
for begin, end, values, _cells in self._parsed_column_batches(
columns, region, batch_size):
yield begin, end, values
def _parsed_column_batches(
self,
columns: dict[str, int],
region: tuple[str, int | None, int | None],
batch_size: int,
extra_columns: frozenset[int] = frozenset(),
) -> Generator[
tuple[np.ndarray, np.ndarray,
dict[str, np.ndarray], dict[int, np.ndarray]],
None, None]:
"""The column read every bulk reader is made of: parse, plus cells.
One statement of the parse loop, because there is more than one
reader over it: :meth:`_value_array_batches` and
:meth:`AlleleScore._allele_array_batches`. Two copies of it is how
the two would come to disagree about a batch's positions, its NA
handling or its dtypes -- the drift ADR 0008 spends its length on.
``extra_columns`` are fetched but NOT parsed, and reach the caller
through the raw ``cells`` alongside the parsed values. That is the
whole of what a reader wanting a non-score column adds: it asks for
the index and reads it out itself, rather than teaching this loop
what the column means.
"""
chrom, pos_begin, pos_end = region
defs = {
score_id: self.score_definitions[score_id]
for score_id in columns
}
wanted = set(columns.values()) | set(extra_columns)
for begin, end, cells in self.table.get_region_value_arrays(
chrom, pos_begin, pos_end, wanted, batch_size):
yield begin, end, {
score_id: defs[score_id].parse_array(cells[column])
for score_id, column in columns.items()
}, cells
def _score_column_indexes(self, scores: list[str]) -> dict[str, int]:
"""Score id -> payload column index, resolved once for a whole scan.
No cast needed: ``score_index`` is an ``int``. A VCF score is
addressed by ``col_name`` and has none, which is how the type already
says the VCF backend does not reach here.
"""
return {
score_id: self.score_definitions[score_id].score_index
for score_id in scores
}
[docs]
def get_all_chromosomes(self) -> list[str]:
"""The chromosome names the score's table holds, in table order.
Raises ``ValueError`` on a score that is not open.
"""
self._require_open()
return self.table.get_chromosomes()
[docs]
def has_chromosome(self, chrom: str) -> bool:
"""Answer whether this score's table carries ``chrom``.
The yes/no half of :meth:`get_all_chromosomes`, and what every
contig screen in the read path asks -- the annotators' pre-read
screens, the shared region-read refusal, the position kind's
absent-contig branch and both allele reads' refusals. All of them
spelled ``chrom not in self.get_all_chromosomes()``, a walk of the
ordered list whose cost grew with the resource's contig count and
with the contig's index in it, and grew WORST for a contig the
resource does not carry, since that walks the whole list before it
can say no (gain#1304).
Raises ``ValueError`` on a score that is not open, exactly as
:meth:`get_all_chromosomes` does and for the same reason: those
screens got that refusal for free from the accessor they used, and a
predicate that answered ``False`` instead would turn "this score was
never opened" into "this score does not carry that contig" at every
one of them.
The ordered list keeps its meaning, its order and its aliasing for
the callers that genuinely need the collection -- a statistics scan
splitting a genome into region tasks, the resource implementation's
contig report. This is for the ones that only ever asked a question.
"""
self._require_open()
return self.table.has_chromosome(chrom)
@property
def chrom_length_source(self) -> ChromLengthSource:
"""What the score's own file can say a contig's length is.
Declared on the backend class, so answered on a closed score --
which is what a caller weighing whether to open the table for its
lengths at all (coverage's second rung, gain#1448) needs. The
lengths themselves come from
:func:`~.chrom_lengths.derive_chrom_lengths`.
"""
return self.table.chrom_length_source
[docs]
def resource_files(self) -> set[str]:
"""The resource's files the score reads: the data file, plus the
index on a backend that reads one.
Answered by the table, which knows how it opens, on a closed score
-- the file set has to be known before any of the files is.
"""
return self.table.resource_files()
[docs]
def region_values_from_records(
self,
records: Iterator[Record],
chrom: str,
pos_begin: int | None = None, # ruff: ignore[unused-method-argument]
pos_end: int | None = None, # ruff: ignore[unused-method-argument]
scores: Sequence[str] | None = None,
) -> Generator[
tuple[int, int, list[ScoreValue]], None, None]:
"""Extract this kind's ``(begin, end, values)`` from raw records.
The region read expressed as a function OF a record stream, which is
what lets the two consumers of a region differ by what they COMPOSE
rather than by a flag: :meth:`fetch_region_segments_scores` is this
applied to :meth:`fetch_records`, and the statistics scan is this
applied to
``validate_records(score, fetch_records(...))``. Neither can quietly
acquire the other's behaviour, and no argument travels down to say
which of the two is reading (ADR 0008).
``chrom``, ``pos_begin`` and ``pos_end`` name the region the records
were asked for. Nothing is fetched here, and nothing is reshaped to
the window either -- what a partial overlap means belongs to the
caller (ADR 0008); a consumer answering a question about the window
clips with :func:`~.records.clip_span`. The positions are what the
guards below are about.
The guards run when this is CALLED rather than on the first
``next()`` -- the pattern :meth:`fetch_records` documents -- which is
why the streaming half lives in ``_score_segments``. They
stay here rather than moving down into ``fetch_records`` because
that is not where the request is: half of what they resolve is the
score ids, which ``fetch_records`` is not even given, and the contig
half would be a THIRD screen of the same contig in one read -- after
this one and before the backend's own, which refuses an unknown
contig from ``get_records_in_region`` regardless. This used to be
argued on cost, from a time when each screen walked the ordered
contig list; since gain#1304 the screen is a set lookup and the
argument is only about where a request is resolved, which is here.
What a kind yields is :meth:`_score_segments`, and not this method:
the resolution above is the same for every kind, the reading below
is not, and since gain#1282 they are split so that a kind states
only the half that is its own. Override THAT to say what a record
means here; overriding this one would take the resolution with it,
and would be skipped by a read that enters below it holding
definitions it has already resolved.
"""
score_defs = self._region_read_defs(chrom, scores)
return self._score_segments(records, score_defs)
def _region_read_defs(
self, chrom: str, scores: Sequence[str] | None,
) -> list[GenomicScoreDef]:
"""Refuse a region request this score cannot serve, before any record.
Shared by every kind's :meth:`region_values_from_records`, so a
closed score, an unknown contig and an unknown score id are refused
alike whatever the kind. The score ids are resolved once for the
whole region rather than per record, and before the first record
rather than on it: a typo answering differently on a populated contig
than on an empty one is the failure that eagerness prevents.
That argument is about a read that MATERIALISES a contig's positions,
and it is why this refusal stays the default. The two aggregating
reads of :class:`~.position.PositionScore` are the deliberate
exception: they answer a question about a WINDOW rather than about
the contig, and a genome-wide fold over a track that skips a
chromosome is the normal case, not a typo (gain#1211). They take
:meth:`_read_defs_for_any_contig` instead and compose the absent
contig as one uncovered run; every other read, this method's own
callers included, keeps the refusal.
"""
self._require_open_and_known_chrom(chrom)
return self._resolve_score_defs(scores)
def _read_defs_for_any_contig(
self, scores: Sequence[str] | None,
) -> list[GenomicScoreDef]:
"""Resolve a read's score ids, without judging any contig.
:meth:`_region_read_defs` minus the contig refusal, for the reads
that treat a contig the score never mentions as uncovered rather
than as an error (gain#1211). Everything else it checks is checked
here, in the same order: a closed score and an unknown score id are
refused on an absent contig exactly as on a populated one, so the
exemption widens what a contig may be and nothing else.
Offered by the base because the guards it reuses are, but it is
:class:`~.position.PositionScore`'s two aggregating reads that may
take it and no others: "an absent contig is uncovered" is a claim on
the logical read plane, and a kind without one -- an allele score --
has nothing for it to mean. A kind that takes this door instead of
:meth:`_region_read_defs` answers "no records" to a typo.
"""
self._require_open()
return self._resolve_score_defs(scores)
def _score_segments(
self,
records: Iterator[Record],
score_defs: list[GenomicScoreDef],
) -> Generator[
tuple[int, int, list[ScoreValue]], None, None]:
"""Stream this kind's segments for an already-resolved request.
**The per-kind hook.** :meth:`region_values_from_records` minus
the resolution -- the half of a region read that differs by kind,
split from the half that does not. A kind states its reading by
overriding THIS, and inherits one resolution of the request rather
than spelling out a second (gain#1282); an allele score is the one
kind that does, reading each record as the point it sits at.
It is also why the two ways into the segment stream cannot drift in
what a segment is: :meth:`region_values_from_records` resolves and
composes this body, and a caller that has resolved already composes
the same body directly. A kind that overrode the resolving entry
instead would be read by one and skipped by the other, which is
drift by another name.
Every record is yielded at its full extent, including one that only
partly overlaps -- or entirely misses -- the region it was fetched
for. What a partial overlap means depends on what the caller is
computing, so deciding it belongs to the window-answering consumers,
each of which clips with :func:`clip_span` (ADR 0008). That is the
BASE body's reading -- what a position score and a fragment score
both mean by a record -- and an overriding kind replaces it.
A record whose end precedes its begin is refused: that is a claim
about the record itself, not about any window. A record outside
the queried region is NOT refused -- a backend answering a region
query with such a record is misconfigured rather than holding bad
data (a table whose index and ``pos_end`` name different columns,
gain#553), which ADR 0008 refuses at ``open()`` and deliberately
not here.
The hottest loop in the read path, so it reads its record slots
directly rather than through the helpers that wrap them (gain#823):
:meth:`get_score_values_from_record` is a method call around a
comprehension over defs already resolved for the whole region, and it
remains, unchanged, for its other callers. The refusal those helpers
carried is kept here, as one comparison raising the shared
:func:`~gain.genomic_resources.resource_errors.inverted_span_error`;
see ``test_segment_path_refuses_a_backwards_record``.
"""
extract = self._extract_value
for record in records:
rec_begin = record[POS_BEGIN]
rec_end = record[POS_END]
if rec_end < rec_begin:
raise inverted_span_error(
record[CHROM], rec_begin, rec_end,
record[REF], record[ALT])
yield (rec_begin, rec_end, [
extract(record, score_def) for score_def in score_defs])
[docs]
def fetch_region_segments_scores(
self,
chrom: str,
pos_begin: int | None = None,
pos_end: int | None = None,
scores: Sequence[str] | None = None,
*,
score_filter: ScoreFilter | None = None,
) -> Generator[
tuple[int, int, list[ScoreValue]], None, None]:
"""Yield ``(begin, end, values)`` per record touching the region.
One tuple per underlying RECORD -- a segment, at that record's own
extent -- not one value per position. The region's records, read as
this kind means them. A record straddling the region's edge is
reported whole: what a partial overlap means depends on what the
caller is computing, so a caller answering a question about the
window composes :func:`~.records.clip_to_region` over this stream, or
calls :func:`~.records.clip_span` per segment (ADR 0008).
``score_filter`` -- from :meth:`compile_filter` -- travels to
:meth:`fetch_records` and nowhere else: the records it rejects
never reach the transform, so a rejected record costs no value
extraction. It reads the RECORD, so it may name any score the
resource defines, including one outside ``scores``. It is not a
second kind of read and it selects nothing the record read would
not: the composition below is the same one either way, with a
filtered stream in place of an unfiltered one.
Its ownership check is the one refusal here that does NOT land on
the call. It rides :meth:`fetch_records`, whose generator body
defers it and which says so, where the request checks
:meth:`region_values_from_records` runs are eager -- so a filter
compiled against a different score is refused on the first
``next()``, not from the call that a closed score, an unknown
contig and an unknown score id are refused from.
A plain read: it checks nothing. The statistics scan reads the same
records through the same transform with
:func:`~gain.genomic_resources.statistics.record_validation.validate_records`
composed in front, and that extra link -- visible at the consumer,
in ``genomic_scores_impl/scan.py`` -- is the whole of the
difference
between the two (ADR 0008).
One body per kind, in :meth:`_score_segments`, rather than one per
kind per consumer: two that had to agree is how the paths drift.
The fragment plane reads through this method for that reason
(gain#1272), and the position plane reaches the same body with an
already-resolved request (gain#1282) -- different entries, one
reading.
"""
return self.region_values_from_records(
self.fetch_records(
chrom, pos_begin, pos_end, score_filter=score_filter),
chrom, pos_begin, pos_end, scores)
[docs]
@classmethod
@abstractmethod
def record_weight(cls, left: int, right: int) -> int:
"""How many times one record's value counts when aggregating.
The rule is a property of the resource TYPE: a position-score
record counts once per base pair of the queried region it covers,
an allele line counts once, a fragment counts once however long it
is. One record, one count, is the answer for everything except a
position score.
**The kind's single statement of that rule**, and every reader that
weighs a RECORD goes through it: :meth:`aggregate_region` folds
with it, the per-record statistics scan calls it, and the bulk scan
broadcasts it over a whole batch through :meth:`record_weights`.
One statement, so a kind cannot weigh its records one way in one
of those and another way in the next. Pinned by
test_the_weight_rule_is_stated_once_per_kind.
A position score's logical plane is deliberately NOT among them.
It tiles POSITIONS rather than folding records --
``get_scores_in_region_agg`` weighs each run by its length -- and
since gain#1131 that is the path annotation takes. The two agree
for this kind, a clipped record's width being exactly the length of
the run it covers, but they agree by arithmetic rather than by
reading one statement, so a change here is not automatically a
change there.
**An implementation must be numpy-elementwise** -- an arithmetic
expression over ``left`` and ``right``, or a constant. It is
declared over scalars because that is what its per-record callers
hand it, but :meth:`record_weights` answers a whole batch by
handing it the position COLUMNS instead, and only an elementwise
body gives the same answers that way.
Most ways of breaking that break loudly -- a body branching on
``left`` raises ``ValueError`` on an array's ambiguous truth, one
calling ``int()`` a ``TypeError``. The contract is written down
for the ones that do not: a body REDUCING its arguments
(``int(np.mean(right - left + 1))``) hands back a plain number,
which is then broadcast as though it were every record's weight,
and the two scan paths disagree with nothing raised. That is what
test_the_weight_rule_is_stated_once_per_kind's broadcast-agreement
assertion is there to catch.
Deriving a weight from the span unconditionally is what this hook
exists to prevent: it would give a fragment its length as a weight
and disagree with the fragment score annotator for every fragment
longer than one base pair.
"""
raise NotImplementedError
[docs]
@classmethod
def record_weights(
cls, begins: np.ndarray, ends: np.ndarray,
) -> np.ndarray:
""":meth:`record_weight` over a whole batch's position columns.
The bulk statistics scan has no record to hand the scalar hook, so
it weighs a batch here instead. This does not restate the rule --
it broadcasts the ONE statement of it, which is why the scan may
not read the weight anywhere else.
The widening is the elementwise contract being spent: the hook is
declared over scalars and its bodies are arithmetic, so the same
expression answers a column. A kind whose weight is a CONSTANT
answers with that constant however it was called, so a 0-d result
is filled out to the batch's shape rather than treated as an
error.
"""
# One marker for one fact -- the scalar-declared hook is being
# called on whole columns, which is the contract its docstring
# states and this method exists to spend.
weights = np.asarray(
cls.record_weight(begins, ends)) # type: ignore[arg-type]
if weights.ndim == 0:
weights = np.full(begins.shape, weights)
return weights.astype(np.int64, copy=False)
[docs]
def aggregate_region(
self,
chrom: str,
pos_begin: int | None = None,
pos_end: int | None = None,
scores: list[str | tuple[str, str]] | None = None,
) -> list[ScoreValue]:
"""Reduce a region to one value per requested score.
The aggregating counterpart of :meth:`fetch_region_segments_scores`,
which it is built on: that method yields one entry per record, this
one folds
those entries into a single value per request.
Each request is either a score id -- aggregated with the resource's
own default, which ``score_def.finish_scoredefs`` resolved from
this class's ``DEFAULT_AGGREGATORS`` -- or a ``(score_id, aggregator)``
pair naming one explicitly. The aggregator string is whatever the
config accepts, parametrized forms (``join(,)``) included.
**Returns a list parallel to ``scores``, not a dict.** One score
may legitimately be requested twice with different aggregators --
``["s", ("s", "max")]``, which is what an annotation config does
when it exposes one source as both a min and a max attribute -- and
a dict keyed by score id would silently drop one of them.
An empty region is not an error: each aggregator answers for
itself. ``list`` returns ``[]``; ``max`` returns ``None``; and so
does ``count``, which chooses to report nothing rather than 0 for an
empty region (see ``CountAggregator.get_final``). This method does
not second-guess any of them. That is deliberately unlike the
per-position reads (``get_scores_at_position``,
``fetch_allele_scores``), which answer ``None`` where there is no
data -- aggregating nothing is a well-defined question, reading a
value where there is none is not.
Values reach the aggregator exactly as the record carried them,
``None`` included, because that is what the annotators do (each
aggregator decides what a null means for it) and the point of this
method is to give the answer they would.
"""
requests = resolve_aggregator_requests(
scores,
score_definitions=self.score_definitions,
all_scores=self.get_all_scores(),
resource_id=self.resource_id,
)
# Built BEFORE the fetch, because the fetch is not lazy: the
# not-open and unknown-contig guards of fetch_region_segments_scores run
# when it is CALLED, not on the first next(). An aggregator built
# afterwards would have a misspelled name reported only for the
# regions a resource happens to cover, so `mediann` would be a
# missing-contig error until someone queried a covered contig.
aggregators = build_region_aggregators(
requests, resource_id=self.resource_id)
score_ids = request_score_ids(requests)
return fold_region_segments(
self._aggregation_segments(
chrom, pos_begin, pos_end, score_ids),
aggregators,
requests,
score_ids=score_ids,
weigh=self.record_weight,
)
def _aggregation_segments(
self,
chrom: str,
pos_begin: int | None = None,
pos_end: int | None = None,
scores: list[str] | None = None,
) -> Iterator[tuple[int, int, list[ScoreValue]]]:
"""The segment stream this kind's aggregating reads consume.
The records as this kind means them, which for everything but a
position score is :meth:`fetch_region_segments_scores` unchanged: a kind
that counts a record ONCE counts it wherever the point it collapses
to falls, window or not (see
test_an_allele_point_outside_the_window_still_aggregates_once).
Clipping the records to the window first is a position-score fact,
so it is stated on :class:`~.position.PositionScore` and nowhere
else. This hook is what lets it be: without it the fold would have
to carry a flag saying which kind it is serving, and that flag
would be a second statement of the weight rule.
Underscored by a criterion rather than by a list of callers, as
:meth:`_score_segments` is: the hooks that are part of the read API
are asked for by name (``fetch_region_segments_scores`` IS
:meth:`region_values_from_records`; the scan calls
:meth:`record_weight` by name, and reads the kind's validation rule
by dispatching on its class), while
this one is never a caller's question -- it is composed, from
inside this hierarchy, by whichever reads aggregate. A kind
overrides it; nothing outside the hierarchy calls it. Stated as
the criterion because the census of those reads is not stable: it
was one until gain#1087 gave the annotators' read the same stream.
"""
return self.fetch_region_segments_scores(
chrom, pos_begin, pos_end, scores)