CN101840467B - 蛋白质组过滤进化分类方法及其系统 - Google Patents

Abstract

本发明公开了一种蛋白质组过滤进化分类方法及其系统，所述方法包括：原始数据获取步骤，获取所要进行分类的物种的基因组数据和蛋白质组数据(S110)；预测步骤，基于所述基因组数据，预测基因组序列中的横向转移基因(S120)；过滤步骤，基于蛋白质组数据和横向转移基因，在蛋白质组中过滤掉与横向转移基因对应的蛋白质序列(S130)；计算步骤，用过滤后的蛋白质组序列计算物种间进化距离(S140)；以及推断步骤，利用物种间的进化距离推断物种间的系统发生关系(S150)。根据本发明，可以排除在推断生物系统发生过程时存在的由横向基因转移所引起的干扰信息，从而可以使推断出的生物系统发生过程更加准确，为生物进化分类提供更加可靠的分类依据。

Description

蛋白质组过滤进化分类方法及其系统

技术领域

本发明涉及一种蛋白质组过滤进化分类方法及其系统。

背景技术

目前在生物进化学研究领域存在有基因组中的横向转移基因预测技术、生物物种间性状距离量化技术和叠加树推断技术。这三类技术目前各有多种，但它们是相对独立的三类技术，在目前的生物系统发育研究以及生物分类学研究中，传统的构建物种间系统发生关系方法都是将后两类技术结合起来，却没有看到第一类技术与后两类技术之间在生物学意义上的内在联系，近年来关于物种基因组间横向基因转移现象的研究表明横向转移基因的存在对生物系统发育关系的重建具有明显的干扰作用，所以目前用传统的方法构建出的物种间的系统发生关系是不准确的，根据这样不准确的系统发生关系建立的生物分类系统不能正确地反映生物进化的历程，本发明正是基于上述考虑，首次将第一类技术与后两类技术结合起来，采用蛋白质组过滤的方法，排除了影响物种间系统发生关系构建的横向转移基因这个因素，完善了传统的系统发育树重建方法，并为生物进化分类研究提供了可靠的分类依据。

发明内容

本发明所要解决的技术问题在于提供了一种通过过滤手段有机地将这三种技术整合起来的蛋白质组过滤进化分类方法及其系统。其实现了将横向基因转移所引起的干扰信息排除，并使生物系统发生关系的推断更加合理，使生物进化分类结果更加完善。

本发明第一方面涉及一种蛋白质组过滤进化分类方法，包括：原始数据获取步骤，获取所要进行分类的物种的基因组数据和蛋白质组数据；预测步骤，基于所述基因组数据，预测基因组序列中的横向转移基因；过滤步骤，基于所述蛋白质组数据和所述横向转移基因，在蛋白质组中过滤掉与所述横向转移基因对应的蛋白质序列；计算步骤，用过滤后的所述蛋白质组序列计算物种间进化距离；推断步骤，利用所述物种间的进化距离推断物种间的系统发生关系。

本发明第二方面涉及一种蛋白质组过滤进化分类系统，包括：原始数据获取模块，用于获取所要进行分类的物种的基因组数据和蛋白质组数据；预测模块，用于基于所述基因组数据，预测基因组序列中的横向转移基因；过滤模块，用于基于所述蛋白质组数据和所述横向转移基因，在蛋白质组中过滤掉与所述横向转移基因对应的蛋白质序列；计算模块，用于使用过滤后的所述蛋白质组序列计算物种间进化距离；推断模块，用于利用所述物种间的进化距离推断物种间的系统发生关系。

本发明的有益效果在于，根据本发明可以实现将横向基因转移所引起的干扰信息排除，并使生物系统发生关系的推断更加合理，使生物进化(生物系统)分类结果更加完善。

应用本发明提供的方法对477个原核生物物种进行系统分类，并将分类结果与国际权威的原核生物分类系统《伯杰氏系统细菌学手册》一一比对，分类效果较传统分类方法明显改善。

附图说明

当结合附图考虑时，通过参照下面的详细描述，能够更完整更好地理解本发明以及容易得知其中许多伴随的优点，但此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定，其中：

图1是蛋白质组过滤进化分类方法的流程图；

图2是过滤步骤的流程图；

图3是蛋白质组过滤进化分类系统2000构成的框图；

图4是过滤模块的框图。

具体实施方式

参见附图对本发明进行说明。

蛋白质组过滤进化分类方法

