<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Follow-up questions - 02/02/2022
Table of contents
Phenotypes
Summary
Uncertain ranges
Spans
Genes
Summary
Step1:  Phenotypes
What are the phenotypes of the structural variant? Are they somatic or germline?
Top of page
Step2:  Summary for phenotypes
somatic slightly enriched in complex variants compared to germline, but not much
complex are 1.77% somatic (vs. 0.878% overall)
Top of page
Uncertain ranges
It would be interesting to know the size of the uncertainty when it is known and compare it to the known inner range (ratio between certain range and uncertain range).
Top of page
Imprecise - known bounds
Ex. NC_000011.8
Step3:  Precise
Ex. NC_000016.10
Step7:  Uncertainty from spans
For this we only include the imprecise, known bounds case.  Note numbers from previous notebook
precise
Step9:  Uncertainty from genes
Top of page
<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 consequence_prediction.vep_mapping_pipeline.consequence_mapping import *
from eva_cttv_pipeline.clinvar_xml_io.clinvar_xml_io import *
from eva_cttv_pipeline.evidence_string_generation.clinvar_to_evidence_strings import convert_allele_origins
%matplotlib inline
import matplotlib.pyplot as plt
Explanation: Follow-up questions - 02/02/2022
Table of contents
Phenotypes
Summary
Uncertain ranges
Spans
Genes
Summary
End of explanation
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 get_somatic_germline_counts(dataset):
    all_allele_origins = [convert_allele_origins(record.valid_allele_origins) for record in dataset]
    
    # Our pipeline's definition for distinguishing somatic & germline
    def is_somatic(allele_origins):
        return allele_origins == ['somatic']
    phenotypes_counts = Counter()
    for allele_origins in all_allele_origins:
        germline = False
        somatic = False
        for ao in allele_origins:
            if is_somatic(ao):
                somatic = True
            else:
                germline = True
        if germline and somatic:
            phenotypes_counts['both'] += 1
        if germline and not somatic:
            phenotypes_counts['germline'] += 1
        if somatic and not germline:
            phenotypes_counts['somatic'] += 1
            
    # flat count of allele origins
    flattened_allele_origins = [x for allele_origins in all_allele_origins for ao in allele_origins for x in ao]
    flat_pheno_counts = Counter(flattened_allele_origins)
            
    return phenotypes_counts, flat_pheno_counts
complex_phenotypes, complex_flat_aos = get_somatic_germline_counts(dataset)
complex_phenotypes
# check if these are enriched in somatic relative to full set
full_dataset = ClinVarDataset(os.path.join(PROJECT_ROOT, 'ClinVarFullRelease_2021-07.xml.gz'))
full_phenotypes, full_flat_aos = get_somatic_germline_counts(full_dataset)
full_phenotypes
def percent_somatic(c):
    return c['somatic'] / sum(c.values()) * 100.0
