基于不同尺度tuple词频的微生物高通量测序数据分析协议

摘要

本发明提供了一种基于不同尺度tuple词频的微生物高通量测序数据分析协议，其包括：步骤1：获取宏基因组样本的2‑10bp的短tuple高通量测序数据，采用插值上下文马尔科夫模型进行建模微生物群落的背景基因组，再采用无监督的聚类方法来比较宏基因组样本，得出宏基因组样本的类别信息；步骤2：基于步骤1)中聚类得出的类别信息，将≥30bp的长tuple作为特征，采用有监督的样本分类方法找出描述宏基因组样本类别的特异性特征长tuple序列。本发明混合不同阶次的马尔科夫模型，由数据本身决定各阶次马尔科夫模型所占的权重，并允许分析上下文不连续的序列之间的关系。

说明书

技术领域

本发明涉及信息技术领域和生物领域，具体来说涉及一种基于不同尺度tuple词频的微生物高通量测序数据分析协议

背景技术

微生物群落是地球上生物多样性最为丰富的资源，广泛存在于各种自然环境中，如土壤、人体皮肤和消化系统中。环境中的微生物传统意义上包括细菌、真菌、病毒以及一些古生菌，不同的环境中微生物的物种丰度组成多样性以及微生物功能多样性都存在着巨大的差异。为了更好地洞悉微生物在不同的微生物环境中所起的功能作用以及更好地理解微生物和环境之间的关系，对环境中所有的微生物基因组研究是极有必要的。

传统的测序方法获取的微生物数量很少，不能从整体上去描述微生物群落之间的结构差异。而高通量测序技术可以获得更完整、更准确的微生物群落结构，因此高通量测序技术逐渐成为研究学者们对微生物群落的比较研究的一个强有力的工具，通过高通量测序技术我们可以直接从环境中获取大量的微生物测序样本，基于这些样本，微生物群落的比较方法大量被提取，其中主要包括基于16S rRNA的方法、基于配准的序列比较方法，如Smith-Waterman算法和Blast算法、基于k-tuple的频度统计方法。然而基于16S rRNA的方法，在微生物群落的分析比较上，存在很大的局限性，能得到的微生物群落构成信息和物种分布范围都很有限。基于高通量测序技术获取的微生物测序数据，许多微生物的基因是未知的，现如今的微生物参考数据库是极其不完备的，而基于配准的方法高度依赖已知数据库或已知基因，这就使得配准的准确性和完整性大大降低。

相比于基于配准的方法，基于无需比对的方法克服了高度依赖参考数据库的不足，为基因间的比较提供了更好的选择。k-tuple方法是最具代表性的无需比对方法，基于k-tuple的频度统计方法在微生物群落分析比较方面主要集中在长度较短的tuple层面上(2-10bp)，结合概率背景统计模型和相异度度量方法，通过非监督的聚类方法，在衡量微生物群落的差异性方面具有优良的性能。然而目前基于这种短k-tuple的方法，只能建立tuple分布的总体统计模型，找出群落与群落之间的关系，度量其整体相异度。但是具体来说是哪些特征序列、哪些微生物/基因序列造成群落间的这种差异和样本类别的分组是k-tuple统计模型无法解决的问题。所以通过无监督的聚类方法在对微生物群落进行比较上就显得不那么完整，而针对无监督聚类获取的样本类别，通过有监督的模式分类能够进一步的识别出不同类别高通量测序数据的特异性tuple，能够为刻画不同类别的微生物群落具体差异以及寻找生物标记物提供重要的参考信息。

众所周知，生物样本是由A、C、G、T四种碱基组成的基因序列，k-tuple是指长度为k的连续字符串序列。所以一个测序样本的k-tuple频度特征向量的维度就是4^k维。之前有研究表明，来自同一个基因组的k-tuple频度相近，能够建立tuple分布的总体统计模型，但不同基因组的k-tuple频度有很大区别。在基于k-tuple频度的无需比对的研究方法中，集中于短tuple(2-10bp)上，且应用于无监督的样本聚类上表现比较优越。因此，基于短k-tuple频度的相异度距离度量方法D₂被提出用来评估比较两个微生物群落样本之间的相异度。此后，在D₂基础上又衍生出和为了更好地应用到高通量测序数据上，通过归一化处理改进的和相继被提出用于比较样本之间的相异度。