图1是蛋白质组过滤进化分类方法的流程图。如图1所示，包括以下步骤：原始数据获取步骤S110，获取所要进行分类的物种的基因组数据和蛋白质组数据；预测步骤S120，预测基因组序列中的横向转移基因；过滤步骤S130，在蛋白质组中过滤掉横向转移基因对应的蛋白质序列；计算步骤S140，用过滤后的蛋白质组序列计算物种间进化距离；推断步骤S150，用物种间的进化距离推断物种间的系统发生关系。

在获取步骤S110中，可以通过internet自由免费从多个在线国际生物信息数据库中获取物种基因组序列数据和蛋白质组序列数据，比如可以从美国国家生物技术信息中心的数据库中获取，具体网址为：ftp://ftp.ncbi.nih.gov/genomes/。

预测步骤S120的过程如下：

(1)对基因组中每一个基因(长度小于300bp的基因除外)，计算如下各项，然后分别求各项的均值u和标准差σ：总的GC含量G+C[T]，密码子在第一位上的GC含量G+C[1]、密码子在第二位上的GC含量G+C[2]、密码子在第三位上的GC含量G+C[3]、除终止密码子以外的其他密码子使用频率和氨基酸组成频率，这里的氨基酸是指密码子经翻译后相应的氨基酸。其中，G代表“鸟嘌呤”，是四种核苷酸之一；C代表“胞嘧啶”，是四种核苷酸之一；T代表英文“total”的首字母，即表示“总的”意思。G+C[T]指的是基因序列中G和C的含量占整个基因序列中核苷酸总量的比例，取值范围“0％-100％”，G+C[1]指的是基因序列中三联体密码子第一位上G和C的个数占整个基因序列中三联体密码子总数的百分比，取值范围“0％-100％”，G+C[2]、G+C[3]与G+C[1]类似，只是位置不同。

(2)开一个包含11个基因的“窗口”，扫描整个基因组，寻找所谓的“外源片段”(外源片段的定义详见下文“第三标准”部分)。

(3)将每个基因和基因组都看作具有61个分量的向量(本方法只统计除终止密码子以外的其他61个密码子)，计算每个基因与所在基因组之间的Mahalanobis距离(马氏距离)，计算公式如下：

先计算61×61协方差矩阵S各元素的值：

$S_{i,j}=\underset{K=1}{\overset{N}{\mathrm{\Σ}}}[X_{k,i}-\stackrel{\‾}{X_i}][X_{k,j}-\stackrel{\‾}{X_j}],(i,j=\mathrm{1,2,3},......61)$

S_i，j：协方差矩阵S第i行第j列元素；

N：基因组所含基因的总数；

X_k，i：第i个密码子在第k个基因中的使用频率；

X_k，j：第j个密码子在第k个基因中的使用频率；

每个密码子在整个基因组中的使用频率；

然后计算每个基因与其所在的基因组间的Mahalanobis距离，计算公式如下：

$d^M(X,\stackrel{\‾}{X})={(X-\stackrel{\‾}{X})}^T{S}^{-1}(X-\stackrel{\‾}{X})$

X：61维向量，每个分量分别代表指定基因组中某个基因其各种密码子的使用频率；

61维向量，每个分量分别代表指定基因组其自身各种密码子使用频率；

T：矩阵转置；

S^-1：矩阵求逆。

采用Montecarlo算法，针对每一个具体的基因组，根据其每种密码子的使用频率，随机产生10000条模拟基因序列，采用与上述Mahalanobis距离的计算方法同样的方法，计算每条模拟基因序列与基因组之间的Mahalanobis距离，这样每一个具体的基因组就可获得10000个Mahalanobis距离，求这10000个Mahalanobis距离的均值u′和标准差σ′。具体计算方法如下：

1)首先计算61个密码子分别在某指定基因组中的使用频率(以千分比表示)。

2)计算累计频率值。

比如说，密码子A，B，C，D在某指定基因组中的使用频率(以千分比表示)为：fA＝10‰，fB＝20‰，fC＝30‰，fD＝40‰，则累计频率值分别为：A＝10‰，B＝30‰，C＝60‰，D＝100‰

3)产生一个介于0-1000之间的随机数，根据这个随机数的值介于上述哪两个累计频率值之间，即可确定一个随机密码子，比如，如果随机数的值是2，则可确定A为随机密码子，如果随机数的值是12，则可确定B为随机密码子，如果随机数的值是32，则可确定C为随机密码子，如果随机数的值是62，则可确定D为随机密码子，以此类推。

