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

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

Info

Publication number
CN106202999B
CN106202999B CN201610577084.4A CN201610577084A CN106202999B CN 106202999 B CN106202999 B CN 106202999B CN 201610577084 A CN201610577084 A CN 201610577084A CN 106202999 B CN106202999 B CN 106202999B
Authority
CN
China
Prior art keywords
sample
tuple
frequency
vector
indicate
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201610577084.4A
Other languages
English (en)
Other versions
CN106202999A (zh
Inventor
王颖
汪顺
刘暾东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xiamen University
Original Assignee
Xiamen University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xiamen University filed Critical Xiamen University
Priority to CN201610577084.4A priority Critical patent/CN106202999B/zh
Publication of CN106202999A publication Critical patent/CN106202999A/zh
Application granted granted Critical
Publication of CN106202999B publication Critical patent/CN106202999B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding

Landscapes

  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Software Systems (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Public Health (AREA)
  • Bioethics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

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

Description

基于不同尺度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个向量与最后一个向量之间的互信息量,互信息量的公式如下:
其中,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);依次继续下去,直到找到最后一个跟当前状态向量关联性最大的向量,整个上下文序列树构建完毕。
优选的,在步骤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转移概率。
优选的,对于计算上述每个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);
其中,*表示乘积,λ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的伪似然度,计算公式为:
|MIMM,k,C|表示模型自由参数的个数,当AICR(C)值最小时算出的C值作为样本的阈值;
所述q表示两个字符串之间差异度的卡方检验值,其计算原理如下:
E(Mr,a)=N(Mr)PICM,r(a|Mr);
其中,N(Mr,a),E(Mr,a)分别表示字符串频度的实际值和理论值,将q=Δr(Mr)作为卡方值,自由度为3,作为卡方检验的指标参数。
优选的,步骤1.3中,应用不同的相异度距离度量方法计算各个宏基因组样本频度向量间的相异度距离,所用到的相异度距离度量方法包括计算公式如下:
其中,都是计算两个样本之间的相异度的一种距离度量方法;表示样本X的频度向量,表示样本Y的频度向量; 称为中心化过程,i=1,2,…,4k;CX,i和CY,i分别表示X和Y样本中第i个tuple出现的频度;nX表示样本X中tuple个数的总和,nY表示样本Y中tuple个数的总和;pX,i和pY,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(x,y)是两个宏基因组样本的相异度距离,|Ci|和|Cj|表示两个群组的大小,也就是群组中样本的个数,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序列特征过滤移除,得到一些类别特异性候选特征,对称不确定性定义如下:
其中,NX表示当前类别组成的X样本集中tuple特征出现的频度;sum(NX)表示由当前类别组成的X样本集中特征出现的频度之和;sum(NY)表示目标类别组成的Y样本集中特征出现的频度之和;n(X)和n(Y)分别表示X和Y样本集中样本的个数;θ表示X和Y之间对称不确定性的阈值;
采用SVM分类器,对样本进行有监督分类,找到能够描述微生物群落内部具有差异性的特异性特征。
优选的,在步骤2.4中,基于步骤2.3获得的特异性特征,用留一法(LOOCV)来验证和评估分类器的准确率P:
其中,P表示分类准确率,D是由有限个以(xi,yi)形式表示的样本组合集合xi是样本中除yi以外的属性列表,yi表示样本中类别标号的属性,g表示分类器模型函数,输出结果为该模型的预测结果,f(g(xi),yi)为判别函数,当g(xi)与yi相等时,输出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中的每一列的字符放在一个向量中,形成A1,A2,…,Ak k个向量。分别计算前面k-1个向量与最后一个向量(即当前状态向量)之间的互信息量,互信息量的公式如下:
其中,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向量矩阵,按照公式(1)中的互信息量公式分别计算互信息量,找出每组中与Y互信息量最大的向量As,其中s=1,2,…,w-1,w+1,…,k-1,将此向量对应的下标位置作为上下文序列树的对应叶子(A,C,G,T)的子节点;依次继续下去,直到找到最后一个跟当前状态向量关联性最大的向量,整个上下文序列树构建完毕。对于上述构建的上下文序列树,每条枝都对应一个tuple,按照与当前状态关联性程度从大到小,把这个tuple中的碱基位置自上而下排列,储存在这条枝的节点中。
步骤1.2.2:基于上下文序列树计算每个tuple的马尔科夫概率。每个tuple的马尔科夫概率公式如下:
P(c1c2…ck)=PICM(c1)PICM(c2|c1)…PICM(ck|c1c2…ck-1) (2)
其中,c1c2…ck表示k-tuple序列,PICM(ck|c1c2…ck-1)表示上下文序列c1c2…ck-1转移到当前状态ck的ICM转移概率。
对于计算上述每个ICM转移概率,针对上述k-tuple序列c1c2…ck,应用ICM马尔科夫模型构建的上下文序列树中找出与当前状态ck关联性程度从大到小排序的重要位置,重新构建其上下文序列。例如,要构建r阶马尔科夫模型(r≤k-1),从上下文序列树中找出与当前状态bk关联性程度从大到小排序的重要位置对应的状态分别为c3,c4,…,cr,组成插值上下文序列Mr,那么构建ICM的概率模型如下:
PICM(ck|c1c2c3…ck-1)=PICM,r(ck|Mr) (3)
PICM,r(ck|Mr)=λr*P(ck|Mr)+(1-λr)*PICM,r-1(ck|Mr-1) (4)
其中,λm表示m阶马尔科夫模型概率所占的权重系数;N(Mr,Ak=ck)表示所有的k-tuple中插值上下文序列是Mr,第k个位置是ck的所有tuple的频度之和。对于上述马尔科夫模型概率所占的权重系数的计算公式为:
其中所述C表示样本阈值,它由赤池信息量准则AICR(C)确定,具体公式如下:
AICR(C)=-2λ(S;Mk)+2|MIMM,k,C| (7)
其中,λ(S;Mk)表示样本S的伪似然度,计算公式为:
|MIMM,k,C|表示模型自由参数的个数。当AICR(C)值最小时算出的C值作为样本的阈值。
所述q表示两个字符串之间差异度的卡方检验值,其计算原理如下:
E(Mr,a)=N(Mr)PICM,r(a|Mr) (10)
其中,N(Mr,a),E(Mr,a)分别表示字符串频度的实际值和理论值。将q=Δr(Mr)作为卡方值,自由度为3,作为卡方检验的指标参数。
步骤1.3:计算各个样本频度向量间的相异度距离,生成一个宏基因组样本间的相异度矩阵。应用不同的相异度距离度量方法计算各个宏基因组样本频度向量间的相异度距离,所用到的相异度距离度量方法包括计算公式如下:
其中,都是计算两个样本之间的相异度的一种距离度量方法;表示样本X的频度向量,表示样本Y的频度向量; 称为中心化过程;CX,i和CY,i分别表示X和Y样本中第i个tuple出现的频度;nX表示样本X中tuple个数的总和,nY表示样本Y中tuple个数的总和;pX,i和pY,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(x,y)是两个样本的相异度距离,|Ci|和|Cj|表示两个群组的大小,也就是群组中样本的个数,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序列特征过滤移除。得到一些类别特异性候选特征。对称不确定性定义如下:
其中,NX表示当前类别组成的X样本集中tuple特征出现的频度;sum(NX)表示由当前类别组成的X样本集中特征出现的频度之和;sum(NY)表示目标类别组成的Y样本集中特征出现的频度之和;n(X)和n(Y)分别表示X和Y样本集中样本的个数;θ表示X和Y之间对称不确定性的阈值。
采用SVM分类器,对样本进行有监督分类,找到能够描述微生物群落内部具有差异性的特异性特征。
步骤2.4:基于步骤2.3获得的特异性特征,用留一法(LOOCV)来验证和评估分类器的准确率P:
其中,P表示分类准确率,D是由有限个以(xi,yi)形式表示的样本组合集合xi是样本中除yi以外的属性列表,yi表示样本中类别标号的属性,g表示分类器模型函数,输出结果为该模型的预测结果,f(g(xi),yi)为判别函数,当g(xi)与yi相等时,输出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 al.,A human gut microbial gene catalogue established by metagenomic sequencing.Nature,2010.464(7285):p.59-65.),将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
上述仅为本发明的具体实施方式,但本发明的设计构思并不局限于此,凡利用此构思对本发明进行非实质性的改动,均应属于侵犯本发明保护范围的行为。

