<jupyter_start><jupyter_text>How to associate consequences derived from HGVS expressions with ClinVar records?<jupyter_code>from collections import defaultdict, Counter
from itertools import zip_longest
import json
import os
import re
import sys
import urllib

import numpy as np
import requests

from eva_cttv_pipeline.clinvar_xml_io.clinvar_xml_io import *

%matplotlib inline
import matplotlib.pyplot as plt
from eva_cttv_pipeline.clinvar_xml_io.clinvar_xml_io.hgvs_variant import HgvsVariant
PROJECT_ROOT = '/home/april/projects/opentargets/'
clinvar_path = os.path.join(PROJECT_ROOT, 'ClinVarFullRelease_00-latest.xml.gz')
# clinvar_path = os.path.join(PROJECT_ROOT, 'clinvar-sample.xml.gz')
dataset = ClinVarDataset(clinvar_path)
def vcv(record):
    if record.measure:
        return find_mandatory_unique_element(record.rcv, './MeasureSet').attrib['Acc']
limit = None

# for each HGVS that structural variants pipeline would process, how many rcvs/vcvs are associated?
#  and would any of them potentially get annotated by the simple VEP pipeline? (assuming repeats override complex)
complex_hgvs_to_complex_rcv = defaultdict(list)
complex_hgvs_to_complex_vcv = defaultdict(list)

complex_hgvs_to_other_rcv = defaultdict(list)
complex_hgvs_to_other_vcv = defaultdict(list)

i = 0
for r in dataset:
    if pipeline.can_process(r):
        complex_hgvs = [h for h in r.measure.current_hgvs if h is not None]
        for h in complex_hgvs:
            complex_hgvs_to_complex_rcv[h].append(r.accession)
            complex_hgvs_to_complex_vcv[h].append(vcv(r))
    else:
        if r.measure and r.measure.current_hgvs:
            other_hgvs = [h for h in r.measure.current_hgvs if h is not None]
            for h in other_hgvs:
                if h in complex_hgvs_to_complex_rcv:
                    complex_hgvs_to_other_rcv[h].append(r.accession)
                    complex_hgvs_to_other_vcv[h].append(vcv(r))
    i += 1
    if limit and i > limit:
        break<jupyter_output><empty_output><jupyter_text>Part 1: "Problem" RCVsRCVs with and without coordinates that are associated with the same genomic HGVS.<jupyter_code>from eva_cttv_pipeline.clinvar_xml_io.clinvar_xml_io.hgvs_variant import SequenceType

problem_rcvs = []
for h in complex_hgvs_to_other_rcv.keys():
    if HgvsVariant(h).sequence_type == SequenceType.GENOMIC:
        problem_rcvs.extend(complex_hgvs_to_other_rcv[h])
        problem_rcvs.extend(complex_hgvs_to_complex_rcv[h])
problem_rcvs = set(problem_rcvs)
problem_rcvs  # includes both complex and "other" rcvs
for r in dataset:
    if r.accession in problem_rcvs:
        print(r.accession)
        print(vcv(r))
        print(r.measure.current_hgvs)
        print(r.measure.vcf_full_coords)
        print('\n=========\n')
for h, vcvs in complex_hgvs_to_complex_vcv.items():
    num_vcvs = len(set(vcvs))
    if num_vcvs > 1 and HgvsVariant(h).sequence_type == SequenceType.GENOMIC:
        print(h)
        print(set(vcvs))
        print('\n========\n')
