import os
os.system('annotate_variation.pl -buildver hg19 -downdb -webfrom annovar clinvar_20150629 humandb/')

"""
Explanation: <center> Reproducible and Shareable Quantifications of Pathogenicity </center>
<center>Manrai, Wang, Patel, Kohane.
Pacific Symposium on Biocomputing 2016 Submission</center>
1. Loading and Cleaning Pathogenicity Assertions for all Diseases from ClinVar
The ClinVar database [http://www.ncbi.nlm.nih.gov/clinvar/] contains pathogenicity assertions for genetic variants across a broad set of diseases. First, let's use the ANNOVAR tool [http://annovar.openbioinformatics.org/en/latest/] to download these annotations locally. After running the following command, a new directory called "humandb" will be created to store the ClinVar annotations, and temporary files will be removed. Exit code 0 indicates success.
End of explanation
"""


with open('humandb/hg19_clinvar_20150629.txt') as infile:
    first_five_lines = [next(infile) for i in range(5)]
print first_five_lines

"""
Explanation: Let's look at a few lines of this file:
End of explanation
"""


import pandas as pd
raw_clinvar = pd.read_table('humandb/hg19_clinvar_20150629.txt', header=None)
raw_clinvar.columns = ['Chromosome','Start','Stop','Ref','Alt','ClinSig']
print 'There are %s rows in this version of ClinVar.' %raw_clinvar.shape[0]

distinct_variants = set(zip(raw_clinvar.Chromosome,raw_clinvar.Start,raw_clinvar.Stop,raw_clinvar.Ref,raw_clinvar.Alt))
num_distinct_variants = len(distinct_variants)
print 'There are %s distinct variants in this version of ClinVar.' %num_distinct_variants

"""
Explanation: We can load this file more cleanly with the pandas library (http://pandas.pydata.org/):
End of explanation
"""


raw_clinvar.head(5)

"""
Explanation: Taking a closer look at raw_clinvar, we see that the clinical significance data is squished together in one field (ClinSig):
End of explanation
"""


assertion = [x.split(';')[0].split('=')[1] for x in raw_clinvar.ClinSig]
disease = [x.split(';')[1].split('=')[1] for x in raw_clinvar.ClinSig]
accession = [x.split(';')[3].split('=')[1] for x in raw_clinvar.ClinSig]
clinvar = raw_clinvar.drop('ClinSig',1)
clinvar['assertion'], clinvar['disease'], clinvar['accession'] = assertion, disease, accession
clinvar.head(25)

"""
Explanation: Separating this field:
End of explanation
"""


import re

accession = clinvar.accession.map(lambda x: re.split(r'[,|]+',x) if ('|' in x or ',' in x) else [x])
multi_rows = [i for i, e in enumerate(accession) if len(e) > 1]
multi = clinvar.ix[multi_rows,:]
multi.head(5)

"""
Explanation: Notice that some of variants contain multiple assertions, separated by the ',' or '|' characters. We can extract just these rows as follows:
End of explanation
"""


assertion = clinvar.assertion.map(lambda x: re.split(r'[,|]+',x) if ('|' in x or ',' in x) else [x])
disease = clinvar.disease.map(lambda x: re.split(r'[,|]+',x) if ('|' in x or ',' in x) else [x])
accession = clinvar.accession.map(lambda x: re.split(r'[,|]+',x) if ('|' in x or ',' in x) else [x])
v0 = clinvar.Chromosome.tolist()
v1 = clinvar.Start.tolist()
v2 = clinvar.Stop.tolist()
v3 = clinvar.Ref.tolist()
v4 = clinvar.Alt.tolist()
variants = zip(v0,v1,v2,v3,v4)
flat = []
for i, v in enumerate(variants):
    num_assertions = len(accession[i])
    v = list(v)
    if num_assertions > 1:  
        flat += [v+[assertion[i][j],disease[i][j],accession[i][j]] for j in range(num_assertions)]
    else:
        flat += [v+assertion[i]+disease[i]+accession[i]]
flat_clinvar = pd.DataFrame(flat)
flat_clinvar.columns = ['Chromosome','Start','Stop','Ref','Alt','Assertion','Disease','Accession']
flat_clinvar.head(5)

"""
Explanation: We can flatten the clinvar data so there is one <b>assertion</b> per line:
End of explanation
"""


flat_clinvar = flat_clinvar.drop_duplicates()
print "There are %s distinct assertions in this version of Clinvar." %flat_clinvar.shape[0]

"""
Explanation: Let's deduplicate:
End of explanation
"""


flat_clinvar.to_csv(path_or_buf='flat_clinvar.tsv',sep='\t',index=False,
                   columns=['Chromosome','Start','Stop','Ref','Alt','Assertion',
                               'Disease','Accession'])

"""
Explanation: and output to file:
End of explanation
"""


