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

一种宏基因组和宏转录组样本相异度的比较方法，涉及信息和生物技术。生成样本的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频度的相异度方法D2被提出用来评估比较两个生物样本之间的距离。此后，在D2基础上改进的和相继被提出用于比较样本之间的距离。用和计算距离时需要用到一个合适的背景模型。在之前的研究中，用到的是定阶次马尔克夫模型。然而由于微生物群落是各种基因组的混合物，很难用几个确定的阶次模拟背景模型。对于定阶次马尔克夫模型，阶次越高模型越准确，然而阶次越高，需要的数据量也越多，一般情况下，我们获取的数据量是很难满足需求的。
技术实现思路
本专利技术的目的是针对宏基因组合宏转录组样本，提供一种宏基因组和宏转录组样本相异度的比较方法。本专利技术包括以下步骤：步骤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时，剪掉相应的子节点，相对熵DKL的计算公式如下： D K L = Σ P ^ ( X | μ ω ) l o g ( P ^ ( X | μ ω ) P ^ ( X | ω ) ) < K - - - ( 8 ) ]]> P ^ ( X | ω ) = N ( ω X ) N ( ω ) , P ^ ( X | ω μ ) = N ( μ ω X ) N ( μ ω ) - - - ( 9 ) ]]>其中，ω表示父节点，μω表示子节点，X表示下一个时刻的状态，表示由μω转移到X的转移概率，表示由ω转移到状态X的转移概率，N(ω)表示字符串ω的频度，N(ωX)表示字符串ωX的频度，N(μω)表示字符串μω的频度，N(μωX)表示字符串μωX的频度；所述阈值K由赤池信息量准则确定，具体公式如下：其中，表示样本的伪似然度，d表示测序深度，表示剪枝后的前缀树的节点个数，表示自由参数的选择范围，表示自由参数的个数，选择使的值最小的K作为剪枝的阈值。在步骤3中，所述计算各个样本频度向量间的距离可采用不同的相异度方法计算各个样本频率向量间的相异度距离，所用到的相异度方法包括和计算公式如下： 本文档来自技高网...

【技术保护点】
一种宏基因组和宏转录组样本相异度的比较方法，其特征在于包括以下步骤：步骤1:生成样本的tuple频度向量，对样本中出现的长度为1～10的tuple的频度进行统计，并生成相应样本的频度向量；步骤2:计算tuple的马尔克夫概率，基于变阶次马尔克夫模型估计频度向量中每一个tuple的马尔克夫概率；步骤3：生成样本间相异度矩阵，计算各个样本频度向量间的距离，生成一个样本间的相异度矩阵；步骤4：生成聚类树，根据相异度矩阵生成一个聚类树。

【技术特征摘要】
1.一种宏基因组和宏转录组样本相异度的比较方法，其特征在于包括以下步骤：
步骤1:生成样本的tuple频度向量，对样本中出现的长度为1～10的tuple的频度进行
统计，并生成相应样本的频度向量；
步骤2:计算tuple的马尔克夫概率，基于变阶次马尔克夫模型估计频度向量中每一个
tuple的马尔克夫概率；
步骤3：生成样本间相异度矩阵，计算各个样本频度向量间的距离，生成一个样本间的相
异度矩阵；
步骤4：生成聚类树，根据相异度矩阵生成一个聚类树。
2.如权利要求1所述一种宏基因组和宏转录组样本相异度的比较方法，其特征在于在步
骤1中，所述样本中可能出现的字符串组合为tuple元素，并选择长度为1～10的字符串组
合作为tuple元素。
3.如权利要求1所述一种宏基因组和宏转录组样本相异度的比较方法，其特征在于在步
骤2中，所述计算tuple的马尔克夫概率的具体方法为：
步骤2-1：基于样本的频度向量构建前缀树；
步骤2-2：基于相对熵对所述前缀树进行剪枝；
步骤2-3：基于剪枝后的前缀树计算tuple的马尔克夫概率。
4.如权利要求3所述一种宏基因组和宏转录组样本相异度的比较方法，其特征在于在步
骤2-1中，所述基于样本的频度向量构建前缀树时，前缀树父节点和子节点的关系是：子节
点表示的tuple包含父节点表示的tuple，并且子节点tuple比父节点tuple多出的一个字
符出现在父节点表示的tuple之前；若父节点tuple为CGT，则子节点tuple可能为ACGT，
CCGT，TCGT或者GCGT。
5.如权利要求3所述一种宏基因组和宏转录组样本相异度的比较方法，其特征在于在步
骤2-2中，所述基于相对熵对所述前缀树进行剪枝时，通过计算父节点表示的tuple与子节
点表示的tuple之间的相对熵判断是否剪去子节点：当相对熵小于一定的阈值K时，剪掉相
应的子节点，相对熵DKL的计算公式如下：
D K L = Σ P ^ ( X | μ ω ) l o g ( P ^ ( X | μ ω ) P ^ ( X | ω ) ) < K - - - ( 1 ) ]]> P ^ ( X | ω ) = N ( ω X ) N ( ω ) , P ^ ( X | ω μ ) = N ( μ ω X ) N ( μ ω ) - - - ( 2 ) ]]>其中，ω表示父节点，μω表示子节点，X表示下一个时刻的状态，表示由μω
转移到X的转移概率，表示由ω转移到状态X的转移概率，N(ω)表示字符串ω
的频度，N(ωX)表示字符串ωX的频度，N(μω)表示字符串μω的频度，N(μωX)表示字
符串μωX的频度；所述阈值K由赤池信息量准则确定，具体公式如下：
其中，表示样本的伪似然度，d表示测序深度，表示剪枝后的前缀树
的节点个数，表示自由参数的选择范围，表示自由参数的个数，选择使的值
最小的K作为剪枝的阈值。
6.如权利要求1所述一种宏基因组和宏转录组样本相异度的比较方法，其特征在于在
步骤3中，所述计算各个样本频度向量间的距离可采用不同的相异度方法计算各个样本频率
向量间的相异度距离，所用到的相异度方法包括和计算公式如下：
D 2 S ( c ~ X , c ~ Y ) = Σ i = 1 4 k c ~ X , i c ~ Y , i c ~ X , i 2 + c ~ Y , i 2 - - - ( 4 ) ]]> D 2 * ( c ~ X , c ~ Y ) = Σ i = 1 4 k c ~ X , i c ~ Y , i n X p X , i n Y p Y , i - - - ( 5 ) ]]> d 2 S ( c ~ X , c ~ Y ) = 1 2 ( 1 - D 2 S ( c ~ X , c ~ Y ...

【专利技术属性】
技术研发人员：王颖，廖伟楠，刘暾东，
申请(专利权)人：厦门大学，
类型：发明
国别省市：福建;35