# for two sets of HGVS identifiers associated with two different VCVs, what's the intersection & set difference?
with_coordinates = {'NM_000080.4:c.1327delG', 'LRG_1254t1:c.1327del', 'LRG_1254:g.9185del', 'NG_028005.1:g.70553del', 'NG_008029.2:g.9185del', 'NC_000017.11:g.4898892del', 'NC_000017.10:g.4802186del', None, None, 'p.Glu443Lysfs*64', 'NP_000071.1:p.Glu443LysfsTer64'}
no_coordinates = {'LRG_1254t1:c.1327del', 'NM_000080.4:c.1327del', 'LRG_1254:g.9185del', 'NG_028005.1:g.70553del', 'NG_008029.2:g.9185del', 'NC_000017.11:g.4898892del', None, 'LRG_1254p1:p.Glu443fs', 'NP_000071.1:p.Glu443fs'}
with_coordinates & no_coordinates
with_coordinates - no_coordinates
no_coordinates - with_coordinates<jupyter_output><empty_output><jupyter_text>Notes* Is the same measure (VCV or HGVS) ever associated with different RCVs?    * Yes - unlike for repeat expansion variants, these structural ones (no coordinates but with HGVS) can have the same measure/VCV associated with multiple RCVs    * hence _cannot_ use RCV to connect records to consequences* Is the same HGVS ever associated with RCVs with and without complete coordinates?    * Yes, surprisingly...    * when this happens it's via different VCVs* If complex HGVS is associated with another RCV that does have complete coordinates, are these aligned?    * Not always exactly...    * e.g. `NC_000017.11:g.4898892del` is associated with an RCV with no coordinates and one with coordinates `17_4898890_TC_T`* Is the same HGVS ever associated with different VCVs?    * Yes, very often actually* Is the same _set_ of HGVS ever associated with different VCVs?    * Not for the above case but probably can't tell for certain (and probably not worth relying on this)Bonus question: how is ClinVar associating HGVS with VCV???  Part 2: Consequences of "problem" RCVsWhat happens if we actually run the two relevant consequence mapping pipelines on these RCVs?<jupyter_code>import pandas as pd
pd.set_option('display.max_colwidth', None)

from consequence_prediction.structural_variants import pipeline as structural_pipeline
from consequence_prediction.vep_mapping_pipeline.consequence_mapping import colon_based_id_to_vep_id, process_variants
problem_path = os.path.join(PROJECT_ROOT, 'complex-events/rcvs_sharing_hgvs.xml.gz')
problem_dataset = ClinVarDataset(problem_path)
# convert VEP pipeline to be more usable...
IUPAC_AMBIGUOUS_SEQUENCE = re.compile(r'[^ACGT]')

def vep_pipeline_main(clinvar_xml):
    variants = []
    for clinvar_record in ClinVarDataset(clinvar_xml):
        if clinvar_record.measure is None or not clinvar_record.measure.has_complete_coordinates:
            continue
        m = clinvar_record.measure
        if IUPAC_AMBIGUOUS_SEQUENCE.search(m.vcf_ref + m.vcf_alt):
            continue
        variants.append(f'{m.chr}:{m.vcf_pos}:{m.vcf_ref}:{m.vcf_alt}')

    variants_to_query = [colon_based_id_to_vep_id(v) for v in variants]
    variant_results = process_variants(variants_to_query)
    variant_data = []
    for variant_id, gene_id, gene_symbol, consequence_term, distance in variant_results:
        variant_data.append((variant_id, '1', gene_id, gene_symbol, consequence_term, distance))
    consequences = pd.DataFrame(variant_data, columns=('VariantID', 'PlaceholderOnes', 'EnsemblGeneID',
                                               'EnsemblGeneName', 'ConsequenceTerm', 'Distance'))
    return consequences
vep_consequences = vep_pipeline_main(problem_path)
vep_consequences
struct_consequences = structural_pipeline.main(problem_path)
# haven't implemented the single base deletion case as it's not a range, but I think we'd get the following
# https://rest.ensembl.org/vep/human/region/NC_000017.11:4898892-4898892:1/DEL?content-type=application/json
struct_consequences = struct_consequences.append(
    pd.DataFrame(
        [['NC_000017.11 4898892 4898892 DEL +', 1, 'ENSG00000108556', 'CHRNE', 'frameshift_variant', 0]],
        columns=('VariantID', 'PlaceholderOnes', 'EnsemblGeneID', 'EnsemblGeneName', 'ConsequenceTerm', 'Distance')
    )
)
struct_consequences<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0913955780875933, 'fraction_numerical': 0.038304464159977465, 'mean_word_length': 4.672523961661342, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 15, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 11}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '255985', 'license': 'Apache-2.0', 'n_tokens_mistral': 2615, 'n_tokens_neox': 2357, 'n_words': 575, 'path': 'data-exploration/complex-events/notebooks/rcv-vcv-hgvs.ipynb', 'repo_name': 'EBIvariation/eva-cttv-pipeline', 'size': '7101'}