甲型流感病毒DNA序列的長記憶ARFIMA模型*

劉 娟 高 潔

(江南大學理學院，無錫 214122)(2010年4月16日收到;2010年8月4日收到修改稿)

甲型流感病毒DNA序列的長記憶ARFIMA模型*

劉 娟 高 潔

(江南大學理學院，無錫 214122)(2010年4月16日收到;2010年8月4日收到修改稿)

流感病毒分為三類:甲型(A型)，乙型(B型)，丙型(C型).在這三種類型中甲型(A型)流感病毒是最致命的流感病毒，對人類引起了嚴重疾病.本文對甲型流感病毒DNA序列建立了一種新的時間序列模型，即CGR(Chaos Game Representation)弧度序列.利用CGR坐標將甲流病毒DNA序列轉(zhuǎn)換成CGR弧度序列，且引入長記憶ARFIMA模型去擬合此類序列，發(fā)現(xiàn)隨機找來的10條H1N1序列，10條H3N2序列都具有長相關(guān)性且擬合很好，并且還發(fā)現(xiàn)這兩種序列可以嘗試用不同的 ARFIMA模型去識別，其中 H1N1可用 ARFIMA(0，d，5)模型去識別，H3N2可用ARFIMA(1，d，1)模型去識別.

甲型流感，時間序列模型，CGR，ARFIMA(p，d，q)模型

PACS:87.10.Vg，02.50.Fz

1.引 言

流感是一種反復出現(xiàn)的傳染病，在全球引起了高發(fā)病率和高死亡率[]1.流感病毒分為三類:甲型(A型)，乙型(B型)，丙型(C型).在這三種類型中甲型(A型)流感(以下簡稱甲流)病毒是最致命的流感病毒，給人類帶來了嚴重的疾病.甲流病毒根據(jù)其表面的血凝素(hemagglutinin，HA)和神經(jīng)氨酸酶(neuraminidase，NA)基因的不同又可分成16個HA亞型(H1-H16)和9個NA亞型(N1—N9)，不同的HA和NA形成了甲流病毒的許多亞型，如H1N1，H3N2，H5N1 等等［2—4］.筆者參看了許多文獻幾乎沒有看到用時間序列模型來挖掘甲型流感病毒的特性的，因而本文采用時間序列模型來分析甲型流感病毒.

1992年，Peng等［5］提出了 DNA一維游走模型.同年Voss等［6］提出了不同的觀點，他們發(fā)現(xiàn) DNA序列的譜密度顯示的1/fβ噪聲無處不在，意味著當0<β<1存在長相關(guān)性，認為不僅在非編碼區(qū)序列中在編碼序列中也存在長相關(guān)性.另一方面Buldyrev 等［7，8］設(shè)計了一個廣義-Lévy 游走模型去生成一個模型序列，使用所有可用的DNA序列發(fā)現(xiàn)主要在非編碼序列中呈現(xiàn)長相關(guān)性.基于該模型，Tai等［9］提出了一個二維修正-Lévy游走模型.為區(qū)分 C和T，A和 G，Lou等建立了二維和三維游走模型［10］，Yu 等則建立了圖譜［11，12］，來研究 DNA 序列的相關(guān)性.2006 年，Lopes和 Nunes［13］引入長記憶ARFIMA(0，d，0)模型去擬合 DNA序列的一維游走序列.2009年，Gao等［14］基于 CGR(chaos game representation)坐標提出了一種DNA序列轉(zhuǎn)換成一個時間序列(CGR-游走序列)的方法，并引入長記憶ARFIMA(p，d，q)模型來分析.

本文對甲流病毒DNA序列提供了一種新的時間序列模型，即CGR弧度序列.利用CGR坐標將甲流病毒DNA序列轉(zhuǎn)換成CGR弧度序列，且引入長記憶ARFIMA模型去擬合此類序列，發(fā)現(xiàn)隨機找來的10條H1N1序列，10條H3N2序列都具有長相關(guān)性且擬合很好，并且還發(fā)現(xiàn)這兩種序列可以嘗試用不同的 ARFIMA模型去識別，其中 H1N1可用ARFIMA(0，d，5)模型去識別，H3N2可用 ARFIMA(1，d，1)模型去識別.

