蔡斌 彭謹 江華 楊浩 孫明偉 CharlesDamienLu 胡衛(wèi)建 曾俊

【摘要】目的 禽流感疫情的爆發(fā)和傳播受到多種自然因素的影響。今欲嘗試將地理信息系統(tǒng)與基因進化樹分析相結合，以建立一種基于基因序列變異追蹤中國禽流感病毒地理傳播的技術。方法 禽流感病毒基因來源于美國國立醫(yī)學圖書館（National Library for Medicine, NLM）數(shù)據(jù)庫，所獲得的基因組數(shù)據(jù)利用E-Utilities軟件包轉化為結構體后，可用Matlab軟件閱讀。結構體主要字段包括PB2、PB1、PA、HA、NP、HA、M1和NS1 8個片段，分別代表流感病毒的8個不同的基因片段?；诮Y構體字段，利用計算生物學的方法比較不同傳播能力禽流感病毒的同義突變/非同義突變基因（Ka/Ks）比例，確定不同選擇壓力之下A型禽流感病毒的基因突變模式。進而選擇Ka/Ks比例最大的基因片段，采用Jukes-Cantor算法估計氨基酸序列變異的進化距離，然后對不同爆發(fā)點的H5N1型禽流感進行進化樹聚類。將聚類信息輸入Google Earth，并利用不同圖層地理信息對影響爆發(fā)點分布的因素做單因素分析。結果 比較分析A型禽流感所有的8個基因序列可以看出，NS1、HA和NA蛋白的Ka/Ks比值較大。三者中，HA基因的Ka/Ks比值最大，可以代表病毒的傳播能力。利用分級聚類的思路對HA基因轉錄的氨基酸相似程度進行比較, 發(fā)現(xiàn)自2003年以來亞洲地區(qū)爆發(fā)的H5N1型禽流感之間的關系可以表示為一個由30個節(jié)點構成的進化樹，其中14個節(jié)點為分支節(jié)點，16個節(jié)點為葉子結點。把分支樹的前三個節(jié)點作為分類標準，可以把所有16個病毒株分為四類。這四類病毒在地理空間的分布呈現(xiàn)一定規(guī)律。計算發(fā)現(xiàn)禽流感爆發(fā)相關地理因素排序分別為：內(nèi)陸水體>主要鐵路交通線>家禽密度。結論 對中國H5N1病毒株基因序列變異的地理分布分析顯示，禽流感病毒爆發(fā)與候鳥遷徙、家禽運輸密切相關。

【關鍵詞】禽流感；病毒基因變異；谷歌地圖；地理信息系統(tǒng)；基因進化樹；同義突變/非同義突變基因比例；Jukes-Cantor算法；中國

Tracking the spread of avian influenza in China: a model based on evolutionary genetics analysis and geographic visualization CAI Bin, PENG Jin, JIANG Hua, YANG Hao, SUN Mingwei, Charles Damien Lu, Hu Wei-jian, ZENG Jun. Computational Biology Team, Metabolomics and Multidisciplinary Laboratory for Trauma Research, Sichuan Provincial Peoples Hospital, Sichuan Academy of Medical Sciences. Chengdu 610101,China

Corresponding author: JIANG Hua, Email: cdjianghua@gmail.com

