周加興,唐旭清,b*

(江南大學(xué)a.理學(xué)院;b.無錫市生物計(jì)算工程技術(shù)研究中心,中國江蘇無錫214122)

甲型流感是由病毒引起的一種高致病率和高死亡率的傳染病,每年造成大約500萬人致病和50萬人死亡[1]。20世紀(jì),暴發(fā)了由不同基因型引起的3次大規(guī)模的流感大流行,分別是1918年的H1N1(“西班牙流感”)、1957年的H2N2(“亞洲流感”)和1968年的H3N2(“香港流感”)。21世紀(jì)初,一種新的甲型流感(“豬流感”)在北美洲暴發(fā),并迅速傳遍世界各地。從2009年4月至2010年8月,官方確認(rèn)致死病例達(dá)18 500例,而且這個(gè)數(shù)字很可能只是死亡病例中的一部分[2,3]。每次流感病毒的暴發(fā),對人類和畜類的生命造成嚴(yán)重的危害,對社會生產(chǎn)和生活帶來諸多不便,所以對流感大流行建立模型進(jìn)行預(yù)測分析迫在眉睫。

目前,已有許多學(xué)者在流感流行預(yù)測方面進(jìn)行了研究。Banerjee等[4]利用堿基間的鍵能大小,由血凝素(hemagglutinin,HA)蛋白質(zhì)發(fā)生抗原漂移的難易程度進(jìn)而預(yù)測流感發(fā)生。Yang等[5]用貝葉斯推理結(jié)合流行病學(xué)模型對香港地區(qū)1998-2013年間的流感暴發(fā)情況進(jìn)行了綜合分析。楊冬紅等[6]分析了流感暴發(fā)與厄爾尼諾現(xiàn)象的關(guān)系。Jacobs等[7]分析了流感暴發(fā)與感染人群年齡之間的關(guān)系。此外,Dugas等[8]基于谷歌流感模型,利用廣義線性回歸分析,結(jié)合流感的時(shí)序數(shù)據(jù),對流感暴發(fā)進(jìn)行了預(yù)測。

隨著生物技術(shù)的發(fā)展,大量蛋白質(zhì)序列通過實(shí)驗(yàn)被測定。而采用實(shí)驗(yàn)的方法去預(yù)測蛋白質(zhì)的結(jié)構(gòu)和功能會耗費(fèi)大量的人力物力。近年來,隨著理論方法的成熟,理論分析方法預(yù)測蛋白質(zhì)結(jié)構(gòu)和功能的應(yīng)用越來越廣泛。薛曉麗等[9]基于H1N1的HA蛋白質(zhì)進(jìn)化樹揭示了H1N1病毒進(jìn)化關(guān)系。靳佩軒等[10]基于流感病毒10種組成蛋白質(zhì)的氨基酸序列,對未來流感病毒的變異和暴發(fā)進(jìn)行了預(yù)測。Deng等[11]針對不同亞型的HA和NA,利用系統(tǒng)進(jìn)化樹對洞庭湖地區(qū)的禽流感進(jìn)行了深入的分析。這些研究顯示,采用HA蛋白質(zhì)序列的相似性來研究流感病毒蛋白質(zhì)功能和變異關(guān)系是一種有效可行的方法。

生物多樣性是指在確定空間內(nèi)生物的豐富程度及變異程度,生物多樣性越大,就可能出現(xiàn)更優(yōu)良、更有生命力的子代。國內(nèi)部分科研人員在通過評價(jià)生物多樣性來預(yù)測生物在環(huán)境中的生存能力和適應(yīng)能力方面進(jìn)行了較為深入的研究。彭麗潭等[12]通過分層聚類和信息融合分析了氣候?qū)Φろ旡Q種群的影響。李桑等[13]通過DGGE技術(shù)研究了鄉(xiāng)村沼氣池污泥微生物的生物多樣性對細(xì)菌微生物群落的影響。本文通過從生物多樣性的角度解釋病毒的變異、進(jìn)化等關(guān)系,從而為流感暴發(fā)研究提供依據(jù)和支撐。