2.ARFIMA模型

如果隨機過程 {xt}是平穩(wěn)的，且滿足方程Φ(B)Δdxt= Θ(B)εt，其中，－ 0.5

則稱{xt}服從 －0.5

因此，{xt}可看作是分數(shù)差分噪聲導出的 ARMA(p，q)過程.當2d－1=－1時，d=0即為短記憶過程;所以當2d－1>－1時，d∈(0，0.5)具有長記憶的特征.

3.基于CGR的時間序列模型

1990年Jeffrey提出了一種 DNA序列可視化的方法即 CGR方法［15］.CGR是一種迭代映射技術(shù)，它把序列中的每個單元，如蛋白質(zhì)序列中氨基酸，DNA中的核苷酸，映射到一個連續(xù)的坐標空間中去.

正方形的四個頂點對應四種核苷酸.在這里，用DNA序列代替隨機數(shù)，每一個堿基的坐標都可以來確定下一個堿基的位置.我們?nèi)?A(0，0)，T(1，0)，G(1，1)，C(0，1)，并且取點(0.5，0.5)為起始點.

下面給出DNA迭代函數(shù)，也可以認為是 CGR算法 的 公 式 化 形 式［15，16］. 對 于 一 個 序 列 S =

其中 gi={(0，0)，(1，0)，(1，1)，(0，1)}，gi和 si相對應.

對于一個DNA序列，定義

其中yn是CGRn的y坐標值，xn是 CGRn的 x坐標值.則得到一個數(shù)據(jù)序列 {Rn:n=1，2，…，N}，我們把它作為一個時間序列，并稱它為“CGR弧度序列”.

以甲流病毒 H1N1序列 CY056890為例，數(shù)據(jù)來自 NCBI網(wǎng)站，其網(wǎng)址:http://www.ncbi.nlm.nih.gov/.

它的CGR弧度序列“游走圖”如下表1.

表1 CY056890序列所選部分前8個ACATGGTA游走結(jié)果

4.甲流 H1N1型病毒 CY056890的數(shù)據(jù)分析

圖1(a)是CY056890序列CGR弧度序列圖(位置380—2170)，樣本容量1791.這些數(shù)據(jù)變動較大，呈現(xiàn)非平穩(wěn)特征.考慮對此過程作d階差分.先對原序列作對數(shù)變換然后再做一階差分結(jié)果如圖1(b)所示，可見除少數(shù)地方呈現(xiàn)異方差外，基本呈現(xiàn)平穩(wěn)性.

圖2(a)(ACF)和圖2(b)(PACF)為樣本取對數(shù)再一階差分后的自相關(guān)函數(shù)圖形和偏自相關(guān)函數(shù)圖形.可見ACF衰減迅速，而PACF衰減緩慢，這意味著原序列具有長記憶特征.

圖1 (a)甲流H1N1型病毒CY056890的弧度序列圖;(b)取對數(shù)再一階差分圖

圖2 (a)取對數(shù)再一階差分的樣本自相關(guān)圖;(b)取對數(shù)再一階差分的樣本偏自相關(guān)圖

圖3給出了方差圖［17］是一個估計長記憶參數(shù)d的有用工具.對于一個長記憶時間序列{Rn}，它的均值珔Rk的方差滿足作log［Var(珔Rk)］關(guān)于log(k)的散點圖，對散點圖線性擬合，可估計得到線性方程的斜率為 －0.6877，令該斜率為2d－1= －0.6877，即可得 d的估計值0.156.

