From c6e56a75e53e4cc0f24f699b90ddbae9927e1c5e Mon Sep 17 00:00:00 2001 From: Peter-J-Freeman Date: Mon, 8 Jun 2026 11:54:10 +0100 Subject: [PATCH 1/5] feat: Align gene to transcript transcript filter Old version does not discriminate between raw and all in select_transcripts field. This does not match VV anf VF processing and leads to slower responses when only latest transcript versions are required. Now, raw returns all transcripts at all versions and all returns all transcripts at latest versions. Closes issue https://github.com/openvar/variantValidator/issues/804 --- VariantValidator/modules/gene2transcripts.py | 68 +++++++++++++++++++- tests/test_gene2transcript.py | 5 ++ 2 files changed, 71 insertions(+), 2 deletions(-) diff --git a/VariantValidator/modules/gene2transcripts.py b/VariantValidator/modules/gene2transcripts.py index 1acfcf0c..e132e02d 100644 --- a/VariantValidator/modules/gene2transcripts.py +++ b/VariantValidator/modules/gene2transcripts.py @@ -283,6 +283,9 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel annotation = g2t.db.get_transcript_annotation(tx[3]) if tx[3] in sel_tx_lst: kept_tx.append(tx) + + # The syntax if x in y is preferred here as it will prevent cases of ["NM_12345.6", "mane_select"] from + # causing issues by defaulting to mane_select (or others in the lists below) elif "mane_select" in sel_tx_lst: if '"mane_select": true' in annotation: kept_tx.append(tx) @@ -293,9 +296,15 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel if '"mane_select": true' in annotation or '"refseq_select": true' in annotation \ or '"ensembl_select": true' in annotation: kept_tx.append(tx) - elif "all" in sel_tx_lst or None in sel_tx_lst: + elif "all" in sel_tx_lst or None in sel_tx_lst or "raw" in sel_tx_lst: kept_tx.append(tx) + if "all" in sel_tx_lst: + logger.info("Set filter to all") + kept_tx = filter_latest_transcripts(kept_tx) + + logger.info(f"Select Transcripts: {sel_tx_lst} retained transcripts {kept_tx}") + tx_for_gene = kept_tx for line in tx_for_gene: @@ -412,7 +421,7 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel error = 'Currently unable to update gene_ids or transcript information records because ' \ 'VariantValidator %s' % str(e) # my_variant.warnings.append(error) - logger.info(error) + logger.warning(error) tx_description = g2t.db.get_transcript_description(tx) # Get annotation @@ -534,6 +543,61 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel return g2d_data + +def get_accession_parts(accession): + """ + Split transcript accession into base accession and version. + + This function is different from the Validator object transcript_filter as it filters the return of tx_for gene and + not tx_for_region which has a different list structure. + + Examples: + NM_000088.4 -> ('NM_000088', 4) + ENST00000225964.10 -> ('ENST00000225964', 10) + ENST00000486572.1/GRCh38 -> ('ENST00000486572', 1) + """ + + # remove optional genome build suffix + accession = accession.split("/")[0] + + # split version + base, version = accession.rsplit(".", 1) + + return base, int(version) + + +def filter_latest_transcripts(rows): + """ + Remove 'blat' rows and keep only latest transcript versions. + """ + + # remove blat rows + rows = [row for row in rows if row[5] != "blat"] + + # find highest version per accession + latest_versions = {} + + for row in rows: + base, version = get_accession_parts(row[3]) + + if ( + base not in latest_versions + or version > latest_versions[base] + ): + latest_versions[base] = version + + # keep only latest versions + filtered_rows = [] + + for row in rows: + base, version = get_accession_parts(row[3]) + + if version == latest_versions[base]: + filtered_rows.append(row) + + return filtered_rows + + def lovd_syntax_check_g2t(query, lovd_syntax_check): # Get additional warnings check_with_lovd = lovd_api.lovd_syntax_check(query, do_lovd_check=lovd_syntax_check, is_a_gene=True) diff --git a/tests/test_gene2transcript.py b/tests/test_gene2transcript.py index 78a0fc93..d0bfc805 100644 --- a/tests/test_gene2transcript.py +++ b/tests/test_gene2transcript.py @@ -169,6 +169,9 @@ def test_select_transcripts(self): output_d = self.vv.gene2transcripts('BRCA2', select_transcripts='flibble') output_e = self.vv.gene2transcripts('BRAF', select_transcripts='mane_select') output_f = self.vv.gene2transcripts('BRAF', select_transcripts='mane') + output_g = self.vv.gene2transcripts('BRCA2', select_transcripts='all') + output_h = self.vv.gene2transcripts('BRCA2', select_transcripts='raw') + print(output) assert len(output['transcripts']) >= 3 assert len(output_b['transcripts']) == 1 @@ -176,6 +179,8 @@ def test_select_transcripts(self): assert len(output_d['transcripts']) == 0 assert len(output_e['transcripts']) == 1 assert len(output_f['transcripts']) >= 2 + assert (len(output_b['transcripts']) < len(output_c['transcripts']) < len(output_g['transcripts']) + < len(output_h['transcripts'])) def test_symbol_valid_hgnc_id(self): symbol = 'HGNC:2197' From e02f93e7eb5b9d46f9700ac1346743ad6cc09210 Mon Sep 17 00:00:00 2001 From: Peter-J-Freeman Date: Mon, 8 Jun 2026 23:40:38 +0100 Subject: [PATCH 2/5] bug: Start list of bad alignments in transcript_map_data.py Certain alignments are too bad to justify keeping. The list in the file transcript_map_data.py can be used to filter mapping options, even if VVTA is updated and new rogue alignments are identified. --- VariantValidator/modules/gapped_mapping.py | 17 ++++++++++--- VariantValidator/modules/mappers.py | 4 +-- .../modules/transcript_map_data.py | 17 +++++++++++++ VariantValidator/modules/vvMixinConverters.py | 11 +++++++- tests/test_inputs.py | 25 ++++++++++++++----- 5 files changed, 62 insertions(+), 12 deletions(-) diff --git a/VariantValidator/modules/gapped_mapping.py b/VariantValidator/modules/gapped_mapping.py index 701b9c3e..fd06225b 100644 --- a/VariantValidator/modules/gapped_mapping.py +++ b/VariantValidator/modules/gapped_mapping.py @@ -1864,6 +1864,9 @@ def g_to_t_compensation(self, ori, hgvs_coding, rec_var): hgvs_genomic.posedit.pos.end.base = start hgvs_genomic = self.variant.hn.normalize(hgvs_genomic) + logger.info(f"g_to_t_compensation returning hgvs_genomic {hgvs_genomic }hgvs_coding {hgvs_coding}, " + f"suppress c normalization {suppress_c_normalization}") + return hgvs_genomic, suppress_c_normalization, hgvs_coding def g_to_t_gapped_mapping_stage2(self, ori, hgvs_coding, hgvs_genomic): @@ -2002,6 +2005,7 @@ def g_to_t_gapped_mapping_stage2(self, ori, hgvs_coding, hgvs_genomic): else: hgvs_coding = copy.deepcopy(hgvs_refreshed_variant) + logger.info(f"g_to_t_gaped_mapping_stage2 returning hgvs_coding {hgvs_coding}") return hgvs_coding def g_to_t_gap_compensation_version3(self, hgvs_alt_genomic, hgvs_coding, ori, alt_chr, rec_var): @@ -2009,7 +2013,9 @@ def g_to_t_gap_compensation_version3(self, hgvs_alt_genomic, hgvs_coding, ori, a self.orientation = int(ori[0]['alt_strand']) hgvs_genomic = copy.deepcopy(hgvs_alt_genomic) - logger.debug('g_to_t gap code 3 active') + logger.debug(f"g_to_t_gap_compensation_version3 triggered with hgvs_alt_genomic {hgvs_alt_genomic}," + f" hgvs_coding {hgvs_coding}, alt_chr {alt_chr}, rec_var {rec_var}, ori {ori}") + rn_hgvs_genomic = self.variant.reverse_normalizer.normalize(hgvs_alt_genomic) self.hgvs_genomic_possibilities.append([rn_hgvs_genomic, ['false', 'false']]) if self.orientation != -1: @@ -2712,6 +2718,8 @@ def g_to_t_gap_compensation_version3(self, hgvs_alt_genomic, hgvs_coding, ori, a except UnboundLocalError: pass + logger.debug(f"g_to_t_gap_compensation_version3 returning hgvs_genomic {hgvs_genomic}, " + f"hgvs_coding {hgvs_coding}") return hgvs_alt_genomic, hgvs_coding def dup_ins_5prime_shift(self, stored_hgvs_not_delins, saved_hgvs_coding): @@ -2984,8 +2992,11 @@ def transcript_disparity(self, reverse_normalized_hgvs_genomic, stored_hgvs_not_ c_tx_hgvs_not_delins = self.validator.vm.n_to_c(self.tx_hgvs_not_delins) except Exception: c_tx_hgvs_not_delins = copy.copy(self.tx_hgvs_not_delins) - genomic_gap_fill_variant_alt = self.validator.vm.t_to_g(c_tx_hgvs_not_delins, self.hgvs_genomic_5pr.ac, - alt_aln_method=self.validator.alt_aln_method) + genomic_gap_fill_variant_alt = self.validator.myvm_t_to_g(c_tx_hgvs_not_delins, + self.hgvs_genomic_5pr.ac, + self.variant.no_norm_evm, + self.variant.hn, + self.variant.map_dat) # Ensure an ALT exists try: diff --git a/VariantValidator/modules/mappers.py b/VariantValidator/modules/mappers.py index b2715c01..437a1db0 100644 --- a/VariantValidator/modules/mappers.py +++ b/VariantValidator/modules/mappers.py @@ -1040,12 +1040,12 @@ def final_tx_to_multiple_genomic(variant, validator, tx_variant, liftover_level= # Loop out gap code under these circumstances! if variant.map_dat.is_gapped_map(variant.hgvs_coding.ac,hgvs_alt_genomic.ac,validator): # warn on gap_compensation for - logger.debug("gap_compensation_3 done for %s when mapped to %s" % - (variant.hgvs_coding.ac, hgvs_alt_genomic.ac)) gap_mapper = gapped_mapping.GapMapper(variant, validator) hgvs_alt_genomic, hgvs_coding = gap_mapper.g_to_t_gap_compensation_version3( hgvs_alt_genomic, variant.hgvs_coding, ori, alt_chr, rec_var) + logger.info(f"gap_compensation_3 done for {variant.hgvs_coding} mapped to {hgvs_alt_genomic}") variant.hgvs_coding = hgvs_coding + logger.info(f"hgvs_coding updated to {hgvs_coding}") # Check for mismatched sequence in dup variants if hgvs_alt_genomic.posedit.edit.type == hgvs_coding.posedit.edit.type and \ diff --git a/VariantValidator/modules/transcript_map_data.py b/VariantValidator/modules/transcript_map_data.py index 58e33450..1280b202 100644 --- a/VariantValidator/modules/transcript_map_data.py +++ b/VariantValidator/modules/transcript_map_data.py @@ -1,4 +1,8 @@ import copy +import logging + +logger = logging.getLogger(__name__) + """ A module for holding the TranscriptMapData specific to an individual transcript in order to avoid having to repeatedly re-pull from the database. @@ -22,6 +26,9 @@ def __init__(self,hdp=None): self.mapped_strands = {} # made using above map data, for dict fetch self.mapping_types = {} # the mapping type for each tx->alt map self.exon_data = {} # the set of exon data per mapping + self._known_bad_alignments = [ + ["NM_001009944.3", "NT_187607.1"] + ] def mapping_options(self,tx_ac,hdp=None): """ @@ -45,6 +52,14 @@ def mapping_options(self,tx_ac,hdp=None): "provider (hdp) for use as a data source") self.mapping_opts[tx_ac] = cur_hdp.get_tx_mapping_options( tx_ac,gap_warn=True) + + self.mapping_opts[tx_ac] = [ + option + for option in self.mapping_opts[tx_ac] + if [option[0], option[1]] not in self._known_bad_alignments + ] + + logger.info(f"Mapping options for {tx_ac}: {self.mapping_opts[tx_ac]}") return self.mapping_opts[tx_ac] def map_strand(self,tx_ac,alt_ac,hdp=None): @@ -157,7 +172,9 @@ def tx_exons(self, tx_ac, alt_ac, alt_aln_method, hdp=None): # If on the reverse strand, reverse the order of elements if tx_exons[0]['alt_strand'] == -1: tx_exons = tx_exons[::-1] + logger.debug(f"Exon data for {tx_ac}: {tx_exons}") return tx_exons else: + logger.debug(f"Exon data for {tx_ac}: {tx_exons}") return tx_exons diff --git a/VariantValidator/modules/vvMixinConverters.py b/VariantValidator/modules/vvMixinConverters.py index 2a84bd24..08a2821d 100644 --- a/VariantValidator/modules/vvMixinConverters.py +++ b/VariantValidator/modules/vvMixinConverters.py @@ -1096,7 +1096,7 @@ def myvm_t_to_g(self, hgvs_c, alt_chr, no_norm_evm, hn, map_dat): hgvs_genomic.posedit.edit.alt = hgvs_genomic.posedit.edit.ref if hgvs_genomic.posedit.edit.type == 'ins' and utilise_gap_code is True: try: - pre_norm_genomic = copy.copy(hgvs_genomic)# can move ins variants (and in doing so break mid base == original bases assumption) + pre_norm_genomic = copy.copy(hgvs_genomic) # can move ins variants (and in doing so break mid base == original bases assumption) hgvs_genomic = hn.normalize(hgvs_genomic) if stored_hgvs_c.posedit.edit.alt and len(stored_hgvs_c.posedit.edit.alt) + 2 == \ len(hgvs_c.posedit.edit.alt) and hgvs_c.posedit.edit.alt == pre_norm_genomic.posedit.edit.alt: @@ -1119,6 +1119,14 @@ def myvm_t_to_g(self, hgvs_c, alt_chr, no_norm_evm, hn, map_dat): hgvs_genomic.posedit.pos.end.base = start hgvs_genomic = hn.normalize(hgvs_genomic) + except AttributeError as e: + if "'Dup' object has no attribute 'alt'" in str(e): + logger.error( + f"Code triggered previously in very poor alignment so not able to fully test, refer to " + f"test_inputs.py tests test_alt_gapping_bug: " + f"hgvs_genomic: {hgvs_genomic}, stored_hgvs_c: {stored_hgvs_c}") + raise + # Statements required to reformat the stored_hgvs_c into a useable synonym if (stored_hgvs_c.posedit.edit.ref == '' or stored_hgvs_c.posedit.edit.ref is None) and expand_out: if stored_hgvs_c.type == 'c': @@ -1396,6 +1404,7 @@ def myvm_t_to_g(self, hgvs_c, alt_chr, no_norm_evm, hn, map_dat): hgvs_genomic = hn.normalize(hgvs_genomic) except: pass + # Correct expansion ref + 2 elif expand_out and ( len(hgvs_genomic.posedit.edit.ref) == (len(stored_hgvs_c.posedit.edit.ref) + 2)): # >= 3: diff --git a/tests/test_inputs.py b/tests/test_inputs.py index 5530dff1..415ae76b 100644 --- a/tests/test_inputs.py +++ b/tests/test_inputs.py @@ -31304,7 +31304,7 @@ def test_issue_733a(self): results = self.vv.validate('chr11:118650341:C:T', 'GRCh37','NM_004397.6').format_as_dict(test=True) assert 'NM_004397.6:c.369G>A' in results - def polyadenylation_a(self): + def test_polyadenylation_a(self): # Test that it fails for genome mismatch results = self.vv.validate('NM_001424184.1:c.438G>A', 'GRCh38', 'all').format_as_dict(test=True) assert 'NM_001424184.1:c.438G>A' in results @@ -31353,17 +31353,30 @@ def test_issue_801(self): assert "NM_014249.4:c.*557del" in results.keys() assert "NM_001281446.1:c.*557del" in results.keys() - def issue_818(self): + def test_issue_818(self): results = self.vv.validate('NC_000022.10:g.19929250_19929251insCCCCGCC', 'GRCh38', 'mane_select', liftover_level=True).format_as_dict(test=True) assert "NM_006440.5:c.70_76dup" in results.keys() assert results["NM_006440.5:c.70_76dup"][ "hgvs_predicted_protein_consequence"] == { - "lrg_slr": "LRG_417p1:p.V26Gfs*132", - "lrg_tlr": "LRG_417p1:p.Val26GlyfsTer132", - "slr": "NP_006431.2:p.V26Gfs*132", - "tlr": "NP_006431.2:p.Val26GlyfsTer132" + "slr": "NP_006431.2:p.(V26Gfs*132)", + "tlr": "NP_006431.2:p.(Val26GlyfsTer132)" } + def test_alt_gapping_bug(self): + results = self.vv.validate('chr16:2089739:G:C', 'GRCh38', 'mane_select', liftover_level=True).format_as_dict(test=True) + assert "NM_001009944.3:c.12900C>G" in results.keys() + assert "NT_187607.1:g.705079dup" not in str(results["NM_001009944.3:c.12900C>G"]["alt_genomic_loci"]) + + def test_alt_gapping_bug_b(self): + results = self.vv.validate('chr16:2089760:C:G', 'GRCh38', 'mane_select', liftover_level=True).format_as_dict(test=True) + assert "NM_001009944.3:c.12879G>C" in results.keys() + assert "NT_187607.1:g.705077_705079dup" not in str(results["NM_001009944.3:c.12879G>C"]["alt_genomic_loci"]) + + def test_alt_gapping_bug_c(self): + results = self.vv.validate('chr16:2090285:C:CT', 'GRCh38', 'mane_select', liftover_level=True).format_as_dict(test=True) + assert "NM_001009944.3:c.12443dup" in results.keys() + assert "NT_187607.1:g.713119_713120insTT" not in str(results["NM_001009944.3:c.12443dup"]["alt_genomic_loci"]) + # # Copyright (C) 2016-2026 VariantValidator Contributors # From 0d9fee463edcf6a650c413a4bf3023b5734dc5d5 Mon Sep 17 00:00:00 2001 From: Peter-J-Freeman Date: Tue, 9 Jun 2026 15:54:59 +0100 Subject: [PATCH 3/5] bug: Adapt g2t transcript filtering for ENS/GRC Ensembl transcripts with the same ID but requiring genome build info require additional handling --- VariantValidator/modules/gene2transcripts.py | 46 ++++++++++++-------- 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/VariantValidator/modules/gene2transcripts.py b/VariantValidator/modules/gene2transcripts.py index e132e02d..69840713 100644 --- a/VariantValidator/modules/gene2transcripts.py +++ b/VariantValidator/modules/gene2transcripts.py @@ -546,24 +546,28 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel def get_accession_parts(accession): """ - Split transcript accession into base accession and version. - - This function is different from the Validator object transcript_filter as it filters the return of tx_for gene and - not tx_for_region which has a different list structure. + Split transcript accession into accession, version, and genome build. Examples: - NM_000088.4 -> ('NM_000088', 4) - ENST00000225964.10 -> ('ENST00000225964', 10) - ENST00000486572.1/GRCh38 -> ('ENST00000486572', 1) + NM_000088.4 + -> ('NM_000088', 4, None) + + ENST00000225964.10 + -> ('ENST00000225964', 10, None) + + ENST00000486572.1/GRCh38 + -> ('ENST00000486572', 1, 'GRCh38') """ - # remove optional genome build suffix - accession = accession.split("/")[0] + accession_part, _, genome_build = accession.partition("/") - # split version - base, version = accession.rsplit(".", 1) + base, _, version = accession_part.rpartition(".") - return base, int(version) + return ( + base, + int(version), + genome_build or None, + ) def filter_latest_transcripts(rows): @@ -574,25 +578,29 @@ def filter_latest_transcripts(rows): # remove blat rows rows = [row for row in rows if row[5] != "blat"] - # find highest version per accession + # find highest version per accession/build latest_versions = {} for row in rows: - base, version = get_accession_parts(row[3]) + base, version, genome_build = get_accession_parts(row[3]) + + key = (base, genome_build) if ( - base not in latest_versions - or version > latest_versions[base] + key not in latest_versions + or version > latest_versions[key] ): - latest_versions[base] = version + latest_versions[key] = version # keep only latest versions filtered_rows = [] for row in rows: - base, version = get_accession_parts(row[3]) + base, version, genome_build = get_accession_parts(row[3]) + + key = (base, genome_build) - if version == latest_versions[base]: + if version == latest_versions[key]: filtered_rows.append(row) return filtered_rows From 53cf22d05f2d544ec2503fe1c0797dbf26a6a7e1 Mon Sep 17 00:00:00 2001 From: Peter-J-Freeman Date: Tue, 9 Jun 2026 17:31:27 +0100 Subject: [PATCH 4/5] bug: central filter for hdp.get_tx_for_gene hdp.get_tx_for_gene has specific structures that need to be filtered out. This can be plugged into all other uses of hdp.get_tx_for_gene and adapted for other similar lists from hdp --- VariantValidator/modules/gene2transcripts.py | 460 ++++++++++--------- 1 file changed, 251 insertions(+), 209 deletions(-) diff --git a/VariantValidator/modules/gene2transcripts.py b/VariantValidator/modules/gene2transcripts.py index 69840713..62f93bb7 100644 --- a/VariantValidator/modules/gene2transcripts.py +++ b/VariantValidator/modules/gene2transcripts.py @@ -14,8 +14,15 @@ logger = logging.getLogger(__name__) -def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, select_transcripts=None, - transcript_set=None, genome_build=None, bypass_genomic_spans=False, lovd_syntax_check=False): +def gene2transcripts(g2t, + query, + validator=False, + bypass_web_searches=False, + select_transcripts=None, + transcript_set="refseq", + genome_build="GRCh38", + bypass_genomic_spans=False, + lovd_syntax_check=False): """ Generates a list of transcript (UTA supported) and transcript names from a gene symbol or RefSeq transcript ID :param g2t: variant object @@ -26,6 +33,7 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel :param transcript_set: String that defines all, refseq or ensembl :param genome_build: String GRCh37 or GRCh38 :param bypass_genomic_spans: bool + :param lovd_syntax_check: bool :return: dictionary of transcript information """ # Set LOVD data @@ -299,9 +307,14 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel elif "all" in sel_tx_lst or None in sel_tx_lst or "raw" in sel_tx_lst: kept_tx.append(tx) + # Clean structures + kept_tx = clean_transcripts(kept_tx, genome_build=genome_build, transcript_set=transcript_set) + if "all" in sel_tx_lst: logger.info("Set filter to all") - kept_tx = filter_latest_transcripts(kept_tx) + kept_tx = filter_latest_transcripts( + kept_tx + ) logger.info(f"Select Transcripts: {sel_tx_lst} retained transcripts {kept_tx}") @@ -316,162 +329,162 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel if re.match("N[MR]_", line[3]): continue - if (line[3].startswith('NM_') or line[3].startswith('NR_') or line[3].startswith('ENST')) and \ - '..' not in line[3] and \ - '_NG_' not in line[3] and \ - "~" not in line[3]: - - # Filter for only requested genome build - if genome_build is not None: - chr_num = seq_data.to_chr_num_refseq(line[4], genome_build) - if chr_num is None: - continue + # Transcript ID + tx = line[3] - # Transcript ID - tx = line[3] + # Protein id + prot_id = g2t.hdp.get_pro_ac_for_tx_ac(tx) - # Protein id - prot_id = g2t.hdp.get_pro_ac_for_tx_ac(tx) - - # Get additional tx_ information - try: - tx_exons = g2t.hdp.get_tx_exons(tx, line[4], line[5]) - except vvhgvs.exceptions.HGVSDataNotAvailableError: - continue - tx_orientation = tx_exons[0]['alt_strand'] - - # Fetch the sequence details, length is the first item - tx_anno = g2t.hdp.get_tx_seq_anno(tx) - tx_len = tx_anno[0] - - # Exon Set - # Collect genomic span for the transcript against known genomic/gene reference sequences - gen_start_pos = None - gen_end_pos = None - exon_set = [] - # get total exons - total_exons = len(tx_exons) - # Set exon counter for current exon + # Get additional tx_ information + try: + tx_exons = g2t.hdp.get_tx_exons(tx, line[4], line[5]) + except vvhgvs.exceptions.HGVSDataNotAvailableError: + continue + tx_orientation = tx_exons[0]['alt_strand'] + + # Fetch the sequence details, length is the first item + tx_anno = g2t.hdp.get_tx_seq_anno(tx) + tx_len = tx_anno[0] + + # Exon Set + # Collect genomic span for the transcript against known genomic/gene reference sequences + gen_start_pos = None + gen_end_pos = None + exon_set = [] + # get total exons + total_exons = len(tx_exons) + # Set exon counter for current exon + if tx_orientation == 1: + current_exon_number = 0 + else: + current_exon_number = total_exons + 1 + for tx_pos in tx_exons: if tx_orientation == 1: - current_exon_number = 0 + current_exon_number = current_exon_number + 1 else: - current_exon_number = total_exons + 1 - for tx_pos in tx_exons: - if tx_orientation == 1: - current_exon_number = current_exon_number + 1 - else: - current_exon_number = current_exon_number - 1 - # Collect the exon_set information - """ - tx_exons have the following attributes:: - { - 'tes_exon_set_id' : 98390 - 'aes_exon_set_id' : 298679 - 'tx_ac' : 'NM_199425.2' - 'alt_ac' : 'NC_000020.10' - 'alt_strand' : -1 - 'alt_aln_method' : 'splign' - 'ord' : 2 - 'tx_exon_id' : 936834 - 'alt_exon_id' : 2999028 - 'tx_start_i' : 786 - 'tx_end_i' : 1196 - 'alt_start_i' : 25059178 - 'alt_end_i' : 25059588 - 'cigar' : '410=' - } - """ - current_exon = {"transcript_start": tx_pos['tx_start_i'] + 1, - "transcript_end": tx_pos['tx_end_i'], - "genomic_start": tx_pos['alt_start_i'] + 1, - "genomic_end": tx_pos['alt_end_i'], - "cigar": tx_pos['cigar'], - "exon_number": current_exon_number - } - exon_set.append(current_exon) - start_pos = tx_pos['alt_start_i'] - end_pos = tx_pos['alt_end_i'] - if gen_start_pos is None: - gen_start_pos = start_pos - else: - if int(start_pos) < int(gen_start_pos): - gen_start_pos = int(start_pos) - if gen_end_pos is None: - gen_end_pos = end_pos - else: - if int(end_pos) > int(gen_end_pos): - gen_end_pos = int(end_pos) - - # reverse the exon_set to maintain gene and not genome orientation if gene is -1 orientated - if tx_orientation == -1: - exon_set.reverse() - if bypass_genomic_spans is True: - gen_span = False - elif ('NG_' in line[4] or 'NC_0' in line[4]) and line[5] != 'blat': - gen_span = True + current_exon_number = current_exon_number - 1 + # Collect the exon_set information + """ + tx_exons have the following attributes:: + { + 'tes_exon_set_id' : 98390 + 'aes_exon_set_id' : 298679 + 'tx_ac' : 'NM_199425.2' + 'alt_ac' : 'NC_000020.10' + 'alt_strand' : -1 + 'alt_aln_method' : 'splign' + 'ord' : 2 + 'tx_exon_id' : 936834 + 'alt_exon_id' : 2999028 + 'tx_start_i' : 786 + 'tx_end_i' : 1196 + 'alt_start_i' : 25059178 + 'alt_end_i' : 25059588 + 'cigar' : '410=' + } + """ + current_exon = {"transcript_start": tx_pos['tx_start_i'] + 1, + "transcript_end": tx_pos['tx_end_i'], + "genomic_start": tx_pos['alt_start_i'] + 1, + "genomic_end": tx_pos['alt_end_i'], + "cigar": tx_pos['cigar'], + "exon_number": current_exon_number + } + exon_set.append(current_exon) + start_pos = tx_pos['alt_start_i'] + end_pos = tx_pos['alt_end_i'] + if gen_start_pos is None: + gen_start_pos = start_pos else: - gen_span = False + if int(start_pos) < int(gen_start_pos): + gen_start_pos = int(start_pos) + if gen_end_pos is None: + gen_end_pos = end_pos + else: + if int(end_pos) > int(gen_end_pos): + gen_end_pos = int(end_pos) + + # reverse the exon_set to maintain gene and not genome orientation if gene is -1 orientated + if tx_orientation == -1: + exon_set.reverse() + if bypass_genomic_spans is True: + gen_span = False + elif ('NG_' in line[4] or 'NC_0' in line[4]) and line[5] != 'blat': + gen_span = True + else: + gen_span = False - tx_description = g2t.db.get_transcript_description(tx) + tx_description = g2t.db.get_transcript_description(tx) - if tx_description == 'none': - try: - g2t.db.update_transcript_info_record(tx, g2t) - except fn.DatabaseConnectionError as e: - error = 'Currently unable to update gene_ids or transcript information records because ' \ - 'VariantValidator %s' % str(e) - # my_variant.warnings.append(error) - logger.warning(error) - tx_description = g2t.db.get_transcript_description(tx) - - # Get annotation + if tx_description == 'none': try: - tx_annotation = g2t.db.get_transcript_annotation(tx) - tx_annotation = json.loads(tx_annotation) - # Missing annotation data - except json.decoder.JSONDecodeError: - continue + g2t.db.update_transcript_info_record(tx, g2t) + except fn.DatabaseConnectionError as e: + error = 'Currently unable to update gene_ids or transcript information records because ' \ + 'VariantValidator %s' % str(e) + # my_variant.warnings.append(error) + logger.warning(error) + tx_description = g2t.db.get_transcript_description(tx) - # Check for duplicates - if tx not in recovered: - recovered.append(tx) - if len(line) >= 3 and isinstance(line[1], int): + # Get annotation + try: + tx_annotation = g2t.db.get_transcript_annotation(tx) + tx_annotation = json.loads(tx_annotation) + # Missing annotation data + except json.decoder.JSONDecodeError: + continue + + # Check for duplicates + if tx not in recovered: + recovered.append(tx) + if len(line) >= 3 and isinstance(line[1], int): + genes_and_tx.append({'reference': tx, + 'description': tx_description, + 'annotations': tx_annotation, + 'translation': prot_id, + 'length': tx_len, + 'coding_start': line[1] + 1, + 'coding_end': line[2], + # 'orientation': tx_orientation, + 'genomic_spans': {} + }) + else: + genes_and_tx.append({'reference': tx, + 'description': tx_description, + 'annotations': tx_annotation, + 'translation': prot_id, + 'length': tx_len, + 'coding_start': None, + 'coding_end': None, + # 'orientation': tx_orientation, + 'genomic_spans': {} + }) + # LRG information + lrg_transcript = g2t.db.get_lrg_transcript_id_from_refseq_transcript_id(tx) + if lrg_transcript != 'none': + if line[1] is None: genes_and_tx.append({'reference': tx, 'description': tx_description, 'annotations': tx_annotation, 'translation': prot_id, 'length': tx_len, - 'coding_start': line[1] + 1, - 'coding_end': line[2], + 'coding_start': None, + 'coding_end': None, # 'orientation': tx_orientation, 'genomic_spans': {} }) - else: - genes_and_tx.append({'reference': tx, + elif sel_tx_lst is False: + genes_and_tx.append({'reference': lrg_transcript, 'description': tx_description, 'annotations': tx_annotation, - 'translation': prot_id, 'length': tx_len, - 'coding_start': None, - 'coding_end': None, - # 'orientation': tx_orientation, + 'translation': lrg_transcript.replace('t', 'p'), + 'coding_start': line[1] + 1, + 'coding_end': line[2], 'genomic_spans': {} }) - # LRG information - lrg_transcript = g2t.db.get_lrg_transcript_id_from_refseq_transcript_id(tx) - if lrg_transcript != 'none': - if line[1] is None: - genes_and_tx.append({'reference': tx, - 'description': tx_description, - 'annotations': tx_annotation, - 'translation': prot_id, - 'length': tx_len, - 'coding_start': None, - 'coding_end': None, - # 'orientation': tx_orientation, - 'genomic_spans': {} - }) - elif sel_tx_lst is False: + else: + if lrg_transcript in sel_tx_lst: genes_and_tx.append({'reference': lrg_transcript, 'description': tx_description, 'annotations': tx_annotation, @@ -481,51 +494,40 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel 'coding_end': line[2], 'genomic_spans': {} }) + + # Add the genomic span information + if gen_span is True: + for check_tx in genes_and_tx: + lrg_transcript = g2t.db.get_lrg_transcript_id_from_refseq_transcript_id(tx) + if check_tx['reference'] == tx: + if gen_start_pos < gen_end_pos: + check_tx['genomic_spans'][line[4]] = {'start_position': gen_start_pos + 1, + 'end_position': gen_end_pos, + 'orientation': tx_orientation, + 'exon_structure': exon_set, + "total_exons": total_exons} else: - if lrg_transcript in sel_tx_lst: - genes_and_tx.append({'reference': lrg_transcript, - 'description': tx_description, - 'annotations': tx_annotation, - 'length': tx_len, - 'translation': lrg_transcript.replace('t', 'p'), - 'coding_start': line[1] + 1, - 'coding_end': line[2], - 'genomic_spans': {} - }) - - # Add the genomic span information - if gen_span is True: - for check_tx in genes_and_tx: - lrg_transcript = g2t.db.get_lrg_transcript_id_from_refseq_transcript_id(tx) - if check_tx['reference'] == tx: - if gen_start_pos < gen_end_pos: - check_tx['genomic_spans'][line[4]] = {'start_position': gen_start_pos + 1, - 'end_position': gen_end_pos, - 'orientation': tx_orientation, - 'exon_structure': exon_set, - "total_exons": total_exons} - else: - check_tx['genomic_spans'][line[4]] = {'start_position': gen_end_pos + 1, - 'end_position': gen_start_pos, - 'orientation': tx_orientation, - 'exon_structure': exon_set, - "total_exons": total_exons} - if lrg_transcript != 'none': - if check_tx['reference'] == lrg_transcript: - if 'NG_' in line[4]: - lrg_id = g2t.db.get_lrg_id_from_refseq_gene_id(line[4]) - if lrg_id[0] in lrg_transcript: - check_tx['genomic_spans'][line[4]] = {'start_position': gen_start_pos + 1, - 'end_position': gen_end_pos, - 'orientation': 1, - 'exon_structure': exon_set, - "total_exons": total_exons} - - check_tx['genomic_spans'][lrg_id[0]] = {'start_position': gen_start_pos + 1, - 'end_position': gen_end_pos, - 'orientation': 1, - 'exon_structure': exon_set, - "total_exons": total_exons} + check_tx['genomic_spans'][line[4]] = {'start_position': gen_end_pos + 1, + 'end_position': gen_start_pos, + 'orientation': tx_orientation, + 'exon_structure': exon_set, + "total_exons": total_exons} + if lrg_transcript != 'none': + if check_tx['reference'] == lrg_transcript: + if 'NG_' in line[4]: + lrg_id = g2t.db.get_lrg_id_from_refseq_gene_id(line[4]) + if lrg_id[0] in lrg_transcript: + check_tx['genomic_spans'][line[4]] = {'start_position': gen_start_pos + 1, + 'end_position': gen_end_pos, + 'orientation': 1, + 'exon_structure': exon_set, + "total_exons": total_exons} + + check_tx['genomic_spans'][lrg_id[0]] = {'start_position': gen_start_pos + 1, + 'end_position': gen_end_pos, + 'orientation': 1, + 'exon_structure': exon_set, + "total_exons": total_exons} # Return data dict if bypass_web_searches is True: @@ -546,28 +548,75 @@ def gene2transcripts(g2t, query, validator=False, bypass_web_searches=False, sel def get_accession_parts(accession): """ - Split transcript accession into accession, version, and genome build. + Split transcript accession into base accession and version. Examples: - NM_000088.4 - -> ('NM_000088', 4, None) + NM_000088.4 -> ('NM_000088', 4) + ENST00000225964.10 -> ('ENST00000225964', 10) + """ - ENST00000225964.10 - -> ('ENST00000225964', 10, None) + base, version = accession.rsplit(".", 1) - ENST00000486572.1/GRCh38 - -> ('ENST00000486572', 1, 'GRCh38') + return base, int(version) + + +def clean_transcripts(rows, genome_build="GRCh38", transcript_set=None): """ + Clean transcripts. + + Accepts the hdp.get_tx_for_gene list format + + :param rows: transcript list from hdp.get_tx_for_gene + :param genome_build: genome build string GRCh37 or GRCh38 + :param transcript_set: transcript set string refseq or ensembl + + return: cleaned transcript list + """ + + # remove blat rows + rows = [row for row in rows if row[5] != "blat"] - accession_part, _, genome_build = accession.partition("/") + # VVTA occasionally contains Ensembl accessions suffixed + # with genome builds (e.g. /GRCh37 or /GRCh38). + # Exclude these duplicated build-specific records. + # It also contains some strange RefSeq characters which + # we also need to filter out + # May need to be modified when we correctly handle such formats currently filtered + + if transcript_set == "ensembl": + rows = [ + row for row in rows + if ( + row[3].startswith("ENST") and + "/" not in row[3] and + ( + genome_build is None or + seq_data.to_chr_num_refseq(row[4], genome_build) is not None + ) + ) + ] + + elif transcript_set == "refseq": + rows = [ + row for row in rows + if ( + ( + row[3].startswith("NM_") or + row[3].startswith("NR_") + ) and + ".." not in row[3] and + "_NG_" not in row[3] and + "~" not in row[3] and + ( + genome_build is None or + seq_data.to_chr_num_refseq(row[4], genome_build) is not None + ) + ) + ] + + return rows - base, _, version = accession_part.rpartition(".") - return ( - base, - int(version), - genome_build or None, - ) def filter_latest_transcripts(rows): @@ -575,32 +624,25 @@ def filter_latest_transcripts(rows): Remove 'blat' rows and keep only latest transcript versions. """ - # remove blat rows - rows = [row for row in rows if row[5] != "blat"] - - # find highest version per accession/build + # find highest version per accession latest_versions = {} for row in rows: - base, version, genome_build = get_accession_parts(row[3]) - - key = (base, genome_build) + base, version = get_accession_parts(row[3]) if ( - key not in latest_versions - or version > latest_versions[key] + base not in latest_versions + or version > latest_versions[base] ): - latest_versions[key] = version + latest_versions[base] = version # keep only latest versions filtered_rows = [] for row in rows: - base, version, genome_build = get_accession_parts(row[3]) - - key = (base, genome_build) + base, version = get_accession_parts(row[3]) - if version == latest_versions[key]: + if version == latest_versions[base]: filtered_rows.append(row) return filtered_rows From 938e826b7cb0ecbd6528a4aa55b8418fc5da40e4 Mon Sep 17 00:00:00 2001 From: Peter-J-Freeman Date: Tue, 9 Jun 2026 17:33:34 +0100 Subject: [PATCH 5/5] feat: Add toggle for shorthand_vcf Allows the API to optionally set the shorthand vcf format with a mod to the interface. The gui will keep the toggle switched on i.e. True. Also, range increased to 500 to allow capture or pretty much all vars in dbSNP --- VariantValidator/modules/variant.py | 1 + VariantValidator/modules/vvMixinCore.py | 11 +++++-- tests/test_inputs.py | 40 ++++++++++++++++++++----- 3 files changed, 42 insertions(+), 10 deletions(-) diff --git a/VariantValidator/modules/variant.py b/VariantValidator/modules/variant.py index 6b04a173..5b882f25 100644 --- a/VariantValidator/modules/variant.py +++ b/VariantValidator/modules/variant.py @@ -23,6 +23,7 @@ def __init__(self, original, quibble=None, warnings=None, write=True, primary_as self.input_parses = None # quibble as hgvs variant object self.transcript_type = None self.lovd_syntax_check = None + self.shorthand_vcf = None self.lovd_messages = None self.lovd_corrections = None diff --git a/VariantValidator/modules/vvMixinCore.py b/VariantValidator/modules/vvMixinCore.py index 9dfcc676..28479f06 100644 --- a/VariantValidator/modules/vvMixinCore.py +++ b/VariantValidator/modules/vvMixinCore.py @@ -49,7 +49,8 @@ def validate(self, select_transcripts, transcript_set=None, liftover_level=False, - lovd_syntax_check=False): + lovd_syntax_check=False, + shorthand_vcf=False): """ This is the main validator function. :param batch_variant: A string containing the variant to be validated @@ -59,6 +60,7 @@ def validate(self, Selecting multiple transcripts will lead to a multiple variant outputs. :param transcript_set: 'refseq' or 'ensembl' :param lovd_syntax_check: True or False + :param shorthand_vcf: True or False :return: """ logger.debug("Running validate with inputs %s and assembly %s", batch_variant, selected_assembly) @@ -78,6 +80,9 @@ def validate(self, self.selected_assembly = selected_assembly self.select_transcripts = select_transcripts + # Set output VCF format + self.shorthand_vcf = shorthand_vcf + # Set LOVD syntax checker self.lovd_syntax_check = lovd_syntax_check @@ -1607,11 +1612,13 @@ def _apply_met_variation(data): else: variant.hgvs_refseqgene_variant = '' hgd = "hgvs_genomic_description" - def _vcf_abrv(hgvs,vcf,max_non_abrv_len=100): + def _vcf_abrv(hgvs,vcf,max_non_abrv_len=500): """ Abbreviate long del/dup/ins type vcf Use start-stop as pos ref N alt as type, not pos as start ref and alt as (extra long seq""" + if self.shorthand_vcf is False: + return vcf try: # only shorten vcf if ref based and long # (and not uncertain which should not have vcf data) diff --git a/tests/test_inputs.py b/tests/test_inputs.py index 415ae76b..e5d7ed99 100644 --- a/tests/test_inputs.py +++ b/tests/test_inputs.py @@ -1159,7 +1159,7 @@ def test_variant20(self): def test_variant21(self): variant = 'NM_000518.4:c.316_*100del' - results = self.vv.validate(variant, 'GRCh37', 'all').format_as_dict(test=True) + results = self.vv.validate(variant, 'GRCh37', 'all', shorthand_vcf=True).format_as_dict(test=True) print(results) assert results['flag'] == 'gene_variant' @@ -1178,13 +1178,37 @@ def test_variant21(self): assert results['NM_000518.4:c.316_*100del']['hgvs_lrg_variant'] == '' self.assertCountEqual(results['NM_000518.4:c.316_*100del']['alt_genomic_loci'], []) assert results['NM_000518.4:c.316_*100del']['primary_assembly_loci']['hg19'] == { - 'hgvs_genomic_description': 'NC_000011.9:g.5246728_5246956del', 'vcf': {'alt': 'DEL', 'chr': 'chr11', 'pos': '5246728-5246956', 'ref': 'N'}} - assert results['NM_000518.4:c.316_*100del']['primary_assembly_loci']['hg38'] == { - 'hgvs_genomic_description': 'NC_000011.10:g.5225498_5225726del', 'vcf': {'alt': 'DEL', 'chr': 'chr11', 'pos': '5225498-5225726', 'ref': 'N'}} - assert results['NM_000518.4:c.316_*100del']['primary_assembly_loci']['grch37'] == { - 'hgvs_genomic_description': 'NC_000011.9:g.5246728_5246956del', 'vcf': {'alt': 'DEL', 'chr': '11', 'pos': '5246728-5246956', 'ref': 'N'}} + "hgvs_genomic_description": "NC_000011.9:g.5246728_5246956del", + "vcf": { + "alt": "A", + "chr": "chr11", + "pos": "5246727", + "ref": "AATCCAGATGCTCAAGGCCCTTCATAATATCCCCCAGTTTAGTAGTTGGACTTAGGGAACAAAGGAACCTTTAATAGAAATTGGACAGCAAGAAAGCGAGCTTAGTGATACTTGTGGGCCAGGGCATTAGCCACACCAGCCACCACTTTCTGATAGGCAGCCTGCACTGGTGGGGTGAATTCTTTGCCAAAGTGATGGGCCAGCACACAGACCAGCACGTTGCCCAGGAG" + }} + assert results['NM_000518.4:c.316_*100del']['primary_assembly_loci']['hg38'] == { + "hgvs_genomic_description": "NC_000011.10:g.5225498_5225726del", + "vcf": { + "alt": "A", + "chr": "chr11", + "pos": "5225497", + "ref": "AATCCAGATGCTCAAGGCCCTTCATAATATCCCCCAGTTTAGTAGTTGGACTTAGGGAACAAAGGAACCTTTAATAGAAATTGGACAGCAAGAAAGCGAGCTTAGTGATACTTGTGGGCCAGGGCATTAGCCACACCAGCCACCACTTTCTGATAGGCAGCCTGCACTGGTGGGGTGAATTCTTTGCCAAAGTGATGGGCCAGCACACAGACCAGCACGTTGCCCAGGAG" + }} + assert results['NM_000518.4:c.316_*100del']['primary_assembly_loci']['grch37'] == { + "hgvs_genomic_description": "NC_000011.9:g.5246728_5246956del", + "vcf": { + "alt": "A", + "chr": "11", + "pos": "5246727", + "ref": "AATCCAGATGCTCAAGGCCCTTCATAATATCCCCCAGTTTAGTAGTTGGACTTAGGGAACAAAGGAACCTTTAATAGAAATTGGACAGCAAGAAAGCGAGCTTAGTGATACTTGTGGGCCAGGGCATTAGCCACACCAGCCACCACTTTCTGATAGGCAGCCTGCACTGGTGGGGTGAATTCTTTGCCAAAGTGATGGGCCAGCACACAGACCAGCACGTTGCCCAGGAG" + }} assert results['NM_000518.4:c.316_*100del']['primary_assembly_loci']['grch38'] == { - 'hgvs_genomic_description': 'NC_000011.10:g.5225498_5225726del', 'vcf': {'alt': 'DEL', 'chr': '11', 'pos': '5225498-5225726', 'ref': 'N'}} + "hgvs_genomic_description": "NC_000011.10:g.5225498_5225726del", + "vcf": { + "alt": "A", + "chr": "11", + "pos": "5225497", + "ref": "AATCCAGATGCTCAAGGCCCTTCATAATATCCCCCAGTTTAGTAGTTGGACTTAGGGAACAAAGGAACCTTTAATAGAAATTGGACAGCAAGAAAGCGAGCTTAGTGATACTTGTGGGCCAGGGCATTAGCCACACCAGCCACCACTTTCTGATAGGCAGCCTGCACTGGTGGGGTGAATTCTTTGCCAAAGTGATGGGCCAGCACACAGACCAGCACGTTGCCCAGGAG" + }} assert results['NM_000518.4:c.316_*100del']['reference_sequence_records'] == { 'transcript': 'https://www.ncbi.nlm.nih.gov/nuccore/NM_000518.4', 'protein': 'https://www.ncbi.nlm.nih.gov/nuccore/NP_000509.1', @@ -30368,7 +30392,7 @@ def test_variant333(self): def test_variant334(self): variant = 'NM_000061.2:c.588_589insCTACATAG' - results = self.vv.validate(variant, 'GRCh37', 'all').format_as_dict(test=True) + results = self.vv.validate(variant, 'GRCh37', 'all', shorthand_vcf=True).format_as_dict(test=True) print(results) assert results['flag'] == 'gene_variant'