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/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/.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..a7b719253 100644
--- a/app/__init__.py
+++ b/app/__init__.py
@@ -1,7 +1,15 @@
+import os
+
import connexion
import flask
+import hgvs
from flask_cors import CORS
-import os
+
+# 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():
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..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(
@@ -831,6 +832,7 @@ 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)
if not variants:
return jsonify({"resourceType": "Parameters"})
@@ -1095,6 +1097,7 @@ 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)
if not variants:
return jsonify({"resourceType": "Parameters"})
diff --git a/app/input_normalization.py b/app/input_normalization.py
new file mode 100644
index 000000000..4a2893562
--- /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.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.
+# 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..e05a82df1 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -9,7 +9,7 @@ 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
streamlit-aggrid==1.0.5
pytest==7.1.1
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..ae65dda90 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..0f59676a0 100644
--- a/utilities/vcf2json.py
+++ b/utilities/vcf2json.py
@@ -1,11 +1,39 @@
-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], ref_seq, pos, ref, alt)
def add_phase_records(record, phased_rec_map, sample_position):
@@ -204,7 +232,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 +242,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)