Perl語言整合5款軟件進行真菌效應蛋白預測

馮世鵬

海南省熱帶生物資源可持續利用重點實驗室，海南大學 農學院，海南 海口 570228

依據植物對病原物的抗病機制不同，將植物抗病性分為被動抗性（passive resistance）和主動抗性（active resistance）。被動抗性（也稱物理化學抗性）是植物自身性狀所決定的抗病性，如植物通過結構屏障（如角質層、蠟質等）和體內的化學成分（如次生代謝產物等）阻止病原物浸入[1]。主動抗性指植物在受到病原物侵染或誘導劑處理時誘導激活一系列防衛反應，抑制病原物浸染，包括一系列生理生化的變化和抗病相關基因的激活表達，如侵染部位的木質化、富含羥脯氨酸的糖蛋白增加、植保素（phytoalexins）積累等[2]。

植物的主動抗病性主要借助于植物的先天免疫反應，可分為2 種。第一種是基于細胞表面的模式識別受體（pattern-recognization receptors，PRR）對病原物相關分子模式（pathogen-associated molecular patterns，PAMP）的識別，該免疫過程被稱為PAMP觸發的免疫（PAMP triggered immunity，PTI），該免疫過程可幫助植物抵抗大部分病原微生物的入侵；第二種免疫發生在細胞內部，主要依靠抗病基因（resistance gene，R gene）編碼的蛋白產物直接或間接地識別病原微生物產生的效應因子（effector）并激發免疫反應，也稱為效應因子觸發的免疫（effector triggered immunity，ETI）。

廣義的效應因子是指病原物可引起宿主細胞結構和功能改變的蛋白或小分子物質，它們促進侵染進程，或激發免疫反應，或兩者都有。廣義的效應因子包括PAMP、毒素、降解酶等[3]，其中蛋白類是效應因子中的重要組成部分，也稱為效應蛋白（effector proteins）[4]。不同物種效應蛋白結構有差異，如細菌通過Ⅲ型分泌系統（type Ⅲsecretion system，T3SS）分泌效應蛋白[5-6]，卵菌的效應蛋白含有RXLR 結構域[7]等。

真菌效應蛋白的特點是：①無典型的結構域。雖然有些研究發現部分效應蛋白有2 個或多個亮氨酸，以形成如二硫鍵這類結構，但整體來說，真菌效應蛋白沒有共有的結構特征[8]，效應蛋白之間及其與其他蛋白間的相似性均不高，因此在進行預測時常與Swiss-prot 數據庫進行blast 比對，無顯著相似結構蛋白作為潛在效應蛋白的條件之一。②屬于胞外分泌蛋白。蛋白分泌至細胞外主要有2 條途徑，一是經典途徑，通過蛋白質N 端信號肽序列將其引導分泌至胞外，約90%的分泌蛋白通過這種途徑分泌至胞外[9]；另一條是非經典途徑，少數分泌蛋白通過該途徑分泌至胞外[10]。在進行效應蛋白預測時經常考慮的是經典途徑分泌的蛋白，主要通過蛋白N 端信號肽來預測分泌蛋白。

有多款軟件可進行分泌蛋白的預測，每款軟件的預測方法及功能各有側重，Caccia等[11]對這些軟件及其作用進行了總結。將各軟件按照預測的結果來分類，可分為如下幾類：①蛋白信號肽及其剪切位點預測，常見的預測軟件如SignalP[12]，其他功能類似的軟件有PrediSI[13]、SigCleave[11]、SpScan[11]、Phobius[14-15]、Philius[16]、RPSP[17]等；②非經典途徑分泌蛋白預測，常見的預測軟件如SecretomeP[18]，其他功能類似的軟件有SOSUI-GramN[19]，NClassG+[20]等；③蛋白亞細胞定位預測，常見的如ProtComp[11]、TargetP[21]、WoLF PSORT[22]，其他有ChloroP[23]等；④跨膜蛋白預測，常見的預測軟件有TMHMM[24]、TMpred[25]。

