甲型流感病毒HA蛋白質序列的預測

張 玲， 高 潔

（江南大學 理學院，江蘇 無錫 214122）

流感是一種傳染病，它具有一定的周期性和反復性，其在全球引起的發病率和死亡率非常高[1-2]。流感病毒基因組由8個獨立的RNA片段組成，分別編碼多個與病毒結構和病毒復制有關的蛋白質分子。其中，最受關注的蛋白分子是血凝素（hemagglutinin，HA）。血凝素蛋白與病毒易感的宿主范圍和宿主對病毒感染產生的免疫反應有直接聯系[3]。且蛋白質是生命的物質基礎，沒有蛋白質就沒有生命[4]。因此，許多學者和專家研究甲型流感病毒并取得了一些成就。

他們從不同的角度研究流感病毒序列，且在尋找著最佳模型。李曉燕等對新型H1N1流感病毒HA基因特性分析[5]。基于聚類算法ModuleFind，梅娟等通過最大化蛋白質網絡的模塊性來尋找具有較強蛋白質結構的劃分[6]。劉娟、任迪用時間序列的方法及整數階差分和分數階差分預測DNA序列[7-8]。基于CGR-游走模型，作者用分數階差分的ARFIMA模型預測HA蛋白質序列。對所選取的1943～2013年同源性相對較高的71條蛋白質序列，我們用 ARFIMA（p，d，q）模型對前 20個位置去擬合并且預測除極個別外由預報區域圖顯示原始數據都在預報區域內，表明模型建立的比較合理，預報效果很好，這對流感病毒的研究和預測有著重要的意義。

1 方法

1.1 基于詳細HP模型的蛋白質序列的CGR-游走模型

CGR是一種迭代映射技術[9]，2004年，喻祖國等人提出了基于詳細HP模型的蛋白質序列的CGR方法[10]。在詳細的HP模型中，將20種氨基酸分為四大類:非極性、負極性、無電荷極性和正極性，分別用 np、nep、up、pp 來表示，np={A，I，L，M，F，P，W，V}，nep={D，E}，up={N，C，Q，G，S，T，Y}，pp={R，H，K}。對于一個給定長為n的蛋白質序列S=s1s2Ksn，其中 si，i=1，2，…，n 是 20 種氨基酸中的一種，定義如下:

則得到一條序列 X（s）=c1c2…cn，其中 ci∈{A0，A1，A2,A3}。

再定義序列X（s）的CGR，類似于DNA序列的CGR定義，通過迭代過程得到該序列CGR CGRi=CGRi-1-0.5·（CGRi-1-ci），i=1，…，n，CGR0=（0.5，0.5）。

對于一個蛋白質序列，定義tk=yk/xk，其中yk是CGRk的y坐標值，xk是CGRk的x坐標值，則得到一個數據序列{tk:k=1，2，…，n}，把它作為一個時間序列，并稱它為“CGR-游走序列”。

1.2 ARFIMA模型

定義 1 {εt}稱為白噪聲（white noise）序列[11]，簡記為Xt～WN（μ，σ2）。 如果時間序列滿足如下性質:

（1）任取 t∈T，有 EXt=u；

定義3 如果隨機過程{Xt}是平穩的，且滿足差分方程 Φ（B）▽dxt=Θ（B）εt，其中{εt}為白噪聲序列，Eεt=0，Eεt2=σ2ε＜∞Θ（B）=1-φ1B-K-φpBp為 p 階自回歸系數多項式；Θ（B）=1-θ1B-K-θqBq，為 q 階移動平均系數多項式，-0.5＜d＜0.5，則稱{Xt}服從-0.5＜d＜0.5 的ARFIMA（p，d，q）模型[4]。

2 數據分析

選取同源性相對較高的71條流感病毒HA蛋白質序列（1943～2013年），數據來自 NCBI網站，其網址:http://www.ncbi.nlm.nih.gov/.

