CN105787296B - 一种宏基因组和宏转录组样本相异度的比较方法 - Google Patents
一种宏基因组和宏转录组样本相异度的比较方法 Download PDFInfo
- Publication number
- CN105787296B CN105787296B CN201610100159.XA CN201610100159A CN105787296B CN 105787296 B CN105787296 B CN 105787296B CN 201610100159 A CN201610100159 A CN 201610100159A CN 105787296 B CN105787296 B CN 105787296B
- Authority
- CN
- China
- Prior art keywords
- sample
- tuple
- dissimilarity
- frequency
- markov
- 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
Links
- 238000013459 approach Methods 0.000 title abstract description 5
- 230000000052 comparative effect Effects 0.000 title abstract 2
- 239000013598 vector Substances 0.000 claims abstract description 26
- 239000011159 matrix material Substances 0.000 claims abstract description 15
- 238000000034 method Methods 0.000 claims description 30
- 238000013138 pruning Methods 0.000 claims description 11
- 235000012571 Ficus glomerata Nutrition 0.000 claims description 6
- 244000153665 Ficus glomerata Species 0.000 claims description 6
- 230000007704 transition Effects 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000012163 sequencing technique Methods 0.000 claims description 2
- 238000000691 measurement method Methods 0.000 claims 1
- 230000000694 effects Effects 0.000 abstract description 4
- 238000009394 selective breeding Methods 0.000 abstract 1
- 239000000523 sample Substances 0.000 description 27
- 230000000813 microbial effect Effects 0.000 description 6
- 108090000623 proteins and genes Proteins 0.000 description 5
- 238000012165 high-throughput sequencing Methods 0.000 description 4
- 239000012472 biological sample Substances 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 1
- 238000007417 hierarchical cluster analysis Methods 0.000 description 1
- 244000005700 microbiome Species 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biophysics (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
一种宏基因组和宏转录组样本相异度的比较方法，涉及信息和生物技术。生成样本的tuple频度向量，对样本中出现的长度为1～10的tuple的频度进行统计，并生成相应样本的频度向量；计算tuple的马尔克夫概率，基于变阶次马尔克夫模型估计频度向量中每一个tuple的马尔克夫概率；生成样本间相异度矩阵，计算各个样本频度向量间的距离，生成一个样本间的相异度矩阵；生成聚类树，根据相异度矩阵生成一个聚类树。无需人工选择马尔克夫阶次，能根据数据特效自动地选择马尔克夫阶次；对宏基因组和宏转录组数据的聚类效果明显优于定阶次马尔克夫模型的聚类效果。
Description
技术领域
本发明涉及信息和生物技术，尤其是涉及一种宏基因组和宏转录组样本相异度的比较方法。
背景技术
微生物群落间的比较对于理解微生物和环境之间的关系至关重要。高通量测序技术已经成为表征微生物群落的一个强有力的工具。对于不同基因间的比较，基于配准的序列比较方法，如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的计算公式如下：
其中，ω表示父节点，μω表示子节点，X表示下一个时刻的状态，表示由μω转移到X的转移概率，表示由ω转移到状态X的转移概率，N(ω)表示字符串ω的频度，N(ωX)表示字符串ωX的频度，N(μω)表示字符串μω的频度，N(μωX)表示字符串μωX的频度；所述阈值K由赤池信息量准则确定，具体公式如下：
其中，表示样本的伪似然度，d表示测序深度，表示剪枝后的前缀树的节点个数，表示自由参数的选择范围，表示自由参数的个数，选择使的值最小的K作为剪枝的阈值。
在步骤3中，所述计算各个样本频度向量间的距离可采用不同的相异度方法计算各个样本频率向量间的相异度距离，所用到的相异度方法包括和计算公式如下：
其中，表示样本X的频度向量，表示样本Y的频度向量，表示样本X第i个tuple的频度，表示样本Y的第i个tuple的频度，nX表示样本X中tuple个数的总和，nY表示样本Y中tuple个数的总和，pX,i表示样本X中第i个tuple的马尔克夫概率，pY,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,即进行层次聚类分析，可最终得到一个聚类树。通过观察聚类树，可以判断聚类情况的好坏。
虽然本发明公开的实施方式如上，但所述的内容只是为了便于理解本发明采样的实施方式，本发明所述的方法还可有多种实施例。在不背离本发明实质的情况下，熟悉本领域的技术人员当可根据本发明做出各种相应的改变或变形，但这些相应的改变或变形都应属于本发明。
Claims (7)
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字符串的左端；例如子节点ACGT的父节点tuple为CGT。
5.如权利要求3所述一种宏基因组和宏转录组样本相异度的比较方法，其特征在于在步骤2-2中，所述基于相对熵对所述前缀树进行剪枝时，通过计算父节点表示的tuple与子节点表示的tuple之间的相对熵判断是否剪去子节点：当相对熵小于一定的阈值K时，剪掉相应的子节点，相对相对熵DKL的计算公式如下：
其中ω表示父节点，μω表示子节点，X表示下一个时刻的状态，表示由μω转移到X的转移概率，表示由ω转移到状态X的转移概率，N(ω)表示字符串ω的频度，N(ωX)表示字符串ωX的频度，N(μω)表示字符串μω的频度，N(μωX)表示字符串μωX的频度；所述阈值K由赤池信息量准则确定，具体公式如下：
其中R表示样本数据集，表示以K为阈值，对样本数据R采用步骤2计算求得的马尔科夫概率，表示对步骤2求得的马尔科夫概率的自然对数，即为样本的伪似然度，
d表示测序深度，表示步骤2-2中剪枝后的前缀树，表示剪枝后的前缀树的节点个数，表示自由参数的选择范围，表示自由参数的个数，选择使的值最小的K作为剪枝的阈值。
6.如权利要求1所述一种宏基因组和宏转录组样本相异度的比较方法，其特征在于在步骤3中，所述计算样本间的相异度，可采用所有基于频度概率期望的距离度量方法计算各个样本频率向量间的相异度距离，例如和计算公式如下：
其中，表示样本X的频率向量，表示样本Y的频度向量，表示样本X第i个tuple的频度，表示样本Y的第i个tuple的频度。nX表示样本X中tuple个数的总和，nY表示样本Y中tuple个数的总和,pX,i表示样本X中第i个tuple的马尔克夫概率，pY,i表示样本Y中第i个样本的马尔克夫概率。和是两种非中心化的相异度距离度量定义，被称为相异度和相异度，而是相异度中心化后的相异度距离度量定义，被称为相异度，是相异度中心化后的相异度距离度量定义，被称为相异度。表示样本X和样本Y的相异度，表示样本X和样本Y的相异度，表示样本X和样本Y的相异度，表示样本X和样本Y的相异度。
7.如权利要求1所述一种宏基因组和宏转录组样本相异度的比较方法，其特征在于在步骤4中，所述生成聚类树，根据相异度矩阵生成一个聚类树是根据层次聚类算法由相异度矩阵得到聚类树。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610100159.XA CN105787296B (zh) | 2016-02-24 | 2016-02-24 | 一种宏基因组和宏转录组样本相异度的比较方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610100159.XA CN105787296B (zh) | 2016-02-24 | 2016-02-24 | 一种宏基因组和宏转录组样本相异度的比较方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105787296A CN105787296A (zh) | 2016-07-20 |
| CN105787296B true CN105787296B (zh) | 2018-07-17 |
Family
ID=56402865
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610100159.XA Expired - Fee Related CN105787296B (zh) | 2016-02-24 | 2016-02-24 | 一种宏基因组和宏转录组样本相异度的比较方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105787296B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106202999B (zh) * | 2016-07-21 | 2018-12-11 | 厦门大学 | 基于不同尺度tuple词频的微生物高通量测序数据分析协议 |
| CN110634538A (zh) * | 2019-08-26 | 2019-12-31 | 上海科技发展有限公司 | 利福平耐药结核菌的检测方法、装置、设备和存储介质 |
| CN111564179B (zh) * | 2020-05-09 | 2022-04-29 | 厦门大学 | 一种基于三元组神经网络的物种生物学分类方法及系统 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2012097152A2 (en) * | 2011-01-13 | 2012-07-19 | Laboratory Corporation Of America Holdings | Methods and systems for predictive modeling of hiv-1 replication capacity |
| CN104616264A (zh) * | 2015-02-12 | 2015-05-13 | 厦门大学 | 基因芯片图像的自动对比度增强方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB0317335D0 (en) * | 2003-07-24 | 2003-08-27 | Sec Dep For The Home Departmen | Improvements in and relating to interpretation |
-
2016
- 2016-02-24 CN CN201610100159.XA patent/CN105787296B/zh not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2012097152A2 (en) * | 2011-01-13 | 2012-07-19 | Laboratory Corporation Of America Holdings | Methods and systems for predictive modeling of hiv-1 replication capacity |
| EP2663943A2 (en) * | 2011-01-13 | 2013-11-20 | Laboratory Corporation of America Holdings | Methods and systems for predictive modeling of hiv-1 replication capacity |
| CN104616264A (zh) * | 2015-02-12 | 2015-05-13 | 厦门大学 | 基因芯片图像的自动对比度增强方法 |
Non-Patent Citations (2)
| Title |
|---|
| "Comparison of Metatranscriptomic Samples Based on k-Tuple Frequencies";王颖 等;《PloS one》;20140930;第9卷(第1期);第1-8页 * |
| "基于k-tuple频度统计的微生物群落测序数据分析";刘麟;《中国优秀硕士论文全文数据库 基础科学辑》;20140815(第8期);第6-7、9-10、19页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105787296A (zh) | 2016-07-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| WO2017157183A1 (zh) | 一种自动多阀值特征过滤方法及装置 | |
| CN106202999B (zh) | 基于不同尺度tuple词频的微生物高通量测序数据分析协议 | |
| CN110516818A (zh) | 一种基于集成学习技术的高维度数据预测方法 | |
| CN111564179B (zh) | 一种基于三元组神经网络的物种生物学分类方法及系统 | |
| CN107145523B (zh) | 基于迭代匹配的大型异构知识库对齐方法 | |
| CN105787296B (zh) | 一种宏基因组和宏转录组样本相异度的比较方法 | |
| CN113377964B (zh) | 知识图谱链接预测方法、装置、设备及存储介质 | |
| CN107229945A (zh) | 一种基于竞争学习的深度聚类方法 | |
| Wang et al. | GAEM: a hybrid algorithm incorporating GA with EM for planted edited motif finding problem | |
| CN117059169A (zh) | 基于参数自适应成长优化器的生物多序列比对方法及系统 | |
| CN110491443B (zh) | 一种基于投影邻域非负矩阵分解的lncRNA蛋白质关联预测方法 | |
| Comin et al. | Fast entropic profiler: An information theoretic approach for the discovery of patterns in genomes | |
| CN113658641B (zh) | 一种噬菌体分类方法、装置、设备及存储介质 | |
| CN113066528B (zh) | 基于主动半监督图神经网络的蛋白质分类方法 | |
| Roos et al. | Analysis of textual variation by latent tree structures | |
| KR101522087B1 (ko) | 미스매치를 고려한 염기 서열 정렬 시스템 및 방법 | |
| CN106156107A (zh) | 一种新闻热点的发现方法 | |
| CN114400043B (zh) | 基于孪生神经网络的半监督宏基因组分箱方法 | |
| CN116504315A (zh) | 一种基于改进的began网络的单细胞rna测序数据缺失的插补方法 | |
| Elishco et al. | The capacity of some Pólya string models | |
| CN104462817A (zh) | 基于蒙特卡洛和非负矩阵因子分解的基因选择和癌症分类方法 | |
| CN115543803A (zh) | 基于改进遗传算法的软件测试用例智能生成方法及系统 | |
| CN105162648A (zh) | 基于骨干网络扩展的社团检测方法 | |
| Liu et al. | Discovery of deep order-preserving submatrix in DNA microarray data based on sequential pattern mining | |
| CN113033768A (zh) | 一种基于图卷积网络的缺失特征重表示方法及系统 |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20180717 |