基于基因进化树和地理数据库追踪禽流感病毒变异

蔡斌 彭谨 江华 杨浩 孙明伟 CharlesDamienLu 胡卫建 曾俊

【摘要】目的 禽流感疫情的爆发和传播受到多种自然因素的影响。今欲尝试将地理信息系统与基因进化树分析相结合，以建立一种基于基因序列变异追踪中国禽流感病毒地理传播的技术。方法 禽流感病毒基因来源于美国国立医学图书馆（National Library for Medicine, NLM）数据库，所获得的基因组数据利用E-Utilities软件包转化为结构体后，可用Matlab软件阅读。结构体主要字段包括PB2、PB1、PA、HA、NP、HA、M1和NS1 8个片段，分别代表流感病毒的8个不同的基因片段。基于结构体字段，利用计算生物学的方法比较不同传播能力禽流感病毒的同义突变/非同义突变基因（Ka/Ks）比例，确定不同选择压力之下A型禽流感病毒的基因突变模式。进而选择Ka/Ks比例最大的基因片段，采用Jukes-Cantor算法估计氨基酸序列变异的进化距离，然后对不同爆发点的H5N1型禽流感进行进化树聚类。将聚类信息输入Google Earth，并利用不同图层地理信息对影响爆发点分布的因素做单因素分析。结果 比较分析A型禽流感所有的8个基因序列可以看出，NS1、HA和NA蛋白的Ka/Ks比值较大。三者中，HA基因的Ka/Ks比值最大，可以代表病毒的传播能力。利用分级聚类的思路对HA基因转录的氨基酸相似程度进行比较, 发现自2003年以来亚洲地区爆发的H5N1型禽流感之间的关系可以表示为一个由30个节点构成的进化树，其中14个节点为分支节点，16个节点为叶子结点。把分支树的前三个节点作为分类标准，可以把所有16个病毒株分为四类。这四类病毒在地理空间的分布呈现一定规律。计算发现禽流感爆发相关地理因素排序分别为：内陆水体>主要铁路交通线>家禽密度。结论 对中国H5N1病毒株基因序列变异的地理分布分析显示，禽流感病毒爆发与候鸟迁徙、家禽运输密切相关。

【关键词】禽流感；病毒基因变异；谷歌地图；地理信息系统；基因进化树；同义突变/非同义突变基因比例；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流感）是由流感病毒引起的急性呼吸道传染病。该病病毒易发生抗原性变异，加之人群对变异株普遍易感，常迅速导致世界性的流感大流行^{［1- 2］}。在流感病毒的发病过程中，禽流感（avian flu）由于其在世界各国此起彼伏的暴发和流行，及可能与人流感病毒进行基因重组，成为研究的热点。虽然迄今并没有高致病性禽流感病毒能在人与人之间直接传播的确切证据，但是由于两种病毒的基因相似性使得人们担心这种可能性正在逐渐增大。 流感病毒基因组主要由8个RNA片段组成。这8个片段为：聚合酶B2（polymerase B2, PB2）、聚合酶B1（polymerase B1）、聚合酶A（PA）、血凝素A（HA）、神经氨酸酶（neuraminidase, NA）核蛋白（nucleoprotein, NP）、基质（M1）以及非结构蛋白（non-structural, NS1）^［3-5］。其中病毒的高度变异能力主要依赖于其病毒颗粒最外层的两种表面抗原：血凝素抗原和神经氨酸酶抗原^［6-7］。在致病力最强的A型禽流感当中，编码HA和NA的基因由于氨基酸序列的不同又被分为HA（H1-16）以及NA（N1-9）几种蛋白质亚型^{［8- 9］}。