用和计算距离时需要结合合适的背景模型来进行建模。在之前的研究中，用到的是定阶次马尔科夫模型和基于变阶次的马尔科夫模型。然而由于微生物群落是由各种各样不同种类的微生物基因组混合组成，很难用几个确定的阶次模拟背景模型，而且需要手动设定模型的阶次，然后去集中对不同阶次模型对聚类结果的优良效果做评估，工作量和计算成本都非常高。对于定阶次马尔科夫模型，阶次越高模型越准确，然而阶次越高，需要的数据量也越多，一般情况下，我们获取的数据量是很难满足需求的。而且基于变阶次的马尔科夫模型在对选择模型阶次时，对其构建的前缀树进行减枝的过程中，需要手动设定阈值，大大提高了模型的不精确性和计算的复杂性。

发明内容

本发明的主要目的在于克服现有技术中的上述缺陷，提出一种通过变尺度tuple频度来区分微生物的群落类别，且基于获取的群落类别，能够找出区分群落类别的特异性信息的基于不同尺度tuple词频的微生物高通量测序数据分析协议。

本发明采用如下技术方案：

基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，包括如下步骤：

步骤1：获取宏基因组样本的2-10bp的短tuple高通量测序数据，采用插值上下文马尔科夫模型进行建模微生物群落的背景基因组，再采用无监督的聚类方法来比较宏基因组样本，得出宏基因组样本的类别信息；

步骤2：基于步骤1)中聚类得出的类别信息，将≥30bp的长tuple作为特征，采用有监督的样本分类方法找出描述宏基因组样本类别的特异性特征长tuple序列。

优选的，所述步骤1)具体包括如下：

步骤1.1:获取宏基因组样本的高通量测序数据，生成该宏基因组样本的短tuple特征频度向量，对每个宏基因组样本中出现的长度为2-10bp的tuple的频度进行统计，并生成相应宏基因组样本的频度向量；

步骤1.2:采用插值上下文马尔科夫模型对微生物群落的背景基因组进行建模，估计频度向量中每一个tuple的马尔科夫概率；

步骤1.3：计算各个宏基因组样本频度向量间的相异度距离，生成一个宏基因组样本间的相异度矩阵；

步骤1.4：根据相异度矩阵生成一个聚类树，用于判断宏基因组样本与样本之间的关系，找出样本的类别信息。

优选的，所述步骤1.1中，tuple特征定义为宏基因组样本中可能出现的字符串组合，并选择长度为2-10bp的字符串组合作为所述tuple特征。

优选的，步骤1.2中，计算tuple的马尔科夫概率具体方法如下：

步骤1.2.1：基于宏基因组样本的频度向量的互信息量构建上下文序列树；

步骤1.2.2：基于上下文序列树计算每个tuple的马尔科夫概率。

优选的，步骤1.2.1具体如下：基于宏基因组样本的频度向量，将k长度tuple中的每一列的字符放在一个向量中，形成A₁,A₂,…,A_k>

$I(A_w,B)=\underset{i}{\Sigma }\underset{j}{\Sigma }P(a_i,b_j)*log\frac{P\left(a_i\right)*P\left(b_j\right)}{P(a_i,b_j)}$

其中,w＝1,2,…,k-1；B＝A_k；a_i,b_j表示向量A,B中的变量；P(a_i,b_j)表示a_i,b_j在对应向量中同时出现的联合概率；P(a_i)表示a_i在对应向量中出现的概率；

找出与B互信息量最大的向量A_w，将此向量对应的下标位置作为上下文序列树的顶点；然后按照该向量中出现的四种不同字符(A,C,G,T)将所有tuple分成四组，；最后对四组tuple向量矩阵，按照公式中的互信息量公式分别计算互信息量，找出每组中与B互信息量最大的向量A_s，s＝1,2,…,w-1,w+1,…,k-1，将此向量对应的下标位置作为上下文序列树的对应叶子的子节点(A,C,G,T)；依次继续下去，直到找到最后一个跟当前状态向量关联性最大的向量，整个上下文序列树构建完毕。

优选的，在步骤1.2.2中，每个tuple的马尔科夫概率公式如下：

P(c₁c₂…c_k)＝P_ICM(c₁)P_ICM(c₂|c₁)…P_ICM(c_k|c₁c₂…c_k-1)

其中，c₁c₂…c_k表示k-tuple序列，P_ICM(c_k|c₁c₂…c_k-1)表示上下文序列c₁c₂…c_k-1转移到当前状态c_k的ICM转移概率。

优选的，对于计算上述每个ICM转移概率，针对所述k-tuple序列c₁c₂…c_k，应用ICM马尔科夫模型构建的上下文序列树中找出与当前状态c_k关联性程度从大到小排序的重要位置，重新构建其上下文序列，具体如下：要构建r阶马尔科夫模型，r≤k-1，从上下文序列树中找出与当前状态c_k关联性程度从大到小排序的重要位置对应的状态分别为c₃,c₄…,c_r，组成插值上下文序列M_r,那么构建ICM的概率模型如下：

P_ICM(c_k|c₁c₂c₃…c_k-1)＝P_ICM,r(c_k|M_r)；

P_ICM,r(c_k|M_r)＝λ_r*P(c_k|M_r)+(1-λ_r)*P_ICM,r-1(c_k|M_r-1)；

$P\left(c_k\right|M_r)=\underset{a\in A,C,G,T}{\Sigma }\frac{N(M_r,A_k=c_k)}{N(M_r,A_k=a)};$

其中，*表示乘积，λ_m表示m阶马尔科夫模型概率所占的权重系数；N(M_r,X_k＝c_k)表示所有的k-tuple中插值上下文序列是M_r，第k个位置是c_k的所有tuple的频度之和，对于上述马尔科夫模型概率所占的权重系数的计算公式为：

其中所述C表示样本阈值，它由赤池信息量准则AIC^R(C)确定，具体公式如下：

AIC^R(C)＝-2λ(S；M_k)+2|M_IMM,k,C|；

其中，λ(S；M_k)表示样本S的伪似然度，计算公式为：

$\lambda (S;M_k)=\underset{c_1,\mathrm{...}{c}_k\in (A,C,G,T)}{\Sigma }N(c_1,\mathrm{...}{c}_{k-1}){\mathrm{logP}}_{ICM}\left(c_k\right|c_1,\mathrm{...}{c}_{k-1});$

|M_IMM,k,C|表示模型自由参数的个数，当AIC^R(C)值最小时算出的C值作为样本的阈值；

所述q表示两个字符串之间差异度的卡方检验值，其计算原理如下：

${\Delta }_r\left(M_r\right)=\underset{a\in A,C,G,T}{\Sigma }\frac{{\left(N\right(M_r,a)-E(M_r,a\left)\right)}^2}{E(M_r,a)};$

E(M_r,a)＝N(M_r)P_ICM，r(a|M_r)；

其中，N(M_r,a)，E(M_r,a)分别表示字符串频度的实际值和理论值，将q＝Δ_r(M_r)作为卡方值，自由度为3，作为卡方检验的指标参数。

优选的，步骤1.3中，应用不同的相异度距离度量方法计算各个宏基因组样本频度向量间的相异度距离，所用到的相异度距离度量方法包括和计算公式如下：

$D_2^S({\tilde{c}}_X,{\tilde{c}}_Y)=\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{X,i}{\tilde{C}}_{Y,i}}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}};$

$d_2^S({\tilde{c}}_X,{\tilde{c}}_Y)=\frac{1}{2}(1-\frac{D_2^S({\tilde{c}}_X,{\tilde{c}}_Y)}{\sqrt{\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{X,i}^2}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}}\sqrt{\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{Y,i}^2}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}}});$

$D_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)=\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{X,i}{\tilde{C}}_{Y,i}}{\sqrt{n_X{p}_{X,i}}\sqrt{n_Y{p}_{Y,i}}}$

$d_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)=\frac{1}{2}(1-\frac{D_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)}{\sqrt{\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{X,i}^2}{n_X{p}_{X,i}}}\sqrt{\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{Y,i}^2}{n_Y{p}_{Y,i}}}});$

其中，和都是计算两个样本之间的相异度的一种距离度量方法；表示样本X的频度向量，表示样本Y的频度向量；称为中心化过程，i＝1,2,…,4^k；C_X,i和C_Y,i分别表示X和Y样本中第i个tuple出现的频度；n_X表示样本X中tuple个数的总和，n_Y表示样本Y中tuple个数的总和；p_X,i和p_Y,i分别表示在插值上下文马尔科夫背景模型下，样本X中第i个tuple的马尔克夫概率和样本Y中第i个tuole的马尔科夫概率；

