"""The ``allele_score_annotator``.
Annotates with scores keyed by allele -- frequencies, pathogenicity
predictions and the like -- read from an ``allele_score`` resource, by
exact allele match or reduced over a region.
"""
import textwrap
from typing import Any
from gain.annotation.annotatable import Annotatable, VCFAllele
from gain.annotation.annotation_config import (
AnnotationConfigurationError,
AnnotatorInfo,
Attribute,
)
from gain.annotation.annotation_pipeline import (
AnnotationPipeline,
Annotator,
AttributeSpec,
)
from gain.annotation.annotator_base import AnnotatedValues
from gain.annotation.genomic_score_annotator_base import (
GenomicScoreAnnotatorBase,
)
from gain.genomic_resources.aggregators import (
ScoreAggregationQuery,
)
from gain.genomic_resources.genomic_scores import (
AlleleAggregate,
AlleleScore,
allele_key,
build_allele_score_from_resource,
)
from gain.genomic_resources.resource_types import (
PREFERRED_ALLELE_SCORE_TYPE,
)
from gain.genomic_resources.score_filter import ScoreFilterError
[docs]
def build_allele_score_annotator(pipeline: AnnotationPipeline,
info: AnnotatorInfo) -> Annotator:
return AlleleScoreAnnotator(pipeline, info)
[docs]
class AlleleScoreAnnotator(GenomicScoreAnnotatorBase):
"""Annotator for allele-level genomic scores (frequencies, pathogenicity…).
Operates in one of two modes, selected by the ``mode`` parameter:
- ``allele`` (**default**): a ``VCFAllele`` the resource keys as a line
of its own gets an exact chrom/pos/ref/alt lookup, exact on the
position too. On a resource declaring ``allele_multiplicity: one``
-- the default -- that is the allele's only line and its scores
(``AlleleScore.get_allele_scores_for_allele``); on a ``many``
resource the allele's lines are reduced, each attribute by its
aggregator, as a region fold reduces them
(``AlleleScore.get_allele_scores_for_allele_agg``). Which alleles
are matched is the resource's
``allele_score_mode``: every ``VCFAllele`` on an ``alleles``
resource, only a substitution on a ``substitutions`` one. Any other
``VCFAllele`` -- an indel or complex allele on a ``substitutions``
resource -- takes the ``region`` path over the bases it covers,
subject to ``region_length_cutoff``, as does every annotatable that
is not a ``VCFAllele``.
- ``region``: the score reduces all allele lines that overlap the
annotatable's span, in one streaming walk
(``AlleleScore.get_allele_scores_in_region_agg``). Works with any
``Annotatable`` (``VCFAllele``, ``Region``, CNV, …), and folds every
``VCFAllele`` whatever the resource's mode, subject to
``region_length_cutoff``. An aggregator
must be defined for every score attribute, either in the attribute
config or as the score's ``aggregator`` default in the resource YAML;
an attribute with neither -- only a ``bool`` score can be in that
position -- is refused when the pipeline loads, in either mode,
because a CNV or a region takes the region path whatever the mode.
Virtual ``allele`` attribute
----------------------------
All annotators expose a virtual attribute ``"allele"``
(``is_default=False``)
that is synthesised rather than read from the data file.
- On an exact match: returns ``["chrom:pos:ref:alt"]`` for the matched
line. On a ``many`` resource, the distinct keys of the allele's
lines that pass ``allele_filter`` -- one, or several when an
``include_attributes`` score differs across them.
- On a region fold, in either mode: returns the distinct
``"chrom:pos:ref:alt"`` strings of the lines that pass the optional
``allele_filter``, in the order the lines were first met -- the
resource's own genomic order.
Optionally append score values to each allele string with
``include_attributes``. The string's format is the score's,
:func:`~gain.genomic_resources.genomic_scores.allele.allele_key`, so
the two modes cannot drift.
An aggregator named on this attribute reduces nothing, in either
mode. It is not a score: its value is the keys the annotator
synthesised, and ``region`` mode has always answered them beside the
reductions rather than as one of them. Exact-match mode used to
differ -- the base folded its one-element list -- and stopped in
gain#1133, so the two modes now agree.
``allele_filter``
-----------------
An optional annotator-level boolean expression evaluated against each
record before it is included in the result. The annotator only resolves
the parameter; the expression language belongs to the score, so see
:meth:`GenomicScore.compile_filter` for the operators it admits and what
a name may contain.
"""
ACCEPTED_RESOURCE_TYPES = (PREFERRED_ALLELE_SCORE_TYPE,)
def __init__(self, pipeline: AnnotationPipeline, info: AnnotatorInfo):
resource = self.resolve_resource(pipeline, info)
self.allele_score = build_allele_score_from_resource(resource)
self.allele_filter = None
allele_filter_str = info.parameters.get("allele_filter")
if allele_filter_str is not None:
assert isinstance(allele_filter_str, str)
try:
self.allele_filter = self.allele_score.compile_filter(
allele_filter_str)
except ScoreFilterError as e:
# Named after the parameter the user wrote: the score knows
# nothing about how the expression reached it (cf. gain#477).
raise AnnotationConfigurationError(
f"Error parsing allele_filter: {e}") from e
mode = info.parameters.get("mode", "allele")
if mode not in {"allele", "region"}:
raise AnnotationConfigurationError(
f"Invalid mode '{mode}' for allele_score_annotator; "
"valid values are 'allele' and 'region'")
self.mode = mode
super().__init__(pipeline, info, self.allele_score)
info.documentation += textwrap.dedent(f"""
Annotator to use with scores that depend on allele like
variant frequencies, etc.
**Mode** (``mode`` parameter, applies to ``VCFAllele`` inputs only):
- ``allele`` (default): routed by the resource's ``allele_score_mode``.
On an ``alleles`` resource every allele is an exact chrom/pos/ref/alt
match. On a ``substitutions`` resource only a substitution is; an
insertion, deletion or complex allele aggregates the allele lines over
the bases it covers.
- ``region``: aggregates scores for all allele lines overlapping the
annotatable's span, whatever the resource's mode.
Non-``VCFAllele`` annotatables always use region aggregation.
<a href="{self.BASE_DOC_URL}#allele-score-annotator" target="_blank">More info</a>
""") # ruff: ignore[line-too-long]
self.allele_attribute = None
self.attrs_to_include: list[str] = []
for attr in self._attributes:
if attr.source == "allele":
attrs_to_include = attr.parameters.get(
"include_attributes", [])
if isinstance(attrs_to_include, str):
attrs_to_include = [attrs_to_include]
self.attrs_to_include = list(attrs_to_include)
self.allele_attribute = attr
continue
self.add_score_aggregator_documentation(
attr, "aggregator", attr.aggregator)
# One query per SCORE attribute, as `PositionScoreAnnotator` builds
# its list and for the reasons its query block gives. The virtual
# `allele` attribute is not a score and asks the read for the keys
# instead.
self._region_queries = [
ScoreAggregationQuery(attr.source, self._query_aggregator(attr))
for attr in self._attributes
if attr is not self.allele_attribute
]
# Resolved now for the refusal alone -- at load, in BOTH modes; see
# the `region` bullet of the class docstring -- and the same for the
# `include_attributes` ids.
self.allele_score.resolve_aggregation_queries(self._region_queries)
self.allele_score.resolve_allele_key_scores(self.attrs_to_include)
# The scores a `one` resource's exact match reads: the score
# attributes' and the ones the allele keys are suffixed with.
self._exact_scores: list[str] = list(dict.fromkeys(
[*self.simple_score_queries, *self.attrs_to_include],
)) or self.allele_score.get_all_scores()
[docs]
def get_attribute_defaults(
self, spec: AttributeSpec,
) -> dict[str, Any]:
defaults = super().get_attribute_defaults(spec)
if "aggregator" not in defaults:
score_def = self.allele_score.get_score_definition(spec.source)
if score_def is not None \
and score_def.aggregator is not None:
defaults["aggregator"] = score_def.aggregator
return defaults
[docs]
def get_attribute_specs(self) -> dict[str, AttributeSpec]:
"""Return score attribute specs plus the virtual ``allele``."""
result = super().get_attribute_specs()
result["allele"] = AttributeSpec(
source="allele",
value_type="list",
description="The allele in the format 'chr:pos:ref:alt'",
is_default=False,
internal_default=False,
)
return result
[docs]
def build_score_aggregator_documentation(
self, attr: Attribute,
) -> list[str]:
"""Collect score aggregator documentation."""
allele_doc = self._build_score_aggregator_documentation(
attr, "aggregator", attr.aggregator,
)
return [allele_doc]
def _annotate_allele(
self, annotatable: VCFAllele,
) -> AnnotatedValues:
"""Answer an exact chrom/pos/ref/alt match off the allele plane.
Read by the resource's ``allele_multiplicity``: the allele's only
row on a ``one`` resource, its rows folded on a ``many`` one.
Either way the match is exact on the position too, so a row
starting before ``pos`` whose span reaches it is not the allele.
"""
if self.allele_score.multiplicity is AlleleScore.Multiplicity.MANY:
return self._annotate_allele_rows(annotatable)
return self._annotate_allele_row(annotatable)
def _annotate_allele_row(
self, annotatable: VCFAllele,
) -> AnnotatedValues:
"""Answer the only row of one allele on a ``one`` resource.
The row's values as they stand, nothing folded; ``None`` under
every attribute when no row holds the allele or the filter rejects
it.
"""
values = self.allele_score.get_allele_scores_for_allele(
annotatable.chrom,
annotatable.position,
annotatable.reference,
annotatable.alternative,
scores=self._exact_scores,
score_filter=self.allele_filter,
)
if values is None:
return self._empty_result()
# Widened, because the virtual `allele` attribute below is a LIST of
# strings and the score's own values are scalars.
scores: dict[str, Any] = dict(
zip(self._exact_scores, values, strict=True))
if self.allele_attribute is not None:
# The same helper the region read builds its keys with, so
# the two paths spell an allele identically.
scores[self.allele_attribute.source] = [allele_key(
annotatable.chromosome, annotatable.position,
annotatable.reference, annotatable.alternative,
[scores[a] for a in self.attrs_to_include])]
# Not ``fold_own_values``: the virtual ``allele`` attribute's
# key list is the answer, not something to reduce.
return self._from_sources(scores)
def _annotate_allele_rows(
self, annotatable: VCFAllele,
) -> AnnotatedValues:
"""Fold the rows of one allele on a ``many`` resource.
Each score attribute by its own aggregator, resolved as on the
region fold; the virtual ``allele`` attribute takes the distinct
keys of the rows. ``None`` under every attribute when no row
holds the allele; rows the filter rejects in full are an empty
selection, which each aggregator answers for.
"""
aggregate = self.allele_score.get_allele_scores_for_allele_agg(
annotatable.chrom, annotatable.position,
annotatable.reference, annotatable.alternative,
queries=self._region_queries,
allele_keys=self._allele_keys_request(),
score_filter=self.allele_filter,
)
return self._from_aggregate(aggregate)
def _allele_keys_request(self) -> list[str] | None:
"""What a folding read is asked for keys: ``None`` builds none."""
return (
self.attrs_to_include
if self.allele_attribute is not None else None)
def _from_aggregate(
self, aggregate: AlleleAggregate | None,
) -> AnnotatedValues:
"""Answer a folding read's aggregate, keyed by attribute name.
``None`` is absent data -- no record was there -- and answers
``None`` for every attribute. An aggregate whose fold saw nothing is
different: records were there and the filter rejected them all,
so each aggregator has answered for an empty selection and the
keys are empty. The values are paired back over the attributes
that built the queries, the ``allele`` attribute taking the keys.
"""
if aggregate is None:
return self._empty_result()
return self._pair_aggregated(
aggregate.values, len(self._region_queries),
resource_id=self.allele_score.resource_id,
reduced=lambda attr: attr is not self.allele_attribute,
otherwise=lambda _attr: list(aggregate.allele_keys or ()))
def _annotate_region(
self, annotatable: Annotatable,
) -> AnnotatedValues:
"""Answer the region already reduced, keyed by attribute name.
The SCORE reduces (gain#1163): one value per query and, when the
virtual ``allele`` attribute is configured, the distinct allele
keys, off a single walk that never materialises the records.
That removes the per-record list this path used to hold beside
the aggregators (gain#834); what an aggregator itself keeps --
``list``, the ``str`` default, keeps every value -- is the
aggregator's property and stays.
"""
aggregate = self.allele_score.get_allele_scores_in_region_agg(
annotatable.chrom, annotatable.position, annotatable.pos_end,
queries=self._region_queries,
allele_keys=self._allele_keys_request(),
score_filter=self.allele_filter,
)
return self._from_aggregate(aggregate)
def _matched_exactly(self, allele: VCFAllele) -> bool:
"""Tell whether the resource keys ``allele`` as a line of its own.
Every allele on an ``alleles`` resource; only a substitution on a
``substitutions`` one. The score's reads are mode-blind, so this
is the one place the resource's mode is honoured.
"""
if self.allele_score.alleles_mode():
return True
return allele.type == Annotatable.Type.SUBSTITUTION
def _do_annotate(
self, annotatable: Annotatable,
context: dict[str, Any], # ruff: ignore[unused-method-argument]
) -> AnnotatedValues:
"""Route the annotatable to the exact match or the region fold.
A ``VCFAllele`` in ``allele`` mode goes where the resource's
``allele_score_mode`` says it can be matched (gain#1748): every
allele on an ``alleles`` resource, only a substitution on a
``substitutions`` one. The rest -- and a ``VCFAllele`` in
``region`` mode, an explicit override -- is region-folded.
Every fold is subject to ``region_length_cutoff``; the exact
match is not a fold and never is.
"""
if not self.allele_score.has_chromosome(annotatable.chromosome):
return self._empty_result()
if (
isinstance(annotatable, VCFAllele)
and self.mode != "region"
and self._matched_exactly(annotatable)
):
return self._annotate_allele(annotatable)
too_long = self._beyond_region_length_cutoff(annotatable)
if too_long is not None:
return too_long
return self._annotate_region(annotatable)