CN101847179A - 通过模型预测流感抗原的方法及应用 - Google Patents

Abstract

本发明公开了一种通过模型预测流感抗原的方法及应用，通过抽提影响流感抗原的12个特征：五个抗原决定簇氨基酸的突变个数、HA蛋白氨基酸的五个理化特性、影响受体结合因素、糖基化位点改变的个数，氨基酸五个理化特性包括疏水性、体积变化、带电性、极性、可积表面积，并对3681对已知的抗原相似的病毒对和1720对抗原变异的病毒对的上述12个特征进行统计，建立一个抗原关系的预测模型，单纯从序列出发，就能给出病毒之间的抗原关系，简单、方便、灵敏度高。通过网络的方式能够形象的展示抗原进化的过程，用于揭示流感传播规律和筛选疫苗候选株等。

Description

通过模型预测流感抗原的方法及应用

技术领域

本发明涉及一种流感抗原预测技术，尤其涉及一种通过模型预测流感抗原的方法及应用。

背景技术

流感病毒是一种全球流行的病毒，它每年感染300～500万人，其中有25～50万人死亡，对人类社会造成巨大危害。流感病毒分A、B、C三个型，A型和B型对人威胁较大，其中A型流感抗原变异频繁，对人类威胁最大。自1968年进入人群后，H3N2亚型流感病毒在人群里占主导地位。H3N2亚型流感病毒基因组包含8个片段，编码11个蛋白，其中HA跟NA是主要的表面抗原蛋白。相对于其它基因，HA变异最快，使得抗原不断发生变化。注射疫苗是目前最有效防治流感的办法之一，由于流感病毒不断发生改变，所以必须不断更新疫苗成分。

目前，世界卫生组织(WHO)通过与全球四个流感参比和研究合作中心及不同国家/地区的流感监测中心紧密合作，实时监测流感的抗原变化情况，并在每年二月(针对北半球)以及九月(针对南半球)通过评估全球流感流行情况推荐下一个流行季使用的疫苗株，指导疫苗的生产。但受人H3N2流感病毒的全球传播规律影响，疫苗株在不同地区的效果存在差异。对于源头地区，例如东亚、东南亚地区，因为新的抗原优势株在这个地区先出现并流行，使得现行推荐疫苗株对这个地区的保护性很差。最理想的情况是各个地区进行抗原监测，针对不同地区的差别分别推荐疫苗株。

目前，使用HI(血凝抑制反应)的方法对抗原进行检测，但这种方法费时、费力，而且有时候不够灵敏。

H3N2病毒无休止地对人体免疫系统的逃避，使得其进化路径沿着一个主干行进，在进化树上表现为一条主干的进化模式，而其主干上的位点变化对其抗原性的进化起主要作用。其中，不同位点对抗原改变的贡献是不尽相同的，有的位点贡献大，而有的位点贡献相对小，但更多时候，抗原的改变是多个位点协同变化的结果。H3N2病毒HA上存在五个抗原表位，是免疫系统抗体的主要识别区域，这些区域上的位点变化对于流感病毒抗原改变有显著的贡献。

现有技术中，先找出跟抗原相关的位点，即所谓正选择位点，然后建立不同的位点模型来模拟和预测抗原变异。这些方法都有一定的预测能力，或多或少抓住了人H3N2流感病毒进化的一些规律。

据目前的研究表明，影响抗原的所谓正选择位点是随时间变化的，即使是同一个位点，其结构背景不同，其变化的效果可能完全不一样。因此，这些基于位点的预测方法的缺点很明显：对应时间段得到的规律只适用于对应时间段的数据，用到其他时间段上效果就会很差。

发明内容

本发明的目的是提供一种简单、方便、灵敏度高的通过模型预测流感抗原的方法及应用。

本发明的目的是通过以下技术方案实现的：

本发明通过模型预测流感抗原的方法，包括通过以下表1和式(1)构建的模型预测病毒对之间抗原相似还是抗原变异：

抽提影响流感抗原的12个特征：五个抗原决定簇氨基酸的突变个数、HA蛋白氨基酸的五个理化特性、影响受体结合因素、糖基化位点改变的个数，所述HA蛋白氨基酸的五个理化特性包括疏水性、体积变化、带电性、极性、可积表面积；

对3681对已知抗原相似的病毒对和1720对抗原变异的病毒对的上述12个特征进行统计，得到：

表1：

表1中的数据分别表示抗原相似的病毒对的特征改变大和小的数量及抗原变异的病毒对的特征改变大和小的数量；