對分泌蛋白的預測使用一款軟件如SignalP[26]或者WoLF PSORT[27]，或者聯合使用多款軟件，多款軟件有不同的組合形式[28-33]。

綜合已有文獻報道的分泌蛋白軟件預測使用情況，結合效應蛋白的特征，本研究擬用Perl語言整合SignalP4.1、ProtComp、TargetP1.1、TMHMM2.0、WoLF POSRT等5款軟件進行真菌效應蛋白預測。

1 材料和方法

1.1 測試數據來源

香蕉枯萎病菌（Fusarium oxysporumf.sp.cubense race 1）為鐮刀菌屬的真菌，基因組大小約50 Mb，其蛋白質組包含15438 條蛋白序列（BioProject Accession：PRJNA174274），在NCBI 下載這些蛋白序列進行測試。

1.2 真菌潛在效應蛋白確定

本文將無典型結構域的胞外分泌蛋白定義為潛在的真菌效應蛋白。

真菌效應蛋白無典型結構域，大多數效應蛋白的功能未知。在NCBI中登錄的蛋白，其功能有實驗驗證的，將其實驗驗證結果列出；暫未進行實驗驗證的，會通過nr、Swiss-prot、KEGG、interpro、COG 等數據庫的blast 搜索，根據相似蛋白的功能來推測未知蛋白的功能（即蛋白注釋）；其中功能既無實驗驗證結果又與其他數據庫蛋白無顯著相似性的蛋白常被命名為“hypothetical protein（預測蛋白，其功能未知）”。根據這種特點，本文擬通過名稱中出現的“hypothetical protein”這一名詞區分與其他蛋白無顯著相似性的蛋白，這樣無須進行相關數據庫的blast 搜索（蛋白注釋過程也完成這一步，因此可省略）。對于研究者自行測序的新物種，先對所預測蛋白進行注釋，注釋完的數據也可利用本文提供的方法進行后續分析。

真菌效應蛋白為分泌蛋白，這類蛋白含有信號肽序列，定位于細胞外，并非跨膜蛋白。通過軟件SignalP 預測蛋白N 端含有信號肽序列，通過3 款亞細胞定位軟件ProtComp、TargetP、WoLF PSORT預測蛋白均定位于胞外，通過軟件TMHMM 預測不是跨膜蛋白。

1.3 預測軟件介紹

信號肽預測軟件SignalP4.1。通過神經網絡法進行N 端信號肽及其剪切位點預測，包含5 個得分參數，即C值（剪切位點原始分值）、S值（信號肽序列分值）、Y 值（C 值幾何均值與S 值的聯合得分值）、mean S（S 值的算術均值）、D 值（mean S 值與最大Y值的權重均值）。其中D 值用于區分信號肽（SP=‘YES’）與非信號肽（SP=‘NO’）。本研究使用參數D值“YES”，預測蛋白含有N 端信號肽序列。SignalP是應用最廣泛的預測軟件之一，Choi 等[26]測試發現其對分泌蛋白的預測準確性達89%，因此可單獨用于蛋白分泌組預測。安裝軟件版本signalp-4.1c.Linux.tar.gz（下載地址：http://www.cbs.dtu.dk/services/SignalP/）。

亞細胞定位預測軟件ProtComp。可預測蛋白定位于何種細胞器或區域，如Nuclear（核）、Plasma membrane（質膜）、Extracellular（胞外）、Cytoplasmic（細胞質）、Mitochondrial（線粒體）、Endoplasm.retic.（內質網）、Peroxisomal（過氧化物酶體）、Lysosomal（溶酶體）、Golgi（高爾基體）、Vacuolar（液泡）。該軟件使用5 種預測方法，即LocDB（與已知定位蛋白同源比對）、PotLocDB（與理論分析定位蛋白同源比對）、NNets（神經網絡法）、Pentamers（查詢序列與數據庫序列的pentamers 分布比較法）和Integral（將上述4 種方法分值進行最終網絡法整合），預測結果供研究者自行判斷。本研究選用Integral 值“Extracellular”預測蛋白定位于胞外。安裝軟件版本protcompan.tar.bz2（下載地址：http://linux1.softberry.com/berry.phtml?topic=index&group=programs&subgroup=proloc）。