Claims (11)

1.基于不同尺度tuple词频的微生物高通量测序数据分析协议,其特征在于,包括如下步骤:
步骤1:获取宏基因组样本的2-10bp的短tuple高通量测序数据,采用插值上下文马尔科夫模型进行建模微生物群落的背景基因组,再采用无监督的聚类方法来比较宏基因组样本,得出宏基因组样本的类别信息;
步骤2:基于步骤1)中聚类得出的类别信息,将≥30bp的长tuple作为特征,采用有监督的样本分类方法找出描述宏基因组样本类别的特异性特征长tuple序列;
所述步骤1)具体包括如下:
步骤1.1:获取宏基因组样本的高通量测序数据,生成该宏基因组样本的短tuple特征频度向量,对每个宏基因组样本中出现的长度为2-10bp的tuple的频度进行统计,并生成相应宏基因组样本的频度向量;所述步骤1.1中,tuple特征定义为宏基因组样本中可能出现的字符串组合,并选择长度为2-10bp的字符串组合作为所述tuple特征;
步骤1.2:采用插值上下文马尔科夫模型对微生物群落的背景基因组进行建模,估计频度向量中每一个tuple的马尔科夫概率;
步骤1.3:计算各个宏基因组样本频度向量间的相异度距离,生成一个宏基因组样本间的相异度矩阵;
步骤1.4:根据相异度矩阵生成一个聚类树,用于判断宏基因组样本与样本之间的关系,找出样本的类别信息。
2.如权利要求1所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议,其特征在于,步骤1.2中,计算tuple的马尔科夫概率具体方法如下:
步骤1.2.1:基于宏基因组样本的频度向量的互信息量构建上下文序列树;
步骤1.2.2:基于上下文序列树计算每个tuple的马尔科夫概率。
3.如权利要求2所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议,其特征在于,步骤1.2.1具体如下:基于宏基因组样本的频度向量,将k长度tuple中的每一列的字符放在一个向量中,形成A1,A2,…,Akk个向量,分别计算前面k-1个向量与最后一个向量之间的互信息量,互信息量的公式如下:
其中,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之一;依次继续下去,直到找到最后一个跟当前状态向量关联性最大的向量,整个上下文序列树构建完毕。
4.如权利要求2所述的基于不同尺度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转移概率。
5.如权利要求4所述的基于不同尺度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);
其中,*表示乘积,λr表示r阶马尔科夫模型概率所占的权重系数;N(Mr,Ak=ck)表示所有的k-tuple中插值上下文序列是Mr,第k个位置是ck的所有tuple的频度之和,对于上述马尔科夫模型概率所占的权重系数的计算公式为:
其中C表示样本阈值,它由赤池信息量准则AICR(C)确定,具体公式如下:
AICR(C)=-2λ(S;Mk)+2|MIMM,k,C|;
其中,λ(S;Mk)表示样本S的伪似然度,计算公式为:
|MIMM,k,C|表示模型自由参数的个数,当AICR(C)值最小时算出的C值作为样本的阈值;
q表示两个字符串之间差异度的卡方检验值,其计算原理如下:
E(Mr,a)=N(Mr)PICM,r(a|Mr);
其中,N(Mr,a),E(Mr,a)分别表示字符串频度的实际值和理论值,将q=Δr(Mr)作为卡方值,自由度为3,作为卡方检验的指标参数,N(Mr)表示插值上下文序列为Mr的所有tuple频度之和,PICM,r(a|Mr)表示上下文序列Mr转移到当前状态a的ICM转移概率。
6.如权利要求1所述基于不同尺度tuple词频的微生物高通量测序数据分析协议,其特征在于,步骤1.3中,应用不同的相异度距离度量方法计算各个宏基因组样本频度向量间的相异度距离,所用到的相异度距离度量方法包括计算公式如下:
其中,都是计算两个样本之间的相异度的一种距离度量方法;表示样本X的频度向量,表示样本Y的频度向量; 称为中心化过程,i=1,2,...,4k;CX,i和CY,i分别表示X和Y样本中第i个tuple出现的频度;nX表示样本X中tuple个数的总和,nY表示样本Y中tuple个数的总和;pX,i和pY,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
其中,d(x,y)是两个宏基因组样本的相异度距离,如果不同样本间的距离越小,那么d(x,y)的值就越小;d(x,x)表示相同样本间的距离为0。
7.如权利要求1所述基于不同尺度tuple词频的微生物高通量测序数据分析协议,其特征在于,在步骤1.4中,在n*n相异度矩阵的基础上,根据非加权平均法层次聚类算法计算出两个群组的相异度距离,其定义如下:
d(x,y)是两个宏基因组样本的相异度距离,|Ci|和|Cj|表示两个群组的大小,也就是群组中样本的个数,i,j=1,2,...,n;由两两群组的相异度距离得到聚类树,从聚类树中可以直观的看出群落中各个样本之间的结构关系,得出样本之间的类别信息。
8.如权利要求1所述基于不同尺度tuple词频的微生物高通量测序数据分析协议,其特征在于,所述步骤2具体包括以下子步骤:
步骤2.1:对样本中出现的长度为40bp的长tuple的频度进行统计,并生成相应样本的频度向量;
步骤2.2:对每个样本的tuple频度向量进行并行处理,生成所有样本的长tuple频度向量矩阵,然后过滤掉冗余的特征;
步骤2.3:基于所述步骤1中获取的样本类别信息,应用过滤好的样本特征对样本进行有监督分类,找到对分类效果具有很强辨识性的特异性tuple特征;
步骤2.4:基于步骤2.3获得的特异性特征,用留一法来验证和评估分类器的准确率。
9.如权利要求8所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议,其特征在于,在步骤2.2中,将需要分类样本的tuple频度向量合并在一起,生成一个tuple频度向量矩阵A,A表示为M×N的频度矩阵,其中N表示样本数量,M表示特征维度。
10.如权利要求8所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议,其特征在于,在步骤2.3中,基于步骤1中获取的样本类别信息,选定训练集和测试集样本,在训练集中选定当前类别和目标类别,然后根据对称不确定性大于设定的阈值时,把冗余的tuple序列特征过滤移除,得到一些类别特异性候选特征,对称不确定性定义如下:
其中,NX表示当前类别组成的X样本集中tuple特征出现的频度;sum(NX)表示由当前类别组成的X样本集中特征出现的频度之和;sum(NY)表示目标类别组成的Y样本集中特征出现的频度之和;n(X)和n(Y)分别表示X和Y样本集中样本的个数;q表示X和Y之间对称不确定性的阈值;
采用SVM分类器,对样本进行有监督分类,找到能够描述微生物群落内部具有差异性的特异性特征。
11.如权利要求8所述的基于不同尺度tuple词频的微生物高通量测序数据分析协议,其特征在于,在步骤2.4中,基于步骤2.3获得的特异性特征,用留一法来验证和评估分类器的准确率P:
其中,P表示分类准确率,D是由有限个以(xi,yi)形式表示的样本组合集合xi是样本中除yi以外的属性列表,yi表示样本中类别标号的属性,g表示分类器模型函数,输出结果为该模型的预测结果,f(g(xi),yi)为判别函数,当g(xi)与yi相等时,输出1,否则,输出0。
CN201610577084.4A 2016-07-21 2016-07-21 基于不同尺度tuple词频的微生物高通量测序数据分析协议 Expired - Fee Related CN106202999B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610577084.4A CN106202999B (zh) 2016-07-21 2016-07-21 基于不同尺度tuple词频的微生物高通量测序数据分析协议

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610577084.4A CN106202999B (zh) 2016-07-21 2016-07-21 基于不同尺度tuple词频的微生物高通量测序数据分析协议

