<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step2:  Detailed stats of HGVS in ClinVar
Looking only at records with no functional consequences and no complete chr_pos_ref_alt coordinates. Based on June consequence predictions and ClinVar dump.
Table of contents
Sequence types
Variant types
Ranges
Span lengths
Intronic numbering
Summary
Step3:  Sequence types
Refer to documentation here.
Top of page
Step4:  Variant types
Documentation here
Compare stats here - distribution over all of ClinVar, using XML attribute not HGVS
Top of page
Step5:  Ranges
Mostly documented here though note the first case is not uncertain.
g.12345_12678del- no uncertainty
g.(123456_234567)_(345678_456789)del - variable endpoints within a known range
g.(?_234567)_(345678_?)del - variable endpoints within an unknown range
Also split into genomic vs. other sequences - the former we can definitely deal with, not sure about the others.
See below for intronic coordinate ranges...
Top of page
Step7:  Span lengths
Lengths of spans - definite or minimum provided.  For now only compute for genomic sequences, as getting span length for coding/noncoding is a bit more complicated.
Top of page
Step8:  Intronic numbering
Documentation here and particularly the diagram.
Only used for coding and non-coding reference sequences, so not relevant if we focus on genomic.
Top of page
<ASSISTANT_TASK:>
Python Code:
import os
import re
import sys
import numpy as np
from eva_cttv_pipeline.clinvar_xml_utils import *
from eva_cttv_pipeline.clinvar_identifier_parsing import *
%matplotlib inline
import matplotlib.pyplot as plt
PROJECT_ROOT = '/home/april/projects/opentargets/complex-events'
# dump of all records with no functional consequences and no complete coordinates
# uses June consequence pred + ClinVar 6/26/2021
no_consequences_path = os.path.join(PROJECT_ROOT, 'no-conseq_no-coords.xml.gz')
dataset = ClinVarDataset(no_consequences_path)
def count_hgvs(dataset, regex_dict, exclusive=False, limit=None, include_no_hgvs=True, include_none=True):
    
    Counts records in dataset with HGVS matching a collection of regexes.
    Can be exclusive or non-exclusive counts (see below).
    If limit is provided, will count at most that many records (useful for testing).
    
    Notes:
      * records with multiple HGVS expressions need at least one matching a given regex to be counted once
      * can also count measures with no HGVS and ones that match none of the regexes (only if not exclusive)
      * non-exclusive => record has an HGVS expression that matches this regex.
                         "If we do support X, how many records could we get?"
      * exclusive => record _only_ has HGVS expressions that match this regex (out of this collection).
                     "If we don't support X, how many records must we lose?"
        
    
    n = 0
    # just use a dict instead of a counter, so we have a predictable key order
    results = {k: 0 for k in regex_dict}
    if include_no_hgvs:
        results['no hgvs'] = 0
    if not exclusive and include_none:
        results['none'] = 0
    for record in dataset:
        if not record.measure:
            continue
        if not record.measure.hgvs:
            if include_no_hgvs:
                results['no hgvs'] += 1
            continue
        hs = [h for h in record.measure.hgvs if h is not None]
        n += 1
        if limit and n > limit:
            break
        temp_results = {
            k: any(r.match(h) for h in hs)
            for k, r in regex_dict.items()
        }
        any_match = False
        for k in regex_dict:
            if exclusive:
                if temp_results[k] and not any(temp_results[j] for j in regex_dict if j != k):
                    results[k] += 1
            else:
                if temp_results[k]:
                    results[k] += 1
                    any_match = True
        if not exclusive and include_none and not any_match:
            results['none'] += 1
            
    return results
def print_example_matches(dataset, regex_dict, size=1, limit=None, include_none=True):
    
    Like count_hgvs but returns (size) example matches for each regex where possible.
    
    n = 0
    all_matches = {k: [] for k in regex_dict}
    if include_none:
        all_matches['none'] = []
    for record in dataset:
        if not record.measure or not record.measure.hgvs:
            continue
        hs = [h for h in record.measure.hgvs if h is not None]
        n += 1
        if limit and n > limit:
            break
            
        for h in hs:
            any_match = False
            for k, r in regex_dict.items():
                if r.match(h):
                    all_matches[k].append(h)
                    any_match = True
            if not any_match and include_none:
                all_matches['none'].append(h)
    result = {
        k: [v[i] for i in np.random.choice(len(v), size=min(len(v), size), replace=False)] if v else []
        for k, v in all_matches.items()
    }
    for k in result:
        print(k)
        for s in result[k]:
            print(f'  {s}')
        print('\n==========\n')