基于进化生物学视角，禽流感变异对于流行病学的意义在于：（1）病毒基因的连续变异可能使得这种病毒获得新宿主；（2）病毒基因的保守序列用以追踪禽流感病毒来源。例如，WHO建立的Influza数据库每年都在发布新的禽流感病毒基因序列变异测序。考虑到这些变异都是发生在一定的时间和空间当中，候鸟迁徙、人类家禽运输等活动都有可能对禽流感的传播和变异起到重要的作用。这些活动很可能使得相隔很远的两个不同的地理位置出现相同基因型禽流感的爆发。另一方面，我们也可以观察到地理上相邻位置的流感病毒基因型可能出现明显差异。由此，以下两类问题——（1）“哪些地理空间中禽流感更易變异？”；（2）“不同基因型的禽流感的是否通过类似的途径进行传播？”——对于研究流感病毒的传播特性和变异规律至为重要。要获得这样的知识，需要引入新的、能够同时在基因变异和地理空间分布上对禽流感的情况进行追踪和标记的研究工具。

1 资料与方法

1.1 研究的基本流程

研究的基本流程如下。（1）筛选出能够正确标记流感病毒抗原变异的基因：笔者首先对比不同流感病毒基因片段的同义突变和非同义突变（Ka/Ks）比值，丛中筛选出非同义突变明显偏高的基因片段，将这些基因片段的突变作为流感病毒获得新的传播能力的标记；（2）基因序列聚类以定义变异程度：利用进化树、主成分分析等方法对这些基因片段进行聚类；（3）把不同聚类的病毒类型通过google earth API数据库叠加到地图上，与我们感兴趣的地理分布特征进行相关性分析，以探索影响禽流感爆发流行的自然环境和地理相关因素。

1.2 禽流感病毒基因序列数据获取与整理

禽流感病毒基因来源于美国国立医学图书馆（National Library for Medicine, NLM）数据库。下载选择参数为：病毒分型=A，宿主=Avian, 国家/地区=Canada & Hong Kong，病毒片段=ALL，病毒基因表型=H5N1 & H2N3。下载1997和2001年两年的数据。另外，还根据基因变异度分析的结果，选取病毒分型=A，宿主=Avian,国家和地区=Asia & African病毒片段=HA，年份=2001-2007年。将所获的基因组数据利用美国国立卫生研究院（National Institute for Health, NIH）开发的E-Utilities软件包转化为可由Matlab （版本7.10.0.499，Mathwork, USA）阅读的结构体^［10］。结构体主要字段包括PB2、PB1、PA、HA、NP、HA、M1和NS1，分别代表流感病毒的8个不同的RNA片段。

1.3 计算Ka/Ks比值

其比较思路是首先对比所下载的基因组数据中的氨基酸序列和核酸序列，确定核酸序列密码子序列。再利用上述密码子序列分别计算同义突变和非同义突变的数量。利用非同义突变的发生概率确定选择压力，同时估算多重突变对结果进行校正（即：多次发生突变重新变回原来的密码子）的概率。最终获得各段基因的Ka/Ks比值。

1.4 HA基因进化树和相似度聚类

Ka/Ks分析发现HA基因变异最能表征病毒的传播能力。据此在上述数据库中选取2001年以来，亚洲和非洲H5N1流感的HA基因片段测序数据。通过测量该基因表达的氨基酸序列变异程度，建立各段氨基酸序列相似程度的度量单位。利用Jukes-Cantor方法构建系统发育树。获得各地不同时间H1N1禽流感亲缘关系数据图。并且利用分级聚类的原理将禽流感分为4类。具体算法为：

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

d为两个氨基酸序列之间的距离， P为两个序列氨基酸的相似程度。P=1意味着两氨基酸序列完全相同，P=0意味着两氨基酸序列构完全不同。

1.5 Google earth 地理信息系统叠加

Google Earth地理信息系统（http://www.google.com/earth/index.html）是一个可以通过API函数进行图层叠加开发的开放式平台，2005年Declan Butler为google earth增加了一个能够标定流感基因分型和爆发人数的独立图层（人及禽流感）^［11-12］。在这个图层基础上，利用H1N1的禽流感分类信息叠加上人口聚集区、铁路交通、养鸡场分布、以及主要内陆水体四条基本信息。为每一个爆发点构建一个向量，如果爆发点周围存在上述任何一种地理单元，则为向量中的对应单元数据赋值1；反之则为0。采用单因素分析对禽流感爆发点与潜在传播影响因素之间的相关程度进行估计。

