一种宏基因组和宏转录组样本相异度的比较方法

摘要

一种宏基因组和宏转录组样本相异度的比较方法，涉及信息和生物技术。生成样本的tuple频度向量，对样本中出现的长度为1～10的tuple的频度进行统计，并生成相应样本的频度向量；计算tuple的马尔克夫概率，基于变阶次马尔克夫模型估计频度向量中每一个tuple的马尔克夫概率；生成样本间相异度矩阵，计算各个样本频度向量间的距离，生成一个样本间的相异度矩阵；生成聚类树，根据相异度矩阵生成一个聚类树。无需人工选择马尔克夫阶次，能根据数据特效自动地选择马尔克夫阶次；对宏基因组和宏转录组数据的聚类效果明显优于定阶次马尔克夫模型的聚类效果。

说明书

技术领域

本发明涉及信息和生物技术，尤其是涉及一种宏基因组和宏转录组样本相异度的比较方 法。

背景技术

微生物群落间的比较对于理解微生物和环境之间的关系至关重要。高通量测序技术已经 成为表征微生物群落的一个强有力的工具。对于不同基因间的比较，基于配准的序列比较方 法，如Smith-Waterman算法和Blast算法已经被广泛应用。然而对于高通量测序数据，基于 配准的方法变得不再适用，主要由于以下原因：首先，基于配准的方法高度依赖已知数据库 或已知基因，然而许多微生物的基因是未知的，这就影响了配准的准确性。其次，基于配准 的方法要对短序列进行组装，这项工程太耗时。因此，免配准的方法为基因间的比较提供了 更好的选择。k-tuple方法是一个经典的免配准方法。生物样本是由A、C、G、T四种碱基组 成的序列，因此可以看成是由A、C、G、T四种字符组成的文本序列。k-tuple是指长度为k 的连续字符串。之前的研究表明，来自同一个基因组的k-tuple频度相近，但不同基因组的 k-tuple频度有很大区别。因此，基于k-tuple频度的相异度方法D₂被提出用来评估比较两 个生物样本之间的距离。此后，在D₂基础上改进的和相继被提出用于比较样本之 间的距离。

用和计算距离时需要用到一个合适的背景模型。在之前的研究中，用到的是定 阶次马尔克夫模型。然而由于微生物群落是各种基因组的混合物，很难用几个确定的阶次模 拟背景模型。对于定阶次马尔克夫模型，阶次越高模型越准确，然而阶次越高，需要的数据 量也越多，一般情况下，我们获取的数据量是很难满足需求的。

发明内容

本发明的目的是针对宏基因组合宏转录组样本，提供一种宏基因组和宏转录组样本相异 度的比较方法。

本发明包括以下步骤：

步骤1:生成样本的tuple频度向量，对样本中出现的长度为1～10的tuple的频度进行 统计，并生成相应样本的频度向量；

步骤2:计算tuple的马尔克夫概率，基于变阶次马尔克夫模型估计频度向量中每一个 tuple的马尔克夫概率；

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

步骤4：生成聚类树，根据相异度矩阵生成一个聚类树。

在步骤1中，所述样本中可能出现的字符串组合为tuple元素，并选择长度为1～10的 字符串组合作为tuple元素。

在步骤2中，所述计算tuple的马尔克夫概率的具体方法可为：

步骤2-1：基于样本的频度向量构建前缀树；

步骤2-2：基于相对熵对所述前缀树进行剪枝；

步骤2-3：基于剪枝后的前缀树计算tuple的马尔克夫概率。

在步骤2-1中，所述基于样本的频度向量构建前缀树时，前缀树父节点和子节点的关系 是：子节点表示的tuple包含父节点表示的tuple，并且子节点tuple比父节点tuple多出 的一个字符出现在父节点表示的tuple之前；例如，父节点tuple为CGT，则子节点tuple 可能为ACGT，CCGT，TCGT或者GCGT。

在步骤2-2中，所述基于相对熵对所述前缀树进行剪枝时，通过计算父节点表示的tuple 与子节点表示的tuple之间的相对熵判断是否剪去子节点：当相对熵小于一定的阈值K时， 剪掉相应的子节点，相对熵D_KL的计算公式如下：

$D_{KL}=\Sigma \hat{P}\left(X\right|\mu \omega )log\left(\frac{\hat{P}\left(X\right|\mu \omega )}{\hat{P}\left(X\right|\omega )}\right)<K---\left(8\right)$