$P_{\mathrm{odds\; ratio}}=\frac{1721}{3682}\×{\left(\frac{3683}{1722}\right)}^{12}\×\underset{j=1}{\overset{12}{\mathrm{\Π}}}\frac{1+\underset{y_{i=1}}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}}{1+\underset{y_i=0}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}},{\tilde{x}}_{\mathrm{ij}}=\left\{\begin{array}{c}{x}_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=1\\ 1-x_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=0\end{array}\right.---\left(1\right),$ 式中：

P_odds ratio表示优胜率，P_odds ratio＜1时被预测的病毒对抗原相似，P_odds ratio＞1时被预测的病毒对抗原变异；

X_j，new表示待预测的病毒对的第j个特征改变情况，X_j，new以及P_odds ratio通过以下方法计算：

首先，分别对所述影响流感抗原的12个特征进行量化，并分别取以下阈值0、2、0、1、0、1.82、54.667、2.493、34.867、0.098、113.607、1；

然后，对被预测病毒对的12个特征分别进行比较，当二者的特征差异小于其对应的阈值时，X_j，new＝0，

对应表1第j个特征在抗原相似的病毒对中特征改变小的数量，

对应表1第j个特征在抗原变异的病毒对中特征改变小的数量；当二者的差异大于或等于其对应的阈值时，X_j，new＝1，

对应表1第j个特征在抗原相似的病毒对中特征改变大的数量，

对应表1第j个特征在抗原变异的病毒对中特征改变大的数量。

本发明的上述通过模型预测流感抗原的方法的应用，用于构建抗原关系网络：

通过预测病毒对之间抗原相似还是抗原变异的关系，建立所有病毒之间的抗原关系网络，把每个病毒作为节点，把抗原相似的病毒之间给一个连线，构成所述抗原关系网络。

由上述本发明提供的技术方案可以看出，本发明所述的通过模型预测流感抗原的方法及应用，通过抽提一些反映抗体抗原相互作用破坏程度的特性，建立一个抗原关系的预测模型，单纯从序列出发，就能给出病毒之间的抗原关系，简单、方便、灵敏度高。通过网络的方式能够形象的展示抗原进化的过程。

附图说明

图1为本发明中模型构建的技术路线示意图；

图2为本发明中抗原关系预测结构模型示意图。

具体实施方式

本发明的通过模型预测流感抗原的方法，其较佳的具体实施方式是，包括：

通过以下表1和式(1)构建的模型预测病毒对之间抗原相似还是抗原变异：

抽提影响流感抗原的12个特征：五个抗原决定簇氨基酸的突变个数、HA蛋白氨基酸的五个理化特性、影响受体结合因素、糖基化位点改变的个数，所述HA蛋白氨基酸的五个理化特性包括疏水性、体积变化、带电性、极性、可积表面积；

对3681对已知抗原相似的病毒对和1720对抗原变异的病毒对的上述12个特征进行统计，得到：

表1：

表1中的数据分别表示抗原相似的病毒对的特征改变大和小的数量及抗原变异的病毒对的特征改变大和小的数量；

$P_{\mathrm{odds\; ratio}}=\frac{1721}{3682}\×{\left(\frac{3683}{1722}\right)}^{12}\×\underset{j=1}{\overset{12}{\mathrm{\Π}}}\frac{1+\underset{y_{i=1}}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}}{1+\underset{y_i=0}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}},{\tilde{x}}_{\mathrm{ij}}=\left\{\begin{array}{c}{x}_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=1\\ 1-x_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=0\end{array}\right.---\left(1\right),$ 式中：

P_odds ratio表示优胜率，P_odds ratio＜1时被预测的病毒对抗原相似，P_odds ratio＞1时被预测的病毒对抗原变异；

X_j，new表示待预测的病毒对的第j个特征改变情况，X_j，new以及P_odds ratio通过以下方法计算：

首先，分别对所述影响流感抗原的12个特征进行量化，并分别取以下阈值0、2、0、1、0、1.82、54.667、2.493、34.867、0.098、113.607、1；

然后，对被预测病毒对的12个特征分别进行比较，当二者的特征差异小于其对应的阈值时，X_j，new＝0，

对应表1第j个特征在抗原相似的病毒对中特征改变小的数量，对应表1第j个特征在抗原变异的病毒对中特征改变小的数量；当二者的差异大于或等于其对应的阈值时，X_j，new＝1，

对应表1第j个特征在抗原相似的病毒对中特征改变大的数量，

对应表1第j个特征在抗原变异的病毒对中特征改变大的数量。

所述的3681对已知的抗原相似的病毒对和1720对抗原变异的病毒对可以通过以下方法得到(也可以通过其它的方法得到)：

已知Smith等人把1968年到2003年间的253株人H3N2流感病毒划分成11个抗原类；

对于这253株病毒，如果两个病毒处于同一抗原类，就认为它们是抗原相似株；如果这两个病毒处于不同的抗原类，就认为他们是抗原变异株，得到31878对两两病毒间的抗原关系；

选取两两病毒HA1蛋白序列差异数介于1-9的病毒对。

所述的阈值通过对所述的3681对已知的抗原相似的病毒对和1720对抗原变异的病毒对进行统计得到。

本发明的上述的通过模型预测流感抗原的方法的应用，用于构建抗原关系网络：

通过预测病毒对之间抗原相似还是抗原变异的关系，建立所有病毒之间的抗原关系网络，把每个病毒作为节点，把抗原相似的病毒之间给一个连线，构成所述抗原关系网络。

还包括对所述抗原关系网络进行聚类。

所述抗原关系网络的聚类包括抽提出所述抗原关系网络中的局部连接密度比较大的区域，作为抗原相似簇；

所述抗原关系网络的聚类可以通过MCL方法(The Markov Cluster Algorithm，马尔可夫聚类算法)，也可以采用其它的方法；

具体可以用于揭示流感传播规律。

还可以用于按以下原则进行疫苗候选株的筛选：

当有新的抗原相似簇出现，并且新的抗原相似簇所占比例不断增加，则选择该抗原相似簇作为疫苗株候选；

如果有多个新抗原相似簇同时满足上面条件，则选择变化更显著的抗原相似簇作为疫苗候选株。

本发明中的模型是通过以下方法得到：

首先，构建训练数据集：

Smith等人，把1968年到2003年间的253株人H3N2流感病毒划分成11个抗原类。通过如下原则构建训练模型需要的训练数据集。对于这253株病毒，如果两个病毒处于同一抗原类，就认为它们是抗原相似株；而如果这两个病毒处于不同的抗原类，就认为他们是抗原变异株，这样可以得到31878对两两病毒间的抗原关系。但这其中包含太多抗原变异病毒对的数据，会影响模型构建，因此选取两两病毒HA1蛋白序列差异数介于1-9的病毒对构建训练数据集，包括3681对抗原相似病毒对，以及1720对抗原变异病毒对。

然后，进行特征选择：

基于流感病毒抗原改变的结构本质抽提了12个特征用于构建抗原关系预测模型。这些特征包括：每个抗原决定簇的氨基酸突变个数(共五个抗原决定簇)，五种氨基酸理化特性(疏水性、体积、带电性、极性、可积表面积)，对受体结合影响，还有就是糖基化位点的改变。

之后，对特征量化，并对特征离散：

给定一个特定的病毒对，就可以通过比较它们的HA1氨基酸序列的差异，计算出上面提到的12个的特性的量化值。

本发明中的模型最终只给出给定的两个病毒抗原是变了还是没变(及两个状态0/1，分别代表抗原相似以及抗原变异)，所有特征也离散成改变大小两个状态(0/1，分别代表特性改变不能导致抗原改变以及能够导致抗原改变)。对于每一个特征，离散的原理就是找到一个阈值，使得以这个阈值为界对训练数据集中的病毒抗原关系对进行划分，划分的结果跟真实的抗原关系匹配最好。通过训练数据集学到的1-12特征的离散化阈值：0、2、0、1、0、1.82、54.667、2.493、34.867、0.098、113.607、1。

最后，进行模型构建：

通过构建12个特征的朴素贝叶斯模型(

Bayes Model)来预测给定病毒对的抗原关系。假定选取的特征满足伯努利模型(Bernoulli Model)，先检验分布满足正态分布，由贝叶斯理论，可以得到给定病毒对的抗原关系P_odds ratio(抗原改变比上抗原不改变的比率)：

$P_{\mathrm{odds\; ratio}}=\frac{\frac{1+\underset{y_i=1}{\mathrm{\Σ}}1}{2+N}}{\frac{1+\underset{y_i=0}{\mathrm{\Σ}}1}{2+N}}\underset{j=1}{\overset{m}{\mathrm{\Π}}}\frac{\frac{1+\underset{y_i=1}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}}{2+\underset{y_i=1}{\mathrm{\Σ}}1}}{\frac{1+\underset{y_i=0}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}}{2+\underset{y_i=0}{\mathrm{\Σ}}1}}=\frac{1+\underset{y_i=1}{\mathrm{\Σ}}1}{1+\underset{y=0}{\mathrm{\Σ}}1}\underset{j=1}{\overset{m}{\mathrm{\Π}}}\frac{\frac{1+\underset{y_i=1}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}}{2+\underset{y_i=1}{\mathrm{\Σ}}1}}{\frac{1+\underset{y_i=0}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}}{2+\underset{y_i=0}{\mathrm{\Σ}}1}},{\tilde{x}}_{\mathrm{ij}}=\left\{\begin{array}{c}{x}_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=1\\ 1-x_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=0\end{array}\right.$

其中y_i表示训练数据集中第i对病毒对的抗原关系(0/1，分别表示抗原相似跟抗原变异)。x_ij表示训练数据集中第i对病毒第j个特性离散值(0/1，分别代表特性改变不能导致抗原改变以及能够导致抗原改变)。m表示我们抽提得到的12个特征(m＝12)。计算针对整个训练数据进行，其实训练数据集给出的就是特性改变跟抗原改变的关系，也即训练集抽提如表1所示。

可以得到

进而得到式(1)；

给定一对病毒，通过12个特征，如果P_odds ratio＞1，抗原变异否则抗原相似。

具体计算实例：

给定一对病毒A/Fujian/411/2002跟A/HongKong/1186/2003；

其HA1的氨基酸差异包括124(S-＞N)，138(A-＞S)，193(S-＞N)，226(V-＞I)，227(S-＞P)，根据每个特征的阈值，可以得到每个特征的改变大小情况，即X_new＝(X_1，new，...，X_12，new)为(1，0，0，0，0，0，0，0，0，0，1，0)

然后根据表1和式(1)计算得到：

$P_{\mathrm{odds\; ratio}}=\frac{1721}{3682}\×{\left(\frac{3683}{1722}\right)}^{12}\×$

$(\frac{1720}{1134}\×\frac{1567}{3664}\×\frac{920}{1686}\×\frac{1464}{2061}\×\frac{991}{1955}\×\frac{35}{771}\×\frac{181}{2341}\×\frac{117}{2469}\×\frac{214}{2465}\×\frac{204}{2122}\×\frac{1556}{1354}\×\frac{1708}{3494});$

$=0.00042565940779$

因为P_odds ratio＜1，因此我们预测A/Fujian/411/2002根A/HongKong/1186/2003抗原相似。

抗原的变化本质上因为位点变化导致抗体抗原的相互作用发生变化，而单纯基于具体位点的模型显然不能反映这个本质。

本发明从结构角度考虑，抽提一些反映抗体抗原相互作用破坏程度的特性，建立一个抗原关系的预测模型，单纯从序列出发，就能给出病毒之间的抗原关系。通过网络的方式能够形象的展示抗原进化的过程。

通过预测抗原关系以及相关分析，可以得到中国大陆人H3N2流感病毒进化的规律，揭示优势抗原由南方到北方的传播规律。还可以更加细致分析亚洲不同地区流感的传播规律，揭示出亚洲热带、亚热带地区的流感传播的源头地位。能够有效监测抗原状态，进行疫苗候选株筛选。如应用到中国大陆地区，筛选出的疫苗候选株能够有效保护这个地区人群。而考虑到人H3N2流感病毒的全球传播趋势，这种基于起源地的抗原监测以及疫苗候选株筛选技术对流感防治意义重大。

下面对本发明的原理和模型构建的过程进行详细的论述：

具体如图1所示，包括：

1、首先以Smith数据建立训练数据集，下载序列并抽提特征及量化，将特征离散化，建立抗原预测模型，验证模型，构建抗原关系网络，最后进行网络聚类。

2、序列数据：

从GenBank下载所有人H3N2流感病毒的HA蛋白质序列，截取其HA1区域，除去较短病毒序列(＜100aa)以及一些特殊序列(相对于所有病毒变化比较显著以及可能年代上标注错误)，共得到7297条病毒蛋白质序列，其中4711条有月份信息，这些病毒序列按照其采集地划分到不同的区域。另外，基于国家流感中心的流感监测网络，收集并测序中国大陆范围的932条病毒，其中506条有月份信息，加上从公共数据库中收集的中国大陆地区的序列数据，共得到中国大陆地区病毒序列1339条，其中705条有月份信息。另外根据病毒分离地点的不同对病毒的爆发区域进行划分，以秦岭淮河一线把中国划分成南北方。

3、训练数据集：

Smith等，对1968年到2003年间的253株人H3N2流感病毒进行抗原测定并最终划分成11个抗原类，代表这段时间内全世界人H3N2流感病毒进化过程。

按照如下的原则从文中抽提训练数据集：

如果两个病毒处于同一抗原类，就认为它们抗原相似；而如果这两个病毒处于不同的抗原类，就认为他们抗原改变。同时，考虑到大于九个位点的氨基酸突变就会产生抗原改变，数据集中包含太多的这类数据会对模型造成影响，因此只保留小于等于九个位点突变的数据。最终从这套数据中抽提得到的训练数据集包含多对病毒的抗原关系，通过比较对应病毒对的HA1序列得到。

4、特征抽提及其量化：

抗原的改变本质上是抗体抗原相互作用的变化，因此按照经验以及数据分析抽提得到四组12个特性来反映抗原关系的改变：

第一组特性包括五个特性，它们分别是流感HA上五个抗原决定簇的位点改变个数。这些特性广泛被人们所接受，主要反映了抗体结合区域的位点变化对抗原的影响；

第二组特性也包含五个特性，这些特性主要从氨基酸变化对物理学直接相互作用破坏程度上来考虑，它们分别是疏水性、体积、带电性、极性和可积表面积；

另外，受体结合区域及其周围区域位点的变化，将影响抗原的变化。这主要是两方面的叠加效应起作用：首先，受体结合区域及其周围区域的位点变化，将导致原来能结合到这个区域的抗体不能再结合，从而产生抗原变化；另外，抗体不能结合也为病毒更有效的结合宿主细胞表面的受体，为病毒的传播提供有利条件，这种适应性的优势将使对应抗原变异株更容易成为抗原优势株。

基于以上分析，把影响受体结合作为第三组特性；

糖基化位点的变化也将影响抗体与流感HA的相互作用，因而也将影响抗原的改变，把其作为第四组特性。

为了表征每个特性对抗原影响的程度，首先对它们进行量化：

第一组的五个特性的变化直接用对应抗原决定簇氨基酸改变个数来表示；

第二组中不同氨基酸物理化学特性的变化，可以从AAindex(氨基酸指数)数据库中抽提。AAindex是一个代表各种理化和氨基酸及氨基酸对生化特性的数值指标数据库，其中需要的特性分别是代表疏水性的FASG890101、代表体积变化的GRAR740103、代表带电性的ZIMJ680104、代表极性的CHAM820101以及代表可积表面变化的JANJ780101；

用距离受体结合区域的远近来度量第三组特性，但这要求首先要确定受体结合区域。应该说，从不同位点到受体结合区域的距离这个角度来讲，人H3N2流感病毒的HA的结构变化不大，因此用早期的H3N2结构(pdb 1MQN)为模板来计算。受体结合区域由三个结构元件组成：135～138的loop(环)、190～198的helix(螺旋)以及221～228的loop。以这三个结构元件出发，结合模板结构，确定131～138、155～160、186～196以及218～228为受体结合区域。任意位点对受体结合的影响用这个位点到受体结合区域的最短距离来表示，但为了体现影响大小与距离大小的关系，用HA上所有位点距离受体结合区域的最远距离减去这个距离来表示，这样距离受体结合区域越近，其影响受体结合的数值就越大。

以上第二组以及第三组特性的计算取两两病毒序列位点变化所导致特性变化的最大的三个值的平均值，之所以这样计算一定程度反映位点变化与抗原变化的一种关联，取平均值避免与第一组特性的重复；

最后一组特性直接用糖基化位点改变个数来计算，糖基化位点的预测用NetNGlyc程序实现，用0.5作为阀值。

5、特征的离散化：

连续变量的过拟和是机器学习中经常遇到的问题，为了避免过拟和，这里对每个特性值进行离散化。前面抽提的训练数据集，其抗原状态已经被离散化：如果训练数据集中包含N对病毒，对于任意一对病毒i(i＝1，...，N)，其抗原关系用y_i表示，如果抗原相似y_i＝0，反之y_i＝1。对于每一个特性j，其量化s_ij的离散化在这里就是找到一个合适的阀值，使得其对抗原关系的区分最好。如果用N₁表示抗原发生改变的病毒对数，用N₀表示抗原未发生改变的病毒对数，则N₀+N₁＝N。对于任意一个特性j，给定阀值c，定义：

$N_{j,11}^{\left(c\right)}=\{i:y_i=1,s_{\mathrm{ij}}>c\},$

$N_{j,10}^{\left(c\right)}=N_1-N_{j,11}^{\left(c\right)},$

$N_{j,01}^{\left(c\right)}=\{i:y_i=0,s_{\mathrm{ij}}>c\},$

$N_{j,00}^{\left(c\right)}=N_0-N_{j,01}^{\left(c\right)}.$

对于特性j最好的阀值通过下面公式得到：

$c_j^{*}=\underset{c}{\arg \max }\underset{p=0}{\overset{1}{\mathrm{\Σ}}}\underset{q=0}{\overset{1}{\mathrm{\Σ}}}\frac{(N_{j,\mathrm{pq}}^{\left(c\right)}-E_{j,\mathrm{pq}}^{\left(c\right)})}{E_{j,\mathrm{pq}}^{\left(c\right)}},$

其中

$E_{j,\mathrm{pq}}^{\left(c\right)}=\frac{N_p(N_{j,0q}^{\left(c\right)}+N_{j,1q}^{\left(c\right)})}{N},p,q\∈\left\{\mathrm{0,1}\right\}.$

理论上，上面的操作是找出N个事例的2×2列联表的最显著卡方检验结果。按照上面方法计算得到的阀值，每个特性将被离散化成0-1，分别表示特性改变没有造成抗原改变以及造成抗原改变，用X表示。

6、抗原预测贝叶斯模型：

如图2所示，

Bayes Model(贝叶斯模型)在统计学习中广泛应用，其基本假设是每个特性之间是独立的。如果用Y表示抗原状态(Y＝0表示抗原相似，Y＝1表示抗原变异，而用X₁，...，X_m(m＝12)表示每一个特性状态，应用Bayes定理：

$P\left(Y\right|X_1,...,X_m)=\frac{P\left(Y\right)\underset{j=1}{\overset{m}{\mathrm{\Π}}}P\left(X_j\right|Y)}{P(X_1,...,X_m)}.$

定义抗原改变的概率比上抗原不改变的概率为优胜率(odds ratio)，它可以通过下面公式计算：

$\frac{P(Y=1|X_1,...,X_m)}{P(Y=0|X_1,...,X_m)}=\frac{P(Y=1)}{P(Y=0)}\underset{j=1}{\overset{m}{\mathrm{\Π}}}\frac{P\left(X_j\right|Y=1)}{P\left(X_j\right|Y=1)}.$

进一步假设Y以及给定Y的每一个X_j满足Bernoulli models，比如：

X_j|p_0j，Y＝0～Bernoulli(p_0j)，

X_j|p_1j，Y＝1～Bernoulli(p_1j)，j＝1，...，m，

Y|p_y～Bernoulli(p_u).

如果认为p_y，p_0j以及p_1j的先验概率为均匀分布，定义训练数据集的抗原状态矢量为y＝(y₁，...，y_N)以及离散化后的特性值矩阵X＝(x_ij)，i＝1，...，N；j＝1，...m，给定训练数据集，p_y，p_0j以及p_1j的后验概率可以很容易计算出来：

$p_{0j}|X,y~\mathrm{Beta}(1+\underset{y_i=0}{\mathrm{\Σ}}{x}_{\mathrm{ij}},1+\underset{y_i=0}{\mathrm{\Σ}}(1-x_{\mathrm{ij}})),$

$p_{1j}|X,y~\mathrm{Beta}(1+\underset{y_i=1}{\mathrm{\Σ}}{x}_{\mathrm{ij}},1+\underset{y_i=1}{\mathrm{\Σ}}(1-x_{\mathrm{ij}})),$

$p_y|X,y~\mathrm{Beta}(1+\underset{y_i=1}{\mathrm{\Σ}}1,1+\underset{y_i=0}{\mathrm{\Σ}}1).$

对于一个给定新的特性的观测量X_new＝(X_1，new，...，X_m，new)，可以得到

$P(X_{j,\mathrm{new}}=1|Y_{\mathrm{new}}=0,X,y)=E\left(P(X_{j,\mathrm{new}}=1|Y_{\mathrm{new}}=0,p_{0j})\right|X,y)$

$=E\left(p_{0j}\right|X,y)=\frac{1+\underset{y_i=0}{\mathrm{\Σ}}{x}_{\mathrm{ij}}}{2+\underset{y_i=0}{\mathrm{\Σ}}1},$

$P(X_{j,\mathrm{new}}=1|Y_{\mathrm{new}}=1,X,y)=E\left(P(X_{j,\mathrm{new}}=1|Y_{\mathrm{new}}=1,p_{1j})\right|X,y)$

$=E\left(p_{1j}\right|X,y)=\frac{1+\underset{y_i=1}{\mathrm{\Σ}}{x}_{\mathrm{ij}}}{2+\underset{y_i=1}{\mathrm{\Σ}}1},$

$P(Y_{\mathrm{new}}=1|X,y)=E\left(P(Y_{\mathrm{new}}=1|p_y)\right|X,y)$

$=E\left(p_y\right|X,y)=\frac{1+\underset{y_i=1}{\mathrm{\Σ}}1}{2+N}.$

而对于给定新的特性的观测量，其预测的优胜率(odds ratio)如下计算：

$\frac{P(Y_{\mathrm{new}}=1|X_{\mathrm{new}},X,y)}{P(Y_{\mathrm{new}}=0|X_{\mathrm{new}},X,y)}=\frac{P(Y_{\mathrm{new}}=1|X,y)}{P(Y_{\mathrm{new}}=0|X,y)}\underset{j=1}{\overset{m}{\mathrm{\Π}}}\frac{P\left(X_{j,\mathrm{new}}\right|Y_{\mathrm{new}}=1,X,y)}{P\left(X_{j,\mathrm{new}}\right|Y_{\mathrm{new}}=0,X,y)}$

$=\frac{1+\underset{y_i=1}{\mathrm{\Σ}}1}{1+\underset{y_i=0}{\mathrm{\Σ}}1}\underset{j=1}{\overset{m}{\mathrm{\Π}}}\frac{2+\underset{y_i=0}{\mathrm{\Σ}}1}{2+\underset{y_i=1}{\mathrm{\Σ}}1}\×\frac{1+\underset{y_i=1}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}}{1+\underset{y_i=0}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}}{\tilde{x}}_{\mathrm{ij}}=\left\{\begin{array}{c}{x}_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=1\\ 1-x_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=0\end{array}\right..$

如果优胜率大于1，认为抗原发生了改变，反之抗原没有发生改变。

7、模型验证：

为了得到上述的抗原关系预测模型对训练数据集本身的预测能力，对训练数据集作10-fold交叉验证。把训练数据集随机分成十份，然后每次留出其中的一份作新的测试数据集，而其余的九份为新的训练数据集，这样重复十次使得每一份都被作为测试数据集被预测一遍，得到预测准确率。同时，为了检验上述的抗原关系预测模型是否反映抗原变化的本质规律，进行前瞻性测试。从Smith文中的数据出发，按照年份信息把数据集分成不同时间段的数据集，分别用时间靠前的数据作训练数据集，而用时间靠后的数据作测试数据集。应该注意，基于每一次用到的新的训练数据集，都要重新进行特性的离散化，重新学习预测模型。

8、构建抗原关系网络：

人H3N2流感病毒的进化过程是一个新抗原替换旧抗原不断反复的过程，本发明能够预测两两病毒的抗原关系，因此可以建立所有病毒之间的抗原关系网络。同时，人H3N2流感病毒的进化表现为从抗原相似的病毒中通过不断突变产生抗原变异株，基于此理解把每个病毒作为节点，而把抗原相似(优胜率odds ratio＜＝1)的病毒之间给一个连线，这样就构成一个抗原相关性网络。因为这个抗原相关性网络展示了人H3N2流感病毒通过不断抗原积累变化的进化过程，可以用这个网络来形象的反映抗原进化。用比较通用的Cytoscape来显示所有网络，并用yFiles Organic Layout来组织网络，这种显示方式能最大程度把网络的模块化给展示出来，反映出人H3N2流感病毒成簇进化的特征。

9、网络聚类：

为了抽提出抗原关系网络中的局部连接密度比较大的区域，也即认为的抗原相似簇，需要对得到的抗原关系网络进行聚类。网络聚类有很多种方法，但基于以下的考虑，选取MCL方法：

MCL方法对每个病毒都进行分类，这符合任何病毒都可以划分抗原状态符合，虽然有一些病毒的抗原状态可能跟其他抗原状态有所区别，但这些抗原状态不是凭空出现，而是由其他一些主要抗原进化而来，从这个角度讲应该对每个病毒都给出抗原状态，这便于分析抗原进化过程。MCL还能够利用边的权重，也就是可以利用病毒与病毒之间的抗原改变的优胜率作为权重来对抗原进行分类。

抗原的进化是一种跟分子进化相比“更不连续”的过程，因此利用好边的权重将对于有效划分抗原类提供帮助。这样，加上权重，对上面得到的抗原相似性网络进行网络聚类，得到不同的网络模块。这些网络模块表现为抗原更相似，就把它们定义为不同抗原类。在分析人H3N2流感病毒整体上抗原进化规律以及传播规律时，我们采用优胜率的负对数作为权重来进行网络聚类得到抗原类；而在疫苗候选株筛选时，因为需要跟踪抗原的细致动态变化，我们用优胜率的倒数作为权重。

本发明可以有以下两方面的应用：

一方面是，揭示流感传播规律：

人H3N2流感病毒的进化是新抗原替换旧抗原的过程，在传播上表现为新抗原优势株由源头向其他区域传播的过程。能够用抗原关系网络来描述抗原的进化过程，并能够通过网络聚类得到不同的抗原类，这样就可以通过分析不同抗原类传播过程来得到传播规律。人H3N2流感病毒存在流行季的概念，即如果按照时间顺序，以月为单位，把每个时间点的病毒数列出来，可以看到在流行季病毒数比较多，而非流行季时间段病毒数相对来说就非常少。因此，考虑到不同流行季所监测病毒数以及测序病毒数的不均匀性，我们采用任意时间点的每个抗原类所包含的病毒数与其前后一个流行季时间段总病毒数的比例来描述病毒爆发以及抗原变化情况。用含有月份信息的数据，通过上面的处理，再比较不同地区新抗原出现的先后就可以确定其传播规律。

中国人H3N2流感病毒进化：

亚洲包括中国在内很可能是全球人H3N2流感病毒优势抗原的起源地，因而加强这个地区的流感抗原监测尤为重要。而很多研究也不约而同的指出，加强对亚洲热带、亚热带地区的流感监测，对全球流感的防治意义重大。中国大陆作为这个地区的大国，在流感的进化中扮演举足轻重的地位。而中国流感中心在中国流感监测中起领导地位，利用其健全的监测网络，可以详细描述人H3N2流感病毒在中国的抗原进化规律，以便更深入有效的理解H3N2在全球的进化与传播。

中国国家流感中心建立了完善有效的人流感监测网络，即广泛分布于全国各个省份直辖市的哨点医院，每天进行采样，并依托分布于全国的流感网络实验室，上报国家流感中心，流感中心对流感的流行情况进行综合的分析与评价。基于这个监测网络，流感中心每年对人H3N2流感病毒的抗原状态进行有效监测，并结合病毒的序列分子进化分析，及时更新和推荐参考株，紧密与WHO沟通推荐疫苗候选株。流感中心的监测网络及其收集的数据对流感区域预防和防治起到十分重要的作用。

为了从整体上有效描述人H3N2流感病毒在中国的进化规律，把WHO历年推荐的疫苗株加到中国大陆的序列库中，通过本发明发展的抗原关系预测模型对两两病毒抗原状态进行预测，并以此建立相应的抗原关系网络，进而对这个网络进行网络聚类，得到不同的抗原类。从进化树上看，中国的人H3N2流感病毒也满足主干式的进化模式，反映在抗原关系网络上，表现为抗原类的不断替换。分析中国流感的进化过程可以看到，中国国家流感中心推荐的参考株以及WHO的疫苗株很好的代表人H3N2流感病毒在中国的抗原进化过程，能够很好的覆盖抗原整个进化过程。

但是，在抗原类WI05(以A/Wisconsin/67/2005疫苗株为代表株)与BR07(以A/Brisbane/10/2007疫苗株为代表株)之间，中国人H3N2流感病毒还存在一个独立的抗原类，而中国流感中心也推荐了相应的参考株用来表征这个抗原类，根据其参考株A/Jiangxidonghu/312/2006，把这个抗原类命名为JX06。JX06在中国只持续了很短的时间(06-07)，并很快被BR07抗原类取代。在中国，JX06却在其流行期间占据H3N2的主导地位，并在抗原相关性网络中形成一个独立一簇。后续的分析表明，JX06在世界其他地方并没有流行，说明JX06在全球范围内并没有成为优势株。人H3N2流感病毒进化树的主干代表优势株的进化过程，如果把抗原流行株的抗原状态描述到进化树的主干上，JX06出现在一个单独的分支上，这可能说明JX06只是一个抗原变异体，不具备很高的适应性，因而没有在全球流行。在人H3N2流感病毒的进化历史上，还存在一个类似的抗原类BE89(以A/Beijing/353/1989疫苗株为代表株)，在进化树上处于一个独立的分支，而在抗原关系图上也表现为游离于其他抗原类之外，但跟JX06不同的是BE89在全球范围内流行过，说明区域变异株的不同命运：如果适应性够好，则会像BE89一样全球流行；而如果适应性不是太好，则只能像JX06一样在局部区域短暂流行，并很快被其他优势抗原替代。JX06的例子同时也说明，人流感抗原进化的复杂性，不同地区可能存在不同，因此对于不同地区，监测了解其抗原状态变化非常重要。

中国人H3N2流感病毒传播规律及其与其他亚洲国家的关系

上面提到了解人H3N2流感病毒区域进化规律对其防治的重要性，那人H3N2流感病毒在中国内部的传播规律如何？有什么特点？中国幅员辽阔，地大物博，地形多变，人口众多，横跨南北的疆域使得其几乎包含所有的气候类型。以秦岭淮河一线，把中国划分成南北方，而南北方在流感的流行规律上却也截然不同：北方属于典型的北半球气候，只有冬季一个流感流行季，在夏天几乎分离不到病毒；而南方则一年有多个流行季，全年流行。

很多文章对不同地区的传播规律有过研究，John Paget et al通过监测临床流感病例的峰值在不同地区的差异来看流感在欧洲的传播规律，而Wladimir J.Alonso et al利用类似的办法来研究流感在巴西的传播，但这些方法都依赖于详细的流行病学监测。而人H3N2流感病毒的进化就是优势抗原不断替代旧抗原的过程，换句话说就是优势抗原不断扩张的过程。既然本发明的方法从给定病毒序列的基础上，就能够有效描述出不同抗原类的进化过程，那本发明就可以直接看出优势抗原的出现在不同区域间有什么特征，从而研究H3N2的传播规律。把序列数据相对较多的2002到2008年这段时间内中国大陆南北方的抗原类进化过程按照月为时间单位描绘出来，可以看到，北方只有一个冬季流行季，而南方在一年则表现出多个流行季：很多时候是春季一个流行季夏季一个流行季。而不同抗原类在所分析的时间段里都表现出由南向北的传播规律，而且一般是优势抗原株在夏季流行季流行，并在随后的冬季流行季传到北方并引起流感流行。

中国南方属于热带、亚热带气候，有很多研究者都指出热带、亚热带地区在流感进化中的特殊地位：因为流感在这些地区全年流行，因而更容易产生优势抗原株。从中国南北方传播的分析中也可以看到，中国南方在中国地区人H3N2流感病毒进化中有着重要地位：优势抗原类都是起源于南方，并在南方先流行，并在随后传到北方。诚然Smith把整个东亚、东南亚地区作为优势抗原的源头，但至少从中国分析的结果可以看出这个地区还是存在差异。

承接对中国传播规律的分析以及亚洲不同地区之间关系的思考，用同样的方法来研究亚洲不同地区间的传播规律。分析可以看到，西亚、亚洲北部(蒙古)以及东亚的韩国、日本有着跟中国北方一致的只在冬季流行季流行的特征，而包括中国香港、澳门以及台湾在内的地区则跟中国南方的流行特点一致：全年流行，一年有多个流感流行季。另外东南亚、南亚地区同样有着跟中国南方一致的流感流行特征。这样可以根据流感的这种流行特征把亚洲地区划分成两个地区：一个是包括中国北方在内的温带地区，这个地区还包括西亚、亚洲北部(蒙古)以及东亚的韩国、日本；另外一个地区是热带、亚热带地区，包括中国南部、中国香港、中国澳门、中国台湾、南亚以及东南亚。而分析抗原类的传播过程可以看到，人H3N2流感病毒在这段时间内都是由亚洲热带、亚热带地区传播到亚洲温带地区，也就是新的优势抗原类一般在亚洲热带、亚热带地区先流行，并在随后的冬季流行季到达亚洲温带地区。而优势抗原株在热带、亚热带地区的传播很复杂，不是由单一的国家或地区向其他地区传播，而是一个复杂的整体，在这个整体内，优势抗原不断出现，并传播到亚洲其他地区。

另一方面是，筛选疫苗候选株：

本发明的方法能有效描述抗原的进化过程，因此本发明可以通过监测抗原状态的变化来及时地筛选疫苗候选株。这里的抗原状态用某一个时间点的对应抗原类所包含的病毒数除以这个时间点所有病毒数的比例来表示，这样能够反映前后抗原比例的变化。中国作为流感新抗原株起源地之一，抗原很多时候都超前，使得现有的疫苗不能很好的保护中国的人群。这里以中国2002年到2008年数据的为例，用本发明的方法来筛选疫苗候选株，看筛选的疫苗候选株的保护效果。考虑到疫苗株制备需要至少6个月的时间，以及中国处于北半球，以对应冬季流行季(10月到来年三月)的抗原状态变化来推荐下一个流行季的疫苗株。考虑到这种流行季相关的推荐方式，以及对应时间段的数据问题，以三个月为单位即季度为单位来分析。这样，以对应冬季流行季前两个季度的抗原状态变化来推荐来年流行季疫苗株，原则是：

1)有新的抗原类出现，并且新的抗原类所占比例不断增加，

