Source code for pysec2pri.parsers.ensembl

"""Ensembl core flat-file parser for secondary-to-primary gene identifier mappings.

Reads the MySQL flat-file dumps Ensembl publishes per release/species:

- ``stable_id_event.txt.gz``: cumulative old (secondary) -> new (primary)
  stable-ID history. Each release's table already holds the *entire* history
  up to that release, so one release's mappings describe the full state of
  Ensembl at that point in time (no cross-release chain resolution needed
  here -- :mod:`pysec2pri.consolidate` handles cross-release temporality).
- ``mapping_session.txt.gz``: maps each ``mapping_session_id`` to the
  ``created`` timestamp used for the per-row ``mapping_date``.
- ``gene.txt.gz`` / ``xref.txt.gz`` / ``external_synonym.txt.gz``: the
  current gene set, its display labels, and symbol synonyms.

This parser extracts:
1. ID-to-ID mappings: retired/old Ensembl gene IDs -> current Ensembl gene IDs
2. Label-to-label mappings: gene symbol synonyms -> current display label

Uses SSSOM-compliant MappingSet classes with cardinality computation.
"""

from __future__ import annotations

from collections.abc import Iterable
from pathlib import Path
from typing import Any

import polars as pl
from sssom_schema import Mapping

from pysec2pri.logging import logger
from pysec2pri.parsers.base import (
    ALL_SPECIES,
    WITHDRAWN_ENTRY,
    WITHDRAWN_ENTRY_LABEL,
    BaseMappingSet,
    BaseParser,
)

# Column layouts verified against the Ensembl release 115 core schema

_STABLE_ID_EVENT_COLUMNS = [
    "old_stable_id",
    "old_version",
    "new_stable_id",
    "new_version",
    "mapping_session_id",
    "type",
    "score",
]
_MAPPING_SESSION_COLUMNS = [
    "mapping_session_id",
    "old_db_name",
    "new_db_name",
    "old_release",
    "new_release",
    "old_assembly",
    "new_assembly",
    "created",
]
_GENE_COLUMNS = [
    "gene_id",
    "biotype",
    "analysis_id",
    "seq_region_id",
    "seq_region_start",
    "seq_region_end",
    "seq_region_strand",
    "display_xref_id",
    "source",
    "description",
    "is_current",
    "canonical_transcript_id",
    "stable_id",
    "version",
    "created_date",
    "modified_date",
]
_XREF_COLUMNS = [
    "xref_id",
    "external_db_id",
    "dbprimary_acc",
    "display_label",
    "version",
    "description",
    "info_type",
    "info_text",
]
_EXTERNAL_SYNONYM_COLUMNS = ["xref_id", "synonym"]

# Explicit dtypes where the natural column type (e.g. stable_id_event.score,
# mostly 0 with rare floats) would otherwise be mis-inferred from a sample.
_DType = type[pl.DataType]
_STABLE_ID_EVENT_DTYPES: dict[str, _DType] = {
    "old_version": pl.Int64,
    "new_version": pl.Int64,
    "mapping_session_id": pl.Int64,
    "score": pl.Float64,
}
_GENE_DTYPES: dict[str, _DType] = {
    "display_xref_id": pl.Int64,
    "is_current": pl.Int64,
}
_XREF_DTYPES: dict[str, _DType] = {"xref_id": pl.Int64}
_EXTERNAL_SYNONYM_DTYPES: dict[str, _DType] = {"xref_id": pl.Int64}


def _scan_ensembl_tsv(
    path: Path, columns: list[str], schema_overrides: dict[str, _DType] | None = None
) -> pl.LazyFrame:
    """Lazily scan a header-less Ensembl flat-file dump."""
    return pl.scan_csv(
        path,
        separator="\t",
        has_header=False,
        new_columns=columns,
        null_values=["\\N"],
        quote_char=None,
        infer_schema_length=10000,
        schema_overrides=schema_overrides,
    )


def _ensembl_date_to_iso(value: object) -> str | None:
    """Convert a ``mapping_session.created`` timestamp to ISO ``YYYY-MM-DD``.

    Args:
        value: Raw ``created`` cell, e.g. ``"2002-09-03 11:19:48"``.

    Returns:
        ISO-8601 date string, or ``None`` if *value* is missing/malformed.
    """
    if value is None:
        return None
    text = str(value).strip()
    if not text:
        return None
    return text.split(" ", 1)[0]