【Abstract】Objective To explore the diverse natural and human factors affect the outbreak and spread of avian influenza. We integrated geographic visualization and evolutionary genetics technique to establish a method to track spread of avian influenza in China. Methods The sequence data of type A avian flu virus were provided by NCBI Nucleotide and Protein Databases. We transformed the original data to readable structures for Matlab using E-Utilities software. These MATLAB readable structures represented 8 genes of the virus, they are: RNA polymerase B2 (PB2), polymerase B1 (PB1), polymerase A (PA), hemagglutinin (HA), nucleoprotein (NP), neural aminidase (NA), matrix (M1), and non-structural (NS1) proteins. Based on these readable structures, we compared Ka/Ks ratio of different virus strains and identified the gene mutation patterns under different selection pressures. Then we selected the gene that exhibited the highest Ka/Ks ratio and performed a phylogenetic analysis by Jukes-Cantor algorithm. Google Earth layer tools were then used to integrate gene variation and geographic transmission information.Results When we compared these 8 virus genes, the NS1, HA and NA were found to exhibit high Ka/Ks ratio and could be seen to represent the transmission capacity of the virus. Among these, the HA gene has the highest Ka/Ks ratio. When we compare the amino acids encoded by the HA gene using clustering analysis, we found that the relationship between H5N1 avian influenza strains since 2003 in Asia made up an evolutionary tree. This evolutionary tree contained 30 nodes (14 branch nodes and 16 leaf nodes). All genes were classified into 4 major groups by the first 3 nodes. And these 4 groups exhibit clear geographic patterns in their spread. The impact of geographic factors on the outbreak of avian influenza in China can be ranked as: inland water bodies (lakes, reservoirs) > major railway paths > density of poultry.Conclusions The analysis on the dominant strainsgene mutations in China s H5N1 found that the outbreaks of avian influenza correlate with avian migration and poultry transportation.

【Key words】Avian flu； Gene Mutation；Google earth；Geographic Visualization；Gene Evolutionary Tree；Ka/Ks；Jukes-Cantor Algorithm； China

流行性感冒（influenza流感）是由流感病毒引起的急性呼吸道傳染病。該病病毒易發(fā)生抗原性變異，加之人群對變異株普遍易感，常迅速導致世界性的流感大流行^{［1- 2］}。在流感病毒的發(fā)病過程中，禽流感（avian flu）由于其在世界各國此起彼伏的暴發(fā)和流行，及可能與人流感病毒進行基因重組，成為研究的熱點。雖然迄今并沒有高致病性禽流感病毒能在人與人之間直接傳播的確切證據(jù)，但是由于兩種病毒的基因相似性使得人們擔心這種可能性正在逐漸增大。 流感病毒基因組主要由8個RNA片段組成。這8個片段為：聚合酶B2（polymerase B2, PB2）、聚合酶B1（polymerase B1）、聚合酶A（PA）、血凝素A（HA）、神經(jīng)氨酸酶（neuraminidase, NA）核蛋白（nucleoprotein, NP）、基質（M1）以及非結構蛋白（non-structural, NS1）^［3-5］。其中病毒的高度變異能力主要依賴于其病毒顆粒最外層的兩種表面抗原：血凝素抗原和神經(jīng)氨酸酶抗原^［6-7］。在致病力最強的A型禽流感當中，編碼HA和NA的基因由于氨基酸序列的不同又被分為HA（H1-16）以及NA（N1-9）幾種蛋白質亞型^{［8- 9］}。

基于進化生物學視角，禽流感變異對于流行病學的意義在于：（1）病毒基因的連續(xù)變異可能使得這種病毒獲得新宿主；（2）病毒基因的保守序列用以追蹤禽流感病毒來源。例如，WHO建立的Influza數(shù)據(jù)庫每年都在發(fā)布新的禽流感病毒基因序列變異測序?？紤]到這些變異都是發(fā)生在一定的時間和空間當中，候鳥遷徙、人類家禽運輸?shù)然顒佣加锌赡軐η萘鞲械膫鞑ズ妥儺惼鸬街匾淖饔?。這些活動很可能使得相隔很遠的兩個不同的地理位置出現(xiàn)相同基因型禽流感的爆發(fā)。另一方面，我們也可以觀察到地理上相鄰位置的流感病毒基因型可能出現(xiàn)明顯差異。由此，以下兩類問題——（1）“哪些地理空間中禽流感更易變異？”；（2）“不同基因型的禽流感的是否通過類似的途徑進行傳播？”——對于研究流感病毒的傳播特性和變異規(guī)律至為重要。要獲得這樣的知識，需要引入新的、能夠同時在基因變異和地理空間分布上對禽流感的情況進行追蹤和標記的研究工具。

1 資料與方法

1.1 研究的基本流程

