生物信息學方法預測病原體抗原蛋白序列中多肽疫苗候選表位

趙靜靜,澹小秀,王廣志,歐陽儉,簡星星,3*,謝 鷺*

（1.上海海洋大學食品學院,中國上海 201306;2.上海市生物醫藥技術研究院上海生物信息技術研究中心,中國上海 201203;3.中南大學湘雅醫院生物信息中心,中國湖南 長沙 410008）

疫苗是當前人類對抗病原生物最經濟有效的方式?，F階段,針對病原體蛋白中能夠引起免疫系統效應功能的表位而設計的多肽疫苗是新型疫苗研發的熱點。多肽疫苗一般是指利用化學合成技術,按照病原體抗原基因中具有免疫原性的一段或多段氨基酸序列構成而制備的疫苗[1]。其因成分簡單、特異性強、安全性高、易于保藏等優點被廣泛應用于傳染病和腫瘤的預防和治療中,例如:預防人乳頭瘤狀病毒的疫苗[2]、預防乙型肝炎病毒的HBsAg/preS疫苗[3]、新型艾滋病疫苗RV144[4]及個性化腫瘤新抗原疫苗[5]等。

近年來,得益于生物信息學技術的迅猛發展,反向疫苗學在疫苗設計領域顯示出廣闊的應用空間。以“序列-結構-功能”思想為依據,借助生物信息學工具對病原體高通量的組學信息（基因組學、轉錄組學、蛋白質組學等）預先進行嚴格的預測和篩選,能極大地提高研究人員發現病原體蛋白質序列中候選表位的效率。隨著候選表位的持續產出和收集,相關表位數據庫應運而生,如IEDB（Immune Epitope Database）[6]及本課題組自行構建的dbPepNeo[7]等。研究人員利用這些數據集開發了諸多表位預測和篩選的工具,但其中的多數均只針對特定的功能進行開發,如:NetMHC-pan EL 4.0[8]主要用于預測與主要組織相容性復合體（major histocompatibility complex,MHC）結合的抗原表位;本課題組開發的INeo-Epp[9]則用于預測具有免疫原性的抗原表位。因此,表位預測工具功能的多樣化,對于候選表位的鑒定以及多肽疫苗的設計具有重要指導意義[10～12]。

本研究基于反向疫苗學策略,綜合利用多種生物信息學工具,構建了一個完整的、易于操作的表位預測和篩選流程,并以當下正在全球蔓延的嚴重急性呼吸綜合征冠狀病毒2（severe acute respiratory syndrome coronavirus 2,SARS-CoV-2）為實例,成功地檢驗了該流程的實用性。文中共篩選到34條關于SARS-CoV-2的T細胞候選表位,其中20條候選表位與數據庫中經過驗證的表位高度同源,且能夠被T細胞受體（T cell receptor,TCR）特異性識別。此外,為了驗證流程的適用性,我們還對鼠類肉瘤病毒癌基因KRAS與大腸桿菌基因OmpC所編碼蛋白質的抗原表位進行了示例性分析。本研究不僅為新冠疫苗設計提供了候選表位,而且還建立了一個廣泛適用于多肽疫苗候選表位預測和篩選的操作流程。

1 材料和方法

1.1 數據下載

本研究涉及的SARS-CoV-2的結構蛋白包括刺突糖蛋白（spike glycoprotein,S）、膜糖蛋白（membrane glycoprotein,M）、包膜蛋白（envelope protein,E）與核衣殼蛋白（nucleocapsid phosphoprotein,N）,它們的蛋白質序列與鼠類肉瘤病毒癌基因KRAS及大腸桿菌基因OmpC相關的蛋白質序列文件均從NCBI（https://www.ncbi.nlm.nih.gov/）數據庫獲取,檢索ID分別為QHR63290.2、QHD43419.1、QHD-43418.1、QHD43423.2、CAC5395073.1 和 ADU34-074.2。

1.2 表位預測識別