pathogenic_SNPs = flat_clinvar[(flat_clinvar.Start == flat_clinvar.Stop) & (flat_clinvar.Assertion == 'pathogenic') &
                              (flat_clinvar.Alt.str.len() == 1) & (flat_clinvar.Ref.str.len() == 1)]
print 'There are %s pathogenic SNP assertions in this version of ClinVar.' % pathogenic_SNPs.shape[0]

"""
Explanation: To reduce to Pathogenic SNPs, we can simply run:
End of explanation
"""


pathogenic_SNPs.to_csv(path_or_buf='flat_clinvar_pathogenic_SNPs.tsv',sep='\t',index=False,
                       columns=['Chromosome','Start','Stop','Ref','Alt','Assertion',
                               'Disease','Accession'])

"""
Explanation: Let's output these to file:
End of explanation
"""


os.system('annotate_variation.pl -downdb -webfrom annovar -build hg19 esp6500si_ea humandb/')
os.system('annotate_variation.pl -downdb -webfrom annovar -build hg19 esp6500si_aa humandb/')
os.system('annotate_variation.pl -downdb -webfrom annovar -build hg19 esp6500si_all humandb/')

"""
Explanation: 2a. Pulling Allele Frequency Data
2a.1. NHLBI ESP
We first retrieve frequency data from the NHLBI ESP [http://evs.gs.washington.edu/EVS/]:
End of explanation
"""


infile_name = 'flat_clinvar_pathogenic_SNPs.tsv'
outfile_name = 'ESP'
os.system('annotate_variation.pl -filter -dbtype esp6500si_all -build hg19 -out %s %s humandb/' % (outfile_name,infile_name))

ESP_data = pd.read_table('ESP.hg19_esp6500si_all_dropped',
                     header=None)
ESP_data = ESP_data.drop(0,1)
ESP_data.columns = ['ESP_Overall_Frequency','Chromosome','Start','Stop','Ref','Alt',
                'Assertion','Disease','Accession']
print "There are %s variant-disease assertions with ESP frequency data available." %ESP_data.shape[0]
ESP_data.head(5)

"""
Explanation: We can use these data to filter clinical variants:
End of explanation
"""


os.system('annotate_variation.pl -downdb -webfrom annovar -build hg19 exac03 humandb/')

infile_name = 'flat_clinvar_pathogenic_SNPs.tsv'
outfile_name = 'ExAc'
os.system('annotate_variation.pl -filter -build hg19 -dbtype exac03 -out %s %s humandb/' % (outfile_name,infile_name))

ExAc_data = pd.read_table('ExAc.hg19_exac03_dropped',
                     header=None)
ExAc_data = ExAc_data.drop(0,1)
ExAc_data.columns = ['ExAc_Overall_Frequency','Chromosome','Start','Stop','Ref','Alt',
                'Assertion','Disease','Accession']
print "There are %s variant-disease assertions with ExAc frequency data available." %ExAc_data.shape[0]
ExAc_data.head(5)

"""
Explanation: 2a.2. Broad ExAc
Now we retrieve frequency data from the Broad ExAc [http://exac.broadinstitute.org]:
End of explanation
"""


merged = pd.merge(ESP_data,ExAc_data,how='outer',on=['Chromosome','Start','Stop','Ref','Alt','Disease',
                                                     'Assertion','Accession'])
print "There are %s variant-disease assertions with frequency data in either ExAC or ESP." %merged.shape[0]

num_SNPs = merged.shape[0]
print "There are %s SNP variant-disease assertions with frequency data in either ExAc or ESP." %num_SNPs

num_variants = len(set(zip(merged.Chromosome,merged.Start,merged.Stop,merged.Ref,merged.Alt)))
print "There are %s distinct variants in the set of variant-disease assertions with frequency data in either ExAc or ESP." %num_variants

num_diseases = len(set(zip(merged.Disease)))
print "There are %s distinct diseases in the set of variant-disease assertions with frequency data in either ExAc or ESP." %num_diseases

merged.to_csv(path_or_buf='merged_output.csv',sep=',',index=False,
                   columns=['Chromosome','Start','Stop','Ref','Alt','Assertion',
                            'Disease','Accession','ExAc_Overall_Frequency','ESP_Overall_Frequency'])

"""
Explanation: 2a.3. Merging
Let's merge the two data structures:
End of explanation
"""


from bs4 import BeautifulSoup
import urllib2

def getPubMedIDs(RCV):
    try:
        response = urllib2.urlopen("http://www.ncbi.nlm.nih.gov/clinvar/" + RCV)
    except urllib2.HTTPError:
        print "Page not found"
    else:
        soup = BeautifulSoup(response.read())
        tag = soup.body.find("div", id="clinvar_rec_pubmed_ids1")
        if tag:
            return tag.text.split(", ")
    return []

"""
Explanation: 2b. Pulling Pubmed Publications
The following function, using the Python module BeautifulSoup, retrieves PubMed IDs directly from the ClinVar website:
End of explanation
"""


