"""Reading a VCF's INFO fields as genomic scores.
Everything the score layer knows about VCF, in one module: how a VCF header's
INFO metadata becomes score definitions, and how one of those scores is read
off a record. Both encode the same thing -- INFO field semantics, and
``Number=1``/``A``/``R``/``.`` in particular -- and they used to sit 362 lines
apart in ``genomic_scores``.
**The VCF table itself is not here and does not belong here.**
``genomic_position_table.table_vcf`` produces records: it owns the payload's
shape, the pysam proxies it carries and the constants that name them. This
module interprets those records as scores. That is the same seam the record
contract draws everywhere else -- a backend yields records, the score layer
says what they mean -- and it is why the table layer still imports nothing
from the score layer.
"""
from __future__ import annotations
from typing import Any
from gain import logging
from gain.genomic_resources.genomic_position_table.record import (
CHROM,
PAYLOAD,
POS_BEGIN,
Record,
)
from gain.genomic_resources.genomic_position_table.table_vcf import (
ALLELE_INDEX,
INFO,
INFO_META,
VARIANT,
)
from gain.genomic_resources.resource_errors import score_configuration_error
from gain.genomic_resources.score_def import GenomicScoreDef, ScoreValue
logger = logging.getLogger(__name__)
VCF_TYPE_CONVERSION_MAP = {
"Integer": "int",
"Float": "float",
"String": "str",
"Flag": "bool",
}
#: The INFO ``Number`` values a header DECLARES for a field whose value
#: reaches a score's parser as a scalar. pysam decodes ``0`` (a ``Flag``) to
#: a ``bool`` and ``1`` to a single value, and :func:`extract_vcf_value`
#: indexes ``A``/``R`` down to one element before parsing. Every other
#: declared shape -- unbounded ``.`` and any fixed arity above one --
#: decodes to a tuple, and the only thing that reads a tuple is the
#: ``|``-joining ``converter``. (The genotype-arity ``G`` is on this side
#: too, but never reaches the converter: pysam will not read a per-genotype
#: INFO field, so :func:`_refuse_genotype_arity` refuses the definition.)
#:
#: It decides THREE things, on both the header side and the config-override
#: side: which fields get that converter, which may take a stated ``type:``
#: as their parser, and which DECLARE the ``str`` the join produces rather
#: than the ``Type=`` of one element (gain#1259). All three are the same
#: question -- does a value reach the parser whole -- so they are answered
#: from one set and cannot drift apart.
#:
#: It is the DECLARED number, not the shape that actually arrives: a row
#: may carry two values for a field its header calls ``Number=1``. This
#: set says what the header promises; :func:`extract_vcf_value` is what
#: holds a row to it (gain#1257).
#:
#: ``0`` is in the set on the strength of the CONFIG side. A ``Flag`` needs
#: no join -- ``converter`` returns a ``bool`` unchanged -- so the header side
#: reads the same whether it is called scalar or not. A config that types one
#: does not: dbSNP's ``GNO`` is ``Number=0,Type=Flag`` in the file and
#: ``type: int`` in the resource, and it is the config's ``int`` that makes it
#: read ``1``/``0`` rather than ``True``/``False``, with a categorical
#: histogram built on those. (``RV``, typed ``bool``, cannot show the
#: difference.)
_SCALAR_VALUED_NUMBERS = (0, 1, "A", "R")
def _check_number_arity(
record: Record, score_def: GenomicScoreDef, number: str | int,
count: int,
) -> None:
"""Warn if an INFO field's value count is not what its ``Number`` fixes.
``Number=A`` declares one value per ALT allele and ``Number=R`` one per
allele including the reference, so on a well-formed record the count
follows from the ALT column; ``Number=1`` fixes it at one. (``0``, a
``Flag``, is accepted for symmetry with the scalar set; pysam never
hands a tuple for one.) No mismatch is rejected anywhere -- htslib
does not enforce ``Number`` on read, and resource load never looks at a
row -- and every shape is silent to a reader: a SHORT per-allele tuple
leaves the alleles past its end with no value (they read null, see
:func:`extract_vcf_value`), an OVER-LONG one has values no allele will
ever select, and a tuple where a scalar was promised has no rule that
picks one of its elements, so the row reads null (gain#1257). Only the
resource's author can fix any of these, and they cannot fix what
nothing reports.
**Once per field per table, not once per line.** The flag is
``score_def.number_mismatch_warned``, and the definition is per-table
(see its comment in ``score_def``, which also says why one flag serves
both shapes). This is the per-record score read -- keeping it cheap is
what #237 was about -- and a field whose arity is wrong on one row is
normally wrong on every row, so per-line logging would bury the run in
identical lines. Testing the flag FIRST is what bounds the cost on a
malformed table: once it has warned, every later record of it stops at
a single attribute read.
A WELL-FORMED per-allele table has no such short circuit and pays the
ALT lookup on every record -- there is no cheaper way to learn how many
values a row ought to carry, and remembering the first row's count would
only be right for a table whose alleles never vary. The cost lands
where it can be afforded: only inside the tuple branch of the read,
which a well-formed ``Number=1`` row never enters (it decodes to a
scalar) and which already builds a fresh pysam ``VariantMetadata`` per
read -- so it is a fraction added to an already-expensive branch, and
nothing at all on the common scalar shape.
The message names the field, its number, the count it carried against
the count it owed, and the row that tripped it. It cannot name the
resource: the read is a pure function of ``(record, score_def)`` and
neither carries a resource id. The locus is the better half of that
trade anyway -- it names the offending row of the offending file, which
is what an author has to open.
"""
if score_def.number_mismatch_warned:
return
if number in ("A", "R"):
alts = record[PAYLOAD][VARIANT].alts
# ``alts`` is None for a record whose ALT is absent ('.'): zero ALT
# alleles, so a Number=A field carries no values at all and a
# Number=R field carries only the reference's.
expected = len(alts) if alts is not None else 0
if number == "R":
expected += 1
owed = (
f"for {expected} allele(s); the VCF row is malformed and an "
f"allele with no value of its own reads as null")
else:
expected = 1
owed = (
"where the header declares one; the VCF row is malformed and "
"reads as null")
if count == expected:
return
score_def.number_mismatch_warned = True
logger.warning(
"INFO field %s (Number=%s) of %s:%s carries %s value(s) %s. "
"Reported once per table.",
score_def.score_id, number, record[CHROM], record[POS_BEGIN],
count, owed)
def _reports_empty_elements(meta: Any) -> bool:
"""Whether an absent element of this field can be called malformed.
Only in a ``String`` field. There pysam decodes the spec's own
missing-value token ``.`` to the literal string ``'.'`` and reserves
``None`` for an element the row declared and did not supply, so a
``None`` is malformed data and can be reported as such.
Every other type collapses the two: ``AF=0.5,.`` and ``AF=0.5,`` both
decode to ``(0.5, None)``, and so does a whole-field ``AF=.``. A report
there would call a spec-legal row malformed and send the resource's
author to a locus with nothing wrong at it, on data where a missing
element is ordinary -- so the empty element is dropped, silently.
"""
return bool(meta.type == "String")
def _report_empty_element(
record: Record, score_def: GenomicScoreDef,
) -> None:
"""Report this field's empty element, once per field per table.
The flag is ``score_def.empty_element_warned`` -- the field's own, not a
share of the arity check's, so an arity report and an empty-element
report cannot silence each other. Why once, and why the message names
the field and the locus, is #289's reasoning unchanged: see
:func:`_check_number_arity`.
"""
if score_def.empty_element_warned:
return
score_def.empty_element_warned = True
logger.warning(
"INFO field %s of %s:%s carries an empty element; the VCF row is "
"malformed and an element with no value of its own contributes "
"nothing to the value read. Reported once per table.",
score_def.score_id, record[CHROM], record[POS_BEGIN])
def _report_no_value(
record: Record, score_def: GenomicScoreDef,
) -> None:
"""Report a field present with no value at all, once per table.
``ORIGIN=`` names the key and supplies nothing after the ``=``, which
pysam decodes to the empty tuple -- no element to be empty, so
:func:`_report_empty_element` never sees it, and the score reads null.
Shares that function's flag: both say one thing about one field -- it
declared a value it does not carry -- and a reader needs to hear it once.
"""
if score_def.empty_element_warned:
return
score_def.empty_element_warned = True
logger.warning(
"INFO field %s of %s:%s is present but carries no value at all; the "
"VCF row is malformed and the score reads as null. Reported once "
"per table.",
score_def.score_id, record[CHROM], record[POS_BEGIN])
def _drop_empty_elements(
record: Record, score_def: GenomicScoreDef, value: tuple, meta: Any,
) -> tuple:
"""Return ``value`` without its empty elements, reporting the first row.
An element of a multi-valued INFO field is EMPTY when the row declares it
and supplies nothing -- ``ORIGIN=1,``, ``ORIGIN=a,,b``, ``ORIGIN=,b`` --
and pysam decodes such an element as ``None``. Nothing rejects it:
not pysam, not resource load. The read answers it by the rule #256 and
#289 settled for the per-allele shapes -- no value, no score -- so the
element contributes nothing to the joined value. (A tuple left with no
elements at all is turned into a null by the caller, the only place that
knows what the join would otherwise have produced.)
``None`` is the whole test, deliberately: this is NOT a falsy filter. A
``String`` field's ``.`` decodes as the literal string ``'.'`` and is a
value dbSNP's ``CAF``/``TOPMED`` carry meaningfully, ``'0'`` is a value,
and ``''`` is a value; only an element pysam decoded as absent is one.
Dropping is unconditional; REPORTING is not, and
:func:`_reports_empty_elements` says which fields have earned it.
The membership test comes FIRST, before both the type test and the flag
-- the reverse of the arity check's order. It is the only work a
WELL-FORMED table pays here (one C-level scan of a tuple that has to be
walked to be joined anyway) and it returns the tuple itself, so the
common shape allocates nothing; testing the flag first would put an
attribute read in front of every well-formed record to save one on the
malformed ones.
"""
if None not in value:
return value
if _reports_empty_elements(meta):
_report_empty_element(record, score_def)
return tuple(element for element in value if element is not None)
# How a score's value is read off a record: chosen once per opened score by
# ``GenomicScore.open``, from the table's type, and called per value.
#
# This replaces the four score-line CLASSES the score layer used to route
# between (#239 had already reduced them to records plus a wrapper; this
# removes the wrapper). A score line existed to hold two things -- which
# payload slot a score lives in, and, for VCF, the per-record pysam proxies --
# and neither needs an object any more: the first is ``score_index`` on the
# definition, the second is in the payload. What is left is a value read that
# is a pure function of ``(record, score_def)``, so the per-line allocation
# goes and the routing stays exactly where it was.
def _refuse_overridden_type(
resource_id: str, score_id: str, number: Any, meta_type: str,
config_type: str | None, *, is_scalar: bool,
) -> None:
"""Refuse a stated ``type:`` the field's join cannot produce.
**The whole rule lives here, and the call is unconditional.** Whether
an entry contradicts its header is one question -- it stated a type, the
field is not scalar, and the type is not already the ``str`` the join
produces -- and splitting it across the call site would leave a reader of
this docstring believing it owns a rule it only half owns.
A ``scores:`` entry over a multi-valued INFO field may describe it
freely -- ``desc``, aggregators and a histogram config are all its own --
but it may not claim a value the field cannot hold. What the field
reads is ``|``-joined text, so its value type is ``str`` by definition;
an entry stating anything else is a contradiction between the config and
the header, and gain#1336 raises on it.
**Stating ``str`` is not a contradiction and is not refused.** It is
the type the join produces, so nothing is being claimed that the field
cannot hold, and it is what every multi-valued entry in the deployed
GRRs states (ClinVar's twenty, dbSNP's ``CAF``/``TOPMED``). Leaving
``type:`` unstated is equally fine -- the definition takes ``str`` from
the header side either way.
The refusal is raised while the definitions are BUILT, so it reaches
every consumer of the resource: ``repo-repair`` reports the resource
failed by name before any region task is planned, and an annotation
pipeline naming it fails to load with the same message.
``resource_id`` is threaded in from
:meth:`GenomicScore._build_scoredefs` for the refusal messages alone --
a parse handed only a header and a config could name the field but not
the resource, which is not enough to find the file to edit in a
repository of thousands.
"""
if config_type is None or is_scalar or config_type == "str":
return
# No "state 'str' instead" advice: that edit is not inert. An entry's
# NA sentinels are normalized against the type it states, so rewriting
# the type also changes which values read null (gain#1284). The message
# says what was claimed and what the field holds; the author chooses
# between stating ``str`` and dropping the ``type:`` line.
raise score_configuration_error(
resource_id, score_id,
f"states 'type: {config_type}', but its ##INFO line declares "
f"Number={number},Type={meta_type}: a field the header declares "
f"multi-valued reads '|'-joined text, so its value type is 'str'. "
f"State 'type: str' or leave 'type:' unstated.")
def _refuse_genotype_arity(
resource_id: str, score_id: str, header_entry: Any,
) -> None:
"""Refuse a field the header declares ``Number=G``.
``G`` is a legal INFO arity -- one value per genotype -- and the header
parses; what pysam will not do is READ the value: ``info.get`` on a row
carrying such a field raises ``ValueError: genotype is only valid as a
format field``. :func:`extract_vcf_value` makes that lookup outside the
``try`` in :meth:`GenomicScoreDef.parse_value` that guards the PARSE, on
purpose, so the error escaped a fetch uncaught, naming neither resource
nor field, from a resource that had opened without complaint
(gain#1258). The shape is visible in the header, so it is refused here
instead, where the resource and the field can both be named and every
consumer sees one attributed error -- through
:func:`score_configuration_error`, which says why a header-only
resource still counts as one.
"""
if header_entry.number != "G":
return
raise score_configuration_error(
resource_id, score_id,
f"is declared Number=G,Type={header_entry.type} in its ##INFO line; "
f"pysam does not read a per-genotype INFO field, so this score can "
f"never be read. List the fields you want in a 'scores:' block that "
f"leaves it out ('merge_vcf_scores' unset or false), or change the "
f"header.")
def _refuse_undeclared_id(
resource_id: str, score_id: str, vcf_header_info: dict[str, Any],
) -> None:
"""Refuse a ``scores:`` entry whose ``id`` no ``##INFO`` line declares.
The entry's ``id`` is the one thing that has to match the header
(:func:`parse_vcf_scoredefs` says why), so one that does not is a
contradiction between the config and the header, like the two rules
above, and is refused the same way: while the definitions are built,
naming the resource and the score.
It used to escape as a bare ``KeyError`` from whichever indexer reached
it first -- pysam's header metadata in filter mode, the header-derived
definitions in merge mode -- and ``KeyError`` is not a fault
``report_resource_failure`` attributes to a resource, so ``repo-repair``
printed an unexpected internal error with a traceback (gain#1489).
The message lists every id the header DOES declare, in header order,
so a typo reads off the line. Not capped: the deployed VCF resources
declare forty-odd fields (ClinVar 39, dbSNP 48), so any cap short
enough to matter would hide the list exactly where it is needed.
"""
if score_id in vcf_header_info:
return
raise score_configuration_error(
resource_id, score_id,
f"names no ##INFO field of the VCF header, which declares: "
f"{', '.join(vcf_header_info)}. A VCF score's 'id' is its INFO "
f"key; fix the 'scores:' entry or the header.")
def _refuse_overridden_address(
resource_id: str, config_scoredef: GenomicScoreDef,
) -> None:
"""Refuse a ``scores:`` entry addressing a column other than its ``id``.
A VCF score has no column address to give: it reads INFO ``<id>``, and
:func:`parse_vcf_scoredefs` takes ``col_name``/``col_index`` from the
header side unconditionally. So a ``column_name:`` (or legacy
``name:``) that differs from the ``id``, or any ``column_index:``
(``index:``), is a line the reader would ignore while the author reads
it as the score's source -- a contradiction between the config and the
header, refused like the rules above. An address EQUAL to the ``id`` is
redundant, not wrong, and is the spelling every deployed VCF resource
uses, so it passes.
"""
score_id = config_scoredef.score_id
if config_scoredef.col_index is not None:
raise score_configuration_error(
resource_id, score_id,
f"states column_index {config_scoredef.col_index}, but a VCF "
f"score has no column index: it reads the INFO field named by "
f"its 'id', so this score reads '{score_id}'. Drop the "
f"'column_index:' (or legacy 'index:') line.")
col_name = config_scoredef.col_name
if col_name is None or col_name == score_id:
return
raise score_configuration_error(
resource_id, score_id,
f"states column_name '{col_name}', but a VCF score reads the INFO "
f"field named by its 'id', so this score reads '{score_id}'. Drop "
f"the 'column_name:' (or legacy 'name:') line, or make it "
f"'{score_id}'.")
[docs]
def parse_vcf_scoredefs(
vcf_header_info: dict[str, Any] | None,
config_scoredefs: dict[str, GenomicScoreDef] | None, *,
resource_id: str,
merge: bool = False,
) -> dict[str, GenomicScoreDef]:
"""Build score definitions from a VCF header's INFO metadata.
Every INFO field the header declares becomes a score, typed through
``VCF_TYPE_CONVERSION_MAP`` and described by the header's own description.
This is why a VCF resource needs no ``scores:`` block to be usable: the
file documents its own scores.
``value_parser`` is set to ``None`` for ``Number`` of 1, ``A`` or ``R``,
because pysam already decodes those to a scalar (or to a tuple that
:func:`extract_vcf_value` indexes by allele). Every other shape keeps
``converter``, which joins a tuple on '|' -- the VCF-local convention for
a field whose arity the header does not fix.
``converter`` joins with ``str`` -- a ``Number=.``/``Type=Integer`` field
therefore reads as text, which is now what such a field DECLARES as well
(gain#1259; it used to declare the header's ``Type=`` and was the one
place the definition and the value disagreed) -- and which is what made
it the SILENT half of #630: an empty element would render as the
four-character string ``'None'``. It is not guarded here. The tuples
that reach it have already had their empty elements dropped by
:func:`extract_vcf_value`, the only route to it and the only layer
holding the record a report has to name; this parser sees a value, not
a row.
``config_scoredefs`` is what the resource's own ``scores:`` block declared,
and overrides the header for the fields it names: description, aggregators
and NA values all take the config's value when it gives one, falling back
to the header's. The value TYPE and the value PARSER are overridable only
together, and only for a field whose header declares a scalar ``Number``
(:data:`_SCALAR_VALUED_NUMBERS`). An entry that leaves ``type:`` unstated
takes neither, so it reads exactly what the header-only resource reads
(gain#1221). A field the header declares MULTI-VALUED takes neither
either, whatever ``type:`` says: it keeps the header's ``converter``,
because that converter IS the field's ``|``-join and a value type cannot
describe a tuple (gain#1233), and it keeps the ``str`` that join produces,
because a type is not merely descriptive -- it selects the histogram, and
a joined field declaring ``int`` aborted its own statistics build in
``np.isnan`` (gain#1259). An entry stating such a type is REFUSED by
:func:`_refuse_overridden_type` -- it contradicts the header, and
gain#1336 raises rather than discarding it. Column addressing is NOT
overridable -- a VCF score is its INFO key, so ``col_name``/``col_index``
always come from the header side -- which is why an entry whose ``id``
the header does not declare is a contradiction too, REFUSED by
:func:`_refuse_undeclared_id` before anything else is asked of it
(gain#1489), and why an entry that states an address OTHER than its
``id`` -- a differing ``column_name:``, or any ``column_index:`` -- is
one as well, REFUSED by :func:`_refuse_overridden_address` (gain#1498);
an address equal to the ``id`` is redundant and passes. The fourth
refusal, :func:`_refuse_genotype_arity`, is the header's own claim
rather than the config's (gain#1258). The order is id, arity, address,
type: the arity check runs BEFORE the two config-side ones because a
per-genotype field is unreadable whatever the entry states and neither
of their advised edits can change that, and the address check precedes
the type check because a wrong address is wrong whatever type is
stated. ``resource_id`` is threaded in for those four messages alone.
``merge`` decides what happens to header fields the config does not
mention: ``False`` (the default) returns only the configured scores, so
the config acts as a filter; ``True`` keeps the rest as the header
defined them. It is the resource's ``merge_vcf_scores`` setting.
"""
def converter(val: Any) -> Any:
try:
if isinstance(val, tuple):
return "|".join(map(str, val))
except TypeError:
pass
return val
vcf_scoredefs = {}
assert vcf_header_info is not None
for key, value in vcf_header_info.items():
# A scalar-valued field needs no parser at all: pysam has already
# decoded it (and a per-allele one is indexed down to one element
# before the parse). Everything else keeps ``converter``, whose
# whole job is the tuple join. The SAME set decides the config
# override below, so the two cannot drift apart.
is_scalar = value.number in _SCALAR_VALUED_NUMBERS
vcf_scoredefs[key] = GenomicScoreDef(
score_id=key,
col_name=key,
col_index=None,
desc=value.description or "",
value_parser=None if is_scalar else converter,
# ``Type=`` describes ONE element; the joined value is text
# whatever that says, so a multi-valued field declares what it
# actually holds rather than what its elements are (gain#1259).
value_type=(
VCF_TYPE_CONVERSION_MAP[value.type] if is_scalar else "str"),
na_values=(),
aggregator=None,
small_values_desc=None,
large_values_desc=None,
hist_conf=None,
)
# A configured id the header does not declare is refused before either
# loop below indexes the header by it: neither the arity nor the type
# can be asked of a field that is not there, and both indexers raised a
# bare KeyError naming neither the resource nor the block (gain#1489).
if config_scoredefs is not None:
for score in config_scoredefs:
_refuse_undeclared_id(resource_id, score, vcf_header_info)
# Every field that becomes a definition -- all of them with no config or
# a merging one, else the ones the config names -- is checked here,
# before the override loop's two config-side refusals (the docstring
# says why that order).
defined = (
vcf_scoredefs if config_scoredefs is None or merge
else config_scoredefs)
for score in defined:
_refuse_genotype_arity(resource_id, score, vcf_header_info[score])
if config_scoredefs is None:
return vcf_scoredefs
# allow overriding of vcf-generated scoredefs
scoredefs = {}
for score, config_scoredef in config_scoredefs.items():
vcf_scoredef = vcf_scoredefs[score]
_refuse_overridden_address(resource_id, config_scoredef)
# ONE rule for both, which is why neither is a ``config.x or vcf.x``
# (a ``None`` on either side means "nothing to parse", not
# "unstated", and ``str``'s falsiness is not the question either):
# the config's type AND its parser apply exactly where the parse is
# handed a scalar, and nowhere else.
#
# A field the header declares multi-valued keeps the header's
# converter, because that converter IS the field's ``|``-join:
# taking the config's parser there fed the raw tuple to
# ``int``/``float`` (a logged non-value per row) or to ``str`` (the
# tuple's repr, silently) -- gain#1233. Taking it whenever the
# config merely stated a ``type:`` is the older, wider form of the
# same mistake, which read a ``Flag`` as ``1.0`` (gain#1221).
#
# The TYPE travels with the parser rather than following the config
# on its own, which is what gain#1233 left it doing. That rested on
# a type being merely descriptive for a value that is joined text
# either way; it is not. It selects the histogram, so a joined
# field declaring ``int`` was answered with a NUMBER histogram and
# aborted its own statistics build in ``np.isnan`` (gain#1259). A
# config may still describe such a field -- ``desc``, aggregators
# and histogram config are all its own -- but not by claiming its
# value is something the join cannot produce.
config_type = config_scoredef.value_type
header_entry = vcf_header_info[score]
number = header_entry.number
is_scalar = number in _SCALAR_VALUED_NUMBERS
takes_config_type = config_type is not None and is_scalar
_refuse_overridden_type(
resource_id, score, number, header_entry.type,
config_type, is_scalar=is_scalar)
value_type = (
config_type if takes_config_type else vcf_scoredef.value_type)
value_parser = (
config_scoredef.value_parser if takes_config_type
else vcf_scoredef.value_parser)
scoredef = GenomicScoreDef(
score_id=vcf_scoredef.score_id,
desc=config_scoredef.desc or vcf_scoredef.desc,
value_type=value_type,
aggregator=config_scoredef.aggregator,
small_values_desc=config_scoredef.small_values_desc,
large_values_desc=config_scoredef.large_values_desc,
col_name=vcf_scoredef.col_name,
col_index=vcf_scoredef.col_index,
hist_conf=config_scoredef.hist_conf,
value_parser=value_parser,
na_values=config_scoredef.na_values or vcf_scoredef.na_values,
)
scoredefs[score] = scoredef
if merge:
for score, vcf_scoredef in vcf_scoredefs.items():
if score in scoredefs:
continue
scoredefs[score] = vcf_scoredef
return scoredefs