亞細胞定位預測軟件TargetP1.1。依據N 端的一定序列長度預測真核生物蛋白質定位（Loc）于葉綠體（C，含葉綠體轉運肽cTP）、線粒體（M，含線粒體轉運肽mTP）、含信號肽分泌蛋白（S，含信號肽序列SP），或者其他地方（-）。通過預測的最高值與第二高值之間的差異分值（diff）進行可信度分區，1 代表diff＞0.8，2 代表0.8＞diff＞0.6，3 代表0.6＞diff＞0.4，4 代表0.4＞diff＞0.2，5代表diff＜0.2。本研究使用參數Loc值“S”，預測蛋白含信號肽序列并分泌至胞外。安裝軟件版本targetp-1.1b.Linux.tar.Z（下載地址：http://signalfind.org/tatfind.html）。

亞細胞定位軟件WoLF PSORT。可預測蛋白定位于extr（extracellular，胞外）、cysk（cytoskeleton，細胞骨架）、E.R.（endoplasmic reticulum，內質網）、golg（Golgi apparatus，高爾基體）、mito（mitochondria，線粒體）、nucl（nucleus，細胞核）、plas（plasma membrane，細胞膜）、pero（peroxisome，過氧化物酶體）、vacu（vacuolar membrane，液泡膜）、chlo（chloroplast，葉綠體）、lyso（lysozyme，溶酶體），并對可能的定位情況進行打分，分值越高定位于該處的可能性越高。本研究使用extr 值＞17，預測蛋白定位于胞外。O'Connell等[27]發現其對胞外分泌蛋白預測準確性達92.7%，因此單獨用其進行真菌分泌組蛋白預測。安裝軟件版本WoLFPSORT_package_v0.2（下載地址：http://wolfpsort.org/）。

跨膜蛋白預測軟件TMHMM2.0。預測蛋白跨膜螺旋的數量（Number of predicted TMHs，PredHel）和跨膜區氨基酸長度（Exp number of AAs in TMHs，ExpAA），同時給出查詢蛋白長度值（len）、前60 個氨基酸所預測跨膜螺旋數（First60）、及拓撲結構（Topology）等參數。ExpAA 長度大于18 則證明有跨膜結構或者有信號肽。本研究使用參數ExpAA值＜18，預測蛋白不是跨膜蛋白。安裝軟件版本tmhmm-2.0.Linux.tar.gz（下載地址：http://www.cbs.dtu.dk/services/TMHMM/）。

所使用電腦為Dell 工作站T5600，內存32G，操作系統CentOS6.5，系統自帶的Perl 版本5.10.1，所有軟件的安裝根據相關軟件說明進行，安裝完畢將相關命令輸出至環境變量，保證該命令可直接調用。

相關軟件及其參數見表1。

1.4 使用Perl語言進行軟件整合

1.4.1 方法一：讓各軟件對蛋白質組進行平行預測再取結果交集 將蛋白質組分為各個蛋白質后，每一款軟件都對任一蛋白進行預測，各軟件的排列順序對該方法最終運行時間及最終預測結果均無影響（圖1）。用該方法可獲得每一款軟件對整個蛋白質組完整的預測結果（*.signalp.out、*.protcomp.out、*.targetp.out、*.tmhmm.out、*.wolfpsort.out）及最后的各軟件預測結果的交集，即潛在的效應蛋白（*.effector.list、*.effector.fa、*.effector.info），其中*.effector.list 文件存放的是潛在效應蛋白名稱清單，*.effector.fa 文件存放fasta格式的潛在效應蛋白序列，*.effector.info文件存放的是5 款軟件對各潛在效應蛋白預測結果的匯總。此方法的優點是可獲得所有軟件預測結果，缺點是運行時間明顯延長。