2)如果有多个新抗原类同时满足上面条件，则选择变化更显著的抗原类作为疫苗候选株。

基于中国监测数据的疫苗候选株筛选

加深对人H3N2流感病毒进化规律的理解，特别是对其抗原进化规律的把握，了解其传播规律，最终的目的是为了能及时有效的对流感进行防治。现在最有效的防治方式还是疫苗，通过及时准确地推荐疫苗株，就能有效减少流感对于人类造成的伤害。WHO通过全球流感监测网络，对全球流感进行抗原监测，再结合流行病学以及分子进化分析，适时推荐疫苗株。WHO每年分两次分别对南北半球不同的流行季推荐疫苗株，指导疫苗的生产。但因为抗原监测方法的滞后性及灵敏度不够，加之人H3N2流感病毒区域传播所造成的抗原状态不同步问题，使得很多时间、很多地区存在疫苗株与流行株不匹配的问题。因为优势抗原更早的在起源地流行，这个问题在亚洲尤为突出。因此，对于起源地抗原状态的监测，以及基于此的疫苗株推荐就意义重大，这样可以一方面有效保护起源地地区的人群，另外对其他地区疫苗推荐以及流感防治具有十分重要的指导意义。这里以中国数据为例来筛选疫苗候选株。

流感抗原进化就是优势抗原不断替代旧抗原的过程，优势抗原类一旦出现并在人群中流行开，那它将迅速成为流行株，替代原有的旧抗原类。基于这些认识，结合中国处于北半球的事实，从序列出发对疫苗候选株进行筛选：在北半球冬季流行季结束之前，如果有新的抗原类出现，并且其所占比例不断增加，那就使用这个新抗原类合适的病毒株为随后的流行季的疫苗候选株。应用到中国的数据上，分别推荐02-03到07-08流行季代表抗原类FU02、FU02、CA04、WI05、WI05以及BR07的疫苗株，而在中国地区，这段时间流行的优势抗原类分别是FU02、FU02、CA04、WI05、JX06以及BR07，其中03-04流行季流行的抗原类跟FU02抗原上非常相似，可以认为是FU02。在这段时间，只有06-07流行季疫苗株跟流行株不匹配，而综观WHO推荐的疫苗株，则没有一个流行季匹配，而其对北半球其它地区也只有两个流行季能有效保护。这一方面说明现有方式推荐的疫苗株中国以及其他地区的保护不好，同时也说明通过本发明筛选疫苗候选株方法的的合理性与有效性。

