diff --git a/README.md b/README.md index 636ae26..4fedfea 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ [![Actions status](https://github.com/genomicmedlab/agct/actions/workflows/checks.yaml/badge.svg)](https://github.com/genomicmedlab/agct/actions/checks.yaml) -Drop-in replacement for the [pyliftover](https://github.com/konstantint/pyliftover) tool, using the St. Jude's [chainfile](https://docs.rs/chainfile/latest/chainfile/) crate. +Lift over positions between genomic reference assemblies, using the St. Jude's [chainfile](https://docs.rs/chainfile/latest/chainfile/) crate. Enables significantly faster chainfile loading from cold start (see `analysis/`). @@ -24,8 +24,9 @@ python3 -m pip install agct Initialize a class instance: ```python3 -from agct import Converter, Genome -c = Converter(Genome.HG38, Genome.HG19) +from agct import Converter, Assembly, Strand + +c = Converter(Assembly.HG19, Assembly.HG38) ``` > If a chainfile is unavailable locally, it's downloaded from UCSC and saved using the `wags-tails` package -- see the [wags-tails configuration instructions](https://github.com/GenomicMedLab/wags-tails?tab=readme-ov-file#configuration) for information on how to designate a non-default storage location. @@ -33,8 +34,7 @@ c = Converter(Genome.HG38, Genome.HG19) Call ``convert_coordinate()``: ```python3 -c.convert_coordinate("chr7", 140453136, "+") -# [['chr7', 140152936, ]] +c.convert_coordinate("chr7", 140453136, Strand.POSITIVE) ``` ## Development diff --git a/src/agct/__init__.py b/src/agct/__init__.py index 37ef9c3..8d6eb91 100644 --- a/src/agct/__init__.py +++ b/src/agct/__init__.py @@ -1,5 +1,6 @@ """Provide fast liftover in Python via the ``chainfile`` crate.""" -from agct.converter import Converter, Genome, Strand +from agct.assembly_registry import Assembly, get_assembly_from_refget_id +from agct.converter import Converter, Strand -__all__ = ["Converter", "Genome", "Strand"] +__all__ = ["Assembly", "Converter", "Strand", "get_assembly_from_refget_id"] diff --git a/src/agct/assembly_registry.py b/src/agct/assembly_registry.py new file mode 100644 index 0000000..7d6b320 --- /dev/null +++ b/src/agct/assembly_registry.py @@ -0,0 +1,160 @@ +"""Assembly registry. + +This module maintains a curated mapping from GA4GH VRS `SequenceReference` +refget accession identifiers (the ``"SQ."`` IDs) to human reference assemblies +and exposes a helper to look up the assembly for a given ID. Assemblies are +represented by the :class:`Assembly` enum (UCSC-style labels: ``hg19``, ``hg38``). + +Use :func:`get_assembly_from_refget_id` when you have a refget accession ID and +need to determine which assembly it belongs to. The mapping is intentionally +conservative and currently reflects internal use cases; extend it as needed. +""" + +import logging +import re +from enum import StrEnum + +_logger = logging.getLogger(__name__) + + +class Assembly(StrEnum): + """Constrain reference genome assembly values. + + We could conceivably support every UCSC chainfile offering, but for now, we'll + stick with internal use cases only. + """ + + HG38 = "hg38" + HG19 = "hg19" + + +ID_TO_REFERENCE_MAP = { + Assembly.HG38: [ + # chr1 + "SQ.Ya6Rs7DHhDeg7YaOSg1EoNi3U_nQ9SvO", + # chr2 + "SQ.pnAqCRBrTsUoBghSD1yp_jXWSmlbdh4g", + # chr3 + "SQ.Zu7h9AggXxhTaGVsy7h_EZSChSZGcmgX", + # chr4 + "SQ.HxuclGHh0XCDuF8x6yQrpHUBL7ZntAHc", + # chr5 + "SQ.aUiQCzCPZ2d0csHbMSbh2NzInhonSXwI", + # chr6 + "SQ.0iKlIQk2oZLoeOG9P1riRU6hvL5Ux8TV", + # chr7 + "SQ.F-LrLMe1SRpfUZHkQmvkVKFEGaoDeHul", + # chr8 + "SQ.209Z7zJ-mFypBEWLk4rNC6S_OxY5p7bs", + # chr9 + "SQ.KEO-4XBcm1cxeo_DIQ8_ofqGUkp4iZhI", + # chr10 + "SQ.ss8r_wB0-b9r44TQTMmVTI92884QvBiB", + # chr11 + "SQ.2NkFm8HK88MqeNkCgj78KidCAXgnsfV1", + # chr12 + "SQ.6wlJpONE3oNb4D69ULmEXhqyDZ4vwNfl", + # chr13 + "SQ._0wi-qoDrvram155UmcSC-zA5ZK4fpLT", + # chr14 + "SQ.eK4D2MosgK_ivBkgi6FVPg5UXs1bYESm", + # chr15 + "SQ.AsXvWL1-2i5U_buw6_niVIxD6zTbAuS6", + # chr16 + "SQ.yC_0RBj3fgBlvgyAuycbzdubtLxq-rE0", + # chr17 + "SQ.dLZ15tNO1Ur0IcGjwc3Sdi_0A6Yf4zm7", + # chr18 + "SQ.vWwFhJ5lQDMhh-czg06YtlWqu0lvFAZV", + # chr19 + "SQ.IIB53T8CNeJJdUqzn9V_JnRtQadwWCbl", + # chr20 + "SQ.-A1QmD_MatoqxvgVxBLZTONHz9-c7nQo", + # chr21 + "SQ.5ZUqxCmDDgN4xTRbaSjN8LwgZironmB8", + # chr22 + "SQ.7B7SHsmchAR0dFcDCuSFjJAo7tX87krQ", + # chr23 + "SQ.w0WZEvgJF0zf_P4yyTzjjv9oW1z61HHP", + # chr24/chrY + "SQ.8_liLu1aycC0tPQPFmUaGXJLDs5SbPZ5", + ], + Assembly.HG19: [ + # chr1 + "SQ.S_KjnFVz-FE7M0W6yoaUDgYxLPc1jyWU", + # chr2 + "SQ.9KdcA9ZpY1Cpvxvg8bMSLYDUpsX6GDLO", + # chr3 + "SQ.VNBualIltAyi2AI_uXcKU7M9XUOuA7MS", + # chr4 + "SQ.iy7Zfceb5_VGtTQzJ-v5JpPbpeifHD_V", + # chr5 + "SQ.vbjOdMfHJvTjK_nqvFvpaSKhZillW0SX", + # chr6 + "SQ.KqaUhJMW3CDjhoVtBetdEKT1n6hM-7Ek", + # chr7 + "SQ.IW78mgV5Cqf6M24hy52hPjyyo5tCCd86", + # chr8 + "SQ.tTm7wmhz0G4lpt8wPspcNkAD_qiminj6", + # chr9 + "SQ.HBckYGQ4wYG9APHLpjoQ9UUe9v7NxExt", + # chr10 + "SQ.-BOZ8Esn8J88qDwNiSEwUr5425UXdiGX", + # chr11 + "SQ.XXi2_O1ly-CCOi3HP5TypAw7LtC6niFG", + # chr12 + "SQ.105bBysLoDFQHhajooTAUyUkNiZ8LJEH", + # chr13 + "SQ.Ewb9qlgTqN6e_XQiRVYpoUfZJHXeiUfH", + # chr14 + "SQ.5Ji6FGEKfejK1U6BMScqrdKJK8GqmIGf", + # chr15 + "SQ.zIMZb3Ft7RdWa5XYq0PxIlezLY2ccCgt", + # chr16 + "SQ.W6wLoIFOn4G7cjopxPxYNk2lcEqhLQFb", + # chr17 + "SQ.AjWXsI7AkTK35XW9pgd3UbjpC3MAevlz", + # chr18 + "SQ.BTj4BDaaHYoPhD3oY2GdwC_l0uqZ92UD", + # chr19 + "SQ.ItRDD47aMoioDCNW_occY5fWKZBKlxCX", + # chr20 + "SQ.iy_UbUrvECxFRX5LPTH_KPojdlT7BKsf", + # chr21 + "SQ.LpTaNW-hwuY_yARP0rtarCnpCQLkgVCg", + # chr22 + "SQ.XOgHwwR3Upfp5sZYk6ZKzvV25a4RBVu8", + # chr23 + "SQ.v7noePfnNpK8ghYXEqZ9NukMXW7YeNsm", + # chr24/chrY + "SQ.BT7QyW5iXaX_1PSX-msSGYsqRdMKqkj-", + ], +} + + +_error_pattern = re.compile(r"^SQ.[0-9A-Za-z_\\-]{32}$") + + +def get_assembly_from_refget_id(refget_accession: str) -> Assembly | None: + """Given a GA4GH SequenceReference refget accession ID, get back its corresponding reference genome + + .. code-block:: pycon + + >>> from agct.assembly_registry import get_assembly_from_refget_id + >>> get_assembly_from_refget_id("SQ.pnAqCRBrTsUoBghSD1yp_jXWSmlbdh4g") + Genome.HG38 + + Use to get an appropriate converter from a referenced GA4GH variation object. + + :param refget_accession: sequence reference (must start with `"SQ."`) + :return: a reference assembly if recognized + :raise ValueError: if input appears to be in an invalid format for a refget accession ID + """ + for reference_genome, chromosome_accession_ids in ID_TO_REFERENCE_MAP.items(): + if refget_accession in chromosome_accession_ids: + return reference_genome + if not re.match(_error_pattern, refget_accession): + msg = f"refget accession ID must be in format 'SQ.ABCDEFGHIJKLMNOPQRSTUVWXYZ123456'; got {refget_accession}" + _logger.error(msg) + raise ValueError(msg) + return None diff --git a/src/agct/converter.py b/src/agct/converter.py index 1252272..747b8fd 100644 --- a/src/agct/converter.py +++ b/src/agct/converter.py @@ -10,6 +10,7 @@ from wags_tails.utils.storage import get_data_dir import agct._core as _core +from agct.assembly_registry import Assembly _logger = logging.getLogger(__name__) @@ -21,17 +22,6 @@ class Strand(str, Enum): NEGATIVE = "-" -class Genome(str, Enum): - """Constrain genome values. - - We could conceivably support every UCSC chainfile offering, but for now, we'll - stick with internal use cases only. - """ - - HG38 = "hg38" - HG19 = "hg19" - - class Converter: """Chainfile-based liftover provider for a single sequence to sequence association. @@ -39,50 +29,48 @@ class Converter: def __init__( self, - from_db: Genome | str | None = None, - to_db: Genome | str | None = None, + from_assembly: Assembly | None = None, + to_assembly: Assembly | None = None, chainfile: str | None = None, ) -> None: """Initialize liftover instance. - :param from_db: database name, e.g. ``"hg19"``. Must be different than ``to_db`` - If ``chainfile`` is provided, will ignore this argument - :param to_db: database name, e.g. ``"hg38"``. Must be different than ``from_db`` - If ``chainfile`` is provided, will ignore this argument + Initialize with either ``from_assembly`` and ``to_assembly``, or directly with a chainfile. + + * If using assembly params, both must be provided, and they must be different + * If using assembly params, the ``wags-tails`` library will be used to locate and, if + necessary, acquire the pertinent chainfile from the UCSC web server. See + the `wags-tails documentation `_ for more info. + * If ``chainfile`` arg is provided, all other args are ignored. + + :param from_assembly: Assembly name, e.g. ```` + :param to_assembly: database name, e.g. ``"hg38"``. :param chainfile: Path to chainfile - If not provided, must provide both ``from_db`` and ``to_db`` so that - ``wags-tails`` can download the corresponding chainfile :raise ValueError: if required arguments are not passed or are invalid :raise FileNotFoundError: if unable to open corresponding chainfile :raise _core.ChainfileError: if unable to read chainfile (i.e. it's invalid) """ if not chainfile: - if from_db is None or to_db is None: - msg = "Must provide both `from_db` and `to_db`" + if from_assembly is None or to_assembly is None: + msg = "Must provide both `from_assembly` and `to_assembly`" raise ValueError(msg) - if from_db == to_db: + if from_assembly == to_assembly: msg = "Liftover must be to/from different sources." raise ValueError(msg) - if isinstance(from_db, str): - try: - from_db = Genome(from_db) - except ValueError as e: - msg = f"Unable to coerce from_db value '{from_db}' to a known reference genome: {list(Genome)}" - raise ValueError(msg) from e - if isinstance(to_db, str): - try: - to_db = Genome(to_db) - except ValueError as e: - msg = f"Unable to coerce to_db value '{to_db}' to a known reference genome: {list(Genome)}" - raise ValueError(msg) from e + try: + file_prefix = f"chainfile_{from_assembly.value}_to_{to_assembly.value}" + except AttributeError as e: + msg = f"Assembly args must be instance of `agct.assembly_registry.Genome`, instead got from_assembly={from_assembly} and to_assembly={to_assembly}" + _logger.exception(msg) + raise ValueError(msg) from e data_handler = CustomData( - f"chainfile_{from_db.value}_to_{to_db.value}", + file_prefix, "chain", lambda: "", - self._download_function_builder(from_db, to_db), + self._download_function_builder(from_assembly, to_assembly), data_dir=get_data_dir() / "ucsc-chainfile", ) file, _ = data_handler.get_latest() @@ -98,15 +86,17 @@ def __init__( raise @staticmethod - def _download_function_builder(from_db: Genome, to_db: Genome) -> Callable: + def _download_function_builder( + from_assembly: Assembly, to_assembly: Assembly + ) -> Callable: """Build downloader function for chainfile corresponding to source/destination params. Wags-Tails' custom data handler takes a downloader callback function. We construct it here, curried with from/to values in the download URL. - :param from_db: genome lifting from - :param to_db: genome lifting to + :param from_assembly: genome lifting from + :param to_assembly: genome lifting to :return: Function that downloads appropriate chainfile from UCSC """ @@ -116,7 +106,7 @@ def _download_data(version: str, file: Path) -> None: # noqa: ARG001 :param version: not used :param file: path to save file to """ - url = f"https://hgdownload.soe.ucsc.edu/goldenPath/{from_db.value}/liftOver/{from_db.value}To{to_db.value.title()}.over.chain.gz" + url = f"https://hgdownload.soe.ucsc.edu/goldenPath/{from_assembly.value}/liftOver/{from_assembly.value}To{to_assembly.value.title()}.over.chain.gz" download_http(url, file, handler=handle_gzip) return _download_data @@ -130,9 +120,9 @@ def convert_coordinate( .. code-block:: python - from agct import Converter, Strand + from agct import Converter, Strand, Assembly - c = Converter("hg19", "hg38") + c = Converter(Assembly.HG19, Assembly.HG38) c.convert_coordinate("chr7", 140453136, Strand.POSITIVE) # returns [['chr7', 140753336, '+']] diff --git a/tests/test_assembly_registry.py b/tests/test_assembly_registry.py new file mode 100644 index 0000000..283e58b --- /dev/null +++ b/tests/test_assembly_registry.py @@ -0,0 +1,39 @@ +import pytest + +from agct.assembly_registry import Assembly, get_assembly_from_refget_id + + +def test_assembly_fetcher(): + assert ( + get_assembly_from_refget_id("SQ.ss8r_wB0-b9r44TQTMmVTI92884QvBiB") + == Assembly.HG38 + ) + assert ( + get_assembly_from_refget_id("SQ.-BOZ8Esn8J88qDwNiSEwUr5425UXdiGX") + == Assembly.HG19 + ) + assert ( + get_assembly_from_refget_id("SQ.0iKlIQk2oZLoeOG9P1riRU6hvL5Ux8TV") + == Assembly.HG38 + ) + assert ( + get_assembly_from_refget_id("SQ.KqaUhJMW3CDjhoVtBetdEKT1n6hM-7Ek") + == Assembly.HG19 + ) + assert ( + get_assembly_from_refget_id("SQ.LpTaNW-hwuY_yARP0rtarCnpCQLkgVCg") + == Assembly.HG19 + ) + + +def test_unknown_sequences(): + assert get_assembly_from_refget_id("SQ.KqX7SGMq4GLY21TCAHPMWfNy-cFdT44h") is None + + +def test_invalid_input(): + input_string = "ga4gh:SQ.KqaUhJMW3CDjhoVtBetdEKT1n6hM-7Ek" + with pytest.raises( + ValueError, + match=f"refget accession ID must be in format 'SQ.ABCDEFGHIJKLMNOPQRSTUVWXYZ123456'; got {input_string}", + ): + get_assembly_from_refget_id(input_string) diff --git a/tests/test_converter.py b/tests/test_converter.py index d35f738..201d986 100644 --- a/tests/test_converter.py +++ b/tests/test_converter.py @@ -5,7 +5,7 @@ import pytest -from agct import Converter, Genome +from agct import Assembly, Converter def test_valid(data_dir: Path): @@ -17,18 +17,20 @@ def test_valid(data_dir: Path): def test_invalid(): """Test invalid initialization""" - with pytest.raises(ValueError, match="Must provide both `from_db` and `to_db`"): + with pytest.raises( + ValueError, match="Must provide both `from_assembly` and `to_assembly`" + ): Converter() with pytest.raises( ValueError, match=re.escape("Liftover must be to/from different sources.") ): - Converter(Genome.HG19, Genome.HG19) + Converter(Assembly.HG19, Assembly.HG19) with pytest.raises( ValueError, match=re.escape( - "Unable to coerce to_db value 'hg18' to a known reference genome: [, ]" + "Assembly args must be instance of `agct.assembly_registry.Genome`, instead got from_assembly=hg19 and to_assembly=hg18" ), ): - Converter(Genome.HG19, "hg18") + Converter(Assembly.HG19, "hg18") diff --git a/tests/test_liftover.py b/tests/test_liftover.py index 4e0d403..566ef09 100644 --- a/tests/test_liftover.py +++ b/tests/test_liftover.py @@ -1,11 +1,11 @@ """Run liftover tests.""" -from agct import Converter, Genome, Strand +from agct import Assembly, Converter, Strand def test_hg19_to_hg38(): """Test hg19 to hg38 lifter.""" - converter = Converter(Genome.HG19, Genome.HG38) + converter = Converter(Assembly.HG19, Assembly.HG38) result = converter.convert_coordinate("chr7", 140439611) assert len(result) == 1 @@ -34,7 +34,7 @@ def test_hg19_to_hg38(): def test_hg38_to_hg19(): """Test hg38 to hg19 lifter.""" - converter = Converter(Genome.HG38, Genome.HG19) + converter = Converter(Assembly.HG38, Assembly.HG19) result = converter.convert_coordinate("chr7", 140739811) assert len(result) == 1