$\hat{P}\left(X\right|\omega )=\frac{N\left(\omega X\right)}{N\left(\omega \right)},\hat{P}\left(X\right|\omega \mu )=\frac{N\left(\mu \omega X\right)}{N\left(\mu \omega \right)}---\left(9\right)$

其中，ω表示父节点，μω表示子节点，X表示下一个时刻的状态，表示由μω 转移到X的转移概率，表示由ω转移到状态X的转移概率，N(ω)表示字符串ω 的频度，N(ωX)表示字符串ωX的频度，N(μω)表示字符串μω的频度，N(μωX)表示字 符串μωX的频度；所述阈值K由赤池信息量准则确定，具体公式如下：

其中，表示样本的伪似然度，d表示测序深度，表示剪枝后的前缀树 的节点个数，表示自由参数的选择范围，表示自由参数的个数，选择使的值 最小的K作为剪枝的阈值。

在步骤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^{*}({\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(12\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(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的频度向量，表示样本X第i个 tuple的频度，表示样本Y的第i个tuple的频度，n_X表示样本X中tuple个数的总和， n_Y表示样本Y中tuple个数的总和，p_X,i表示样本X中第i个tuple的马尔克夫概率，p_Y,i表示样本Y中第i个样本的马尔克夫概率。

在步骤4中，所述生成聚类树，根据相异度矩阵生成一个聚类树可根据层次聚类算法由 相异度矩阵得到聚类树。

由于定阶次马尔克夫模型存在以上局限性，本发明提出了应用变阶次马尔克夫模型的基 于k-tuple频度的微生物群落的比较方法。定阶次马尔克夫模型是变阶次马尔克夫模型的一 种特殊情况，在变阶次马尔克夫模型中，阶次根据数据性质可以为任意阶次，无需人为选择 阶次。变阶次马尔克夫模型最大的优势是其在实际应用中的灵活性和适应性。

与现有技术相比，本发明具有如下优点：本发明使用的方法无需人工选择马尔克夫阶次， 能根据数据特效自动地选择马尔克夫阶次；本发明对宏基因组和宏转录组数据的聚类效果明 显优于定阶次马尔克夫模型的聚类效果。

具体实施方式

以下将详细说明本发明的实施方式，借此本发明的实施人员可充分理解本发明如何利用 技术手段来解决技术问题。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各 个实施例的各个特征可以相互结合，所形成的技术方案均在本发明的技术保护范围之内。

利用高通量测序技术能从微生物群落中获取大量的宏基因组合宏转录组样本，通过比较 这些宏基因组或者宏转录组样本，能够更深入地了解微生物和环境之间的关系。

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

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

为了计算tuple元素的变阶次马尔克夫概率在本实施例中需要执行步骤2。在步骤2中， 首先执行步骤2-1：根据样本的所有tuple建立一个前缀树。前缀树父节点和子节点的关系 如下，父节点表示的tuple包含在子节点表示的tuple中，并且子节点tuple比父节点tuple 多出的一个字符出现在父节点tuple之前。

在本实施例步骤2中，为了对实施例中的前缀树进行剪枝操作，需要执行步骤2-2,依次 判断前缀树中每一个叶子节点是否可以被剪去。剪枝策略是通过计算叶子节点所代表的tuple 与其父节点所代表的tuple之间的相对熵，当它们之间的相对熵小于某一个特定的阈值K时， 叶子节点将被剪去。

在步骤2中，按照步骤2-2所示的方法循环检查叶子节点是否符合剪枝条件，直到没有 叶子节点可以被剪去时，剪枝完成，进而进行步骤2中的2-3子步骤操作。由剪枝的过程可 知，被剪去叶子转移到下一个状态的概率可以用与其最接近的未被剪去的祖先节点的转移概 率替代。根据这一原则，可以估算出每一个tuple的马尔可夫概率。

为了计算实施例中k-tuple向量之间的距离，接下来实施步骤3。对k-tuple向量分别采 取和的相异度方法计算距离。其中用到的tuple的变阶次马尔克夫概率，在步骤2中 已求得。

实施例步骤3可以得到一个相异度矩阵，由这个相异度矩阵进行步骤4,即进行层次聚类 分析，可最终得到一个聚类树。通过观察聚类树，可以判断聚类情况的好坏。

虽然本发明公开的实施方式如上，但所述的内容只是为了便于理解本发明采样的实施方 式，本发明所述的方法还可有多种实施例。在不背离本发明实质的情况下，熟悉本领域的技 术人员当可根据本发明做出各种相应的改变或变形，但这些相应的改变或变形都应属于本发 明。