Publications (2)

Publication Number Publication Date
CN106202999A CN106202999A (zh) 2016-12-07
CN106202999B true CN106202999B (zh) 2018-12-11

Family

ID=57491188

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610577084.4A Expired - Fee Related CN106202999B (zh) 2016-07-21 2016-07-21 基于不同尺度tuple词频的微生物高通量测序数据分析协议

Country Status (1)

Country Link
CN (1) CN106202999B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108268753B (zh) * 2018-01-25 2021-12-03 清华大学 一种微生物组识别方法和装置、设备
CN110797088B (zh) * 2019-10-17 2020-09-15 南京医基云医疗数据研究院有限公司 全基因组重测序分析及用于全基因组重测序分析的方法
CN110782949A (zh) * 2019-10-22 2020-02-11 王文婷 一种基于最大最小化序列搜索的多层基因加权分组方法
CN111326215B (zh) * 2020-02-07 2022-04-29 厦门大学 一种基于k-tuple频度的核酸序列搜索方法及系统
CN111564179B (zh) * 2020-05-09 2022-04-29 厦门大学 一种基于三元组神经网络的物种生物学分类方法及系统
CN112863593B (zh) * 2021-02-05 2024-02-20 厦门大学 基于皮肤宏基因组数据的身份鉴定特征提取方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102203788A (zh) * 2008-10-31 2011-09-28 雅培制药有限公司 用于装配成小组的癌细胞系以用于测试一种或多种药物组合物的功效的方法
WO2014200991A1 (en) * 2013-06-10 2014-12-18 University Of Virginia Patent Foundation System, method and computer readable medium for rapid dna identification
CN105787296A (zh) * 2016-02-24 2016-07-20 厦门大学 一种宏基因组和宏转录组样本相异度的比较方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102203788A (zh) * 2008-10-31 2011-09-28 雅培制药有限公司 用于装配成小组的癌细胞系以用于测试一种或多种药物组合物的功效的方法
JP5391279B2 (ja) * 2008-10-31 2014-01-15 アッヴィ・インコーポレイテッド 1種以上の医薬組成物の有効性を試験することに使用する癌細胞系のパネルを構築するための方法
WO2014200991A1 (en) * 2013-06-10 2014-12-18 University Of Virginia Patent Foundation System, method and computer readable medium for rapid dna identification
EP3008028A1 (en) * 2013-06-10 2016-04-20 University Of Virginia Patent Foundation System, method and computer readable medium for rapid dna identification
CN105787296A (zh) * 2016-02-24 2016-07-20 厦门大学 一种宏基因组和宏转录组样本相异度的比较方法