4)重复上述第三步，直到模拟序列的长度达到500个密码子为止。

5)重复上述第三步及第四步，直到产生10000条模拟序列为止。

6)计算每条模拟序列与基因组间的Mahalanobis距离。

7)计算所有10000个Mahalanobis距离的均值u′和标准差σ′。

8)对每一个基因组皆重复上述第一至第七步，获得每一个基因组各自的均值u′和标准差σ′。

其中，横向转移基因判定标准：

(1)某基因的总的GC含量G+C[T]与所在基因组的总的GC含量的差大于1.5σ，则判定此基因为横向转移基因(第一标准)。

(2)某基因的密码子第一和第三位上的GC含量G+C[1]、G+C[3]同时大于或小于所在基因组的密码子第一和第三位上的GC含量G+C[1]，G+C[3]，并且至少其中一个与所在基因组相应值的差大于1.5σ，则判定此基因为横向转移基因(第二标准)。

(3)如果发现某个包含11个基因的“窗口”中有5个或5个以上基因符合上述第一或第二标准，则称这个“窗口”为一个“外源片段”，在“外源片段”中如果任何一个基因的G+C[T]和“外源片段”的G+C[T]同时大于或小于所在基因组的G+C[T]，则即使不符合所述第一标准和所述第二标准也判定此基因为横向转移基因(第三标准)。

(4)如果某基因与所在基因组之间的Mahalanobis距离与此基因组的u′值之差大于或小于2σ′(此处的u′，σ′均由Montecarlo算法得出)，则判定此基因为横向转移基因(第四标准)。

(5)氨基酸组成上的差异是作为限制性条件使用的，当某个基因符合上述第一标准，或者符合第二标准、或者符合第四标准时，只有在其每种氨基酸频率与整个基因组的相应氨基酸频率之差不大于3σ时才可接受为横向转移基因，但符合第三标准的“外源片段”中的基因不受这个条件的限制(第五标准)。

图2是过滤步骤的流程图，如图2所示，在过滤步骤S130中，包括以下步骤：步骤S132，针对每个待分类的物种读取其蛋白质组序列数据，同时读取预测出的横向转移基因对应的蛋白质版本号；步骤S134，跳过蛋白质组序列数据中蛋白质版本号与预测出的横向转移基因对应的蛋白质版本号相同的所有蛋白质序列，将其余的蛋白质组序列重新写入一个新文件，即可获得过滤后的蛋白质组序列数据。

上述步骤可以利用计算机程序进行，例如可以使用C语言编写一个消息传递并行程序来实现，具体地可以是：

步骤一：开辟N个进程，N等于待分类的物种的个数。

步骤二：每个进程读取一个物种的蛋白质组序列数据，同时读取预测出的横向转移基因对应的蛋白质版本号。

步骤三：每个进程皆执行写文件操作，跳过全蛋白质组序列数据中蛋白质版本号与预测出的横向转移基因对应的蛋白质版本号相同的所有蛋白质序列，将其余的蛋白质组序列重新写入一个新文件。重新写出的文件即是过滤后的蛋白质组。

步骤四：关闭所有开辟的进程退出程序。

计算步骤S140包括以下过程：

给定一条长度为L的氨基酸序列(在具体应用中将过滤后的全蛋白质组序列看做是一条连续的长度为L的氨基酸序列)，计数这条序列中长度固定为K的各种短串的出现频数，氨基酸共有20种，所有短串的种类数为N＝20K。

一条指定长度为L的氨基酸序列，用f(a₁a₂…a_k)来表示K串a₁a₂…a_k出现的频数，其中每一个a_i(i＝1，2，3，......k)都是20种氨基酸中的一种，这个频数用给定的长度为L的氨基酸序列中K串的总数(L-K+1)来除，根据大数定理，商可以作为该短串a₁a₂…a_k在该蛋白质序列中出现的概率p(a₁a₂…a_k)，计算公式如下：

$p(a_1{a}_2...a_k)=\frac{f(a_1{a}_2...a_k)}{(L-K+1)}$

随机突变产生短串a₁a₂…a_k的概率由马尔可夫概率模型来预测得出，详情如下所述：

首先如上所述计算此指定长度为L的氨基酸序列中所有长度为(K-1)和(K-2)的短串的出现概率，那么任意一个长度为K的短串a₁a₂…a_k通过随机突变而产生的概率p⁰(a₁a₂…a_k)都可以通过一个马尔可夫模型来预测，预测公式如下：

