<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 T cell epitopes of SARS-CoV2
Methods
Predict MHC-I binders for sars-cov2 reference sequences (S and N important)
Align with sars-cov and get conserved epitopes.
Best alleles to use?
Multiple sequence alignment of each protein to reference
find conservation of binders with closest peptide in each HCov sequence and determine identity
References
J. Mateus et al., “Selective and cross-reactive SARS-CoV-2 T cell epitopes in unexposed humans,” Science (80-. )., vol. 3871, no. August, p. eabd3871, Aug. 2020.
S. F. Ahmed, A. A. Quadeer, and M. R. McKay, “Preliminary Identification of Potential Vaccine Targets for the COVID-19 Coronavirus (SARS-CoV-2) Based on SARS-CoV Immunological Studies.,” Viruses, vol. 12, no. 3, 2020.
A. Grifoni et al., “A sequence homology and bioinformatic approach can predict candidate targets for immune responses to SARS-CoV-2,” Cell Host Microbe, pp. 1–10, 2020.
V. Baruah and S. Bose, “Immunoinformatics-aided identification of T cell and B cell epitopes in the surface glycoprotein of 2019-nCoV,” J. Med. Virol., no. February, pp. 495–500, 2020.
Common coronoviruses
https
Step1:  load ref genomes
Step2:  find orthologs in each genome
blast the genomes to find corresponding protein as names are ambigious
Step3:  predict MHC-I and MHC-II epitopes
Step4:  conservation
Step5:  Find conserved predicted epitopes in all proteins
Step6:  Compare predictions to mateus exp results
Step7:  check epitope selection method
Promiscuous binders are those high scoring above some threshold in multiple alleles. There are several ways to select them that can give different results. By default epitopepredict selects those in each allele above a percentile score cutoff and then counts how many alleles each peptide is present in. We can also limit our set in each protein across a genome to prevent large proteins dominating the list. We can also select by score and protein rank. The overlap is shown in the venn diagram.
<ASSISTANT_TASK:>
Python Code:
import os, math, time, pickle, subprocess
from importlib import reload
from collections import OrderedDict, defaultdict
import numpy as np
import pandas as pd
pd.set_option('display.width', 180)
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
import pathogenie
from Bio import SeqIO,AlignIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
Explanation: T cell epitopes of SARS-CoV2
Methods
Predict MHC-I binders for sars-cov2 reference sequences (S and N important)
Align with sars-cov and get conserved epitopes.
Best alleles to use?
Multiple sequence alignment of each protein to reference
find conservation of binders with closest peptide in each HCov sequence and determine identity
References
J. Mateus et al., “Selective and cross-reactive SARS-CoV-2 T cell epitopes in unexposed humans,” Science (80-. )., vol. 3871, no. August, p. eabd3871, Aug. 2020.
S. F. Ahmed, A. A. Quadeer, and M. R. McKay, “Preliminary Identification of Potential Vaccine Targets for the COVID-19 Coronavirus (SARS-CoV-2) Based on SARS-CoV Immunological Studies.,” Viruses, vol. 12, no. 3, 2020.
A. Grifoni et al., “A sequence homology and bioinformatic approach can predict candidate targets for immune responses to SARS-CoV-2,” Cell Host Microbe, pp. 1–10, 2020.
V. Baruah and S. Bose, “Immunoinformatics-aided identification of T cell and B cell epitopes in the surface glycoprotein of 2019-nCoV,” J. Med. Virol., no. February, pp. 495–500, 2020.
Common coronoviruses
https://www.cdc.gov/coronavirus/types.html
End of explanation
labels = {'sars':'NC_004718.3','scov2':'NC_045512.2','229E':'NC_002645.1','NL63':'NC_005831.2','OC43':'NC_006213.1','HKU1':'NC_006577.2'}
genomes = []
for l in labels:
    df = ep.genbank_to_dataframe(labels[l]+'.gb',cds=True)
    df['label'] = l
    genomes.append(df)
genomes = pd.concat(genomes)
scov2_df = genomes[genomes.label=='scov2']
scov2_df = scov2_df.drop_duplicates('locus_tag')
#print (genomes[['label','gene','product','length']])
def get_seqs(gene):
    seqs = []
    sub = genomes[genomes['gene']==gene]
    for i,r in sub.iterrows():
        s=SeqRecord(Seq(r.translation),id=r.label)
        seqs.append(s)
    return seqs
seqs=get_seqs('S')
Explanation: load ref genomes
End of explanation
pathogenie.tools.dataframe_to_fasta(genomes, idkey='locus_tag', descrkey='product', outfile='proteins.fa')
pathogenie.tools.make_blast_database('proteins.fa', dbtype='prot')
def get_orthologs(gene):
    sub = scov2_df[scov2_df['gene']==gene].iloc[0]    
    rec = SeqRecord(Seq(sub.translation),id=sub.gene)
    bl = pathogenie.tools.blast_sequences('proteins.fa', [rec], maxseqs=10, evalue=1e-4,
                              cmd='blastp', threads=4)
    bl = bl.drop_duplicates('sseqid')
    #print (bl.sseqid)
    found = genomes[genomes.locus_tag.isin(bl.sseqid)].drop_duplicates('locus_tag')
    #print (found)
    recs = pathogenie.tools.dataframe_to_seqrecords(found,
                            seqkey='translation',idkey='label',desckey='product')
    return recs