print('percent somatic for complex:', percent_somatic(complex_phenotypes))
print('percent somatic for all:', percent_somatic(full_phenotypes))
Explanation: Phenotypes
What are the phenotypes of the structural variant? Are they somatic or germline?
Top of page
End of explanation
(67974756 - 67967551) / (67974774 - 67967534)
(67967551 - 67967534) + (67974774 - 67974756)
Explanation: Summary for phenotypes
somatic slightly enriched in complex variants compared to germline, but not much
complex are 1.77% somatic (vs. 0.878% overall)
Top of page
Uncertain ranges
It would be interesting to know the size of the uncertainty when it is known and compare it to the known inner range (ratio between certain range and uncertain range).
Top of page
Imprecise - known bounds
Ex. NC_000011.8:g.(67967534_67967551)_(67974756_67974774)del
Certainty ratio = smallest possible span divided by largest possible span
Also interested in absolute size of uncertain bounds region = range of possible start points + range of possible end points
End of explanation
(26775295 - 26547773) / 101991189  # length of chr 15
Explanation: Precise
Ex. NC_000016.10:g.12595039_12636793del
Certainty of 1 (no uncertainty)
Imprecise - unknown bounds
Ex. NC_000015.10:g.(?_26547773)_(26775295_?)del
In theory we could compute the same numbers using the full length of the reference sequence, but not sure it's meaningful.
End of explanation
sequence_identifier = r'[a-zA-Z0-9_.]+'
genomic_sequence = f'^({sequence_identifier}):g\.'
# only INS, DEL, DUP supported by VEP
variant_type_regex = {
    re.compile(f'{genomic_sequence}.*?del(?!ins).*?') : 'DEL',
    re.compile(f'{genomic_sequence}.*?dup.*?') : 'DUP',
    re.compile(f'{genomic_sequence}.*?(?<!del)ins.*?') : 'INS',
}
# for this we EXCLUDE unknown bounds, and capture all numeric bounds on endpoints
def_range = r'([0-9]+)_([0-9]+)'
var_range = r'\(([0-9]+)_([0-9]+)\)_\(([0-9]+)_([0-9]+)\)'
ch = r'[^?_+-]'
def_span_regex = re.compile(f'{genomic_sequence}{ch}*?{def_range}{ch}*?$')
var_span_regex = re.compile(f'{genomic_sequence}{ch}*?{var_range}{ch}*?$')
def endpoint_bounds(dataset, include_precise=False, limit=None):
    Returns inner and outer bounds on endpoints (duplicating inner/outer if precise).
    n = 0
    all_bounds = []
    all_hgvs = []
    
    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:
            # NC_000011.8:g.(67967534_67967551)_(67974756_67974774)del
            var_match = var_span_regex.match(h)
            if var_match and all(var_match.group(i) for i in range(2,6)):
                # use terminology from dbVar data model
                # see https://www.ncbi.nlm.nih.gov/core/assets/dbvar/files/dbVar_VCF_Submission.pdf
                outer_start = int(var_match.group(2))
                inner_start = int(var_match.group(3))
                inner_stop = int(var_match.group(4))
                outer_stop = int(var_match.group(5))   
                all_bounds.append(((outer_start, inner_start), (inner_stop, outer_stop)))
                all_hgvs.append(h)
                # presumably all hgvs expressions for one record have the same span, don't double count
                break
                
            elif include_precise:
                # NC_000016.10:g.12595039_12636793del
                def_match = def_span_regex.match(h)
                if def_match and def_match.group(2) and def_match.group(3):
                    outer_start = inner_start = int(def_match.group(2))
                    inner_stop = outer_stop = int(def_match.group(3))
                    all_bounds.append(((outer_start, inner_start), (inner_stop, outer_stop)))
                    all_hgvs.append(h)
                    break
                
    return all_hgvs, all_bounds
all_hgvs, all_bounds = endpoint_bounds(dataset)
def is_valid(bounds):
    # invalid if any range is negative
    return (bounds[0][1] >= bounds[0][0]
           and bounds[1][1] >= bounds[1][0]
           and bounds[1][0] >= bounds[0][1])
def certainty_ratio(bounds):
    For an HGVS range (A_B)_(C_D), this computes (C-B) / (D-A)
    return (bounds[1][0] - bounds[0][1]) / (bounds[1][1] - bounds[0][0])
def uncertain_bounds_region(bounds):
    For an HGVS range (A_B)_(C_D), this computes (A-B) + (D-C)
    return (bounds[0][1] - bounds[0][0]) + (bounds[1][1] - bounds[1][0])
