diff --git a/VariantValidator/modules/format_converters.py b/VariantValidator/modules/format_converters.py index d3426a5c..42d352c7 100644 --- a/VariantValidator/modules/format_converters.py +++ b/VariantValidator/modules/format_converters.py @@ -15,10 +15,10 @@ logger = logging.getLogger(__name__) -def initial_format_conversions(variant, validator, select_transcripts_dict_plus_version): +def initial_format_conversions(variant, validator, select_transcripts_dict_plus_version, batch_list): # VCF type 1 - toskip = vcf2hgvs_stage1(variant, validator) + toskip = vcf2hgvs_stage1(variant, batch_list) if toskip: return True @@ -28,7 +28,7 @@ def initial_format_conversions(variant, validator, select_transcripts_dict_plus_ if toskip: return True - toskip = gene_symbol_catch(variant, validator, select_transcripts_dict_plus_version) + toskip = gene_symbol_catch(variant, validator, select_transcripts_dict_plus_version, batch_list) if toskip: return True @@ -38,13 +38,13 @@ def initial_format_conversions(variant, validator, select_transcripts_dict_plus_ return True # Find not_sub type in input e.g. GGGG>G - toskip = vcf2hgvs_stage4(variant, validator) + toskip = vcf2hgvs_stage4(variant, batch_list) if toskip: return True # Extract variants from HGVS allele descriptions # http://varnomen.hgvs.org/recommendations/DNA/variant/alleles/ - toskip = allele_parser(variant, validator, validator) + toskip = allele_parser(variant, validator, validator, batch_list) if toskip: return True @@ -153,7 +153,7 @@ def final_hgvs_convert(variant,validator): return False -def vcf2hgvs_stage1(variant, validator): +def vcf2hgvs_stage1(variant, batch_list): """ VCF2HGVS stage 1. converts chr-pos-ref-alt into chr:posRef>Alt The output format is a common mistake caused by inaccurate conversion of @@ -218,8 +218,8 @@ def vcf2hgvs_stage1(variant, validator): primary_assembly=variant.primary_assembly, order=variant.order) query_b = Variant(variant.original, quibble=input_b, warnings=variant.warnings, primary_assembly=variant.primary_assembly, order=variant.order) - validator.batch_list.append(query_a) - validator.batch_list.append(query_b) + batch_list.append(query_a) + batch_list.append(query_b) logger.info("Submitting new variant with format %s", input_a) skipvar = True elif vcf_data[3]: @@ -392,7 +392,7 @@ def vcf2hgvs_stage2(variant, validator): return skipvar -def gene_symbol_catch(variant, validator, select_transcripts_dict_plus_version): +def gene_symbol_catch(variant, validator, select_transcripts_dict_plus_version, batch_list): """ Searches for gene symbols that have been used as reference sequence identifiers. Provides a sufficiently repremanding warning, but also provides @@ -450,7 +450,7 @@ def gene_symbol_catch(variant, validator, select_transcripts_dict_plus_version): query = Variant(variant.original, quibble=refreshed_description, warnings=variant.warnings, primary_assembly=variant.primary_assembly, order=variant.order) - validator.batch_list.append(query) + batch_list.append(query) logger.info('HGVS variant nomenclature does not allow the use of a gene symbol (' + query_a_symbol + ') in place of a valid reference sequence') logger.info("Submitting new variant with format %s", refreshed_description) @@ -557,7 +557,7 @@ def refseq_catch(variant, validator, select_transcripts_dict_plus_version): logger.info('NG_:c.PositionVariation descriptions should not be used unless a transcript ' 'reference sequence has also been provided e.g. NG_(NM_):c.PositionVariation. ' 'Resubmitting corrected version.') - validator.batch_list.append(query) + batch_list.append(query) logger.info("Submitting new variant with format %s", refreshed_description) else: variant.warnings.append('A transcript reference sequence has not been provided e.g. ' @@ -590,7 +590,7 @@ def refseq_catch(variant, validator, select_transcripts_dict_plus_version): return skipvar -def vcf2hgvs_stage4(variant, validator): +def vcf2hgvs_stage4(variant, batch_list): """ VCF2HGVS conversion step 4 has two purposes 1. VCF is frequently inappropriately converted into HGVS like descriptions @@ -636,7 +636,7 @@ def vcf2hgvs_stage4(variant, validator): query = Variant(variant.original, quibble=refreshed_description, warnings=variant.warnings, primary_assembly=variant.primary_assembly, order=variant.order) - validator.batch_list.append(query) + batch_list.append(query) logger.info('Multiple ALT sequences detected. Auto-submitting all possible combinations.') logger.info("Submitting new variant with format %s", refreshed_description) skipvar = True @@ -952,7 +952,7 @@ def remap_intronic(hgvs_transy, hgvs_genomic, variant, validator): except AttributeError: pass -def allele_parser(variant, validation, validator): +def allele_parser(variant, validation, validator, batch_list): """ HGVS allele string parsing function Occurance #1 Takes a single HGVS allele description and separates each allele into a @@ -1049,7 +1049,7 @@ def allele_parser(variant, validation, validator): for allele in alleles: query = Variant(variant.original, quibble=allele, warnings=variant.warnings, write=True, primary_assembly=variant.primary_assembly, order=variant.order) - validation.batch_list.append(query) + batch_list.append(query) logger.info("Submitting new variant with format %s", allele) variant.write = False return True diff --git a/VariantValidator/modules/gapped_mapping.py b/VariantValidator/modules/gapped_mapping.py index fd06225b..9a9b0572 100644 --- a/VariantValidator/modules/gapped_mapping.py +++ b/VariantValidator/modules/gapped_mapping.py @@ -8,7 +8,6 @@ from VariantValidator.modules.hgvs_utils import hgvs_delins_parts_to_hgvs_obj, hgvs_dup_to_delins from VariantValidator.modules.variant import TranscriptMapData from VariantValidator.modules.utils import simple_dna_revcomp -import traceback logger = logging.getLogger(__name__) diff --git a/VariantValidator/modules/mappers.py b/VariantValidator/modules/mappers.py index 1633e566..44278613 100644 --- a/VariantValidator/modules/mappers.py +++ b/VariantValidator/modules/mappers.py @@ -21,7 +21,7 @@ class MappersError(Exception): class TranscriptMappingError(Exception): pass -def gene_to_transcripts(variant, validator, select_transcripts_dict): +def gene_to_transcripts(variant, validator, select_transcripts_dict, batch_list): logger.info(f"Mapping {variant.hgvs_formatted} to transcripts") g_query = variant.hgvs_formatted # set hdp for exon mapping fetch before first use @@ -170,7 +170,7 @@ def gene_to_transcripts(variant, validator, select_transcripts_dict): query = Variant(variant.original, quibble=genomic_input, warnings=variant.warnings, primary_assembly=variant.primary_assembly, order=variant.order, selected_assembly=variant.selected_assembly) - validator.batch_list.append(query) + batch_list.append(query) logger.info('Submitting new variant with format %s', genomic_input) else: error = 'TranscriptIdentificationWarning: Mapping unavailable for RefSeqGene ' + str(variant.hgvs_formatted) + \ @@ -253,7 +253,7 @@ def gene_to_transcripts(variant, validator, select_transcripts_dict): expanded_repeat=variant.expanded_repeat) # since we already fetched the exon mappings etc and python uses just a pointer for this set map_dat query.map_dat = variant.map_dat - validator.batch_list.append(query) + batch_list.append(query) logger.info("Submitting new variant with format %s", str(c_description)) # Call next description diff --git a/VariantValidator/modules/transcript_map_data.py b/VariantValidator/modules/transcript_map_data.py index 7db95ce2..7ada466b 100644 --- a/VariantValidator/modules/transcript_map_data.py +++ b/VariantValidator/modules/transcript_map_data.py @@ -1,5 +1,6 @@ import copy import logging +import time logger = logging.getLogger(__name__) @@ -141,15 +142,49 @@ def mapped_exons(self,tx_ac,alt_ac,alt_aln_method=None,hdp=None): "provider (hdp) for use as a data source") if tx_ac not in self.exon_data: self.exon_data[tx_ac] = {} + if alt_ac not in self.exon_data[tx_ac]: if alt_aln_method: aln_method = alt_aln_method else: - aln_method = self.map_type(tx_ac,alt_ac) - self.exon_data[tx_ac][alt_ac] = cur_hdp.get_tx_exons( - tx_ac,alt_ac,aln_method) + aln_method = self.map_type(tx_ac, alt_ac) + + max_retries = 3 + retry_delay = 0.2 + + for attempt in range(1, max_retries + 1): + try: + self.exon_data[tx_ac][alt_ac] = cur_hdp.get_tx_exons( + tx_ac, + alt_ac, + aln_method, + ) + break + + except KeyError as e: + logger.warning( + f"Attempt {attempt}/{max_retries}: " + f"Failed get_tx_exons " + f"tx_ac={tx_ac} " + f"alt_ac={alt_ac} " + f"aln_method={aln_method} " + f"key_error={e}" + ) + + if attempt == max_retries: + logger.exception( + f"Failed get_tx_exons after {max_retries} retries " + f"tx_ac={tx_ac} " + f"alt_ac={alt_ac} " + f"aln_method={aln_method}" + ) + raise + + time.sleep(retry_delay) + if alt_ac not in self.exon_data[tx_ac]: return [] + return self.exon_data[tx_ac][alt_ac] def tx_exons(self, tx_ac, alt_ac, alt_aln_method, hdp=None): diff --git a/VariantValidator/modules/vvMixinCore.py b/VariantValidator/modules/vvMixinCore.py index a68a0f1a..746f09be 100644 --- a/VariantValidator/modules/vvMixinCore.py +++ b/VariantValidator/modules/vvMixinCore.py @@ -123,15 +123,16 @@ def validate(self, batch_queries = [batch_variant] if isinstance(batch_queries, int): batch_queries = [str(batch_queries)] + # Turn each variant into a dictionary. The dictionary will be compiled during validation - self.batch_list = [] + batch_list = [] for queries in batch_queries: if isinstance(queries, int): queries = str(queries) queries = str(queries) queries = queries.strip() query = Variant(queries) - self.batch_list.append(query) + batch_list.append(query) logger.info("Submitting variant with format %s", queries) # Create List to carry batch data output @@ -151,8 +152,8 @@ def validate(self, flag : mitochondrial """ - logger.debug("Batch list length " + str(len(self.batch_list))) - for my_variant in self.batch_list: + logger.debug("Batch list length " + str(len(batch_list))) + for my_variant in batch_list: # Create Normalizers my_variant.hn = vvhgvs.normalizer.Normalizer(self.hdp, @@ -389,8 +390,10 @@ def validate(self, my_variant.warnings.append("Reference sequence type o. should only be used for circular " "reference sequences that are not mitochondrial. Instead use m.") try: - toskip = format_converters.initial_format_conversions(my_variant, self, - select_transcripts_dict_plus_version) + toskip = format_converters.initial_format_conversions(my_variant, + self, + select_transcripts_dict_plus_version, + batch_list) except vvhgvs.exceptions.HGVSError as e: # import traceback @@ -689,7 +692,7 @@ def validate(self, # Now start mapping from genome to transcripts if my_variant.reftype == ':g.': try: - toskip = mappers.gene_to_transcripts(my_variant, self, select_transcripts_dict) + toskip = mappers.gene_to_transcripts(my_variant, self, select_transcripts_dict, batch_list) except IndexError: my_variant.output_type_flag = 'warning' error = '%s cannot be validated in the context of genome build %s, ' \ @@ -737,6 +740,8 @@ def validate(self, my_variant.warnings.append(error) exc_type, exc_value, last_traceback = sys.exc_info() logger.error(str(exc_type) + " " + str(exc_value)) + import traceback + traceback.print_exc() continue # Outside the for loop @@ -744,7 +749,7 @@ def validate(self, logger.debug("End of 1st for loop - Finalising formatting") # order the rows - by_order = sorted(self.batch_list, key=lambda x: x.order) + by_order = sorted(batch_list, key=lambda x: x.order) for variant in by_order: ############################################################### # Runtime information and errors at warning and above only!!! # diff --git a/VariantValidator/modules/vvMixinInit.py b/VariantValidator/modules/vvMixinInit.py index c143e64f..1a53d35d 100644 --- a/VariantValidator/modules/vvMixinInit.py +++ b/VariantValidator/modules/vvMixinInit.py @@ -192,7 +192,6 @@ def __init__(self): self.selected_assembly = None self.select_transcripts = None self.alt_aln_method = None - self.batch_list = [] # Create additional normalizers def create_additional_normalizers_and_mappers(self): diff --git a/tests/test_transcript_map_data.py b/tests/test_transcript_map_data.py index 53dd00e4..28ea30d6 100644 --- a/tests/test_transcript_map_data.py +++ b/tests/test_transcript_map_data.py @@ -1,5 +1,6 @@ import pytest from unittest.mock import MagicMock +from unittest.mock import patch from VariantValidator.modules.transcript_map_data import TranscriptMapData @@ -246,6 +247,102 @@ def test_mapped_exons_requires_hdp(): "NM_000001.1", "NC_000001.11" ) + + +@patch("VariantValidator.modules.transcript_map_data.time.sleep") +def test_mapped_exons_keyerror_retry(mock_sleep, mock_hdp): + mock_hdp.get_tx_exons.side_effect = [ + KeyError("temporary failure"), + KeyError("temporary failure"), + [ + {"alt_strand": 1, "exon": 1}, + {"alt_strand": 1, "exon": 2}, + ], + ] + + tmd = TranscriptMapData(mock_hdp) + + result = tmd.mapped_exons( + "NM_000001.1", + "NC_000001.11", + ) + + assert len(result) == 2 + assert mock_hdp.get_tx_exons.call_count == 3 + assert mock_sleep.call_count == 2 + + +@patch("VariantValidator.modules.transcript_map_data.time.sleep") +def test_mapped_exons_keyerror_retry_failure(mock_sleep, mock_hdp): + mock_hdp.get_tx_exons.side_effect = KeyError("permanent failure") + + tmd = TranscriptMapData(mock_hdp) + + with pytest.raises(KeyError): + tmd.mapped_exons( + "NM_000001.1", + "NC_000001.11", + ) + + assert mock_hdp.get_tx_exons.call_count == 3 + assert mock_sleep.call_count == 2 + +def test_mapping_options_with_explicit_hdp(mock_hdp): + tmd = TranscriptMapData() + + result = tmd.mapping_options( + "NM_000001.1", + hdp=mock_hdp, + ) + + assert len(result) == 2 + assert tmd.hdp is mock_hdp + + +def test_mapped_exons_with_explicit_hdp(mock_hdp): + tmd = TranscriptMapData() + + result = tmd.mapped_exons( + "NM_000001.1", + "NC_000001.11", + hdp=mock_hdp, + ) + + assert len(result) == 2 + assert tmd.hdp is mock_hdp + + +def test_map_type_prefers_non_blat(mock_hdp): + mock_hdp.get_tx_mapping_options.return_value = [ + ["NM_000001.1", "NC_000001.11", "splign", True, 1], + ["NM_000001.1", "NC_000001.11", "blat", True, 1], + ] + + tmd = TranscriptMapData(mock_hdp) + + assert ( + tmd.map_type( + "NM_000001.1", + "NC_000001.11", + ) + == "splign" + ) + + +def test_is_gapped_map_cache_after_fetch(mock_hdp): + tmd = TranscriptMapData(mock_hdp) + + tmd.is_gapped_map( + "NM_000001.1", + "NC_000001.11", + ) + + assert ( + tmd.gap_status["NM_000001.1"]["NC_000001.11"] + is True + ) + + # # Copyright (C) 2016-2026 VariantValidator Contributors #