研究的基本流程如下。（1）篩選出能夠正確標記流感病毒抗原變異的基因：筆者首先對比不同流感病毒基因片段的同義突變和非同義突變（Ka/Ks）比值，叢中篩選出非同義突變明顯偏高的基因片段，將這些基因片段的突變作為流感病毒獲得新的傳播能力的標記；（2）基因序列聚類以定義變異程度：利用進化樹、主成分分析等方法對這些基因片段進行聚類；（3）把不同聚類的病毒類型通過google earth API數(shù)據(jù)庫疊加到地圖上，與我們感興趣的地理分布特征進行相關性分析，以探索影響禽流感爆發(fā)流行的自然環(huán)境和地理相關因素。

1.2 禽流感病毒基因序列數(shù)據(jù)獲取與整理

禽流感病毒基因來源于美國國立醫(yī)學圖書館（National Library for Medicine, NLM）數(shù)據(jù)庫。下載選擇參數(shù)為：病毒分型=A，宿主=Avian, 國家/地區(qū)=Canada & Hong Kong，病毒片段=ALL，病毒基因表型=H5N1 & H2N3。下載1997和2001年兩年的數(shù)據(jù)。另外，還根據(jù)基因變異度分析的結果，選取病毒分型=A，宿主=Avian,國家和地區(qū)=Asia & African病毒片段=HA，年份=2001-2007年。將所獲的基因組數(shù)據(jù)利用美國國立衛(wèi)生研究院（National Institute for Health, NIH）開發(fā)的E-Utilities軟件包轉化為可由Matlab （版本7.10.0.499，Mathwork, USA）閱讀的結構體^［10］。結構體主要字段包括PB2、PB1、PA、HA、NP、HA、M1和NS1，分別代表流感病毒的8個不同的RNA片段。

1.3 計算Ka/Ks比值

其比較思路是首先對比所下載的基因組數(shù)據(jù)中的氨基酸序列和核酸序列，確定核酸序列密碼子序列。再利用上述密碼子序列分別計算同義突變和非同義突變的數(shù)量。利用非同義突變的發(fā)生概率確定選擇壓力，同時估算多重突變對結果進行校正（即：多次發(fā)生突變重新變回原來的密碼子）的概率。最終獲得各段基因的Ka/Ks比值。

1.4 HA基因進化樹和相似度聚類

Ka/Ks分析發(fā)現(xiàn)HA基因變異最能表征病毒的傳播能力。據(jù)此在上述數(shù)據(jù)庫中選取2001年以來，亞洲和非洲H5N1流感的HA基因片段測序數(shù)據(jù)。通過測量該基因表達的氨基酸序列變異程度，建立各段氨基酸序列相似程度的度量單位。利用Jukes-Cantor方法構建系統(tǒng)發(fā)育樹。獲得各地不同時間H1N1禽流感親緣關系數(shù)據(jù)圖。并且利用分級聚類的原理將禽流感分為4類。具體算法為：

d = －19/20 log(1－p ×20/19)

d為兩個氨基酸序列之間的距離， P為兩個序列氨基酸的相似程度。P=1意味著兩氨基酸序列完全相同，P=0意味著兩氨基酸序列構完全不同。

1.5 Google earth 地理信息系統(tǒng)疊加

Google Earth地理信息系統(tǒng)（http://www.google.com/earth/index.html）是一個可以通過API函數(shù)進行圖層疊加開發(fā)的開放式平臺，2005年Declan Butler為google earth增加了一個能夠標定流感基因分型和爆發(fā)人數(shù)的獨立圖層（人及禽流感）^［11-12］。在這個圖層基礎上，利用H1N1的禽流感分類信息疊加上人口聚集區(qū)、鐵路交通、養(yǎng)雞場分布、以及主要內(nèi)陸水體四條基本信息。為每一個爆發(fā)點構建一個向量，如果爆發(fā)點周圍存在上述任何一種地理單元，則為向量中的對應單元數(shù)據(jù)賦值1；反之則為0。采用單因素分析對禽流感爆發(fā)點與潛在傳播影響因素之間的相關程度進行估計。

2 結果

2.1 追尋禽流感病毒基因的變異程度

