<jupyter_start><jupyter_text>Data Structuring and Pruning<jupyter_code># Load datasets

import json
import pathlib
import importlib
from collections import defaultdict, Counter
import pyupset as pyu
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import re
from matplotlib.gridspec import GridSpec
import csv
import requests
import pickle
import os


from vicckb import model as viccdb
from vicckb.definitions import DATA_ROOT, PROJECT_ROOT
%matplotlib inline

REPOPATH = PROJECT_ROOT.parent
OUTPATH = REPOPATH / 'out'
FIGPATH = OUTPATH / 'figures'

os.makedirs(FIGPATH, exist_ok=True)
# reload module and load data
importlib.reload(viccdb)
vdb = viccdb.ViccDb()
vdb.report_groups()
# for now, omit brca from analysis
brca = vdb.select(lambda x: x['source'] == 'brca')
core_vdb = vdb - brca
core_vdb.report_groups()
# remove biological associations
oncokb_biological = core_vdb.select(lambda x: x['source'] == 'oncokb' and 'biological' in x['raw'])
oncokb_biological.report_groups(core_vdb)
core_vdb = core_vdb - oncokb_biological<jupyter_output>oncokb: 3801 (93.9% of superset)
Total: 3801 (22.5% of superset)<jupyter_text>Evidence UniquenessThis section deals with non-unique entries from the database. This is a temporary measure until the importers are fixed. As such, it uses private variables and non-standard methods to hack around the built-in uniqueness assumptions that are violated by these data. Remove this entire section once the source hash checks pass.<jupyter_code>core_vdb.report_groups()
# Non-unique raw entries
raw_duplicates = core_vdb.select(lambda x: len(core_vdb._hashed[hash(x)]) > 1)
raw_duplicates.report_groups(core_vdb)
cgi_dups = raw_duplicates.by_source('cgi')
cgi_clean = [x for x in cgi_dups if x['raw']['Drug status']]
test = viccdb.ViccDb([x for x in core_vdb if x not in cgi_dups] + list(cgi_clean))
test.report_groups(core_vdb)
# Test matches expectation, moving to core
core_vdb = test
pmkb_dups = raw_duplicates.by_source('pmkb')
len(pmkb_dups._hashed)
len(pmkb_dups)
merged_associations = list()
for hash_key, equivalent_associations in pmkb_dups._hashed.items():
    root_association = equivalent_associations.pop()
    for other_association in equivalent_associations:
        root_association['features'].append(other_association['features'][0])
    merged_associations.append(root_association)
test = viccdb.ViccDb([x for x in core_vdb if x not in merged_associations] + list(merged_associations))
test.report_groups(core_vdb)
x = len(core_vdb.by_source('pmkb')) - len(pmkb_dups) + len(pmkb_dups._hashed)
print("Expecting {} associations for PMKB".format(x))
# Test matches expectation, moving to core
core_vdb = test<jupyter_output><empty_output><jupyter_text>Evidence filteringRemoval of all evidence without associated publications, followed by removal of all associations without evidence.<jupyter_code>def clean_refs(association):
    evidences = association['association']['evidence']
    evidence_indices_to_delete = list()
    for i, evidence in enumerate(evidences):
        assert isinstance(publications, list)
        publications = [x for x in evidence['info']['publications'] if x]
        evidence['info']['publications'] = publications
        if not publications:
            evidence_indices_to_delete.append(i)
    for index in sorted(evidence_indices_to_delete, reverse=True):
        del association['association']['evidence'][index]

map(clean_refs, core_vdb)
core_missing_ref = core_vdb.select(lambda x: not any(x.publications))
core_missing_ref.report_groups(core_vdb)
core_vdb = core_vdb - core_missing_ref
core_vdb.report_groups()
# All associations should have an evidence level

core_vdb.select(lambda x: not x.evidence_level).report_groups(core_vdb)<jupyter_output>Total: 0 (0.0% of superset)<jupyter_text>Feature coordinatesWhat follows is a detailed look at associations without start and end coordinates after normalization, and a set of regular expression filters to separate out these associations into chunks that can be annotated with gene- or exon-level coordinates, as appropriate.<jupyter_code># Associations with more than 1 feature
x = core_vdb.select(lambda x: len(x.features) > 1)
x.report_groups(vdb)
# Associations without at least 1 complete and valid feature
no_features = core_vdb.select(lambda x: len(x.features) == 0)
no_features.report_groups(vdb)
vdb[0]['association']['phenotype']
# Associations with coordinate features
coord_featured = core_vdb - no_features
coord_featured.report_groups(core_vdb)<jupyter_output>cgi: 1063 (99.2% of superset)
civic: 3323 (99.5% of superset)
jax: 5736 (99.8% of superset)
molecularmatch: 2063 (99.2% of superset)
oncokb: 245 (99.2% of superset)
pmkb: 369 (99.5% of superset)
Total: 12799 (99.6% of superset)<jupyter_text>Fix PMKB features<jupyter_code>with open('vicckb/data/gene_strand.pkl', 'rb') as f:
    gene_strand = pickle.load(f)
COMPLEMENT = {
    'A': 'T',
    'T': 'A',
    'G': 'C',
    'C': 'G',
    '-': '-'
}

complement_map = str.maketrans(COMPLEMENT)

def get_gene_strand(gene, trx_id):
    strand = gene_strand.get((gene, trx_id), None)
    if strand is None:
        cmd = '''wget -q -O - 'http://grch37.ensembl.org/biomart/martservice?query=<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE Query>
<Query  virtualSchemaName = "default" formatter = "CSV" header = "0" uniqueRows = "0" count = "" datasetConfigVersion = "0.6" >
<Dataset name = "hsapiens_gene_ensembl" interface = "default" >
<Filter name = "hgnc_symbol" value = "{}"/>
<Attribute name = "ensembl_transcript_id" />
<Attribute name = "strand" />
</Dataset>
</Query>' | grep '{}'
'''.format(gene, trx_id)
        result = !{cmd}
        strand = result[0].split(',')[-1]
        gene_strand[(gene, trx_id)] = strand
        return strand
    return strand

for association in core_vdb.by_source('pmkb'):
    del(association._features)
    for feature in association['features']:
        gene = feature['geneSymbol']
        trx_id = feature['attributes']['transcript']['string_value']
        strand = get_gene_strand(gene, trx_id)
        if strand == '-1':
            # feature['ref'] = feature['ref'][::-1].translate(complement_map)
            try:
                feature['alt'] = feature['alt'][::-1].translate(complement_map)
            except KeyError:
                continue

with open('vicckb/data/gene_strand.pkl', 'wb') as f:
    pickle.dump(gene_strand, f)<jupyter_output><empty_output><jupyter_text>Remainder of section is inactivated code for identifying associations without coordinatesimport redef feature_filter(re_obj, associations):     report matches and return non-matches    found = list(filter(lambda x: re_obj.search(x['feature_names']) is not None, associations))    not_found = list(filter(lambda x: re_obj.search(x['feature_names']) is None, associations))    report_groups(found)    return(not_found)amp_re = re.compile(r'(amplification)|(loss)|(amp)', re.IGNORECASE)remainder = feature_filter(amp_re, no_partial_coord_featured_with_feature_names) fusion_re = re.compile(r'(\w{2,}-\w{2,})|(fusion)', re.IGNORECASE)r2 = feature_filter(fusion_re, remainder) ppm_re = re.compile(r'\w+(:| )[a-z]\d+[a-z]?(fs\*?)?$', re.IGNORECASE)r3 = feature_filter(ppm_re, r2) indel_re = re.compile(r'\w+(:| )\w+(ins\w+)|(del($|ins\w+))|(dup$)')r4 = feature_filter(indel_re, r3) bucket_re = re.compile(r'[A-Z0-9]+( (in)?act)?( oncogenic)? mut((ant)|(ation))?$')r5 = feature_filter(bucket_re, r4) exon_re = re.compile(r'exon', re.IGNORECASE)r6 = feature_filter(exon_re, r5) expression_re = re.compile(r'(exp)|(^\w+ (pos(itive)?)|(neg(ative)?)|(biallelic inactivation)$)|(truncating)|(deletion)', re.IGNORECASE)r7 = feature_filter(expression_re, r6) report_groups(r7) get_feature_names([x for x in r7 if x['source'] == 'cgi'])  Diseases<jupyter_code>disease_missing = core_vdb.select(lambda x: x.disease is None)
disease_missing.report_groups(core_vdb)
disease_missing[0]['association']
# Fix DOID

for association in core_vdb.select(lambda x: x.disease and x.disease.name == 'CNS Cancer'):
    association['association']['phenotype']['type']['term'] = 'central nervous system cancer'
    association['association']['phenotype']['type']['id'] = 'DOID:3620'
mismatched_do_cgi = core_vdb.select(lambda x: x.disease and \
                                x.disease.source.lower().endswith('doid') and \
                                not x.disease.id.startswith('DOID') and \
                                x.source == 'cgi'
                               )
mismatched_do_cgi.report_groups(core_vdb)
mismatched_do_cgi[0]['association']['phenotype']
mismatched_do_cgi[0]['raw']
from vicckb.definitions import DATA_ROOT
from vicckb.harmonizers import DiseaseHarmonizer
adh = DiseaseHarmonizer(map_file=(DATA_ROOT / 'disease_alias.tsv'),
                             disease_ontology='DOID')
disease = mismatched_do_cgi[0]['raw']['Primary Tumor type']
adh.harmonize(disease)
for cgi_association in mismatched_do_cgi:
    disease = cgi_association['raw']['Primary Tumor type']
    harmonized = adh.harmonize(disease)
    cgi_association['association']['phenotype'] = {
        'description': disease,
        'type': {
            'id': harmonized['id'],
            'source': harmonized['ontology'],
            'term': harmonized['term']
        }
    }
mismatched_do_cgi = core_vdb.select(lambda x: x.disease and \
                                x.disease.source.lower().endswith('doid') and \
                                not x.disease.id.startswith('DOID') and \
                                x.source == 'cgi'
                               )
mismatched_do_cgi.report_groups(core_vdb)
core_vdb.select(lambda x: not x.disease and x.source == 'cgi').report_groups()
other_do_cgi = core_vdb.select(lambda x: x.disease and \
                x.source == 'cgi' and \
                not x.disease.source.lower().endswith('doid'))
other_do_cgi.report_groups()
for cgi_association in other_do_cgi:
    disease = cgi_association['raw']['Primary Tumor type']
    harmonized = adh.harmonize(disease)
    cgi_association['association']['phenotype'] = {
        'description': disease,
        'type': {
            'id': harmonized['id'],
            'source': harmonized['ontology'],
            'term': harmonized['term']
        }
    }