本研究使用生物信息學工具NetMHCpan EL 4.0[8]（http://tools.immuneepitope.org/mhci/）,對SARSCoV-2的S蛋白、M蛋白、E蛋白與N蛋白及KRAS和OmpC蛋白序列中潛在的T細胞表位分別進行預測。然后,依據“抗原肽-MHC”復合物間的親和力（通常,%rank<0.5被定義為強結合,%rank<2為弱結合,其他為不結合）進行篩選,保留%rank在0～2的短肽作為結合表位,進一步評估其作為多肽疫苗候選表位的潛力。

1.3 表位免疫原性與抗原性預測

為了篩選出具有免疫原性的候選表位,我們使用IEDB中提供的免疫原性預測工具[6,13]（http://tools.iedb.org/immunogenicity/）及本課題組自行開發的INeo-Epp[9]（http://www.biostatistics.online/INeo-Epp/）兩個生物信息學工具,對S蛋白、M蛋白、E蛋白與N蛋白及KRAS和OmpC蛋白結合表位的免疫原性分別進行預測。IEDB與INeo-Epp工具是通過不同的算法對“抗原肽-MHC”復合物錨定位置處氨基酸的理化特性進行建模,預測效果表現良好[前者的曲線下面積（area under the curve,AUC）=0.69,后者AUC=0.81]。對于免疫原性評分（score）,一般認為評分越高表示該表位引發免疫反應的可能性越大[14]。為了降低假陽性率,通常以0.2為標準進行篩選。因此,本研究僅保留IEDB中score大于0.2和INeo-Epp中輸出結果為強陽性（positive-high,PH）的表位。

隨后,我們使用Vaxijen v2.0[15]（http://www.ddgpharmfac.net/vaxijen/）工具預測結合表位的抗原性。Vaxijen是第一個利用非序列比對策略預測表位抗原性的生物信息學工具,也是反向疫苗學中常用的工具。該模型設置了0～1的閾值,以預測病毒、細菌和腫瘤蛋白質序列中的保護性抗原。當閾值為0.4時,該模型在病毒測試數據集的效果較好（AUC=0.74）。因此,本研究在0.4的閾值下獲取S蛋白、M蛋白、E蛋白、N蛋白、KRAS蛋白及OmpC蛋白結合表位的抗原性分值。

1.4 表位理化特性分析

為了降低多肽疫苗輸注后患者發生不良反應的概率,本研究使用AllerTOP v2.0[16]（https://www.ddg-pharmfac.net/AllerTOP/）、ClanTox[17]（http://www.clantox.cs.huji.ac.il/）工具對結合表位的理化特性進行分析。AllerTOP和ClanTox擁有非常簡潔的操作界面,不需要額外設置閾值即可得到預測結果。此外,AllerTOP還是當前較被認可的預測蛋白質致敏性的工具,它通過自協方差和交叉協方差變換的方法對氨基酸進行編碼,綜合利用邏輯回歸、決策樹、隨機森林、樸素貝葉斯、多層感知器、K-近鄰等算法進行建模,最后選擇結果最好的K-近鄰算法（AUC=0.85）完成工具的構建。而ClanTox是一種蛋白質毒素分類器,可以根據蛋白質的一級結構計算出該蛋白質是否為毒性蛋白質。在本研究中,我們僅保留AllerTOP和ClanTox輸出結果為非致敏（non-allergen）與非毒性（nontoxin）的表位作為SARS-CoV-2多肽疫苗研發的候選表位。

1.5 人群覆蓋率評估

