基于不同尺度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频度特征向量的维度就是4k维。之前有研究表明，来自同一个基因组的k-tuple频度相近，能够建立tuple分布的总体统计模型，但不同基因组的k-tuple频度有很大区别。在基于k-tuple频度的无需比对的研究方法中，集中于短tuple(2-10bp)上，且应用于无监督的样本聚类上表现比较优越。因此，基于短k-tuple频度的相异度距离度量方法D2被提出用来评估比较两个微生物群落样本之间的相异度。此后，在D2基础上又衍生出和为了更好地应用到高通量测序数据上，通过归一化处理改进的和相继被提出用于比较样本之间的相异度。用和计算距离时需要结合合适的背景模型来进行建模。在之前的研究中，用到的是定阶次马尔科夫模型和基于变阶次的马尔科夫模型。然而由于微生物群落是由各种各样不同种类的微生物基因组混合组成，很难用几个确定的阶次模拟背景模型，而且需要手动设定模型的阶次，然后去集中对不同阶次模型对聚类结果的优良效果做评估，工作量和计算成本都非常高。对于定阶次马尔科夫模型，阶次越高模型越准确，然而阶次越高，需要的数据量也越多，一般情况下，我们获取的数据量是很难满足需求的。而且基于变阶次的马尔科夫模型在对选择模型阶次时，对其构建的前缀树进行减枝的过程中，需要手动设定阈值，大大提高了模型的不精确性和计算的复杂性。
技术实现思路
本专利技术的主要目的在于克服现有技术中的上述缺陷，提出一种通过变尺度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中的每一列的字符放在一个向量中，形成A1,A2,…,Ak k个向量，分别计算前面k-1个向量与最后一个向量之间的互信息量，互信息量的公式如下： I ( A w , B ) = Σ i Σ j P ( a i , b j ) * l o g P ( a i ) * P ( b本文档来自技高网...

【技术保护点】
基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，包括如下步骤：步骤1：获取宏基因组样本的2‑10bp的短tuple高通量测序数据，采用插值上下文马尔科夫模型进行建模微生物群落的背景基因组，再采用无监督的聚类方法来比较宏基因组样本，得出宏基因组样本的类别信息；步骤2：基于步骤1)中聚类得出的类别信息，将≥30bp的长tuple作为特征，采用有监督的样本分类方法找出描述宏基因组样本类别的特异性特征长tuple序列。

