<jupyter_start><jupyter_text>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* <NAME> 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.* <NAME>, <NAME>, and <NAME>, “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.* <NAME> 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.* <NAME> and <NAME>, “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. Epitope Loss in Mutations* https://www.biorxiv.org/content/10.1101/2020.03.27.012013* https://www.biorxiv.org/content/10.1101/2020.04.10.029454v1?ct=* https://www.biorxiv.org/content/10.1101/2020.04.07.030924v1 Common coronoviruses* https://www.cdc.gov/coronavirus/types.html<jupyter_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 pybioviz
from bokeh.io import show, output_notebook
output_notebook()
import pathogenie
from Bio import SeqIO,AlignIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord<jupyter_output><empty_output><jupyter_text>load ref genomes<jupyter_code>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')<jupyter_output><empty_output><jupyter_text>find orthologs in each genome blast the genomes to find corresponding protein as names are ambigious<jupyter_code>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)
import Levenshtein
for a in aln:    
    r=Levenshtein.ratio(str(aln[1].seq).replace('-',''), str(a.seq).replace('-',''))
    print (round(r,2), a.id)
spikesars = SeqIO.to_dict(seqs)['sars'].seq
spikesars
p = pybioviz.plotters.plot_sequence_alignment(aln, annot = {'polybasic cleavage site':690,'RBD contact residues':480})
#output_file('alignment.html')
show(p)
sc2 = ep.genbank_to_dataframe('NC_045512.2.gb',cds=True)
sc2 = sc2.drop_duplicates('gene')<jupyter_output><empty_output><jupyter_text>predict MHC-I and MHC-II epitopes<jupyter_code>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,cpus=10,path='netmhcpan',length=9,overwrite=False,verbose=True)
P1.load(path='netmhcpan')
pb1 = P1.promiscuous_binders(n=3, cutoff=.95)
P2 = base.get_predictor('netmhciipan') 
P2.predict_sequences(sc2, alleles=m2_alleles,cpus=10,path='netmhciipan',length=15,overwrite=False,verbose=True)
P3 = base.get_predictor('tepitope') 
P3.predict_sequences(sc2, alleles=m2_alleles,cpus=10,path='tepitope',length=15,overwrite=False)
P3.load(path='tepitope')
P2.load(path='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()<jupyter_output><empty_output><jupyter_text>conservation: find identity to closest peptide in each HCoV sequence<jupyter_code>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')<jupyter_output><empty_output><jupyter_text>Find conserved predicted epitopes in all proteins<jupyter_code>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]<jupyter_output><empty_output><jupyter_text>Compare predictions to mateus exp results<jupyter_code>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))<jupyter_output><empty_output><jupyter_text>check epitope selection methodPromiscuous 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.<jupyter_code>reload(base)
P = base.get_predictor('tepitope') 
P.predict_sequences(sc2, alleles=m2_alleles[:4],names=['GU280_gp01','GU280_gp02','GU280_gp03','GU280_gp04'],cpus=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()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14664291281594477, 'fraction_numerical': 0.03941654604164347, 'mean_word_length': 5.527616279069767, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 42, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '50401', 'license': 'MIT', 'n_tokens_mistral': 3650, 'n_tokens_neox': 3285, 'n_words': 654, 'path': 'sarscov2/epitopes.ipynb', 'repo_name': 'vishalbelsare/teaching', 'size': '9009'}