根據先前的研究[18～19],我們構建了12種在中國人群中常見的人類白細胞抗原（human leukocyte antigen,HLA）的等位基因（allele）分型。它們分別為HLA-A*01:01、HLA-A*02:01、HLA-A*03:01、HLA-A*11:01、HLA-A*23:01、HLA-A*24:02、HLA-B*07:02、HLA-B*08:01、HLA-B*35:01、HLA-B*40:01、HLA-B*44:02、HLA-B*44:03（表1）。此外,我們使用IEDB提供的人口覆蓋率計算工具Population Coverage[20]（http://tools.immuneepitope.org/population/）評估了每條候選表位在中國及世界人群中的覆蓋率。需要注意的是,使用Population Coverage工具評估表位在特定人群中的覆蓋率時,需要提交一份包含該表位所對應HLA等位基因分布頻率的信息文件。

表1 12種HLA等位基因的分布頻率Table 1 Distribution frequency of 12 HLA alleles

1.6 同源比對及特異性TCR檢索

IEDB[21]（http://www.iedb.org/）、VDJdb[22]（https://vdjdb.cdr3.net/）、McAPS-TCR[23]（http://friedmanlab.weizmann.ac.il/McPAS-TCR/）是當前表位及相關TCR數據儲存最多的數據庫。首先,本研究收集和整理了3個數據庫中有關human、CD8+T、TCRβ的表位信息及其相關的TCR序列信息。然后,應用BLAST（Basic Local Alignment Search Tool）[24]（https://blast.ncbi.nlm.nih.gov/Blast.cgi/）工具,將本研究篩選的候選表位與數據庫中的表位進行同源比對,記錄比對結果中所有同源表位的序列信息及E value和identity值。最后,在已整理的TCR數據中手動檢索與該同源表位匹配的TCR序列信息。

2 結果

2.1 SARS-CoV-2結構蛋白中T細胞識別的候選表位

本研究綜合多個生物信息學工具（圖1A）,系統地篩選了SARS-CoV-2結構蛋白序列中可激發T細胞免疫應答的表位,即T細胞表位。首先,我們構建了12種在中國人群普遍存在的HLA等位基因分型:HLA-A*01:01、HLA-A*02:01、HLAA*03:01、HLA-A*11:01、HLA-A*23:01、HLA-A*24:02、HLA-B*07:02、HLA-B*08:01、HLAB*35:01、HLA-B*40:01、HLA-B*44:02、HLAB*44:03;然后,預測了SARS-CoV-2結構蛋白中長度為8～12個氨基酸的T細胞表位。通過“抗原肽-MHC”結合預測工具NetMHCpan EL,共得到117 180條預測表位。之后,根據%rank分值篩選出0～2的結合表位,僅有2 593條表位顯示能與MHC分子結合。

免疫原性是指抗原誘導機體發生特異性免疫應答的性能,抗原性則是指抗原與其所誘導產生的抗體或致敏淋巴細胞特異性結合的能力。簡而言之,表位的免疫原性與抗原性越強,表明其被T細胞捕獲并引發免疫應答的潛力越高?？紤]到表位的免疫原性和抗原性均是決定該抗原能否被T淋巴細胞識別的關鍵特征,我們進一步對結合表位的免疫原性及抗原性進行了分析。首先,使用IEDB和INeo-Epp兩個工具分別對2 593條SARS-CoV-2結合表位的免疫原性進行預測,保留在兩個工具中結果都較好的表位（IEDB,score>0.2;INeo-Epp,score>0.5即為PH）,共獲得131條;隨后,使用Vaxijen工具分析這131條表位的抗原性,最終僅獲得71條既有免疫原性又有抗原性的表位。

疫苗設計的前提是對該病毒候選表位有充分的理解。為此,我們使用蛋白質理化特性分析工具AllerTOP和ClanTox評估了這71條表位的致敏性與毒性,發現其中非致敏和非毒性的候選表位有49條。此外,由于HLA等位基因的多態性,我們發現候選表位中存在同一表位與多個HLA等位基因結合的情況,例如:表位M135～144和表位S1110～1121以不同的%rank 分值與 HLA-B*40:01、HLA-B*44:02及HLA-B44:03三種分型結合。候選表位與HLA等位基因分型的結合率如圖1B所示。