從NCBI下載了1902-2016年99 861條流感病毒的HA,進(jìn)行40維蛋白質(zhì)特征提取,將蛋白質(zhì)序列轉(zhuǎn)換成數(shù)值信息。同時(shí),對每年的數(shù)據(jù)進(jìn)行層次聚類,引入最優(yōu)層次聚類評價(jià)指標(biāo),獲取每一年數(shù)據(jù)的最優(yōu)聚類數(shù)。最后,通過計(jì)算每一年的種群熵評價(jià)指標(biāo)來度量每一年流感病毒的生物多樣性,并計(jì)算種群熵變化率以衡量流感病毒的變異速率。

1 材料與方法

1.1 數(shù)據(jù)來源

本文所需數(shù)據(jù)均從NCBI(http://www.ncbi.nlm.nih.gov/gennomes/FLU/Database/nph-select.cgi)網(wǎng)站中Molecular Databases的Protein Sequence上下載。數(shù)據(jù)包括1902-2016年的99 861個(gè)流感病毒 HA蛋白質(zhì)序列,其中 1903-1917、1920、1921、1922、1923、1924、1926、1928、1929、1932、1941、1944等年份HA蛋白質(zhì)序列的數(shù)據(jù)缺失。

1.2 特征提取

特征提取是數(shù)據(jù)處理與序列比對的關(guān)鍵步驟,蛋白質(zhì)序列比對現(xiàn)行的方法主要基于氨基酸序列的組成頻率與位置關(guān)系。本文采用融合蛋白質(zhì)組成、氨基酸理化屬性及耦合信息的40維特征向量。如給定長度為n的蛋白質(zhì)序列記為S=s1s2…si,si(i=1,2,…,n)表示20種基本氨基酸中的一種,對每一條蛋白質(zhì)序列構(gòu)造40維的特征向量V(C),即V(C)=(V1(G),V2(F),V3(H))[14]。其中,V1(G)表示HA蛋白質(zhì)序列中各種氨基酸出現(xiàn)的頻數(shù)構(gòu)成的20維向量;V2(F)表示氨基酸的理化性質(zhì)分成W1=(R,D,E,N,Q,K,H)非極性且親水、W2=(L,I,V,A,M,F)非極性且疏水、W3=(S,T,Y,W)極性且親水性、W4=(P,G,C)極性且疏水性等四類,每類氨基酸在HA蛋白質(zhì)序列中出現(xiàn)的頻數(shù)構(gòu)成的4維向量;V3(H)表示HA蛋白質(zhì)序列中按理化屬性分為四類之后,氨基酸兩兩相連出現(xiàn)的頻數(shù)構(gòu)成的16維向量。

1.3 最優(yōu)聚類評價(jià)指標(biāo)

聚類分析旨在使得相似度較高的聚成同類,相似度較小的劃為不同類,即同類之間的相似度高,不同類之間的相似度低。本文采用層次聚類處理流感病毒的HA蛋白質(zhì)序列信息,定義蛋白質(zhì)序列x1,x2,…，xm間的差異[15]為:

其中,xij為樣本xi的第j個(gè)屬性。Vm和Vn是HA蛋白質(zhì)序列的兩類,如果Vm和Vn的相似度最大,則將這兩類合并。合并后新的聚類表示為:

對于層次聚類建立的分層遞階結(jié)構(gòu),確定合適的粒度是關(guān)鍵問題[16,17]。假設(shè)d是可數(shù)集X={x1,x2,…,xn}的度量標(biāo)準(zhǔn)(如歐氏距離等)。存在一個(gè)粒度空間X(λ)??Td(X),記X(λ)={a1,a2,…,acλ},ak={xk1,聚類中心ak且為樣本中心。引入類內(nèi)偏差Sin(X(λ))和類間偏差Sbetween(X(λ)),如下:

最優(yōu)層次評價(jià)指標(biāo)(HEI)為:

聚類最優(yōu)評價(jià)指標(biāo)被用來確定合理的分類和粒度。當(dāng)Sbetween(X(λ))＜Sin(X(λ))時(shí),HEI(X(λ))取最小值,取得合理的粒度,其目標(biāo)優(yōu)化模型如下:

1.4 生物多樣性

種群熵值是評價(jià)生物多樣性的指標(biāo)之一。假設(shè)P(t)為第t代種群,種群中個(gè)體的數(shù)量為N,根據(jù)個(gè)體之間的差異可將種群劃分為m個(gè)部分[18]:P1(t),P2(t),…,Pm(t)。顯然,并且對于均有。設(shè)k1,k2,…,km分別為P1(t),P2(t),…,Pm(t)中個(gè)體數(shù)目,定義第t代種群的熵其中Pi=ki/N。

1.5 算法流程

Step 1 對流感病毒數(shù)據(jù)進(jìn)行預(yù)處理,并剔除掉部分流感病毒。

Step 2 將蛋白質(zhì)序列數(shù)值化轉(zhuǎn)化成40維特征向量,并計(jì)算每一年與第一年數(shù)據(jù)的歐氏距離,然后繪出基于歐氏距離的流感病毒變異進(jìn)化分布圖。

Step 3 對所得40維特征向量進(jìn)行基于歐氏距離的層次聚類,根據(jù)HEI(X(λ),確定每一年的最優(yōu)聚類數(shù)。

Step 4 基于最優(yōu)聚類數(shù),計(jì)算每一年種群熵值和種群熵變化率并繪出種群熵的變化圖。

2 結(jié)果與分析

首先對流感病毒的HA數(shù)據(jù)進(jìn)行預(yù)處理,按年份對其進(jìn)行40維特征向量提取。因?yàn)閿?shù)據(jù)過少,不利于評價(jià)當(dāng)年的生物多樣性分析,所以剔除掉 1919、1925、1936、1937、1952、1955 等年份的數(shù)據(jù)。然后,利用算法計(jì)算每一年的種群熵值(表1),并繪出種群熵變化圖(圖1)。同時(shí),計(jì)算每一年的數(shù)據(jù)與第一年(1902年)數(shù)據(jù)的歐氏距離并繪出病毒的變異進(jìn)化分布圖(圖2)。最后,計(jì)算出20世紀(jì)種群熵變化率top10(表2)和21世紀(jì)種群熵變化率 top5(表3)。

表1 1902-2016年的種群熵值Table1 The population entropy from 1902 to 2016

表2 1902-2000年的種群熵變化率top10Table2 Top 10 years in change rate of the population entropy from 1902 to 2000

圖1 1902-2016年種群熵變化圖Fig.1 The population entropy variation from 1902 to 2016

從圖1中歷年信息熵值的變化可以宏觀地看出,信息熵整體呈現(xiàn)上升趨勢,表明了流感病毒生物多樣性越來越豐富。這現(xiàn)象符合實(shí)際中生物發(fā)展的規(guī)律,因?yàn)椴《镜牟粩嘧儺惡椭亟M產(chǎn)生了新的基因型和亞種,使得流感病毒的生物多樣性變得越來越豐富,即種群熵的值越來越大。20世紀(jì)以來,1918-1919年暴發(fā)了由H1N1引發(fā)的著名的“西班牙流感”,據(jù)初步統(tǒng)計(jì),導(dǎo)致2 000～5 000萬人死亡[19,20]。從表1可以看出,1918年流感病毒的種群熵值為0.702 9,明顯比1902年與1927年等鄰近年份的病毒種群熵值高。1955-1957年暴發(fā)了“亞洲流感”,導(dǎo)致全球大約280萬人死亡[21]。表1顯示,1956年的種群熵值為1.258 0,1957年的種群熵值為1.550 2,均也明顯高于與它鄰近的年份,同樣表明1956年和1957年的流感病毒生物多樣性豐富。1968年暴發(fā)了由甲型流感病毒H3N2所致的“香港流感”,1965年的種群熵為1.319 7,1968年的種群熵為1.305 1,也明顯高于其他鄰近的年份,說明流感病毒變異重組幾年前就已經(jīng)開始,一直處于潛伏期,直到1968年暴發(fā)[22]。1977年、1987年和1997年都發(fā)生了不同規(guī)模的流感暴發(fā)事件[23,24]。進(jìn)入21世紀(jì),2009年暴發(fā)了由新型甲型H1N1引發(fā)的“墨西哥流感”,迅速傳遍全球?qū)е轮辽?8 500人死亡。從表1中可看出,2005-2009年的種群熵值分別為2.328 6、2.438 5、2.456 5、2.453 0、2.376 3,都高于鄰近的年份,說明流感病毒從2005年開始出現(xiàn)較大的變異,直到2009年暴發(fā)成流感大流行[25]。同樣2012-2014年的種群熵也高于鄰近年份,查閱文獻(xiàn)顯示在2013-2014年暴發(fā)了H1N1流感季節(jié)性流行[26]。