對于這71條蛋白質序列，選每一條序列的第三個位置為例來如下分析，首先把HA蛋白質序列轉換成數值，定義:

得到相對應的HA蛋白質數值序列為:[0 0 0 1 0 0 3 0 0 0 3 0 2 2 0 3 2 0 3 0 0 0 2 3 3 0 0 2 2 2 0 3 3 0 0 0 0 0 3 2 2 2 0 3 0 0 0 0 3 0 0 2 0 0 0 0 0 3 0 0 3 2 0 2 0 3 0 3 3 3 2]，再 CGR 混沌游走，經過計算得到分數階差分d=0.151，再取對數及0.151階差分，可以得到一組新的數值序列:[0 0 0 0-2.833 -2.833 -2.833 2.026 6 2.026 6 2.026 6 2.026 6 4.852 5 4.852 5 0.233 6 0.084 16 0.084 16 1.322 5 0.361 1 0.361 1 1.569 1.569 1.569 1.569 2 0.067 0 1.109 5 1.940 23 1.940 23 1.940 23 0.168 8 0.060 7 0.026 64 0.026 64 1.194 6 2.060 9 2.060 9 2.060 9 2.060 9 2.060 9 2.060 9 5.700 8 0.412 71 0.157 0.070 55 0.070 55 1.211 4 1.211 4 1.211 4 1.211 4 1.211 4 4.337 4.337 4.337 0.121 16 0.121 16 0.121 16 0.121 16 0.121 16 0.121 16 4.174 7 4.174 7 4.174 7 6.356 4 0.446 03 0.446 03 0.106 61 0.106 61 1.461 4 1.461 4 2.839 6 3.754 3 4.542 2]，這 71 條序列的各個位置都要經過這樣的預處理。

然后把這一組新的HA蛋白質序列轉化成如圖1所示的時序圖 （CGR游走和取對數0.151階差分），圖中顯示蛋白質序列在零上下波動，呈現出基本的平穩性。

圖1 第三個位置對數差分后的時序圖Fig.1 Time series figure of the third position sequence

圖 2（ACF）和（PACF）為第三位置序列取對數再0.151階差分后的自相關函數圖和偏自相關函數圖。序列的自相關圖衰減迅速，而序列的偏自相關圖衰減緩慢，這意味著原序列具有長記憶特征。

圖2 對數0.151階差分序列的自相關函數圖和偏自相關函數圖Fig.2 Sample ACF and PACF of the 0.151 order differenced log data of protein sequences

我們以純粹的隨機測試的分數差分序列，結果見表1。P＜0.000 1＜0.05，所以經取對數分數階差分后此序列不是白噪聲序列。因此，我們可以用ARFIMA（p，d，q）模型擬合該序列。

根據Akaike信息判別準則[12]，可選擇ARFIMA（4，0.151，0）模型來擬合此每條序列的第三個位置序列。表2給出了模型的參數估計。

表1 白噪聲檢驗Table 1 Autocorrelation check for white noise

表2 參數最小二乘估計Table 2 Conditional least squares estimation

從表中可以看出參數的P值都小于0.1，這表明 ARFIMA（4，0.151，0）模型能有效地擬合該序列。為檢驗該模型的合理性，選擇了一個合適的檢驗統計量LB檢驗統計量[13]:

其中，rk是滯后的樣本自相關函數，n是樣本容量，M是一個取定的比n小的正整數。

表3顯示了對于各滯后階數，LB統計量的P值均顯著大于0.1，意味著擬合模型的殘差序列應為白噪聲 （純隨機），因而可認為ARFIMA（4，0.151，0）模型能合理有效地擬合該序列。

表3 殘差的自相關檢驗Table 3 Autocorrelation check of residuals

表 4則利用ARFIMA （4，0.151，0） 模型對 HA序列進行短期預測，即預測后十年2014～2023年第三個位置的預報值。（注:這些預報數據是混沌游走后取對數差分得到的，所以我們可以推出相應的數值。）

表4 2014-2023年的第三個位置的預報值Table 4 Forecast value of the third position for 2014-2023 year