1.4.2 方法二：讓各軟件按順序進行預測，邊運行邊取交集 將蛋白質組分為各個蛋白質后，讓各軟件按順序對所取各蛋白進行預測，邊預測邊取前面各軟件預測結果的交集，如果其中任一款預測軟件所預測結果不符合參數設定，則終止其后的軟件對該蛋白的進一步預測分項，這樣可快速獲得最終的交集結果（圖2），其中*.signalp.final.out、*.protcomp.out、*.targetp.out、*.tmhmm.out、*.wolfpsort.out 這5 個文件存放的是各軟件對潛在效應蛋白的預測結果，這些結果與方法一的結果的區別是：前者保存的各軟件對整個蛋白質組完整的預測結果，后者保存的是各軟件對最終潛在效應蛋白的預測結果。*.effector.list、*.effector.fa、*.effector.info 這3 個文件存放的內容與方法一同名文件存放相同的內容（表2）。在相關軟件的順序安排上，由于各軟件運行速度有差異，會導致最終運行時間的差異，因此將運行速度快的軟件或方法盡可能排在前面（表2），由于SignalP是預測信號肽蛋白最準確，也是目前預測分泌蛋白使用最多的軟件，因此將其排列在除Linux系統命令之外的其他各軟件之前，WoLF PSORT 是另一款可單獨用于分泌蛋白預測的軟件，將其放在最后，這樣最終軟件排列順序依次為Linux 系統命令、SignalP、ProtComp、TMHMM、TargetP、WoLF PSORT。

圖1 方法一預測流程圖（各軟件平行預測）

2 結果與討論

用香蕉枯萎病菌蛋白質組（共有蛋白序列15438 seqs）進行2 種方法的測試，2 種方法均可預測最終潛在效應蛋白348 seqs。

圖2 方法二預測流程圖（各軟件先后預測）

方法一在獲得最終潛在效應蛋白的同時，可以獲得其他軟件單獨運行的結果，預測流程及結果見圖1，預測過程耗時8 h 25 min。具體結果如下：①用Linux 命令找出“hypothetical protein”7672 seqs；②用SignalP 軟件預測含信號肽的蛋白1494 seqs；③用ProtComp 軟件預測胞外分泌蛋白1668 seqs；④用TMHMM 軟件預測非跨膜蛋白12387 seqs；⑤用TargetP 軟件預測胞外蛋白2595 seqs；⑥用WoLF PSORT 軟件預測胞外蛋白1027 seqs；⑦各軟件預測結果交集即為最終潛在效應蛋白348 seqs。

方法二進行各軟件的先后順序預測，前一款軟件預測結果不符合參數設置則不再進行后續軟件的預測，對香蕉枯萎病菌蛋白質組15438 seqs 進行預測，流程圖見圖2，整個過程耗時40 min，相對于方法一大大縮短了預測時間，而所預測的最終潛在效應蛋白一致（348 seqs）。預測過程是這樣的：①首先用Linux 命令找出其中的“hypothetical protein”7672 seqs，去掉非“hypothetical protein”7766 seqs；②然后用SignalP 軟件預測“hypothetical protein”7672 seqs 中含信號肽序列的蛋白828 seqs，去掉不含信號肽序列的蛋白6844 seqs；③用ProtComp 軟件預測828 seqs 中定位于胞外的蛋白539 seqs，去掉定位于其他地方的蛋白289 seqs；④用TMHMM 軟件預測539 seqs 中不含跨膜結構的蛋白505 seqs，去掉含跨膜結構域的蛋白34 seq；⑤用TargetP 軟件預測505 seqs 中定位與胞外的蛋白501 seqs，去掉定位于其他地方的蛋白4 seqs；⑥最后用WoLF PSORT 軟件預測501 seqs 中定位于胞外的蛋白348 seqs，去掉153 seqs，這樣獲得最終潛在效應蛋白348 seqs。