最后,在忽略同一表位與不同HLA等位基因結合的情況下,我們仍在49條表位中找到了34條獨特的表位（表2）,其中15條來自S蛋白序列、9條來自M蛋白序列、5條來自E蛋白序列、5條來自N蛋白序列（圖1C）?？傮w來看,這34條候選表位經過了系統和嚴格的篩選,可以與MHC分子結合,并由抗原提呈細胞加工后呈遞到細胞表面,供TCR識別,繼而誘導T細胞發揮效應功能。

表2 SARS-CoV-2的候選表位信息Table 2 Information about candidate epitopes of SARS-CoV-2

圖1 多肽疫苗候選表位預測流程及其在SARS-CoV-2中的應用（A）SARS-CoV-2 T細胞候選表位的預測流程圖;（B）候選表位與HLA等位基因的結合率;（C）不同預測過程中得到的SARSCoV-2表位數量。Fig.1 The workflow for predicting candidate epitopes for peptide vaccines and its application in SARS-CoV-2（A）Screening process of SARS-CoV-2 T cell candidate epitopes;（B）Binding rates of candidate epitopes and HLA alleles;（C）The number of SARS-CoV-2 epitopes retained during the screening process.

2.2 候選表位的人群覆蓋率

由于與表位結合的HLA等位基因在不同國家的分布頻率差異明顯,我們評估了SARS-CoV-2的34條候選表位在中國及世界范圍的覆蓋率。首先,我們獲取了HLA等位基因在中國及世界范圍內的分布頻率[18～19]。然后,基于HLA等位基因的分布頻率,使用Population Coverage工具計算了每條候選表位在中國和世界人群的覆蓋率。結果顯示,本研究篩選的SARS-CoV-2候選表位在中國人群廣泛覆蓋（1%～91%）（表2）。其中,表位S1110～1121與表位 M135～144在中國人群的適用范圍高達91%,同時在世界人群的適用范圍也達到71%。

2.3 同源表位及其誘導的TCR序列

為了檢驗本研究所構建預測流程的實用性,我們使用BLAST工具將本文篩選的34條SARSCoV-2表位與數據庫中經過驗證的表位進行同源比對。結果顯示,有20條SARS-CoV-2候選表位與 IEDB、VDJdb、McAPS-TCR 3個數據庫中經過驗證的表位同源（表3）。另外,當比對結果的參數E value調整為 0.05 時,E29～38和 E30～37這兩條表位與數據庫中LLAILTYYV（blast epitope）表位的identity高達100%。與此同時,我們還在數據庫中檢索到6條與該同源表位相關的特異性TCR序列（CASSLVRDRHTEAFF、CASSPTGTGGSDTQYF、CASSQAGEQYF、CASSWVGGADTQYF、CASTVRQGSNQPQHF 和 CSASFHNGFWGGTEAFF）。

表3 SARS-CoV-2的同源表位及相關TCR序列信息Table 3 Information about homologous epitopes and relevant TCR sequences of SARS-CoV-2

2.4 其他病原體蛋白質序列中候選表位的預測與篩選

為了進一步驗證該流程的適用性,我們在NCBI中檢索了鼠類肉瘤病毒癌基因KRAS及大腸桿菌基因OmpC相關的蛋白質序列,并對蛋白質序列中適合多肽疫苗研發的表位進行了預測和篩選。從表4可知,本研究共預測到5條有關KRAS的表位,即 KLVVVGAGGVGK、VVVGAGGVGK、LVVVGAGGVGK、QYMRTGEGFL和YMRTGEGF?；谶@5條表位的%rank分值、免疫原性、抗原性、致敏性及毒性,我們發現7～16區段的表位VVVGAGGVGK有潛力成為KRAS多肽疫苗研發的候選表位。另外,從表5的同源比對結果可以看出,VVVGAGGVGK與數據庫中YMDDVVLGA表位的identity達到80%,同時還有一條特異性TCR序列CASSYLTGEGDYGYTF與該同源表位相關,這進一步表明表位VVVGAGGVGK具有免疫原性,能夠誘導T細胞免疫應答。

