import os, math, time, pickle, subprocess
from importlib import reload
from collections import OrderedDict
import numpy as np
import pandas as pd
pd.set_option('display.width', 100)
import epitopepredict as ep
from epitopepredict import base, sequtils, plotting, peptutils, analysis
from IPython.display import display, HTML, Image
%matplotlib inline
import matplotlib as mpl
import pylab as plt

"""
Explanation: Predicting cattle T cell Epitopes in M.bovis
This is a simplified version of an analysis we used to predict cattle T cell epitopes in the proteome of M.bovis. We apply three binding prediction methods to the M. bovis proteome using a subset of human MHC-II alleles to approximate the bovine immune response. Two different strategies are then applied to filter the resulting set of binders: global binder ranking and epitope density. Both these sets were used to make 20mer peptides for synthesis. 
Several other metrics are applied to produce a final list of peptide candidates.
References

Integrated computational prediction and experimental validation identifies promiscuous T cell epitopes in the proteome of pathogenic mycobacteria. Damien Farrell, Gareth Jones, Chris Pirson, Kerri M Malone, Kevin Rue-Albrecht, Anthony J Chubb, H Martin Vordermeier, Stephen V. Gordon. 2016
M. Pandya, “Definition of Bovine Leukocyte Antigen Diversity and Peptide Binding Profiles for Epitope Discovery,” 2016.
O. T. Schubert et al., “Absolute Proteome Composition and Dynamics during Dormancy and Resuscitation of Mycobacterium tuberculosis,” Cell Host Microbe, 2015.
End of explanation
"""


mhc2alleles = ep.get_preset_alleles('bovine_like_mhc2')
mhc1alleles = ['BoLA-1:01901', 'BoLA-2:00801','BoLA-2:01201','BoLA-4:02401','BoLA-2:07001',
               'BoLA-3:01701','BoLA-1:02301','BoLA-2:01801','BoLA-6:01302','BoLA-3:00201']

"""
Explanation: Alleles

8 bovine similar HLA MHC-II alleles were used to approximate cattle responses
BoLA class I alleles used here are those most frequent in 81 Holstein cattle from the UVM Dairy Center of Excellence research herd assayed by PCR-SSP. (Different alleles were used in the original study).
End of explanation
"""


#srm data
srm = pd.read_csv('srm_mtb.csv')
#srm = srm[srm.length<=400]
proteome = sequtils.genbank_to_dataframe('MTB-H37Rv.gb',cds=True)
#filter proteins on srm data
proteome = proteome.merge(srm[['locus_tag','concentration']],on='locus_tag',how='inner')
proteome = proteome[proteome.concentration>0]
print (len(proteome))
proteome.head()

"""
Explanation: Filtering with MTB srm data
Here we load the protein sequences from a genbank file. MTB is used since it is 99% similar to Bovis. You can find this file at https://www.ncbi.nlm.nih.gov/genome/166?genome_assembly_id=159857
The proteins are filtered using absolute abundance levels of proteins in an unfractionated mixed lysate of M. tuberculosis H37Rv cultures identified by Schubert et al. All proteins undetected using selected reaction monitoring (SRM) (concentration = 0 and <2 peptides detected) in this study were filtered out. This leaves 1870 out of ~4000 proteins.
End of explanation
"""


P1 = base.get_predictor('tepitope')
b = P1.predict_sequences(proteome, alleles=mhc2alleles, threads=8, path='mtb_tepitope')

P2 = base.get_predictor('netmhciipan')
b = P2.predict_sequences(proteome, alleles=mhc2alleles, threads=8, path='mtb_netmhciipan',overwrite=False)

reload(base)
P3 = base.get_predictor('netmhcpan')
b = P3.predict_sequences(proteome[:200], alleles=mhc1alleles, threads=8, path='mtb_netmhcpan',overwrite=False)

"""
Explanation: Predict binders
Run binding predictions for each of the three and save to disk. This enables us to reload the results to re-analyse.
End of explanation
"""


reload(analysis)
predictors = ['netmhcpan','tepitope','netmhciipan']
clusters = {}
for name in predictors:
    P = base.get_predictor(name)
    path = 'mtb_'+name
    P.load(path)   
    pb=P.promiscuous_binders(n=3)
    #print (pb[:20])
    cl = analysis.find_clusters(pb, min_binders=3)
    clusters[name] = cl
 

"""
Explanation: Compute clusters of promiscuous binders
Here we have detected clusters of MHC II binders. These clusters vary in length and for peptides synthesis we must create a list of 20-mers covering shorter sequences or splitting longer sequences into 2. For the 20mers hydrophobicity and net charge are calculated.
End of explanation
"""


c1=clusters['netmhciipan']
c2=clusters['tepitope']
x = analysis.get_overlaps(c1,c2,how='inside')
print (x)

final = ep.analysis.create_nmers(cl, proteome, key='20mer', length=20, margin=2)     
final = analysis.peptide_properties(final, colname='20mer')
print (final.head())   

reload(plotting)
#print (P.get_names())
ax = plotting.plot_tracks([P],name='Bla-g-5',cutoff=.95, cutoff_method='default',n=2, legend=True, figsize=(12,5),regions=cl)


"""
Explanation: overlapping clusters?
End of explanation
"""