seqs = get_orthologs('S')
aln = pathogenie.clustal_alignment(seqs=seqs)
print (aln)
spikesars = SeqIO.to_dict(seqs)['sars'].seq
spikesars
sc2 = ep.genbank_to_dataframe('NC_045512.2.gb',cds=True)
sc2 = sc2.drop_duplicates('gene')
Explanation: find orthologs in each genome
blast the genomes to find corresponding protein as names are ambigious
End of explanation
m1_alleles = ep.get_preset_alleles('broad_coverage_mhc1')
m2_alleles = ep.get_preset_alleles('mhc2_supertypes')
P1 = base.get_predictor('netmhcpan') 
P1.predict_sequences(sc2, alleles=m1_alleles,threads=10,path='scov2_netmhcpan',length=9,overwrite=False)
P1.load(path='scov2_netmhcpan')
pb1 = P1.promiscuous_binders(n=3, cutoff=.95)
P2 = base.get_predictor('netmhciipan') 
P2.predict_sequences(sc2, alleles=m2_alleles,threads=10,path='scov2_netmhciipan',length=15,overwrite=False,verbose=True)
P3 = base.get_predictor('tepitope') 
P3.predict_sequences(sc2, alleles=m2_alleles,threads=10,path='scov2_tepitope',length=15,overwrite=False)
P3.load(path='scov2_tepitope')
P2.load(path='scov2_netmhciipan')
pb2 = P2.promiscuous_binders(n=3, cutoff=.95, limit=70)
#rb2 = P2.promiscuous_binders(n=3, cutoff_method='rank', cutoff=40)
pb2.name.value_counts()
#pb3 = P3.promiscuous_binders(n=3, cutoff=50, cutoff_method='rank', limit=50)
#pb3.name.value_counts()
Explanation: predict MHC-I and MHC-II epitopes
End of explanation
import difflib
def get_conservation(x, w):    
    m = difflib.get_close_matches(x, w, n=1, cutoff=.67)
    if len(m)==0:
        return 0
    else:
        m=m[0]
        s = difflib.SequenceMatcher(None, x, m)
        return s.ratio()
def find_epitopes_conserved(pb,gene,locus_tag):
    
    seqs = get_orthologs(gene)
    df = pb[pb.name==locus_tag]
    #print (df)
    print (len(seqs),len(df))
    s=seqs[0]
    for s in seqs:
        if s.id == 'scov2': 
            continue
        w,ss = peptutils.create_fragments(seq=str(s.seq), length=11)
        df.loc[:,s.id] = df.peptide.apply(lambda x: get_conservation(x, w),1)
    df.loc[:,'total'] = df[df.columns[8:]].sum(1)
    df = df.sort_values('total',ascending=False)
    df = df[df.total>0]
    df = df.round(2)
    return df
df = find_epitopes_conserved(pb2, 'S','GU280_gp02')
#df.to_csv('S_netmhciipan_conserved.csv')
Explanation: conservation: find identity to closest peptide in each HCoV sequence
End of explanation
res=[]
for i,r in scov2_df.iterrows():
    print (r.locus_tag,r.gene)    
    df = find_epitopes_conserved(pb2,r.gene,r.locus_tag)
    df['gene'] = r.gene
    res.append(df)
res = pd.concat(res).sort_values('total',ascending=False).dropna().reset_index()
print (len(res),len(pb2))
res.to_csv('scov2_netmhciipan_conserved.csv')
cols = ['gene','peptide','pos','alleles','sars','229E','NL63','OC43','HKU1']
h=res[:30][cols].style.background_gradient(cmap="ocean_r",subset=['sars','229E','NL63','OC43','HKU1']).set_precision(2)
#res[:30][cols]
Explanation: Find conserved predicted epitopes in all proteins
End of explanation
s1 = pd.read_csv('mateus_hcov_reactive.csv')
hits=[]
w = list(res.peptide)
for i,r in s1.iterrows():    
    m = difflib.get_close_matches(r.Sequence, w, n=2, cutoff=.6)
    #print (r.Sequence,m,r.Protein)
    if len(m)>0:
        hits.append(m)
    else:
        hits.append(None)
        
s1['hit'] = hits
display(s1)
print (len(s1.hit.dropna())/len(s1))
Explanation: Compare predictions to mateus exp results
End of explanation
reload(base)
P = base.get_predictor('tepitope') 
P.predict_sequences(sc2, alleles=m2_alleles[:4],names=['GU280_gp01','GU280_gp02','GU280_gp03','GU280_gp04'],threads=10,length=9)
pb= P.promiscuous_binders(n=2, cutoff=.98, limit=20)
pb.name.value_counts()
rb= P.promiscuous_binders(n=3, cutoff_method='rank',cutoff=30,limit=20)
rb.name.value_counts()
sb= P.promiscuous_binders(n=2, cutoff_method='score',cutoff=3.5,limit=20)
sb.name.value_counts()
from matplotlib_venn import venn3
ax = venn3((set(pb.peptide),set(rb.peptide),set(sb.peptide)), set_labels = ('default', 'ranked', 'score'))
b=P.get_binders(cutoff=10, cutoff_method='rank')
func = max
s=b.groupby(['peptide','pos','name']).agg({'allele': pd.Series.count,
             'core': base.first, P.scorekey:[func,np.mean],
             'rank': np.median})
s.columns = s.columns.get_level_values(1)
s.rename(columns={'max': P.scorekey, 'count': 'alleles','median':'median_rank',
                 'first':'core'}, inplace=True)
s = s.reset_index()
s
s.name.value_counts()
s=s.groupby('name').head(10)
s.name.value_counts()
Explanation: check epitope selection method
Promiscuous binders are those high scoring above some threshold in multiple alleles. There are several ways to select them that can give different results. By default epitopepredict selects those in each allele above a percentile score cutoff and then counts how many alleles each peptide is present in. We can also limit our set in each protein across a genome to prevent large proteins dominating the list. We can also select by score and protein rank. The overlap is shown in the venn diagram.
End of explanation
<END_TASK>