$P^0(a_1{a}_2...a_k)=\frac{p(a_1{a}_2...a_{k-1})p(a_2{a}_3...a_k)}{p(a_2{a}_3...a_{k-1})}$

p⁰(a₁a₂…a_k)代表马尔可夫模型预测结果，也就是因随机突变而产生短串a₁a₂…a_k的概率。

计算短串a₁a₂…a_k的修正概率a(a₁a₂…a_k)，计算公式如下：

$a(a_1{a}_2...a_k)=\left\{\begin{array}{cc}\frac{p(a_1{a}_2...a_k)-p^0(a_1{a}_2...a_k)}{p^0(a_1{a}_2...a_k)}& (\mathrm{when}{p}^0(a_1{a}_2...a_k)\≠0)\\ 0& (\mathrm{when}{p}^0(a_1{a}_2...a_k)=0)\end{array}\right.$

以此修正概率来构造物种X的组分矢量A，第i种短串的修正概率a_i(a₁a₂…a_k)作为组分矢量A的第i个元素，其中i的取值范围从1到N＝20^K，把这些元素按照固定顺序排列，就得到了物种X的组分矢量A：

A＝(a₁，a₂，...，a_N)

同样，对于物种Y我们也可以构造相应的组分矢量B：

B＝(b₁，b₂，...，b_N)

物种X和Y之间的关联C(X，Y)可定义为N＝20^K维空间中两个矢量的夹角余弦值，计算公式如下：

$C(X,Y)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{a}_i\×b_i}{{(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{a}_i^2\×\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{b}_i^2)}^{\frac{1}{2}}}$

两个物种之间的距离D(X，Y)定义为：

$D(X,Y)=\frac{1-C(X,Y)}{2}$

所有物种两两之间的距离都被按上述方法确定之后，就形成了用于构建系统发生树的距离矩阵。

推断步骤S150可以采用国际著名系统发生树构建软件包PHYLIP-3.68中的Neighber.exe程序来完成。

蛋白质组过滤进化分类系统

图3是蛋白质组过滤进化分类系统2000构成的框图，如图2所示，蛋白质组过滤进化分类系统2000包括：原始数据获取模块210，获取所要进行分类的物种的基因组数据和蛋白质组数据；预测模块220，预测基因组序列中的横向转移基因；过滤模块230，在蛋白质组中过滤掉横向转移基因对应的蛋白质序列；计算模块240，用过滤后的蛋白质组序列计算物种间进化距离；推断模块250，用物种间的进化距离推断物种间的系统发生关系。

预测模块220的预测过程如下：对基因组中每一个其长度大于等于300bp的基因，计算以下各项：总的GC含量G+C[T]、密码子在第一位上的GC含量G+C[1]、密码子在第二位上的GC含量G+C[2]、密码子在第三位上的GC含量G+C[3]、除终止密码子以外的其他密码子使用频率和氨基酸(密码子经翻译后相应的氨基酸)组成频率，然后求各项的均值u和标准差σ；

开一个包含11个基因的“窗口”，扫描整个基因组，寻找所谓的“外源片段”(外源片段的定义详见下文“第三标准”部分)；

将每个基因和基因组都看作具有61个分量的向量(本方法只统计除终止密码子以外的其他61个密码子)，计算每个基因与所在基因组之间的Mahalanobis距离，所述Mahalanobis距离的计算包括：

先计算61×61协方差矩阵S各元素的值：

$S_{i,j}=\underset{K=1}{\overset{N}{\mathrm{\Σ}}}[X_{k,i}-\stackrel{\‾}{X_i}][X_{k,j}-\stackrel{\‾}{X_j}],(i,j=\mathrm{1,2,3},......61)$

其中，S_i，j：协方差矩阵S第i行第j列元素，N：基因组所含基因的总数，X_k，i：第i个密码子在第k个基因中的使用频率，

第i个密码子在整个基因组中平均使用频率；

X_k，j：第j个密码子在第k个基因中的使用频率，

第j个密码子在整个基因组中平均使用频率；

然后计算所述每个基因与其所在的基因组间的Mahalanobis距离，计算公式如下：

$d^M(X,\stackrel{\‾}{X})={(X-\stackrel{\‾}{X})}^T{S}^{-1}(X-\stackrel{\‾}{X})$