上述分析顯示,20世紀(jì)發(fā)生流感大流行的間隔時(shí)間約在10年左右,而進(jìn)入21世紀(jì)后流感大流行發(fā)生的間隔時(shí)間縮短至4～5年左右。因?yàn)椴《緩淖儺惖奖┌l(fā)具有一定的潛伏期,本文選擇對20世紀(jì)每一年的種群熵與5年前的種群熵值作差,21世紀(jì)每一年的種群熵與3年前的種群熵值作差。再將其與時(shí)間間隔的比值記為病毒變異的變化率。

圖2 流感病毒變異進(jìn)化分布圖Fig.2 Mutation distribution of influenza viruses

表3 2000-2016年的種群熵變化率top5Table3 Top 5 years in change rate of the population entropy from 2000 to 2016

表2顯示在1956年、1963年、1976年和1998年等年份的種群熵變化率較大,即病毒的變異速率快,變異明顯。相對,在病毒較明顯的變異之后發(fā)生了1957年、1968年、1977年和1997年的流感大流行。表3顯示在2006年、2007年和2014年等年份的種群熵變化率較大,即病毒變異速率快,變異明顯。同樣,在病毒較明顯的變異之后發(fā)生了2009年和2014年流感的大流行。

計(jì)算出每一年的數(shù)據(jù)與第一年的數(shù)據(jù)的歐氏距離,并繪出相對第一年的變異進(jìn)化圖(圖2),可以更容易看出流感病毒的變異情況。結(jié)合表1和圖2發(fā)現(xiàn),1902-1970年的種群熵值不高且變異程度也比較低,而1970-2016年的種群熵值較高且變異程度很活躍。但實(shí)際上,1918年、1957年以及1968年的3次流感大流行的波及范圍廣、死亡率高、破壞力強(qiáng);相比之下,1977年、1997年以及2009年的大流行死亡率較低、破壞力也較弱。1918年“西班牙流感”的高死亡率主要集中在15～35歲的青年之間,原因在于1889-1890年暴發(fā)了由不同基因型引發(fā)的“俄羅斯流感”,導(dǎo)致這個(gè)年齡段的人失去了對由H1N1病毒引發(fā)的流感的免疫能力[27,28]。1957年由H2N2引發(fā)的“亞洲流感”病毒攜帶的HA抗原和NA抗原從未在人體中出現(xiàn)過,它是由人H1N1亞型與禽流感病毒的3個(gè)基因片段重組而來。從圖3流感病毒系統(tǒng)發(fā)育樹也可以看出,1918年和1957年的流感病毒同源性較低,所以1957年暴發(fā)流感時(shí)人類體內(nèi)沒有相關(guān)的抗體,從而導(dǎo)致全球各地280萬人死亡[29]。而1968年由H3N2引發(fā)的“香港流感”,其流感病毒的N2亞型與1957年的H3N2一樣,這樣多數(shù)人體內(nèi)存在N2表面蛋白抗體;此外,由于該病毒擁有與1889年大流行的H3N8病株同樣的H3抗原,這樣部分老年人也擁有對H3N2的免疫能力,從而使死亡人數(shù)降低,約750 000人死亡[22,30]。