print "Expecting 5 PubMedIDs for RCV000019428:", getPubMedIDs("RCV000019428")
print "Expecting 2 PubMedIDs for RCV000019429:", getPubMedIDs("RCV000019429")
print "Expecting 0 PubMedIDs for RCV000116253.2:", getPubMedIDs("RCV000116253.2")

"""
Explanation: Let's test this function on a few examples:
End of explanation
"""


import os
os.system('mkdir plots')

a = flat_clinvar.groupby(['Chromosome','Start','Stop','Ref','Alt'])

var_counts = list(a.size().values)
assertions_per_variant = [var_counts.count(i) for i in range(1,20)]

%matplotlib inline

import matplotlib
import matplotlib.pyplot as plt

plt.figure(figsize = (7,7))
plt.bar(range(1,len(assertions_per_variant)+1),assertions_per_variant,
        align='center',color='teal')
plt.xlabel('Assertions per Variant',fontsize=14)
plt.ylabel('Count',fontsize=14)
#plt.title('Pathogenicity Assertions per Variant in ClinVar')
plt.xlim([0,15])
plt.xticks(range(len(assertions_per_variant)),fontsize=12)
plt.savefig('plots/assertions_per_variant.pdf')

assertions_per_variant[0]/(len(a.groups))
len(a.groups)

"""
Explanation: 3. Penetrance
We are interested in the quantitative disease risk associated with pathogenic variation. We define the <i>penetrance</i> of genotype $G=j$ for disease $D$ as:
\begin{eqnarray}
\text{Penetrance} = P(D|G=j) = \frac{k \times P(G=j|D)}{P(G=j)} \
= \frac{k \times P(G=j|D)}{P(G=j|D)\times k + P(G=j|\overline{D})\times (1-k)}
\end{eqnarray}
where $k$ is the prevalence of the disease, $D$ signifies affected and $\overline{D}$ signifies unaffected.
4. Summarizing the Data
4.1 Diseasome-wide analyses with ClinVar variants
First, let's create a directory to store plots:
End of explanation
"""


# parameters above as Python variables
import numpy as np

prev = 1.0/500
het_range = np.linspace((0.001),(0.1),10)
pen = {} # to be computed below
genetic_model = 'AD'

"""
Explanation: 4.2 Uncertainty in the disease risk conveyed by pathogenic variation
Here we perform a sensitivity analysis for penetrance using the disease hypertrophic cardiomyopathy (HCM). Making our assumptions explicit:
<ul>
<li>The prevalence of HCM is $1/500$ [Maron et al., <i>Circulation</i> 1995]</li>
<li>There are more than a thousand causal variants for HCM [Maron et al., <i>J Am. Coll. Cardiol.</i>, 2012]. A conservative estimate for $h$, the heterogeneity parameter, is therefore $h=0.1$. A more realistic estimate is $h=0.001$</li>
<li>HCM is autosomal dominant: a single copy of a highly-penetrant causal allele is sufficient to cause disease</li>
</ul>
End of explanation
"""


HCM = merged[merged.Disease.str.contains('hypertrophic_cardiomyopathy')]
HCM = HCM.drop_duplicates(subset = ['Chromosome','Start','Stop','Ref','Alt']) # 81 distinct variants
raw_names = zip(HCM.Chromosome,HCM.Start,HCM.Ref,HCM.Alt)
clean_names = [str(chrom)+':'+str(int(pos))+str(ref)+'>'+alt for (chrom,pos,ref,alt) in raw_names]
HCM[['Chromosome','Start','Ref','Alt','ExAc_Overall_Frequency','ESP_Overall_Frequency','Accession','Disease']].head(5)

import numpy as np

max_freqs = [np.nanmax(np.array([x,y])) for (x,y) in zip(HCM.ESP_Overall_Frequency,HCM.ExAc_Overall_Frequency)]
max_freqs = [x**2 + 2*x*(1-x) for x in max_freqs]

for i,var in enumerate(clean_names):
    these_pens = []
    for h in het_range:
        f = max_freqs[i]
        these_pens += [min(prev*h/f,1)]
    pen[var] = these_pens

# now export to R
import csv
outfile = open('sens_analysis','w')
for k,v in pen.items():
    outfile.write(k)
    for i in v:
        outfile.write(','+str(i))
    outfile.write('\n')
outfile.close()

"""
Explanation: Recall:
\begin{eqnarray}
\text{Penetrance} = P(D|G=j) = \frac{k \times P(G=j|D)}{P(G=j)} \
= \frac{k \times P(G=j|D)}{P(G=j|D)\times k + P(G=j|\overline{D})\times (1-k)}
\end{eqnarray}
Thus, grabbing all HCM variants and computing penetrance can be accomplished by:
End of explanation
"""


HCM_PMIDs = {}
for acc in set(HCM.Accession):
    HCM_PMIDs[acc] = getPubMedIDs(acc)
HCM_PMIDs

"""
Explanation: Let's retrieve all PubMed IDs for these accessions:
End of explanation
"""


