CN102750461A - 一种可得到完全解的生物序列局部比对方法 - Google Patents
一种可得到完全解的生物序列局部比对方法 Download PDFInfo
- Publication number
- CN102750461A CN102750461A CN2012101966689A CN201210196668A CN102750461A CN 102750461 A CN102750461 A CN 102750461A CN 2012101966689 A CN2012101966689 A CN 2012101966689A CN 201210196668 A CN201210196668 A CN 201210196668A CN 102750461 A CN102750461 A CN 102750461A
- Authority
- CN
- China
- Prior art keywords
- sequence
- score
- comparison
- result
- biological
- 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.)
- Granted
Links
Images
Landscapes
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
Abstract
一种可得到完全解的生物序列局部比对方法，包含以下步骤：步骤1：采用一种生物序列作为基准序列，另一种生物序列作查询序列，设定匹配得分Sa，不匹配得分Sb，起始罚分Sg，扩展罚分Ss，分数阈值H；步骤2：进行基准序列的后缀树分支与查询序列的比对，步骤如下：步骤3：整合各分支比对得分结果，取最大值作为两个生物序列的最终比对得分结果。步骤4：根据最终比对得分结果，寻找查询序列和基准序列中具有相似功能的片段或判断查询序列和基准序列之间的同源性关系。本发明采用BWT索引，结合过滤和重用技术，进行基准序列的后缀树分支与查询序列的比对，得出生物序列比对的完全解，弥补现有方法准确度不够或效率低下的问题。
Description
技术领域
本发明属于数据库和生物信息学领域，具体涉及一种可得到完全解的生物序列局部比对方法。
背景技术
在生物信息学研究中，经常需要将获得的基因或蛋白质序列（设为P）与已知的生物序列（设为T）进行比对。在很多时候，T和P从整体来看也许并不相似，然而二者却可能包含非常相似的子序列。局部比对的目的就是要找出这类具有高度相似性的子序列。局部比对技术在生物信息学研究中有重要的应用，例如可用于基因和蛋白功能研究、物种同源性研究等。将两条不同的基因序列进行局部比对，通过分析二者相似的子序列，从而找出两条基因序列中具有相似功能的基因片段。通过把新发现的蛋白序列同功能已知的蛋白序列作比对，则可以推测新蛋白的功能，指导新药的开发。在不同的物种中，特定基因的序列变异可用于研究物种之间的同源性。将两个物种的基因进行局部比对时，错配与突变相应，空位与插入或缺失对应，比对的结果可用于判断基因的相似度，此外也可以在基因组层面比较序列的相同与差异之处，将结果用于构建进化树。因此，探索生物序列的局部相似性具有非常重要的意义。如何能够准确快速的进行局部比对人们提出了挑战。目前已有的较为经典的算法包括Smith-Waterman、FASTA、BLAST等。
Smith-Waterman基于动态规划的思想，算法考虑两个序列中任意长度的子串，在计算得分的过程中允许匹配、不匹配和插入空格的操作。用这种方法得到的比对也是局部比对中得分最高者。该方法的时间和空间复杂度都是O（mn），其中m和n分别是P和T的长度。尽管该方法能够找到所有符合条件的结果，但是时间和空间上消耗都太大，以至于很少在实际中应用。
FASTA是一种经典的、基于启发式算法的生物序列局部比对工具，其基本思路是首先在T中精确匹配很短的序列片段K-tuple，之后采用启发式算法将构成的动态规划矩阵中同一或相邻对角线中位置相近的片段连接起来，构成局部比对结果。该方法虽有较高的效率，但并不能保证得到所有符合条件的结果。
BLAST是一种经典且非常流行的生物序列局部比对工具。与FASTA相似，BLAST也是基于启发式算法。首先在T和P中定位匹配程度超过一定阈值的短片段对segment pair，然后从这些位置开始向左右扩展得到满足给定阈值的局部比对结果。该方法虽然具有很高的效率，但同样不能保证找到所有符合条件的结果。
发明内容
针对现有技术存在的不足，本发明提供一种可得到完全解的生物序列局部比对方法，利用BWT索引，结合过滤和重用技术，得到生物序列局部比对的完全解。
本发明采用一组广泛应用的得分模式，在这组得分模式中，匹配得分是指若两个对应字符相同则为一个匹配（match），每一个匹配（match）得Sa分，不匹配得分是指若两个对应字符不相同，即需进行替换操作，则为一个不匹配（mismatch），不匹配得Sb分，若需进行插入或删除操作，则插入一个gap（连续插入r个空格）得分为Sg+r×Ss，其中Sg是gap起始罚分（gap opening penalty），即每插入一个gap需罚相应的分数，Ss是gap扩展罚分（gapextension penalty），也就是每插入一个空格罚Ss分。Sa为正分，Sb，Sg和Ss均为负分。
本发明方法包含以下步骤：
步骤1：采用一种生物序列作为基准序列T，另一种生物序列作查询序列P；
步骤2：进行基准序列的后缀树分支与查询序列的比对，步骤如下：
步骤2.1：设定匹配得分Sa，不匹配得分Sb，起始罚分Sg，扩展罚分Ss，分数阈值H；
步骤2.2：对基准序列T的逆序列T-1构建BWT索引；
BWT最早应用在数据压缩方面，对一个字符串进行BWT变换后并不改变字符串中字符的值，只是将它们的位置进行了改变。本方法通过BWT索引模拟后缀树遍历，构建BWT索引，步骤如下：
步骤2.2.1：在T-1的末尾增加一特殊字符$，使该字符小于T-1序列中所有字符；
步骤2.2.2：对T-1的后缀数组按字典序进行排序；
步骤2.2.3：建立数组SA，使其代表排序后的后缀数组中第i个位置的子序列在T-1中出现的开头位置；
步骤2.2.4：经过上述BWT变换后得到的序列的第i个字符BWT[i]=T-1[SA[i]-1]，若SA[i]-1=0，则BWT[i]=$；
步骤2.3：按基准序列的后缀树分支进行局部比对，计算各分支比对得分结果；
设X代表T的后缀树中从根节点到叶子节点的任意一支的字符串，后缀树分支可以定义成后缀树中从根节点到叶子节点的一条路径；
计算所有满足如下条件的比对：
sim(X[1，i]，P[y，j])≥H(1≤i≤|X|)(1≤y≤j≤|P|)
其中sim(X[1，i]，P[y，j])代表X[1，i]与P[y，j]比对的分数，H代表给定的分数阈值。P[y，j]代表P中从y到j位置的子串；
按基准序列的后缀树分支进行局部比对，计算各分支比对得分结果，具体按如下步骤进行：
步骤2.3.1：过滤，包括长度过滤、分数过滤、区域过滤、前缀过滤；
具体如下；
1）长度过滤：比对过程中，计算基准序列的子序列与查询序列构成的矩阵时，限定在一定长度范围内，过滤掉不必要的计算；
当计算矩阵MX时，只有当i满足如下条件时才需要计算MX（i，j）和其后续值：
其中i，j分别为MX矩阵中的横纵坐标位置，MX为采用动态规划方法计算基准序列子序列X与查询序列P的比对得分时构成的矩阵，MX（i，j）代表MX矩阵中第i行，第j列的位置的值，其后续值是矩阵中通过该位置的值计算出来的值，
代表对
向上取整，
代表对
向下取整。
2）分数过滤：比对过程中，将分数小于一定值的比对得分结果及时过滤掉；
对于从T中任意位置πt(1≤πt≤n)开始的任意子串X[1，i]，如果MX（i，j）满足如下条件则不需计算MX（i，j）及其后续值：
m代表P的长度，n代表T的长度；
3）区域过滤：将比对过程中的计算限定在以完全匹配开始的区域中；
可能作为比对结果的分数分布在特定的叉子型的区域中，因此该区域外的数据不需要计算。如图2所示，EMR代表完全匹配区域，即在叉子开头须有q长的子串与X中从1到q位置构成的子串X[1，q]完全匹配，其中
NGR代表无gap区域，在该区域中只有匹配与不匹配，不能插入gap。GR代表有gap区域，该区域中可能插入gap。FGOE为第一次插入gap的位置，也是NGR和GR的交界点，当满足score>|Sg+Ss|时出现。Extension Entry为叉子中可能允许出现gap的边界位置，当score>|Ss|时出现。
4）前缀过滤：若存在前缀支配关系，则直接过滤掉矩阵中对应的以完全匹配开始的区域；
若T的后缀树中的一条路径的前缀与另一条路径的前缀在T中相邻出现，则这两个前缀具有支配关系；
通过本发明提出的前缀支配关系，矩阵中的某些特定的叉子型区域可被直接过滤掉。设p和p1，…，pk为T的后缀树中从根节点到叶子节点的路径。使Xq和X1 q，…，Xk q分别代表p和p1，…，pk的q长前缀。若对于每个出现在t位置的Xq，总能在{X1 q，…，Xk q}中找到Xi q出现在t-1位置，则表示Xi q支配Xq。因此首先对T构建q长子串支配关系索引，可以在常量时间内判断某两个q长子串是否具有支配关系。如图2所示，在计算叉子型区域之前，首先检查P[πp-1，πp+q-2]是否在T中支配P[πp，πp+q-1]，若存在这种支配关系，则该叉子型区域可不需计算。
步骤2.3.2：计算结果重用；
由于P存在相同的子序列，因此在计算过程中一部分计算结果是相同的，通过本发明提出的计算结果重用技术可复重用这部分计算结果。如图3所示，可重用部分为叉子型区域间从FGOE开始的相同子串对应部分。设有两个叉子型区域f1和f2，并且他们的FGOE分别为（i，j1）和（i，j2）。设Ps为P[j1，m]和P[j2，m]的共同前缀，则对和i′≥i，MX(i′，j1+s)=MX(i′，j2+s)。在计算过程中在线建立共同前缀树用于匹配Ps，在计算叉子型区域时匹配该共同前缀树，匹配部分可直接复制之前计算的结果。
步骤3：整合各分支比对得分结果，取最大值作为两个生物序列的最终比对得分结果；
设p为T的后缀树中从根节点到叶子节点的一支，X代表p构成的子串，本发明需找到所有比对分数大于H的X[1，i]和P[y，j]。对于每一个X[1，i]对应T中t1，…，tk不同的位置。对于所有以T中t+i(t1≤t≤tk)和P中j结尾的子串，选择分数最大值作为最终比对结果。
计算完T的后缀树的一支后，保存对应局部矩阵的结果，即完成了P与T中某一子串的比对。当完成P与T中所有子串的比对后，将具有相同结尾位置的子串的结果进行整合，即取其中的最大值作为比对结果。
步骤4：根据最终比对得分结果，寻找查询序列和基准序列中具有相似功能的片段或判断查询序列和基准序列之间的同源性关系。
有益效果：本发明采用BWT索引，结合过滤和重用技术，能够得出生物序列比对的完全解，弥补了现有方法准确度不够或效率低下的问题，在保证找到全部解的前提下拥有很高的效率。对用于进行功能性研究的短查询和进行同源性研究的长查询均表现出了极高的查询效率。
附图说明
图1是本发明实施例一种可得到完全解的生物序列局部比对方法流程图；
图2是本发明实施例的区域过滤示意图；
图3是本发明实施例的计算结果重用示意图；
图4是本发明实施例的过滤计算后的局部矩阵示意图。
具体实施方式
下面结合附图对本发明具体实施方式做进一步说明。
本发明的实施方式使用GNU C++实现，测试硬件环境为Intel 2.93GHz四核i7处理器，8G内存，500G硬盘，装载Ubuntu（Linux）操作系统。
本实施方式的一种可得到完全解的生物序列局部比对方法流程如图1所示，包含以下步骤：
步骤1：采用一种生物序列作为基准序列T，另一种生物序列作查询序列P；
本发明的实施例1，采用两组DNA序列分别构成T和P，如下所示：
T=GCTAACTGCTAGCTGCGAGTTACC
P=GCTACCTGCTAGCTGCTAGCTGTG
步骤2：进行基准序列的后缀树分支与查询序列的比对，步骤如下：
步骤2.1：用户自行设定Sa=1Sb=-3Sg=-5Ss=-2H=7；
步骤2.2：对基准序列T的逆序列T-1构建BWT索引；
基准序列的逆序列T-1=CCATTGAGCGTCGATCGTCAATCG
通过BWT索引模拟后缀树遍历，构建BWT索引，步骤如下：
步骤2.2.1：在T-1的末尾增加一特殊字符$，使该字符小于T-1序列中所有字符，其形式如下：
CCATTGAGCGTCGATCGTCAATCG$
步骤2.2.2：对T-1的后缀数组按字典序进行排序；
其中字符串后对应的数字为对应的子串在T-1出现的开头位置。
步骤2.2.3：建立数组SA，使SA[i]代表排序后的后缀树中第i个位置的子序列在T-1中出现的开头位置；
SA[]={25，20，7，21，14，3，19，2，1，23，12，16，9，24，6，13，8，17，10，18，22，11，15，5，4}
步骤2.2.4：经过上述BWT变换后得到的序列的第i个字符为T-1中第SA[i]-1个字符，即BWT[i]=T-1[SA[i]-1]，若SA[i]-1=0，则BWT[i]=$；
构建后的索引如下：
GCGAGCTC$TTTGCTCACCGAGATA
本发明方法建立的BWT索引与后缀树相比极大的节省了空间，同时可实现如下本发明需要的模拟遍历后缀树的3个功能：
1）给定P中一个q长的短串Sq，检查Sq是否在T中出现。由于给定一个短串X的SA范围后，可以通过T的BWT索引在常量时间内确定在X之前添加一字符c构成的字符串cX的SA范围，若该SA范围合法则表示存在，反之表示不存在。同理，通过已经建立T-1的BWT索引可判断Xc对应的SA范围，从而确定Xc是否存在。所以循环此过程q次可判断Sq是否在T中出现。
2）给定一字符串X，找到X在T中出现的所有位置。利用T-1的BWT索引在判断X是否在T中出现的同时可得到X对应的SA范围，若X的SA范围为[x，y]，则X在T中出现的每个位置对应|T|-SA[i]-|X|+2（x≤i≤y）。
3）给定一字符串Xq，遍历T的后缀树中所有以Xq为前缀的分支。由于字符集已知，在得到Xq的SA范围后，在从字符集中顺序选取字符c构成Xqc判断是否存在，若存在继续遍历，反之则回溯。
步骤2.3：按基准序列的后缀树分支进行局部比对，计算各分支比对得分结果；
设X代表T的后缀树即T的子序列中从根节点到叶子节点的任意一支的字符串，本实施例以后缀树的一支为例，设
X=GCTAACTGCTAGCTGCGAGTTACC
计算所有满足如下条件的比对：
sim(X[1，i]，P[y，j])≥H(1≤i≤|X|)(1≤y≤j≤|P|)
其中sim(X[1，i]，P[y，j])代表X[1，i]与P[y，j]比对的分数，H代表给定的分数阈值。
按基准序列的后缀树分支进行局部比对，计算各分支比对得分结果，具体按如下步骤进行：
步骤2.3.1：过滤；
过滤包括长度过滤、分数过滤、区域过滤、前缀过滤，具体如下；
1）长度过滤；计算矩阵MX及其后续值时，限定在一定长度范围内，过滤掉不必要的计算，MX为采用动态规划方法计算基准序列子序列X与查询序列P的比对得分时构成的矩阵；
当计算MX时，只有当i满足如下条件时才需要计算MX（i，j）和其后续值：
其中i，j分别为MX矩阵中的横纵坐标位置。本例中，分数阈值H=7，匹配得分Sa=1，P的长度m=24，起始罚分Sg=-5，扩展罚分Ss=-2，计算
计算只需考虑T中长度在7和30之间的子串。
2）分数过滤；比对过程中，将分数小于一定值的比对得分结果及时过滤掉；
对于从T中πt(1≤πt≤n)位置开始的任意子串X[1，i]，本例中，πt=1，如果MX（i，j）满足如下条件则不需计算MX（i，j）及其后续值：
T的长度n=24；
如计算MX[5，19]，计算H-(m-j)×Sa-1=1，πt=1，计算
由于MX[5，19]≤1，因此MX[5，19]及其后续值均不需计算。
3）区域过滤：将比对过程中的计算限定在特定区域中；
可能作为比对结果的分数分布在特定的叉子型的区域中，如图2所示，因此区域外的数据不需要计算。代表完全匹配区域，即在叉子开头须有q长的子串与X[1，q]完全匹配，其中
计算得出q=4；NGR代表无gap区域，在该区域中只有匹配与不匹配，不能插入gap。GR代表有gap区域，该区域中可能插入gap。FGOE为第一次插入gap的位置，也是NGR和GR的交界点，当满足score>|Sg+Ss|时出现。Extension Entry为叉子中可能允许出现gap的边界位置，当score>|Ss|时出现。
EMR区域、NGR区域和GR区域共同构成叉子型区域；
如MX[1，1]为EMR区域的过滤，MX[6，7]为NGR区域的过滤，MX[12，14]为GR区域的过滤。
4）前缀过滤：若存在前缀支配关系，则直接过滤掉局部矩阵中的特定区域；
若T的后缀树中的一条路径的前缀与另一条路径的前缀在T中相邻出现，则这两个前缀具有支配关系；
通过本实施方式提出的前缀支配关系，矩阵中的某些特定的叉子型区域可被直接过滤掉。设p和p1，…，pk为T的后缀树中从根节点到叶子节点的路径。使Xq和X1 q，…，Xk q分别代表p和p1，…，pk的q长前缀。若对于每个出现在t位置的Xq，总能在{X1 q，…，Xk q}中找到Xi q出现在t-1位置，则表示Xi q支配Xq。因此首先对T构建q长子串支配关系索引，可以在常量时间内判断某两个q长子串是否具有支配关系。如图2所示，在计算叉子型区域之前，首先检查P[πp-1，πp+q-2]是否在T中支配P[πp，πp+q-1]，若存在这种支配关系，则该叉子型区域可不需计算。
本例中，考虑T的后缀树的另一支X′=CTAACTGCTAGCTGCGAGTTACC，，因为GCTA支配CTAG，因此在计算以P[9，12]或P[16，19]为对应EMR区域的叉子型区域时，可直接过滤掉不需计算。
步骤2.3.2：计算结果重用；
由于查询序列中存在相同的子序列，因此在比对过程中一部分计算结果是相同的，可重复使用这部分计算结果。
由于P存在相同的子串，因此在计算过程中一部分计算结果是相同的，通过本实施方式提出的计算结果重用技术可复重用这部分计算结果。如图3所示，可重用部分为叉子型区域间从FGOE开始的相同子串对应部分。设有两个叉子型区域f1和f2，并且他们的FGOE分别为（i，j1）和（1，j2）。设Ps为P[j1，m]和P[j2，m]的共同前缀，则对
和i′≥i，MX(i′，j1+s)=MX(i′，j2+s)。在计算过程中在线建立共同前缀树用于匹配Ps，在计算叉子型区域时匹配该共同前缀树，匹配部分可直接复制之前计算的结果。
如图4所示，P[13，15]与P[20，22]是相同的子串，P[20，22]=P[13，15]，因此在计算过程中的计算结果也是相同的，因此对应部分可直接复用之前计算的结果。
计算后的矩阵MX如图4所示，图中的FF为通过区域过滤过滤掉的分数，SF为通过分数过滤过滤掉的分数，RU为重用的分数。
重复上述过程对各分支进行局部比对，找到所有分数大于H，即大于7的结果。
步骤3：整合各分支结果，得出两个生物序列的最终比对结果；
对所有分数大于H，即大于7的结果进行整合，最终局部比对结果如下：
（1）Query:GCTACCTGCTA
Subject:GCTAACTGCTA
Score:7
（2）Query:GCTACCTGCTAG
Subject:GCTAACTGCTAG
Score:8
（3）Query:GCTACCTGCTAGC
Subject:GCTAACTGCTAGC
Score:9
（4）Query:GCTACCTGCTAGCT
Subject:GCTAACTGCTAGCT
Score:10
（5）Query:GCTACCTGCTAGCTG
Subject:GCTAACTGCTAGCTG
Score:11
（6）Query:GCTACCTGCTAGCTGC
Subject:GCTAACTGCTAGCTGC
Score:12
（7）Query:GCTACCTGCTAGCTGCT
Subject:GCTAACTGCTAGCTGCG
Score:9
（8）Query:GCTACCTGCTAGCTGCTA
Subject:GCTAACTGCTAGCTGCGA
Score:10
（9）Query:GCTACCTGCTAGCTGCTAG
Subject:GCTAACTGCTAGCTGCGAG
Score:11
（10）Query:GCTACCTGCTAGCTGCTAGC
Subject:GCTAACTGCTAGCTGCGAGT
Score:8
（11）Query:GCTACCTGCTAGCTGCTAGCT
Subject:GCTAACTGCTAGCTGCGAGTT
Score:9
（12）Query:GCTAGCTGCTA
Subject:GCTAACTGCTA
Score:7
（13）Query:GCTAGCTGCTAG
Subject:GCTAACTGCTAG
Score:8
（14）Query:GCTAGCTGCTAGC
Subject:GCTAACTGCTAGC
Score:9
（16）Query:GCTAGCTGCTAGCT
Subject:GCTAACTGCTAGCT
Score:10
（17）Query:GCTAGCTGCTAGCTG
Subject:GCTAACTGCTAGCTG
Score:11
其中Query后为P中子串，Subject为T中子串，Score为比对分数。选择分数最大值作为最终比对结果。
步骤4：根据最终比对得分结果，寻找查询序列和基准序列中具有相似功能的片段或判断查询序列和基准序列之间的同源性关系。
比对得分结果中，包含相似性大于用户给定分数阈值的子序列对，可用于功能性和同源性（homology）研究等许多方面。例如，通过比对得到的分数大于分数阈值的子序列对，可以判断查询序列和基准序列相应的基因片段可能具有相似的功能。通过分析两个序列中相似性大于用户给定分数阈值的子序列占整体的比例，或通过分析局部子序列的得分情况，可判断两个基因的同源性。一般来说，当相似程度高于50％时，比较容易推测查询序列和基准序列可能是同源序列；而当相似性程度低于20％时，就难以确定或者根本无法确定其是否具有同源性。其中相似程度可用序列比对过程中查询序列和基准序列之间相似性大于用户给定分数阈值的子序列占整体的比例。
实施例2：
基准序列T从人类基因序列（GRCh37）中提取，大小为1Gb，查询序列P从老鼠基因的第一条染色体（MGSCv37chr1）中提取，按不同长度，每种长度从随机位置提取100条。
对上述两个序列执行本发明方法，步骤如下：
步骤1：采用一种生物序列作为基准序列T，另一种生物序列作查询序列P；
由于数据量太大，因此取T的一个后缀X，即后缀树的一个分支作为实例来说明实施过程。
X=ATGCCTGATGCATGATACAGGCTT
P=ATGCTTGATGCATGATGCATGAGA
步骤2：进行基准序列的后缀树分支与查询序列的比对，步骤如下：
步骤2.1：用户自行设定Sa=1Sb=-3Sg=-5Ss=-2H=7；
步骤2.2：对基准序列T的逆序列T-1构建BWT索引；
同样由于为T-1构建的BWT索引太大，因此省略该步骤的描述，具体执行过程参见实施例1。
步骤2.3：按基准序列的后缀树分支进行局部比对，计算各分支比对得分结果；
设X代表T的后缀树即T的子序列中从根节点到叶子节点的任意一支的字符串，本实施例以后缀树的一支为例，设
X=ATGCCTGATGCATGATACAGGCTT
计算所有满足如下条件的比对：
sim(X[1，i]，P[y，j])≥H(1≤i≤|X|)(1≤y≤j≤|P|)
其中sim(X[1，i]，P[y，j])代表X[1，i]与P[y，j]比对的分数，H代表给定的分数阈值。
按基准序列的后缀树分支进行局部比对，计算各分支比对得分结果，具体按如下步骤进行：
步骤2.3.1：过滤；
过滤包括长度过滤、分数过滤、区域过滤、前缀过滤，具体如下；
1）长度过滤；计算矩阵MX及其后续值时，限定在一定长度范围内，过滤掉不必要的计算，MX为采用动态规划方法计算基准序列子序列X与查询序列P的比对得分时构成的矩阵；
当计算MX时，只有当i满足如下条件时才需要计算MX（i，j）和其后续值：
其中i，j分别为MX矩阵中的横纵坐标位置。本例中，分数阈值H=7，匹配得分Sa=1，P的长度m=24，起始罚分Sg=-5，扩展罚分Ss=-2，计算
计算只需考虑T中长度在7和30之间的子串。
2）分数过滤；比对过程中，将分数小于一定值的比对得分结果及时过滤掉；
对于从T中πt(1≤πt≤n)位置开始的任意子串X[1，i]，本例中，πt=1，如果MX（i，j）满足如下条件则不需计算MX（i，j）及其后续值：
T的长度n=24；
如计算MX[22，21]，计算H-(m-j)×Sa-1=3，πt=1，计算
由于
MX[22，21]≤3，因此MX[22，21]及其后续值均不需计算。
3）区域过滤：将比对过程中的计算限定在特定区域中；
可能作为比对结果的分数分布在特定的叉子型的区域中，因此区域外的数据不需要计算。代表完全匹配区域，即在叉子开头须有q长的子串与X[1，q]完全匹配，其中
计算得出q=4；NGR代表无gap区域，在该区域中只有匹配与不匹配，不能插入gap。GR代表有gap区域，该区域中可能插入gap。FGOE为第一次插入gap的位置，也是NGR和GR的交界点，当满足score>|Sg+Ss|时出现。Extension Entry为叉子中可能允许出现gap的边界位置，当score>|Ss|时出现。
EMR区域、NGR区域和GR区域共同构成叉子型区域；
如MX[4，18]为EMR区域的过滤，MX[7，15]为NGR区域的过滤，MX[13，19]为GR区域的过滤。
4）前缀过滤：若存在前缀支配关系，则直接过滤掉局部矩阵中的特定区域；
若T的后缀树中的一条路径的前缀与另一条路径的前缀在T中相邻出现，则这两个前缀具有支配关系；
通过本实施方式提出的前缀支配关系，矩阵中的某些特定的叉子型区域可被直接过滤掉。设p和p1，…，pk为T的后缀树中从根节点到叶子节点的路径。使Xq和X1 q，…，Xk q分别代表p和p1，…，pk的q长前缀。若对于每个出现在t位置的Xq，总能在{X1 q，…，Xk q}中找到Xi q出现在t-1位置，则表示Xi q支配Xq。因此首先对T构建q长子串支配关系索引，可以在常量时间内判断某两个q长子串是否具有支配关系。如图2所示，在计算叉子型区域之前，首先检查P[πp-1，πp+q-2]是否在T中支配P[πp，πp+q-1]，若存在这种支配关系，则该叉子型区域可不需计算。
本例中，考虑T的后缀树的另一支X′=TGCCTGATGCATGATACAGGCTT，因为ATGC支配TGCA，因此在计算以P[9，12]或P[16，19]为对应EMR区域的叉子型区域时，可直接过滤掉不需计算。
步骤2.3.2：计算结果重用；
由于查询序列中存在相同的子序列，因此在比对过程中一部分计算结果是相同的，可重复使用这部分计算结果。
由于P存在相同的子串，因此在计算过程中一部分计算结果是相同的，通过本实施方式提出的计算结果重用技术可复重用这部分计算结果。如图3所示，可重用部分为叉子型区域间从FGOE开始的相同子串对应部分。设有两个叉子型区域f1和f2，并且他们的FGOE分别为（i，j1）和（1，j2）。设Ps为P[j1，m]和P[j2，m]的共同前缀，则对
和i′≥i，MX(i′，j1+s)=MX(i′，j2+s)。在计算过程中在线建立共同前缀树用于匹配Ps，在计算叉子型区域时匹配该共同前缀树，匹配部分可直接复制之前计算的结果。
P[13，15]与P[20，22]是相同的子串，P[20，22]=P[13，15]，因此在计算过程中的计算结果也是相同的，因此对应部分可直接复用之前计算的结果。
重复上述过程对各分支进行局部比对，找到所有分数大于H，即大于7的结果。
步骤3：整合各分支结果，得出两个生物序列的最终比对结果；
对所有分数大于H，即大于7的结果进行整合，最终局部比对结果如下：
（1）Query:ATGCTTGATGC
Subject:ATGCCTGATGC
Score:7
（2）Query:ATGCTTGATGCA
Subject:ATGCCTGATGCA
Score:8
（3）Query:ATGCTTGATGCAT
Subject:ATGCCTGATGCAT
Score:9
（4）Query:ATGCTTGATGCATG
Subject:ATGCCTGATGCATG
Score:10
（5）Query:ATGCTTGATGCATGA
Subject:ATGCCTGATGCATGA
Score:11
（6）Query:ATGCTTGATGCATGAT
Subject:ATGCCTGATGCATGAT
Score:12
（7）Query:ATGCTTGATGCATGATG
Subject:ATGCCTGATGCATGATA
Score:9
（8）Query:ATGCTTGATGCATGATGC
Subject:ATGCCTGATGCATGATAC
Score:10
（9）Query:ATGCTTGATGCATGATGCA
Subject:ATGCCTGATGCATGATACA
Score:11
（10）Query:ATGCTTGATGCATGATGCAT
Subject:ATGCCTGATGCATGATACAG
Score:8
（11）Query:ATGCTTGATGCATGATGCATG
Subject:ATGCCTGATGCATGATACAGG
Score:9
（12）Query:ATGCATGATGC
Subject:ATGCCTGATGC
Score:7
（13）Query:ATGCATGATGCA
Subject:ATGCCTGATGCA
Score:8
（14）Query:ATGCATGATGCAT
Subject:ATGCCTGATGCAT
Score:9
（16）Query:ATGCATGATGCATG
Subject:ATGCCTGATGCATG
Score:10
（17）Query:ATGCATGATGCATGA
Subject:ATGCCTGATGCATGA
Score:11
其中Query后为P中子串，Subject为T中子串，Score为比对分数。选择分数最大值作为最终比对结果。
采用本发明的方法，不同长度的P得到完全解的平均时间下表所示：
其中ALAE为采用本实施方式的方法，BLAST为目前流行的局部比对工具，m为查询序列P的长度。表中列举了查询长度从1kb到1Mb变化时两种方法的平均比对时间。可以看出，ALAE总体优于BLAST。需要注意的是，BLAST不能保证找到完全解，而ALAE能够找到完全解。
实施例3：
本实施例采用三种链霉菌（Streptomyces）的基因组进行局部比对，分别为天蓝色链霉菌（S.coelicolor）基因组，全长8,667,507bp；阿维链霉菌（S.avermitilis）基因组，长度为9.02Mb；灰色链霉菌（S.griseus）的线性染色体，全长为8,545,929bp。由于数据量太大，因此取计算过程中的一个小片段作为实例说明。
步骤1到步骤2与实施例1和实施例2类似，因此不做赘述。
取计算过程中的小片段作为实例说明，即
X=TGACCGATGACTGATGTCTAACGG
P=TGACGGATGACTGATGACTGATAT
步骤3：整合各分支结果，得出两个生物序列的最终比对结果；
（1）Query:TGACGGATGAC
Subject:TGACCGATGAC
Score:7
（2）Query:TGACGGATGACT
Subject:TGACCGATGACT
Score:8
（3）Query:TGACGGATGACTG
Subject:TGACCGATGACTG
Score:9
（4）Query:TGACGGATGACTGA
Subject:TGACCGATGACTGA
Score:10
（5）Query:TGACGGATGACTGAT
Subject:TGACCGATGACTGAT
Score:11
（6）Query:TGACGGATGACTGATG
Subject:TGACCGATGACTGATG
Score:12
（7）Query:TGACGGATGACTGATGA
Subject:TGACCGATGACTGATGT
Score:9
（8）Query:TGACGGATGACTGATGAC
Subject:TGACCGATGACTGATGTC
Score:10
（9）Query:TGACGGATGACTGATGACT
Subject:TGACCGATGACTGATGTCT
Score:11
（10）Query:TGACGGATGACTGATGACTG
Subject:TGACCGATGACTGATGTCTA
Score:8
（11）Query:TGACGGATGACTGATGACTGA
Subject:TGACCGATGACTGATGTCTAA
Score:9
（12）Query:TGACTGATGAC
Subject:TGACCGATGAC
Score:7
（13）Query:TGACTGATGACT
Subject:TGACCGATGACT
Score:8
（14）Query:TGACTGATGACTG
Subject:TGACCGATGACTG
Score:9
（16）Query:TGACTGATGACTGA
Subject:TGACCGATGACTGA
Score:10
（17）Query:TGACTGATGACTGAT
Subject:TGACCGATGACTGAT
Score:11
步骤4：根据最终比对得分结果，寻找查询序列和基准序列中具有相似功能的片段，或判断查询序列和基准序列之间的同源性关系。
经过对局部比对结果的分析可知，天蓝色链霉菌、阿维链霉菌和灰色链霉菌DNA序列相似性高达80％以上，而且对应片段编码的蛋白功能也相同或相近，证实三种链霉菌之间在DNA序列水平上具有很高的同源性。
Claims (3)
1.一种可得到完全解的生物序列局部比对方法，其特征在于：包含以下步骤：
步骤1：采用一种生物序列作为基准序列，另一种生物序列作查询序列；
步骤2：进行基准序列的后缀树分支与查询序列的比对，步骤如下：
步骤2.1：设定匹配得分Sa，不匹配得分Sb，起始罚分Sg，扩展罚分Ss，分数阈值H；
步骤2.2：对基准序列的逆序列T-1构建BWT索引；
步骤2.3：按基准序列的后缀树分支进行局部比对，计算各分支比对得分结果；
步骤3：整合各分支比对得分结果，取最大值作为两个生物序列的最终比对得分结果；
步骤4：根据最终比对得分结果，寻找查询序列和基准序列中具有相似功能的片段或判断查询序列和基准序列之间的同源性关系。
2.根据权利要求1所述的可得到完全解的生物序列局部比对方法，其特征在于：步骤2.2对基准序列的逆序列T-1构建BWT索引，包括如下步骤：
步骤2.2.1：在基准序列的逆序列的末尾增加一字符，使该字符小于基准序列的逆序列中所有字符；
步骤2.2.2：对基准序列的逆序列的后缀树分支按字典序进行排序；
步骤2.2.3：建立数组SA，使其代表排序后的后缀数组中第i个位置的子序列在T-1中出现的开头位置；
步骤2.2.4：经过BWT变换后得到的序列的第i个字符为T-1中第SA[i]-1个字符。
3.根据权利要求1所述的可得到完全解的生物序列局部比对方法，其特征在于：所述步骤2.3按基准序列的后缀树分支进行局部比对，计算各分支比对得分结果，按如下步骤进行：
步骤2.3.1：过滤；
过滤包括长度过滤、分数过滤、区域过滤、前缀过滤，具体如下：
1)长度过滤：比对过程中，只有当i满足如下条件时才需要计算基准序列子序列X与查询序列P的比对得分构成的矩阵MX(i，j)和MX(i，j)后续值：
其中i，j分别为MX矩阵中的横纵坐标位置，m为P的长度；
2)分数过滤：比对过程中，对于从T中πt(1≤πt≤n)位置开始的任意子串X[1，i]，如果MX(i，j)满足如下条件则不需计算MX(i，j)及MX(i，j)后续值：
其中，n为T的长度；
3)区域过滤：将比对过程中的计算限定在以完全匹配开始的区域中；
4)前缀过滤：若存在前缀支配关系，则直接过滤掉矩阵中对应的以完全匹配开始的区域；
所述的支配关系是指若T的后缀树中的一条路径的前缀与另一条路径的前缀在T中相邻出现，则这两个前缀具有支配关系；
步骤2.3.2：计算结果重用；
由于查询序列中存在相同的子序列，因此在比对过程中一部分计算结果是相同的，可重复使用这部分计算结果。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201210196668.9A CN102750461B (zh) | 2012-06-14 | 2012-06-14 | 一种可得到完全解的生物序列局部比对方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201210196668.9A CN102750461B (zh) | 2012-06-14 | 2012-06-14 | 一种可得到完全解的生物序列局部比对方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN102750461A true CN102750461A (zh) | 2012-10-24 |
| CN102750461B CN102750461B (zh) | 2015-04-22 |
Family
ID=47030637
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201210196668.9A Active CN102750461B (zh) | 2012-06-14 | 2012-06-14 | 一种可得到完全解的生物序列局部比对方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN102750461B (zh) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103793626A (zh) * | 2012-10-29 | 2014-05-14 | 三星Sds株式会社 | 碱基序列比对系统及方法 |
| CN103902599A (zh) * | 2012-12-27 | 2014-07-02 | 北京新媒传信科技有限公司 | 模糊查找的方法和装置 |
| CN104395900A (zh) * | 2013-03-15 | 2015-03-04 | 北京未名博思生物智能科技开发有限公司 | 序列比对的空间计数运算方法 |
| CN104881592A (zh) * | 2015-02-11 | 2015-09-02 | 哈尔滨工业大学深圳研究生院 | 一种dna序列比对中的打分方法 |
| CN104899476A (zh) * | 2015-06-15 | 2015-09-09 | 中国人民解放军国防科学技术大学 | 一种对多序列bwt索引构建进行并行加速的方法 |
| CN105051741A (zh) * | 2012-12-17 | 2015-11-11 | 微软技术许可有限责任公司 | 并行局部序列对齐 |
| CN105956417A (zh) * | 2016-05-04 | 2016-09-21 | 西安电子科技大学 | 云环境下基于编辑距离的相似碱基序列查询方法 |
| CN106529212A (zh) * | 2016-10-19 | 2017-03-22 | 哈尔滨工业大学深圳研究生院 | 基于序列依赖频率矩阵的生物序列进化信息提取方法及其应用 |
| CN111445952A (zh) * | 2020-03-25 | 2020-07-24 | 山东大学 | 超长基因序列的相似性快速比对方法及系统 |
| CN111916153A (zh) * | 2020-06-17 | 2020-11-10 | 电子科技大学 | 一种并行多重序列比对方法 |
| CN112614544A (zh) * | 2020-12-28 | 2021-04-06 | 杭州瑞普基因科技有限公司 | Kraken2软件输出结果的优化方法及鉴定样本中物种类型的方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101158952A (zh) * | 2007-11-22 | 2008-04-09 | 中国人民解放军国防科学技术大学 | 基于流处理的生物序列数据库搜索多层次加速方法 |
| CN101984445A (zh) * | 2010-03-04 | 2011-03-09 | 深圳华大基因科技有限公司 | 一种基于聚合酶链式反应产物测序序列分型的实现方法和系统 |
-
2012
- 2012-06-14 CN CN201210196668.9A patent/CN102750461B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101158952A (zh) * | 2007-11-22 | 2008-04-09 | 中国人民解放军国防科学技术大学 | 基于流处理的生物序列数据库搜索多层次加速方法 |
| CN101984445A (zh) * | 2010-03-04 | 2011-03-09 | 深圳华大基因科技有限公司 | 一种基于聚合酶链式反应产物测序序列分型的实现方法和系统 |
Non-Patent Citations (1)
| Title |
|---|
| 熊文林: "面向基因组重测序的BWT索引压缩算法", 《中国优秀硕士论文全文数据库（信息科技辑）2012年》, no. 5, 15 May 2012 (2012-05-15) * |
Cited By (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103793626A (zh) * | 2012-10-29 | 2014-05-14 | 三星Sds株式会社 | 碱基序列比对系统及方法 |
| CN105051741A (zh) * | 2012-12-17 | 2015-11-11 | 微软技术许可有限责任公司 | 并行局部序列对齐 |
| CN105051741B (zh) * | 2012-12-17 | 2018-04-10 | 微软技术许可有限责任公司 | 并行局部序列对齐 |
| CN103902599A (zh) * | 2012-12-27 | 2014-07-02 | 北京新媒传信科技有限公司 | 模糊查找的方法和装置 |
| CN103902599B (zh) * | 2012-12-27 | 2017-04-05 | 北京新媒传信科技有限公司 | 模糊查找的方法和装置 |
| CN104395900B (zh) * | 2013-03-15 | 2017-08-25 | 北京未名博思生物智能科技开发有限公司 | 序列比对的空间计数运算方法 |
| CN104395900A (zh) * | 2013-03-15 | 2015-03-04 | 北京未名博思生物智能科技开发有限公司 | 序列比对的空间计数运算方法 |
| CN104881592A (zh) * | 2015-02-11 | 2015-09-02 | 哈尔滨工业大学深圳研究生院 | 一种dna序列比对中的打分方法 |
| CN104881592B (zh) * | 2015-02-11 | 2018-04-06 | 哈尔滨工业大学深圳研究生院 | 一种dna序列比对中的打分方法 |
| CN104899476A (zh) * | 2015-06-15 | 2015-09-09 | 中国人民解放军国防科学技术大学 | 一种对多序列bwt索引构建进行并行加速的方法 |
| CN105956417A (zh) * | 2016-05-04 | 2016-09-21 | 西安电子科技大学 | 云环境下基于编辑距离的相似碱基序列查询方法 |
| CN106529212A (zh) * | 2016-10-19 | 2017-03-22 | 哈尔滨工业大学深圳研究生院 | 基于序列依赖频率矩阵的生物序列进化信息提取方法及其应用 |
| CN106529212B (zh) * | 2016-10-19 | 2019-01-25 | 哈尔滨工业大学深圳研究生院 | 基于序列依赖频率矩阵的生物序列进化信息提取方法 |
| CN111445952A (zh) * | 2020-03-25 | 2020-07-24 | 山东大学 | 超长基因序列的相似性快速比对方法及系统 |
| CN111445952B (zh) * | 2020-03-25 | 2024-01-26 | 山东大学 | 超长基因序列的相似性快速比对方法及系统 |
| CN111916153A (zh) * | 2020-06-17 | 2020-11-10 | 电子科技大学 | 一种并行多重序列比对方法 |
| CN112614544A (zh) * | 2020-12-28 | 2021-04-06 | 杭州瑞普基因科技有限公司 | Kraken2软件输出结果的优化方法及鉴定样本中物种类型的方法 |
| CN112614544B (zh) * | 2020-12-28 | 2024-05-17 | 杭州瑞普基因科技有限公司 | Kraken2软件输出结果的优化方法及鉴定样本中物种类型的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN102750461B (zh) | 2015-04-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN102750461A (zh) | 一种可得到完全解的生物序列局部比对方法 | |
| Heather et al. | High-throughput sequencing of the T-cell receptor repertoire: pitfalls and opportunities | |
| Stranneheim et al. | Classification of DNA sequences using Bloom filters | |
| Belacel et al. | Clustering methods for microarray gene expression data | |
| CN103778349A (zh) | 一种基于功能模块的生物分子网络分析的方法 | |
| Torshizi et al. | A new cluster validity measure based on general type-2 fuzzy sets: application in gene expression data clustering | |
| Gerstein | Measurement of the effectiveness of transitive sequence comparison, through a third'intermediate'sequence. | |
| CN112735537A (zh) | 一种新的脑胶质瘤分子分型方法 | |
| CN112750502A (zh) | 一种基于二维分布结构判定的单细胞转录组测序数据聚类推荐方法 | |
| CN110910950A (zh) | 一种联合分析单细胞scRNA-seq和scATAC-seq的流程方法 | |
| CN110400602A (zh) | 一种基于测序数据的abo血型系统分型方法及其应用 | |
| CN109545283A (zh) | 一种基于序列模式挖掘算法的系统发生树构建方法 | |
| CN106156541B (zh) | 分析个体两类状态的免疫差异的方法和装置 | |
| Pirim et al. | Clustering high throughput biological data with B-MST, a minimum spanning tree based heuristic | |
| CN105426700A (zh) | 一种批量计算基因组直系同源基因进化速率的方法 | |
| KR101839572B1 (ko) | 질병 관련 유전자 관계 분석 장치 및 방법 | |
| Casey et al. | A machine learning approach to prostate cancer risk classification through use of RNA sequencing data | |
| Jothi | A betweenness centrality guided clustering algorithm and its applications to cancer diagnosis | |
| Dinu et al. | Clustering based on Rank Distance with Applications on DNA | |
| CN110534158B (zh) | 一种基因序列比对方法、装置、服务器及介质 | |
| Ryan et al. | Multi-Omic Graph Diagnosis (MOGDx): A data integration tool to perform classification tasks for heterogeneous diseases | |
| KR101482010B1 (ko) | 전체 유전체 서열분석을 위한 초고속 범용 검색장치 및 방법 | |
| KR101701168B1 (ko) | 유전자 패스웨이 활성지수의 세부적 정량화를 위한 유전자 프로파일 방법 | |
| Zhang et al. | nSEA: n-Node Subnetwork Enumeration Algorithm Identifies Lower Grade Glioma Subtypes with Altered Subnetworks and Distinct Prognostics | |
| CN107038350B (zh) | 一种药物的长非编码rna靶点预测方法和系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant |