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/gene2transcripts.py b/VariantValidator/modules/gene2transcripts.py index 1acfcf0c..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 @@ -283,6 +291,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 +304,20 @@ 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) + # 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 + ) + + logger.info(f"Select Transcripts: {sel_tx_lst} retained transcripts {kept_tx}") + tx_for_gene = kept_tx for line in tx_for_gene: @@ -307,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: + 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: - gen_span = False + 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.info(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, @@ -472,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: @@ -534,6 +545,109 @@ 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. + + Examples: + NM_000088.4 -> ('NM_000088', 4) + ENST00000225964.10 -> ('ENST00000225964', 10) + """ + + base, version = accession.rsplit(".", 1) + + 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"] + + # 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 + + + + +def filter_latest_transcripts(rows): + """ + Remove 'blat' rows and keep only latest transcript versions. + """ + + # 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/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/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/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/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_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' diff --git a/tests/test_inputs.py b/tests/test_inputs.py index 5530dff1..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' @@ -31304,7 +31328,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 +31377,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 #