【技术特征摘要】
1.基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，包括如下步骤：步骤1：获取宏基因组样本的2-10bp的短tuple高通量测序数据，采用插值上下文马尔科夫模型进行建模微生物群落的背景基因组，再采用无监督的聚类方法来比较宏基因组样本，得出宏基因组样本的类别信息；步骤2：基于步骤1)中聚类得出的类别信息，将≥30bp的长tuple作为特征，采用有监督的样本分类方法找出描述宏基因组样本类别的特异性特征长tuple序列。2.如权利要求1所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，所述步骤1)具体包括如下：步骤1.1:获取宏基因组样本的高通量测序数据，生成该宏基因组样本的短tuple特征频度向量，对每个宏基因组样本中出现的长度为2-10bp的tuple的频度进行统计，并生成相应宏基因组样本的频度向量；步骤1.2:采用插值上下文马尔科夫模型对微生物群落的背景基因组进行建模，估计频度向量中每一个tuple的马尔科夫概率；步骤1.3：计算各个宏基因组样本频度向量间的相异度距离，生成一个宏基因组样本间的相异度矩阵；步骤1.4：根据相异度矩阵生成一个聚类树，用于判断宏基因组样本与样本之间的关系，找出样本的类别信息。3.如权利要求2所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，所述步骤1.1中，tuple特征定义为宏基因组样本中可能出现的字符串组合，并选择长度为2-10bp的字符串组合作为所述tuple特征。4.如权利要求2所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，步骤1.2中，计算tuple的马尔科夫概率具体方法如下：步骤1.2.1：基于宏基因组样本的频度向量的互信息量构建上下文序列树；步骤1.2.2：基于上下文序列树计算每个tuple的马尔科夫概率。5.如权利要求4所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，步骤1.2.1具体如下：基于宏基因组样本的频度向量，将k长度tuple中的每一列的字符放在一个向量中，形成A1,A2,…,Akk个向量，分别计算前面k-1个向量与最后一个向量之间的互信息量，互信息量的公式如下： I ( A w , B ) = Σ i Σ j P ( a i , b j ) * l o g P ( a i ) * P ( b j ) P ( a i , b j ) ]]>其中,w＝1,2,...,k-1；B＝Ak；ai,bj表示向量A,B中的变量；P(ai,bj)表示ai,bj在对应向量中同时出现的联合概率；P(ai)表示ai在对应向量中出现的概率；找出与B互信息量最大的向量Aw，将此向量对应的下标位置作为上下文序列树的顶点；然后按照该向量中出现的四种不同字符(A,C,G,T)将所有tuple分成四组，；最后对四组tuple向量矩阵，按照公式中的互信息量公式分别计算互信息量，找出每组中与B互信息量最大的向量As，s＝1,2,...,w-1,w+1,...,k-1，将此向量对应的下标位置作为上下文序列树的对应叶子的子节点(A,C,G,T)；依次继续下去，直到找到最后一个跟当前状态向量关联性最大的向量，整个上下文序列树构建完毕。6.如权利要求4所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，在步骤1.2.2中，每个tuple的马尔科夫概率公式如下：P(c1c2…ck)＝PICM(c1)PICM(c2|c1)…PICM(ck|c1c2…ck-1)其中，c1c2…ck表示k-tuple序列，PICM(ck|c1c2…ck-1)表示上下文序列c1c2…ck-1转移到当前状态ck的ICM转移概率。7.如权利要求6所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，对于计算上述每个ICM转移概率，针对所述k-tuple序列c1c2…ck，应用ICM马尔科夫模型构建的上下文序列树中找出与当前状态ck关联性程度从大到小排序的重要位置，重新构建其上下文序列，具体如下：要构建r阶马尔科夫模型，r≤k-1，从上下文序列树中找出与当前状态ck关联性程度从大到小排序的重要位置对应的状态分别为c3,c4…,cr，组成插值上下文序列Mr,那么构建ICM的概率模型如下：PICM(ck|c1c2c3…ck-1)＝PICM,r(ck|Mr)；PICM,r(ck|Mr)＝λr*P(ck|Mr)+(1-λr)*PICM,r-1(ck|Mr-1)； P ( c k | M r ) = Σ a ∈ A , C , G , T N ( M r , A k = c k ) N ( M r , A k = a ) ; ]]>其中，*表示乘积，λm表示m阶马尔科夫模型概率所占的权重系数；N(Mr,Xk＝ck)表示所有的k-tuple中插值上下文序列是Mr，第k个位置是ck的所有tuple的频度之和，对于上述马尔科夫模型概率所占的权重系数的计算公式为：其中所述C表示样本阈值，它由赤池信息量准则AICR(C)确定，具体公式如下：AICR(C)＝-2λ(S；Mk)+2|MIMM,k,C|；其中，λ(S；Mk)表示样本S的伪似然度，计算公式为： λ ( S ; M k ) = Σ c 1 , ... c k ∈ ( A , C , G , T ) N ( c 1 , ... c k - 1 ) logP I C M ( c k | c 1 , ... c k - 1 ) ; ]]>|MIMM,k,C|表示模型自由参数的个数，当AICR(C)值最小时算出的C值作为样本的阈值；所述q表示两个字符串之间差异度的卡方检验值，其计算原理如下： Δ r ( M r ) = Σ a ∈ A , C , G , T ( N ( M r , a ) - E ( M r , a ) ) 2 E ( M r , a ) ; ]]>E(Mr,a)＝N(Mr)PICM,r(a|Mr)；其中，N(Mr,a)，E(Mr,a)分别表示字符串频度的实际值和理论值，将q＝Δr(Mr)作为卡方值，自由度为3，作为卡方检验的指标参数。8.如权利要求2所述基于不同尺度tuple词频的微生物高通量测序数据分析协议，其特征在于，步骤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 ; ]]> d 2 S ( c ~ X , c ~ Y ) = 1 2 ( 1 - D 2 S ( c ~ X , c ~ Y ) Σ i = 1 4 k C ~ X , i 2 C ~ X , i 2 ...

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