Explanation: Detailed stats of HGVS in ClinVar
Looking only at records with no functional consequences and no complete chr_pos_ref_alt coordinates. Based on June consequence predictions and ClinVar dump.
Table of contents
Sequence types
Variant types
Ranges
Span lengths
Intronic numbering
Summary
End of explanation
# be more lenient than what we currently in identifier_parsing
# for example this allows things like `chr11` or `LRG_199p1`
sequence_identifier = r'[a-zA-Z0-9_.]+:'
seq_type_dict = {
    'coding': re.compile(sequence_identifier + r'c\.'),
    'genomic': re.compile(sequence_identifier + r'g\.'),
    'non-coding': re.compile(sequence_identifier + r'n\.'),  # transcript but not coding for a protein
    'protein': re.compile(sequence_identifier + r'p\.'),
    'mitochondrial': re.compile(sequence_identifier + r'm\.'),
    'circular': re.compile(sequence_identifier + r'o\.'),
    'RNA': re.compile(sequence_identifier + r'r\.'),
}
print_example_matches(dataset, seq_type_dict, size=5)
seq_type_counts = count_hgvs(dataset, seq_type_dict, exclusive=False)
plt.figure(figsize=(15,7))
plt.title('Sequence type (non-exclusive)')
plt.bar(seq_type_counts.keys(), seq_type_counts.values())
seq_type_counts
# coding or non-coding
2192 + 225
# have hgvs in general
17649 - 4030
seq_type_counts_exclusive = count_hgvs(dataset, seq_type_dict, exclusive=True)
plt.figure(figsize=(15,7))
plt.title('Sequence type (exclusive)')
plt.bar(seq_type_counts_exclusive.keys(), seq_type_counts_exclusive.values())
seq_type_counts_exclusive
Explanation: Sequence types
Refer to documentation here.
Top of page
End of explanation
genomic_sequence = f'^{sequence_identifier}g\.'
all_other_sequence = f'^{sequence_identifier}[a-fh-z]\.'
# double-counts hybrid things, e.g.
#  * NC_000013.9:g.93703239_93802554del99316insCTA
#  * NC_000016.9:g.2155486_2155487ins2145304_2155487inv
variant_regex = {
    'substitution (genomic)': re.compile(f'{genomic_sequence}.*?>.*?'),
    'deletion (genomic)': re.compile(f'{genomic_sequence}.*?del(?!ins).*?'),
    'duplication (genomic)': re.compile(f'{genomic_sequence}.*?dup.*?'),
    'insertion (genomic)': re.compile(f'{genomic_sequence}.*?(?<!del)ins.*?'),
    'inversion (genomic)': re.compile(f'{genomic_sequence}.*?inv.*?'),
    'delins (genomic)': re.compile(f'{genomic_sequence}.*?delins.*?'),
    
    'substitution (other)': re.compile(f'{all_other_sequence}.*?>.*?'),
    'deletion (other)': re.compile(f'{all_other_sequence}.*?del(?!ins).*?'),
    'duplication (other)': re.compile(f'{all_other_sequence}.*?dup.*?'),
    'insertion (other)': re.compile(f'{all_other_sequence}.*?(?<!del)ins.*?'),
    'inversion (other)': re.compile(f'{all_other_sequence}.*?inv.*?'),
    'delins (other)': re.compile(f'{all_other_sequence}.*?delins.*?'),
}
print_example_matches(dataset, variant_regex, size=5)
variant_counts = count_hgvs(dataset, variant_regex, include_no_hgvs=False, exclusive=False)
plt.figure(figsize=(15,7))
plt.title('Variant type')
plt.xticks(rotation='vertical')
plt.bar(variant_counts.keys(), variant_counts.values())
variant_counts
Explanation: Variant types
Documentation here
Compare stats here - distribution over all of ClinVar, using XML attribute not HGVS
Top of page
End of explanation
genomic_sequence = f'^{sequence_identifier}g\.'
coding_sequence = f'^{sequence_identifier}c\.'
noncoding_sequence = f'^{sequence_identifier}n\.'
other_sequence = f'^{sequence_identifier}[abdefh-mo-z]\.' # r'^' + sequence_identifier + r'[a-fh-z]\.'
num_range = r'[0-9]+_[0-9]+'
unk_range = r'[0-9?]+_[0-9?]+'
ch = r'[^?_+-]'  # we allow characters on either side of the range, but none of this guff
# g.12345_12678del
def definite_range(sequence_type):
    return re.compile(f'{sequence_type}{ch}*?{num_range}{ch}*?$')
# g.(123456_234567)_(345678_456789)del
def variable_range(sequence_type):
    return re.compile(f'{sequence_type}{ch}*?\({num_range}\)_\({num_range}\){ch}*?$')
# g.(?_234567)_(345678_?)del
def unknown_range(sequence_type):
    return re.compile(f'{sequence_type}{ch}*?(?=.*?\?.*?)\({unk_range}\)_\({unk_range}\){ch}*?$')