def _clean_assembly(value: object) -> str | None:
    r"""Return a ``mapping_session`` assembly cell stripped, or ``None`` if empty.

    Args:
        value: Raw ``old_assembly``/``new_assembly`` cell, e.g. ``"GRCh38"``
            (``\N`` nulls already arrive as ``None`` from the scan).

    Returns:
        The assembly name, or ``None`` when missing/blank.
    """
    if value is None:
        return None
    text = str(value).strip()
    return text or None


def _load_session_meta(
    mapping_session_path: Path | str | None,
) -> tuple[dict[str, str | None], dict[str, tuple[str | None, str | None]]]:
    """Load per-session mapping dates and genome builds from ``mapping_session``.

    Args:
        mapping_session_path: Optional path to ``mapping_session.txt``.

    Returns:
        ``(date_by_session, assembly_by_session)``. Each assembly value is an
        ``(old_assembly, new_assembly)`` pair -- the true builds the secondary
        (old) and primary (current) IDs of that session's mappings belonged to
        (issue #51). Both dicts are empty when no path is given.
    """
    date_by_session: dict[str, str | None] = {}
    assembly_by_session: dict[str, tuple[str | None, str | None]] = {}
    if mapping_session_path is None:
        return date_by_session, assembly_by_session
    sessions = _scan_ensembl_tsv(Path(mapping_session_path), _MAPPING_SESSION_COLUMNS).collect()
    for row in sessions.iter_rows(named=True):
        session_id = str(row["mapping_session_id"])
        date_by_session[session_id] = _ensembl_date_to_iso(row["created"])
        assembly_by_session[session_id] = (
            _clean_assembly(row.get("old_assembly")),
            _clean_assembly(row.get("new_assembly")),
        )
    return date_by_session, assembly_by_session


