From faa468de747424f2f8648b27803f57a1e7f850c5 Mon Sep 17 00:00:00 2001 From: Mihai Todor Date: Thu, 21 Sep 2023 04:37:19 +0100 Subject: [PATCH 1/2] Add SPDI normalization with refseq files --- .env | 1 + .github/workflows/cicd.yml | 1 + .gitignore | 6 +- .vscode/settings.json | 3 - README.md | 155 +++++++++++++- app/__init__.py | 7 + app/api_spec.yml | 58 +++++ app/common.py | 179 ++-------------- app/endpoints.py | 52 ++++- app/input_normalization.py | 201 ++++++++++++++++++ app/utilities_endpoints.py | 41 ++-- fetch_utilities_data.sh | 19 ++ requirements.txt | 8 + tests/conftest.py | 21 +- .../find_population_specific_variants/3.json | 6 +- .../test_normalize_variant_utility.py | 26 +++ .../test_population_genotype_operations.py | 17 +- tests/unit_tests/test_seqrepo_data_packing.py | 46 ++++ tests/utilities.py | 7 + utilities/SPDI_Normalization.py | 163 -------------- utilities/common.py | 5 +- utilities/gene_ref_seq.py | 2 +- utilities/{pyard.py => hla.py} | 6 +- utilities/pack_utilities_data.py | 187 ++++++++++++++++ utilities/requirements.txt | 10 + utilities/spdi.py | 12 ++ utilities/spdi_refseq.py | 109 ++++++++++ utilities/vcf2json.py | 43 +++- 28 files changed, 1005 insertions(+), 386 deletions(-) create mode 100644 app/input_normalization.py create mode 100644 tests/integration_tests/test_normalize_variant_utility.py create mode 100644 tests/unit_tests/test_seqrepo_data_packing.py delete mode 100644 utilities/SPDI_Normalization.py rename utilities/{pyard.py => hla.py} (64%) create mode 100644 utilities/pack_utilities_data.py create mode 100644 utilities/requirements.txt create mode 100644 utilities/spdi.py create mode 100644 utilities/spdi_refseq.py diff --git a/.env b/.env index c0ed5abaa..c7c606c6f 100644 --- a/.env +++ b/.env @@ -1,2 +1,3 @@ UTILITIES_DATA_VERSION=113c119 +UTA_DATABASE_SCHEMA=uta_20240523b PYARD_DATABASE_VERSION=3580 diff --git a/.github/workflows/cicd.yml b/.github/workflows/cicd.yml index ad6fbf3c6..4ac95ba11 100644 --- a/.github/workflows/cicd.yml +++ b/.github/workflows/cicd.yml @@ -33,6 +33,7 @@ jobs: run: ./fetch_utilities_data.sh && python -m pytest env: MONGODB_READONLY_PASSWORD: ${{ secrets.MONGODB_READONLY_PASSWORD }} + UTA_DATABASE_URL: ${{ secrets.UTA_DATABASE_URL }} deploy: name: Deploy to dev diff --git a/.gitignore b/.gitignore index a20a98cb0..cad092af9 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,8 @@ .pytest_cache __pycache__ .venv -utilities/FASTA -/data secrets.env +/data +/seqrepo +/tmp +/utilities/FASTA diff --git a/.vscode/settings.json b/.vscode/settings.json index 451c92c0d..61327f206 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -10,9 +10,6 @@ "python.testing.pytestArgs": [ "." ], - "[python]": { - "editor.defaultFormatter": "ms-python.autopep8", - }, "autopep8.args": [ "--max-line-length=200" ], diff --git a/README.md b/README.md index 75964c677..07bbc74ee 100644 --- a/README.md +++ b/README.md @@ -42,17 +42,156 @@ The operations return the following status codes: ## Testing -To run the [integration tests](https://github.com/FHIR/genomics-operations/tree/main/tests), you can use the VS Code Testing functionality which should discover them automatically. You can also -run `python3 -m pytest` from the terminal to execute them all. +For local development, you will have to create a `secrets.env` file in the root of the repo and add in it the MongoDB +password and the UTA Postgres database connection string (see the UTA section below for details): + +``` +MONGODB_READONLY_PASSWORD=... +UTA_DATABASE_URL=... +``` + +Then, you will need to run `fetch_utilities_data.sh` in a terminal to fetch the required data files: + +```shell +$ ./fetch_utilities_data.sh +``` + +To run the [integration tests](https://github.com/FHIR/genomics-operations/tree/main/tests), you can use the VS Code +Testing functionality which should discover them automatically. You can also 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. -## Update py-ard database +## Heroku Deployment + +Currently, there are two environments running in Heroku: +- Dev: +- Prod: + +Pull requests will trigger a deployment to the dev environment automatically after being merged. + +The ["Manual Deployment"](https://github.com/FHIR/genomics-operations/actions/workflows/manual_deployment.yml) workflow +can be used to deploy code to either the `dev` or `prod` environments. To do so, please select "Run workflow", ignore +the "Use workflow from" dropdown which lists the branches in the current repo (I can't disable / remove it) and then +select the environment, the branch and the repository. By default, the `https://github.com/FHIR/genomics-operations` +repo is specified, but you can replace it with any any fork. + +Deployments to the prod environment can only be triggered manually from the `main` branch of the repo using the Manual +Deployment. + +### Heroku Stack + +Make sure that the Python version under [`runtime.txt`](./runtime.txt) is +[supported](https://devcenter.heroku.com/articles/python-support#supported-runtimes) by the +[Heroku stack](https://devcenter.heroku.com/articles/stack) that is currently running in each environment. + +## UTA Database + +The Biocommons [hgvs](https://github.com/biocommons/hgvs) library which is used for variant parsing, validation and +normalisation requires access to a copy of the [UTA](https://github.com/biocommons/uta) Postgres database. + +We have provisioned a Heroku Postgres instance in the Prod environment which contains the imported data from a database +dump as described [here](https://github.com/biocommons/uta#installing-from-database-dumps). + +We define a `UTA_DATABASE_SCHEMA` environment variable in the [`.env`](.env) file which contains the name of the +currently imported database schema. + +### Database import procedure (it will take about 30 minutes): + +- Go to the UTA dump download site (http://dl.biocommons.org/uta/) and get the latest `.pgd.gz` file. +- Go to https://dashboard.heroku.com/apps/fhir-gen-ops/resources and click on the "Heroku Postgres" instance (it will +open a new window) +- Go to the Settings tab +- Click "View Credentials" +- Use the fields from this window to fill in the variables below + +```shell +$ POSTGRES_HOST="" +$ POSTGRES_DATABASE="" +$ POSTGRES_USER="" +$ PGPASSWORD="" +$ UTA_SCHEMA="" # Specify the UTA schema of the UTA dump you downloaded (example: uta_20240523b) +$ gzip -cdq ${UTA_SCHEMA}.pgd.gz | grep -v '^GRANT USAGE ON SCHEMA .* TO anonymous;$' | grep -v '^ALTER .* OWNER TO uta_admin;$' | psql -U ${POSTGRES_USER} -1 -v ON_ERROR_STOP=1 -d ${POSTGRES_DATABASE} -h ${POSTGRES_HOST} -Eae +``` + +Note: The `grep -v` commands are required because the Heroku Postgres instance doesn't allow us to create a new role. + +Once complete, make sure you update the `UTA_DATABASE_SCHEMA` environment variable in the [`.env`](.env) file and commit +it. + +### Connection string + +The connection string for this database can be found in the same Heroku Postgres Settings tab under "View Credentials". +It is pre-populated in the Heroku runtime under the `UTA_DATABASE_URL` environment variable. Additionally, we set the +same `UTA_DATABASE_URL` environment variable in GitHub so the CI can can use this database when running the tests. + +For local development, set `UTA_DATABASE_URL` to the Heroku Postgres connection string in the `secrets.env` file. +Alternatively, you can set it to `postgresql://anonymous:anonymous@uta.biocommons.org/uta` if you'd like to use the HGVS +public instance. + +### Testing the database + +```shell +$ source secrets.env +$ pgcli "${UTA_DATABASE_URL}" +> set schema ''; # Specify the UTA schema of the UTA dump you downloaded (example: uta_20240523b) +> select count(*) from alembic_version + union select count(*) from associated_accessions + union select count(*) from exon + union select count(*) from exon_aln + union select count(*) from exon_set + union select count(*) from gene + union select count(*) from meta + union select count(*) from origin + union select count(*) from seq + union select count(*) from seq_anno + union select count(*) from transcript + union select count(*) from translation_exception; +``` + +### Update utilities data + +The RefSeq metadata from the UTA database needs to be in sync with the RefSeq data which is available for the Seqfetcher +Utility endpoint. Currently, this is stored in GitHub as release artifacts. Similarly, the PyARD SQLite database is also +stored as a release artifact. + +To update the RefSeq data and PyARD database, you will have to run `./utilities/pack_seqrepo_data.py`. Here is a +step-by-step guide on how to do this: + +```shell +$ mkdir seqrepo +$ cd seqrepo +$ python3 -m venv .venv +$ . .venv/bin/activate +$ pip install setuptools==75.7.0 +$ pip install biocommons.seqrepo==0.6.9 +$ # See https://github.com/biocommons/biocommons.seqrepo/issues/171 for a bug that's causing issues with the builtin +$ # rsync on OSX. +# # This OSX-specific. Guess the standard package managers have it available on Linux. +$ brew install rsync +$ # Fetch seqrepo data (should take about 16 minutes) +$ seqrepo --rsync-exe /opt/homebrew/bin/rsync -r . pull --update-latest +$ # If you'll get a "Permission denied" error, then you can run the following command (using the temp directory which +$ # got created): +$ # > chmod +w 2024-02-20.r4521u5y && mv 2024-02-20.r4521u5y 2024-02-20 && ln -s 2024-02-20 latest +$ +$ # Exit venv and cd to genomics-operations repo. +$ +$ # Pack the utilities data (should take about 25 minutes) +$ python ./utilities/pack_utilities_data.py +``` +You should see a warning in the output log if the current `PYARD_DATABASE_VERSION` is outdated and you can change +`PYARD_DATABASE_VERSION` in the `.env` file if you wish to switch to the latest version that is printed in this log. + +Now you should set a new value for `UTILITIES_DATA_VERSION` in the `.env` file, create a new branch and commit this +change in it. Then also create a git tag for this commit with the `UTILITIES_DATA_VERSION` value and push it to GitHub +along with the branch. Now you can use this tag to create a new [release](https://github.com/FHIR/genomics-operations/releases). +Inside this release, you need to attach all the `*.tar.gz` files from the `./tmp` folder which was created after +`pack_utilities_data.py` ran successfully. + +Once the release is published, create PR from this new branch and merge it. -- Run `pyard.init(data_dir='./data/pyard', imgt_version=)` to download the new version -- Run `cd data/pyard && tar -czf pyard.sqlite3.tar.gz pyard-.sqlite3` -- Upload `pyard.sqlite3.tar.gz` in a new release on GitHub -- Update `PYARD_DATABASE_VERSION` in `.env` -- Update `UTILITIES_DATA_VERSION` in `.env` with the new tag ID (short git sha) +Finally, in order to validate the new release locally, run `fetch_utilities_data.sh` locally to recreate the `data` +directory (delete it first if you have it already). diff --git a/app/__init__.py b/app/__init__.py index d806ce99e..621c1b14b 100644 --- a/app/__init__.py +++ b/app/__init__.py @@ -3,6 +3,13 @@ from flask_cors import CORS import os +import hgvs +# Disable the hgvs LRU cache to avoid blowing up memory +# TODO: Revisit this, since this caching might not use a ton of memory. +hgvs.global_config.lru_cache.maxsize = 0 +# Disable HGVS strict bounds checks as a workaround for liftover failures: https://github.com/biocommons/hgvs/issues/717 +hgvs.global_config.mapping.strict_bounds = False + def create_app(): # App and API diff --git a/app/api_spec.yml b/app/api_spec.yml index bf244c351..199e9ebc7 100644 --- a/app/api_spec.yml +++ b/app/api_spec.yml @@ -1330,6 +1330,64 @@ paths: type: string example: "NM_001127510.3:c.145A>T" + /utilities/normalize-variant-hgvs: + get: + summary: "Normalize Variant HGVS" + operationId: "app.utilities_endpoints.normalize_variant_hgvs" + tags: + - "Operations Utilities (not part of balloted HL7 Operations)" + responses: + "200": + description: "Returns a normalized variant in both GRCh37 and GRCh38." + content: + application/json: + schema: + type: object + parameters: + - name: variant + in: query + required: true + description: "Variant." + schema: + type: string + example: "NM_021960.4:c.740C>T" + + /utilities/seqfetcher/1/sequence/{acc}: + get: + summary: "Seqfetcher" + operationId: "app.utilities_endpoints.seqfetcher" + tags: + - "Operations Utilities (not part of balloted HL7 Operations)" + responses: + "200": + description: "Returns RefSeq subsequence" + content: + text/plain: + schema: + type: string + parameters: + - name: acc + in: path + required: true + description: Accession + schema: + type: string + example: "NC_000001.10" + - name: start + in: query + required: true + description: Subsequence start index + schema: + type: integer + example: 10000 + - name: end + in: query + required: true + description: Subsequence end index + schema: + type: integer + example: 10010 + /utilities/normalize-hla: get: description: > diff --git a/app/common.py b/app/common.py index 36cd85522..2799e5ae8 100644 --- a/app/common.py +++ b/app/common.py @@ -8,14 +8,15 @@ import pyliftover import pymongo -import requests from flask import abort -from utilities.pyard import redux +from utilities import hla + +from .input_normalization import normalize # MongoDB Client URIs -FHIR_genomics_data_client_uri = f"mongodb+srv://readonly:{os.getenv('MONGODB_READONLY_PASSWORD')}@cluster0.8ianr.mongodb.net/FHIRGenomicsData" -utilities_data_client_uri = f"mongodb+srv://readonly:{os.getenv('MONGODB_READONLY_PASSWORD')}@cluster0.8ianr.mongodb.net/UtilitiesData" +FHIR_genomics_data_client_uri = f"mongodb+srv://readonly:{os.environ['MONGODB_READONLY_PASSWORD']}@cluster0.8ianr.mongodb.net/FHIRGenomicsData" +utilities_data_client_uri = f"mongodb+srv://readonly:{os.environ['MONGODB_READONLY_PASSWORD']}@cluster0.8ianr.mongodb.net/UtilitiesData" # MongoDB Clients client = pymongo.MongoClient(FHIR_genomics_data_client_uri) @@ -122,10 +123,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 @@ -169,26 +166,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.") @@ -202,22 +179,22 @@ 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"]: + elif not SPDIs["GRCh37SPDI"] and not SPDIs["GRCh38SPDI"]: abort(400, f'Cannot normalize variant: {variant}') else: - normalized_variant = {"variant": variant, "GRCh37": SPDIs["GRCh37"], "GRCh38": SPDIs["GRCh38"]} + normalized_variant = {"variant": variant, "GRCh37": SPDIs["GRCh37SPDI"], "GRCh38": SPDIs["GRCh38SPDI"]} 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"]: + elif not SPDIs["GRCh37SPDI"] and not SPDIs["GRCh38SPDI"]: abort(400, f'Cannot normalize variant: {variant}') else: - normalized_variant = {"variant": variant, "GRCh37": SPDIs["GRCh37"], "GRCh38": SPDIs["GRCh38"]} + normalized_variant = {"variant": variant, "GRCh37": SPDIs["GRCh37SPDI"], "GRCh38": SPDIs["GRCh38SPDI"]} else: abort(400, f'variant ({variant}) is not in the correct format(SPDI|HGVS)') @@ -261,7 +238,7 @@ def get_haplotype(haplotype): haplotype = haplotype.strip() haplotype_return = {'isSystem': False, 'haplotype': haplotype, 'system': None, 'lgxHaplotype': None} try: - haplotype_return['lgxHaplotype'] = redux(haplotype, "lgx") + haplotype_return['lgxHaplotype'] = hla.redux(haplotype, "lgx") except Exception: haplotype_return['lgxHaplotype'] = None if "|" in haplotype: @@ -272,7 +249,7 @@ def get_haplotype(haplotype): haplotype_return['haplotype'] = haplotype haplotype_return['system'] = haplotype_system_url try: - haplotype_return['lgxHaplotype'] = redux(haplotype, "lgx") + haplotype_return['lgxHaplotype'] = hla.redux(haplotype, "lgx") except Exception: haplotype_return['lgxHaplotype'] = None @@ -1053,136 +1030,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 5baef65c8..eaa1fea08 100644 --- a/app/endpoints.py +++ b/app/endpoints.py @@ -180,7 +180,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 = {} @@ -831,13 +835,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() @@ -1095,13 +1109,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() @@ -1464,7 +1488,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 = {} @@ -1973,7 +2001,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: @@ -2185,7 +2217,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..b946fbbd5 --- /dev/null +++ b/app/input_normalization.py @@ -0,0 +1,201 @@ +import os + +import hgvs.assemblymapper +import hgvs.dataproviders.uta +import hgvs.parser + +from utilities.spdi_refseq import normalize as spdi_normalize + +from . import common + +# Set the HGVS_SEQREPO_URL env var so the hgvs library will use the local `utilities/seqfetcher` endpoint instead of +# making NCBI API calls. +port = os.environ('PORT', 5000) # The localhost debugger starts the app on port 5000 +os.environ('HGVS_SEQREPO_URL') = f"http://127.0.0.1:{port}/utilities/seqfetcher" + +database_schema = os.environ['UTA_DATABASE_SCHEMA'] +# Use the biocommons UTA database if we don't specify a custom one. +# Also, make sure the URL uses `postgresql` instead of `postgres` as schema +database_url = f"{os.environ['UTA_DATABASE_URL']}/{database_schema}".replace('postgres://', 'postgresql://') + +hgvsParser = hgvs.parser.Parser() +hgvsDataProvider = hgvs.dataproviders.uta.connect(db_url=database_url) + +# Note: One can set `replace_reference=False` and `prevalidation_level=None` to skip data validation +# Also, until https://github.com/biocommons/hgvs/issues/704 is addressed, the following config settings need to also be set. +# hgvs.config.global_config.normalizer.validate = False +# hgvs.global_config.mapping.prevalidation_level = None # TODO: Open issue +# hgvs.global_config.mapping.replace_reference = False # TODO: Open issue +# Until https://github.com/biocommons/hgvs/issues/705 is addressed, validation cannot be disabled completely. + +b37hgvsAssemblyMapper = hgvs.assemblymapper.AssemblyMapper( + hgvsDataProvider, assembly_name='GRCh37', alt_aln_method='splign') +b38hgvsAssemblyMapper = hgvs.assemblymapper.AssemblyMapper( + hgvsDataProvider, assembly_name='GRCh38', alt_aln_method='splign') + +# ------------------ 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): + 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 spdi_normalize(parsed_variant.ac, pos, ref, alt) + +# ---------------- CONVERT NM_HGVS to canonical SPDIs --------------- + + +def process_NM_HGVS(NM_HGVS): + parsed_variant = hgvsParser.parse_hgvs_variant(NM_HGVS) + + projected_variant_dict = project_variant(parsed_variant) + print( + f"b37projected: {projected_variant_dict['b37projected']}; b38projected: {projected_variant_dict['b38projected']}") + + b37SPDI = normalize_variant(projected_variant_dict['b37projected']) + b38SPDI = normalize_variant(projected_variant_dict['b38projected']) + print(f"b37normalized: {b37SPDI}; b38normalized: {b38SPDI}") + + return {'GRCh37SPDI': b37SPDI, 'GRCh38SPDI': b38SPDI, 'GRCh37HGVS': str(projected_variant_dict['b37projected']), 'GRCh38HGVS': str(projected_variant_dict['b38projected'])} + + +# ---------------- CONVERT NC_HGVS to canonical SPDIs --------------- + + +def process_NC_HGVS(NC_HGVS): + parsed_variant = hgvsParser.parse_hgvs_variant(NC_HGVS) + + try: + transcripts = b38hgvsAssemblyMapper.relevant_transcripts(parsed_variant) + # Since we might not have access to all the transcripts that UTA contains, we select the "Longest Compatible + # Remaining" transcript as described here: https://github.com/GenomicMedLab/cool-seq-tool/blob/main/docs/TranscriptSelectionPriority.md + # - If there is a tie, choose the first-published transcript (lowest-numbered accession for RefSeq/Ensembl) among + # those transcripts meeting this criterion. + # Note: We want the most recent version of a transcript associated with an assembly. + # Eventually, hgvs should include pyliftover for this purpose: https://github.com/biocommons/hgvs/issues/711 + nm_transcripts = [t for t in transcripts if 'NM_' in t] + # Since a transcript looks like 'NM_006015.6', we sort them based on the number following 'NM_' + ordered_transcripts = sorted(nm_transcripts, key=lambda t: int(t[3:t.find('.')])) + # Pick the first transcript accession + transcript_acc = ordered_transcripts[0].split('.')[0] + # Search for its most recent version among nm_transcripts + relevant_transcript = max((t for t in nm_transcripts if t.startswith(transcript_acc)), key=lambda t: int(t.split('.')[-1])) + + var_c = b38hgvsAssemblyMapper.g_to_c( + parsed_variant, relevant_transcript) + + 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']) + b38SPDI = normalize_variant(projected_variant_dict['b38projected']) + print(f"b37normalized: {b37SPDI}; b38normalized: {b38SPDI}") + + except Exception: + provided_genomic_build = common.get_build_and_chrom_by_ref_seq(parsed_variant.ac) + liftover = common.lift_over(parsed_variant.ac, str(parsed_variant.posedit.pos.start), str(parsed_variant.posedit.pos.end)) + iv = hgvs.location.Interval(start=liftover['start'], end=liftover['end']) + posedit = hgvs.posedit.PosEdit(pos=iv, edit=parsed_variant.posedit.edit) + lifted_variant = hgvs.sequencevariant.SequenceVariant(ac=liftover['refSeq'], type="g", posedit=posedit) + if provided_genomic_build['build'] == 'GRCh37': + b37SPDI = normalize_variant(parsed_variant) + b38SPDI = normalize_variant(lifted_variant) + else: + b37SPDI = normalize_variant(lifted_variant) + b38SPDI = normalize_variant(parsed_variant) + + return {'GRCh37SPDI': b37SPDI, 'GRCh38SPDI': b38SPDI, 'GRCh37HGVS': str(projected_variant_dict['b37projected']), 'GRCh38HGVS': str(projected_variant_dict['b38projected'])} + +# ---------------- CONVERT NM_SPDI to canonical SPDIs --------------- + + +def process_NM_SPDI(NM_SPDI): + # convert SPDI into NM_HGVS then use NM_HGVS pipeline + ref_seq = 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( + ref_seq+":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( + ref_seq+":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( + ref_seq+":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( + ref_seq+":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 + ref_seq = 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 = (ref_seq+":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 = (ref_seq+":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 = (ref_seq+":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 = (ref_seq+":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/utilities_endpoints.py b/app/utilities_endpoints.py index cfd1b13cd..701e2d813 100644 --- a/app/utilities_endpoints.py +++ b/app/utilities_endpoints.py @@ -1,12 +1,12 @@ - import json from collections import OrderedDict import requests from flask import abort, jsonify -from app import common -from utilities.pyard import redux +from app import common, input_normalization +from utilities import hla +from utilities.spdi_refseq import get_ref_seq_subseq def fetch_concept_map(mapID): @@ -172,7 +172,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})") @@ -356,18 +357,34 @@ def translate_terminology(codeSystem, code): return response +def normalize_variant_hgvs(variant): + try: + return input_normalization.normalize(variant) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(422, 'Failed LiftOver') + + +def seqfetcher(acc, start, end): + try: + return get_ref_seq_subseq(acc, start, end) + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + abort(404, 'Not Found') + + def normalize_hla(allele): try: return { allele: { - "G": redux(allele, "G"), - "P": redux(allele, "P"), - "lg": redux(allele, "lg"), - "lgx": redux(allele, "lgx"), - "W": redux(allele, "W"), - "exon": redux(allele, "exon"), - "U2": redux(allele, "U2"), - "S": redux(allele, "S") + "G": hla.redux(allele, "G"), + "P": hla.redux(allele, "P"), + "lg": hla.redux(allele, "lg"), + "lgx": hla.redux(allele, "lgx"), + "W": hla.redux(allele, "W"), + "exon": hla.redux(allele, "exon"), + "U2": hla.redux(allele, "U2"), + "S": hla.redux(allele, "S") } } except Exception as err: diff --git a/fetch_utilities_data.sh b/fetch_utilities_data.sh index b215dd81b..e31be189d 100755 --- a/fetch_utilities_data.sh +++ b/fetch_utilities_data.sh @@ -9,6 +9,25 @@ if [ -d ./data ]; then exit 0 fi +mkdir -p ./data/refseq +( + cd ./data/refseq + + echo "Downloading refseq files..." + + curl -sLO https://github.com/FHIR/genomics-operations/releases/download/${UTILITIES_DATA_VERSION}/GRCh37_refseq.tar.gz + curl -sLO https://github.com/FHIR/genomics-operations/releases/download/${UTILITIES_DATA_VERSION}/GRCh38_refseq.tar.gz + curl -sLO https://github.com/FHIR/genomics-operations/releases/download/${UTILITIES_DATA_VERSION}/rna_refseq.tar.gz + + echo "Extracting refseq files..." + + tar -xzf GRCh37_refseq.tar.gz + tar -xzf GRCh38_refseq.tar.gz + tar -xzf rna_refseq.tar.gz + + echo "Finished extracting refseq files." +) + mkdir -p ./data/pyard ( cd ./data/pyard diff --git a/requirements.txt b/requirements.txt index 9c07a31b7..6216fe975 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,8 +11,13 @@ pyranges==0.0.120 python_dateutil==2.8.2 requests==2.28.1 streamlit==1.19.0 +<<<<<<< HEAD streamlit-aggrid==1.0.5 pytest==7.1.1 +======= +pytest==7.1.1 +deepdiff==6.7.1 +>>>>>>> 1b4630a (Add SPDI normalization with refseq files) python-dotenv==1.0.1 toml==0.10.2 pandas==2.2.2 @@ -20,6 +25,9 @@ numpy==2.0.2 pyliftover==0.4 hgvs==1.5.4 py-ard==1.5.3 +<<<<<<< HEAD openpyxl==3.1.5 pyvcf3==1.0.3 tqdm==4.65.1 +======= +>>>>>>> 1b4630a (Add SPDI normalization with refseq files) diff --git a/tests/conftest.py b/tests/conftest.py index 2eafd8e38..074274d02 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,3 +1,5 @@ +from threading import Thread + import pytest from dotenv import load_dotenv @@ -8,7 +10,20 @@ load_dotenv("secrets.env") +@pytest.fixture(scope='session') +def app(): + # Create app with shutdown hook + app = create_app() + + # Hack: Start Flask Werkzeug server as a daemon thread and let it dangle (skip `thread.join()`). When the tests + # finish, it will be the only thread left running and the Python process exits when only daemon threads are left. + # Details here: https://docs.python.org/3/library/threading.html#thread-objects + thread = Thread(target=app.run, daemon=True, kwargs=dict(port=5000)) + thread.start() + + yield app + + @pytest.fixture(scope='module') -def client(): - with create_app().test_client() as c: - yield c +def client(app): + return app.test_client() diff --git a/tests/expected_outputs/find_population_specific_variants/3.json b/tests/expected_outputs/find_population_specific_variants/3.json index dfaab8c5a..48d11ee97 100644 --- a/tests/expected_outputs/find_population_specific_variants/3.json +++ b/tests/expected_outputs/find_population_specific_variants/3.json @@ -17,7 +17,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1090 + "value": 1115 } } ] @@ -38,7 +38,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1090 + "value": 1115 } } ] @@ -59,7 +59,7 @@ { "name": "denominator", "valueQuantity": { - "value": 1090 + "value": 1115 } } ] diff --git a/tests/integration_tests/test_normalize_variant_utility.py b/tests/integration_tests/test_normalize_variant_utility.py new file mode 100644 index 000000000..261786e03 --- /dev/null +++ b/tests/integration_tests/test_normalize_variant_utility.py @@ -0,0 +1,26 @@ +import tests.utilities as tu + + +""" +Normalize Variant Utility Tests +-------------------------------- +""" + + +def test_normalize_variant_utility(client): + url = tu.normalize_variant_query('variant=NM_021960.4%3Ac.740C>T') + response = client.get(url) + + assert response.status_code == 200 + + +def test_normalize_variant_utility_transcript_ref_disagree_position_variant(client): + """ + This is a specific corner case where liftover fails if HGVS strict bounds checks are enforced + Details here: https://github.com/biocommons/hgvs/issues/717 + """ + url = tu.normalize_variant_query('variant=NC_000001.10:g.145592073A>T') + + response = client.get(url) + + assert response.status_code == 200 diff --git a/tests/integration_tests/test_population_genotype_operations.py b/tests/integration_tests/test_population_genotype_operations.py index c3327d55b..6dcfbe7d5 100644 --- a/tests/integration_tests/test_population_genotype_operations.py +++ b/tests/integration_tests/test_population_genotype_operations.py @@ -22,6 +22,12 @@ 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 returned + "HGVSDataNotAvailableError: No alignments for NM_001002811.1 in GRCh38 using splign" before the relevant_transcript + selection for liftover was adjusted to select the newest version of the first-published transcript in + `process_NC_HGVS()` in input_normalization.py. + """ 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) @@ -29,11 +35,14 @@ def test_population_specific_variants_2(client): tu.compare_actual_and_expected_output(f'{tu.FIND_POPULATION_SPECIFIC_VARIANTS_OUTPUT_DIR}2.json', response.json) -# def test_population_specific_variants_3(client): -# url = tu.find_population_specific_variants_query('variants=NC_000001.10:144931726:G:A, NC_000001.10:145532548:T:C, NC_000001.10:145592072:A:T') -# response = client.get(url) +def test_population_specific_variants_3(client): + """ + Same as test 2 above, but without `subject` & `includePatientList`. + """ + url = tu.find_population_specific_variants_query('variants=NC_000001.10:144931726:G:A, NC_000001.10:145532548:T:C, NC_000001.10:145592072:A:T') + response = client.get(url) -# tu.compare_actual_and_expected_output(f'{tu.FIND_POPULATION_SPECIFIC_VARIANTS_OUTPUT_DIR}3.json', response.json) + tu.compare_actual_and_expected_output(f'{tu.FIND_POPULATION_SPECIFIC_VARIANTS_OUTPUT_DIR}3.json', response.json) """ diff --git a/tests/unit_tests/test_seqrepo_data_packing.py b/tests/unit_tests/test_seqrepo_data_packing.py new file mode 100644 index 000000000..e3056a78c --- /dev/null +++ b/tests/unit_tests/test_seqrepo_data_packing.py @@ -0,0 +1,46 @@ +import pytest + +from utilities.pack_utilities_data import pack_ref_seq +from utilities.spdi_refseq import RefSeq + + +@pytest.mark.parametrize( + "name, ref_seq, exp_packed_ref_seq", + [ + ('happy flow', 'ACGT', b'\x88\xf6'), + ('empty RefSeq', '', b''), + ('2 nucleotides -> one byte', 'GT', b'\xda'), + ('3 nucleotides -> two bytes', 'GTA', b'\x1a\xfe'), + ('5 nucleotides -> two bytes', 'GTACT', b'\x1a\xb2'), + ('6 nucleotides -> three bytes', 'GTACTG', b'\x1a\x32\xfd'), + ('8 nucleotides -> three bytes', 'GTACTGAT', b'\x1a\x32\x61'), + ('9 nucleotides -> four bytes', 'GTACTGATC', b'\x1a\x32\x61\xf9'), + ('all other codes get mapped to N (including invalid ones)', 'FO0b4R', b'\xff\xff\xff') + ], +) +def test_pack_ref_seq(name, ref_seq, exp_packed_ref_seq): + assert pack_ref_seq(ref_seq) == exp_packed_ref_seq, name + + +@pytest.mark.parametrize( + "name, packed_ref_seq, exp_ref_seq", + [ + ('happy flow', b'\x88\xf6', 'ACGT'), + ('empty RefSeq', b'', ''), + ('2 nucleotides -> one byte', b'\xda', 'GT'), + ('3 nucleotides -> two bytes', b'\x1a\xfe', 'GTA'), + ('5 nucleotides -> two bytes', b'\x1a\xb2', 'GTACT'), + ('6 nucleotides -> three bytes', b'\x1a\x32\xfd', 'GTACTG'), + ('8 nucleotides -> three bytes', b'\x1a\x32\x61', 'GTACTGAT'), + ('9 nucleotides -> four bytes', b'\x1a\x32\x61\xf9', 'GTACTGATC'), + ('Ns are decoded correctly', b'\xff\xff\xff', 'NNNNNN') + ], +) +def test_unpack_ref_seq(name, packed_ref_seq, exp_ref_seq): + assert RefSeq(packed_ref_seq, len(exp_ref_seq))[:] == exp_ref_seq, name + + +def test_unpack_ref_seq_error(): + with pytest.raises(NotImplementedError) as exc_info: + RefSeq(b'\x06', 1)[:] + assert str(exc_info.value) == "unexpected hex-encoded code: '6'" diff --git a/tests/utilities.py b/tests/utilities.py index 2cff8a27e..e0f176d33 100644 --- a/tests/utilities.py +++ b/tests/utilities.py @@ -63,6 +63,9 @@ NORMALIZE_HLA_URL = "/utilities/normalize-hla" NORMALIZE_HLA_OUTPUT_DIR = "tests/expected_outputs/normalize_hla/" +NORMALIZE_VARIANT_URL = "/utilities/normalize-variant" +NORMALIZE_VARIANT_OUTPUT_DIR = "tests/expected_outputs/normalize_variant/" + SORT_JSON_DATA_OUTPUT_DIR = "tests/expected_outputs/sort_json_data/" @@ -146,6 +149,10 @@ def normalize_hla_utility_query(query): return f"{NORMALIZE_HLA_URL}?{query}" +def normalize_variant_query(query): + return f"{NORMALIZE_VARIANT_URL}?{query}" + + def sort_json_data(data): """ Recursively sorts the keys of a JSON object or the elements of a JSON list. diff --git a/utilities/SPDI_Normalization.py b/utilities/SPDI_Normalization.py deleted file mode 100644 index 1424281d9..000000000 --- a/utilities/SPDI_Normalization.py +++ /dev/null @@ -1,163 +0,0 @@ -import pyfastx -from threading import Lock - - -# 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): - 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] - - # Step 0 - start = pos - end = pos + len(ref) - - # Step 1 - # a. Trim Common Suffix - - small = len(alt) if len(alt) <= len(ref) else len(ref) - - suffix_length = 0 - - last = -1 - while small != 0: - if ref[last] == alt[last]: - suffix_length += 1 - end -= 1 - last -= 1 - small -= 1 - else: - break - - if suffix_length > 0: - ref = ref[0:-suffix_length] - alt = alt[0:-suffix_length] - - # b. Trim Common Prefix - - small = len(alt) if len(alt) <= len(ref) else len(ref) - - prefix_length = 0 - - first = 0 - while small != 0: - if ref[first] == alt[first]: - prefix_length += 1 - start += 1 - first += 1 - small -= 1 - else: - break - - if prefix_length > 0: - ref = ref[prefix_length:] - alt = alt[prefix_length:] - - # Step 2 - # a. Check if both are empty - # Not Possible in this case - - # b. Check if both are non-empty - - if ref and alt: - return f"{ref_seq}:{start}:{ref}:{alt}" - - # c. Check if one is empty - - if alt: - allele = alt - type_ = 'I' - else: - allele = ref - type_ = 'D' - - # Step 3 - # a. Left Roll - - if type_ == 'I': - left_roll_bound = start - - temp_allele = allele - - while (str(ref_seq_fasta[left_roll_bound-1]).upper() == temp_allele[-1].upper()): - temp_allele = str(ref_seq_fasta[left_roll_bound-1]).upper() + temp_allele[:-1] - left_roll_bound -= 1 - - # . b. Right Roll - - right_roll_bound = start - - temp_allele = allele - - while (str(ref_seq_fasta[right_roll_bound]).upper() == temp_allele[0].upper()): - temp_allele = temp_allele[1:] + str(ref_seq_fasta[right_roll_bound]).upper() - right_roll_bound += 1 - - # Step 4 - # a. Prepend - if ref_seq_fasta[left_roll_bound:start]: - ref = str(ref_seq_fasta[left_roll_bound:start]) + ref - alt = str(ref_seq_fasta[left_roll_bound:start]) + alt - - # b. append - if ref_seq_fasta[start:right_roll_bound]: - ref = ref + str(ref_seq_fasta[start:right_roll_bound]) - alt = alt + str(ref_seq_fasta[start:right_roll_bound]) - - # c. Modify Start and End, and Return SPDI - start = left_roll_bound - end = right_roll_bound - - return (f"{ref_seq}:{start}:{ref.upper()}:{alt.upper()}") - - else: - left_roll_bound = start - - while (str(ref_seq_fasta[left_roll_bound-1]).upper() == str(ref_seq_fasta[left_roll_bound + (len(allele)-1)]).upper()): - left_roll_bound -= 1 - - # . b. Right Roll - - right_roll_bound = start - - while (str(ref_seq_fasta[right_roll_bound]).upper() == (ref_seq_fasta[right_roll_bound + (len(allele))]).upper()): - right_roll_bound += 1 - - # Step 4 - # a. Prepend - if ref_seq_fasta[left_roll_bound:start]: - ref = str(ref_seq_fasta[left_roll_bound:start]) + ref - alt = str(ref_seq_fasta[left_roll_bound:start]) + alt - - # b. append - if ref_seq_fasta[start:right_roll_bound]: - ref = ref + str(ref_seq_fasta[start:right_roll_bound]) - alt = alt + str(ref_seq_fasta[start:right_roll_bound]) - - # c. Modify Start and End, and Return SPDI - start = left_roll_bound - end = right_roll_bound - - return (f"{ref_seq}:{start}:{ref.upper()}:{alt.upper()}") diff --git a/utilities/common.py b/utilities/common.py index 025bd9229..207e83f23 100644 --- a/utilities/common.py +++ b/utilities/common.py @@ -1,18 +1,17 @@ import os import re from enum import Enum +from pathlib import Path import pandas as pd import pymongo - -from pathlib import Path from dotenv import load_dotenv BASE_DIR = Path(__file__).resolve().parent PROJECT_ROOT = BASE_DIR.parent load_dotenv(PROJECT_ROOT / "secrets.env") -utilities_data_client_uri = f"mongodb+srv://readonly:{os.getenv('MONGODB_READONLY_PASSWORD')}@cluster0.8ianr.mongodb.net/UtilitiesData" +utilities_data_client_uri = f"mongodb+srv://readonly:{os.environ('MONGODB_READONLY_PASSWORD')}@cluster0.8ianr.mongodb.net/UtilitiesData" utilities_client = pymongo.MongoClient(utilities_data_client_uri) utilities_db = utilities_client.UtilitiesData transcript_data = utilities_db.Transcripts diff --git a/utilities/gene_ref_seq.py b/utilities/gene_ref_seq.py index b9e9d0c63..ec806df15 100644 --- a/utilities/gene_ref_seq.py +++ b/utilities/gene_ref_seq.py @@ -56,7 +56,7 @@ } -def _get_ref_seq_by_chrom(build, chrom): +def get_ref_seq_by_chrom(build, chrom): if build in ["GRCh37"]: return (b37[chrom]) diff --git a/utilities/pyard.py b/utilities/hla.py similarity index 64% rename from utilities/pyard.py rename to utilities/hla.py index 9847209bb..87fec860f 100644 --- a/utilities/pyard.py +++ b/utilities/hla.py @@ -1,14 +1,14 @@ import os from os.path import isdir -import pyard +from pyard import init as pyard_init # Make sure the pyard folder exists locally if not isdir('./data/pyard'): exit("Missing pyard folder. Please run fetch_utilities_data.sh!") -pyard_database_version = os.getenv('PYARD_DATABASE_VERSION', '3580') -ard = pyard.init(data_dir='./data/pyard', cache_size=1, imgt_version=pyard_database_version) +pyard_database_version = os.environ['PYARD_DATABASE_VERSION'] +ard = pyard_init(data_dir='./data/pyard', cache_size=1, imgt_version=pyard_database_version) def redux(glstring, redux_type="lgx"): diff --git a/utilities/pack_utilities_data.py b/utilities/pack_utilities_data.py new file mode 100644 index 000000000..8e514ac83 --- /dev/null +++ b/utilities/pack_utilities_data.py @@ -0,0 +1,187 @@ +import argparse +import os +import pathlib +import tarfile +from math import ceil + +import hgvs.assemblymapper +import hgvs.dataproviders.uta +import hgvs.parser +from biocommons.seqrepo import SeqRepo +from dotenv import load_dotenv +from pyard import init as pyard_init +from pyard.load import load_latest_version as pyard_load_latest_version + +load_dotenv() +# Load secrets from secrets.env file if available +load_dotenv("secrets.env") + + +def code_to_hex(code): + """ + Convert sequence code to a 3 bit integer. + Supports only 6 codes: A, C, G, T, U and a generic placeholder (designated N) for everything else. + Details here: https://en.wikipedia.org/wiki/FASTA_format#Sequence_representation) + + Nucleic Acid Code | Meaning | Mnemonic + A | A | Adenine + C | C | Cytosine + G | G | Guanine + T | T | Thymine + U | U | Uracil + N | A C G T U | Nucleic acid + """ + match code.upper(): + case 'A': return 0x0 + case 'C': return 0x1 + case 'G': return 0x2 + case 'T': return 0x3 + case 'U': return 0x4 + # All other codes including `N`s + case _: return 0x7 + + +def encode_8_codes_to_3_bytes(codes, byte_count, ba): + # Pad the codes list with extra `N`s if its length is less than 8. This can happen at the end of a RefSeq and these + # extra codes get discarded by the decoder which receives the original length of the RefSeq, thus knowing when to + # stop decoding. + input_codes_len = len(codes) + if input_codes_len < 8: + codes = codes + 'N' * (8 - len(codes)) + + # We need 3 bits to represent each code and, since we need to pack them into bytes, some of them will end up being + # split across 2 bytes. The pattern repeats after LCM(3,8) = 24 (3 bytes), when the next code will start at index 0 + # in the next byte. + + c1 = code_to_hex(codes[0]) + c2 = code_to_hex(codes[1]) << 3 + # Chop off the top bit of codes[2] since it will be stored in the next byte. + c3 = (code_to_hex(codes[2]) << 6) & 0xff + ba[byte_count + 0] = c3 | c2 | c1 + + # Bail out if there are no more codes left for the next byte. + if input_codes_len < 3: + return + + c1 = code_to_hex(codes[2]) >> 2 + c2 = code_to_hex(codes[3]) << 1 + c3 = code_to_hex(codes[4]) << 4 + # Chop off the top 2 bits of codes[5] since they'll be stored in the next byte. + c4 = (code_to_hex(codes[5]) << 7) & 0xff + ba[byte_count + 1] = c4 | c3 | c2 | c1 + + # Bail out if there are no more codes left for the next byte. + if input_codes_len < 6: + return + + c1 = code_to_hex(codes[5]) >> 1 + c2 = code_to_hex(codes[6]) << 2 + c3 = code_to_hex(codes[7]) << 5 + ba[byte_count + 2] = c3 | c2 | c1 + + +def pack_ref_seq(ref_seq): + # Packing the RefSeq using 3 bits per code will yield a packed RefSeq that's 3/8 times shorter. + packed_len = ceil(len(ref_seq) * 3 / 8) + ba = bytearray(packed_len) + + # Pack each sequence of 8 codes into 3 bytes. + byte_count = 0 + for i in range(0, len(ref_seq), 8): + encode_8_codes_to_3_bytes(ref_seq[i:i + 8], byte_count, ba) + byte_count += 3 + + return ba + + +def compress_ref_seq(path): + with tarfile.open(f'{path}_refseq.tar.gz', 'w:gz') as tar: + tar.add(path, arcname=os.path.basename(path)) + + +def pack_ref_seq_to_file(seq_repo, acc, output_dir): + # Coalesce the entire sequence into a string to load it all in memory. + # Otherwise, the data packing will be very slow. + seq = str(seq_repo[acc]) + + with open(os.path.join(output_dir, f'{acc}_{len(seq)}.refseq'), 'wb') as binary_file: + binary_file.write(pack_ref_seq(seq)) + + +def pack_rnas(seq_repo, output_dir): + ref_seq_dir = os.path.join(output_dir, 'rna') + pathlib.Path(ref_seq_dir).mkdir(parents=True, exist_ok=True) + + uta_rnas = hgvsDataProvider._fetchall(r"""select distinct(split_part(ac, '.', 1) || '.' || regexp_replace(split_part(ac, '.', 2), '/\d+$', '')) from transcript where ac ~ '^NM_'""") + for rna in [rna[0] for rna in uta_rnas]: + pack_ref_seq_to_file(seq_repo, rna, ref_seq_dir) + + compress_ref_seq(ref_seq_dir) + + +def pack_chromosomes(seq_repo, output_dir, build): + ref_seq_dir = os.path.join(output_dir, build) + pathlib.Path(ref_seq_dir).mkdir(parents=True, exist_ok=True) + + for acc in [acc for acc in hgvsDataProvider.get_assembly_map(build) if acc.startswith('NC_')]: + pack_ref_seq_to_file(seq_repo, acc, ref_seq_dir) + + compress_ref_seq(ref_seq_dir) + + +def compress_pyard_data(output_dir): + pyard_dir = os.path.join(output_dir, 'pyard') + pathlib.Path(pyard_dir).mkdir(parents=True, exist_ok=True) + + db_version = os.environ['PYARD_DATABASE_VERSION'] + + latest_db_version = pyard_load_latest_version() + if db_version != latest_db_version: + print(f'Warning: Latest PyARD database version is {latest_db_version}, but PYARD_DATABASE_VERSION={db_version}') + + pyard_init(data_dir=pyard_dir, imgt_version=db_version) + + with tarfile.open(os.path.join(output_dir, 'pyard.sqlite3.tar.gz'), 'w:gz') as tar: + file = f'pyard-{db_version}.sqlite3' + tar.add(os.path.join(pyard_dir, file), arcname=file) + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Pack genome and RNA refseq data from seqrepo') + + parser.add_argument('--uta_database_schema', + help='UTA database schema', + default=os.environ['UTA_DATABASE_SCHEMA']) + parser.add_argument('--uta_database_url', help='UTA database URL', + # Use the biocommons UTA database if we don't specify a custom one. + default=os.environ['UTA_DATABASE_URL']) + parser.add_argument('--seqrepo_dir', + help='Seqrepo directory', + default='./seqrepo/latest') + parser.add_argument('--output_dir', + help='Output directory', + default='tmp') + args = parser.parse_args() + + # Also, make sure the URL uses `postgresql` instead of `postgres` as schema + database_url = f"{args.uta_database_url}/{args.uta_database_schema}".replace('postgres://', 'postgresql://') + + hgvsDataProvider = hgvs.dataproviders.uta.connect(db_url=database_url) + + seq_repo = SeqRepo(args.seqrepo_dir) + + pathlib.Path(args.output_dir).mkdir(parents=True, exist_ok=True) + + print('Packing GRCh37 refseq...') + pack_chromosomes(seq_repo, args.output_dir, 'GRCh37') + + print('Packing GRCh38 refseq...') + pack_chromosomes(seq_repo, args.output_dir, 'GRCh38') + + print('Packing RNA refseq...') + pack_rnas(seq_repo, args.output_dir) + + print('Packing pyard data...') + compress_pyard_data(args.output_dir) + + print('All done!') diff --git a/utilities/requirements.txt b/utilities/requirements.txt new file mode 100644 index 000000000..a8d982266 --- /dev/null +++ b/utilities/requirements.txt @@ -0,0 +1,10 @@ +pyfastx==0.8.4 +pymongo[srv]==4.2.0 +pyranges==0.0.120 +python_dateutil==2.8.2 +pyvcf3==1.0.3 +requests==2.28.1 +tqdm==4.65.1 +pandas==2.2.2 +numpy==2.0.2 +bioutils==0.6.1 diff --git a/utilities/spdi.py b/utilities/spdi.py new file mode 100644 index 000000000..7c2103f7f --- /dev/null +++ b/utilities/spdi.py @@ -0,0 +1,12 @@ +from bioutils.normalize import NormalizationMode +from bioutils.normalize import normalize as spdi_normalize + + +def normalize(ref_seq, acc, pos, ref, alt): + new_ival, new_alleles = spdi_normalize(sequence=ref_seq, + interval=(pos, pos + len(ref)), + alleles=(None, alt), + mode=NormalizationMode.EXPAND, + anchor_length=0, + trim=ref != alt) + return (f"{acc}:{new_ival[0]}:{new_alleles[0]}:{new_alleles[1]}") diff --git a/utilities/spdi_refseq.py b/utilities/spdi_refseq.py new file mode 100644 index 000000000..b901290d6 --- /dev/null +++ b/utilities/spdi_refseq.py @@ -0,0 +1,109 @@ +from glob import glob +from math import floor +from os.path import isdir +from pathlib import Path +from sys import exit +from threading import Lock + +from .spdi import normalize as spdi_normalize + +# Make sure the refseq folder exists locally +if not isdir('./data/refseq'): + exit("Missing refseq folder. Please run fetch_utilities_data.sh!") + + +def hex_to_code(hex): + """ + Convert a 3 bit integer to a sequence code. See `code_to_hex()` for details. + """ + match hex: + case 0x0: return 'A' + case 0x1: return 'C' + case 0x2: return 'G' + case 0x3: return 'T' + case 0x4: return 'U' + case 0x7: return 'N' + # Error out if we get any unexpected hex-encoded code + case _: + raise NotImplementedError(f"unexpected hex-encoded code: '{hex}'") + + +class RefSeq: + 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 __len__(self): + return self.__len + + def __getitem__(self, subscript): + if isinstance(subscript, slice): + codes = [] + # 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) + + 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) & 0x7) + + +def get_ref_seq(acc): + try: + ref_seq_file_pattern = f'./data/refseq/**/{acc}_*.refseq' + found_files = glob(ref_seq_file_pattern) + # TODO: Don't store NM accessions for each build since they're redundant + if not found_files: + raise ValueError(f'failed to find expected refseq file for accession "{acc}"') + + ref_seq_file = found_files[0] + ref_seq_length = int(Path(ref_seq_file).stem.rpartition('_')[-1]) + + with open(ref_seq_file, mode='rb') as file: + ref_seq = RefSeq(file.read(), ref_seq_length) + except Exception as err: + print(f"failed to read refseq file: {err=}, {type(err)=}") + raise + + return ref_seq + + +ref_seq_lock = Lock() + + +def get_ref_seq_subseq(acc, 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(acc)[start:end] + + +def normalize(acc, pos, ref, alt): + # Need to serialise this if we can't keep all the RefSeq data in memory + with ref_seq_lock: + ref_seq_fasta = get_ref_seq(acc) + return spdi_normalize(ref_seq=ref_seq_fasta, acc=acc, pos=pos, ref=ref, alt=alt) diff --git a/utilities/vcf2json.py b/utilities/vcf2json.py index 38a435a51..ac3d9fe34 100644 --- a/utilities/vcf2json.py +++ b/utilities/vcf2json.py @@ -1,11 +1,40 @@ -import vcf import json -from collections import OrderedDict -from gene_ref_seq import _get_ref_seq_by_chrom -from SPDI_Normalization import get_normalized_spdi -import common import re import uuid +from collections import OrderedDict +from threading import Lock + +import common +import pyfastx +import vcf + +from .gene_ref_seq import get_ref_seq_by_chrom +from .spdi import normalize + +# 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 normalize_spdi(ref_seq, pos, ref, alt, build): + fasta = get_fasta(BUILD37_FILE) if build == 'GRCh37' else get_fasta(BUILD38_FILE) + + return normalize(fasta[ref_seq], pos, ref, alt) def add_phase_records(record, phased_rec_map, sample_position): @@ -204,7 +233,7 @@ def vcf2json(vcf_filename=None, ref_build=None, patient_id=None, output_json["GT"] = GT output_json["genomicSourceClass"] = genomic_source_class # populate SPDI - ref_seq = _get_ref_seq_by_chrom(ref_build, common.extract_chrom_identifier(record.CHROM)) + ref_seq = get_ref_seq_by_chrom(ref_build, common.extract_chrom_identifier(record.CHROM)) if not record.is_sv and output_json["ALT"] is not None: spdi = (f'{ref_seq}:{record.POS - 1}:{record.REF}:{output_json["ALT"]}') @@ -214,7 +243,7 @@ def vcf2json(vcf_filename=None, ref_build=None, patient_id=None, # Calculate SPDI using SPDI_Normalization logic if len(record.REF) != len(output_json["ALT"]): - output_json["SPDI"] = get_normalized_spdi(ref_seq, (record.POS - 1), record.REF, output_json["ALT"], ref_build) + output_json["SPDI"] = normalize_spdi(ref_seq, (record.POS - 1), record.REF, output_json["ALT"], ref_build) # populate allelicState alleles = common.get_allelic_state(record, ratio_ad_dp, sample_position) From 305e0638ba39ce5e3c68dce44c9e6f64eb5b741b Mon Sep 17 00:00:00 2001 From: Mihai Todor Date: Thu, 20 Feb 2025 16:52:29 +0000 Subject: [PATCH 2/2] WIP Signed-off-by: Mihai Todor --- .vscode/launch.json | 1 + app/__init__.py | 5 ++-- app/endpoints.py | 55 ++++++++------------------------------ app/input_normalization.py | 4 +-- requirements.txt | 10 +------ utilities/common.py | 2 +- utilities/vcf2json.py | 7 +++-- 7 files changed, 22 insertions(+), 62 deletions(-) diff --git a/.vscode/launch.json b/.vscode/launch.json index 8d87e770a..f8d7e2015 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -23,6 +23,7 @@ "type": "debugpy", "request": "launch", "program": "${file}", + "cwd": "${fileDirname}", "console": "integratedTerminal", "justMyCode": false } diff --git a/app/__init__.py b/app/__init__.py index 621c1b14b..a7b719253 100644 --- a/app/__init__.py +++ b/app/__init__.py @@ -1,9 +1,10 @@ +import os + import connexion import flask +import hgvs from flask_cors import CORS -import os -import hgvs # Disable the hgvs LRU cache to avoid blowing up memory # TODO: Revisit this, since this caching might not use a ton of memory. hgvs.global_config.lru_cache.maxsize = 0 diff --git a/app/endpoints.py b/app/endpoints.py index eaa1fea08..fb59c5847 100644 --- a/app/endpoints.py +++ b/app/endpoints.py @@ -1,7 +1,8 @@ from collections import OrderedDict + from flask import abort, jsonify -from app import common -from app import utilities_endpoints + +from app import common, utilities_endpoints def find_subject_variants( @@ -180,11 +181,7 @@ def find_subject_specific_variants( subject = subject.strip() common.validate_subject(subject) - try: - variants = list(map(common.get_variant, variants)) - except Exception as err: - print(f"Unexpected {err=}, {type(err)=}") - abort(422, 'Failed LiftOver') + variants = list(map(common.get_variant, variants)) # Query query = {} @@ -836,22 +833,13 @@ def find_subject_tx_implications( ranges = list(map(common.get_range, ranges)) common.get_lift_over_range(ranges) - try: - variants = common.get_variants(ranges, query) - except Exception as err: - print(f"Unexpected {err=}, {type(err)=}") - abort(422, 'Failed LiftOver') - + variants = common.get_variants(ranges, query) if not variants: return jsonify({"resourceType": "Parameters"}) normalized_variants = [{variant["BUILD"]: variant["SPDI"]} for variant in variants] if variants and not ranges: - try: - normalized_variants = list(map(common.get_variant, variants)) - except Exception as err: - print(f"Unexpected {err=}, {type(err)=}") - abort(422, 'Failed LiftOver') + normalized_variants = list(map(common.get_variant, variants)) # Result Object result = OrderedDict() @@ -1110,22 +1098,13 @@ def find_subject_dx_implications( ranges = list(map(common.get_range, ranges)) common.get_lift_over_range(ranges) - try: - variants = common.get_variants(ranges, query) - except Exception as err: - print(f"Unexpected {err=}, {type(err)=}") - abort(422, 'Failed LiftOver') - + variants = common.get_variants(ranges, query) if not variants: return jsonify({"resourceType": "Parameters"}) normalized_variants = [{variant["BUILD"]: variant["SPDI"]} for variant in variants] if variants and not ranges: - try: - normalized_variants = list(map(common.get_variant, variants)) - except Exception as err: - print(f"Unexpected {err=}, {type(err)=}") - abort(422, 'Failed LiftOver') + normalized_variants = list(map(common.get_variant, variants)) # Result Object result = OrderedDict() @@ -1488,11 +1467,7 @@ def find_population_specific_variants( # Parameters variants = list(map(lambda x: x.strip().split(","), variants)) for i in range(len(variants)): - try: - variants[i] = list(map(common.get_variant, variants[i])) - except Exception as err: - print(f"Unexpected {err=}, {type(err)=}") - abort(422, 'Failed LiftOver') + variants[i] = list(map(common.get_variant, variants[i])) # Query query = {} @@ -2001,11 +1976,7 @@ def find_population_tx_implications( return jsonify({"resourceType": "Parameters"}) if variants: - try: - variants = list(map(common.get_variant, variants)) - except Exception as err: - print(f"Unexpected {err=}, {type(err)=}") - abort(422, 'Failed LiftOver') + variants = list(map(common.get_variant, variants)) condition_code_list = [] if conditions: @@ -2217,11 +2188,7 @@ def find_population_dx_implications( return jsonify({"resourceType": "Parameters"}) if variants: - try: - variants = list(map(common.get_variant, variants)) - except Exception as err: - print(f"Unexpected {err=}, {type(err)=}") - abort(422, 'Failed LiftOver') + variants = list(map(common.get_variant, variants)) condition_code_list = [] if conditions: diff --git a/app/input_normalization.py b/app/input_normalization.py index b946fbbd5..4a2893562 100644 --- a/app/input_normalization.py +++ b/app/input_normalization.py @@ -10,8 +10,8 @@ # Set the HGVS_SEQREPO_URL env var so the hgvs library will use the local `utilities/seqfetcher` endpoint instead of # making NCBI API calls. -port = os.environ('PORT', 5000) # The localhost debugger starts the app on port 5000 -os.environ('HGVS_SEQREPO_URL') = f"http://127.0.0.1:{port}/utilities/seqfetcher" +port = os.environ.get('PORT', 5000) # The localhost debugger starts the app on port 5000 +os.environ['HGVS_SEQREPO_URL'] = f"http://127.0.0.1:{port}/utilities/seqfetcher" database_schema = os.environ['UTA_DATABASE_SCHEMA'] # Use the biocommons UTA database if we don't specify a custom one. diff --git a/requirements.txt b/requirements.txt index 6216fe975..e05a82df1 100644 --- a/requirements.txt +++ b/requirements.txt @@ -9,15 +9,10 @@ pyfastx==2.0.1 pymongo==4.2.0 pyranges==0.0.120 python_dateutil==2.8.2 -requests==2.28.1 +requests==2.32.5 streamlit==1.19.0 -<<<<<<< HEAD streamlit-aggrid==1.0.5 pytest==7.1.1 -======= -pytest==7.1.1 -deepdiff==6.7.1 ->>>>>>> 1b4630a (Add SPDI normalization with refseq files) python-dotenv==1.0.1 toml==0.10.2 pandas==2.2.2 @@ -25,9 +20,6 @@ numpy==2.0.2 pyliftover==0.4 hgvs==1.5.4 py-ard==1.5.3 -<<<<<<< HEAD openpyxl==3.1.5 pyvcf3==1.0.3 tqdm==4.65.1 -======= ->>>>>>> 1b4630a (Add SPDI normalization with refseq files) diff --git a/utilities/common.py b/utilities/common.py index 207e83f23..ae65dda90 100644 --- a/utilities/common.py +++ b/utilities/common.py @@ -11,7 +11,7 @@ PROJECT_ROOT = BASE_DIR.parent load_dotenv(PROJECT_ROOT / "secrets.env") -utilities_data_client_uri = f"mongodb+srv://readonly:{os.environ('MONGODB_READONLY_PASSWORD')}@cluster0.8ianr.mongodb.net/UtilitiesData" +utilities_data_client_uri = f"mongodb+srv://readonly:{os.environ['MONGODB_READONLY_PASSWORD']}@cluster0.8ianr.mongodb.net/UtilitiesData" utilities_client = pymongo.MongoClient(utilities_data_client_uri) utilities_db = utilities_client.UtilitiesData transcript_data = utilities_db.Transcripts diff --git a/utilities/vcf2json.py b/utilities/vcf2json.py index ac3d9fe34..0f59676a0 100644 --- a/utilities/vcf2json.py +++ b/utilities/vcf2json.py @@ -7,9 +7,8 @@ import common import pyfastx import vcf - -from .gene_ref_seq import get_ref_seq_by_chrom -from .spdi import normalize +from gene_ref_seq import get_ref_seq_by_chrom +from spdi import normalize # Fasta file handles cache fasta_cache = {} @@ -34,7 +33,7 @@ def get_fasta(file): def normalize_spdi(ref_seq, pos, ref, alt, build): fasta = get_fasta(BUILD37_FILE) if build == 'GRCh37' else get_fasta(BUILD38_FILE) - return normalize(fasta[ref_seq], pos, ref, alt) + return normalize(fasta[ref_seq], ref_seq, pos, ref, alt) def add_phase_records(record, phased_rec_map, sample_position):