如果一个数据集中有n个样本，根据相异度距离度量公式计算出的两两样本间的相异度，生成一个n*n维的相异度距离矩阵，这个矩阵定义如下：

N(n,n)＝(d(x,y))_n×n,d(x,y)＝d(y,x),d(x,x)＝0

其中，d(x,y)是两个宏基因组样本的相异度距离，如果不同样本间的距离越小，那么d(x,y)的值就越小；d(x,x)表示相同样本间的距离为0。

优选的，在步骤1.4中，在n*n相异度矩阵的基础上，根据非加权平均法层次聚类算法计算出两个群组的相异度距离，其定义如下：

$d(C_i,C_j)=\frac{1}{\left|C_i\right|·\left|C_j\right|}\underset{x\in C_i}{\Sigma }\underset{y\in C_j}{\Sigma }d(x,y)$

d(x,y)是两个宏基因组样本的相异度距离，|C_i|和|C_j|表示两个群组的大小，也就是群组中样本的个数，i,j＝1,2,…,n；由两两群组的相异度距离得到聚类树，从聚类树中可以直观的看出群落中各个样本之间的结构关系，得出样本之间的类别信息。

优选的，所述步骤2具体包括以下子步骤：

步骤2.1：对样本中出现的长度为40bp的长tuple的频度进行统计，并生成相应样本的频度向量；

步骤2.2：对每个样本的tuple频度向量进行并行处理，生成所有样本的长tuple频度向量矩阵，然后过滤掉冗余的特征；

步骤2.3：基于所述步骤1中获取的样本类别信息，应用过滤好的样本特征对样本进行有监督分类，找到对分类效果具有很强辨识性的特异性tuple特征；

步骤2.4：基于步骤2.3获得的特异性特征，用留一法(LOOCV)来验证和评估分类器的准确率。

优选的，在步骤2.2中，将需要分类样本的tuple频度向量合并在一起，生成一个tuple频度向量矩阵A，A表示为M×N的频度矩阵，其中N表示样本数量，M表示特征维度。

优选的，在步骤2.3中，基于步骤1中获取的样本类别信息，选定训练集和测试集样本，在训练集中选定当前类别和目标类别，然后根据对称不确定性大于设定的阈值时，把冗余的tuple序列特征过滤移除，得到一些类别特异性候选特征，对称不确定性定义如下：

其中，N_X表示当前类别组成的X样本集中tuple特征出现的频度；sum(N_X)表示由当前类别组成的X样本集中特征出现的频度之和；sum(N_Y)表示目标类别组成的Y样本集中特征出现的频度之和；n(X)和n(Y)分别表示X和Y样本集中样本的个数；θ表示X和Y之间对称不确定性的阈值；

采用SVM分类器，对样本进行有监督分类，找到能够描述微生物群落内部具有差异性的特异性特征。

优选的，在步骤2.4中，基于步骤2.3获得的特异性特征，用留一法(LOOCV)来验证和评估分类器的准确率P：

$P=\frac{\underset{i=1}{\overset{\left|D\right|}{\Sigma }}f\left(g\right(x_i),y_i)}{\left|D\right|}$

其中，P表示分类准确率，D是由有限个以(x_i,y_i)形式表示的样本组合集合x_i是样本中除y_i以外的属性列表，y_i表示样本中类别标号的属性，g表示分类器模型函数，输出结果为该模型的预测结果，f(g(x_i),y_i)为判别函数，当g(x_i)与y_i相等时，输出1，否则，输出0。

由上述对本发明的描述可知，与现有技术相比，本发明具有如下有益效果：

1、本发明基于高通量测序数据，通过变尺度tuple频度，不但能够清楚的区分微生物的群落类别，而且基于获取的群落类别，能够找出区分群落类别的特异性信息。

2、本发明使用的方法无需人工选择马尔科夫阶次，能根据数据本身特点自动地选择马尔科夫阶次，而且马尔科夫模型中对应的上下文序列可以是连续的也可以是不连续的；

3、本发明对宏基因组数据的聚类效果明显优于定阶次马尔科夫模型的聚类效果；

