From 8321c7d5f72c1da953488435e6169e2b65b3cea0 Mon Sep 17 00:00:00 2001 From: Mihai Todor Date: Thu, 5 Oct 2023 18:08:42 +0300 Subject: [PATCH 1/5] Update tests --- .../expected_outputs/find_population_dx_implications/1.json | 2 +- .../expected_outputs/find_population_dx_implications/2.json | 2 +- .../expected_outputs/find_population_dx_implications/3.json | 2 +- .../expected_outputs/find_population_dx_implications/4.json | 2 +- .../find_population_specific_haplotypes/1.json | 2 +- .../find_population_specific_haplotypes/2.json | 2 +- .../find_population_specific_haplotypes/3.json | 2 +- .../find_population_specific_haplotypes/4.json | 6 +++--- .../find_population_specific_variants/1.json | 2 +- .../find_population_specific_variants/2.json | 2 +- .../find_population_structural_intersecting_variants/1.json | 2 +- .../find_population_structural_intersecting_variants/2.json | 4 ++-- .../find_population_structural_intersecting_variants/3.json | 2 +- .../find_population_structural_subsuming_variants/1.json | 2 +- .../find_population_structural_subsuming_variants/2.json | 4 ++-- .../expected_outputs/find_population_tx_implications/2.json | 2 +- .../expected_outputs/find_population_tx_implications/3.json | 2 +- .../expected_outputs/find_population_tx_implications/4.json | 2 +- .../expected_outputs/find_population_tx_implications/5.json | 2 +- 19 files changed, 23 insertions(+), 23 deletions(-) diff --git a/tests/expected_outputs/find_population_dx_implications/1.json b/tests/expected_outputs/find_population_dx_implications/1.json index 90b8d1d06..c524713e8 100644 --- a/tests/expected_outputs/find_population_dx_implications/1.json +++ b/tests/expected_outputs/find_population_dx_implications/1.json @@ -13,7 +13,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_dx_implications/2.json b/tests/expected_outputs/find_population_dx_implications/2.json index b2b592014..b4823e5bd 100644 --- a/tests/expected_outputs/find_population_dx_implications/2.json +++ b/tests/expected_outputs/find_population_dx_implications/2.json @@ -13,7 +13,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } } ] diff --git a/tests/expected_outputs/find_population_dx_implications/3.json b/tests/expected_outputs/find_population_dx_implications/3.json index 90b8d1d06..c524713e8 100644 --- a/tests/expected_outputs/find_population_dx_implications/3.json +++ b/tests/expected_outputs/find_population_dx_implications/3.json @@ -13,7 +13,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_dx_implications/4.json b/tests/expected_outputs/find_population_dx_implications/4.json index 9dea191bd..98e3d8ed1 100644 --- a/tests/expected_outputs/find_population_dx_implications/4.json +++ b/tests/expected_outputs/find_population_dx_implications/4.json @@ -13,7 +13,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_specific_haplotypes/1.json b/tests/expected_outputs/find_population_specific_haplotypes/1.json index 5a0089e83..7a30e6d86 100644 --- a/tests/expected_outputs/find_population_specific_haplotypes/1.json +++ b/tests/expected_outputs/find_population_specific_haplotypes/1.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_specific_haplotypes/2.json b/tests/expected_outputs/find_population_specific_haplotypes/2.json index eb17abda0..d0e3ead2d 100644 --- a/tests/expected_outputs/find_population_specific_haplotypes/2.json +++ b/tests/expected_outputs/find_population_specific_haplotypes/2.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_specific_haplotypes/3.json b/tests/expected_outputs/find_population_specific_haplotypes/3.json index fb743c1d9..62f4d9b37 100644 --- a/tests/expected_outputs/find_population_specific_haplotypes/3.json +++ b/tests/expected_outputs/find_population_specific_haplotypes/3.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_specific_haplotypes/4.json b/tests/expected_outputs/find_population_specific_haplotypes/4.json index e1ae8333a..28be67cab 100644 --- a/tests/expected_outputs/find_population_specific_haplotypes/4.json +++ b/tests/expected_outputs/find_population_specific_haplotypes/4.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { @@ -42,7 +42,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { @@ -67,7 +67,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_specific_variants/1.json b/tests/expected_outputs/find_population_specific_variants/1.json index 99a6fa779..c5f1386bd 100644 --- a/tests/expected_outputs/find_population_specific_variants/1.json +++ b/tests/expected_outputs/find_population_specific_variants/1.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } } ] diff --git a/tests/expected_outputs/find_population_specific_variants/2.json b/tests/expected_outputs/find_population_specific_variants/2.json index 2b9cda4f0..c2a300357 100644 --- a/tests/expected_outputs/find_population_specific_variants/2.json +++ b/tests/expected_outputs/find_population_specific_variants/2.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_structural_intersecting_variants/1.json b/tests/expected_outputs/find_population_structural_intersecting_variants/1.json index a7b5a52ce..dc4769d1c 100644 --- a/tests/expected_outputs/find_population_structural_intersecting_variants/1.json +++ b/tests/expected_outputs/find_population_structural_intersecting_variants/1.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_structural_intersecting_variants/2.json b/tests/expected_outputs/find_population_structural_intersecting_variants/2.json index 215a90aaf..387ea9191 100644 --- a/tests/expected_outputs/find_population_structural_intersecting_variants/2.json +++ b/tests/expected_outputs/find_population_structural_intersecting_variants/2.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } } ] @@ -38,7 +38,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_structural_intersecting_variants/3.json b/tests/expected_outputs/find_population_structural_intersecting_variants/3.json index a7b5a52ce..dc4769d1c 100644 --- a/tests/expected_outputs/find_population_structural_intersecting_variants/3.json +++ b/tests/expected_outputs/find_population_structural_intersecting_variants/3.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_structural_subsuming_variants/1.json b/tests/expected_outputs/find_population_structural_subsuming_variants/1.json index 8d1692f8c..fa120a928 100644 --- a/tests/expected_outputs/find_population_structural_subsuming_variants/1.json +++ b/tests/expected_outputs/find_population_structural_subsuming_variants/1.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } } ] diff --git a/tests/expected_outputs/find_population_structural_subsuming_variants/2.json b/tests/expected_outputs/find_population_structural_subsuming_variants/2.json index 2d61a2459..546f1bed3 100644 --- a/tests/expected_outputs/find_population_structural_subsuming_variants/2.json +++ b/tests/expected_outputs/find_population_structural_subsuming_variants/2.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } } ] @@ -38,7 +38,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_tx_implications/2.json b/tests/expected_outputs/find_population_tx_implications/2.json index c5ea09a1e..bc99cb176 100644 --- a/tests/expected_outputs/find_population_tx_implications/2.json +++ b/tests/expected_outputs/find_population_tx_implications/2.json @@ -13,7 +13,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_tx_implications/3.json b/tests/expected_outputs/find_population_tx_implications/3.json index fd3248467..113c88a65 100644 --- a/tests/expected_outputs/find_population_tx_implications/3.json +++ b/tests/expected_outputs/find_population_tx_implications/3.json @@ -13,7 +13,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_tx_implications/4.json b/tests/expected_outputs/find_population_tx_implications/4.json index fd3e1cd00..e18257812 100644 --- a/tests/expected_outputs/find_population_tx_implications/4.json +++ b/tests/expected_outputs/find_population_tx_implications/4.json @@ -13,7 +13,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { diff --git a/tests/expected_outputs/find_population_tx_implications/5.json b/tests/expected_outputs/find_population_tx_implications/5.json index 40b7b5ecd..664a06555 100644 --- a/tests/expected_outputs/find_population_tx_implications/5.json +++ b/tests/expected_outputs/find_population_tx_implications/5.json @@ -13,7 +13,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1114 + "value": 1115 } }, { From cc91547bd325c01c73a67075f53e7942217a6e4d Mon Sep 17 00:00:00 2001 From: Mihai Todor Date: Thu, 7 Sep 2023 16:50:00 +0100 Subject: [PATCH 2/5] Update python version --- runtime.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/runtime.txt b/runtime.txt index fadb07024..5fa6f4127 100644 --- a/runtime.txt +++ b/runtime.txt @@ -1 +1 @@ -python-3.10.0 +python-3.10.13 From f3411a0e0b960a1454cad7b1d49d9c6677bdd9f5 Mon Sep 17 00:00:00 2001 From: Mihai Todor Date: Thu, 21 Sep 2023 04:37:19 +0100 Subject: [PATCH 3/5] Add SPDI normalization with refseq files --- Procfile | 2 +- app/__init__.py | 4 + app/common.py | 160 +---------------- app/endpoints.py | 52 +++++- app/input_normalization.py | 170 ++++++++++++++++++ app/refseq.py | 28 +++ fetch_refseq.sh | 21 +++ requirements.txt | 3 +- .../test_population_genotype_operations.py | 6 +- utilities/SPDI_Normalization.py | 139 +++++++++++--- 10 files changed, 392 insertions(+), 193 deletions(-) create mode 100644 app/input_normalization.py create mode 100644 app/refseq.py create mode 100755 fetch_refseq.sh diff --git a/Procfile b/Procfile index 40b0792d0..3cd2d20c7 100644 --- a/Procfile +++ b/Procfile @@ -1 +1 @@ -web: gunicorn run:app +web: ./fetch_refseq.sh && gunicorn run:app diff --git a/app/__init__.py b/app/__init__.py index d806ce99e..ea4314107 100644 --- a/app/__init__.py +++ b/app/__init__.py @@ -2,9 +2,13 @@ import flask from flask_cors import CORS import os +# from .refseq import download_refseq_files def create_app(): + # First ensure we have the refseq files locally + # download_refseq_files() + # App and API options = { 'swagger_url': '/', diff --git a/app/common.py b/app/common.py index fca128f5a..c83e6354b 100644 --- a/app/common.py +++ b/app/common.py @@ -2,12 +2,12 @@ from threading import Lock from uuid import uuid4 import pyliftover -import requests from datetime import datetime import pymongo from flask import abort from itertools import groupby import re +from .input_normalization import normalize # MongoDB Client URIs FHIR_genomics_data_client_uri = "mongodb+srv://download:download@cluster0.8ianr.mongodb.net/FHIRGenomicsData" @@ -116,10 +116,6 @@ def get_liftover(from_db, to_db): SUPPORTED_GENOMIC_SOURCE_CLASSES = ['germline', 'somatic'] -NCBI_VARIATION_SERVICES_BASE_URL = 'https://api.ncbi.nlm.nih.gov/variation/v0/' - -CHROMOSOME_CSV_FILE = 'app/_Dict_Chromosome.csv' - # Utility Functions @@ -163,26 +159,6 @@ def merge_ranges(ranges): return merged_ranges -def get_hgvs_contextuals_url(hgvs): - return f"{NCBI_VARIATION_SERVICES_BASE_URL}hgvs/{hgvs}/contextuals" - - -def get_spdi_all_equivalent_contextual_url(contextual_SPDI): - return f'{NCBI_VARIATION_SERVICES_BASE_URL}spdi/{contextual_SPDI}/all_equivalent_contextual' - - -def get_spdi_canonical_representative_url(contextual_SPDI): - return f'{NCBI_VARIATION_SERVICES_BASE_URL}spdi/{contextual_SPDI}/canonical_representative' - - -def build_spdi(seq_id, position, deleted_sequence, inserted_sequence): - return f"{seq_id}:{position}:{deleted_sequence}:{inserted_sequence}" - - -def get_spdi_elements(response_object): - return (response_object['seq_id'], response_object['position'], response_object['deleted_sequence'], response_object['inserted_sequence']) - - def validate_subject(patient_id): if not patients_db.find_one({"patientID": patient_id}): abort(400, f"Patient ({patient_id}) not found.") @@ -196,7 +172,7 @@ def get_variant(variant): variant = variant.lstrip() if variant.count(":") == 1: # HGVS expression - SPDIs = hgvs_2_contextual_SPDIs(variant) + SPDIs = normalize(variant) if not SPDIs: abort(400, f'Cannot normalize variant: {variant}') elif not SPDIs["GRCh37"] and not SPDIs["GRCh38"]: @@ -205,7 +181,7 @@ def get_variant(variant): normalized_variant = {"variant": variant, "GRCh37": SPDIs["GRCh37"], "GRCh38": SPDIs["GRCh38"]} elif variant.count(":") == 3: # SPDI expression - SPDIs = SPDI_2_contextual_SPDIs(variant) + SPDIs = normalize(variant) if not SPDIs: abort(400, f'Cannot normalize variant: {variant}') elif not SPDIs["GRCh37"] and not SPDIs["GRCh38"]: @@ -1002,136 +978,6 @@ def get_intersected_regions(bed_id, build, chrom, start, end, intersected_region intersected_regions.append(f'{ref_seq}:{max(start, csePair["Start"])}-{min(end, csePair["End"])}') -def hgvs_2_contextual_SPDIs(hgvs): - - # convert hgvs to contextualSPDI - url = get_hgvs_contextuals_url(hgvs) - headers = {'Accept': 'application/json'} - - r = requests.get(url, headers=headers) - if r.status_code != 200: - return False - - response = r.json() - raw_data = response['data'] - raw_SPDI = raw_data['spdis'][0] - - seq_id, position, deleted_sequence, inserted_sequence = get_spdi_elements(raw_SPDI) - - contextual_SPDI = build_spdi(seq_id, position, deleted_sequence, inserted_sequence) - - # convert contextualSPDI to build37 and build38 contextual SPDIs - url = get_spdi_all_equivalent_contextual_url(contextual_SPDI) - headers = {'Accept': 'application/json'} - - r = requests.get(url, headers=headers) - if r.status_code != 200: - return False - - response = r.json() - raw_SPDI_List = response['data']['spdis'] - - b37SPDI = None - b38SPDI = None - for item in raw_SPDI_List: - if item['seq_id'].startswith("NC_"): - temp = get_build_and_chrom_by_ref_seq(item['seq_id']) - if temp: - seq_id, position, deleted_sequence, inserted_sequence = get_spdi_elements(item) - - if temp['build'] == 'GRCh37': - b37SPDI = build_spdi(seq_id, position, deleted_sequence, inserted_sequence) - elif temp['build'] == 'GRCh38': - b38SPDI = build_spdi(seq_id, position, deleted_sequence, inserted_sequence) - else: - return False - - return {"GRCh37": b37SPDI, "GRCh38": b38SPDI} - - -def hgvs_2_canonical_SPDI(hgvs): - - # convert hgvs to contextualSPDI - url = get_hgvs_contextuals_url(hgvs) - headers = {'Accept': 'application/json'} - - r = requests.get(url, headers=headers) - if r.status_code != 200: - return False - - response = r.json() - raw_data = response['data'] - raw_SPDI = raw_data['spdis'][0] - - seq_id, position, deleted_sequence, inserted_sequence = get_spdi_elements(raw_SPDI) - - contextual_SPDI = build_spdi(seq_id, position, deleted_sequence, inserted_sequence) - - # convert contextualSPDI to canonical SPDI - url = get_spdi_canonical_representative_url(contextual_SPDI) - headers = {'Accept': 'application/json'} - - r = requests.get(url, headers=headers) - if r.status_code != 200: - return False - - response = r.json() - raw_SPDI = response['data'] - - seq_id, position, deleted_sequence, inserted_sequence = get_spdi_elements(raw_SPDI) - - canonical_SPDI = build_spdi(seq_id, position, deleted_sequence, inserted_sequence) - - return {"canonicalSPDI": canonical_SPDI} - - -def SPDI_2_contextual_SPDIs(spdi): - url = get_spdi_all_equivalent_contextual_url(spdi) - headers = {'Accept': 'application/json'} - - r = requests.get(url, headers=headers) - if r.status_code != 200: - return False - - response = r.json() - raw_SPDI_List = response['data']['spdis'] - - b37SPDI = None - b38SPDI = None - for item in raw_SPDI_List: - if item['seq_id'].startswith("NC_"): - temp = get_build_and_chrom_by_ref_seq(item['seq_id']) - if temp: - seq_id, position, deleted_sequence, inserted_sequence = get_spdi_elements(item) - - if temp['build'] == 'GRCh37': - b37SPDI = build_spdi(seq_id, position, deleted_sequence, inserted_sequence) - elif temp['build'] == 'GRCh38': - b38SPDI = build_spdi(seq_id, position, deleted_sequence, inserted_sequence) - else: - return False - - return {"GRCh37": b37SPDI, "GRCh38": b38SPDI} - - -def SPDI_2_canonical_SPDI(spdi): - url = get_spdi_canonical_representative_url(spdi) - headers = {'Accept': 'application/json'} - - r = requests.get(url, headers=headers) - if r.status_code != 200: - return False - - response = r.json() - raw_SPDI = response['data'] - - seq_id, position, deleted_sequence, inserted_sequence = get_spdi_elements(raw_SPDI) - - canonical_SPDI = build_spdi(seq_id, position, deleted_sequence, inserted_sequence) - - return {"canonicalSPDI": canonical_SPDI} - - def query_clinvar_by_variants(normalized_variant_list, code_list, query, population=False): variant_list = [] for item in normalized_variant_list: diff --git a/app/endpoints.py b/app/endpoints.py index b04c53abf..3c62ba98e 100644 --- a/app/endpoints.py +++ b/app/endpoints.py @@ -179,7 +179,11 @@ def find_subject_specific_variants( subject = subject.strip() common.validate_subject(subject) - variants = list(map(common.get_variant, variants)) + try: + variants = list(map(common.get_variant, variants)) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(422, 'Failed LiftOver') # Query query = {} @@ -833,13 +837,23 @@ def find_subject_tx_implications( if ranges: ranges = list(map(common.get_range, ranges)) common.get_lift_over_range(ranges) - variants = common.get_variants(ranges, query) + + try: + variants = common.get_variants(ranges, query) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(422, 'Failed LiftOver') + if not variants: return jsonify({"resourceType": "Parameters"}) normalized_variants = [{variant["BUILD"]: variant["SPDI"]} for variant in variants] if variants and not ranges: - normalized_variants = list(map(common.get_variant, variants)) + try: + normalized_variants = list(map(common.get_variant, variants)) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(422, 'Failed LiftOver') # Result Object result = OrderedDict() @@ -1105,13 +1119,23 @@ def find_subject_dx_implications( if ranges: ranges = list(map(common.get_range, ranges)) common.get_lift_over_range(ranges) - variants = common.get_variants(ranges, query) + + try: + variants = common.get_variants(ranges, query) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(422, 'Failed LiftOver') + if not variants: return jsonify({"resourceType": "Parameters"}) normalized_variants = [{variant["BUILD"]: variant["SPDI"]} for variant in variants] if variants and not ranges: - normalized_variants = list(map(common.get_variant, variants)) + try: + normalized_variants = list(map(common.get_variant, variants)) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(422, 'Failed LiftOver') # Result Object result = OrderedDict() @@ -1373,7 +1397,11 @@ def find_population_specific_variants( # Parameters variants = list(map(lambda x: x.strip().split(","), variants)) for i in range(len(variants)): - variants[i] = list(map(common.get_variant, variants[i])) + try: + variants[i] = list(map(common.get_variant, variants[i])) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(422, 'Failed LiftOver') # Query query = {} @@ -1884,7 +1912,11 @@ def find_population_tx_implications( return jsonify({"resourceType": "Parameters"}) if variants: - variants = list(map(common.get_variant, variants)) + try: + variants = list(map(common.get_variant, variants)) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(422, 'Failed LiftOver') condition_code_list = [] if conditions: @@ -2090,7 +2122,11 @@ def find_population_dx_implications( return jsonify({"resourceType": "Parameters"}) if variants: - variants = list(map(common.get_variant, variants)) + try: + variants = list(map(common.get_variant, variants)) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(422, 'Failed LiftOver') condition_code_list = [] if conditions: diff --git a/app/input_normalization.py b/app/input_normalization.py new file mode 100644 index 000000000..3221f153c --- /dev/null +++ b/app/input_normalization.py @@ -0,0 +1,170 @@ +import hgvs.parser +import hgvs.dataproviders.uta +import hgvs.assemblymapper +from utilities.SPDI_Normalization import get_normalized_spdi + +hgvsParser = hgvs.parser.Parser() +hgvsDataProvider = hgvs.dataproviders.uta.connect( + db_url="postgresql://anonymous:anonymous@uta.biocommons.org/uta/uta_20210129") +b37hgvsAssemblyMapper = hgvs.assemblymapper.AssemblyMapper( + hgvsDataProvider, assembly_name='GRCh37', alt_aln_method='splign', replace_reference=True) +b38hgvsAssemblyMapper = hgvs.assemblymapper.AssemblyMapper( + hgvsDataProvider, assembly_name='GRCh38', alt_aln_method='splign', replace_reference=True) + +# ------------- point to latest data source ------------------------ +# at unix command line: export UTA_DB_URL=postgresql://anonymous:anonymous@uta.biocommons.org/uta/uta_20210129 + +# ------------------ PARSE ------------- + + +def parse_variant(variant): + parsed_variant_dict = dict() + parsed_variant_dict['parsed'] = hgvsParser.parse_hgvs_variant(variant) + return parsed_variant_dict + +# ------------------ PROJECT ------------- + + +def project_variant(parsed_variant): + projected_variant_dict = dict() + projected_variant_dict['b37projected'] = b37hgvsAssemblyMapper.c_to_g( + parsed_variant) + projected_variant_dict['b38projected'] = b38hgvsAssemblyMapper.c_to_g( + parsed_variant) + return projected_variant_dict + +# ---------------- NORMALIZE to canonical SPDIs --------------- + + +def normalize_variant(parsed_variant, build): + pos = parsed_variant.posedit.pos.start.base-1 + if parsed_variant.posedit.edit.ref: + ref = parsed_variant.posedit.edit.ref + else: # ref is blank for insertions + ref = '' + pos = pos+1 + if str(parsed_variant.posedit.edit) == 'dup': + alt = ref+ref + elif parsed_variant.posedit.edit.alt: + alt = parsed_variant.posedit.edit.alt + else: # alt is blank for deletions + alt = '' + return get_normalized_spdi(parsed_variant.ac, pos, ref, alt, build) + +# ---------------- CONVERT NM_HGVS to canonical SPDIs --------------- + + +def process_NM_HGVS(NM_HGVS): + parsed_variant_dict = parse_variant(NM_HGVS) + print(f"parsed: {parsed_variant_dict['parsed']}") + + projected_variant_dict = project_variant(parsed_variant_dict['parsed']) + print( + f"b37projected: {projected_variant_dict['b37projected']}; b38projected: {projected_variant_dict['b38projected']}") + + b37SPDI = normalize_variant( + projected_variant_dict['b37projected'], 'GRCh37') + b38SPDI = normalize_variant( + projected_variant_dict['b38projected'], 'GRCh38') + print(f"b37normalized: {b37SPDI}; b38normalized: {b38SPDI}") + + return {"GRCh37": b37SPDI, "GRCh38": b38SPDI} + + +# ---------------- CONVERT NC_HGVS to canonical SPDIs --------------- + + +def process_NC_HGVS(NC_HGVS): + parsed_variant_dict = parse_variant(NC_HGVS) + parsed_variant = parsed_variant_dict['parsed'] + print(f"parsed: {parsed_variant_dict['parsed']}") + + transcripts = b38hgvsAssemblyMapper.relevant_transcripts(parsed_variant) + relevantTranscript = next((tr for tr in transcripts if tr.startswith("NM_"))) + var_c = b38hgvsAssemblyMapper.g_to_c( + parsed_variant, relevantTranscript) + + projected_variant_dict = project_variant(var_c) + print( + f"b37projected: {projected_variant_dict['b37projected']}; b38projected: {projected_variant_dict['b38projected']}") + + b37SPDI = normalize_variant( + projected_variant_dict['b37projected'], 'GRCh37') + b38SPDI = normalize_variant( + projected_variant_dict['b38projected'], 'GRCh38') + print(f"b37normalized: {b37SPDI}; b38normalized: {b38SPDI}") + + return {"GRCh37": b37SPDI, "GRCh38": b38SPDI} + +# ---------------- CONVERT NM_SPDI to canonical SPDIs --------------- + + +def process_NM_SPDI(NM_SPDI): + # convert SPDI into NM_HGVS then use NM_HGVS pipeline + refSeq = NM_SPDI.split(":")[0] + pos = int(NM_SPDI.split(":")[1])+1 + ref = NM_SPDI.split(":")[2] + alt = NM_SPDI.split(":")[3] + + if len(ref) == len(alt) == 1: # SNV + var_n = hgvsParser.parse_hgvs_variant( + refSeq+":n."+str(pos)+ref+">"+alt) + elif len(ref) == 0: # INS (e.g. NM_007294.3:c.5533_5534insG) + start = pos-1 + end = start+1 + var_n = hgvsParser.parse_hgvs_variant( + refSeq+":n."+str(start)+"_"+str(end)+'ins'+alt) + elif len(alt) == 0: # DEL (e.g. NM_000527.5:c.1350_1355del) + start = pos + end = start+len(ref)-1 + var_n = hgvsParser.parse_hgvs_variant( + refSeq+":n."+str(start)+"_"+str(end)+'del') + elif len(alt) != 0 and len(ref) != 0: # DELINS (e.g. NM_007294.3:c.5359_5363delinsAGTGA) + start = pos + end = start+len(ref)-1 + var_n = hgvsParser.parse_hgvs_variant( + refSeq+":n."+str(start)+"_"+str(end)+'delins'+alt) + NM_HGVS = b38hgvsAssemblyMapper.n_to_c(var_n) + + return process_NM_HGVS(str(NM_HGVS)) + +# ---------------- CONVERT NC_SPDI to canonical SPDIs --------------- + + +def process_NC_SPDI(NC_SPDI): + # convert SPDI into NC_HGVS then use NC_HGVS pipeline + refSeq = NC_SPDI.split(":")[0] + pos = int(NC_SPDI.split(":")[1])+1 + ref = NC_SPDI.split(":")[2] + alt = NC_SPDI.split(":")[3] + + if len(ref) == len(alt) == 1: # SNV + NC_HGVS = (refSeq+":g."+str(pos)+ref+">"+alt) + elif len(ref) == 0: # INS (e.g. NM_007294.3:c.5533_5534insG) + start = pos-1 + end = start+1 + NC_HGVS = (refSeq+":g."+str(start)+"_"+str(end)+'ins'+alt) + elif len(alt) == 0: # DEL (e.g. NM_000527.5:c.1350_1355del) + start = pos + end = start+len(ref)-1 + NC_HGVS = (refSeq+":g."+str(start)+"_"+str(end)+'del') + elif len(alt) != 0 and len(ref) != 0: # DELINS (e.g. NM_007294.3:c.5359_5363delinsAGTGA) + start = pos + end = start+len(ref)-1 + NC_HGVS = (refSeq+":g."+str(start)+"_"+str(end)+'delins'+alt) + + return process_NC_HGVS(str(NC_HGVS)) + + +def normalize(variant): + print(f"submitted: {variant}") + if variant.upper().startswith('NM'): + if variant.count(':') == 3: + return process_NM_SPDI(variant) + else: + return process_NM_HGVS(variant) + elif variant.upper().startswith('NC'): + if variant.count(':') == 3: + return process_NC_SPDI(variant) + else: + return process_NC_HGVS(variant) diff --git a/app/refseq.py b/app/refseq.py new file mode 100644 index 000000000..50dad1983 --- /dev/null +++ b/app/refseq.py @@ -0,0 +1,28 @@ +from pathlib import Path +from urllib.request import urlretrieve +import tarfile +import sys + +SEQ_FILES_RELEASE_TAG = '113c119' + + +def download_refseq_files(): + try: + # Make sure the parent folder exists + Path('refseq').mkdir(exist_ok=True) + + for build in ['GRCh37', 'GRCh38']: + filename = build + 'seq.tar.gz' + filepath = f'refseq/{filename}' + + if not Path(filepath).is_file(): + # Download seq files + urlretrieve('https://github.com/FHIR/genomics-operations/releases/download/' + SEQ_FILES_RELEASE_TAG + '/' + + filename, filepath) + + # Unarchive seq files each time on startup, just in case something went wrong the previous time + with tarfile.open(filepath, "r|gz") as f: + f.extractall(path='refseq') + + except Exception as err: + sys.exit(f"Failed to download refseq files: {err=}, {type(err)=}") diff --git a/fetch_refseq.sh b/fetch_refseq.sh new file mode 100755 index 000000000..9fc41c78f --- /dev/null +++ b/fetch_refseq.sh @@ -0,0 +1,21 @@ +#!/bin/sh + +if [ -d ./refseq ]; then + echo "refseq files already fetched." + exit 0 +fi + +mkdir -p ./refseq +cd ./refseq + +echo "Downloading refseq files..." + +curl -sLO https://github.com/FHIR/genomics-operations/releases/download/113c119/GRCh37seq.tar.gz +curl -sLO https://github.com/FHIR/genomics-operations/releases/download/113c119/GRCh38seq.tar.gz + +echo "Extracting refseq files..." + +tar -xzf GRCh37seq.tar.gz +tar -xzf GRCh38seq.tar.gz + +echo "Finished extracting refseq files." diff --git a/requirements.txt b/requirements.txt index 776b6ae64..04821d920 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,8 +3,9 @@ dnspython==2.2.1 Flask==2.2.2 flask_cors==4.0.0 gunicorn==20.1.0 +hgvs==1.5.4 pandas==1.3.5 -pyfastx==0.9.1 +pyfastx==2.0.1 pyliftover==0.4 pymongo==4.2.0 pyranges==0.0.120 diff --git a/tests/integration_tests/test_population_genotype_operations.py b/tests/integration_tests/test_population_genotype_operations.py index c3327d55b..d2f5e2343 100644 --- a/tests/integration_tests/test_population_genotype_operations.py +++ b/tests/integration_tests/test_population_genotype_operations.py @@ -22,11 +22,15 @@ def test_population_specific_variants_1(client): def test_population_specific_variants_2(client): + # PDE4DIP gene is on minus strand in build 37, but is on positive strand in build 38. + # The hgvs library returns "HGVSDataNotAvailableError: No alignments for NM_001002811.1 in GRCh38 using splign" in + # this case. + url = tu.find_population_specific_variants_query( 'subject=HG00403&variants=NC_000001.10:144931726:G:A&variants=NC_000001.10:145532548:T:C&variants=NC_000001.10:145592072:A:T&includePatientList=true') response = client.get(url) - tu.compare_actual_and_expected_output(f'{tu.FIND_POPULATION_SPECIFIC_VARIANTS_OUTPUT_DIR}2.json', response.json) + assert response.status_code == 422 # def test_population_specific_variants_3(client): diff --git a/utilities/SPDI_Normalization.py b/utilities/SPDI_Normalization.py index 1424281d9..035befc23 100644 --- a/utilities/SPDI_Normalization.py +++ b/utilities/SPDI_Normalization.py @@ -1,34 +1,123 @@ -import pyfastx +from glob import glob +from math import ceil +from pathlib import Path from threading import Lock +MAX_REFSEQ_CACHE_SIZE = 300 * 1024 * 1024 # 450MB + + +def hex2Code(code): + """ + Convert a 4 byte integer to a sequence code. See `code2Hex()` for details + """ + match code: + case 0x0: return 'A' + case 0x1: return 'C' + case 0x2: return 'G' + case 0x3: return 'T' + case 0x4: return 'U' + case 0x5: return 'R' + case 0x6: return 'Y' + case 0x7: return 'K' + case 0x8: return 'M' + case 0x9: return 'S' + case 0xA: return 'W' + case 0xB: return 'B' + case 0xC: return 'D' + case 0xD: return 'H' + case 0xE: return 'V' + case 0xF: return 'N' + case _: + raise NotImplementedError(f"unsupported code '{code}'") + + +class RefSeq: + def __init__(self, refSeq, length): + """ + Store the length so we can easily tell when the index operator receives an out of bounds subscript because we + want to represent each symbol using 4 bits and we have 16 symbols, so there's no room left to represent an + invalid symbol which denotes end-of-sequence. + """ + self.__packedRefSeq = refSeq + self.__length = length + + # Size in bytes + self.size = ceil(length / 2) + + def __getitem__(self, subscript): + if isinstance(subscript, slice): + codes = [] + for i in range(*subscript.indices(self.__length)): + codes.append(self[i]) + return "".join(codes) + + byte = self.__packedRefSeq[subscript // 2] + + if subscript % 2 == 0: + return hex2Code((byte >> 4) & 0xF) + + return hex2Code(byte & 0xF) + + +class RefSeqCache: + def __init__(self): + self.cache = {} + self.size = 0 + + def __getKey(self, build, refSeqName): + return f'{build}/{refSeqName}' + + def contains(self, build, refSeqName): + return self.__getKey(build, refSeqName) in self.cache + + def get(self, build, refSeqName): + return self.cache[self.__getKey(build, refSeqName)] + + def add(self, build, refSeqName, refSeq): + """ + Evacuates enough items from the cache until we have room for the new one and then adds it to the cache. + """ + + while self.size + refSeq.size > MAX_REFSEQ_CACHE_SIZE: + key, cachedRefSeq = next(iter(self.cache.items())) + self.size -= cachedRefSeq.size + del self.cache[key] + + self.cache[self.__getKey(build, refSeqName)] = refSeq + self.size += refSeq.size + + +refSeqCache = RefSeqCache() +refSeqCacheLock = Lock() + + +def get_ref_seq(build, refSeqName): + if not refSeqCache.contains(build, refSeqName): + try: + refSeqFilePattern = f'refseq/{build}/{refSeqName}_*.seq' + foundFiles = glob(refSeqFilePattern) + if len(foundFiles) != 1: + raise ValueError(f'failed to find expected refseq file for build "{build}" and RefSeq "{refSeqName}" ') + + refSeqFile = foundFiles[0] + refSeqLength = Path(refSeqFile).stem.rpartition('_')[-1] + with open(refSeqFile, mode='rb') as file: + refSeq = RefSeq(file.read(), int(refSeqLength)) + except Exception as err: + print(f"failed to read refseq file: {err=}, {type(err)=}") + raise + refSeqCache.add(build, refSeqName, refSeq) + return refSeqCache.get(build, refSeqName) -# Fasta file handles cache -fasta_cache = {} -fasta_lock = Lock() -BUILD37_FILE = 'FASTA/GRCh37_latest_genomic.fna.gz' -BUILD38_FILE = 'FASTA/GRCh38_latest_genomic.fna.gz' - - -def get_fasta(file): - with fasta_lock: - if file not in fasta_cache: - try: - fasta = pyfastx.Fasta(file) - except Exception as err: - print(f"Unexpected {err=}, {type(err)=}") - raise - fasta_cache[file] = fasta - return fasta_cache[file] +def get_normalized_spdi(ref_seq, pos, ref, alt, build): + # Need to serialise this if we can't keep all the ref seq data in memory + with refSeqCacheLock: + return get_normalized_spdi_impl(ref_seq, pos, ref, alt, build) -def get_normalized_spdi(ref_seq, pos, ref, alt, build): - if build == 'GRCh37': - GRCh37_ref_seq_fasta = get_fasta(BUILD37_FILE) - ref_seq_fasta = GRCh37_ref_seq_fasta[ref_seq] - elif build == 'GRCh38': - GRCh38_ref_seq_fasta = get_fasta(BUILD38_FILE) - ref_seq_fasta = GRCh38_ref_seq_fasta[ref_seq] +def get_normalized_spdi_impl(ref_seq, pos, ref, alt, build): + ref_seq_fasta = get_ref_seq(build, ref_seq) # Step 0 start = pos From 5c21e72ce014dfe891544ebdacb63bdf327296f9 Mon Sep 17 00:00:00 2001 From: Mihai Todor Date: Wed, 27 Sep 2023 03:19:33 +0100 Subject: [PATCH 4/5] WIP --- utilities/SPDI_Normalization.py | 57 +++++++++++---------------------- 1 file changed, 19 insertions(+), 38 deletions(-) diff --git a/utilities/SPDI_Normalization.py b/utilities/SPDI_Normalization.py index 035befc23..3969270ab 100644 --- a/utilities/SPDI_Normalization.py +++ b/utilities/SPDI_Normalization.py @@ -2,11 +2,12 @@ from math import ceil from pathlib import Path from threading import Lock +from memory_profiler import profile -MAX_REFSEQ_CACHE_SIZE = 300 * 1024 * 1024 # 450MB +MAX_REFSEQ_CACHE_SIZE = 200 * 1024 * 1024 # 450MB -def hex2Code(code): +def hex_to_code(code): """ Convert a 4 byte integer to a sequence code. See `code2Hex()` for details """ @@ -31,6 +32,11 @@ def hex2Code(code): raise NotImplementedError(f"unsupported code '{code}'") +def refseq_length_to_size(length): + # Each byte contains at most 2 codes + return ceil(length / 2) + + class RefSeq: def __init__(self, refSeq, length): """ @@ -42,7 +48,7 @@ def __init__(self, refSeq, length): self.__length = length # Size in bytes - self.size = ceil(length / 2) + self.size = refseq_length_to_size(length) def __getitem__(self, subscript): if isinstance(subscript, slice): @@ -54,41 +60,12 @@ def __getitem__(self, subscript): byte = self.__packedRefSeq[subscript // 2] if subscript % 2 == 0: - return hex2Code((byte >> 4) & 0xF) - - return hex2Code(byte & 0xF) - - -class RefSeqCache: - def __init__(self): - self.cache = {} - self.size = 0 - - def __getKey(self, build, refSeqName): - return f'{build}/{refSeqName}' - - def contains(self, build, refSeqName): - return self.__getKey(build, refSeqName) in self.cache + return hex_to_code((byte >> 4) & 0xF) - def get(self, build, refSeqName): - return self.cache[self.__getKey(build, refSeqName)] + return hex_to_code(byte & 0xF) - def add(self, build, refSeqName, refSeq): - """ - Evacuates enough items from the cache until we have room for the new one and then adds it to the cache. - """ - - while self.size + refSeq.size > MAX_REFSEQ_CACHE_SIZE: - key, cachedRefSeq = next(iter(self.cache.items())) - self.size -= cachedRefSeq.size - del self.cache[key] - - self.cache[self.__getKey(build, refSeqName)] = refSeq - self.size += refSeq.size - -refSeqCache = RefSeqCache() -refSeqCacheLock = Lock() +refSeqLock = Lock() def get_ref_seq(build, refSeqName): @@ -100,9 +77,13 @@ def get_ref_seq(build, refSeqName): raise ValueError(f'failed to find expected refseq file for build "{build}" and RefSeq "{refSeqName}" ') refSeqFile = foundFiles[0] - refSeqLength = Path(refSeqFile).stem.rpartition('_')[-1] + refSeqLength = int(Path(refSeqFile).stem.rpartition('_')[-1]) + + # Make room for this refseq in the cache + refSeqCache.reserve_space(refSeqLength) + with open(refSeqFile, mode='rb') as file: - refSeq = RefSeq(file.read(), int(refSeqLength)) + refSeq = RefSeq(file.read(), refSeqLength) except Exception as err: print(f"failed to read refseq file: {err=}, {type(err)=}") raise @@ -112,7 +93,7 @@ def get_ref_seq(build, refSeqName): def get_normalized_spdi(ref_seq, pos, ref, alt, build): # Need to serialise this if we can't keep all the ref seq data in memory - with refSeqCacheLock: + with refSeqLock: return get_normalized_spdi_impl(ref_seq, pos, ref, alt, build) From 3f80058802cb75d53efbad99d818bbb575dc586c Mon Sep 17 00:00:00 2001 From: Mihai Todor Date: Thu, 5 Oct 2023 17:14:43 +0300 Subject: [PATCH 5/5] Test SPDI 3 bit packing RefSeq scheme --- .github/workflows/main.yml | 12 ++-- README.md | 6 ++ app/__init__.py | 4 -- app/api_spec.yml | 55 +++++++++++--- app/refseq.py | 28 -------- app/utilities_endpoints.py | 12 +++- fetch_refseq.sh | 8 +-- utilities/SPDI_Normalization.py | 124 ++++++++++++++++---------------- 8 files changed, 137 insertions(+), 112 deletions(-) delete mode 100644 app/refseq.py diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 318b9bafe..580e6cfd7 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -28,20 +28,24 @@ jobs: args: --extend-ignore E501,E741 - name: Run Tests - run: python -m pytest + # run: ./fetch_refseq.sh python -m pytest + run: "echo ${heroku_app_name} && false" + with: + heroku_app_name: ${{secrets.HEROKU_DEV_APP_NAME}} + # TODO: Add a way to deploy to Prod manually deploy: name: Deploy runs-on: ubuntu-latest - if: ${{ github.event_name == 'push' && github.ref == 'refs/heads/main' && always() && contains(join(needs.*.result, ','), 'success') }} + if: ${{ contains(join(needs.*.result, ','), 'success') }} needs: [test] steps: - uses: actions/checkout@v2 - - uses: akhileshns/heroku-deploy@v3.12.12 + - uses: akhileshns/heroku-deploy@v3.12.14 with: heroku_api_key: ${{secrets.HEROKU_API_KEY}} - heroku_app_name: ${{secrets.HEROKU_APP_NAME}} + heroku_app_name: ${{secrets.HEROKU_DEV_APP_NAME}} heroku_email: ${{secrets.HEROKU_EMAIL}} diff --git a/README.md b/README.md index 0ef80cf03..f4785b1e3 100644 --- a/README.md +++ b/README.md @@ -48,3 +48,9 @@ run `python3 -m pytest` from the terminal to execute them all. Additionally, since the tests run against the Mongo DB database, if you need to update the test data in this repo, you can run `OVERWRITE_TEST_EXPECTED_DATA=true python3 -m pytest` from the terminal and then create a pull request with the changes. + +## Development environment on Heroku + +Pull requests will trigger a deployment to this environment automatically which is accessible at the following URL: + +https://fhir-gen-ops-dev-ca42373833b6.herokuapp.com/ diff --git a/app/__init__.py b/app/__init__.py index ea4314107..d806ce99e 100644 --- a/app/__init__.py +++ b/app/__init__.py @@ -2,13 +2,9 @@ import flask from flask_cors import CORS import os -# from .refseq import download_refseq_files def create_app(): - # First ensure we have the refseq files locally - # download_refseq_files() - # App and API options = { 'swagger_url': '/', diff --git a/app/api_spec.yml b/app/api_spec.yml index 73ba7c996..23403ce7f 100644 --- a/app/api_spec.yml +++ b/app/api_spec.yml @@ -102,7 +102,7 @@ paths: type: boolean default: false description: Include sequence phase relationships in response if set to true. - + /subject-operations/genotype-operations/$find-subject-specific-variants: get: description: > @@ -177,7 +177,7 @@ paths: - "germline" - "somatic" description: Enables an App to limit results to those that are 'germline' or 'somatic'. Default is to include variants irrespective of genomic source class. - + /subject-operations/genotype-operations/$find-subject-structural-intersecting-variants: get: description: > @@ -262,7 +262,7 @@ paths: type: boolean default: false description: Include variants in response if set to true. - + /subject-operations/genotype-operations/$find-subject-structural-subsuming-variants: get: description: > @@ -346,7 +346,7 @@ paths: type: boolean default: false description: Include variants in response if set to true. - + /subject-operations/genotype-operations/$find-subject-haplotypes: get: description: > @@ -422,7 +422,7 @@ paths: - "germline" - "somatic" description: Enables an App to limit results to those that are 'germline' or 'somatic'. Default is to include variants irrespective of genomic source class. - + /subject-operations/genotype-operations/$find-subject-specific-haplotypes: get: description: > @@ -497,7 +497,7 @@ paths: - "germline" - "somatic" description: Enables an App to limit results to those that are 'germline' or 'somatic'. Default is to include variants irrespective of genomic source class. - + /subject-operations/phenotype-operations/$find-subject-tx-implications: get: description: |- @@ -614,7 +614,7 @@ paths: - "germline" - "somatic" description: Enables an App to limit results to those that are 'germline' or 'somatic'. Default is to include variants irrespective of genomic source class. - + /subject-operations/phenotype-operations/$find-subject-dx-implications: get: description: |- @@ -713,7 +713,7 @@ paths: - "germline" - "somatic" description: Enables an App to limit results to those that are 'germline' or 'somatic'. Default is to include variants irrespective of genomic source class. - + /subject-operations/metadata-operations/$find-study-metadata: get: description: |- @@ -1147,6 +1147,7 @@ paths: type: string pattern: '^\s*[Nn][Pp]_\d{4,10}(\.)?(\d{1,2})?\s*$' example: "NP_000005.3" + /utilities/find-the-gene: get: summary: "Find The Gene" @@ -1170,6 +1171,42 @@ paths: pattern: '^\s*[Nn][Cc]_\d{4,10}(\.)(\d{1,2}):\d{1,10}-\d{1,10}\s*$' example: "NC_000001.11:11794399-11794400" + /utilities/seqfetcher/1/sequence/{ref_seq}: + get: + summary: "Seqfetcher" + operationId: "app.utilities_endpoints.seqfetcher" + tags: + - "Seqfetcher Utility" + responses: + '200': + description: "Returns RefSeq subsequence" + content: + text/plain: + schema: + type: string + parameters: + - name: ref_seq + in: path + required: true + description: RefSeq + schema: + type: string + example: "NC_000001.10" + - name: start + in: query + required: true + description: Subsequence start index + schema: + type: integer + example: 1 + - name: end + in: query + required: true + description: Subsequence end index + schema: + type: integer + example: 2 + tags: - name: Subject Genotype Operations - name: Subject Phenotype Operations @@ -1178,7 +1215,7 @@ tags: - name: Population Phenotype Operations - name: Feature Coordinates Utility description: This utility returns genomic feature coordinates and other annotations. All data are from NCBI Human Genome Resources. For chromosomes, build 37 and build 38 reference sequences are returned. For genes, genomic coordinates are returned, along with a list of transcripts. MANE transcript is flagged. For transcripts, genomic coordinates are returned, along with the gene name and composite exons, along with exon coordinates. For proteins, the corresponding transcript is returned. - + - name: Find The Gene Utility description: This utility returns all genes that intersect with a provided genomic region. diff --git a/app/refseq.py b/app/refseq.py deleted file mode 100644 index 50dad1983..000000000 --- a/app/refseq.py +++ /dev/null @@ -1,28 +0,0 @@ -from pathlib import Path -from urllib.request import urlretrieve -import tarfile -import sys - -SEQ_FILES_RELEASE_TAG = '113c119' - - -def download_refseq_files(): - try: - # Make sure the parent folder exists - Path('refseq').mkdir(exist_ok=True) - - for build in ['GRCh37', 'GRCh38']: - filename = build + 'seq.tar.gz' - filepath = f'refseq/{filename}' - - if not Path(filepath).is_file(): - # Download seq files - urlretrieve('https://github.com/FHIR/genomics-operations/releases/download/' + SEQ_FILES_RELEASE_TAG + '/' + - filename, filepath) - - # Unarchive seq files each time on startup, just in case something went wrong the previous time - with tarfile.open(filepath, "r|gz") as f: - f.extractall(path='refseq') - - except Exception as err: - sys.exit(f"Failed to download refseq files: {err=}, {type(err)=}") diff --git a/app/utilities_endpoints.py b/app/utilities_endpoints.py index c2f278ff3..03fdc0857 100644 --- a/app/utilities_endpoints.py +++ b/app/utilities_endpoints.py @@ -1,6 +1,7 @@ from flask import abort, jsonify from collections import OrderedDict from app import common +from utilities import SPDI_Normalization def get_feature_coordinates( @@ -133,7 +134,8 @@ def get_feature_coordinates( protein = protein.split('.')[0] try: - result = common.proteins_data.aggregate([{"$match": {"proteinRefSeq": {'$regex': ".*"+str(protein).replace('*', r'\*')+".*"}}}]) + result = common.proteins_data.aggregate( + [{"$match": {"proteinRefSeq": {'$regex': ".*"+str(protein).replace('*', r'\*')+".*"}}}]) result = list(result) except Exception as e: print(f"DEBUG: Error({e}) under get_feature_coordinates(protein={protein})") @@ -189,3 +191,11 @@ def find_the_gene(range=None): output.append(ord_dict) return (jsonify(output)) + + +def seqfetcher(ref_seq, start, end): + try: + subseq = SPDI_Normalization.get_ref_seq_subseq('GRCh37', ref_seq, start, end) + except Exception: + subseq = SPDI_Normalization.get_ref_seq_subseq('GRCh38', ref_seq, start, end) + return f'>{ref_seq}:{start}-{end} Homo sapiens chromosome 1, GRCh37.p13 Primary Assembly\n{subseq}\n\n' diff --git a/fetch_refseq.sh b/fetch_refseq.sh index 9fc41c78f..f409c6c49 100755 --- a/fetch_refseq.sh +++ b/fetch_refseq.sh @@ -10,12 +10,12 @@ cd ./refseq echo "Downloading refseq files..." -curl -sLO https://github.com/FHIR/genomics-operations/releases/download/113c119/GRCh37seq.tar.gz -curl -sLO https://github.com/FHIR/genomics-operations/releases/download/113c119/GRCh38seq.tar.gz +curl -sLO https://github.com/FHIR/genomics-operations/releases/download/113c119/GRCh37_refseq.tar.gz +curl -sLO https://github.com/FHIR/genomics-operations/releases/download/113c119/GRCh38_refseq.tar.gz echo "Extracting refseq files..." -tar -xzf GRCh37seq.tar.gz -tar -xzf GRCh38seq.tar.gz +tar -xzf GRCh37_refseq.tar.gz +tar -xzf GRCh38_refseq.tar.gz echo "Finished extracting refseq files." diff --git a/utilities/SPDI_Normalization.py b/utilities/SPDI_Normalization.py index 3969270ab..dc2dd890d 100644 --- a/utilities/SPDI_Normalization.py +++ b/utilities/SPDI_Normalization.py @@ -1,99 +1,99 @@ from glob import glob -from math import ceil +from math import floor from pathlib import Path from threading import Lock -from memory_profiler import profile -MAX_REFSEQ_CACHE_SIZE = 200 * 1024 * 1024 # 450MB - -def hex_to_code(code): +def hex_to_code(hex): """ - Convert a 4 byte integer to a sequence code. See `code2Hex()` for details + Convert a 3 bit integer to a sequence code. See `code_to_hex()` in pack_refseq.py for details. """ - match code: + match hex: case 0x0: return 'A' case 0x1: return 'C' case 0x2: return 'G' case 0x3: return 'T' - case 0x4: return 'U' - case 0x5: return 'R' - case 0x6: return 'Y' - case 0x7: return 'K' - case 0x8: return 'M' - case 0x9: return 'S' - case 0xA: return 'W' - case 0xB: return 'B' - case 0xC: return 'D' - case 0xD: return 'H' - case 0xE: return 'V' - case 0xF: return 'N' + case 0x4: return 'N' + # Blow up if we get any unexpected hex-encoded code case _: - raise NotImplementedError(f"unsupported code '{code}'") - - -def refseq_length_to_size(length): - # Each byte contains at most 2 codes - return ceil(length / 2) + raise NotImplementedError(f"unexpected hex-encoded code: '{hex}'") class RefSeq: - def __init__(self, refSeq, length): - """ - Store the length so we can easily tell when the index operator receives an out of bounds subscript because we - want to represent each symbol using 4 bits and we have 16 symbols, so there's no room left to represent an - invalid symbol which denotes end-of-sequence. - """ - self.__packedRefSeq = refSeq - self.__length = length - - # Size in bytes - self.size = refseq_length_to_size(length) + # TODO: Consider adding `__str__()` method (and maybe `__repr__()` too?) + + def __init__(self, packed_ref_seq, len): + self.__packed_ref_seq = packed_ref_seq + # Store the RefSeq length so we can easily tell when the index operator receives an out of bounds subscript + # because we want to represent each code using 3 bits and we don't want to introduce an end-of-sequence code, in + # case we'll want to leverage the unused codes for something else. + self.__len = len def __getitem__(self, subscript): if isinstance(subscript, slice): codes = [] - for i in range(*subscript.indices(self.__length)): + # TODO: It should be more efficient to process 3 bytes at a time for long ranges. + for i in range(*subscript.indices(self.__len)): codes.append(self[i]) return "".join(codes) - byte = self.__packedRefSeq[subscript // 2] + if subscript >= self.__len: + raise IndexError(f"list index out of range: '{subscript}' must be less than '{self.__len}'") + + # Select the byte(s) and the offset of the code we wish to extract from the packed RefSeq data. + packed_subscript = floor(subscript * 3 / 8) + data = self.__packed_ref_seq[packed_subscript] + match subscript % 8: + case 0: offset = 0 + case 1: offset = 3 + case 2: + offset = 6 + # This code is packed across two bytes. + data += self.__packed_ref_seq[packed_subscript + 1] << 8 + case 3: offset = 1 + case 4: offset = 4 + case 5: + offset = 7 + # This code is packed across two bytes. + data += self.__packed_ref_seq[packed_subscript + 1] << 8 + case 6: offset = 2 + case 7: offset = 5 + + # Move the relevant three bits to the bottom, extract them using a bit mask and convert them back to a code. + return hex_to_code((data >> offset) & 0x3) - if subscript % 2 == 0: - return hex_to_code((byte >> 4) & 0xF) - return hex_to_code(byte & 0xF) +def get_ref_seq(build, ref_seq_name): + try: + ref_seq_file_pattern = f'refseq/{build}/{ref_seq_name}_*.refseq' + found_files = glob(ref_seq_file_pattern) + if len(found_files) != 1: + raise ValueError(f'failed to find expected refseq file for build "{build}" and RefSeq "{ref_seq_name}" ') + ref_seq_file = found_files[0] + ref_seq_length = int(Path(ref_seq_file).stem.rpartition('_')[-1]) -refSeqLock = Lock() + with open(ref_seq_file, mode='rb') as file: + refSeq = RefSeq(file.read(), ref_seq_length) + except Exception as err: + print(f"failed to read refseq file: {err=}, {type(err)=}") + raise + return refSeq -def get_ref_seq(build, refSeqName): - if not refSeqCache.contains(build, refSeqName): - try: - refSeqFilePattern = f'refseq/{build}/{refSeqName}_*.seq' - foundFiles = glob(refSeqFilePattern) - if len(foundFiles) != 1: - raise ValueError(f'failed to find expected refseq file for build "{build}" and RefSeq "{refSeqName}" ') - refSeqFile = foundFiles[0] - refSeqLength = int(Path(refSeqFile).stem.rpartition('_')[-1]) +ref_seq_lock = Lock() - # Make room for this refseq in the cache - refSeqCache.reserve_space(refSeqLength) - with open(refSeqFile, mode='rb') as file: - refSeq = RefSeq(file.read(), refSeqLength) - except Exception as err: - print(f"failed to read refseq file: {err=}, {type(err)=}") - raise - refSeqCache.add(build, refSeqName, refSeq) - return refSeqCache.get(build, refSeqName) +def get_ref_seq_subseq(build, ref_seq, start, end): + # Need to serialise this if we can't keep all the RefSeq data in memory + with ref_seq_lock: + return get_ref_seq(build, ref_seq)[start:end] def get_normalized_spdi(ref_seq, pos, ref, alt, build): - # Need to serialise this if we can't keep all the ref seq data in memory - with refSeqLock: + # Need to serialise this if we can't keep all the RefSeq data in memory + with ref_seq_lock: return get_normalized_spdi_impl(ref_seq, pos, ref, alt, build)