<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 How to associate consequences derived from HGVS expressions with ClinVar records?
Step1:  Part 1
Step2:  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
<ASSISTANT_TASK:>
Python 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
Explanation: How to associate consequences derived from HGVS expressions with ClinVar records?
End of explanation
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
Explanation: Part 1: "Problem" RCVs
RCVs with and without coordinates that are associated with the same genomic HGVS.
End of explanation
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
Explanation: 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" RCVs
What happens if we actually run the two relevant consequence mapping pipelines on these RCVs?
End of explanation
<END_TASK>