<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 epitopepredict API basic usage
Methodology
Predictors for each method inherit from the Predictor class and all implement a predict method for scoring a single sequence. This may wrap methods from other modules and/or call command line predictors. For example the TepitopePredictor uses the mhcpredict.tepitope module. This method should return a Pandas DataFrame. The predictProteins method is used for multiple proteins contained in a dataframe of sequences in a standard format. This is created from a genbank or fasta file (see examples below). For large numbers of sequences predictProteins you should provide a path so that the results are saved as each protein is completed to avoid memory issues, since many alleles might be called for each protein. Results are saved with one file per protein in csv format. Results can be loaded into the predictor individually or all together using the load method.
Step1:  train models for basicmhc1 predictor before first use
this will be done automatically when you call the prediction methods if not run first
Step2:  predict random peptides
Step3:  results are returned as dataframes, sorted by score
Step4:  predict from a sequence
Step5:  predict n-mers from multiple protein sequence
This example loads protein sequences from mycobacterium tuberculosis from a genbank file
Step6:  get promiscuous binders - those peptides above the cutoff for at least n alleles
Step7:  run other predictors
Step8:  plot binders in a protein sequence
<ASSISTANT_TASK:>
Python Code:
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, tepitope, plotting, utilities, peptutils, mhclearn
from IPython.display import display, HTML, Image
%matplotlib inline
import matplotlib as mpl
import pylab as plt
#get preset alleles
m2_alleles = ep.get_preset_alleles('mhc2_supertypes')
m1_alleles = ep.get_preset_alleles('mhc1_supertypes')
print (m1_alleles)
print (m2_alleles)
data = mhclearn.get_training_set()
data.length.value_counts()
e = mhclearn.get_evaluation_set1()
print (e.length.value_counts())
Explanation: epitopepredict API basic usage
Methodology
Predictors for each method inherit from the Predictor class and all implement a predict method for scoring a single sequence. This may wrap methods from other modules and/or call command line predictors. For example the TepitopePredictor uses the mhcpredict.tepitope module. This method should return a Pandas DataFrame. The predictProteins method is used for multiple proteins contained in a dataframe of sequences in a standard format. This is created from a genbank or fasta file (see examples below). For large numbers of sequences predictProteins you should provide a path so that the results are saved as each protein is completed to avoid memory issues, since many alleles might be called for each protein. Results are saved with one file per protein in csv format. Results can be loaded into the predictor individually or all together using the load method.
End of explanation
reload(mhclearn)
mhclearn.train_models(overwrite=True, alleles=['H-2-Kb'])#,'HLA-A*01:01','HLA-A*02:01','HLA-A*02:02','HLA-A*02:03','HLA-A*02:06'])
a = mhclearn.get_allele_names()
print (len(a))
Explanation: train models for basicmhc1 predictor before first use
this will be done automatically when you call the prediction methods if not run first
End of explanation
reload(base)
reload(mhclearn)
seqs = peptutils.create_random_sequences(1000)
print (seqs[:5])
df = pd.DataFrame(seqs,columns=['peptide'])
P = base.get_predictor('basicmhc1')
b = P.predict_peptides(df.peptide, alleles=m1_alleles, threads=8)
Explanation: predict random peptides
End of explanation
print (b[:5])
Explanation: results are returned as dataframes, sorted by score
End of explanation
seq = 'MTDDPGSGFTTVWNAVVSELNGDPKVDDGP'
b = P.predict_sequences(seq, alleles=m1_alleles, length=9)
b
Explanation: predict from a sequence
End of explanation
#load protein sequences into a dataframe
prots = ep.genbank_to_dataframe(base.mtb_genome, cds=True)
prots[:5]
proteins = ['Rv3615c','Rv3875','Rv1886c']
P1 = base.get_predictor('basicmhc1')
binders = P1.predict_sequences(prots, names=proteins, alleles=m1_alleles, length=9, threads=8)
binders.head()
Explanation: predict n-mers from multiple protein sequence
This example loads protein sequences from mycobacterium tuberculosis from a genbank file
End of explanation
pb1 = P1.promiscuous_binders(n=2, cutoff=5, cutoff_method='rank')
pb1[:5]
Explanation: get promiscuous binders - those peptides above the cutoff for at least n alleles
End of explanation
reload(base)
P2 = base.get_predictor('netmhcpan')
P2.predict_sequences(prots, names=proteins, alleles=m1_alleles, length=9, threads=8)
reload(base)
P3 = base.get_predictor('mhcflurry')
P3.predict_sequences(prots, names=proteins, alleles=m1_alleles, length=9, verbose=True)
prot='Rv3615c'
b1 = P1.get_binders(cutoff=5,name=prot,cutoff_method='rank')
b2 = P2.get_binders(cutoff=5,name=prot,cutoff_method='rank')
b3 = P3.get_binders(cutoff=5,name=prot,cutoff_method='rank')
from matplotlib_venn import venn3
ax = venn3((set(b1.peptide),set(b2.peptide),set(b3.peptide)), set_labels = ('basicmhc1','netmhcpan','mhcflurry'))
plt.title('top 5 ranked binders in compared, Rv00010c')
plt.savefig('basicmhc1_binders_compared.jpg',dpi=150)
pb1 = P1.promiscuous_binders(n=2,cutoff=10,limit=10,cutoff_method='rank')
pb2 = P2.promiscuous_binders(n=2,cutoff=10,limit=10,cutoff_method='rank')
pb3 = P3.promiscuous_binders(n=2,cutoff=10,limit=10,cutoff_method='rank')
ax = venn3((set(pb1.peptide),set(pb2.peptide),set(pb3.peptide)), set_labels = ('basicmhc1','netmhcpan','mhcflurry'))
plt.title('promiscuous binders compared, 3 methods')
plt.savefig('basicmhc1_prom_binders_compared.jpg',dpi=150)
Explanation: run other predictors
End of explanation
reload(plotting)
#get names of proteins stored in results of predictor
print (P2.get_names())
ax = plotting.plot_tracks([P1,P2,P3],name='Rv3875',cutoff=5,cutoff_method='rank',n=2,legend=True)
plt.tight_layout()
plt.savefig('comparison_rv3875.png',dpi=150)
#ax = plotting.plot_binder_map(P1,name='Rv0011c',cutoff=10)
#plt.savefig('mhc_rv0011c_map.png',dpi=150)
reload(plotting)
from bokeh.io import show, output_notebook
output_notebook()
p = plotting.bokeh_plot_tracks([P1,P2,P3],name='Rv0008c',cutoff=10,cutoff_method='rank',n=2,width=800)
show(p)
Explanation: plot binders in a protein sequence
End of explanation
<END_TASK>