Also Published As

Publication number Publication date
CN106202999A (zh) 2016-12-07

Similar Documents

Publication Publication Date Title
CN106202999B (zh) 基于不同尺度tuple词频的微生物高通量测序数据分析协议
CN110222745A (zh) 一种基于相似性学习及其增强的细胞类型鉴定方法
CN104040561B (zh) 通过质谱术和分数规整识别微生物的方法
Wang et al. Structured subcomposition selection in regression and its application to microbiome data analysis
CN105825078B (zh) 基于基因大数据的小样本基因表达数据分类方法
CN106682454B (zh) 一种宏基因组数据分类方法和装置
CN110659378B (zh) 基于对比相似性损失函数的细粒度图像检索方法
Liao et al. A new unsupervised binning approach for metagenomic sequences based on n-grams and automatic feature weighting
Rasheed et al. Metagenomic taxonomic classification using extreme learning machines
CN109002859B (zh) 基于主成分分析的传感器阵列特征选择和阵列优化方法
CN104966105A (zh) 一种鲁棒机器错误检索方法与系统
CN105469108B (zh) 基于生物学数据的聚类方法及系统、聚类结果评价方法及系统
CN102254033A (zh) 基于熵权重的全局k-均值聚类方法
CN101923604A (zh) 基于邻域粗糙集的加权knn肿瘤基因表达谱分类方法
CN106548041A (zh) 一种基于先验信息和并行二进制微粒群算法的肿瘤关键基因识别方法
Pouyan et al. Clustering single-cell expression data using random forest graphs
Dash et al. Performance analysis of clustering techniques over microarray data: A case study
CN107423697A (zh) 基于非线性融合深度3d卷积描述子的行为识别方法
CN107392249A (zh) 一种k近邻相似度优化的密度峰聚类方法
CN105139037A (zh) 基于最小生成树的集成多目标进化自动聚类方法
CN110070070B (zh) 一种动作识别方法
Wu On biological validity indices for soft clustering algorithms for gene expression data
KR102376212B1 (ko) 신경망 기반의 유전자 선택 알고리즘을 이용한 유전자 발현 마커 선별 방법
CN108090514B (zh) 基于两阶段密度聚类的红外图像识别方法
CN112801197A (zh) 一种基于用户数据分布的K-means方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20181211