其中，X：61维向量，每个分量分别代表指定基因组中某个基因其各种密码子的使用频率，

61维向量，每个分量分别代表指定基因组其自身各种密码子使用频率，T：矩阵转置，S^-1：矩阵求逆。

采用Montecarlo算法，针对每一个具体的基因组，根据其每种密码子的使用频率，随机产生10000条模拟基因序列，采用与上述计算Mahalanobis距离同样的计算方法，计算每条模拟基因序列与基因组之间的Mahalanobis距离，这样每一个具体的基因组就可获得10000个Mahalanobis距离，求这10000个Mahalanobis距离的均值u′和标准差σ′。

在上述的蛋白质组过滤进化分类系统2000中，所述横向转移基因的判定标准包括：

第一标准：所述某基因总的GC含量G+C[T]与所在基因组总的GC含量的差大于1.5σ，则判定此基因为横向转移基因；

第二标准：某基因的密码子第一和第三位上的GC含量G+C[1]、G+C[3]同时大于或小于所在基因组的密码子第一和第三位上的GC含量G+C[1]、G+C[3]，并且至少其中一个与所在基因组相应值的差大于1.5σ，则判定此基因为横向转移基因；

第三标准：如果发现某个包含11个基因的“窗口”中有5个或5个以上基因符合上述第一或第二标准，则称这个“窗口”为一个“外源片段”，在“外源片段”中如果任何一个基因的G+C[T]和“外源片段”的G+C[T]同时大于或小于所在基因组的G+C[T]，则即使不符合所述第一标准和所述第二标准也判定此基因为横向转移基因；

第四标准：如果某基因与所在基因组之间的Mahalanobis距离与此基因组的u′值之差大于或小于2σ′(此处的u′，σ′均由Montecarlo算法得出)，则判定此基因为横向转移基因；

第五标准：氨基酸组成上的差异是作为限制性条件使用的，当某个基因符合所述第一标准，或符合所述第二标准、或者符合所述第四标准时，只有在其每种氨基酸频率与整个基因组的相应氨基酸频率之差不大于3σ时才可接受为横向转移基因，但是在“外源片段”中的基因不受这个条件的限制。

图4是过滤模块的框图，如图4所示，在上述的蛋白质组过滤进化分类系统2000中，所述过滤模块230包括：读取模块232，针对每个待分类的物种读取其全蛋白质组序列数据，同时读取预测出的横向转移基因对应的蛋白质版本号；获取模块234，跳过所述蛋白质组序列数据中蛋白质版本号与预测出的横向转移基因对应的蛋白质版本号相同的所有蛋白质序列，获得过滤后的蛋白质组。

根据本发明，本发明所涉及的蛋白质组过滤进化分类系统及其方法可用于生物系统发育分析及生物分类，但本方法与其他同类方法相比较，首次增加了过滤由横向转移基因引起的干扰信息的功能，从而使生物系统发生过程的推断更加准确可靠，使生物分类更加合理。

如上所述，对本发明的实施例进行了详细地说明，但是只要实质上没有脱离本发明的发明点及效果可以有很多的变形，这对本领域的技术人员来说是显而易见的。因此，这样的变形例也全部包含在本发明的保护范围之内。

Claims (4)

1.一种蛋白质组过滤进化分类方法，包括：

原始数据获取步骤，获取所要进行分类的物种的基因组数据和蛋白质组数据；

预测步骤，基于所述基因组数据，预测基因组序列中的横向转移基因；

过滤步骤，基于所述蛋白质组数据和所述横向转移基因，在蛋白质组中过滤掉与所述横向转移基因对应的蛋白质序列；

计算步骤，用过滤后的所述蛋白质组序列计算物种间进化距离；以及

推断步骤，利用所述物种间的进化距离推断物种间的系统发生关系，

其中，所述预测步骤包括以下过程：针对每一个基因组中每一个其长度大于等于300bp的基因，计算以下各项：总的GC含量、密码子在第一位上的GC含量、密码子在第二位上的GC含量、密码子在第三位上的GC含量、除终止密码子以外的其他密码子使用频率和氨基酸组成频率，然后求各项指标各自的均值u和标准差σ；

开一个包含11个基因的“窗口”，扫描整个基因组，寻找所谓的“外源片段”；

将每个基因和基因组都看作具有61个分量的向量，计算每个基因与所在基因组之间的Mahalanobis距离，所述Mahalanobis距离的计算方法包括：