4、本发明为了更好更完整的比较分析微生物群落的群落结构和找出微生物群落的类别差异性，针对不同尺度k-tuple词频的微生物高通量测序数据，采用基于无监督聚类和基于聚类得到的样本类别进行有监督分类相结合的微生物群落比较分析方法，把微生物群落的比较分析从统计分布层面扩展到物种和基因分析层面。

5、本发明为了更好地刻画微生物群落间的特异性信息，基于无监督聚类获得的样本类别，我们首次将≥30bp的长tuple作为特征，应用于微生物高通量测序数据比较分析协议的有监督的样本分类方法中，用来找出区分样本类别的特异性tuple序列特征。实例实验表明，当k-tuple长度等于40bp时最能代表两类数据存在的差异性。

附图说明

图1为采用固定阶次马尔科夫模型方法进行聚类的结果；

图2为插值上下文马尔科夫模型方法进行聚类的结果。

具体实施方式

以下通过具体实施方式对本发明作进一步的描述。

本发明提供一种基于不同尺度tuple词频的微生物高通道量测序数据分析协议。基于2-10bp的短tuple高通量测序数据，应用所述插值上下文马尔科夫模型进行建模微生物群落的背景基因组，来比较宏基因组样本，得出宏基因组样本的类别信息。并且找出样本分类的特异性特征，所述方法包括以下步骤：

步骤1：获取宏基因组样本的2-10bp的短tuple高通量测序数据，采用插值上下文马尔科夫模型进行建模微生物群落的背景基因组，再采用无监督的聚类方法来比较宏基因组样本，得出宏基于组样本的类别信息。具体包括如下

步骤1.1:获取宏基因组样本高通量测序数据，生成该宏基因组样本的短tuple特征频度向量，对每个宏基因组样本中出现的长度为2-10bp的tuple的频度进行统计，并生成相应宏基因组样本的频度向量；其中tuple特征定义为为宏基因组样本中可能出现的字符串组合，并选择长度为2-10bp的字符串组合作为所述tuple特征。

步骤1.2:基于插值上下文马尔科夫模型对微生物群落的背景基因组进行建模，估计频度向量中每一个tuple的马尔科夫概率；计算tuple的马尔科夫概率具体方法如下：

步骤1.2.1：基于宏基因组样本的频度向量的互信息量构建上下文序列树；。该步骤中，基于样本的频度向量，将k长度tuple中的每一列的字符放在一个向量中，形成A₁,A₂,…,A_k>

$I(A_w,B)=\underset{i}{\Sigma }\underset{j}{\Sigma }P(a_i,b_j)*log\frac{P\left(a_i\right)*P\left(b_j\right)}{P(a_i,b_j)}---\left(1\right)$

其中,w＝1,2,…,k-1；B＝A_k；a_i,b_j表示向量A,B中的变量；P(a_i,b_j)表示a_i,b_j在对应向量中同时出现的联合概率；P(a_i)表示a_i在对应向量中出现的概率。

找出与B互信息量最大的向量A_w，将此向量对应的下标位置作为上下文序列树的顶点；然后按照该向量中出现的四种不同字符(A,C,G,T)将所有tuple分成四组；最后对四组tuple向量矩阵，按照公式(1)中的互信息量公式分别计算互信息量，找出每组中与Y互信息量最大的向量A_s，其中s＝1,2,…,w-1,w+1,…,k-1，将此向量对应的下标位置作为上下文序列树的对应叶子(A,C,G,T)的子节点；依次继续下去，直到找到最后一个跟当前状态向量关联性最大的向量，整个上下文序列树构建完毕。对于上述构建的上下文序列树，每条枝都对应一个tuple，按照与当前状态关联性程度从大到小，把这个tuple中的碱基位置自上而下排列,储存在这条枝的节点中。

步骤1.2.2：基于上下文序列树计算每个tuple的马尔科夫概率。每个tuple的马尔科夫概率公式如下：

P(c₁c₂…c_k)＝P_ICM(c₁)P_ICM(c₂|c₁)…P_ICM(c_k|c₁c₂…c_k-1)(2)

其中，c₁c₂…c_k表示k-tuple序列，P_ICM(c_k|c₁c₂…c_k-1)表示上下文序列c₁c₂…c_k-1转移到当前状态c_k的ICM转移概率。