從圖3可以看出，除極個別外由預報區域圖顯示原始數據都在預報區域內，表明該模型建立的比較合理，預報效果很好。

根據上面對蛋白質序列的第三個位置的選擇模型、參數估計、LB統計量檢驗，我們可以預測2014-2023年第三個位置的HA蛋白質序列的預報值為:3 2 1 0 1 1 1 1 2 2，進而可以預知2014-2023年第三個位置HA蛋白質是非極性、負極性、無電荷極性還是正極性。

圖3 預報區域圖Fig.3 Forecast figure

3 結語

基于CGR混沌游走模型和分數階差分模型，作者采取是一種縱向的預測方法，即用ARFIMA（p，d，q）模型預測未來年甲型流感病毒HA蛋白質序列。以1943～2013年這71條蛋白質序列的第三個位置為例，我們得到用 ARFIMA（p，d，q）模型對其前20個位置去擬合并且預測，發現除極個別外原始數據都在預報區域內除極個別外，表明模型建立的比較合理，預報效果很好。

我們可以用此方法分析和研究預測未來年的流感病毒蛋白質序列，這樣節省大量的精力和財力。當然這種方法有一定的缺點，如選定的位置不足夠多等，會影響預測值的精確性，還有待改進。

[1]Morens D，Folkers G，Fauci A.The challenge of emerging and re-emerging infectious diseases[J].Nature，2004，430:242-249.

[2]Muzaffar S B，Ydenberg R C，Jones I L.Avian influenza:an ecological and evolutionary perspective for waterbirdd scientists[J].Waterbirds，2006，29:243-257.

[3]Kobasa D，Takada A，Shinya K，et al.Enhanced virulence of influenza A viruses with the hae-magglutinin of the 1918 pandemic virus[J].Nature，2004，431（7017）:703-707.

[4]GAO Jie，XU Zhen-yuan.Chao game representation （CGR）-walk model for DNA sequences[J].Chinese Physics B，2009，18（11）:370-376.

[5]李曉燕，孔梅，陳錦英，等.新型H1N1流感病毒HA基因特性分析[J].中國衛生檢驗雜志，2010，20（12）:3121-3124.LI Xiao-yan，KONG Mei，CHEN Jin-ying.Analysis on the HA gene characteristics of novel influenza A （1H1N1） virus[J].Chinese Journal of Health Laboratory Technology，2010，20（12）:3121-3124.（in Chinese）

[6]梅娟，何正，王正祥，等.基于網絡模塊性的蛋白質序列聚類[J].食品與生物技術學報，2010，29（1）:123-127.MEI Juan，HE Zheng，WANG Zheng-xiang.Clustering protein sequences through modularity optimization[J].Journal of Food Science and Biotechnology，2010，29（1）:123-127.（in Chinese）

[7]劉娟，高潔.甲型H1N1流感病毒DNA序列堿基的預測[J].生物信息學，2011，9（3）:259-262.LIU Juan，GAO Jie.Forecasting bases for DNA sequences of influenza virus A/H1N1[J].China Journal of Bioinformatics，2011，9（3）:259-262.（in Chinese）

[8]REN Di，GAO jie.Early-warning signals for an outbreak of the influenza pandemic[J].Chin Phys B，2011，20（12）:128701-4.

[9]Jeffrey H J.Chaos game representation of gene structure[J].Nucleic Acid Res，1990，18:2163-2170.

[10]YU Zu-guo，Anh V V，Lau Ka-sing.Fractal analysis of measure representation of large proteins based on the detailed HP model[J].Physicas A，2004，337:171-184.

[11]王燕.應用時間序列分析[M].北京:中國人民大學出版社，2008.

[12]Akaike H.A new look at statistical model identification[J].IEEE transaction on Automatics Control，1974，19:9-14.

[13]Ljung G M，Box G E P.On a measure of lack of fit in time series models[J].Biomertrica，1978，30A:9-14.