Ka/Ks比值分析是一種有效跟蹤和比較病毒基因組中有效和無效突變的技術。同義突變增加表明生物選擇壓力增大，但是同義突變本身并不會改變翻譯蛋白質的氨基酸序列。非同義突變則表明蛋白質氨基酸序列會發(fā)生改變，病毒可能獲得新的抗原特征從而造成爆發(fā)流行。本研究中，比較分析A型禽流感所有的8個基因序列可以看出，NS1、HA和NA蛋白的Ka/Ks比值最大（圖1）。這說明該三個基因所受到的選擇壓力最大，是決定病毒流行的主要因素,該發(fā)現(xiàn)與已有研究一致。而三者中的HA基因所發(fā)生的改變可能使病毒具有了在除鳥類以外的其他物種中傳播的能力。因此，HA基因的變異可作為表征病毒傳播能力的生物學變化。

2.2HA蛋白的系統(tǒng)發(fā)育學分析

利用Jukes-Cantor算法估計氨基酸序列之間的距離，獲得16個基因兩兩間的歐氏距離，總共可以獲得120個距離。利用分級聚類的思想對氨基酸之間的相似程度進行比較?？梢园l(fā)現(xiàn)自2003年以來亞洲地區(qū)爆發(fā)的H5N1型禽流感之間的關系可以表示為一個由30個節(jié)點構成的進化樹，其中14個節(jié)點為分支節(jié)點，16個節(jié)點為葉子結點。把分支樹的前三個節(jié)點作為分類標準，可以把所有16種基因分為4類。類似的，利用氨基酸序列相關性的120個距離也可以在三維空間中獲得16個基因的位置關系（圖2）。

2.3 地圖數(shù)據(jù)庫功能疊加

利用Google Earth選取高密度居民點、內(nèi)陸水體、主要鐵路交通線、養(yǎng)雞場密度等四個指標在地圖上進行標記，同時對進化樹上有親緣關系的禽流感病毒株與上述區(qū)域的重疊度進行相關性分析。結果發(fā)現(xiàn)各HA基因組分型禽流感病毒株與地理空間分布的相關性見圖3。

3 討論

禽流感爆發(fā)受到包括候鳥遷徙、家禽養(yǎng)殖和運輸?shù)榷喾N自然與人類活動因素的影響^［4-5］。如何對這些影響進行定量分析，一直以來都是禽流感流行病學研究中的難點。傳統(tǒng)的方法很難對相同血清型

的禽流感病毒進行進一步分類，也難于對病毒在地理空間中的傳播和變異進行追蹤。上述難題的解決，需要引入新的研究思路?，F(xiàn)代計算生物學技術在近十年取得的長足發(fā)展，已經(jīng)使得對禽流感病毒基因型的細微變異進行精確統(tǒng)計分析和聚類成為可能。同時，地理信息數(shù)據(jù)庫系統(tǒng)及其可視化技術的發(fā)展，使得有可能將病毒的空間分布和傳播信息與其生物學信息相結合。本研究正是在這一思路下，就禽流感在中國的傳播規(guī)律，建立了首個基于計算生物學-地理信息系統(tǒng)方法的模型。

1997年和2001年，甲型流感病毒的亞型（H5N1）在香港爆發(fā)，同其他禽流感不一樣的是，這次爆發(fā)的流感病毒具備了初步的禽-人的傳染能力，并導致了6人死亡。香港特區(qū)政府對病毒進行了全基因組測序，并向世界衛(wèi)生組織提供了病毒基因組全序列^［13］。為了比較和跟蹤禽流感的變異情況，筆者利用另外一種變異相對穩(wěn)定的禽流感病毒：加拿大衛(wèi)生部門保存的1985年和1977年的阿爾伯塔鴨（A/H2N3，Alberta duck）禽流感病毒^［14］來進行對照。通過對比兩種禽流感的基因序列的變異情況，探尋H5N1更加容易感染人類的原因。

