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