[docs] class EnsemblParser(BaseParser): """Parser for Ensembl core flat-file dumps using Polars. Each release's ``stable_id_event`` table is cumulative, so a single parse describes the *whole state* of Ensembl gene IDs at that release (no chain-walking across releases). Returns: - IdMappingSet for ID-to-ID mappings (retired Ensembl gene IDs) - LabelMappingSet for symbol mappings (external gene synonyms) """ datasource_name = "ensembl" def __init__( self, version: str | None = None, show_progress: bool = True, species: str | int = 9606, ) -> None: """Initialize the Ensembl parser. Args: version: Ensembl release number (e.g. ``"115"``). show_progress: Whether to show progress bars during parsing. species: Canonical NCBI taxon ID of the species being parsed. Unlike NCBI, this never filters rows (the input files are already species-specific downloads) -- it only labels the output (see :meth:`_product_slug`). """ super().__init__(version=version, show_progress=show_progress) self.species = species def _product_slug(self, **overrides: Any) -> str | None: """Species slug, suffixed with ``consolidate`` for label-history runs. Label history is a distinct data product at the same release/species, so it gets an IRI segment. Args: **overrides: Forwarded to the base implementation. """ base = super()._product_slug(**overrides) species_slug = str(base) if base is not None else None if getattr(self, "_is_label_history", False): return f"{species_slug}/consolidate" if species_slug else "consolidate" return species_slug
[docs] def parse( self, stable_id_event_path: Path | str | None = None, mapping_session_path: Path | str | None = None, gene_path: Path | str | None = None, ) -> BaseMappingSet: """Parse ``stable_id_event`` (+ ``mapping_session``) into an IdMappingSet. Args: stable_id_event_path: Path to ``stable_id_event.txt`` (can be ``.gz`` compressed). mapping_session_path: Optional path to ``mapping_session.txt``, used to resolve each row's ``mapping_date``. When omitted, rows fall back to the set-level release date. gene_path: Optional path to ``gene.txt``. When supplied, ``_primary_ids`` is populated with every current gene ID. Returns: IdMappingSet with computed cardinalities based on IDs. """ if stable_id_event_path is None: raise ValueError("stable_id_event_path must not be None") stable_id_event_path = Path(stable_id_event_path) self._resolve_version(stable_id_event_path) mappings = self._parse_stable_id_event(stable_id_event_path, mapping_session_path) primary_ids = self._extract_primary_ids(Path(gene_path)) if gene_path is not None else None return self.create_mapping_set(mappings, mapping_type="id", primary_ids=primary_ids)
[docs] def parse_labels( self, gene_path: Path | str | None = None, xref_path: Path | str | None = None, external_synonym_path: Path | str | None = None, ) -> BaseMappingSet: """Parse external gene synonyms into a LabelMappingSet. A release's files state a gene's *current* symbol and its synonyms, but not the symbols it previously had, so a single release yields alias mappings only. Renames are recovered by consolidating across releases. Args: gene_path: Path to ``gene.txt``. xref_path: Path to ``xref.txt``. external_synonym_path: Optional path to ``external_synonym.txt``. When omitted, the returned set carries only the full primary label set (no synonym mappings). Returns: LabelMappingSet with computed cardinalities based on labels. """ if gene_path is None or xref_path is None: raise ValueError("gene_path and xref_path must not be None") gene_path = Path(gene_path) xref_path = Path(xref_path) self._resolve_version(gene_path) logger.warning( "%s releases carry no previous-symbol field, so this label mapping set holds " "alias mappings only. Consolidate to also recover symbol renames.", self.datasource_name.upper(), ) mappings: list[Mapping] = [] if external_synonym_path is not None: mappings = self._parse_external_synonyms( gene_path, xref_path, Path(external_synonym_path) ) return self.create_mapping_set( mappings, mapping_type="label", primary_labels=self._extract_primary_labels(gene_path, xref_path), primary_ids=self._extract_primary_ids(gene_path), )
[docs] def parse_primary_ids(self, gene_path: Path | str | None = None) -> BaseMappingSet: """Return a mapping set containing the full list of current Ensembl gene IDs. Args: gene_path: Path to ``gene.txt`` (can be ``.gz`` compressed). Returns: :class:`~pysec2pri.parsers.base.IdMappingSet` with ``_primary_ids`` populated. """ if gene_path is None: raise ValueError("gene_path must not be None") gene_path = Path(gene_path) self._resolve_version(gene_path) return self.create_mapping_set( [], mapping_type="id", primary_ids=self._extract_primary_ids(gene_path) )
[docs] def parse_primary_labels( self, gene_path: Path | str | None = None, xref_path: Path | str | None = None, ) -> BaseMappingSet: """Return a mapping set containing the full list of current Ensembl gene labels. Args: gene_path: Path to ``gene.txt``. xref_path: Path to ``xref.txt``. Returns: :class:`~pysec2pri.parsers.base.LabelMappingSet` with ``_primary_labels`` populated. """ if gene_path is None or xref_path is None: raise ValueError("gene_path and xref_path must not be None") gene_path = Path(gene_path) xref_path = Path(xref_path) self._resolve_version(gene_path) return self.create_mapping_set( [], mapping_type="label", primary_labels=self._extract_primary_labels(gene_path, xref_path), )
[docs] def parse_all( self, stable_id_event_path: Path | str | None, mapping_session_path: Path | str | None, gene_path: Path | str | None, xref_path: Path | str | None, external_synonym_path: Path | str | None, ) -> tuple[BaseMappingSet, BaseMappingSet]: """Parse the full set of Ensembl core flat files. Args: stable_id_event_path: Path to ``stable_id_event.txt``. mapping_session_path: Path to ``mapping_session.txt``. gene_path: Path to ``gene.txt``. xref_path: Path to ``xref.txt``. external_synonym_path: Path to ``external_synonym.txt``. Returns: Tuple of (IdMappingSet, LabelMappingSet). """ id_mappings = self.parse(stable_id_event_path, mapping_session_path, gene_path=gene_path) label_mappings = self.parse_labels(gene_path, xref_path, external_synonym_path) return id_mappings, label_mappings
def _parse_stable_id_event( self, file_path: Path, mapping_session_path: Path | str | None, ) -> list[Mapping]: """Parse ``stable_id_event`` for gene ID-to-ID mappings. Args: file_path: Path to ``stable_id_event.txt``. mapping_session_path: Optional path to ``mapping_session.txt``. Returns: List of SSSOM Mapping objects. """ try: df = ( _scan_ensembl_tsv(file_path, _STABLE_ID_EVENT_COLUMNS, _STABLE_ID_EVENT_DTYPES) .filter(pl.col("type") == "gene") .collect() ) except pl.exceptions.NoDataError: return [] if df.is_empty(): return [] date_by_session, assembly_by_session = _load_session_meta(mapping_session_path) m_meta = self.get_mapping_metadata() fixed_base = self._fixed_mapping_fields() replaced: set[str] = set( df.filter( pl.col("old_stable_id").is_not_null() & pl.col("new_stable_id").is_not_null() & (pl.col("old_stable_id") != pl.col("new_stable_id")) )["old_stable_id"].to_list() ) rows_data: list[dict[str, Any]] = [] for row in df.iter_rows(named=True): old_id = row.get("old_stable_id") new_id = row.get("new_stable_id") # No old (secondary) ID: a brand-new gene, not a sec->pri mapping. # old == new: not an ID change. if not old_id or old_id == new_id: continue if not new_id and old_id in replaced: continue session_id = str(row.get("mapping_session_id") or "") old_assembly, new_assembly = assembly_by_session.get(session_id, (None, None)) rows_data.append( self._stable_id_event_row( old_id=old_id, new_id=new_id, score=row.get("score"), mapping_date=date_by_session.get(session_id), old_assembly=old_assembly, new_assembly=new_assembly, m_meta=m_meta, ) ) return self._build_mappings( rows_data, fixed_base, desc="Processing stable_id_event", total=len(rows_data) ) def _stable_id_event_row( self, *, old_id: str, new_id: str | None, score: object, mapping_date: str | None, old_assembly: str | None, new_assembly: str | None, m_meta: dict[str, Any], ) -> dict[str, Any]: """Build one ``stable_id_event`` mapping row (a retirement or a rename). ``subject_source_version``/``object_source_version`` carry the analyzed release (set at the mapping-set level, same as every other datasource); they don't carry genome build. When a rename spans an assembly change, that's noted in ``comment`` instead, so it doesn't compete with the harmonized source-version semantics. Args: old_id: Bare old stable ID (the secondary). new_id: Bare new stable ID, or falsy for a retirement. score: ``stable_id_event.score`` cell (used as confidence). mapping_date: Resolved per-row mapping date, if any. old_assembly: Build the old ID belonged to, if known. new_assembly: Build the new ID belongs to, if known. m_meta: Mapping-level config metadata. Returns: A row dict for :meth:`_build_mappings` (which fills in ``record_id`` from ``subject_id``/``object_id``). """ if not new_id: return { "subject_id": f"ENSEMBL:{old_id}", "object_id": WITHDRAWN_ENTRY, "object_label": WITHDRAWN_ENTRY_LABEL, "predicate_id": "oboInOwl:consider", "mapping_date": mapping_date, } fields: dict[str, Any] = { "subject_id": f"ENSEMBL:{old_id}", "object_id": f"ENSEMBL:{new_id}", "predicate_id": m_meta["predicate_id"], "predicate_label": m_meta.get("predicate_label"), "mapping_date": mapping_date, } if score and float(score) > 0: # type: ignore fields["confidence"] = float(score) # type: ignore if old_assembly and new_assembly and old_assembly != new_assembly: fields["comment"] = f"Assembly changed from {old_assembly} to {new_assembly}." return fields def _parse_external_synonyms( self, gene_path: Path, xref_path: Path, external_synonym_path: Path, ) -> list[Mapping]: """Parse ``external_synonym`` joined through genes' display xref. Only synonyms attached to the *same* xref a gene displays as its canonical label are picked up (there is no ``object_xref`` table in the required input list to resolve arbitrary gene<->xref links). Args: gene_path: Path to ``gene.txt``. xref_path: Path to ``xref.txt``. external_synonym_path: Path to ``external_synonym.txt``. Returns: List of SSSOM Mapping objects for synonym mappings. """ genes = ( _scan_ensembl_tsv(gene_path, _GENE_COLUMNS, _GENE_DTYPES) .filter(pl.col("is_current") == 1) .select(["display_xref_id", "stable_id"]) ) xrefs = _scan_ensembl_tsv(xref_path, _XREF_COLUMNS, _XREF_DTYPES).select( ["xref_id", "display_label"] ) synonyms = _scan_ensembl_tsv( external_synonym_path, _EXTERNAL_SYNONYM_COLUMNS, _EXTERNAL_SYNONYM_DTYPES ) df = ( genes.join(xrefs, left_on="display_xref_id", right_on="xref_id", how="inner") .join(synonyms, left_on="display_xref_id", right_on="xref_id", how="inner") .select(["stable_id", "display_label", "synonym"]) .drop_nulls() .collect() ) if df.is_empty(): return [] fixed = self._fixed_mapping_fields() rows_data: list[dict[str, Any]] = [] for stable_id, display_label, synonym in df.rows(): curie_id = f"ENSEMBL:{stable_id}" rows_data.append( { "object_id": curie_id, "subject_label": synonym, "subject_type": "rdfs literal", "object_label": display_label, "_label_type": "alias", "comment": "Ensembl gene external synonym.", } ) return self._build_mappings( rows_data, fixed, desc="Processing external_synonym", total=len(rows_data) ) def _extract_primary_ids(self, file_path: Path) -> set[str]: """Extract all current Ensembl gene IDs from ``gene.txt``. Args: file_path: Path to ``gene.txt``. Returns: Set of ``ENSEMBL:<stable_id>`` CURIEs. """ df = ( _scan_ensembl_tsv(file_path, _GENE_COLUMNS, _GENE_DTYPES) .filter(pl.col("is_current") == 1) .select("stable_id") .collect() ) return {f"ENSEMBL:{v}" for v in df["stable_id"].drop_nulls().to_list()} def _extract_primary_labels(self, gene_path: Path, xref_path: Path) -> dict[str, set[str]]: """Extract all current gene display labels from ``gene.txt`` + ``xref.txt``. Returns a ``dict`` mapping each label text to the set of primary ``ENSEMBL:<stable_id>`` IDs that carry it. Args: gene_path: Path to ``gene.txt``. xref_path: Path to ``xref.txt``. Returns: ``dict[label, set[ENSEMBL:<stable_id>]]`` """ genes = ( _scan_ensembl_tsv(gene_path, _GENE_COLUMNS, _GENE_DTYPES) .filter(pl.col("is_current") == 1) .select(["display_xref_id", "stable_id"]) ) xrefs = _scan_ensembl_tsv(xref_path, _XREF_COLUMNS, _XREF_DTYPES).select( ["xref_id", "display_label"] ) df = ( genes.join(xrefs, left_on="display_xref_id", right_on="xref_id", how="inner") .select(["stable_id", "display_label"]) .drop_nulls() .collect() ) result: dict[str, set[str]] = {} for stable_id, label in df.rows(): result.setdefault(str(label), set()).add(f"ENSEMBL:{stable_id}") return result
[docs] def current_label_snapshot( self, gene_path: Path | str, xref_path: Path | str ) -> dict[str, str]: """Return ``{bare stable_id -> current display label}`` for every current gene. Unlike :meth:`_extract_primary_labels` (keyed by label, for ambiguity detection and ``to_pri_labels()``), this is keyed by gene -- the natural shape for diffing one release's label snapshot against the next (see :func:`pysec2pri.consolidate.build_label_history`). Args: gene_path: Path to ``gene.txt``. xref_path: Path to ``xref.txt``. Returns: ``dict[stable_id, label]`` (bare stable IDs, no ``ENSEMBL:`` prefix). """ genes = ( _scan_ensembl_tsv(Path(gene_path), _GENE_COLUMNS, _GENE_DTYPES) .filter(pl.col("is_current") == 1) .select(["display_xref_id", "stable_id"]) ) xrefs = _scan_ensembl_tsv(Path(xref_path), _XREF_COLUMNS, _XREF_DTYPES).select( ["xref_id", "display_label"] ) df = ( genes.join(xrefs, left_on="display_xref_id", right_on="xref_id", how="inner") .select(["stable_id", "display_label"]) .drop_nulls() .collect() ) return {str(stable_id): str(label) for stable_id, label in df.rows()}
[docs] def parse_label_history( self, transitions: Iterable[tuple[str, str, str, str | None]] ) -> BaseMappingSet: """Build a LabelMappingSet from precomputed previous->current label transitions. Args: transitions: Iterable of ``(stable_id, prev_label, curr_label, mapping_date)`` tuples, one per gene whose display label changed between two releases (see :func:`pysec2pri.consolidate.build_label_history`). Returns: LabelMappingSet with ``IAO:0100001`` ("term replaced by") mappings. """ self._is_label_history = True try: fixed = self._fixed_mapping_fields() rows_data: list[dict[str, Any]] = [] for stable_id, prev_label, curr_label, mapping_date in transitions: curie_id = f"ENSEMBL:{stable_id}" rows_data.append( { "object_id": curie_id, "subject_label": prev_label, "subject_type": "rdfs literal", "object_label": curr_label, "_label_type": "previous", "mapping_date": mapping_date, "comment": "Derived from snapshots across releases.", } ) mappings = self._build_mappings( rows_data, fixed, desc="Building label history", total=len(rows_data) ) return self.create_mapping_set(mappings, mapping_type="label") finally: self._is_label_history = False
__all__ = ["ALL_SPECIES", "EnsemblParser"]