len(all_bounds)
all_valid_bounds = [bounds for bounds in all_bounds if is_valid(bounds)]
len(all_valid_bounds)
all_certainty_ratios = [certainty_ratio(bounds) for bounds in all_valid_bounds]
all_uncertain_ranges = [uncertain_bounds_region(bounds) for bounds in all_valid_bounds]
# 1.0 is the most certain
print(min(all_certainty_ratios))
print(max(all_certainty_ratios))
plt.figure(figsize=(15,10))
plt.grid(visible=True)
plt.title(f'Variants per certainty ratio (imprecise, known bounds)')
plt.hist(all_certainty_ratios, bins=100)
print(min(all_uncertain_ranges))
print(max(all_uncertain_ranges))
# exclude the max
i = all_uncertain_ranges.index(max(all_uncertain_ranges))
plt.figure(figsize=(15,10))
plt.grid(visible=True)
plt.title('Variants per total size of uncertain bounds region')
plt.hist(all_uncertain_ranges[:i] + all_uncertain_ranges[i+1:], bins=100)
# the max is screwing up all my plots, get rid of it
i = all_uncertain_ranges.index(max(all_uncertain_ranges))
xs = all_certainty_ratios[:i] + all_certainty_ratios[i+1:]
ys = all_uncertain_ranges[:i] + all_uncertain_ranges[i+1:]
print(all_uncertain_ranges[i])
print(all_certainty_ratios[i])
plt.figure(figsize=(12,10))
plt.grid(visible=True)
plt.title('Certainty ratio vs. size of uncertain bounds region')
plt.xlabel('Certainty ratio')
plt.ylabel('Size of uncertain bounds region')
plt.scatter(xs, ys, marker='.')
plt.figure(figsize=(12,10))
plt.grid(visible=True)
plt.title('Certainty ratio vs. size of uncertain bounds region (log scale)')
plt.xlabel('Certainty ratio')
plt.ylabel('Size of uncertain bounds region')
plt.yscale('log')
plt.scatter(xs, ys, marker='.')
Explanation: Uncertainty from spans
For this we only include the imprecise, known bounds case.  Note numbers from previous notebook
precise: 1735
uncertainty is 0
imprecise (known bounds): 559
measures make sense, numbers below
imprecise (unknown bounds, mostly unknown outer bounds): 9311
measures might not make sense
Top of page
End of explanation
def hgvs_and_bounds_to_vep_identifier(all_hgvs, all_bounds):
    for hgvs, bounds in zip(all_hgvs, all_bounds):
        m = def_span_regex.match(hgvs)
        if not m:
            m = var_span_regex.match(hgvs)
            if not m:
                continue
        seq = m.group(1)
        # not everything accepted by VEP, for now we'll be lazy
        if not (seq.startswith('NC') or seq.startswith('LRG') or seq.startswith('NW') or seq.startswith('AC')):
            continue
        variant_type = None
        for r, s in variant_type_regex.items():
            if r.match(hgvs):
                variant_type = s
                break
        if not variant_type:
            continue
        # yield both inner and outer bounds
        # include inner/outer in the identifier so we can connect them later
        yield f'{seq} {bounds[0][1]} {bounds[1][0]} {variant_type} + {hgvs}###INNER'
        yield f'{seq} {bounds[0][0]} {bounds[1][1]} {variant_type} + {hgvs}###OUTER'
# modified from previous notebook...
def grouper(iterable, n):
    args = [iter(iterable)] * n
    return [x for x in zip_longest(*args, fillvalue=None) if x is not None]
def get_vep_results(all_hgvs, all_bounds):   
    variants = [v for v in hgvs_and_bounds_to_vep_identifier(all_hgvs, all_bounds) if v]
    print(f'{len(variants)} parsed into chrom/start/end/type')
    # VEP only accepts batches of 200
    vep_results = []
    for group in grouper(variants, n=200):
        vep_results.extend(query_vep(variants=group, search_distance=VEP_SHORT_QUERY_DISTANCE))
        
    return vep_results
def extract_genes(vep_results):
    results_by_variant = defaultdict(list)
    for result in vep_results:
        variant_identifier = result['id']
        consequences = result.get('transcript_consequences', [])
        results_by_variant[variant_identifier].extend({c['gene_id'] for c in consequences})
    return results_by_variant
def gene_counts(all_hgvs, all_bounds, limit=None):
    Return a map: hgvs -> (num affected genes inner, num affected genes outer)
    if limit:
        vep_results = get_vep_results(all_hgvs[:limit], all_bounds[:limit])
    else:
        vep_results = get_vep_results(all_hgvs, all_bounds)
        
    identifiers_to_genes = extract_genes(vep_results)
    print(f'{len(identifiers_to_genes)} successfully mapped by VEP')
    
    result = defaultdict(dict)
    for identifier, genes in identifiers_to_genes.items():
        hgvs, inner_or_outer = identifier.split('###')
        result[hgvs][inner_or_outer] = len(genes)
    
    return result
result = gene_counts(all_hgvs, all_bounds)
def certainty_ratio_genes(genes_dict):
    return genes_dict['INNER'] / genes_dict['OUTER']
# don't think the UBR size measurement makes sense
all_genes_ratios = [certainty_ratio_genes(x) for x in result.values()]
print(len(all_genes_ratios))
print(min(all_genes_ratios))
print(max(all_genes_ratios))
plt.figure(figsize=(15,10))
plt.grid(visible=True)
plt.title(f'Variants per certainty ratio (target genes version)')
plt.hist(all_genes_ratios, bins=100)
Explanation: Uncertainty from genes
Top of page
End of explanation
<END_TASK>