先计算61×61协方差矩阵S各元素的值：

$S_{i,j}=\underset{K=1}{\overset{N}{\mathrm{\Σ}}}[X_{k,i}-\stackrel{\‾}{X_i}][X_{k,j}-\stackrel{\‾}{X_j}]$ (i，j＝1，2，3，......61)

其中，S_i，j：协方差矩阵S第i行第j列元素，N：基因组所含基因的总数，X_k，i：第i个密码子在第k个基因中的使用频率，第i个密码子在整个基因组中平均使用频率；

X_k，j：第j个密码子在第k个基因中的使用频率，

第j个密码子在整个基因组中平均使用频率；

然后计算所述每个基因与其所在的基因组间的Mahalanobis距离，计算公式如下：

$d^M(X,\stackrel{\‾}{X})={(X-\stackrel{\‾}{X})}^T{S}^{-1}(X-\stackrel{\‾}{X})$

其中，X：61维向量，每个分量分别代表指定基因组中某个基因其各种密码子的使用频率，

61维向量，每个分量分别代表指定基因组其自身各种密码子使用频率，T：矩阵转置，S^-1：矩阵求逆；

采用Montecarlo算法，针对每一个具体的基因组，根据其每种密码子的使用频率，随机产生10000条模拟基因序列，采用所述Mahalanobis距离的计算方法，计算每条模拟基因序列与基因组之间的Mahalanobis距离，这样每一个具体的基因组就可获得10000个Mahalanobis距离，求这10000个Mahalanobis距离的均值u′和标准差σ′，

其中，判定某基因组中某指定基因是否为横向转移基因的判定标准包括：

第一标准：某基因的总的GC含量与所在基因组的总的GC含量的差大于1.5σ，则判定此基因为横向转移基因；

第二标准：基因的密码子在第一位上的GC含量大于所在基因组的密码子在第一位上的GC含量，并且基因的密码子在第三位上的GC含量大于所在基因组的密码子在第三位上的GC含量，并且至少其中一个与所在基因组相应的值的差大于1.5σ；或者基因的密码子在第一位上的GC含量小于所在基因组的密码子在第一位上的GC含量，并且基因的密码子在第三位上的GC含量小于所在基因组的密码子在第三位上的GC含量，并且至少其中一个与所在基因组相应的值的差大于1.5σ，则判定此基因为横向转移基因；

第三标准：如果发现某个包含11个基因的“窗口”中有5个或5个以上基因符合上述第一或第二标准，则称这个“窗口”为一个“外源片段”，在“外源片段”中如果任何一个基因的总的GC含量和“外源片段”的总的GC含量同时大于或小于所在基因组的总的GC含量，则即使不符合所述第一标准和所述第二标准也判定此基因为横向转移基因；

第四标准：如果某基因与所在基因组之间的Mahalanobis距离与u′的差值大于2σ′，则判定此基因为横向转移基因，其中，此处的u′，σ′均由Montecarlo算法得出；

其中，氨基酸组成上的差异是作为限制性条件使用的，当某个基因符合所述第一标准或所述第二标准、或者符合所述第四标准时，只有在其每种氨基酸频率与整个基因组的相应氨基酸频率之差不大于3σ时才可接受为横向转移基因，但是在“外源片段”中的基因不受这个条件的限制。

2.根据权利要求1所述的蛋白质组过滤进化分类方法，其特征在于，所述过滤步骤包括如下过程：

针对每个待分类的物种读取其全蛋白质组序列数据，同时读取预测出的横向转移基因对应的蛋白质版本号；

跳过所述全蛋白质组序列数据中蛋白质版本号与预测出的横向转移基因对应的蛋白质版本号相同的所有蛋白质序列，重新读取此物种全蛋白质组序列数据，即可获得过滤后的蛋白质组序列数据。

3.一种蛋白质组过滤进化分类系统，其特征在于包括：

原始数据获取模块，用于获取所要进行分类的物种的基因组数据和蛋白质组数据；

预测模块，用于基于所述基因组数据，预测基因组序列中的横向转移基因；

过滤模块，用于基于所述蛋白质组数据和所述横向转移基因，在蛋白质组中过滤掉与所述横向转移基因对应的蛋白质序列；

计算模块，用于使用过滤后的所述蛋白质组序列计算物种间进化距离；

推断模块，用于利用所述物种间的进化距离推断物种间的系统发生关系，