core_vdb.select(lambda x: x.disease and \
                x.source == 'cgi' and \
                not x.disease.source.lower().endswith('doid')).report_groups()<jupyter_output>0 total associations<jupyter_text>Drugs<jupyter_code>drugs_missing = core_vdb.select(lambda x: len(x.drugs) == 0)
drugs_missing.report_groups(core_vdb)<jupyter_output>cgi: 112 (10.4% of superset)
civic: 1261 (37.8% of superset)
jax: 457 (8.0% of superset)
molecularmatch: 120 (5.8% of superset)
oncokb: 8 (3.2% of superset)
pmkb: 371 (100.0% of superset)
Total: 2329 (18.1% of superset)<jupyter_text>Genes<jupyter_code>ambiguous = list()
for a in core_vdb:
    a.genes

# assert len(ambiguous) == 0  # Ensure there are no ambiguous genes from knowledgebases<jupyter_output>/Users/awagner/Workspace/git/vicckb/vicckb/model.py:257: UserWarning: Ambiguous gene symbol MLL2 in assertion 235677252030682
  warn('Ambiguous gene symbol {} in assertion {}'.format(g, self))<jupyter_text>CacheSaving core_vdb to cache for testing.<jupyter_code>core_vdb.cache_data()<jupyter_output><empty_output><jupyter_text>Knowledgebase Comparison  Publications All publications<jupyter_code>x = core_vdb.plot_element_by_source('publications', min_bound=4)
f = x['figure']
# f.savefig('out/publications.pdf')
# x['input_data']
## For Sidi

ebs = x['input_data']
df_dict = dict()
for source in ebs:
    fe = list(filter(lambda x: bool(x), ebs[source]))
    df_dict[source] = pd.DataFrame(fe, columns=['attribute'])

with open('example.pkl', 'wb') as f:
    pickle.dump(df_dict, f)
with open('example.pkl', 'rb') as f:
    df_dict = pickle.load(f)
!open .
# Publications uniquely cited
data = x['input_data']
total = 0
for source in data:
    publications_from_elsewhere = set()
    for source2 in data:
        if source == source2:
            continue
        publications_from_elsewhere.update(data[source2])
    unique = data[source] - publications_from_elsewhere
    print("{}: {} resource-specific publications".format(source, len(unique)))
    total += len(unique)
print("{} ({:.2%}) total resource-specific publications".format(total, total / len(set.union(*(data.values())))))
p_sets = core_vdb.get_element_by_source('publications')

len((p_sets['civic'] & p_sets['jax']) - p_sets['pmkb'] - p_sets['oncokb'] - p_sets['molecularmatch'] - p_sets['cgi'])
x = core_vdb.element_by_source_stats('publications')
x['ubiquitous'] # Bose et al. Cancer Discovery 2013<jupyter_output>3696 / 4354 (84.89%) of publications are represented in only 1 resource.
203 / 4354 (4.66%) of publications are represented in the majority of (3) resources.
1 / 4354 (0.02%) of publications are represented across all resources.<jupyter_text>PMIDs<jupyter_code>x = core_vdb.plot_element_by_source('publications', lambda x: x.pmid, min_bound=3)
# f = x['figure']
# f.savefig('out/pmids.pdf')
x = core_vdb.element_by_source_stats('publications', lambda x: x.pmid)
x['ubiquitous'] # Bose et al. Cancer Discovery 2013<jupyter_output>3146 / 3800 (82.79%) of publications are represented in only 1 resource.
203 / 3800 (5.34%) of publications are represented in the majority of (3) resources.
1 / 3800 (0.03%) of publications are represented across all resources.<jupyter_text>Genes<jupyter_code>g_set = core_vdb.get_element_by_source('genes')
x = g_set['civic']
for n, s in g_set.items():
    if n == 'civic':
        continue
    x = x - s
len(x)
no_genes = core_vdb.select(lambda x: not x.genes)
no_genes.report_groups(core_vdb)
with_genes = core_vdb - no_genes
x = with_genes.plot_element_by_source('genes')
# f = x['figure']
# f.savefig('out/genes.pdf')
x = with_genes.element_by_source_stats('genes')
x['ubiquitous']
stats.fisher_exact([
        [203, 4151],
        [97, 318]
    ])<jupyter_output><empty_output><jupyter_text>Features<jupyter_code># suddenly stopped working? maybe a dependency error?

# x = core_vdb.plot_element_by_source('features', min_bound=5)
# f = x['figure']
# f.savefig(str(FIGPATH / 'misc_figures' / 'feature_upset.pdf'))
x = core_vdb.select(lambda x: x.evidence_level == 'A').element_by_source_stats('features')
f_sets = core_vdb.get_element_by_source('features')
# all_features = set.update()
count = Counter()
for s in f_sets.values():
    count.update(s)
# non-scientific-notation calculations

sources = list(f_sets)
for source in sources:
    s = set()
    for other in sources:
        if source == other:
            continue
        s.update(f_sets[other])
    print(f'{source}: {len(f_sets[source] - s)}')
    
cgi_and_okb = f_sets['oncokb'] & f_sets['cgi']
others = f_sets['pmkb'] | f_sets['molecularmatch'] | f_sets['civic'] | f_sets['jax']
print(f'cgi and oncokb: {len(cgi_and_okb - others)}')
count_of_counts = Counter()
count_of_counts.update(count.values())
np.asarray(list(count_of_counts.values()))
labels = sorted(count_of_counts.keys())
values = [count_of_counts[x] for x in labels]
values
fig1, ax1 = plt.subplots()

pie, _, _ = ax1.pie(values, radius=1,
        labels=labels, autopct='%1.1f%%', pctdistance=2)
ax1.axis('equal')
plt.setp(pie, edgecolor='white')

# plt.savefig(str(FIGPATH / 'misc_figures' / 'feature_overlap.pdf'))
plt.show()
ubiquitous_features = list(x['ubiquitous'])
sorted([x.name for x in ubiquitous_features])
def element_uniqueness_across_kbs(element, as_proportion=False):
    e_sets = core_vdb.get_element_by_source(element)
    count = Counter()
    for s in e_sets.values():
        count.update(s)
    count_of_counts = Counter()
    count_of_counts.update(count.values())
    labels = sorted(count_of_counts.keys())
    if as_proportion:
        d = sum(count_of_counts.values())
        values = [count_of_counts[x]/d for x in labels]
    else:
        values = [count_of_counts[x] for x in labels]
    return dict(zip(labels, values))
element_uniqueness_across_kbs('features', as_proportion=True)
g_prop = element_uniqueness_across_kbs('genes', as_proportion=True)
f_prop = element_uniqueness_across_kbs('features', as_proportion=True)
di_prop = element_uniqueness_across_kbs('disease', as_proportion=True)
dr_prop = element_uniqueness_across_kbs('drugs', as_proportion=True)
pub_prop = element_uniqueness_across_kbs('publications', as_proportion=True)
labels = ['Genes', 'Features', 'Diseases', 'Drugs*', 'Publications']
N = len(labels)
ind = np.arange(N)
value_sets = [g_prop, f_prop, di_prop, dr_prop, pub_prop]
width = 0.5
plot_sets = [[x.get(i, 0) for x in value_sets] for i in range(1,7)]
plots = list()
b_sums = np.zeros(N)
for plot_set in plot_sets:
    p = plt.bar(ind, plot_set, width, bottom=b_sums)
    b_sums += np.array(plot_set)
    plots.append(p)

plt.ylabel('Proportion')
plt.xticks(ind, labels)
plt.legend([p[0] for p in plots], range(1,7), bbox_to_anchor=(1.2, 0.5, 0, 0.5))

# plt.savefig(str(FIGPATH / 'misc_figures' / 'elements_overlap.pdf'))
plt.show()<jupyter_output><empty_output><jupyter_text>Tier 1 variants<jupyter_code>tier1 = core_vdb.select(lambda x: x.evidence_level in ['A', 'B'])
genes = set()
features = set()
for a in tier1:
    genes.update(a.genes)
    features.update(a.features)
print(len(genes))
print(len(features))<jupyter_output>236
1512<jupyter_text>Hierarchical searchExisting method is to find an exact match of any features for an association. Below we demonstrate gains through hierarchical search of GenomicFeatures, a core result of this effort.<jupyter_code>hits = core_vdb.search_by_feature(chromosome=7, start=140453136, end=140453136, reference_name='GRCh37')

v600k = [x['best_match']['p'] for x in hits if x['best_match']['feature'].name.endswith('V600K')]
len(v600k)
unique_features = set()
x = [x.features for x in core_vdb]
for fset in x:
    unique_features.update(fset)
unique_features = list(unique_features)
len(unique_features)
# This is a computationally expensive operation (~1 minute for the 2800 searches). Could be sped up through indexed searching.

# feature_hits = dict()
# for feature in unique_features:
#     feature_hits[feature] = core_vdb.search_by_feature(genomic_feature=feature)

# New method

feature_hits = core_vdb.search_by_features(genomic_features=unique_features)

ranking = viccdb.ViccDb.MATCH_RANKING
hits_by_type = Counter()
sources_by_type = defaultdict(Counter)
for match_type in ranking:
    typed_hits = [x for x in feature_hits if ranking.index(x['best_match']['type']) <= ranking.index(match_type)]
    hits_by_type[match_type] = len(typed_hits)
    typed_associations_by_query = defaultdict(set)
    for hit in typed_hits:
        typed_associations_by_query[hit['query']].add(hit['association'])
    for associations in typed_associations_by_query.values():
        sources = {association.source for association in associations}
        sources_by_type[len(sources)][match_type] += 1
hits_by_type
sources_by_type
fig, ax = plt.subplots()

source_counts = sorted(sources_by_type)
groups = viccdb.ViccDb.MATCH_RANKING
width = 0.15
plot_elements = list()
ind = np.arange(len(groups))
for i, source_count in enumerate(source_counts):
    type_counts = sources_by_type[source_count]
    x = [type_counts[k] for k in groups]
    p = ax.bar(ind + width*(i-2.5), x, width, label=source_count)
    plot_elements.append(p)
ax.set_xticks(ind)
ax.set_xticklabels(groups)
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels, title='Sources')
plt.ylabel('Queries')
plt.show()
# fig.savefig(str(FIGPATH / 'misc_figures' / 'search_performance.pdf'), format='pdf')
# feature_hits = core_vdb.search_by_features(genomic_features=ubiquitous_features)
q = {hit['query'] for hit in feature_hits}
len(q)
c = Counter()
for f in unique_features:
    if len(f) > 9 and len(f) < 100:
        l = '10-99'
    elif len(f) >= 100:
        l = '100+'
    else:
        l = len(f)
    c[l] += 1
