[发明专利]一种用数字标识生物序列和推断物种亲缘关系的方法

权利要求书

1.一种用数字标识生物序列和推断物种亲缘关系的方法，所述用数字标识生物序列和推断物种亲缘关系的方法具体包括：

I信息关联和偏信息关联

以给定DNA序列为例，构成序列的元素是碱基A、G、C、T，根据统计学方法：碱基i出现的概率为p_i，其中i＝A、G、C、T；碱基j出现的概率为p_j，其中j＝A、G、C、T；相距k的碱基对i(k)j出现的联合概率记为p_i(k)j，再根据信息学理论可以得出整条序列的碱基关联信息量：

Dk+2=-2Σipilog2pi+Σijpi(k)jlog2pi(k)j(k=0,1,2...)---(1)]]>

称D_k+2为信息关联；

计算p_i(k)j时，为了避免有限长度N产生的边缘效应，可以引进周期性边界条件，即将序列前面的k+1个碱基接在序列的尾部，形成长为N+k+1的序列，然后将相距k的碱基对i(k)j出现的次数记为N_i(k)j，求出联合概率p_i(k)j，其中加之泰勒展开式，公式(1)可以写成

Dk+2=-2Σipilog2pi+ΣijNi(k)jNlog2Ni(k)jN=12ln2Σijpipj(xi(k)j2-13xi(k)j3+16xi(k)j4-...)]]>

此时x_i(k)j＝(N_i(k)j-Np_ip_j)/Np_ip_j；当N大时，信息关联可以写成

Dk+2≅1ln2Σij(pi(k)j-pipj)2pipj]]>

求和跨越16种碱基关联，为描述特定碱基关联，引入偏信息关联(PIC)

F_i(k)j＝(p_i(k)j-p_ip_j)²(2)

II生物序列识别码

生物序列识别码可以表示成矩阵、向量：

FA(k0)AFA(k0)T...FG(k0)GDk0+2FA(k0+1)AFA(k0+1)T...FG(k0+1)GDk0+3...............FA(k0+d-1)AFA(k0+d-1)T...FG(k0+d-1)GDk0+d+1......(3)]]>

FA(k0)A/Dk0+2FA(k0)T/Dk0+2...FG(k0)G/Dk0+2Dk0+2FA(k0+1)A/Dk0+3FA(k0+1)T/Dk0+3...FG(k0+1)G/Dk0+3Dk0+3...............FA(k0+d-1)A/Dk0+d+1FA(k0+d-1)T/Dk0+d+1...FG(k0+d-1)G/Dk0+d+1Dk0+d+1......(4)]]>

[F_A(k0)A,F_A(k0)T,…,F_G(k0)G,D_k0+2,F_A(k0+1)A,…,F_G(k0+1)G,D_k0+3,…D_k0+d+1]……(5)

[F_A(k0)A/D_k0+2,…,F_G(k0)G/D_k0+2,D_k0+2,F_A(k0+1)A/D_k0+3,…,D_k0+3,…D_k0+d+1]……(6)

无论以何种形式表示，生物序列识别码的核心元素都是信息关联和16种偏信息关联数据，且参数k₀＝0，参数d的范围可以根据需要确定，共有d×17个数据；将生物序列识别码以矩阵形式表示或表述；

III重构生物序列亲缘关系

构建物种进化树，步骤如下：

(1)筛选信息关联

构建形如[D_k0+2,D_k0+3,…,D_k0+d+1]或[F_α(k0)β,F_α(k0+1)β,…,F_α(k0+d-1)β]的矢量X，其中α,β∈{A,G,C,T}；对矢量X的元素做方差分析(ANOVA)和多重比较检验(MCT)；在多重比较时，对于给定物种对，只要矢量X的任意元素通过多重比较检验，就认为该矢量X可以成功区分此物种对；将不可识别的物种对数占总物种对数的比例归一化到100后称为矢量X的失败得分，记为W_X(k₀,d)；失败得分低的矢量X对应的信息参数物种特异性强；

(2)评价亲缘关系

进化距离D需要满足以下三个公理：

(i)D_x,y≥0；当且仅当x＝y时D_x,y＝0；

(ii)D_x,y＝D_y,x；

(iii)对于任意物种x,y和z，D_x,z≤D_x,y+D_y,z恒成立；

可以用马氏距离、欧式距离算法计算碱基关联矩阵之间的距离，并以此作为物种序列进化距离；以欧氏距离公式

Dx,y=Σj=1fΣi=1d(Mx(i,j)-My(i,j))2---(7)]]>

其中M代表碱基关联矩阵，f代表矩阵的列向量；d代表构成列向量的元素；当两两物种之间的距离都被确定之后，就得到了距离矩阵，进而用于绘制进化树；

计算进化距离时使用偏信息关联与信息关联的比值，因为由公式(1)和(2)知，信息关联与偏信息关联的数量级分别为～10^-3和～10^-6，以偏信息关联和信息关联的比值作为参数，可使各个列向量数量级相当；

(3)进化树的统计检验

进化树可使用通用的邻接法(Neighbor-Joining，简称NJ)和算术平均不加权组队法(Unweighted Pair-Group Method withArithmetic，简称UPGMA)来构造系统发生树；对于检验所生成的进化树的鲁棒性(Robustness)，提出一种相当于反向解靴带法或Jack-knife检验的新方法进行检验：建立Z行的碱基关联矩阵，令d为从0到Z-1的可变参数；给定一个d值，可以得到一颗树；通过将生成的树与物种分类信息比较确定最佳d值范围；最后将最佳d值范围内的树整合成一致树，即最终的系统发生树；由此得到的系统发生树的每一分支都会有一个统计值，称为Bootstrap值；d值范围越大说明所得系统发生树越稳定。