對于1970年后暴發(fā)的流感基因型,如1977年的H1N1、1997年的H5N1、2009年的 H1N1以及2014年的H1N1等。綜合前文分析及圖3可知,由于前面流感的暴發(fā)使得人體內(nèi)已經(jīng)擁有部分重組或變異流感基因型的抗體,加上預(yù)測系統(tǒng)的成熟建立和醫(yī)療條件的改善,使得1970年后流感引起的死亡率和破壞力沒有前面的幾次高。

3 結(jié)束語

本文通過對99 861個(gè)流感病毒的HA蛋白質(zhì)序列進(jìn)行特征提取,采用層次聚類方法并引入最優(yōu)層次評價(jià)指標(biāo)計(jì)算出每一年的最優(yōu)聚類數(shù)。利用每一年的種群熵值刻畫流感病毒的生物多樣性,從宏觀角度能夠較好吻合歷史數(shù)據(jù)。進(jìn)一步,通過變異進(jìn)化分布圖和種群熵變化率對流感病毒的變異進(jìn)行深入分析。數(shù)據(jù)分析發(fā)現(xiàn),種群熵值能很好地反應(yīng)流感病毒的生物多樣性,種群熵變化率也能很好地反應(yīng)流感病毒的變異和不穩(wěn)定情況。通過流感病毒種群熵值和種群熵變化率對發(fā)生流感大流行的時(shí)間進(jìn)行粗略的估計(jì),表明在21世紀(jì)流感病毒大流行的發(fā)生時(shí)間間隔會縮短且規(guī)模和破壞性都將降低。這些研究可為流感的預(yù)測提供依據(jù)和支撐。