对于计算上述每个ICM转移概率，针对上述k-tuple序列c₁c₂…c_k，应用ICM马尔科夫模型构建的上下文序列树中找出与当前状态c_k关联性程度从大到小排序的重要位置，重新构建其上下文序列。例如，要构建r阶马尔科夫模型(r≤k-1)，从上下文序列树中找出与当前状态b_k关联性程度从大到小排序的重要位置对应的状态分别为c₃,c₄,…,c_r，组成插值上下文序列M_r,那么构建ICM的概率模型如下：

P_ICM(c_k|c₁c₂c₃…c_k-1)＝P_ICM,r(c_k|M_r)(3)

P_ICM,r(c_k|M_r)＝λ_r*P(c_k|M_r)+(1-λ_r)*P_ICM，r-1(c_k|M_r-1)(4)

$P\left(c_k\right|M_r)=\underset{a\in A,C,G,T}{\Sigma }\frac{N(M_r,A_k=c_k)}{N(M_r,A_k=a)}---\left(5\right)$

其中，λ_m表示m阶马尔科夫模型概率所占的权重系数；N(M_r,A_k＝c_k)表示所有的k-tuple中插值上下文序列是M_r，第k个位置是c_k的所有tuple的频度之和。对于上述马尔科夫模型概率所占的权重系数的计算公式为：

其中所述C表示样本阈值，它由赤池信息量准则AIC^R(C)确定，具体公式如下：

AIC^R(C)＝-2λ(S；M_k)+2|M_IMM,k,C|(7)

其中，λ(S；M_k)表示样本S的伪似然度，计算公式为：

$\lambda (S;M_k)=\underset{c_1,\mathrm{...}{c}_k\in (A,C,G,T)}{\Sigma }N(c_1,\mathrm{...}{c}_{k-1}){\mathrm{logP}}_{ICM}\left(c_k\right|c_1,\mathrm{...}{c}_{k-1})---\left(8\right)$

|M_IMM,k,C|表示模型自由参数的个数。当AIC^R(C)值最小时算出的C值作为样本的阈值。

所述q表示两个字符串之间差异度的卡方检验值，其计算原理如下：

${\Delta }_r\left(M_r\right)=\underset{a\in A,C,G,T}{\Sigma }\frac{{\left(N\right(M_r,a)-E(M_r,a\left)\right)}^2}{E(M_r,a)}---\left(9\right)$

E(M_r,a)＝N(M_r)P_ICM,r(a|M_r)(10)

其中，N(M_r,a)，E(M_r,a)分别表示字符串频度的实际值和理论值。将q＝Δ_r(M_r)作为卡方值，自由度为3，作为卡方检验的指标参数。

步骤1.3：计算各个样本频度向量间的相异度距离，生成一个宏基因组样本间的相异度矩阵。应用不同的相异度距离度量方法计算各个宏基因组样本频度向量间的相异度距离，所用到的相异度距离度量方法包括和计算公式如下：

$D_2^S({\tilde{c}}_X,{\tilde{c}}_Y)=\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{X,i}{\tilde{C}}_{Y,i}}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}---\left(11\right)$

$d_2^S({\tilde{c}}_X,{\tilde{c}}_Y)=\frac{1}{2}(1-\frac{D_2^S({\tilde{c}}_X,{\tilde{c}}_Y)}{\sqrt{\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{X,i}^2}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}}\sqrt{\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{Y,i}^2}{\sqrt{{\tilde{C}}_{X,i}^2+{\tilde{C}}_{Y,i}^2}}}})---\left(12\right)$

$D_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)=\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{X,i}{\tilde{C}}_{Y,i}}{\sqrt{n_X{p}_{X,i}}\sqrt{n_Y{p}_{Y,i}}}---\left(13\right)$

$d_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)=\frac{1}{2}(1-\frac{D_2^{*}({\tilde{c}}_X,{\tilde{c}}_Y)}{\sqrt{\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{X,i}^2}{n_X{p}_{X,i}}}\sqrt{\underset{i=1}{\overset{4^k}{\Sigma }}\frac{{\tilde{C}}_{Y,i}^2}{n_Y{p}_{Y,i}}}})---\left(14\right)$

其中，和都是计算两个样本之间的相异度的一种距离度量方法；表示样本X的频度向量，表示样本Y的频度向量；称为中心化过程；C_X,i和C_Y,i分别表示X和Y样本中第i个tuple出现的频度；n_X表示样本X中tuple个数的总和，n_Y表示样本Y中tuple个数的总和；p_X,i和p_Y,i分别表示在插值上下文马尔科夫背景模型下，样本X中第i个tuple的马尔克夫概率和样本Y中第i个tuple的马尔科夫概率。