以中国的数据出发来筛选疫苗候选株，筛选的疫苗候选株能有效的保护中国地区的人群。但通过研究人H3N2流感病毒的传播规律知道，亚洲热带、亚热带地区作为优势抗原的起源地，优势抗原在这个地区流行比其他地区要更早流行，例如比澳洲、北美以及欧洲要早半个流行季甚至更长时间。这就使得针对优势抗原起源地的抗原监测以及疫苗候选株的筛选具有特殊的意义，因为基于这个地区筛选的疫苗候选株，不仅能有效保护这个地区的人群，对世界其他地区流感防治同样有指导意义。通过本发明的分析可以看到，如果优势抗原在亚洲的冬季流行季起源，那么在随后的夏季流行季传到澳洲(澳洲的冬季流行季)，而下一个冬季流行季到达北美以及欧洲；而如果优势抗原是在亚洲夏季流行季起源，那么在同一个或再下一个夏季流行季(澳洲冬季流行)就有可能传到澳洲，并在接下来的冬季流行季到达北美以及欧洲使得可以有至少半个流行季(3-6个月)的预警期，可以根据起源地分别推荐不同的疫苗株，能对其它地区人H3N2流感病毒进行有效防治。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (8)

1.一种通过模型预测流感抗原的方法，其特征在于，包括通过以下表1和式(1)构建的模型预测两两病毒之间抗原相似还是抗原变异：