c<jupyter_output><empty_output><jupyter_text>Sequence Ontology<jupyter_code>count = Counter()
for association in core_vdb:
    for feature in association.features:
        count[feature.so.get('name', 'Uncategorized')] += 1
count.most_common(10)
(5052 + 3263) / sum(count.values())
# Associations lacking any sequence ontology ID

def no_soid(association):
    x = [feature.so.get('name', 'Uncategorized') == 'Uncategorized' for feature in association.features]
    return all(x)

no_soid_group = core_vdb.select(no_soid)
no_soid_group.report_groups(core_vdb)
# Associations lacking at least one sequence ontology ID

def missing_soid(association):
    x = [feature.so.get('name', 'Uncategorized') == 'Uncategorized' for feature in association.features]
    return any(x)

missing_soid_group = core_vdb.select(missing_soid)
missing_soid_group.report_groups(core_vdb)
partial_soid_group = missing_soid_group - no_soid_group
partial_soid_group.report_groups(core_vdb)
partial_soid_group[0].features[2].so<jupyter_output><empty_output><jupyter_text>Project GENIE feature match<jupyter_code># Loading Disease graph here for use in GENIE disease matching

import obonet
import networkx

url = 'https://raw.githubusercontent.com/DiseaseOntology/HumanDiseaseOntology/v2018-05-11/src/ontology/HumanDO.obo'
graph = obonet.read_obo(url)

assert networkx.is_directed_acyclic_graph(graph)
so_by_name = {
    "3'Flank": {'name': 'downstream_transcript_variant', 'soid': 'SO:0001987'},
    "3'UTR": {'name': '3_prime_UTR_variant', 'soid': 'SO:0001624'},
    "5'Flank": {'name': 'upstream_transcript_variant', 'soid': 'SO:0001986'},
    "5'UTR": {'name': '5_prime_UTR_variant', 'soid': 'SO:0001623'},
    "Frame_Shift_Del": {'name': 'frameshift_truncation', 'soid': 'SO:0001910'},
    "Frame_Shift_Ins": {'name': 'frameshift_elongation', 'soid': 'SO:0001909'},
    "In_Frame_Del": {'name': 'inframe_deletion', 'soid': 'SO:0001822'},
    "In_Frame_Ins": {'name': 'inframe_insertion', 'soid': 'SO:0001821'},
    "Intron": {'name': 'intron_variant', 'soid': 'SO:0001627'},
    "Missense_Mutation": {'name': 'missense_variant', 'soid': 'SO:0001583'},
    "Nonsense_Mutation": {'name': 'stop_gained', 'soid': 'SO:0001587'},
    "Nonstop_Mutation": {'name': 'stop_lost', 'soid': 'SO:0001578'},
    "Silent": {'name': 'synonymous_variant', 'soid': 'SO:0001819'},
    "Splice_Region": {'name': 'splice_region_variant', 'soid': 'SO:0001630'},
    "Splice_Site": {'name': 'splice_site_variant', 'soid': 'SO:0001629'},
    "Translation_Start_Site": {'name': 'initiator_codon_variant', 'soid': 'SO:0001582'}
}
alias_to_doids = defaultdict(list)
def map_to_doid(graph, doid):
    for _id in graph.predecessors(doid):
        map_to_doid(graph, _id)
    xrefs = graph.node[doid].get('xref', [])
    for xref in xrefs:
        source, xref_id = xref.split(':')
        alias_to_doids[(source, xref_id)].append(doid)

map_to_doid(graph, 'DOID:162') 
oncotree_to_aliases = dict()
oncotree_types_url = 'http://oncotree.mskcc.org/api/tumorTypes'
resp = requests.get(oncotree_types_url, params={'version': 'oncotree_2018_05_01'})
resp.raise_for_status()
oncotree_types = resp.json()
for o_type in oncotree_types:
    oncotree_to_aliases[o_type['code']] = []
    for source, terms in o_type['externalReferences'].items():
        for term in terms:
            oncotree_to_aliases[o_type['code']].append((source, term))
oncotree_to_doids = dict()
for o_term, aliases in oncotree_to_aliases.items():
    if not aliases:
        oncotree_to_doids[o_term] = None
        continue
    doids = set()
    for alias in aliases:
        alias_doids = alias_to_doids.get(alias, False)
        if alias_doids:
            doids.update(alias_doids)
    if doids:
        oncotree_to_doids[o_term] = list(doids)
    else:
        oncotree_to_doids[o_term] = None
patched_doids = {
#     'MAAP': 'DOID:3608',
#     'SCCNOS': 'DOID:1749',
#     'MACR': 'DOID:0050861',
#     'OCSC': 'DOID:0050866',
#     'UDMN': 'DOID:162',
#     'CUP': 'DOID:162',
#     'CUPNOS': 'DOID:162',
#     'MYF': 'DOID:4971',
#     'HGSOC': 'DOID:0050933',
#     'LGSOC': 'DOID:0050933',
#     'SOC': 'DOID:0050933',
#     'PANET': 'DOID:1798',
#     'IMT': 'DOID:0050905',
#     'OPHSC': 'DOID:0050921',
#     'MDS': 'DOID:0050908',
    'ACYC': 'DOID:0080202',
    'HL': 'DOID:8567',
    'SEM': 'DOID:4440'
}
sample_to_patient = dict()
patient_to_samples = defaultdict(list)
sample_to_doid = dict()
sample_oncotree_code = dict()
required_patches = set()
optional_patches = set()
no_patch = set()

with open(DATA_ROOT / 'GENIE_v3' / 'data_clinical_sample_3.0.0.txt') as f:
    for _ in range(4):
        f.readline()    # get past info lines
    reader = csv.DictReader(f, delimiter="\t")
    for row in reader:
        patient = row['PATIENT_ID']
        sample = row['SAMPLE_ID']
        sample_to_patient[sample] = patient
        patient_to_samples[patient].append(sample)
        oncotree_code = row['ONCOTREE_CODE'].upper()
        doids = oncotree_to_doids[oncotree_code]
        if not doids:
            required_patches.add(oncotree_code)
            doid = patched_doids.get(oncotree_code, None)
        elif len(doids) > 1:
            optional_patches.add(oncotree_code)
            doid = patched_doids[oncotree_code]
        else:
            no_patch.add(oncotree_code)
            doid = doids[0]
        sample_to_doid[sample] = doid
        sample_oncotree_code[sample] = oncotree_code
        
patient_to_samples = dict(patient_to_samples)
# Load genie variants
classification_counter = Counter()
genie_features = list()
genie_features_by_patient = defaultdict(list)
genie_features_by_variant = defaultdict(list)
genie_features_by_sample = defaultdict(list)
unfiltered_patients_with_variants = set()
EXCLUDED_CLASSIFICATIONS = [
    'Silent',
    "3'Flank",
    "3'UTR",
    "5'Flank",
    "5'UTR",
    'Intron',
    'Splice_Region'
]
with open(DATA_ROOT / 'GENIE_v3' / 'data_mutations_extended_3.0.0.txt', 'r') as maf:
    sample_list = maf.readline().strip().split(' ')[1:]
    maf_reader = csv.DictReader(maf, delimiter="\t")
    for row in maf_reader:
        start = row['Start_Position']
        end = row['End_Position']
        chromosome = row['Chromosome']
        patient = sample_to_patient[row['Tumor_Sample_Barcode']]
        unfiltered_patients_with_variants.add(patient)
        if row['Variant_Classification'] in EXCLUDED_CLASSIFICATIONS:
            continue
        if row['Reference_Allele'] != row['Tumor_Seq_Allele1']:
            alt = row['Tumor_Seq_Allele1']
        else:
            alt = row['Tumor_Seq_Allele2']
        reference = row['NCBI_Build']
        feature = viccdb.GenomicFeature(
            chromosome=chromosome,
            start=start,
            end=end,
            referenceName=reference,
            name=':'.join([row['Tumor_Sample_Barcode'], row['HGVSp_Short']]),
            geneSymbol=row['Hugo_Symbol'],
            sequence_ontology=so_by_name[row['Variant_Classification']],
            alt=alt
        )
        genie_features.append(feature)
        genie_features_by_patient[patient].append(feature)
        genie_features_by_variant[(reference, chromosome, start, end, alt)].append(feature)
        genie_features_by_sample[row['Tumor_Sample_Barcode']].append(feature)
len(genie_features)
# Do a GENIE feature search across knowledgebase. Huge search operation takes ~5.5 min to complete
from timeit import default_timer
tick = default_timer()
genie_search_results = core_vdb.search_by_features(genie_features)
tock = default_timer()
print(tock-tick)
featured_patients = set(genie_features_by_patient)
print(f'Avg. queries / second: {len(genie_features)/(tock-tick)}')
print(f'Search results: {len(genie_search_results)}')
print(f'Avg. search results / query: {len(genie_search_results)/len(genie_features)}')
genie_feature_lengths = Counter()
for feature in genie_features:
    genie_feature_lengths[len(feature)] += 1
result_size = defaultdict(Counter)
for result in genie_search_results:
    for match_type in ranking:
        if ranking.index(result['best_match']['type']) <= ranking.index(match_type):
            result_size[result['query']][match_type] += 1
print('Percentage of queries with results: {:.1%}'.format(len(result_size) / len(genie_features)))
exact_match_features = [x for x in result_size 
                        if result_size[x]['exact'] + 
                        result_size[x]['positional'] > 0]
len(exact_match_features)
print('Percentage of queries with exact results: {:.1%}'.format(len(exact_match_features) / len(genie_features)))
data = defaultdict(list)
for feature in genie_features:
    length = len(feature)
    for match_type in ranking:
        data[match_type].append([length, result_size[feature][match_type]])
colors = dict(zip(ranking, ['blue', 'green', 'orange', 'red']))
pcs = list()
labels = list()
plt.figure()
plt.xscale('log')
plt.yscale('log')
for t in reversed(ranking):
    coord_pairs = data[t]
    a = np.array(coord_pairs) + 1
    pc = plt.scatter(a[:,0], a[:,1], marker='.', color=colors[t], alpha=0.05)
    pcs.append(pc)
    labels.append(t)