已有研究發(fā)現(xiàn)，禽流感的7個基因在流感傳播過程中具有不同的功能。HA蛋白可以幫助病毒與宿主細胞粘附在一起，并且進入到宿主細胞內(nèi)部，NA蛋白可以剪切新生成的病毒，幫助病毒從宿主細胞中逸出進入新的健康細胞^［15］。不同的HA/NA蛋白組合對于流感病毒對不同物種呼吸道上皮細胞的侵襲性有不同影響。病毒學中，一個已知的重要事實是：H5N1雞型禽流感和N2N3鴨型禽流感之間最重要的差異，在于前者已經(jīng)初步具有跨種屬間的傳播能力^［16］。本研究從分子進化的角度印證了上述結論：傳播力最強的H5N1雞禽流感的每一段基因的Ka/Ks比值都較傳播力弱的H2N3野鴨禽流感要大，說明H5N1處在一個選擇壓力更高的環(huán)境。對比兩種禽流感病毒全部基因的Ka/Ks比值，可以發(fā)現(xiàn)，HA蛋白最大。這說明HA蛋白的變異最能說明兩種病毒之間的差異。所以，無論是就蛋白質的功能還是基因的變異程度而言，HA蛋白序列的差異都是追蹤不同病毒株之間基因差異的最好標記。

為了追蹤不同時間和地點有親緣關系的病毒株之間的關系。筆者將亞洲和非洲不同地區(qū)不同時期從雞身上分離的H5N1禽流感病毒的HA蛋白進行系統(tǒng)發(fā)育學的分析，分析它們之間相互的關系，并且利用人口聚居點、內(nèi)陸水體、主要鐵路交通線、養(yǎng)雞場的地理參數(shù)與禽流感是否爆發(fā)進行相關性分析。通過系統(tǒng)發(fā)育學比較，可以把H5N1型禽流感病毒株分為以下四種地理類型：1）河北—香港型，2）吉林—朝鮮—橫濱型，3）庫爾干（西伯利亞）—阿富汗—尼日爾型，以及4）河南—廣東—越南型。

對所有禽流感類型來說，在上述四個地理學參數(shù)中，與其爆發(fā)相關的排序分別為：內(nèi)陸水體>主要鐵路交通線>家禽密度。如果分別考慮各病毒株爆發(fā)點與地理參數(shù)的相關性，則與河北—香港型最密切相關的地理參數(shù)為主要鐵路交通線，與吉林—朝鮮—橫濱以及庫爾干—阿富汗—尼日爾相關最為密切的地理參數(shù)為內(nèi)陸水體。河南—廣東—越南型的禽流感爆發(fā)點分布與上述幾個參數(shù)關系均比較相似。

上述結果所提示的最重要的新信息是：河北-香港型的禽流感病毒株是伴隨著鐵路交通線進行擴散傳播的?？紤]到我國每天有三趟由京廣、京九線運輸鮮活禽類的列車在這條鐵路線上運行，對于監(jiān)測禽流感的變異和控制其爆發(fā)，可能有必要對上述運輸過程加以足夠重視。此外，吉林—朝鮮—橫濱，以及庫爾干—阿富汗—尼日爾型的病毒爆發(fā)點與內(nèi)陸水體分布密切相關。候鳥遷徙的路線，補充淡水的水庫湖泊等內(nèi)陸水體對于候鳥的遷徙路線選擇至關重要。甚至有時候少數(shù)的水體和山脈的地理分布走向會讓大量候鳥在特定季節(jié)聚集于某一個相對比較狹窄的擁有水體的區(qū)域形成“鳥道”^［17］。而吉林—朝鮮—橫濱，以及庫爾干—阿富汗—尼日爾型的H5N1禽流感往往在重要內(nèi)陸水體附近爆發(fā)。在這兩條途徑上，由于大洋、荒漠等天然地理因素的影響，使得候鳥在遷徙過程中往往會聚集于一些特定的內(nèi)陸水體。因此這兩條線路及其周圍的禽流感爆發(fā)和傳播監(jiān)測應重視在候鳥遷徙途中可能聚集的水體。

參考文獻

［1］Webster RG, Bean WJ, Gorman OT, et al. Evolution and ecology of influenza A viruses［J］. Microbiol Rev, 1992,56(1):152-179.