抽提影响流感抗原的12个特征：五个抗原决定簇氨基酸的突变个数、HA蛋白氨基酸的五个理化特性、影响受体结合因素、糖基化位点改变的个数，所述HA蛋白氨基酸的五个理化特性包括疏水性、体积变化、带电性、极性、可积表面积；

对3681对已知抗原相似的病毒对和1720对抗原变异的病毒对的上述12个特征进行统计，得到：

表1：

表1中的数据分别表示抗原相似病毒对的特征改变大和小的数量及抗原变异病毒对的特征改变大和小的数量；

$P_{\mathrm{oddsratio}}=\frac{172}{3682}{\left(\frac{3683}{1722}\right)}^{12}\×\underset{j=1}{\overset{12}{\mathrm{\Π}}}\frac{1+\underset{y_i=1}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}}{1+\underset{y_i=0}{\mathrm{\Σ}}{\tilde{x}}_{\mathrm{ij}}},{\tilde{x}}_{\mathrm{ij}}=\left\{\begin{array}{c}{x}_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=1\\ 1-x_{\mathrm{ij}},\mathrm{if}{X}_{j,\mathrm{new}}=0\end{array}\right.---\left(1\right),$ 式中：

P_odds ratio表示优胜率，P_odds ratio＜1时被预测的病毒对抗原相似，P_odds ratio＞1时被预测的病毒对抗原变异；

X_j，new表示待预测的病毒对的第j个特征改变情况，X_j，new以及P_odds ratio通过以下方法计算：

首先，分别对所述影响流感抗原的12个特征进行量化，并分别取以下阈值0、2、0、1、0、1.82、54.667、2.493、34.867、0.098、113.607、1；

然后，对被预测病毒对的12个特征分别进行比较计算，当二者的特征差异小于其对应的阈值时，X_j，new＝0，

对应表1第j个特征在抗原相似的病毒对中特征改变小的数量，

对应表1第j个特征在抗原变异的病毒对中特征改变小的数量；当二者的差异大于或等于其对应的阈值时，X_j，new＝1，

对应表1第j个特征在抗原相似的病毒对中特征改变大的数量，

对应表1第j个特征在抗原变异的病毒对中特征改变大的值。

2.根据权利要求1所述的通过模型预测流感抗原的方法，其特征在于，所述的3681对已知的抗原相似的病毒对和1720对抗原变异的病毒对通过以下方法得到：

已知Smith等人把1968年到2003年间的253株人H3N2流感病毒划分成11个抗原类；

对于这253株病毒，如果两个病毒处于同一抗原类，就认为它们是抗原相似株；如果这两个病毒处于不同的抗原类，就认为他们是抗原变异株，得到31878对两两病毒间的抗原关系；

选取两两病毒HA1蛋白序列差异数介于1-9的病毒对。

3.根据权利要求2所述的通过模型预测流感抗原的方法，其特征在于，所述的阈值通过对所述的3681对已知的抗原相似的病毒对和1720对抗原变异的病毒对进行统计得到。

4.一种权利要求1、2或3所述的通过模型预测流感抗原的方法的应用，其特征在于，用于构建抗原关系网络：

通过预测病毒对之间抗原相似还是抗原变异的关系，建立所有病毒之间的抗原关系网络，把每个病毒作为节点，把抗原相似的病毒之间给一个连线，构成所述抗原关系网络。

5.根据权利要求4所述的通过模型预测流感抗原的方法的应用，其特征在于，包括对所述抗原关系网络进行聚类。

所述抗原关系网络的聚类包括抽提出所述抗原关系网络中的局部连接密度比较大的区域，作为抗原相似簇。

6.根据权利要求5所述的通过模型预测流感抗原的方法的应用，其特征在于，所述抗原关系网络的聚类通过MCL方法。

7.根据权利要求6所述的通过模型预测流感抗原的方法的应用，其特征在于，用于揭示流感的传播规律。

8.根据权利要求6所述的通过模型预测流感抗原的方法的应用，其特征在于，用于按以下原则进行筛选疫苗候选株：

当有新的抗原相似簇出现，并且新的抗原相似簇所占比例不断增加，则选择该抗原相似簇作为疫苗候选株；

如果有多个新抗原相似簇同时满足上面条件，则选择变化更显著的抗原相似簇作为疫苗候选株。

Images