CN101497924A - 一种基于间隙谱的生物序列分析方法 - Google Patents
一种基于间隙谱的生物序列分析方法 Download PDFInfo
- Publication number
- CN101497924A CN101497924A CNA2008100572005A CN200810057200A CN101497924A CN 101497924 A CN101497924 A CN 101497924A CN A2008100572005 A CNA2008100572005 A CN A2008100572005A CN 200810057200 A CN200810057200 A CN 200810057200A CN 101497924 A CN101497924 A CN 101497924A
- Authority
- CN
- China
- Prior art keywords
- character
- biological
- sequence
- frequency
- distance
- 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.)
- Pending
Links
Images
Landscapes
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
Abstract
本发明提供了一种生物序列的分析方法，其包括如下步骤：(1)计算生物序列的间隙谱：计算生物序列中字符之间的距离，分别统计字符之间相同距离的出现频率，构成间隙谱；(2)计算不同生物序列间的相似度；(3)推导不同生物序列的同源性或生物学功能。本发明生物序列的分析方法，具有高速度、高灵敏度和高准确性等优点。
Description
技术领域
本发明属于生物信息学领域，具体涉及一种生物序列的非比对的分析方法。
背景技术
众所周知，生物序列包括核酸、氨基酸序列，含有大量生命信息。目前，生物序列测序已经不是一件难事。在国内外的各数据库中，已积聚了海量的生物序列数据。为了使用好这些海量数据，揭示出生物序列数据背后更深层次的结构、功能信息，产生了计算机化的生物序列分析方法。传统的计算机化的生物序列分析方法的基本思想是当两个分子具有相似的序列时，它们很可能具有相似的三维结构和功能。因此，从数据库中浩瀚的生物序列资源里搜索目的生物序列的同源序列，寻找保守的生物序列模式成为传统生物序列分析的核心内容。根据同时进行比对的生物序列数目，序列比对分为双序列比对和多序列比对。序列比对也可分为全局比对和局部比对，全局比对考虑序列的全局相似性，局部比对考虑序列片段之间的相似性。
用于双序列比对的算法开发于19世纪70年代，开始于由Needleman和Wunsch提出的全局序列比对Needleman-Wunsch算法。Needleman-Wunsch算法的主要思想是利用动态规划的方法计算两条序列之间的最佳比对。根据动态规划算法的复杂性分析，动态规划算法的运算速度与待检序列的长度和数据库大小密切相关。因此对于大规模生物序列的分析，Needleman-Wunsch算法的时间复杂性和空间复杂性都很高。至19世纪80年代，Smith和Waterman提出了局部序列比对的Smith-Waterman算法。Smith-Waterman算法与Needleman-Wunsch算法相似，也存在计算速度慢的问题。
目前，FASTA算法和BLAST算法是替代Smith-Waterman算法的两个流行的局部序列比对算法。与Smith-Waterman算法不同，虽然它们的计算速度比较快，但不能保证找到最佳比对，即灵敏度低。
虽然双序列比对是传统生物序列分析的基础，但对于成组序列，必须进行多序列比对。在目前发展的众多多序列比对方法中，最常用的是来自Da-Fei Feng和Russell Doolittle的Clustal算法。这种方法需要先对所有的序列计算两两比对的分值，然后从关系最近的一对序列开始，逐步加入其他序列。ClustalW是最常用的多序列比对程序之一。多序列比对方法同样不能同时保证计算的高速度和高灵敏度。
人们应用基于序列相似性的传统序列比对分析方法成功地预测出一部分基因的结构和功能。但是随着研究的深入，发现相当一部分功能相似的分子之间并不存在保守序列或共同的功能结构域，如RNA沉默抑制子。这就决定了无法利用基于序列相似性的传统方法来研究这些分子的结构和功能。并且，传统的计算机化生物序列分析是以牺牲灵敏度来换取速度提高的。在处理海量数据时，也不能同时满足高速度和高灵敏度。
发明内容
本发明的目的是提供一种高速度、高灵敏度和高准确性的生物序列的分析方法。
为了实现本发明的目的，本发明的一种生物序列的分析方法，其包括如下步骤：
(1)、计算生物序列的间隙谱：计算生物序列中字符之间的距离，分别统计字符之间相同距离的出现频率，构成间隙谱；
(2)、计算不同生物序列间的相似度；
(3)、推导不同生物序列的同源性或生物学功能：根据步骤(2)计算得到的相似度，如果相似性高，则推导这些生物序列之间可能具有同源性，或可能具有相似的生物学功能。
所述的步骤(1)后还包括如下步骤(1)′：①将间隙谱中的频率数据进行归一化，得到归一化后的频率数据；②计算间隙谱中频率最大值、最小值、均值、中位值、方差中的一种或多种；③将步骤①②得到的数据依次排列组合成一个特征向量，表示一条生物序列，再用线性相关系数或距离法计算不同生物序列间的相似度。
所述的步骤(1)′还包括如下步骤统计归一化后的间隙谱中字符相同距离出现的频率值的高低或者差别，如果在不同生物序列的间隙谱中，某一字符对的某一频率值都较高，则这一字符对是这些生物序列的一种相似(保守)序列模式；如果这一频率值在不同生物序列中差别较大，则这一频率值对应的字符对是这些生物序列的一种差别序列模式。
所述的步骤②还包括计算间隙谱中出现频率最大值、最小值时字符之间的距离。
所述的计算生物序列中字符之间的距离的方法包括如下步骤：
在用一维坐标标识的长度为n的生物序列中，沿正链方向或反链方向找到第一次出现某特定字符的坐标；
沿着该方向找到第二次、第三次、直至第p(p≤n)次出现某特定字符坐标；
不重复计算任意两坐标间的坐标差，即可得到
个字符与字符之间的距离数据。
所述的统计字符之间相同距离的出现频率的方法包括如下步骤：
设两字符i与j之间的距离为dij，dij＝1，2，...，max；把两字符之间的距离dij从最小(即为1)到最大(即为max)进行升序或从最大(即为max)到最小(即为1)进行降序排列，分别统计相同距离的出现频率，得max维向量；对每个max维向量，选其S(S≤max)维子向量，构成间隙谱，再进行上述步骤(1)′。
所述的距离是指在用一维坐标标识的生物序列中两个字符之间的坐标差。
所述的计算生物序列中字符之间的距离包括计算生物序列中一种或多种相同字符之间的距离和/或一种或多种不同字符之间的距离。
所述的统计字符之间相同距离出现的频率包括统计一种或多种相同字符之间相同距离出现的频率和/或统计一种或多种不同字符之间相同距离出现的频率。
对于某些生物序列，可以只统计几种字符之间相同距离的出现频率，构成间隙谱。根据这些生物序列的间隙谱，①计算这些生物序列的相似(保守)序列模式和差别序列模式，推导可能与功能直接相关的序列信息，如不连续的保守序列和不连续的差别序列。或者计算这些生物序列的相似度，推导不同生物序列的同源性或生物学功能，如果相似性高，则推导这些生物序列之间可能具有同源性，或可能具有相似的生物学功能
本发明所述的归一化法可以采用本领域常用的各种归一化法，例如均值归一化法、长度归一化法等。
其中max维向量是指由相同距离的出现频率组成的max维向量。对于16种不同的两字符(即：A-A，A-C，A-G，A-T，C-A，C-C，C-G，C-T，G-A，G-C，G-G，G-T，T-A，T-C，T-G，T-T)，由于两字符之间距离的最大值可能不尽相同，统计后得到的相同距离出现频率组成的向量维数可能也不尽相同，所以用max来笼统表示向量的维数。
其中所述的生物序列对于DNA序列而言，可以是全基因组序列及完整基因序列，也可以是具有不同生物学意义的若干子序列，如内含子、外显子、启动子等。其中可以计算全部16种字符之间的距离(A-A，A-C，A-G，A-T，C-A，C-C，C-G，C-T，G-A，G-C，G-G，G-T，T-A，T-C，T-G，T-T)，也可以根据需要选择计算其中几种字符之间的距离。
本发明的一种基于间隙谱的生物序列分析方法的有益效果，这种方法首先针对生物序列建立一个比较有生物学意义的、比较容易理解的、有效的数学形式描述***。在此基础上，提取合适的生物序列特征，对生物序列进行高速度、高灵敏度和高准确性的序列分析，从而揭示出生物序列数据背后更深层次的功能信息。该方法解决了传统生物序列分析方法所无法克服的低相似性序列分析即灵敏度低的问题，序列分析***的性能显著提高。
附图说明
图1是两DNA序列字符A-A的间隙谱；
图2是两DNA序列字符G-G的间隙谱；
图3是12个冠状病毒的进化树；
图4是12个冠状病毒和12个SARS病毒的进化树；
图5是27条水稻cDNA序列的进化树。
具体实施方式
以下实施例用于说明本发明，但不用来限制本发明的范围。
实施例1 基于间隙谱的DNA序列分析
将全长为20的DNA序列(SEQ1)的子序列(seq1)进行一维坐标标识。
A A A A G G C C T A C A A T T A C T C T 全长为20的DNA序列(负链方向)
将此条全长为20的DNA序列(SEQ1)沿其负链方向进行一维坐标标识。坐标标识起始于第4个碱基(标识为1)，终止于倒数第5个碱基(标识为13)，最终得到具有一维坐标标识的长度为13的DNA子序列(seq1)。标识结果如下：
A G G C C T A C A A T T A 长度为13的DNA子序列
1 2 3 4 5 6 7 8 9 10 11 12 13 一维坐标标识的长度为13的DNA子序列
以字符A为例，沿着该子序列的坐标标识方向找到第二次、第三次、直至最后一次(第五次)出现字符A的坐标，如下所示：
A A A A A DNA子序列中的字符A
1 7 9 10 13 一维坐标标识
不重复计算任意两字符A之间的坐标差，得到字符A与字符A之间的距离数据，见表1。
表1 字符A与字符A之间的10个距离数据
| 第2个字符A | 第3个字符A | 第4个字符A | 第5个字符A | |
| 第1个字符A | 6 | 8 | 9 | 12 |
| 第2个字符A | 2 | 3 | 6 | |
| 第3个字符A | 1 | 4 | ||
| 第4个字符A | 3 |
同字符A与字符A之间的距离计算方法，不重复计算任意两坐标间的坐标差，得到子序列中所有相同字符之间的距离和不同字符之间的距离，结果见表2、3、4、5、6、7、8、9、10、11、12、13、14。
表2 字符A与字符C之间的4个距离数据
| 第1个字符C | 第2个字符C | 第3个字符C | |
| 第1个字符A | 3 | 4 | 7 |
| 第2个字符A | 1 |
表3 字符A与字符G之间的2个距离数据
| 第1个字符G | 第2个字符G | |
| 第1个字符A | 1 | 2 |
表4 字符A与字符T之间的9个距离数据
| 第1个字符T | 第2个字符T | 第3个字符T | |
| 第1个字符A | 5 | 10 | 11 |
| 第2个字符A | 4 | 5 | |
| 第3个字符A | 2 | 3 | |
| 第4个字符A | 1 | 2 |
表5 字符C与字符A之间的11个距离数据
| 第1个字符A | 第2个字符A | 第3个字符A | 第4个字符A | |
| 第1个字符C | 3 | 5 | 6 | 9 |
| 第2个字符C | 2 | 4 | 5 | 8 |
| 第3个字符C | 1 | 2 | 5 |
表6 字符C与字符C之间的3个距离数据
| 第2个字符C | 第3个字符C | |
| 第1个字符C | 1 | 4 |
| 第2个字符C | 3 |
表7 字符C与字符T之间的8个距离数据
| 第1个字符T | 第2个字符T | 第3个字符T | |
| 第1个字符C | 2 | 7 | 8 |
| 第2个字符C | 1 | 6 | 7 |
| 第3个字符C | 3 | 4 |
表8 字符G与字符A之间的8个距离数据
| 第1个字符A | 第2个字符A | 第3个字符A | 第4个字符A | |
| 第1个字符G | 5 | 7 | 8 | 11 |
| 第2个字符G | 4 | 6 | 7 | 10 |
表9 字符G与字符C之间的6个距离数据
| 第1个字符C | 第2个字符C | 第3个字符C | |
| 第1个字符G | 2 | 3 | 6 |
| 第2个字符G | 1 | 2 | 5 |
表10 字符G与字符G之间的1个距离数据
| 第2个字符G | |
| 第1个字符G | 1 |
表11 字符G与字符T之间的6个距离数据
| 第1个字符T | 第2个字符T | 第3个字符T | |
| 第1个字符G | 4 | 9 | 10 |
| 第2个字符G | 3 | 8 | 9 |
表12 字符T与字符A之间的6个距离数据
| 第1个字符A | 第2个字符A | 第3个字符A | 第4个字符A | |
| 第1个字符T | 1 | 3 | 4 | 7 |
| 第2个字符T | 2 | |||
| 第3个字符T | 1 |
表13 字符T与字符C之间的1个距离数据
| 第2个字符C | |
| 第1个字符T | 2 |
表14 字符T与字符T之间的3个距离数据
| 第2个字符T | 第3个字符T | |
| 第1个字符T | 5 | 6 |
| 第2个字符T | 1 |
没有字符C与字符G之间、字符T与字符G之间的距离数据，因此两种距离用0表示。
设两字符i与j(其中i与j可以相同也可以不同)之间的距离为dij，dij＝1，2，...，max；把两字符之间的距离dij从最小(即为1)到最大(即为max)进行升序排列，分别统计相同距离的出现频率，得12维向量，见表15。
表15 两字符之间相同距离出现频率统计
对每个12维向量，都选取从第2维到第11维(共10维)的子向量，共得到16个10维向量，构成16个间隙谱，见表16。对每个10维向量进行如下所述的数据处理：
计算间隙谱中频率最大值、最小值、均值、中位值、方差，结果见表17。
表16 维数选取后两字符之间相同距离出现频率统计
表17 维数选取后的特征分析结果
| 频率统计值 | A-A | A-C | A-G | A-T | C-A | C-C | C-G | C-T | G-A | G-C | G-G | G-T | T-A | T-C | T-G | T-T |
| 最大值 | 2 | 1 | 1 | 2 | 3 | 1 | 0 | 2 | 2 | 2 | 0 | 2 | 1 | 1 | 0 | 1 |
| 最小值 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 均值 | 0.8 | 0.3 | 0.1 | 0.8 | 1 | 0.2 | 0 | 0.7 | 0.8 | 0.5 | 0 | 0.6 | 0.4 | 0.1 | 0 | 0.2 |
| 中位值 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0.5 | 0 | 0 | 0 | 0 |
| 方差 | 0.62 | 0.23 | 0.1 | 0.62 | 0.89 | 0.18 | 0 | 0.46 | 0.4 | 0.5 | 0 | 0.49 | 0.27 | 0.1 | 0 | 0.18 |
将间隙谱中的频率数据，采用均值归一化法进行归一化。
均值归一化的方法为：
其中，fij和
分别为未归一化和归一化后字符i与字符j之间相同距离出现的频率值；uij为字符i与字符j之间相同距离出现频率的均值。
经过均值归一化后所得到的频率数据见表18。
表18 采用均值归一化方法进行归一化后的频率数据
| 字符间距离 | A-A | A-C | A-G | A-T | C-A | C-C | C-G | C-T | G-A | G-C | G-G | G-T | T-A | T-C | T-G | T-T |
| 2 | 1.25 | 0.00 | 10.0 | 2.50 | 2.00 | 0.00 | 0.00 | 1.43 | 0.00 | 4.00 | 0.00 | 0.00 | 2.50 | 10.0 | 0.00 | 0.00 |
| 3 | 2.50 | 3.33 | 0.00 | 1.25 | 1.00 | 5.00 | 0.00 | 1.43 | 0.00 | 2.00 | 0.00 | 167 | 2.50 | 0.00 | 0.00 | 0.00 |
| 4 | 1.25 | 3.33 | 0.00 | 1.25 | 1.00 | 5.00 | 0.00 | 1.43 | 1.25 | 0.00 | 0.00 | 167 | 2.50 | 0.00 | 0.00 | 0.00 |
| 5 | 0.00 | 0.00 | 0.00 | 2.50 | 3.00 | 0.00 | 0.00 | 0.00 | 1.25 | 2.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 5.00 |
| 6 | 2.50 | 0.00 | 0.00 | 0.00 | 1.00 | 0.00 | 0.00 | 1.43 | 1.25 | 2.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 5.00 |
| 7 | 0.00 | 3.33 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 2.86 | 2.50 | 0.00 | 0.00 | 0.00 | 2.50 | 0.00 | 0.00 | 0.00 |
| 8 | 1.25 | 0.00 | 0.00 | 0.00 | 1.00 | 0.00 | 0.00 | 1.43 | 1.25 | 0.00 | 0.00 | 1.67 | 0.00 | 0.00 | 0.00 | 0.00 |
| 9 | 1.25 | 0.00 | 0.00 | 0.00 | 1.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 3.33 | 0.00 | 0.00 | 0.00 | 0.00 |
| 10 | 0.00 | 0.00 | 0.00 | 1.25 | 0.00 | 0.00 | 0.00 | 0.00 | 1.25 | 0.00 | 0.00 | 1.67 | 0.00 | 0.00 | 0.00 | 0.00 |
| 11 | 0.00 | 0.00 | 0.00 | 1.25 | 0.00 | 0.00 | 0.00 | 0.00 | 1.25 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
备注：表中的10.0为10.00
将间隙谱中采用均值归一化方法进行归一化后的频率数据(共160个频率值)以及间隙谱中频率最大值、均值、中位值、方差(共64个值)作为SEQ1的特征，按照A-A，A-C，A-G，A-T，C-A，C-C，C-G，C-T，G-A，G-C，G-G，G-T，T-A，T-C，T-G，T-T的次序组合成一个1×224的特征向量，这样SEQ1即可用这个特征向量表示，见表19。
表19 SEQ1的特征向量
特征向量
1-10 1.25 2.50 1.25 0.00 2.50 0.00 1.25 1.25 0.00 0.00
11-20 0.00 3.33 3.33 0.00 0.00 3.33 0.00 0.00 0.00 0.00
21-30 10.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
31-40 2.50 1.25 1.25 2.50 0.00 0.00 0.00 0.00 1.25 1.25
41-50 2.00 1.00 1.00 3.00 1.00 0.00 1.00 1.00 0.00 0.00
51-60 0.00 5.00 5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
61-70 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
71-80 1.43 1.43 1.43 0.00 1.43 2.86 1.43 0.00 0.00 0.00
81-90 0.00 0.00 1.25 1.25 1.25 2.50 1.25 0.00 1.25 1.25
91-100 4.00 2.00 0.00 2.00 2.00 0.00 0.00 0.00 0.00 0.00
101-110 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
111-120 0.00 1.67 1.67 0.00 0.00 0.00 1.67 3.33 1.67 0.00
121-130 2.50 2.50 2.50 0.00 0.00 2.50 0.00 0.00 0.00 0.00
131-140 10.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
141-150 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
151-160 0.00 0.00 0.00 5.00 5.00 0.00 0.00 0.00 0.00 0.00
161-170 2.00 0.80 1.00 0.62 1.00 0.30 0.00 0.23 1.00 0.10
171-180 0.00 0.10 2.00 0.80 1.00 0.62 3.00 1.00 1.00 0.89
181-190 1.00 0.20 0.00 0.18 0.00 0.00 0.00 0.00 2.00 0.70
191-200 1.00 0.46 2.00 0.80 1.00 0.40 2.00 0.50 0.00 0.50
201-210 0.00 0.00 0.00 0.00 2.00 0.60 0.50 0.49 1.00 0.40
211-220 0.00 0.27 1.00 0.10 0.00 0.10 0.00 0.00 0.00 0.00
221-224 1.00 0.20 0.00 0.18
根据上述方法，对另外一条全长为20的DNA序列(SEQ2)的子序列(seq2)进行分析，得到SEQ2的特征向量，见表20。
长度为16的子序列seq2为：AGCAACATGCGCGCTC
表20 SEQ2的特征向量
特征向量
1-10 2.00 4.00 2.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00
11-20 1.25 0.63 0.00 1.88 0.63 1.25 0.63 1.88 0.63 1.25
21-30 0.91 0.00 1.82 0.91 1.82 0.91 1.82 0.91 0.91 0.00
31-40 0.00 1.67 1.67 0.00 0.00 1.67 1.67 0.00 1.67 1.67
41-50 5.00 0.00 5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
51-60 2.14 0.71 2.14 0.00 1.43 0.71 0.71 0.71 0.71 0.71
61-70 0.00 2.86 0.00 1.43 1.43 1.43 1.43 0.00 1.43 0.00
71-80 2.00 2.00 0.00 4.00 0.00 0.00 0.00 2.00 0.00 0.00
81-90 3.33 3.33 0.00 3.33 0.00 0.00 0.00 0.00 0.00 0.00
91-100 0.00 3.33 1.11 2.22 0.00 1.11 1.11 0.00 1.11 0.00
101-110 3.33 0.00 1.67 0.00 0.00 1.67 0.00 1.67 0.00 1.67
111-120 2.50 0.00 2.50 0.00 5.00 0.00 0.00 0.00 0.00 0.00
121-130 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
131-140 2.50 0.00 2.50 0.00 2.50 0.00 2.50 0.00 0.00 0.00
141-150 0.00 5.00 0.00 5.00 0.00 0.00 0.00 0.00 0.00 0.00
151-160 0.00 0.00 0.00 0.00 0.00 10.00 0.00 0.00 0.00 0.00
161-170 2.00 0.50 0.00 0.50 3.00 1.60 1.50 0.93 2.00 1.10
171-180 1.00 0.54 1.00 0.60 1.00 0.27 1.00 0.20 0.00 0.18
181-190 3.00 1.40 1.00 0.93 2.00 0.70 1.00 0.46 2.00 0.50
191-200 0.00 0.50 1.00 0.30 0.00 0.23 3.00 0.90 1.00 0.99
201-210 2.00 0.60 0.50 0.49 2.00 0.40 0.00 0.49 0.00 0.00
211-220 0.00 0.00 1.00 0.40 0.00 0.27 1.00 0.20 0.00 0.18
221-224 1.00 0.10 0.00 0.10
利用线性相关系数计算用间隙谱表示的上述两条DNA序列的相似度，得到相似度为4％，说明这两条DNA序列很可能不具有同源性。
应用上述方法，可对大量的不同生物序列进行分析，推导不同生物序列的同源性或生物学功能。
在这两条DNA序列的间隙谱中，字符对A-A在距离1、2、3时频率值都较高，见图1，则字符对A-A在距离1、2、3时是这两条DNA序列的一种相似(保守)序列模式；字符对G-G的频率值在间隙谱中差别较大，见图2，则字符对G-G是这两条DNA序列的一种差别序列模式。
实施例2
本实施例是在实施例1的基础上，对两条DNA序列引入一个新的特征，即：针对实施例1中间隙谱(表16)的频率数据，计算间隙谱中对应频率最大时字符之间的距离，计算结果见表21。
表21 频率最大时字符之间的距离
| 频率统计值 | A-A | A-C | A-G | A-T | C-A | C-C | C-G | C-T | G-A | G-C | G-G | G-T | T-A | T-C | T-G | T-T |
| seq1中频率最大时字符间距 | 3，6 | 3，4，7 | 2 | 2，5 | 5 | 3，4 | / | 7 | 7 | 2 | / | 9 | 2，3，4，7 | 2 | / | 5，6 |
| seq2中频率最大时字符间距 | 3 | 5，9 | 4，6，8 | 3，4，7，8，10，11 | 2，4 | 2，4 | 3 | 5 | 2，3，5 | 3 | 2 | 6 | / | 2，4，6，8 | 3，5 | 7 |
将间隙谱中频率最大时字符A-A、C-A、G-C之间的距离也作为这两条DNA序列的特征，按照A-A，A-C，A-G，A-T，C-A，C-C，C-G，C-T，G-A，G-C，G-G，G-T，T-A，T-C，T-G，T-T的次序组合成1×228的特征向量，见表22、表23。
表22 SEQ1的特征向量
特征向量
1-10 1.25 2.50 1.25 0.00 2.50 0.00 1.25 1.25 0.00 0.00
11-20 0.00 3.33 3.33 0.00 0.00 3.33 0.00 0.00 0.00 0.00
21-30 10.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
31-40 2.50 1.25 1.25 2.50 0.00 0.00 0.00 0.00 1.25 1.25
41-50 2.00 1.00 1.00 3.00 1.00 0.00 1.00 1.00 0.00 0.00
51-60 0.00 5.00 5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
61-70 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
71-80 1.43 1.43 1.43 0.00 1.43 2.86 1.43 0.00 0.00 0.00
81-90 0.00 0.00 1.25 1.25 1.25 2.50 1.25 0.00 1.25 1.25
91-100 4.00 2.00 0.00 2.00 2.00 0.00 0.00 0.00 0.00 0.00
101-110 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
111-120 0.00 1.67 1.67 0.00 0.00 0.00 1.67 3.33 1.67 0.00
121-130 2.50 2.50 2.50 0.00 0.00 2.50 0.00 0.00 0.00 0.00
131-140 10.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
141-150 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
151-160 0.00 0.00 0.00 5.00 5.00 0.00 0.00 0.00 0.00 0.00
161-170 2.00 0.80 1.00 0.62 1.00 0.30 0.00 0.23 1.00 0.10
171-180 0.00 0.10 2.00 0.80 1.00 0.62 3.00 1.00 1.00 0.89
181-190 1.00 0.20 0.00 0.18 0.00 0.00 0.00 0.00 2.00 0.70
191-200 1.00 0.46 2.00 0.80 1.00 0.40 2.00 0.50 0.00 0.50
201-210 0.00 0.00 0.00 0.00 2.00 0.60 0.50 0.49 1.00 0.40
211-220 0.00 0.27 1.00 0.10 0.00 0.10 0.00 0.00 0.00 0.00
221-228 1.00 0.20 0.00 0.18 3.00 6.00 5.00 2.00
表23 SEQ2的特征向量
特征向量
1-10 2.00 4.00 2.00 0.00 2.00 0.00 0.00 0.00 0.00 0.00
11-20 1.25 0.63 0.00 1.88 0.63 1.25 0.63 1.88 0.63 1.25
21-30 0.91 0.00 1.82 0.91 1.82 0.91 1.82 0.91 0.91 0.00
31-40 0.00 1.67 1.67 0.00 0.00 1.67 1.67 0.00 1.67 1.67
41-50 5.00 0.00 5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
51-60 2.14 0.71 2.14 0.00 1.43 0.71 0.71 0.71 0.71 0.71
61-70 0.00 2.86 0.00 1.43 1.43 1.43 1.43 0.00 1.43 0.00
71-80 2.00 2.00 0.00 4.00 0.00 0.00 0.00 2.00 0.00 0.00
81-90 3.33 3.33 0.00 3.33 0.00 0.00 0.00 0.00 0.00 0.00
91-100 0.00 3.33 1.11 2.22 0.00 1.11 1.11 0.00 1.11 0.00
101-110 3.33 0.00 1.67 0.00 0.00 1.67 0.00 1.67 0.00 1.67
111-120 2.50 0.00 2.50 0.00 5.00 0.00 0.00 0.00 0.00 0.00
121-130 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
131-140 2.50 0.00 2.50 0.00 2.50 0.00 2.50 0.00 0.00 0.00
141-150 0.00 5.00 0.00 5.00 0.00 0.00 0.00 0.00 0.00 0.00
151-160 0.00 0.00 0.00 0.00 0.00 10.00 0.00 0.00 0.00 0.00
161-170 2.00 0.50 0.00 0.50 3.00 1.60 1.50 0.93 2.00 1.10
171-180 1.00 0.54 1.00 0.60 1.00 0.27 1.00 0.20 0.00 0.18
181-190 3.00 1.40 1.00 0.93 2.00 0.70 1.00 0.46 2.00 0.50
191-200 0.00 0.50 1.00 0.30 0.00 0.23 3.00 0.90 1.00 0.99
201-210 2.00 0.60 0.50 0.49 2.00 0.40 0.00 0.49 0.00 0.00
211-220 0.00 0.00 1.00 0.40 0.00 0.27 1.00 0.20 0.00 0.18
221-228 1.00 0.10 0.00 0.10 3.00 2.00 4.00 3.00
利用线性相关系数计算用间隙谱表示的上述两条DNA序列的相似度，得到相似度为9％，说明这两条DNA序列很可能不具有同源性。
试验例1 利用本发明的生物序列的分析方法研究SARS病毒的进化情况
我们从GenBank中选取了3类共12条有代表性的冠状病毒全基因组序列和12条典型的SARS病毒全基因组序列，见表24。应用实施例1所述的生物序列分析方法对这24条病毒全基因组序列进行了序列分析。
对于每条病毒全基因组序列，我们只对两相同字符(即：A-A，C-C，G-G，T-T)之间相同距离的出现频率进行统计，构成4个间隙谱。对于每条病毒全基因组序列的每个间隙谱，我们都选取字符间距离从1到50并进行均值归一化后的频率值作为序列特征。将选取的4×50个频率数据按照A-A、C-C、G-G、T-T的顺序依次组合成一个特征向量，表示相应的病毒全基因组序列。
通常认为，两序列之间距离越小，相似度越大。根据用特征向量表示的24条病毒全基因组序列，利用欧氏距离计算序列之间的相似度(距离)，见表25。并根据已知的生物学知识，推导未知病毒全基因组序列之间的同源性，结果用进化树表示，见图3、图4。
表24 12个冠状病毒和12个SARS病毒
| 病毒名称 | 类别 |
| cPig 2 | 1 |
| cPig 1 | 1 |
| cMurine 4 | 2 |
| cMurine 3 | 2 |
| cMurine 2 | 2 |
| cMurine 1 | 2 |
| cHuman | 1 |
| cBovine 4 | 2 |
| cBovine 3 | 2 |
| cBovine 2 | 2 |
| cBovine 1 | 2 |
| cAvian | 3 |
| SARS TW1 | 未知 |
| SARS SG5 | 未知 |
| SARS SG4 | 未知 |
| SARS SG3 | 未知 |
| SARS SG2 | 未知 |
| SARS SG1 | 未知 |
| SARS HK3 | 未知 |
| SARS HK2 | 未知 |
| SARS HK1 | 未知 |
| SARS TOR2 | 未知 |
| SARS BJ01 | 未知 |
| SARS Urbani | 未知 |
从图3中可以看到，利用基于间隙谱的生物序列分析方法可以由计算机自动完成3类共12个冠状病毒的分类，分类结果准确无误。从图4中我们可以得出这样的假说：从病毒进化的角度上看，SARS病毒很可能是从第1类cPig系冠状病毒进化而来的。经文献检索，这一假说也是实验生物学家提出的一种主流假说。
另外，在计算环境相同的情况下，该方法的执行速度较传统的生物序列比对方法如BLAST、C1ustalW快1倍以上。
试验例2
利用基于间隙谱的生物序列分析方法对基因按照功能进行分类
我们从GenBank中下载了27条水稻的cDNA序列，按照其编码蛋白行使的功能，将这些序列可分为3类，见表26。应用如实施例1所述的生物序列的分析方法，我们对这27条水稻的cDNA序列的相似度进行了序列分析。
对于每条水稻的cDNA序列，我们只对两相同字符(即A-A，C-C，G-G，T-T)之间相同距离的出现频率进行统计，构成4个间隙谱。对于每条水稻的cDNA序列的每个间隙谱，我们都选取字符间距离从51到100并进行均值归一化后的频率值作为序列特征。将选取的4×50个频率数据按照A-A、C-C、G-G、T-T的顺序依次组合成一个特征向量，表示相应的水稻cDNA序列。
根据用特征向量表示的27条水稻的cDNA序列的特征，利用欧氏距离计算序列之间的相似度，见表27。结果用进化树表示，见图5。
表26 27条水稻的cDNA序列
从图3中我们发现，应用基于间隙谱的生物序列分析方法，我们可以完全准确的将这27条水稻cDNA序列按照其对应的生物学功能进行计算机自动分类，这一技术将在很大程度上有助于未知功能的基因生物学功能预测。
Claims (10)
1、一种生物序列的分析方法，其特征在于，该方法包括如下步骤：
(1)、计算生物序列的间隙谱：计算生物序列中字符之间的距离，分别统计字符之间相同距离的出现频率，构成间隙谱；
(2)、计算不同生物序列间的相似度；
(3)、推导不同生物序列的同源性或生物学功能：根据步骤(2)计算得到的相似度，如果相似性高，则推导这些生物序列之间可能具有同源性，或可能具有相似的生物学功能。
2、如权利要求1所述的分析方法，其特征在于，所述的步骤(1)后还包括如下步骤：①将间隙谱中的频率数据进行归一化，得到归一化后的频率数据；②计算间隙谱中频率最大值、最小值、均值、中位值、方差中的一种或多种；③将步骤①②得到的数据依次排列组合成一个特征向量，表示一条生物序列，再用线性相关系数或距离法计算不同生物序列间的相似度。
3、如权利要求2所述的分析方法，其特征在于，所述的步骤(1)还包括统计归一化后的间隙谱中字符相同距离出现的频率值的高低或者差别，如果在不同生物序列的间隙谱中，某一字符对的某一频率值都较高，则这一字符对是这些生物序列的一种相似序列模式；如果这一频率值在不同生物序列中差别较大，则这一频率值对应的字符对是这些生物序列的一种差别序列模式。
4、如权利要求3所述的分析方法，其特征在于，所述的步骤②还包括计算间隙谱中出现频率最大值、最小值时字符之间的距离。
5、如权利要求1-4任一所述的分析方法，其特征在于，所述的计算生物序列中字符之间的距离的方法包括如下步骤：
在用一维坐标标识的长度为n的生物序列中，沿正链方向或反链方向找到第一次出现某特定字符的坐标；
沿着该方向找到第二次、第三次、直至第p(p≤n)次出现某特定字符的坐标；
不重复计算任意两坐标间的坐标差，即可得到个字符与字符之间的距离数据。
6、如权利要求1-5任一所述的分析方法，其特征在于，所述的统计字符之间相同距离的出现频率的方法包括如下步骤：
设两字符i与j之间的距离为dij，dij＝1，2，...，max；把两字符之间的距离dij进行升序或降序排列，分别统计相同距离的出现频率，得max维向量；对每个max维向量，选其S(S≤max)维子向量，构成间隙谱。
7、如权利要求1-6任一所述的分析方法，其特征在于，所述的距离是指在用一维坐标标识的生物序列中两个字符之间的坐标差。
8、如权利要求1-7任一所述的分析方法，其特征在于，所述的计算生物序列中字符之间的距离包括计算生物序列中一种或多种相同字符之间的距离和/或一种或多种不同字符之间的距离。
9、如权利要求8所述的分析方法，其特征在于，所述的统计字符之间相同距离出现的频率包括统计一种或多种相同字符之间相同距离出现的频率和/或统计一种或多种不同字符之间相同距离出现的频率。
10、权利要求1-9任一所述的分析方法在推导不同生物序列间的生物学关系的应用。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CNA2008100572005A CN101497924A (zh) | 2008-01-30 | 2008-01-30 | 一种基于间隙谱的生物序列分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CNA2008100572005A CN101497924A (zh) | 2008-01-30 | 2008-01-30 | 一种基于间隙谱的生物序列分析方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN101497924A true CN101497924A (zh) | 2009-08-05 |
Family
ID=40945195
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CNA2008100572005A Pending CN101497924A (zh) | 2008-01-30 | 2008-01-30 | 一种基于间隙谱的生物序列分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN101497924A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101840467A (zh) * | 2010-04-20 | 2010-09-22 | 中国科学院研究生院 | 蛋白质组过滤进化分类方法及其*** |
| CN106529212A (zh) * | 2016-10-19 | 2017-03-22 | 哈尔滨工业大学深圳研究生院 | 基于序列依赖频率矩阵的生物序列进化信息提取方法及其应用 |
| CN108846262A (zh) * | 2018-05-31 | 2018-11-20 | 广西大学 | 基于dft的rna二级结构距离计算构建***发育树的方法 |
| CN109937426A (zh) * | 2016-04-11 | 2019-06-25 | 量子生物有限公司 | 用于生物数据管理的***和方法 |
-
2008
- 2008-01-30 CN CNA2008100572005A patent/CN101497924A/zh active Pending
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101840467A (zh) * | 2010-04-20 | 2010-09-22 | 中国科学院研究生院 | 蛋白质组过滤进化分类方法及其*** |
| CN101840467B (zh) * | 2010-04-20 | 2012-07-04 | 中国科学院研究生院 | 蛋白质组过滤进化分类方法及其*** |
| CN109937426A (zh) * | 2016-04-11 | 2019-06-25 | 量子生物有限公司 | 用于生物数据管理的***和方法 |
| CN106529212A (zh) * | 2016-10-19 | 2017-03-22 | 哈尔滨工业大学深圳研究生院 | 基于序列依赖频率矩阵的生物序列进化信息提取方法及其应用 |
| CN106529212B (zh) * | 2016-10-19 | 2019-01-25 | 哈尔滨工业大学深圳研究生院 | 基于序列依赖频率矩阵的生物序列进化信息提取方法 |
| CN108846262A (zh) * | 2018-05-31 | 2018-11-20 | 广西大学 | 基于dft的rna二级结构距离计算构建***发育树的方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Lücking et al. | Multiple ITS haplotypes in the genome of the lichenized basidiomycete Cora inversa (Hygrophoraceae): fact or artifact? | |
| CN101497924A (zh) | 一种基于间隙谱的生物序列分析方法 | |
| CN106228035A (zh) | 基于局部敏感哈希和非参数化贝叶斯方法的高效聚类方法 | |
| WO2017210102A1 (en) | Methods and system for generating and comparing reduced genome data sets | |
| Boden et al. | Alignment-free sequence comparison with spaced k-mers | |
| US20150310168A1 (en) | Method for predicting gene cluster including secondary metabolism-related genes, prediction program, and prediction device | |
| Lin et al. | SSAW: A new sequence similarity analysis method based on the stationary discrete wavelet transform | |
| Rasheed et al. | LSH-Div: Species diversity estimation using locality sensitive hashing | |
| Wei et al. | A new classification method for human gene splice site prediction | |
| CN106557668B (zh) | 基于lf熵的dna序列相似性检验方法 | |
| Endo | Probabilistic nucleotide assembling method for sequencing by hybridization | |
| Xie et al. | Similarity evaluation of DNA sequences based on frequent patterns and entropy | |
| Maier et al. | Entropy predicts sensitivity of pseudorandom seeds | |
| Khaled et al. | Accelerating pairwise DNA Sequence Alignment using the CUDA compatible GPU | |
| Nguyen et al. | Efficient and accurate OTU clustering with GPU-based sequence alignment and dynamic dendrogram cutting | |
| EP2390811B1 (en) | Identification of ribosomal DNA sequences | |
| Böer | Multiple alignment using hidden Markov models | |
| CN101320404B (zh) | 一种生物病毒的计算机自动分类方法 | |
| Weitschek et al. | Classifying bacterial genomes with compact logic formulas on k-Mer frequencies | |
| Bandyopadhyay et al. | A parallel pairwise local sequence alignment algorithm | |
| Das et al. | A novel SFLA based method for gene expression biclustering | |
| CN114334004B (zh) | 一种病原微生物快速比对鉴定方法及其应用 | |
| Song et al. | A new graph theoretic approach for protein threading | |
| Sun et al. | Genome-scale NCRNA homology search using a Hamming distance-based filtration strategy | |
| Kakazu et al. | GPU acceleration for availability scoring of short constituent amino acid sequences |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C12 | Rejection of a patent application after its publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20090805 |