其中，所述预测模块的预测步骤包括以下过程：

对基因组中每一个其长度大于等于300bp的基因，计算以下各项：总的GC含量、密码子在第一位上的GC含量、密码子在第二位上的GC含量、密码子在第三位上的GC含量、除终止密码子以外的其他密码子使用频率和氨基酸组成频率，然后求各项的均值u和标准差σ；

开一个包含11个基因的“窗口”，扫描整个基因组，寻找所谓的“外源片段”；

将每个基因和基因组都看作具有61个分量的向量，计算每个基因与所在基因组之间的Mahalanobis距离，所述Mahalanobis距离的计算方法包括：

先计算61×61协方差矩阵S各元素的值：

$S_{i,j}=\underset{K=1}{\overset{N}{\mathrm{\Σ}}}[X_{k,i}-\stackrel{\‾}{X_i}][X_{k,j}-\stackrel{\‾}{X_j}]$ (i，j＝1，2，3，......61)

其中，S_i，j：协方差矩阵S第i行第j列元素，N：基因组所含基因的总数，X_k，i：第i个密码子在第k个基因中的使用频率，

第i个密码子在整个基因组中平均使用频率；

X_k，j：第j个密码子在第k个基因中的使用频率，

第j个密码子在整个基因组中平均使用频率；

然后计算所述每个基因与其所在的基因组间的Mahalanobis距离，计算公式如下：

$d^M(X,\stackrel{\‾}{X})={(X-\stackrel{\‾}{X})}^T{S}^{-1}(X-\stackrel{\‾}{X})$

其中，X：61维向量，每个分量分别代表指定基因组中某个基因其各种密码子的使用频率，

61维向量，每个分量分别代表指定基因组其自身各种密码子使用频率，T：矩阵转置，S^-1：矩阵求逆；

采用Montecarlo算法，针对每一个具体的基因组，根据其每种密码子的使用频率，随机产生10000条模拟基因序列，采用所述Mahalanobis距离的计算方法，计算每条模拟基因序列与基因组之间的Mahalanobis距离，这样每一个具体的基因组就可获得10000个Mahalanobis距离，求这10000个Mahalanobis距离的均值u′和标准差σ′，

其中，判定某基因组中某指定基因是否为横向转移基因的判定标准包括：

第一标准：某基因的总的GC含量与所在基因组的总的GC含量的差大于1.5σ，则判定此基因为横向转移基因；

第二标准：基因的密码子在第一位上的GC含量大于所在基因组的密码子在第一位上的GC含量，并且密码子在第三位上的GC含量大于所在基因组的密码子在第三位上的GC含量，并且至少其中一个与所在基因组相应的值的差大于1.5σ；或者基因的密码子在第一位上的GC含量小于所在基因组的密码子在第一位上的GC含量，并且基因的密码子在第三位上的GC含量小于所在基因组的密码子在第三位上的GC含量，并且至少其中一个与所在基因组相应的值的差大于1.5σ，则判定此基因为横向转移基因；

第三标准：如果发现某个包含11个基因的“窗口”中有5个或5个以上基因符合上述第一或第二标准，则称这个“窗口”为一个“外源片段”，在“外源片段”中如果任何一个基因的总的GC含量和“外源片段”的总的GC含量同时大于或小于所在基因组的总的GC含量，则即使不符合所述第一标准和所述第二标准也判定此基因为横向转移基因；

第四标准：如果某基因与所在基因组之间的Mahalanobis距离与u′的差值大于2σ′，则判定此基因为横向转移基因，其中，此处的u′，σ′均由Montecarlo算法得出；

其中，氨基酸组成上的差异是作为限制性条件使用的，当某个基因符合所述第一标准或所述第二标准、或者符合所述第四标准时，只有在其每种氨基酸频率与整个基因组的相应氨基酸频率之差不大于3σ时才可接受为横向转移基因，但是在“外源片段”中的基因不受这个条件的限制。

4.根据权利要求3所述的蛋白质组过滤进化分类系统，其特征在于，所述过滤模块包括：

读取模块，针对每个待分类的物种读取其全蛋白质组序列数据，同时读取预测出的横向转移基因对应的蛋白质版本号；

获取模块，跳过所述全蛋白质组序列数据中蛋白质版本号与预测出的横向转移基因对应的蛋白质版本号相同的所有蛋白质序列，重新读取此物种全蛋白质组序列数据，即可获得过滤后的蛋白质组序列数据。

Images