圖3 流感大暴發(fā)年份的流感病毒系統(tǒng)發(fā)育樹Fig.3 Phylogenetic tree of influenza viruses in epidemic years

參考文獻(xiàn)(References):

[1]Xu R,Ekiert D C,Krause J C,et al.Structural basis of preexisting immunity to the 2009 H1N1 pandemic influenza virus[J].Science,2010,328(5976):357-360.

[2]Dawood F S,Iuliano A D,Reed C,et al.Preliminary estimates of global 2009 H1N1 influenza mortality[C]//Infectious Diseases Society of America 2011 Annual Meeting.Arlington:Infectious Diseases Society of America,2011.

[3]NovelSwine-OriginInfluenzaA(H1N1)VirusInvestigationTeam,Dawood F S,Jain S.Emergence of a novel swine-origin influenza a(H1N1)virus in humans[J].The New England Journal of Medicine,2009,360(25):2605-2615.

[4]Banerjee R,Roy A,Das S,et al.Similarity of currently circulating H1N1 virus with the 2009 pandemic clone:viability of an imminent pandemic[J].Infection,Genetics&Evolution,2015,32:107-112.

[5]Yang W,Cowling B J,Lau E H,et al.Forecasting influenza epidemics in Hong Kong[J].PLoS Computational Biology,2015,11(7):e1004383.