根據(jù)上述理由我們可選擇CGR弧度序列顯示長記憶特征.目的是利用上述特點為序列建立一個合適的模型.因此，可以考慮長記憶 ARFIMA(p，d，q)模型(d∈ (0，0.5))，p，q定階時為考慮實用性，僅考慮 p，q均小于等于5的 ARFIMA(p，d，q)模型.由 Akaike 信 息 判 別 準 則［18，19］，可 選 ARFIMA(0，0.156，5)模型來擬合.

為檢驗該模型的合理性，選擇了一個合適的檢驗統(tǒng)計量 LB 檢驗統(tǒng)計量［20，21］

其中rk是滯后k的樣本自相關(guān)函數(shù)，n是樣本容量，M是一個取定的比n小的正整數(shù).

表2顯示了對于各滯后階數(shù)，LB統(tǒng)計量的p值均顯著大于0.1，意味著擬合模型的殘差序列應為白噪聲(純隨機)，因而可以認為 ARFIMA(0，0.156，5)模型能很合理地擬合 CY056890序列的CGR-游走序列.

圖3 CY056890DNA序列的CGR弧度序列方差圖

表2 殘差的自相關(guān)檢驗

表3給出了被選擇的 ARFIMA(0，0.156，5)模型的參數(shù)估計，5個參數(shù)的 T檢驗統(tǒng)計量的p值均顯著小于 0.005.這意味著 ARFIMA(0，0.156，5)模型能有效地擬合這個CGR弧度序列.

表3 條件最小二乘估計

5.其余9條H1N1序列和10條H3N2序列數(shù)據(jù)分析

表4和表5分別給出了隨機選的9條H1N1序列和10條 H3N2序列的數(shù)據(jù)信息、被選擇的ARFIMA(p，d，q)模型及參數(shù)估計.從計算結(jié)果可得d均位于(0，0.5);對于各滯后階數(shù)，LB統(tǒng)計量的 p值除極個別外其余均顯著大于0.1;且每個被選擇的模型中各參數(shù)的T檢驗統(tǒng)計量的p值均顯著小于0.01.所有這些結(jié)果都顯示ARFIMA(p，d，q)模型能很合理很有效地擬合這些不同的CGR弧度序列且還發(fā)現(xiàn)所選H1N1序列均為ARFIMA(0，d，5)模型，所選 H3N2序列均為 ARFIMA(1，d，1)模型.所以我們可嘗試用 ARFIMA(0，d，5)模型，ARFIMA(1，d，1)模型分別去識別H1N1序列，H3N2序列.

表4 9條H1N1序列的數(shù)據(jù)信息、被選擇的ARFIMA模型和參數(shù)估計

表5 10條H3N2序列的數(shù)據(jù)信息、被選擇的ARFIMA模型和參數(shù)估計

6.結(jié) 論

本文基于CGR坐標提出了一種將甲流病毒DNA序列轉(zhuǎn)換成時間序列(CGR弧度序列)的方法，并引入長記憶模型ARFIMA模型來分析，首先分析了甲流H1N1型病毒CY056890序列，從圖1到圖3可知弧度序列顯示長記憶特征，并選擇了ARFIMA(0，0.156，5)模型去擬合它，從表 2到表 3發(fā)現(xiàn)擬合合理有效.

然后又分析了隨機找來的19條序列的CGR弧度序列，從表4和表5可知所有 ARFIMA(p，d，q)模型都有效合理.并且從表4中還發(fā)現(xiàn)所選H1N1序列均為 ARFIMA(0，d，5)模型，表5所選 H3N2序列均為 ARFIMA(1，d，1)模型.

由此可見，DNA序列的CGR弧度序列能由長記憶ARFIMA(p，d，q)模型有效合理地擬合，并且還可嘗試用 ARFIMA(0，d，5)模型，ARFIMA(1，d，1)模型分別去識別H1N1序列、H3N2序列.作為具有完善算法的經(jīng)典時間序列模型，不僅可以幫助我們得到甲流病毒DNA序列清晰的結(jié)構(gòu)，而且還可幫助我們有效識別甲流中的兩種亞型.