如果一个数据集中有n个样本，根据相异度距离度量公式计算出的两两样本间的相异度，生成一个n*n维的相异度距离矩阵，这个矩阵定义如下：

N(n,n)＝(d(x,y))_n×n,d(x,y)＝d(y,x),d(x,x)＝0(15)

其中，d(x,y)是两个宏基因组样本的相异度距离，如果不同样本间的距离越小，那么d(x,y)的值就越小；d(x,x)表示相同样本间的距离为0。

步骤1.4：根据n*n相异度矩阵生成一个聚类树。由此判断宏基因组样本与样本之间的关系。在相异度矩阵的基础上，根据非加权平均法层次聚类算法计算出两个群组的相异度距离。其定义如下：

$d(C_i,C_j)=\frac{1}{\left|C_i\right|·\left|C_j\right|}\underset{x\in C_i}{\Sigma }\underset{y\in C_j}{\Sigma }d(x,y)---\left(16\right)$

d(x,y)是两个样本的相异度距离，|C_i|和|C_j|表示两个群组的大小，也就是群组中样本的个数，i,j＝1,2,…,n。由两两群组的相异度距离得到聚类树，从聚类树中可以直观的看出群落中各个样本之间的结构关系，得出样本之间的类别信息。

步骤2：基于步骤1)中聚类得出的类别信息，将≥30bp的长tuple作为特征，采用有监督的样本分类方法找出描述样本类别的特异性特征tuple序列。具体包括以下子步骤：

步骤2.1：对宏基因组样本中出现的长度为40bp的长tuple的频度进行统计，并生成相应样本的频度向量。该步骤可参照步骤1.1，只是在统计过程中将tuple的长度变长为40bp。

步骤2.2：对每个宏基因组样本的tuple频度向量进行并行处理，生成所有样本的长tuple频度向量矩阵，然后过滤掉冗余的特征。将需要分类样本的tuple频度向量合并在一起，生成一个tuple频度向量矩阵A，A表示为M×N的频度矩阵，其中N表示样本数量，M表示特征维度。

步骤2.3：基于步骤1中获取的样本类别信息，应用过滤好的样本特征对样本进行有监督分类，找到对分类效果具有很强辨识性的特异性tuple特征。具体的：基于步骤1中获取的样本类别信息，选定训练集和测试集样本，在训练集中选定当前类别和目标类别。然后根据对称不确定性大于设定的阈值时，把冗余的tuple序列特征过滤移除。得到一些类别特异性候选特征。对称不确定性定义如下：

其中，N_X表示当前类别组成的X样本集中tuple特征出现的频度；sum(N_X)表示由当前类别组成的X样本集中特征出现的频度之和；sum(N_Y)表示目标类别组成的Y样本集中特征出现的频度之和；n(X)和n(Y)分别表示X和Y样本集中样本的个数；θ表示X和Y之间对称不确定性的阈值。

采用SVM分类器，对样本进行有监督分类，找到能够描述微生物群落内部具有差异性的特异性特征。

步骤2.4：基于步骤2.3获得的特异性特征，用留一法(LOOCV)来验证和评估分类器的准确率P：

$P=\frac{\underset{i=1}{\overset{\left|D\right|}{\Sigma }}f\left(g\right(x_i),y_i)}{\left|D\right|}---\left(18\right)$

其中，P表示分类准确率，D是由有限个以(x_i,y_i)形式表示的样本组合集合x_i是样本中除y_i以外的属性列表，y_i表示样本中类别标号的属性，g表示分类器模型函数，输出结果为该模型的预测结果，f(g(x_i),y_i)为判别函数，当g(x_i)与y_i相等时，输出1，否则，输出0。

本发明针对由高通量测序获得的宏基因组样本，对微生物群落进行比较和分析。接下来详细描述本发明的方法的实施过程。虽然在以下内容中展示了执行步骤的逻辑过程，但在某些情况下，可以不同此处的顺序执行。