range_regex = {
    'definite (genomic)': definite_range(genomic_sequence),
    'variable (genomic)': variable_range(genomic_sequence),
    'unknown (genomic)': unknown_range(genomic_sequence),
    
    'definite (coding)': definite_range(coding_sequence),
    'variable (coding)': variable_range(coding_sequence),
    'unknown (coding)': unknown_range(coding_sequence),
    
    'definite (noncoding)': definite_range(noncoding_sequence),
    'variable (noncoding)': variable_range(noncoding_sequence),
    'unknown (noncoding)': unknown_range(noncoding_sequence),
    
    'definite (other)': definite_range(other_sequence),
    'variable (other)': variable_range(other_sequence),
    'unknown (other)': unknown_range(other_sequence),
}
print_example_matches(dataset, range_regex, size=5)
range_counts = count_hgvs(dataset, range_regex, include_no_hgvs=False)
plt.figure(figsize=(15,7))
plt.xticks(rotation='vertical')
plt.title('Ranges')
plt.bar(range_counts.keys(), range_counts.values())
range_counts
# genomic ranges
1735 + 559 + 9311
# coding / noncoding ranges
264 + 79 + 58
Explanation: Ranges
Mostly documented here though note the first case is not uncertain.
g.12345_12678del- no uncertainty
g.(123456_234567)_(345678_456789)del - variable endpoints within a known range
g.(?_234567)_(345678_?)del - variable endpoints within an unknown range
Also split into genomic vs. other sequences - the former we can definitely deal with, not sure about the others.
See below for intronic coordinate ranges...
Top of page
End of explanation
def span_lengths(dataset, regex, limit=None):
    
    Returns all span lengths for a given regex.
    This will take the first two captured groups of the regex, convert to integers, and subtract the two.
    It will NOT be smart.
    
    n = 0
    all_spans = []
    
    for record in dataset:
        if not record.measure or not record.measure.hgvs:
            continue
        hs = [h for h in record.measure.hgvs if h is not None]
        n += 1
        if limit and n > limit:
            break
            
        for h in hs:
            m = regex.match(h)
            if m and m.group(1) and m.group(2):
                span = int(m.group(2)) - int(m.group(1)) + 1
                if span < 0:
                    print('negative span!!!', h)
                else:
                    all_spans.append(span)
                # presumably all hgvs expressions for one record have the same span, don't double count
                break
                
    return all_spans
# same as previous but with capturing groups added
def_range = r'([0-9]+)_([0-9]+)'
var_range = r'\([0-9?]+_([0-9]+)\)_\(([0-9]+)_[0-9?]+\)'
def_span_regex = re.compile(f'{genomic_sequence}{ch}*?{def_range}{ch}*?$')
var_span_regex = re.compile(f'{genomic_sequence}{ch}*?{var_range}{ch}*?$')
spans = span_lengths(dataset, def_span_regex) + span_lengths(dataset, var_span_regex)
# This is everything with a known minimum span - genomic reference sequence, X_Y or (?_X)_(Y_?)
print(len(spans))
print('Mean:', np.mean(spans))
print('Median:', np.median(spans))
print('Min:', np.min(spans))
print('Max:', np.max(spans))
# actually reasonable spans...
MAX_REASONABLE_SPAN = 20000 #100000
smaller_spans = [x for x in spans if x < MAX_REASONABLE_SPAN]
print(len(smaller_spans))
plt.figure(figsize=(15,10))
plt.grid(visible=True)
plt.title(f'Minimum Spans (less than {MAX_REASONABLE_SPAN})')
# first array is counts per bin
# second array is left edges of bins, plus last right edge
plt.hist(smaller_spans, bins=100)
# VEP acceptable spans
vep_spans = [x for x in spans if x < 5000]
print(len(vep_spans))
Explanation: Span lengths
Lengths of spans - definite or minimum provided.  For now only compute for genomic sequences, as getting span length for coding/noncoding is a bit more complicated.
Top of page
End of explanation
coding_sequence = r'^' + sequence_identifier + r'c\.'
other_sequence = r'^' + sequence_identifier + r'[abd-z]\.'
pivot = r'[-*]?[0-9]+'
offset = r'[+-][0-9]+'
endpoint = pivot + offset
num_range = f'{endpoint}_{endpoint}'
unk_range = f'(?:{endpoint}|\?)_(?:{endpoint}|\?)'
ch = r'[^?_+-]'  # we allow characters on either side of the range, but none of this guff
irange_regex = {
    'definite intron (coding)': re.compile(coding_sequence + f'{ch}*?{num_range}{ch}*?$'),
    'variable intron (coding)': re.compile(coding_sequence + f'{ch}*?\({num_range}\)_\({num_range}\){ch}*?$'),
    'unknown intron (coding)': re.compile(coding_sequence + f'{ch}*?(?=.*?\?.*?)\({unk_range}\)_\({unk_range}\){ch}*?$'),
    'definite intron (other)': re.compile(other_sequence + f'{ch}*?{num_range}{ch}*?$'),
    'variable intron (other)': re.compile(other_sequence + f'{ch}*?\({num_range}\)_\({num_range}\){ch}*?$'),
    'unknown intron (other)': re.compile(other_sequence + f'{ch}*?(?=.*?\?.*?)\({unk_range}\)_\({unk_range}\){ch}*?$'),
}
print_example_matches(dataset, irange_regex, size=5, include_none=False)
irange_counts = count_hgvs(dataset, irange_regex, include_no_hgvs=False, include_none=False)
plt.figure(figsize=(15,7))
plt.title('Ranges')
plt.bar(irange_counts.keys(), irange_counts.values())
irange_counts
sum(irange_counts.values())
Explanation: Intronic numbering
Documentation here and particularly the diagram.
Only used for coding and non-coding reference sequences, so not relevant if we focus on genomic.
Top of page
End of explanation
<END_TASK>