本文僅對甲流中的兩種亞型進行了研究分析，后面我們將研究分析甲流中的其他亞型以及乙型丙型流感病毒.

［1］ Morens D，F(xiàn)olkers G，F(xiàn)auci A 2004 Nature 430 242

［2］ Chen J M，Sun Y X，Liu S 2009 Chinese Science Bulletin 54 1657(in Chinese)［陳繼明、孫映雪、劉 朔 2009科學通報54 1657］

［3］ Webster R G，Bean W J，Gorman O T 1992 Microbiol Rev.56 152

［4］ Shi X M，Shi L，Zhang J F 2010 Chin.Phys.B 19 038701

［5］ Peng C K，Buldyrev S，Goldberg A L，Havlin S，Sciortino F，Simons M，Stanley H E 1992 Nature 356 168

［6］ Voss R F 1992 Phys.Rev.Lett.68 3805

［7］ Buldyrev S V，Goldberger A L，Havlin S，Peng C K，Simon M，Stanley H E 1993 Phys.Rev.E 47 4514

［8］ Buldyrev S V，Goldberger A L，Havlin S，Mantegna R N，Matsa ME，Peng C K，Simon M，Stanley H E 1995 Phys.Rev.E 51 5084

［9］ Tai Y Y，Li P C，Tseng H C 2006 Physica A 369 688

［10］ Luo L F，Lee W J，Jia L J，Ji F M，Tsai L 1998 Phys.Rev.E 58 861

［11］ Yu Z G，Chen G Y 2000 Theor.Phys.33 673

［12］ Yu Z G，Anh V，Gong Z M，Long S C 2002 Chin.Phys.11 1313

［13］ Lopes S R C，Nunes M A 2006 Physica A 361 569

［14］ Gao J，Xu Z Y 2009 Chin.Phys.B 18 370

［15］ Jeffrey H J 1990 Nucleic Acid Res 18 2163

［16］ AlmeidaJonas， carrico Joao A， Maretzek António 2001 Bioinformatics 17 429

［17］ Beran J 1994 Statistics for long-memory Processes(New Work:Chapman Hall)

［18］ Hosking J R M 1984 Water Resour.Res.20 1898

［19］ Crato N，Ray B K 1996 Journal of Forcasting 15 107

［20］ Ljung G M，Box G E P 1978 Biometrika 65 297

［21］ Li W K，Mcleod A I 1986 Biometrika 73 217

Long-memory ARFIMA model for DNA sequences of influenza A virus*

Liu Juan Gao Jie
(School of Science，Jiangnan University，Wuxi 214122，China)(Received 16 April 2010;revised manuscript received 4 August 2010)

Influenza viruses are divided into three types:A，B and C.Among them，type A virus is the most virulent human pathogen and causes the most severe disease.In this paper，we propose a new time series model for influenza A virus DNA sequence，i.e.chaos game representation(CGR)radians series.The CGR coordinates are converted into a time series model，and a long-memory ARFIMA(p，d，q)model is introduced to simulate the time series model.We select randomly 10 H1N1 sequences and 10 H3N2 sequences in analysis.we find in these data a remarkably long-range correlation and fit the model reasonably by ARFIMA(p，d，q)model，and also find that we can use different ARFIMA models to identify the two kinds of sequences，i.e.ARFIMA(0，d，5)model and ARFIMA(1，d，1)model that can identify H1N1 and H3N2 respectively.

influenza A virus，time series model，chaos game representation(CGR)，ARFIMA(p，d，q)model

.E-mail:ezhun6669@sina.com

*江南大學創(chuàng)新團隊發(fā)展計劃(批準號:2008CX002)中央高?；究蒲袠I(yè)務(wù)經(jīng)費專項資金(批準號:JUSRP21117)資助的課題.

.E-mail:ezhun6669@sina.com

*Project supported by the Innovative Research Team of Jiangnan University(Grant No.2008CX002)the Foundamental Research Founds for the Central Universities(Grant No.JUSRP21117).

PACS:87.10.Vg，02.50.Fz