由于真菌經實驗驗證的效應蛋白太少，無法進行本方法的準確性測試，但是通過比較現有的蛋白分泌組所使用的方法[28-33]，本方法與其類似；O'Connell 等[27]通過blast 比對及單獨使用WoLF PSORT 進行2 種炭疽病菌的潛在效應蛋白預測，結果分別為177（禾生炭疽菌，Colletotrichum graminicola）和365 seqs（白菜炭疽菌，C.higginsianum），與本方法的預測結果處于同一數量級。

表2 2種方法運行結果比較表

綜上，我們整合5 款軟件建立了真菌效應蛋白的2 種預測方法，方法一在獲得最終潛在效應蛋白的同時，可獲得各軟件單獨的預測結果，但耗時長；方法二可快速獲得最終潛在效應蛋白預測結果。2種方法均可用于真菌效應蛋白預測，可由研究者選擇使用（2種方法的Perl程序附在文后）。

[1]Hückelhoven R.Transport and secretion in plant-microbe interactions [J].Curr Opin Plant Biol,2007,10(6):573-579.

[2]van Loon L C,Bakker P A,Pieterse C M.Systemic resistance induced by rhizosphere bacteria[J].Annu Rev Phytopathol,2006,36:453-483.

[3]Hogenhout S A,van der Hoorn R A,Terauchi R,et al.Emerging concepts in effector biology of plant associated organisms[J].Mol Plant Microbe Interact,2009,22(2):115-122.

[4]Ellis J,Rafiqi M,Gan P,et al.Recent progress in discovery and functional analysis of effector proteins of fungal and oomycete plant pathogens[J].Curr Opin Plant Biol,2009,12:399-405.

[5]McCann H C,Guttman D S.Evolution of the type III secretion system and its effectors in plant-microbe interactions[J].New Phytol,2008,177(1),33-47.

[6]Block A,Li G,Fu Z Q,et al.Phytopathogen type III effector weaponry and their plant targets[J].Curr Opin Plant Biol,2008,11,396-403.

[7]Dou D,Kale S D,Wang X,et al.RXLR-mediated entry of Phytophthora sojae effector Avr1b into soybean cells does not require pathogen-encoded machinery[J].Plant Cell,2008,20(7),1930-1947.

[8]Houterman P M,Speijer D,Dekker H L,et al.The mixed xylem sap proteome of Fusarium oxysporum-infected tomato plants[J].Mol Plant Pathol,2007,8:215-221.

[9]Tsang A,Butler G,Powlowski J,et al.Analytical and computational approaches to define the Aspergillus niger secretome[J].Fungal Genet Biol,2009,46(Suppl1):S153-S160.

[10]Hirai Y,Nelson C M,Yamazaki K,et al.Non-classical export of epimorphin and its adhesion to avintegrin in regulation of epithelial morphogenesis[J].J Cell Sci,2007,120(12):2032-2043.

[11]Caccia D,Dugo M,Callari M,et al.Bioinformatics tools for secretome analysis[J].Biochim Biophys Acta,2013,183(11):2442-2453.

[12]Petersen T N,Brunak S,Heijne G,et al.SignalP 4.0:discriminating signal peptides from transmembrane regions[J].Nat Methods,2011,8:785-786.

[13]Hiller K,Grote A,Scheer M,et al.PrediSi:prediction of signal peptides and their cleavage positions[J].Nucleic Acids Res,2004,32:W375-W379.

[14]K?ll L,Krogh A,Sonnhammer E L L.A Combined transmembrane topology and signal peptide prediction method[J].J Mol Biol,2004,338(5):1027-1036.

[15]K?ll L,Krogh A,Sonnhammer E L L.Advantages of combined transmembrane topology and signal peptide prediction--the Phobius web server[J].Nucleic Acids Res,2007,35:W429-W32.