表4 KRAS的候選表位信息Table 4 Information about candidate epitopes of KRAS

表5 KRAS的同源表位及相關TCR序列信息Table 5 Information about homologous epitopes and relevant TCR sequences of KRAS

我們在基因OmpC相關的蛋白質序列中共預測到12條表位（表6）。從預測表位的各預測分值及同源比對結果可知,表位YEGFGIGGAI（205～214）、VLPEFGGDTY（127～136）、TDVLPEFGGDTY（125～136）、DVLPEFGGDTY（126～136）可以作為OmpC多肽疫苗研發的候選表位（表6～7）。而且,對比分析 125～136、126～136、127～136 區段可以發現,125～136區段是OmpC蛋白中適用于多肽疫苗研發的優勢區段。由此,我們認為本研究所提出的預測流程可適用于病原體抗原蛋白序列中多肽疫苗候選表位的預測和篩選。

表6 OmpC的候選表位信息Table 6 Information about candidate epitopes of OmpC

表7 OmpC的同源表位及相關TCR序列信息Table 7 Information about homologous epitopes and relevant TCR sequences of OmpC

3 討論

快速、安全、靈活、高效且低成本的新型疫苗研發路線是未來人類對抗病原體最有效的手段[25]。近年來,隨著生物信息學技術的發展和應用,研究人員已經可以從病原體高通量組學信息中挖掘出大量的免疫原性表位序列,同時實現對T/B細胞靶點處表位的精準識別,從而為新型疫苗的研發提供理論指導[26]?；诖?借助生物信息學工具初篩病原體蛋白質序列中的候選表位已經成為疫苗研究中關鍵的一步。在此次新型冠狀病毒肺炎（COVID-19）發生之初,Grifoni團隊[19]就使用Net-MHCpan EL web工具對SARS-CoV-2蛋白質序列中潛在的T、B細胞表位進行了預測,經同源比對發現,這些表位與IEDB中記錄的SARS表位有較高的同源性（B細胞表位為69%～100%;T細胞表位為17%～100%）。由此,我們綜合利用當前的生物信息學工具,對SARS-CoV-2結構蛋白中潛在的T細胞候選表位進行了預測和篩選,共獲得34條符合MHCⅠ類分子的T細胞候選表位信息,其中候選表位 S1110～1121和 M135～144覆蓋高達91%的中國人群和71%的世界人群,可為COVID-19的疫苗設計提供參考。

目前,利用機器學習算法,相關人員已經開發了系列針對抗原提呈過程的生物信息學工具。比如,利用人工神經網絡算法開發的NetMHCpan系列工具,既擁有操作簡潔的web頁面,也有適用于Linux系統的版本。另外,本課題組利用隨機森林算法開發的INeo-Epp工具可以預測抗原肽的免疫原性,其AUC達到0.81,而IEDB中提供的同類型工具的AUC僅為0.69。除本研究中涉及的工具外,還有諸多類似功能的工具可用于疫苗候選表位的預測和篩選,例如:NetCTL 1.2（http://www.cbs.dtu.dk/services/NetCTL/）可用于預測蛋白質序列中潛在的細胞毒性T淋巴細胞（cytotoxic T lymphocyte,CTL）表位[27];MixMHC2pred[28]可用于預測HLAⅡ類分子呈遞的抗原肽。這些生物信息學預測工具在抗原特定功能的預測中均顯示出良好的性能,但在實現疫苗候選表位的完整預測和分析方面仍存在一定難度。

綜上所述,本研究綜合利用多種生物信息學工具,結合自主研發的算法工具,整理并提出了一個可廣泛應用于病原體多肽疫苗候選表位預測和分析的流程。通過該工作流程可以預測出不同抗原蛋白質序列中高可信度的疫苗候選表位,從而為預防和控制感染性疾病與惡性腫瘤的新型疫苗的研制提供參考。