2 结果

2.1 追寻禽流感病毒基因的变异程度

Ka/Ks比值分析是一种有效跟踪和比较病毒基因组中有效和无效突变的技术。同义突变增加表明生物选择压力增大，但是同义突变本身并不会改变翻译蛋白质的氨基酸序列。非同义突变则表明蛋白质氨基酸序列会发生改变，病毒可能获得新的抗原特征从而造成爆发流行。本研究中，比较分析A型禽流感所有的8个基因序列可以看出，NS1、HA和NA蛋白的Ka/Ks比值最大（图1）。这说明该三个基因所受到的选择压力最大，是决定病毒流行的主要因素,该发现与已有研究一致。而三者中的HA基因所发生的改变可能使病毒具有了在除鸟类以外的其他物种中传播的能力。因此，HA基因的变异可作为表征病毒传播能力的生物学变化。

2.2HA蛋白的系统发育学分析

利用Jukes-Cantor算法估计氨基酸序列之间的距离，获得16个基因两两间的欧氏距离，总共可以获得120个距离。利用分级聚类的思想对氨基酸之间的相似程度进行比较。可以发现自2003年以来亚洲地区爆发的H5N1型禽流感之间的关系可以表示为一个由30个节点构成的进化树，其中14个节点为分支节点，16个节点为叶子结点。把分支树的前三个节点作为分类标准，可以把所有16种基因分为4类。类似的，利用氨基酸序列相关性的120个距离也可以在三维空间中获得16个基因的位置关系（图2）。

2.3 地图数据库功能叠加

利用Google Earth选取高密度居民点、内陆水体、主要铁路交通线、养鸡场密度等四个指标在地图上进行标记，同时对进化树上有亲缘关系的禽流感病毒株与上述区域的重叠度进行相关性分析。结果发现各HA基因组分型禽流感病毒株与地理空间分布的相关性见图3。

3 讨论

禽流感爆发受到包括候鸟迁徙、家禽养殖和运输等多种自然与人类活动因素的影响^［4-5］。如何对这些影响进行定量分析，一直以来都是禽流感流行病学研究中的难点。传统的方法很难对相同血清型

的禽流感病毒进行进一步分类，也难于对病毒在地理空间中的传播和变异进行追踪。上述难题的解决，需要引入新的研究思路。现代计算生物学技术在近十年取得的长足发展，已经使得对禽流感病毒基因型的细微变异进行精确统计分析和聚类成为可能。同时，地理信息数据库系统及其可视化技术的发展，使得有可能将病毒的空间分布和传播信息与其生物学信息相结合。本研究正是在这一思路下，就禽流感在中国的传播规律，建立了首个基于计算生物学-地理信息系统方法的模型。

1997年和2001年，甲型流感病毒的亚型（H5N1）在香港爆发，同其他禽流感不一样的是，这次爆发的流感病毒具备了初步的禽-人的传染能力，并导致了6人死亡。香港特区政府对病毒进行了全基因组测序，并向世界卫生组织提供了病毒基因组全序列^［13］。为了比较和跟踪禽流感的变异情况，笔者利用另外一种变异相对稳定的禽流感病毒：加拿大卫生部门保存的1985年和1977年的阿尔伯塔鸭（A/H2N3，Alberta duck）禽流感病毒^［14］来进行对照。通过对比两种禽流感的基因序列的变异情况，探寻H5N1更加容易感染人类的原因。