[16]Reynolds S M,K?ll L,Riffle M E,et al.Transmembrane topology and signal peptide prediction using dynamic bayesian networks[J].PLoS Comput Biol,2004,4(11):e1000213.

[17]Plewczynski D,Slabinski L,Ginalski K,et al.Prediction of signal peptides in protein sequences by neural networks[J].Acta Biochim Pol,2008,55(2):261-267.

[18]Bendtsen J D,Jensen L J,Blom N,et al.Feature based prediction of non-classical and leaderless protein secretion[J].Protein Eng Des Sel,2004,17(4):349-356.

[19]Imai K,Asakawa N,Tsuji T,et al.SOSUI-GramN:high performance prediction for sub-cellular localization of proteins in gram-negative bacteria[J].Bioinformation,2008,2(9):417-421.

[20]Restrepo-Montoya D,Pino C,Nino L F,et al.NClassG+:a classifier for non-classically secreted Gram-positive bacterial proteins[J].BMC Bioinformatics,2011,12:21.

[21]Emanuelsson O,Nielsen H,Brunak S,et al.Predicting subcellular localization of proteins based on their N-terminal amino acid sequence[J].J Mol Biol,2000,300(4):1005-1016.

[22]Horton P,Park K J,Obayashi T,et al.WoLF PSORT:protein localization predictor[J].Nucleic Acids Res,2007,35:W585-W587.

[23]Emanuelsson O,Nielsen H,von Heijne G.ChloroP,a neural network-based method for predicting chloroplast transit peptides and their cleavage sites[J].Protein Sci,1999,8(5):978-984.

[24]Sonnhammer E L,von Heijne G,Krogh A.A hidden Markov model for predicting transmembrane helices in protein sequences[J].Proc Int Conf Intell Syst Mol Biol,1998,6:175-182.

[25]Hofmann K,Stoffel W.TMbase -a database of membranespanning proteins segments[J].Biol Chem Hoppe Seyler,1993,374:166.

[26]Choi J,Park J,Kim D,et al.Fungal secretome database:integrated platform for annotation of fungal secretomes[J].BMC Genomics,2010,11:105.

[27]O'Connell R J,Thon M R,Hacquard S,et al.Lifestyle transitions in plant pathogenic Colletotrichum fungi deciphered by genome and transcriptome analyses[J].Nat Genet,2012,44(9):1060-1065.

[28]王黎,吉愛國.大腸埃希氏桿菌UTI89 分泌蛋白組的預測分析[J].中國生物工程雜志,2007,27(7):100-105.

[29]周曉罡,丁玉梅,張紹松,等.馬鈴薯青枯菌(Ralstonia solanacearum)染色體與質粒分泌蛋白組比對分析[J].中國食用菌,2008,27(suppl.):12-15.

[30]李午佼，顏瑾，王運生.南方和北方根結線蟲分泌蛋白組預測與比較分析[J].湖南農業大學學報(自然科學版),2011,37(4):376-380.

[31]Lum G,Min X J.FunSecKB:the Fungal Secretome Knowledge Base[J].Database,2011,2011:bar001.

[32]Druzhinina I S,Shelest E,Kubicek C P.Novel traits of Trichoderma predicted through the analysis of its secretome[J].FEMS Microbiol Lett,2012,337(1):1-9.

[33]陳相永,陳捷胤,肖紅利,等.植物病原真菌寄生性與分泌蛋白組CAZymes 的比較分析[J].植物病理學報,2014,44(2):163-172.

[34]Emanuelsson O,Nielsen H,Brunak S,et al.Predicting subcellular localization of proteins based on their N-terminal amino acid sequence[J].J Mol Biol,2000,300:1005-1016.

[35]Nielsen H,Engelbrecht J,Brunak S,et al.Identification of prokaryotic and eukaryotic signal peptides and prediction of their cleavage sites[J].Protein Eng,1997,10:1-6.

附件1：方法一所用Perl程序

附件2：方法二所用Perl程序