［2］ WHO. Influenza A (H1N1)-update 10. 2009; Available from: http://dse.healthrepository.org/handle/123456789/175.

［3］ Ahn I, Jeong BJ, Bae SE, et al. Genomic analysis of influenza A viruses, including avian flu (H5N1) strains［J］. Eur J Epidemiol, 2006,21(7):511.

［4］ 劉克洲.高致病性禽流感的流行及防治［J］.中華急診醫(yī)學雜志,2004,13(2):143-144.

［5］ 顧雪峰,沈宏韜,邵傳利，等.人高致病性H5N1亞型禽流感一例臨床報告［J］.中華急診醫(yī)學雜志,2007,16(10):1077-1080.

［6］ Janies D, Hill AW, Guralnick R, et al. Genomic analysis and geographic visualization of the spread of avian influenza (H5N1)［J］. Syst Biol, 2007,56(2):321-329.

［7］ Dugan VG, Chen R, Spiro DJ, et al. The evolutionary genetics and emergence of avian influenza viruses in wild birds［J］. PLoS Path, 2008,4(5):e1000076.

［8］ Lacy DB, Tepp W, Cohen AC, et al. Crystal structure of botulinum neurotoxin type A and implications for toxicity［J］. Nat Struct Mol Biol, 1998,5(10):898-902.

［9］ Schuller DJ, Wilks A, de Montellano PRO, et al. Crystal structure of human heme oxygenase-1［J］. Nat Struct Mol Biol, 1999,6(9):860-867.

［10］Sayers E, Wheeler D. Building customized data pipelines using the entrez programming utilities (eUtils)2004［EB/OL］. http://www.ncbi.nlm.nih.gov/books/NBK1058/.

［11］ Butler D. Virtual globes: The web-wide world［J］. Nature, 2006,439(7078):776-778.

［12］ Butler D. Mashups mix data into global service［J］. Nature, 2006,439(7072):6-7.

［13］ Hatta M, Gao P, Halfmann P, et al. Molecular basis for high virulence of Hong Kong H5N1 influenza A viruses［J］. Science, 2001,293(5536):1840.

［14］ Hinshaw V, Wood J, Webster R, et al. Circulation of influenza viruses and paramyxoviruses in waterfowl originating from two different areas of North America［J］. Bull World Health Organ, 1985,63(4):711.

［15］ Kaverin NV, Matrosovich MN, Gambaryan AS, et al. Intergenic HA-NA interactions in influenza A virus: postreassortment substitutions of charged amino acid in the hemagglutinin of different subtypes［J］. Virus Res, 2000,66(2):123-129.

［16］ Li K, Guan Y, Wang J, et al. Genesis of a highly pathogenic and potentially pandemic H5N1 influenza virus in eastern Asia［J］. Nature, 2004,430(6996):209-213.

［17］ Si Y, Skidmore AK, Wang T, et al. Spatio-temporal dynamics of global H5N1 outbreaks match bird migration patterns［J］. Geos Health,2009, 4(1): 65-78.

（收稿日期：2012-04-01）

（本文編輯：何小軍）

DOI：10.3760/cma.j.issn.1671-0282.2012.08.021

基金項目：四川省科技廳（2011SZ0139）及衛(wèi)生廳科研基金支持 (090442；100552；100553)

作者單位： 四川省醫(yī)學科學院 四川省人民醫(yī)院 創(chuàng)傷代謝組多學科實驗室 計算生物學研究組(蔡斌、彭謹、江華、楊浩、孫明偉、Charles Damien Lu、 曾俊); 四川省醫(yī)學科學院 四川省人民醫(yī)院急救中心,急診外科(蔡斌，胡衛(wèi)建); 四川省人民醫(yī)院城東病區(qū)創(chuàng)傷外科(江華、孫明偉)

通信作者：江華, Email：cdjianghua@gmail.com

中華急診醫(yī)學雜志2012年8月第21卷第8期Chin J Emerg Med,August 2012,Vol.21,No.8

P887-891