[6]楊冬紅,楊學(xué)祥.流感世界大流行的氣候特征[J].沙漠與綠洲氣象(Yang Dong-hong,Yang Xue-xiang.The climatic characteristic of pandemic influenza[J].Desert and Oasis Meteorology,2007,1(3):1-8.

[7]Jacobs J H,Archer B N,Baker M G,et al.Searching for sharp drops in the incidence of pandemic A/H1N1 influenza by single year of age[J].PLoS One,2012,7(8):e42328.

[8]Dugas A F,Jalalpour M,Gel Y,et al.Influenza forecasting with Google flu trends[J].PLoS One,2013,8(2):e56176.

[9]薛曉麗,李陽,唐旭清.基于H1N1型禽流感病毒的HA蛋白序列進(jìn)化樹研究[J].計(jì)算機(jī)應(yīng)用研究(Xue Xiao-li,Li Yang,Tang Xu-qing.Research on evolutionary tree for H1N1 flu virus based on HA protein sequences[J].Application Research of Computers),2015,32(9):2634-2638.

[10]靳佩軒,高潔.流感病毒組成蛋白質(zhì)序列的分析與預(yù)測[J].食品與生物技術(shù)學(xué)報(bào)(Jin Pei-xuan,Gao Jie.Sequence analysis and prediction of the influenza virus protein[J].Journal of Food Science and Biotechnology),2016,35(4):393-398.

[11]Deng G,Tan D,Shi J,et al.Complex reassortment of multiple subtypes of avian influenza viruses in domestic ducks at the Dongting lake region of China[J].Journal of Virology,2013,87(17):9452-9462.

[12]彭麗潭,吳軍,唐旭清.氣候變化對丹頂鶴種群在繁殖地逗留時(shí)間的影響分析[J].生態(tài)與農(nóng)村環(huán)境學(xué)報(bào)(Peng Li-tan,Wu Jun,Tang Xu-qing.Impact of climate change on stay of redcrowned cranes in their breeding habitat[J].Journal of Ecology and Rural Environment),2014,30(3):280-288.

[13]李桑,張琳,邱義蘭,等.鄉(xiāng)村沼氣池污泥微生物多樣性的研究[J].生命科學(xué)研究(Li Sang,Zhang Lin,Qiu Yi-lan,et al.Study on microbial diversities in rural biogas digesters[J].Life Science Research),2015,19(4):321-327.

[14]李巍巍,李陽,唐旭清.不同特征描述下H1N1病毒蛋白序列的比較[J].生命科學(xué)研究(Li Wei-wei,Li Yang,Tang Xuqing.Comparing the H1N1 flu virus protein sequences by different feature vectors[J].Life Science Research),2016,20(2):119-124.

[15]李陽,唐旭清.基于粗?；牧鞲胁《镜鞍走M(jìn)化樹構(gòu)建[J].模式識別與人工智能(Li Yang,Tang Xu-qing.Construction of phylogenetic treeof flu virus proteins based on coarse graining[J].Pattern Recognition and Artificial Intelligence),2016,29(10):936-942.

[16]Tang X Q,Zhu P.Hierarchical clustering problems and analysis of fuzzy proximity relation on granular space[J].IEEE Transactions on Fuzzy Systems,2013,21(5):814-824.

[17]Tang X Q,Zhu P,Cheng J X.The structural clustering and analysis of metric based on granular space[J].Pattern Recognition,2010,43(11):3768-3786.

[18]申元霞,張翠芳.一種新型保持種群多樣性的遺傳算法[J].系統(tǒng)仿真學(xué)報(bào)(Shen Yuan-xia,Zhang Cui-fang.A modified genetic algorithm with maintaining diversity[J].Journal of System Simulation),2005,17(5):1052-1053.

[19]周劍芳,楊磊,藍(lán)雨,等.1918/1919年西班牙流感(H1N1)病原學(xué)概述[J].病毒學(xué)報(bào)(Zhou Jian-fang,Yang Lei,Lan Yu,et al.Epidemiological overview of 1918/1919 Spanish influenza[J].Chinese Journal of Virology),2009,25(suppl.):8-11.

[20]Worobey M,Han G Z,Rambaut A.Genesis and pathogenesis of the 1918 pandemic H1N1 influenza a virus[J].Proceedings of the National Academy of Sciences USA,2014,111(22):8107-8112.

[21]Duff F L.Pandemic influenza in 1957:review of international spread of new Asian strain[J].Journal of the American Medical Association,1958,166(10):1140-1148.

[22]袁帆,藍(lán)雨,郭俊峰,等.1968年流感大流行的流行病學(xué)概述[J].病毒學(xué)報(bào)(Yuan Fan,Lan Yu,Guo Jun-feng,et al.Epidemiological overview of the 1968 influenza pandemic[J].Chinese Journal of Virology),2009,25(suppl.):33-35.

[23]Singh G,Oberoi M S,Kwatra M S,et al.Isolation of influenza virus from horses in the equine influenza outbreak of 1987[J].Current Science,1987,56:1285-1286.

[24]Rocchi G,De Felici A,Ragona G,et al.Influenza activity in metropolitan Rome,Italy,during the cold-weather months of 1976-1977[J].Developments in Biological Standardization,1977,39:425-428.

[25]Dawood F S,Iuliano A D,Reed C,et al.Estimated global mortality associated with the first 12 months of 2009 pandemic influenza a H1N1 virus circulation:a modelling study[J].The Lancet Infectious Diseases,2012,12(9):687-695.

[26]Linderman S L,Chambers B S,Zost S J,et al.Potential antigenic explanation for atypical H1N1 infections among middleaged adults during the 2013-2014 influenza season[J].Proceedings of the National Academy of Sciences USA,2014,111(44):15798-15803.

[27]Brownlee G G,Foder E.The predicted antigenicity of the haemagglutinin of the 1918 Spanish influenza pandemic sug gests an avian origin[J].Philosophical Transactions of the Royal Society B-Biological Sciences,2001,356(1416):1871-1876.

[28]Reid A H,Taubenberger J K,Fanning T G.Evidence of an absence:the genetic origins of the 1918 pandemic influenza virus[J].Nature Reviews Microbiology,2004,2(11):909-914.

[29]Kawaoka Y,Krauss S,Webster R G.Avian-to-human transmission of the PB1 gene of influenza a viruses in the 1957 and 1968 pandemics[J].Journal of Virology,1989,63(11):4603-4608.

[30]Cox N J,Subbarao K.Global epidemiology of influenza:past and present[J].Annual Review of Medicine,2000,51(1):407-421.