plt.xlabel('Feature Size')
plt.ylabel('Interpretations')
plt.title('Feature v Result Size of GENIE Variants')
plt.legend(pcs, labels, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.tight_layout() # drops legend :-/
# plt.savefig(str(FIGPATH / 'feature_v_result_size_genie.pdf'))
# plt.savefig(str(FIGPATH / 'feature_v_result_size_genie.png'))
plt.show()
for t in reversed(ranking):
    data = [result_size[q][t] for q in genie_features
            if result_size[q][t] > 0
           ]
    plt.hist(data, alpha=1, label=t, color=colors[t], bins=150)
plt.xlabel('Interpretations Matched to Variant Query')
plt.ylabel('Queries')
plt.title('Counts of Interpretations per Variant Query by Search Strategy')
plt.legend()
# plt.savefig(str(FIGPATH / 'genie_interpretation_v_query_all.pdf'))
plt.show()
notables = Counter()
for t in reversed(ranking[:3]):
    data = [result_size[q][t] for q in genie_features
            if result_size[q][t] > 0
           ]
    for q in genie_features:
        if result_size[q][t] > 90:
            notables[(q,t)] += 1
    plt.hist(data, alpha=1, label=t, color=colors[t], bins=150)
plt.title('Counts of Interpretations per Variant Query by Search Strategy')
plt.xlabel('Interpretations Matched to Variant Query')
plt.ylabel('Queries')
plt.legend()
# plt.savefig(str(FIGPATH / 'genie_interpretation_v_query_no_regional.pdf'))
plt.show()
coord_q = dict()
for notable in notables.most_common(len(notables)):
    if notable[1] > 950:
        feature = notable[0][0]
        t = notable[0][1]
        queries = notable[1]
        HGVSp = feature.name.split(':')[1]
        coord_q[(result_size[feature][t], queries)] = (feature, t)
    else:
        break
for coord in sorted(coord_q):
    feature, t = coord_q[coord]
    HGVSp = feature.name.split(':')[1]
    queries = coord[1]
    print(f'({result_size[feature][t]}, {queries}): {feature.gene_symbol} {HGVSp} ({t})')<jupyter_output>(92, 1425): KRAS p.G12V (exact)
(94, 998): KRAS p.G12C (exact)
(111, 1062): PIK3CA p.E545K (exact)
(113, 1062): PIK3CA p.E545K (positional)
(117, 1062): PIK3CA p.E545K (focal)
(160, 1080): PIK3CA p.H1047R (exact)
(161, 1720): KRAS p.G12D (exact)
(174, 1080): PIK3CA p.H1047R (positional)
(178, 1080): PIK3CA p.H1047R (focal)
(205, 998): KRAS p.G12C (positional)
(308, 1425): KRAS p.G12V (positional)
(308, 1720): KRAS p.G12D (positional)
(343, 998): KRAS p.G12C (focal)
(457, 1425): KRAS p.G12V (focal)
(457, 1720): KRAS p.G12D (focal)
(565, 1439): BRAF p.V600E (exact)
(584, 1439): BRAF p.V600E (positional)
(655, 1439): BRAF p.V600E (focal)<jupyter_text>Disease<jupyter_code>oncotree_codes = {x['code'] for x in oncotree_types}
print(f'There are {len(oncotree_codes)} diseases in oncotree.')
no_xrefs = {x for x, v in oncotree_to_aliases.items() if not v}
print(f'Of these, {len(no_xrefs)} have no xrefs.')
print(f'{len(required_patches)} of {len(required_patches | no_patch)} oncotree diseases from GENIE do not automatically map to doid.')
print(f'Of these, {len(required_patches & no_xrefs)} have no xrefs.')
a = np.array([
        [len(required_patches & no_xrefs), len(required_patches - no_xrefs)],
        [len(oncotree_codes & no_xrefs), len(oncotree_codes - no_xrefs)]
         ])
print(a)
stats.fisher_exact(a)
c = Counter()
for x in core_vdb:
    try:
        c[x.disease.source] += 1
    except AttributeError:
        c[None] += 1
c
doid = core_vdb.select(lambda x: x.disease is not None and x.disease.id.startswith('DOID:'))
# Problem with number of entries from jax reporting "cancer" as doid type
doid.select(lambda x: x.disease.id == 'DOID:162').report_groups(core_vdb)
id_to_name = {id_: data['name'] for id_, data in graph.nodes(data=True)}
c = Counter()
for x in doid:
    try:
        c[x.disease.id] += 1
    except AttributeError:
        c[None] += 1
for k, v in c.most_common(20):
    print(f'{id_to_name[k]}: {v}')
organ_system_ids = graph.predecessors('DOID:0050686')
benign_and_premalignant_ids = ['DOID:0060072', 'DOID:0060071']
cell_type_cancer_id = 'DOID:0050687'
cancer_id = 'DOID:162'
organ_system = dict()
benign_and_premalignant = dict()

def assign_to_id(ids, id_, d):
    if len(ids) == 0:
        return
    for i in ids:
        d[i] = id_
        assign_to_id(graph.predecessors(i), id_, d)
    return 

for organ_id in (organ_system_ids + [cell_type_cancer_id]):
    assign_to_id(graph.predecessors(organ_id), organ_id, organ_system)
    organ_system[organ_id] = organ_id

for id_ in benign_and_premalignant_ids:
    assign_to_id(graph.predecessors(id_), id_, benign_and_premalignant)
    benign_and_premalignant[id_] = id_
    

return_id = organ_system['DOID:0050615']
graph.node[return_id]
normalized_disease = core_vdb.select(lambda x: x.disease is not None)
normalized_disease.report_groups(core_vdb)
do_sourced = normalized_disease.select(lambda x: x.disease.source in ['DOID', 'http://purl.obolibrary.org/obo/doid'])
do_sourced.report_groups(core_vdb)
cancer_organ_interpretations = do_sourced.select(lambda x: organ_system.get(x.disease.id, False))
benign_premalignant_interpretations = do_sourced.select(lambda x: benign_and_premalignant.get(x.disease.id, False))
cancer_interpretations = do_sourced.select(lambda x: x.disease.id == cancer_id)
other_interpretations = do_sourced - cancer_organ_interpretations - benign_premalignant_interpretations - cancer_interpretations
cancer_organ_interpretations.report_groups(do_sourced)
benign_premalignant_interpretations.report_groups(do_sourced)
cancer_interpretations.report_groups(do_sourced)
other_interpretations.report_groups(do_sourced)
x = other_interpretations.select(lambda x: x.disease.id.split(':')[0] != 'DOID')
x.report_groups(do_sourced)
c = Counter()
for i in x:
    c[i.disease.id.split(':')[0]] += 1
c
entry = x.by_source('civic')[0]
print(entry.disease.source)
print(entry.disease.id)
print(entry.disease.term)
print(entry['feature_names'])
c = Counter()
for association in cancer_organ_interpretations:
    disease_id = association.disease.id
    organ_id = organ_system[disease_id]
    c[organ_id] += 1
for k, v in c.most_common(13):
    print(f'{id_to_name[k]}: {v}')
p_cancer = cancer_organ_interpretations.select(lambda x: id_to_name[organ_system[x.disease.id]] == 'peritoneum cancer')[0]
print(p_cancer.disease)
print(p_cancer.source)
print(p_cancer.description)
len(set(organ_system.values()))
with open(DATA_ROOT / 'TopNodes_DOcancerslim_3_18.json', 'r') as f:
    result = json.load(f)
    nodes = result['graphs'][0]['nodes']
nodes[0]['id'].split('/')[-1].replace('_', ':')
topnodes_docancerslim = list()
doid_re = re.compile(r'DOID:\d+')

for node in nodes:
    doid = node['id'].split('/')[-1].replace('_', ':')
    if doid_re.match(doid):
        topnodes_docancerslim.append(doid)
len(topnodes_docancerslim)
def assign_to_nearest_id(ids, id_, d, terminals):
    if len(ids) == 0:
        return
    for i in ids:
        if i in terminals:
            assignment = i
        else:
            assignment = id_
        d[i] = assignment
        assign_to_nearest_id(graph.predecessors(i), assignment, d, terminals)
    return

topnode_map = dict()
assign_to_nearest_id(['DOID:162'], 'DOID:162', topnode_map, topnodes_docancerslim)
cancer_counts = Counter()

cancer_associations = do_sourced.select(
    lambda x: topnode_map.get(x.disease.id, False))

for association in cancer_associations:
    disease_id = association.disease.id
    topnode_id = topnode_map[disease_id]
    cancer_counts[topnode_id] += 1

cancer_associations.report_groups(do_sourced)
other_associations = do_sourced - cancer_associations
other_associations.report_groups(do_sourced)
len(cancer_counts)
for k, v in cancer_counts.most_common(48):
    print(f'{id_to_name[k]}: {v}')
    
def write_disease_counts(file_handle, disease_counter):
    writer = csv.writer(file_handle)
    s = sum(disease_counter.values())
    writer.writerow(['DOID', 'Disease Name', 'Interpretations', 'Percentage'])
    for k, v in disease_counter.most_common(len(disease_counter)):
        writer.writerow([k, id_to_name[k], v, '{:.2%}'.format(v/s)])
with open('out/interpretation_disease_topnode_counts.csv', 'w') as f:
    write_disease_counts(f, cancer_counts)
# Benign

benign_id = 'DOID:0060072'
benign = dict()

assign_to_id(graph.predecessors(benign_id), benign_id, benign)
benign[benign_id] = benign_id

benign_associations = other_associations.select(
    lambda x: benign.get(x.disease.id, False))

other_associations = other_associations - benign_associations
benign_associations.report_groups(do_sourced)
# pre-malignant

premalignant_id = 'DOID:0060071'
premalignant = dict()

assign_to_id(graph.predecessors(premalignant_id), premalignant_id, premalignant)
premalignant[premalignant_id] = premalignant_id

premalignant_associations = other_associations.select(
    lambda x: premalignant.get(x.disease.id, False))

premalignant_associations.report_groups(do_sourced)
# Make data

common_cancers = list(filter(lambda x: x[0] != "DOID:162", cancer_counts.most_common(6)))

interpretation_group_names=[id_to_name[x[0]] for x in common_cancers] + ['other cancers'] + ['benign', 'other disease']
common_cancer_values = [x[1] for x in common_cancers]
interpretation_group_sizes=common_cancer_values + [sum(cancer_counts.values()) - sum(common_cancer_values),
                                      len(benign_associations),
                                      len(other_associations)]

incidence_by_topnode = Counter()
mortality_by_topnode = Counter()
p = (len(benign_associations) + len(other_associations)) / len(do_sourced)
with open(DATA_ROOT / 'Cancer Incidence and Mortality 2018.csv', 'r') as f:
    reader = csv.DictReader(f)
    for row in reader:
        incidence_by_topnode[topnode_map[row['DOID']]] += int(row['New Cases'])
        mortality_by_topnode[topnode_map[row['DOID']]] += int(row['Estimated Deaths'])
def select_by_percent(counter, percent=5):
    assert percent <= 100
    s = sum(counter.values())
    out = Counter()
    for k, v in counter.most_common(len(counter)):
        if k == "DOID:162":
            continue
        p = v/s
        if p*100 < percent:
            break
        out[k] = v
        print(f'{id_to_name[k]}: {p}')
    return out
with open('out/NCI_disease_topnode_counts.csv', 'w') as f:
    writer = csv.writer(f)
    s1 = sum(incidence_by_topnode.values())
    s2 = sum(mortality_by_topnode.values())
    writer.writerow(
        ['DOID', 'Disease Name', 
         'Estimated New Cases, 2018, US', 'Percentage',
         'Estimated Deaths, 2018, US', 'Percentage'
        ])
    for k, v in incidence_by_topnode.most_common(len(incidence_by_topnode)):
        v2 = mortality_by_topnode[k]
        writer.writerow(
            [k, id_to_name[k], 
             v, '{:.2%}'.format(v/s1),
             v2, '{:.2%}'.format(v2/s2)
            ])
prevalent_incidence = select_by_percent(incidence_by_topnode)
prevalent_mortality = select_by_percent(mortality_by_topnode)
len(incidence_by_topnode)
incidence_group_names = [id_to_name[x] for x in prevalent_incidence] + ['other cancers']
incidence_group_sizes = list(prevalent_incidence.values()) + [
                            sum(incidence_by_topnode.values()) - sum(prevalent_incidence.values())
                        ]
mortality_group_names = [id_to_name[x] for x in prevalent_mortality] + ['other cancers']
mortality_group_sizes = list(prevalent_mortality.values()) + [
                            sum(mortality_by_topnode.values()) - sum(prevalent_mortality.values())
                        ]
RADIUS=1

greens = plt.cm.Greens
grays = plt.cm.Greys
blues = plt.cm.Blues
reds = plt.cm.Reds

colors = [
    greens(.85),
    greens(.7),
    greens(.55),
    greens(.4),
    greens(.25),
    greens(.1),
    grays(.5),
    grays(.25)
]
fig1, ax1 = plt.subplots()
pie, _ = ax1.pie(interpretation_group_sizes, colors=colors, radius=RADIUS,
        labels=interpretation_group_names)
ax1.axis('equal')
plt.setp(pie, edgecolor='white')

# plt.savefig(str(FIGPATH / 'disease_interpretations.pdf'))
plt.show()
colors = [
    blues(.85),
    blues(.75 * 5/6 + .1),
    blues(.75 * 4/6 + .1),
    blues(.75 * 3/6 + .1),
    blues(.75 * 2/6 + .1),
    blues(.75 * 1/6 + .1),
    blues(.1)
]

fig1, ax1 = plt.subplots()
pie, _ = ax1.pie(incidence_group_sizes, colors=colors, radius=RADIUS,
        labels=incidence_group_names)

ax1.axis('equal')
plt.setp(pie, edgecolor='white')

# plt.savefig(str(FIGPATH / 'disease_incidence.pdf'))
plt.show()
colors = [
    reds(.85),
    reds(.7),
    reds(.55),
    reds(.4),
    reds(.25),
    reds(.1)
]

fig1, ax1 = plt.subplots()
pie, _ = ax1.pie(mortality_group_sizes, colors=colors, radius=RADIUS,
        labels=mortality_group_names)

ax1.axis('equal')
plt.setp(pie, edgecolor='white')

# plt.savefig(str(FIGPATH / 'disease_mortality.pdf'))
plt.show()
b = benign_associations[0]
print(b.description)
print(b.disease)
print(b.disease.id)
print(b.source)
# MAKE A DISEASE/GENE PLOT: heatmap of ubiquitous gene x disease, heat = tier 1 evidence<jupyter_output><empty_output><jupyter_text>Export disease counts for supplementary table<jupyter_code># d = cancer_associations.get_element_by_source('disease')
# d2 = benign_associations.get_element_by_source('disease')
d = defaultdict(Counter)
d2 = defaultdict(Counter)
d3 = defaultdict(Counter)
for association in cancer_associations:
    disease = association.disease
    source = association.source
    d[source][disease] += 1
for association in benign_associations:
    disease = association.disease
    source = association.source
    d2[source][disease] += 1
for association in other_associations:
    disease = association.disease
    source = association.source
    d3[source][disease] += 1
sorted(d.keys())
diseases = set()
cancer_diseases = set()
benign_diseases = set()
other_diseases = set()
for x in d.values():
    cancer_diseases.update(x)
    diseases.update(x)
for x in d2.values():
    benign_diseases.update(x)
    diseases.update(x)
for x in d3.values():
    other_diseases.update(x)
    diseases.update(x)<jupyter_output><empty_output><jupyter_text>with open(FIGPATH / 'Data' / 'disease_counts.csv', 'w') as f:    header = ['disease', 'doid', 'TopNode_disease', 'TopNode_doid'] + sorted(d.keys())    writer = csv.DictWriter(f, fieldnames=header)    writer.writeheader()    counts = dict()    for s, v in d.items():        counts[s] = Counter(v)    for s, v in d2.items():        counts[s].update(Counter(v))    for s, v in d3.items():        counts[s].update(Counter(v))    for disease in diseases:        if disease in cancer_diseases:            tn_id = topnode_map[disease.id]            tn_name = id_to_name[tn_id]        elif disease in benign_diseases:            tn_id = benign[disease.id]            tn_name = id_to_name[tn_id]        elif disease in other_diseases:            tn_id = None            tn_name = 'other'        else:            raise ValueError        if not disease.id.startswith('DOID'):            continue        try:            row = {                'disease': id_to_name[disease.id],                'doid': disease.id,                'TopNode_disease': tn_name,                'TopNode_doid': tn_id            }        except KeyError:            print(f'Failed to find name for {disease.id}: {disease.name}')            continue        for s in d:            row[s] = counts[s][disease]        writer.writerow(row)  Clinical actionability improvement<jupyter_code># Stacked bar (Actionability type):
#  Group 1: Average actionability, variant only (SD whiskers?)
#  Group 2: Aggregate actionability, variant only (p-value bar?; with narrow search)
#  Group 3: Aggregate actionability, variant only (p-value bar?; with broad search)
#  Group 4: Average actionability, variant + disease
#  Group 5: Aggregate actionability, variant + disease (p-value bar?; with narrow search)
#  Group 6: Aggregate actionability, variant + disease (p-value bar?; with broad search)
# Actionability by disease type
from collections import defaultdict
genie_search_results[0]['query']
genie_features_by_patient['GENIE-NKI-01CH']
genie_search_results_by_query = defaultdict(list)
for result in genie_search_results:
    genie_search_results_by_query[result['query']].append(result)<jupyter_output><empty_output><jupyter_text>Interpretations and Actionability<jupyter_code>tier1_disease_gene = Counter()
all_disease_gene = Counter()
for association in core_vdb:
    try:
        disease = id_to_name[topnode_map[association.disease.id]]
    except (KeyError, AttributeError):
        continue
    level = association.evidence_level
    for gene in association.genes:
        k = (gene.gene_symbol, disease)
        if level in ['A', 'B']:
            tier1_disease_gene[k] += 1
        all_disease_gene[k] += 1
(gene, disease), count = tier1_disease_gene.most_common(30)[0]
import seaborn as sns
f, ax = plt.subplots(figsize=(3.75, 1.5))
genes = Counter()
diseases = interpretation_group_names[:5]
i = 0
for (gene, disease), count in all_disease_gene.most_common(len(all_disease_gene)):
    i += 1
    if disease not in diseases:
        continue
    genes[gene] += count
    
genes = [x[0] for x in genes.most_common(15)]
        
heat_array_all = np.zeros((len(diseases),len(genes)))
for i in range(len(genes)):
    for j in range(len(diseases)):
        heat_array_all[j][i] = all_disease_gene.get((genes[i], diseases[j]), 0)

sns.heatmap(heat_array_all, xticklabels=genes, yticklabels=diseases, robust=True, cmap='Greens')

# f.savefig(str(FIGPATH / 'misc_figures' / 'all_disease_gene_heatmap.pdf'), format='pdf')
# genes = [x.gene_symbol for x in with_genes.element_by_source_stats('genes')['ubiquitous']]
f, ax = plt.subplots(figsize=(3.75, 1.5))
genes = Counter()
diseases = interpretation_group_names[:5]
i = 0
for (gene, disease), count in tier1_disease_gene.most_common(len(tier1_disease_gene)):
    i += 1
    if disease not in diseases:
        continue
    genes[gene] += count
    
genes = [x[0] for x in genes.most_common(15)]
        
heat_array_tier1 = np.zeros((len(diseases),len(genes)))
for i in range(len(genes)):
    for j in range(len(diseases)):
        heat_array_tier1[j][i] = tier1_disease_gene.get((genes[i], diseases[j]), 0)


sns.heatmap(heat_array_tier1, xticklabels=genes, yticklabels=diseases, robust=True, cmap='Greens')

# f.savefig(str(FIGPATH / 'misc_figures' / 'tier1_disease_gene_heatmap.pdf'), format='pdf')
all_df = pd.DataFrame.from_dict(all_disease_gene, orient='index')
tier1_df = pd.DataFrame.from_dict(tier1_disease_gene, orient='index')
merged_df = all_df.merge(tier1_df, left_index=True, right_index=True, how='outer').fillna(0)
merged_df.columns = ['all_interpretations', 'tier1_interpretations']
sns.set_style("ticks")
f, ax = plt.subplots(figsize=(5, 4))

ax = sns.regplot(x='all_interpretations', y='tier1_interpretations', fit_reg=False, data=np.log2(merged_df + 1),
                 marker='.'
                )
tdf = np.log2(merged_df + 1)
tdf.loc[(tdf['all_interpretations'] > 8) & (tdf['tier1_interpretations'] < 4)]
with open(DATA_ROOT / 'GENIE_v3' / 'data_clinical_sample_3.0.0.txt') as f, open(OUTPATH / 'genie_mapping.csv', 'w') as out_f:
    for _ in range(4):
        f.readline()    # get past info lines
    reader = csv.DictReader(f, delimiter="\t")
    writer = csv.writer(out_f)
    header = ['PATIENT_ID',
              'SAMPLE_ID',
              'ONCOTREE_CODE',
              'SPECIFIC_CANCER_TYPE',
              'CANCER_TYPE',
              'DOID',
              'DO_NAME',
              'TOPNODE_DOID',
              'TOPNODE_NAME'
             ]
    writer.writerow(header)
    for row in reader:
        patient = row['PATIENT_ID']
        sample = row['SAMPLE_ID']
        oncotree_code = row['ONCOTREE_CODE'].upper()
        cancer_type_detailed = row['CANCER_TYPE_DETAILED']
        cancer_type = row['CANCER_TYPE']
        doid = sample_to_doid[sample]
        do_name = id_to_name.get(doid, None)
        topnode_id = topnode_map.get(doid, None)
        topnode_name = id_to_name.get(topnode_id, None)
        writer.writerow([
                patient,
                sample,
                oncotree_code,
                cancer_type_detailed,
                cancer_type,
                doid,
                do_name,
                topnode_id,
                topnode_name
            ])
dd_score = dict()

def disease_dist(doid_1, doid_2):
    try:
        key = tuple(sorted([doid_1, doid_2]))
    except TypeError:
        return -1
    if key in dd_score:
        return dd_score[key]
    best_score = None
    queue = [(doid_1, doid_2, 0, False), (doid_2, doid_1, 0, False)]
    while queue:
        current, target, distance, topnode_hit = queue.pop()
        if current == target:
            if best_score is None or distance < best_score:
                best_score = distance
        else:
            if topnode_hit:
                distance += 1
                topnode_hit = False
            if current in topnode_map and current == topnode_map[current]:
                topnode_hit = True
            try:
                for successor in graph.successors(current):
                    queue.append((successor, target, distance, topnode_hit))
            except networkx.NetworkXError:
                pass
    if best_score is None:
        best_score = -1
    dd_score[key] = best_score
    return best_score
SOURCES = tuple(sorted(core_vdb.sources))
patient_actionability = dict()
disease_actionability_count = Counter()
sample_disease_count = Counter()
disease_feature_count = Counter()
patients_with_topnode = set()
cgi_diseases_with_patient_match = defaultdict(set)
ALL_PATIENT_COUNT = len(patient_to_samples)
for patient in featured_patients:
    samples = patient_to_samples[patient]
    actionable = np.zeros(40) # [0..5]: sources (exact), 6: combined, 7: deprecated. 
                              # +8 for +disease +16 for +dis/tier1
                              # +24 for broad variant, +32 for broad disease
    for sample in samples:
        sample_disease = sample_to_doid[sample]
        features = genie_features_by_sample[sample]
        evidence_level = None
        for feature in features:
            results = genie_search_results_by_query.get(feature, [])
            for result in results:
                best = result['best_match']
                dd = disease_dist(result['association'].disease.id, sample_disease)
                level = result['association']['association']['evidence_label']
                tier1 = level in ['A', 'B']
                idx = SOURCES.index(result['association']['source'])
                if best['type'] in ['exact', 'positional']:
                    actionable[idx] = 1
                    actionable[6] = 1
                    if result['association'].source == 'cgi' and sample_disease:
                            cgi_diseases_with_patient_match[patient].add((result['association'].disease.id, sample_disease))
                    if dd == 0:
                        actionable[14] = 1
                        actionable[idx + 8] = 1
                        if tier1:
                            actionable[22] = 1
                            actionable[idx + 16] = 1
                        if evidence_level is None or level < evidence_level:
                            evidence_level = level
                actionable[idx + 24] = 1
                actionable[30] = 1
                if dd >= 0:
                    actionable[idx + 32] = 1
                    actionable[38] = 1
        try:
            top_disease = id_to_name[topnode_map[sample_disease]]
            sample_disease_count[top_disease] += 1
            disease_feature_count[top_disease] += len(features)
        except KeyError:
            top_disease = None
        if evidence_level is not None and top_disease is not None:
            disease_actionability_count[(top_disease, evidence_level)] += 1
        if top_disease is not None:
            patients_with_topnode.add(sample_to_patient[sample])
    patient_actionability[patient] = actionable
    
actionability_grid = np.array(list(patient_actionability.values()))
actionability_sum = actionability_grid.sum(axis=0)<jupyter_output><empty_output><jupyter_text>investigating CGI discrepancy<jupyter_code># len(cgi_diseases_with_patient_match)
patient_disease_match = dict()
for patient, matches in cgi_diseases_with_patient_match.items():
    best_score = -1
    for cgi_diseases in matches:
        score = disease_dist(*cgi_diseases)
        if score >= 0:
            if best_score == -1:
                best_score = score
            elif score < best_score:
                best_score = score
    if best_score < 0:
        patient_disease_match[patient] = {'total': 0}
        patient_diseases = set([x[1] for x in matches])
        for disease in patient_diseases:
            patient_disease_match[patient][disease] = Counter()
        for cgi_diseases in matches:
            patient_disease_match[patient][cgi_diseases[1]][cgi_diseases[0]] += 1
            patient_disease_match[patient]['total'] += 1
# for patient in sorted(patient_disease_match, key=lambda x: patient_disease_match[x]['total'], reverse=True):
#     print(patient, patient_disease_match[patient]['total'])
disease_counter = Counter()
for patient in patient_disease_match:
    for disease in patient_disease_match[patient]:
        disease_counter[disease] += 1
disease_counter.most_common(30)
cgi_disease_counter = Counter()
for patient in patient_disease_match:
    if 'DOID:3008' in patient_disease_match[patient]:
        for cgi_disease in patient_disease_match[patient]['DOID:3008']:
            cgi_disease_counter[cgi_disease] += 1
cgi_disease_counter.most_common(30)<jupyter_output><empty_output><jupyter_text>disease match analysis<jupyter_code>x = (0, 2, 4, 6, 8, 10) * 4
y = list()
for i in (6, 4, 2, 0):
    y.extend((i,) * 6)
    
levels = ['A','B','C','D']

proportion = np.zeros((4,6))
for i, level in enumerate(levels):
    for j, disease in enumerate(diseases):
        proportion[i][j] = disease_actionability_count[(disease, level)] / sample_disease_count[disease]
proportion[:,5] = (.5,.25,.1,.01)

f, ax = plt.subplots(figsize=(5, 4))
ax.set_xlim(-1,11)
ax.set_ylim(-1,7)
plt.scatter(x, y, s=proportion*4000)
# f.savefig(str(FIGPATH / 'misc_figures' / 'disease_actionability.pdf'), format='pdf')
diseases
interpretable_disease_counts = np.zeros((4,2))
total_counts = np.zeros(2)
genie_topnode_diseases = set([
        x[0] for x in disease_actionability_count
    ])
for disease in genie_topnode_diseases:
    for i, level in enumerate(levels):
        if disease in diseases:
            j = 0
        else:
            j = 1
        interpretable_disease_counts[i][j]   += disease_actionability_count[(disease, level)]
    total_counts[j] += sample_disease_count[disease]

test_matrix = np.zeros((2,2))
test_matrix[0,:] = interpretable_disease_counts.sum(axis=0)
test_matrix[1,:] = total_counts - test_matrix[0,:]
print(test_matrix[0,:] / total_counts)
print(test_matrix)

stats.fisher_exact(test_matrix)
interpretable_disease_counts.sum(axis=0)
test_matrix = np.zeros((2,2))
test_matrix[0,:] = interpretable_disease_counts[:2,:].sum(axis=0)
test_matrix[1,:] = total_counts - test_matrix[0,:]
print(test_matrix[0,:] / total_counts)
print(test_matrix)

result = stats.fisher_exact(test_matrix)
print(result)
fig, ax = plt.subplots(figsize=(7.5,3))

counts = actionability_sum
groups = ['Variant', 'Variant + Disease', 'Variant + Disease + Tier I']
subs = SOURCES + ('aggregate',)
width = 0.15
plot_elements = list()
resource_ind = np.arange(3)

for i, source in enumerate(subs):
    x = np.array([actionability_sum[i + 24], actionability_sum[i + 32], 0]) / ALL_PATIENT_COUNT
    ind = np.array((i, i+12, i+24)) * width
    p = ax.bar(ind, x, width, label=source, color='black')
    plot_elements.append(p)

for i, source in enumerate(subs):
    x = np.array([actionability_sum[i], actionability_sum[i + 8], actionability_sum[i + 16]]) / ALL_PATIENT_COUNT
    ind = np.array((i, i+12, i+24)) * width
    p = ax.bar(ind, x, width, label=source)
    plot_elements.append(p)

ax.set_xticks(ind - 3*width)
# ax.set_xticklabels(groups, rotation='vertical')
handles, labels = ax.get_legend_handles_labels()
half_idx = len(labels) // 2
ax.legend(handles[half_idx:], labels[half_idx:], title='Search Type')
plt.ylabel('% Cohort with Interpretations')

plt.show()
# fig.savefig(str(FIGPATH / 'misc_figures' / 'genie_actionability.pdf'), format='pdf')
(actionability_sum / ALL_PATIENT_COUNT).reshape((5,8))
# Average individual KB variant matching, exact searching
(actionability_sum[:6] / ALL_PATIENT_COUNT).mean()
# Average individual KB variant+disease matching, exact searching
(actionability_sum[8:14] / ALL_PATIENT_COUNT).mean()
# Average individual KB variant+disease+tier I matching, exact searching
(actionability_sum[16:22] / ALL_PATIENT_COUNT).mean()
# Average individual KB variant matching, broad searching
(actionability_sum[24:30] / ALL_PATIENT_COUNT).mean()
# Average individual KB variant+disease matching, broad searching
(actionability_sum[32:38] / ALL_PATIENT_COUNT).mean()
c = defaultdict(Counter)
for a in core_vdb:
    c[a.source][a.evidence_level] += 1
c
for source, counts in c.items():
    tier1 = (counts['A'] + counts['B']) / sum(counts.values())
    print("{}: Tier 1 is {:.1%} of total.".format(source, tier1))
s = len(patients_with_topnode)
s / ALL_PATIENT_COUNT
len(sample_oncotree_code)
code_patient_count = Counter()
for patient, samples in patient_to_samples.items():
    patient_codes = set()
    for sample in samples:
        patient_codes.add(sample_oncotree_code[sample])
    for code in patient_codes:
        code_patient_count[code] += 1
patient_count = len(patient_to_samples)
for code, count in code_patient_count.most_common():
    f = count / patient_count
    print('{}: {:.1%}, {}'.format(code, f, count))
cancer_rank = {pair[0]: rank for rank, pair in enumerate(
        code_patient_count.most_common()
    ) if pair[0] not in patched_doids} # Exclude manual patching for analysis
mapped = list()
unmapped = list()
for code, rank in cancer_rank.items():
    doids = oncotree_to_doids.get(code, False)
    if doids:
        mapped.append(rank)
    else:
        unmapped.append(rank)
print(mapped)
print(unmapped)
print(len(mapped) / (len(unmapped) + len(mapped)))
stats.mannwhitneyu(mapped, unmapped, alternative='two-sided')
cancer_rank
len(unmapped)
np.mean(mapped)
np.mean(unmapped)<jupyter_output><empty_output><jupyter_text>Interpretation gene intersection search<jupyter_code>a = core_vdb[0]
f = genie_features[0]
e_level = defaultdict(dict)
gene_diseases = defaultdict(lambda: defaultdict(set))
for association in core_vdb:
    source = association.source
    for gene in association.genes:
        try:
            disease = association.disease.id
        except AttributeError:
            continue
        key = (gene.symbol, disease)
        gene_diseases[source][gene.symbol].add(disease)
        gene_diseases['aggregate'][gene.symbol].add(disease)
        current = e_level[source].get(key, None)
        new = association.evidence_level
        if current is None or new < current:
            e_level[source][key] = new
        current = e_level['aggregate'].get(key, None)
        if current is None or new < current:
            e_level['aggregate'][key] = new
patient_gene_actionability = dict()
for patient in featured_patients:
    samples = patient_to_samples[patient]
    actionable = np.zeros(24) # [0..5]: sources (exact), 6: combined, 7: deprecated. 
                              # +8 for +disease +16 for +dis/tier1
    aggregate_gene_diseases = gene_diseases['aggregate']
    aggregate_e_level = e_level['aggregate']
    for sample in samples:
        sample_disease = sample_to_doid[sample]
        features = genie_features_by_sample[sample]
        evidence_level = None
        for feature in features:
            feature_gene = feature.gene_symbol
            if feature_gene not in aggregate_gene_diseases:
                continue
            actionable[6] = 1
            for interpretation_disease in aggregate_gene_diseases[feature_gene]:
                if disease_dist(interpretation_disease, sample_disease) >= 0:
                    actionable[14] = 1
                    if aggregate_e_level[(feature_gene, interpretation_disease)] in ['A', 'B']:
                        actionable[22] = 1
                        break
            for i, source in enumerate(SOURCES):
                source_e_level = e_level[source]
                source_gene_diseases = gene_diseases[source]
                if feature_gene not in source_gene_diseases:
                    continue
                actionable[i] = 1
                if not actionable[14]:
                    continue
                for interpretation_disease in source_gene_diseases[feature_gene]:
                    if disease_dist(interpretation_disease, sample_disease) >= 0:
                        actionable[i + 8] = 1
                        if source_e_level[(feature_gene, interpretation_disease)] in ['A', 'B']:
                            actionable[i + 16] = 1
                            break
    patient_gene_actionability[patient] = actionable
    
gene_actionability_grid = np.array(list(patient_gene_actionability.values()))
gene_actionability_sum = gene_actionability_grid.sum(axis=0)
(gene_actionability_sum / ALL_PATIENT_COUNT).reshape((3,8))
fig, ax = plt.subplots(figsize=(7.5,3))

counts = gene_actionability_sum
groups = ['Variant', 'Variant + Disease', 'Variant + Disease + Tier I']
subs = SOURCES + ('aggregate',)
width = 0.15
plot_elements = list()
resource_ind = np.arange(3)

# for i, source in enumerate(subs):
#     x = np.array([gene_actionability_sum[i + 24], gene_actionability_sum[i + 32], 0]) / gene_actionability_grid.shape[0]
#     ind = np.array((i, i+12, i+24)) * width
#     p = ax.bar(ind, x, width, label=source, color='black')
#     plot_elements.append(p)

for i, source in enumerate(subs):
    x = np.array([gene_actionability_sum[i], gene_actionability_sum[i + 8], gene_actionability_sum[i + 16]]) / ALL_PATIENT_COUNT
    ind = np.array((i, i+12, i+24)) * width
    p = ax.bar(ind, x, width, label=source)
    plot_elements.append(p)

ax.set_xticks(ind - 3*width)
# ax.set_xticklabels(groups, rotation='vertical')
handles, labels = ax.get_legend_handles_labels()
# half_idx = len(labels) // 2
half_idx = 0
lgd = ax.legend(handles[half_idx:], labels[half_idx:], title='Search Type',
          bbox_to_anchor=(1,1)
         )
plt.ylabel('% Cohort with Interpretations')

plt.show()
# fig.savefig(str(FIGPATH / 'misc_figures' / 'genie_gene_actionability.pdf'), format='pdf',
#             bbox_extra_artists=(lgd, ), bbox_inches='tight'
#            )<jupyter_output><empty_output><jupyter_text>Quantitative value of harmonization  Genes<jupyter_code>raw_genes = list()<jupyter_output><empty_output><jupyter_text>CGI<jupyter_code>def create_raw_cgi_genes(vdb):
    cgi = vdb.by_source('cgi')
    genes = list()
    for interpretation in cgi:
        for gene in interpretation['raw']['Gene'].split(';'):
            genes.append(gene)
    return genes
cgi_genes = create_raw_cgi_genes(core_vdb)
raw_genes.extend(set(cgi_genes))<jupyter_output><empty_output><jupyter_text>CIViC<jupyter_code>def create_raw_civic_genes(vdb):
    civic = vdb.by_source('civic')
    genes = list()
    for interpretation in civic:
        genes.append(interpretation['raw']['entrez_name'])
    return genes
civic_genes = create_raw_civic_genes(core_vdb)
raw_genes.extend(set(civic_genes))<jupyter_output><empty_output><jupyter_text>JAX-CKB<jupyter_code>def create_raw_jax_genes(vdb):
    jax = vdb.by_source('jax')
    genes = list()
    for interpretation in jax:
        for gene in interpretation['genes']:
            genes.append(gene.strip())
    return genes
jax_genes = create_raw_jax_genes(core_vdb)
raw_genes.extend(set(jax_genes))<jupyter_output><empty_output><jupyter_text>MolecularMatch<jupyter_code>def create_raw_molecularmatch_genes(vdb):
    mm = vdb.by_source('molecularmatch')
    genes = list()
    for interpretation in mm:
        genes2 = interpretation['raw'].get('includeGene0', None)
        if not genes2:
            genes2 = interpretation['raw'].get('includeGene1', None)
        if genes2:
            genes.extend(genes2)
    return genes
mm_genes = create_raw_molecularmatch_genes(core_vdb)
raw_genes.extend(set(mm_genes))<jupyter_output><empty_output><jupyter_text>OncoKB<jupyter_code>def create_raw_okb_genes(vdb):
    okb = vdb.by_source('oncokb')
    genes = list()
    for interpretation in okb:
        gene = interpretation['raw']['clinical']['gene']
        genes.append(gene)
    return genes
okb_genes = create_raw_okb_genes(core_vdb)
raw_genes.extend(set(okb_genes))<jupyter_output><empty_output><jupyter_text>PMKB<jupyter_code>def create_raw_pmkb_genes(vdb):
    pmkb = vdb.by_source('pmkb')
    genes = list()
    for interpretation in pmkb:
        gene = interpretation['raw']['variant']['gene']['name']
        genes.append(gene)
    return genes
pmkb_genes = create_raw_pmkb_genes(core_vdb)
raw_genes.extend(set(pmkb_genes))<jupyter_output><empty_output><jupyter_text>Harmonized genes<jupyter_code>h_genes = list()

for source in core_vdb.sources:
    unique_source_genes = set()
    for interpretation in core_vdb.by_source(source):
        for gene in interpretation.genes:
            unique_source_genes.add(gene.symbol)
    print(f'{source}: {len(unique_source_genes)}')
    h_genes.extend(unique_source_genes)<jupyter_output>molecularmatch: 109
civic: 296
pmkb: 42
oncokb: 44
jax: 107
cgi: 182<jupyter_text>Comparison<jupyter_code>raw_start = len(raw_genes)
raw_uniq = len(set(raw_genes))
h_start = len(h_genes)
h_uniq = len(set(h_genes))
1 - (h_uniq / h_start)
1 - (raw_uniq / raw_start)<jupyter_output><empty_output><jupyter_text>Variants<jupyter_code>raw_features = list()<jupyter_output><empty_output><jupyter_text>CGI<jupyter_code>hgvs_re = re.compile(r'(.*):g.(\d+)(\w+)>(\w+)')

def hgvs_to_coords(hgvs):
    match = hgvs_re.match(hgvs)
    if not match:
        return None
    groups = match.groups()
    return (str(groups[0]), int(groups[1]), int(groups[1]) + len(groups[2]) - 1)

def create_raw_cgi_features(vdb):
    cgi = vdb.by_source('cgi')
    features = list()
    for interpretation in cgi:
        for f in interpretation['raw']['gDNA']:
            if not f:
                continue
            coords = hgvs_to_coords(f)
            if coords:
                features.append(coords)
    return features
            
cgi_features = create_raw_cgi_features(core_vdb)
len(set(cgi_features))
raw_features.extend(set(cgi_features))<jupyter_output><empty_output><jupyter_text>CIViC<jupyter_code>def create_raw_civic_features(vdb):
    civic = vdb.by_source('civic')
    features = list()
    harmonized_features = list()
    for interpretation in civic:
        coordinates = interpretation['raw']['coordinates']
        gf = (
            coordinates['chromosome'],
            coordinates['start'],
            coordinates['stop'],
#             coordinates['variant_bases']
        )
        if not all(gf[:3]):
            continue
        features.append(gf)
    return features
civic_features = create_raw_civic_features(core_vdb)
len(set(civic_features))
raw_features.extend(set(civic_features))<jupyter_output><empty_output><jupyter_text>JAX-CKB There is nothing to do for this resource. Without harmonization or inference routines, cannot ascribe variant names to coordinates.  MolecularMatch<jupyter_code>def create_raw_molecularmatch_features(vdb):
    features = list()
    mm = vdb.by_source('molecularmatch')
    for interpretation in mm:
        for mutation in interpretation['raw']['mutations']:
            try:
                coords = mutation['GRCh37_location'][0] # Take first coord only
            except IndexError:
                continue
            try:
                start = int(coords['start'])
                stop = int(coords['stop'])
            except TypeError:
                continue
            assert start <= stop
            chromosome = coords['chr']
            alt = coords['alt']
            if not all([chromosome, start, stop]):
                continue
            f = (chromosome, start, stop)
            features.append(f)
    return features
    
mm_features = create_raw_molecularmatch_features(core_vdb)
len(set(mm_features))
raw_features.extend(set(mm_features))<jupyter_output><empty_output><jupyter_text>OncoKB All mutations at protein level.  PMKB<jupyter_code>def create_raw_pmkb_features(vdb):
    features = list()
    pmkb = vdb.by_source('pmkb')
    for interpretation in pmkb:
        coordinates = interpretation['raw']['variant']['coordinates']
        for coordinate in coordinates.split(', '):
            chromosome, r = coordinate.split(':')
            start, stop = r.split('-')
            features.append((chromosome, int(start), int(stop)))
    return features
pmkb_features = create_raw_pmkb_features(core_vdb)
len(set(pmkb_features))
raw_features.extend(set(pmkb_features))<jupyter_output><empty_output><jupyter_text>Harmonized features<jupyter_code>h_features = list()

for source in core_vdb.sources:
    unique_source_features = set([x[0] for x in core_vdb.by_source(source).features])
    print(f'{source}: {len(unique_source_features)}')
    h_features.extend(unique_source_features)
len(h_features)<jupyter_output><empty_output><jupyter_text>Comparison<jupyter_code>raw_start = len(raw_features)
raw_uniq = len(set(raw_features))
h_start = len(h_features)
h_uniq = len(set(h_features))
1 - (h_uniq / h_start)
1 - (raw_uniq / raw_start)<jupyter_output><empty_output><jupyter_text>Diseases<jupyter_code>raw_diseases = list()<jupyter_output><empty_output><jupyter_text>CGI<jupyter_code>def create_raw_cgi_diseases(vdb):
    cgi = vdb.by_source('cgi')
    diseases = list()
    for interpretation in cgi:
        for d in interpretation['raw']['Primary Tumor type'].split(';'):
            diseases.append(d.lower())
    return diseases
cgi_diseases = create_raw_cgi_diseases(core_vdb)
raw_diseases.extend(set(cgi_diseases))<jupyter_output><empty_output><jupyter_text>CIViC<jupyter_code>def create_raw_civic_diseases(vdb):
    civic = vdb.by_source('civic')
    diseases = list()
    for interpretation in civic:
        disease = interpretation['raw']['evidence_items'][0]['disease']['display_name']
        diseases.append(disease.lower())
    return diseases
civic_diseases = create_raw_civic_diseases(core_vdb)
raw_diseases.extend(set(civic_diseases))<jupyter_output><empty_output><jupyter_text>JAX-CKB<jupyter_code>def create_raw_jax_diseases(vdb):
    jax = vdb.by_source('jax')
    diseases = list()
    for interpretation in jax:
        disease = interpretation['raw']['indication']['name']
        diseases.append(disease.lower())
    return diseases
jax_diseases = create_raw_jax_diseases(core_vdb)
raw_diseases.extend(set(jax_diseases))<jupyter_output><empty_output><jupyter_text>MolecularMatch<jupyter_code>def create_raw_mm_diseases(vdb):
    mm = vdb.by_source('molecularmatch')
    diseases = list()
    for interpretation in mm:
        disease = interpretation['raw']['includeCondition1'][0]
        diseases.append(disease.lower())
    return diseases
mm_diseases = create_raw_mm_diseases(core_vdb)
raw_diseases.extend(set(mm_diseases))<jupyter_output><empty_output><jupyter_text>OncoKB<jupyter_code>def create_raw_okb_diseases(vdb):
    okb = vdb.by_source('oncokb')
    diseases = list()
    for interpretation in okb:
        disease = interpretation['raw']['clinical']['cancerType']
        diseases.append(disease.lower())
    return diseases
okb_diseases = create_raw_okb_diseases(core_vdb)
raw_diseases.extend(set(okb_diseases))<jupyter_output><empty_output><jupyter_text>PMKB<jupyter_code>def create_raw_pmkb_diseases(vdb):
    pmkb = vdb.by_source('pmkb')
    diseases = list()
    for interpretation in pmkb:
        for tissue in interpretation['raw']['tissues']:
            disease = ' '.join([tissue['name'], interpretation['raw']['tumor']['name']])
            diseases.append(disease.lower())
    return diseases
pmkb_diseases = create_raw_pmkb_diseases(core_vdb)
raw_diseases.extend(set(pmkb_diseases))<jupyter_output><empty_output><jupyter_text>Harmonized diseases<jupyter_code>h_diseases = list()

for source in core_vdb.sources:
    unique_source_diseases = set([x.disease.term for x in core_vdb.by_source(source) if x.disease])
    print(f'{source}: {len(unique_source_diseases)}')
    h_diseases.extend(unique_source_diseases)
len(h_diseases)<jupyter_output><empty_output><jupyter_text>Comparison<jupyter_code>raw_start = len(raw_diseases)
raw_uniq = len(set(raw_diseases))
h_start = len(h_diseases)
h_uniq = len(set(h_diseases))
1 - (h_uniq / h_start)
1 - (raw_uniq / raw_start)<jupyter_output><empty_output><jupyter_text>Drugs<jupyter_code>raw_drugs = list()<jupyter_output><empty_output><jupyter_text>CGI<jupyter_code>def create_raw_cgi_drugs(vdb):
    cgi = vdb.by_source('cgi')
    drugs = list()
    for interpretation in cgi:
        for drug in interpretation['raw']['Drug full name'].split('+'):
            drugs.append(drug.strip().lower())
    return drugs
cgi_drugs = create_raw_cgi_drugs(core_vdb)
raw_drugs.extend(set(cgi_drugs))<jupyter_output><empty_output><jupyter_text>CIViC<jupyter_code>def create_raw_civic_drugs(vdb):
    civic = vdb.by_source('civic')
    drugs = list()
    for interpretation in civic:
        for drug in interpretation['raw']['evidence_items'][0]['drugs']:
            drugs.append(drug['name'].lower())
    return drugs
civic_drugs = create_raw_civic_drugs(core_vdb)
raw_drugs.extend(set(civic_drugs))<jupyter_output><empty_output><jupyter_text>JAX-CKB<jupyter_code>def create_raw_jax_drugs(vdb):
    jax = vdb.by_source('jax')
    drugs = list()
    for interpretation in jax:
        for drug in interpretation['raw']['therapy']['therapyName'].split('+'):
            drugs.append(drug.strip().lower())
    return drugs
jax_drugs = create_raw_jax_drugs(core_vdb)
raw_drugs.extend(set(jax_drugs))<jupyter_output><empty_output><jupyter_text>MolecularMatch<jupyter_code>def create_raw_mm_drugs(vdb):
    mm = vdb.by_source('molecularmatch')
    drugs = list()
    for interpretation in mm:
        try:
            for drug in interpretation['raw']['includeDrug1']:
                for drug2 in drug.split('+'):
                    drugs.append(drug2.strip().lower())
        except KeyError:
            continue
    return drugs
mm_drugs = create_raw_mm_drugs(core_vdb)
raw_drugs.extend(set(mm_drugs))<jupyter_output><empty_output><jupyter_text>OncoKB<jupyter_code>def create_raw_okb_drugs(vdb):
    okb = vdb.by_source('oncokb')
    drugs = list()
    for interpretation in okb:
        for drug in interpretation['raw']['clinical']['drug'].split(','):
            for drug2 in drug.strip().split('+'):
                drugs.append(drug2.strip().lower())
    return drugs
okb_drugs = create_raw_okb_drugs(core_vdb)
raw_drugs.extend(set(okb_drugs))<jupyter_output><empty_output><jupyter_text>PMKB PMKB does not provide drug fields  Harmonized drugs<jupyter_code>h_drugs = list()

for source in core_vdb.sources:
    unique_source_drugs = set()
    for interpretation in core_vdb.by_source(source):
        for drug in interpretation.drugs:
            unique_source_drugs.add(drug.term)
    print(f'{source}: {len(unique_source_drugs)}')
    h_drugs.extend(unique_source_drugs)<jupyter_output>molecularmatch: 110
civic: 313
pmkb: 0
oncokb: 77
jax: 542
cgi: 200<jupyter_text>Comparison<jupyter_code>raw_start = len(raw_drugs)
raw_uniq = len(set(raw_drugs))
h_start = len(h_drugs)
h_uniq = len(set(h_drugs))
1 - (h_uniq / h_start)
1 - (raw_uniq / raw_start)
Counter([x['raw']['clinicalSignificance'] for x in core_vdb.by_source('molecularmatch')])
Counter([x['raw']['clinicalSignificance'] for x in vdb.by_source('molecularmatch')])<jupyter_output><empty_output><jupyter_text>Export for Somatic Reference Sample project<jupyter_code>import csv
with open('out/srs_export.csv', 'w') as f:
    writer = csv.writer(f)
    writer.writerow([
            'chromosome',
            'start',
            'stop',
            'ref',
            'alt',
            'feature_label',
            'evidence_level',
            'drugs',
            'disease_context',
            'disease_id',
            'interpretation_source',
            'pmids',
        ])
    for a in core_vdb:
        pmids = '|'.join([str(x.pmid) for x in a.publications if x.pmid])
        drugs = '|'.join([str(x) for x in a.drugs if x])
        if a.evidence_level not in ['A', 'B']:
            continue
        for drug in a.drugs:
            assert '|' not in str(drug)
        for feature in a.features:
            try:
                feature_name = feature.name
            except AttributeError:
                feature_name = ''
            try:
                disease_name = a.disease.name
            except AttributeError:
                disease_name = a['association']['phenotype'].get('description', '')
            try:
                disease_id = a.disease.id
            except AttributeError:
                continue
            out = [feature.chromosome, feature.start, feature.end, 
                   feature.ref, feature.alt, feature_name, 
                   a.evidence_level, drugs,
                   disease_name, disease_id, a.source, pmids]
            if not (feature.ref and feature.alt):
                continue
            writer.writerow(out)<jupyter_output><empty_output><jupyter_text>Troubleshooting<jupyter_code>x = core_vdb[0]
f = x.features[0]
x.disease.id
len(unfiltered_patients_with_variants)
len(unfiltered_patients_with_variants) / len(patient_to_samples)
brca.report_groups()<jupyter_output>0 total associations
{'fraction_non_alphanumeric': 0.12846325279846021, 'fraction_numerical': 0.0181076837360077, 'mean_word_length': 3.9556350401606424, 'pattern_counts': {'":': 17, '<': 246, '<?xml version=': 1, '>': 249, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '59', 'dir': 'jupyter-structured-clean-dedup', 'id': '165897', 'license': 'MIT', 'n_tokens_mistral': 29246, 'n_tokens_neox': 27285, 'n_words': 5360, 'path': 'analysis.ipynb', 'repo_name': 'ahwagner/vicckb', 'size': '78972'}