已有研究发现，禽流感的7个基因在流感传播过程中具有不同的功能。HA蛋白可以帮助病毒与宿主细胞粘附在一起，并且进入到宿主细胞内部，NA蛋白可以剪切新生成的病毒，帮助病毒从宿主细胞中逸出进入新的健康细胞^［15］。不同的HA/NA蛋白组合对于流感病毒对不同物种呼吸道上皮细胞的侵袭性有不同影响。病毒学中，一个已知的重要事实是：H5N1鸡型禽流感和N2N3鸭型禽流感之间最重要的差异，在于前者已经初步具有跨种属间的传播能力^［16］。本研究从分子进化的角度印证了上述结论：传播力最强的H5N1鸡禽流感的每一段基因的Ka/Ks比值都较传播力弱的H2N3野鸭禽流感要大，说明H5N1处在一个选择压力更高的环境。对比两种禽流感病毒全部基因的Ka/Ks比值，可以发现，HA蛋白最大。这说明HA蛋白的变异最能说明两种病毒之间的差异。所以，无论是就蛋白质的功能还是基因的变异程度而言，HA蛋白序列的差异都是追踪不同病毒株之间基因差异的最好标记。

为了追踪不同时间和地点有亲缘关系的病毒株之间的关系。笔者将亚洲和非洲不同地区不同时期从鸡身上分离的H5N1禽流感病毒的HA蛋白进行系统发育学的分析，分析它们之间相互的关系，并且利用人口聚居点、内陆水体、主要铁路交通线、养鸡场的地理参数与禽流感是否爆发进行相关性分析。通过系统发育学比较，可以把H5N1型禽流感病毒株分为以下四种地理类型：1）河北—香港型，2）吉林—朝鲜—横滨型，3）库尔干（西伯利亚）—阿富汗—尼日尔型，以及4）河南—广东—越南型。

对所有禽流感类型来说，在上述四个地理学参数中，与其爆发相关的排序分别为：内陆水体>主要铁路交通线>家禽密度。如果分别考虑各病毒株爆发点与地理参数的相关性，则与河北—香港型最密切相关的地理参数为主要铁路交通线，与吉林—朝鲜—横滨以及库尔干—阿富汗—尼日尔相关最为密切的地理参数为内陆水体。河南—广东—越南型的禽流感爆发点分布与上述几个参数关系均比较相似。

上述结果所提示的最重要的新信息是：河北-香港型的禽流感病毒株是伴随着铁路交通线进行扩散传播的。考虑到我国每天有三趟由京广、京九线运输鲜活禽类的列车在这条铁路线上运行，对于监测禽流感的变异和控制其爆发，可能有必要对上述运输过程加以足够重视。此外，吉林—朝鲜—横滨，以及库尔干—阿富汗—尼日尔型的病毒爆发点与内陆水体分布密切相关。候鸟迁徙的路线，补充淡水的水库湖泊等内陆水体对于候鸟的迁徙路线选择至关重要。甚至有时候少数的水體和山脉的地理分布走向会让大量候鸟在特定季节聚集于某一个相对比较狭窄的拥有水体的区域形成“鸟道”^［17］。而吉林—朝鲜—横滨，以及库尔干—阿富汗—尼日尔型的H5N1禽流感往往在重要内陆水体附近爆发。在这两条途径上，由于大洋、荒漠等天然地理因素的影响，使得候鸟在迁徙过程中往往会聚集于一些特定的内陆水体。因此这两条线路及其周围的禽流感爆发和传播监测应重视在候鸟迁徙途中可能聚集的水体。

参考文献

［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］.中华急诊医学杂志,2004,13(2):143-144.

［5］ 顾雪峰,沈宏韬,邵传利，等.人高致病性H5N1亚型禽流感一例临床报告［J］.中华急诊医学杂志,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）及卫生厅科研基金支持 (090442；100552；100553)

作者单位： 四川省医学科学院 四川省人民医院 创伤代谢组多学科实验室 计算生物学研究组(蔡斌、彭谨、江华、杨浩、孙明伟、Charles Damien Lu、 曾俊); 四川省医学科学院 四川省人民医院急救中心,急诊外科(蔡斌，胡卫建); 四川省人民医院城东病区创伤外科(江华、孙明伟)

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

中华急诊医学杂志2012年8月第21卷第8期Chin J Emerg Med,August 2012,Vol.21,No.8

P887-891