首先执行步骤1中的步骤1.1,获取宏基因组样本的k-tuple频度向量。k-tuple是指长度为k的连续字符串。在本发明中，通过统计这些字符串在样本中出现的频度，并将这些频度组合成一个k-tuple频度向量，以此来代表整个样本的特征。在本发明中，先选择的tuple尺度标准为长度为2-10bp的字符串作为k-tuple的tuple特征。

为了计算tuple序列特征的插值上下文马尔科夫概率，在本实施例中需要执行步骤1.2。在步骤1.2中，首先执行步骤1.2.1：根据样本的所有tuple建立一个上下文序列树，在构建的过程中需根据互信息量最大的准则，依次找到与当前状态关联性最大的点，然后添加到上下文序列的节点中，每个子节点再按A,C,G,T作为叶子向下分枝，在每个分枝下，再依照互信息量最大的准则，向下添加子节点。上下文序列树中，父节点表示的tuple字符位置包含在子节点表示的tuple字符位置中。

在步骤1.2中，当整个上下文序列树构建好之后，这时候每个tuple的上下文序列按照和当前状态的关联性大小从大到小都会依次存储在树的每个节点中。进而进行步骤1.2中的1.2.2子步骤操作。由构建上下文序列树的过程可知，原先tuple的上下文序列转移到下一个状态的概率可以用经过排好序的上下文序列转移到下一个状态的转移概率替代。根据这一原则，可以估算出每一个tuple的马尔科夫概率。

为了计算实施例中k-tuple向量之间的距离，接下来实施步骤1.3。对k-tuple向量分别采取和的相异度方法计算距离。其中距离度量方法中用到的tuple的插值上下文马尔科夫概率，在步骤1.2中已求得。

实施例步骤1.3可以得到一个相异度矩阵，由这个相异度矩阵进行步骤1.4，即进行无监督的层次聚类分析，最终可得到一个聚类树。通过观察聚类树，可以判断聚类情况的好坏，找出样本的类别信息。

实施例步骤2.1和步骤1.1类似。通过实施例中步骤2.2对得到的tuple频度向量进行合并，得到tuple频度向量矩阵。然后实施该步骤中的特征过滤，将样本中的tuple特征出现频度不为零的特征频度归一化为1，然后利用对称不确定性来计算当前类别和目标类别的相关性熵，将相关性熵大于某一设定的阈值的特征留下来，这些特征就是类别特异性候选特征。

利用这些类别特异性候选特征对基于步骤1所获得的类别信息进行有监督的样本分类，需要执行步骤2.3，应用SVM分类器，选定训练集和测试集样本，在训练集中选定当前类别和目标类别，通过学习建立分类模型，找出能分开当前类别和目标类别的特异性tuple特征。为刻画不同类别的微生物群落具体差异以及寻找生物标记物提供重要的参考信息。最后执行步骤2.4，应用留一法对分类器分类准确性进行评估。为刻画不同类别的微生物群落具体差异以及寻找生物标记物提供重要的参考信息。

我们选取24个人体皮肤微生物群落样本(NCBI基因数据库http://www.ncbi.nlm.nih.gov/)进行无监督聚类实验，分别采用了固定阶次马尔科夫模型和插值上下文马尔科夫模型的方法，结果显示样本按人体左右两个不同位置分别进行聚类，而且插值上下文马尔科夫模型的方法(参见图1，)得出的结果要好于固定阶次马尔科夫模型(参见图2，)。

在有监督聚类过程中，选取了99个健康成年人和25个患有肠胃炎(IBD)病人的粪便样本(Qin,J.,et>metagenomic>)，将25个IBD病人样本和25个健康人样本作为训练集，用SVM分类器建立分类模型；将剩下的74个健康人样本作为测试集，并进行LOOCV实验来评估分类器性能，最终结果显示，分别将40-tuple和k-tuple(k＝2-10)作为特征，评估分类器分类性能，基于40-tuple作为特征构建的分类器只要用一个特征就能平均获得100％的准确率，而基于2-10tuple作为特征时，其最好的分类准确率为88％(k＝7)，且需要200个特征。实验表明，长tuple序列比短tuple包含更多的显著的类别特异性信息。一个最直观的表现就是分类性能的准确性。参见表1

表1

上述仅为本发明的具体实施方式，但本发明的设计构思并不局限于此，凡利用此构思对本发明进行非实质性的改动，均应属于侵犯本发明保护范围的行为。

去获取专利，查看全文>