CN110875084A - 一种核酸序列比对的方法 - Google Patents
一种核酸序列比对的方法 Download PDFInfo
- Publication number
- CN110875084A CN110875084A CN201810914779.6A CN201810914779A CN110875084A CN 110875084 A CN110875084 A CN 110875084A CN 201810914779 A CN201810914779 A CN 201810914779A CN 110875084 A CN110875084 A CN 110875084A
- Authority
- CN
- China
- Prior art keywords
- sequence
- alignment
- input
- penalty
- sequences
- 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
Landscapes
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
Abstract
本发明公开了一种核酸序列比对的方法。该方法包括如下步骤：(1)获取与输入序列类似性积分最高的参考序列；(2)从所述参考序列中截取与输入序列等长的若干个比对序列；(3)将输入序列分别和比对序列比对，获取输入序列与各个比对序列的匹配信息；(4)根据步骤(3)获取的匹配信息，获取比对序列的比对得分；比对得分越高，则输入序列和该比对序列的相似度越高。本发明提供的核酸比对方法可以进行序列的比对，尤其是可以进行小序列的比对和一些存在突变(如插入突变、缺失突变)的小序列的比对。本发明具有重大的应用价值。
Description
技术领域
本发明属于生物信息学领域，具体涉及一种核酸序列比对的方法，尤其涉及较短的核酸序列比对的方法。
背景技术
在待比对序列(如microRNA)较短的情况下，采用序列比对工具比对的准确性较低，尤其是待比对序列存在indel(即插入和/或缺失)时，基本很难比对上。例如：blast软件和bowtie2软件都可以进行小序列的比对且比对速度非常快，但对于一些存在突变(如插入突变、缺失突变)的小片段，则很难比对上。而进行小RNA分析时定量结果容易存在偏差，原因有二：一是部分物种公认的小RNA序列不一定完全准确；二是测序过程中引入的误差(如illumina仪器在高GC区域时测序不准确)导致序列存在一些小的突变，因而难以准确定位。由此可见，现有比对工具不能满足所有类型序列的精确定位，比对结果存在一定程度的不准确性，是造成靶标预测错误的原因之一。
发明内容
本发明的目的是提供一种核酸序列比对的方法，尤其是对于一些存在突变(如插入突变、缺失突变)的小片段进行比对的方法。
本发明首先保护一种核酸序列比对的方法，可包括如下步骤：
(1)获取与输入序列类似性积分最高的参考序列；
(2)从所述参考序列中截取与输入序列等长的若干个比对序列；
(3)将输入序列分别和比对序列比对，获取输入序列与各个比对序列的匹配信息；
(4)根据步骤(3)获取的匹配信息，获取比对序列的比对得分；
比对得分越高，则输入序列和该比对序列的相似度越高。
上述方法中，输入序列即待比对的核酸序列。待比对的核酸序列中可存在一些突变(如插入突变、缺失突变)。输入序列可为较短的待比对的核酸序列。所述较短的待比对的核酸序列的长度可为15-50bp(如15-30bp、30-50bp、15bp、30bp或50bp)。
所述步骤(1)可包括如下步骤：
(1-1)连续不重叠k-mer切分参考序列集中的各个参考序列，得到k-mer组合甲；将k-mer组合甲中的各个k-mer及在参考序列集中相应的参考序列信息进行存储，获得种子库；
(1-2)读取输入序列，连续重叠k-mer切分输入序列，得到N个k-mer组成的k-mer组合乙；
(1-3)将k-mer组合乙中的N个k-mer和其反向互补序列分别在所述种子库中检索，获得每个k-mer和其反向互补序列对应的参考序列，然后进行类似性积分：
如果N个k-mer对应的某参考序列类似性积分最高，则记录该参考序列及N个k-mer在该参考序列上的起始位置和终止位置；
如果N个k-mer的反向互补序列对应的某参考序列类似性积分最高，则记录该参考序列及N个k-mer的反向互补序列在该参考序列上的起始位置和终止位置；
步骤(1-3)记录的参考序列即为获取与输入序列类似性积分最高的参考序列。
所述步骤(1-1)中，k-mer组合甲中的各个k-mer的核苷酸序列可按大写字母存储；特殊嘌呤或嘧啶(即非ATCG)可按大写字母N存储。
所述步骤(1-3)中，检索时k-mer组合乙中各个k-me和其反向互补序列可按大写字母存储；特殊嘌呤或嘧啶(即非ATCG)可按大写字母N存储。
所述步骤(1-1)中，所述参考序列集可由人类参考序列集和/或动物参考序列集和/或植物参考序列集和/或微生物参考序列集组成。在本发明的一个实施例中，所述参考序列集具体可为由表1所示的30条RNA序列组成的人类参考序列集。
所述步骤(1-1)中，连续不重叠k-mer切分的长度根据输入序列的长度进行设定。例如输入序列为较短的待比对的核酸序列(15-50bp)，连续不重叠k-mer切分的长度可为4-8，优选地，k-mer切分的长度为5。
所述步骤(1-2)中，连续重叠k-mer切分的长度和步骤(1-1)中的k-mer切分长度相同。
所述步骤(1-1)中，参考序列信息可包括参考序列的名称和在参考序列上的起始位置。
所述步骤(1-1)中，参考序列信息具体由参考序列的名称和在参考序列上的起始位置组成。
所述步骤(1-3)中，类似性积分原则可为：在同一条参考序列上，如果只存在一个片段或者不连续的多个片段，计1分；在同一条参考序列上，如果存在连续的多个片段，则以连续的多个片段中片段数目最多的进行计分，每个片段计1分。例如在同一条参考序列上，分别存在N个连续的片段和M个连续的片段，N和M均为大于1的自然数且M>N，计M分。所述片段即为k-mer。
所述步骤(2)中，如果所述参考序列的长度小于输入序列，则需在5’末端和/或3’末端的相应位置补“N”。
所述步骤(3)依次可包括如下步骤：
(3-1)将输入序列和各个比对序列转化为二进制字符；
(3-2)进行位运算的异或操作；
(3-3)获取输入序列与各个比对序列的匹配信息。
所述步骤(3-1)中，二进制字符转化方式可为：A可转化为110，T可转化为101，C可转化为010，G可转化为001，N可转化为000。
所述步骤(3-2)中，进行位运算的异或操作时，相同的位为0，不同的位为1。
所述步骤(3-3)中，匹配信息中每三位代表一个碱基的匹配结果，如果三位全为0，表示该碱基匹配，否则表示该碱基不匹配。
所述步骤(4)中每个输入序列的比对得分，要分别按下述(4-1)、(4-3)和(4-4)获得。
所述步骤(4)可包括如下步骤：
(4-1)根据步骤(3)获取的匹配信息，按照罚分规则1获得各个比对序列的罚分值；然后选择罚分值为3以下的比对序列进行如下步骤，获得比对得分：
(4-1-1)根据输入序列和各个罚分值为3以下的比对序列的匹配结果，获得相应的匹配分值；每个匹配碱基计1分；
(4-1-2)获得各个罚分值为3以下的比对序列的比对得分；比对得分＝匹配分值-罚分值；
(4-2)根据步骤(3)获取的匹配信息，将输入序列分为四个区段，分别为5’末端的第1-2个碱基组成的输入区段S1、A个碱基组成的输入区段A、B个碱基组成的输入区段B和3’末端的第1-2个碱基组成的输入区段S2；A和B均为1以上的自然数；
输入区段A与输入区段B相邻的一端有C个碱基与步骤(1)获取的参考序列完全匹配；输入区段B与输入区段A相邻的一端有D个碱基与步骤(1)获取的参考序列完全不匹配；C为自然数且1＜C≤A；D为自然数且1＜D≤B；
(4-3)按照如下步骤进行插入比对；
(4-3-1)从步骤(1)获取的参考序列中截取并组装与输入序列等长的若干个比对序列；每个比对序列由第1-2个碱基组成的比对区段V1、A个碱基组成的比对区段A、一个空格、B-1个碱基组成的比对区段C和第1-2个碱基组成的比对区段V2组成；其中空格在比对区段A和比对区段C之间；A和B均为1以上的自然数；
比对区段A与空格相邻的一端的C个碱基与所述输入区段A中与输入区段B相邻的一端的C个碱基完全匹配；
(4-3-2)将输入序列分别和步骤(4-3-1)获得的比对序列比对，获取输入序列与各个比对序列的匹配信息；
(4-3-3)根据步骤(4-3-2)获取的匹配信息，按照罚分规则2获得各个比对序列的罚分值；选择罚分值为3以下的比对序列进行如下步骤，获得比对得分：
①根据输入序列和各个罚分值为3以下的比对序列的匹配结果，获得相应的匹配分值；每个匹配碱基计1分；
②获得各个罚分值为3以下的比对序列的比对得分；比对得分＝匹配分值-罚分值；
(4-4)按照如下步骤进行缺失比对；
(4-4-1)向输入序列的输入区段A和输入区段B的中间添加一个空格，作为新的输入序列；
(4-4-2)步骤(1)获取的参考序列中截取与新的输入序列等长的若干个比对序列；每个比对序列由第1-2个碱基组成的比对区段V3、A个碱基组成的比对区段A、B+1个碱基组成的比对区段D和第1-2个碱基组成的比对区段V4组成；A和B均为1以上的自然数；
比对区段A与比对区段D相邻的一端的C个碱基与所述输入区段A中与输入区段B相邻的一端的C个碱基完全匹配；
(4-4-3)将新的输入序列分别和步骤(4-4-2)获得的比对序列比对，获取新的输入序列与各个比对序列的匹配信息；
(4-4-4)根据步骤(4-4-3)获取的匹配信息，按照罚分规则2获得各个比对序列的罚分值；选择罚分值为3以下的比对序列进行如下步骤，获得比对得分：
③根据新的输入序列和各个罚分值为3以下的比对序列的匹配结果，获得相应的匹配分值；每个匹配碱基计1分；
④获得各个罚分值为3以下的比对序列的比对得分；比对得分＝匹配分值-罚分值。
步骤(4-3-1)中，如果步骤(1)获取的参考序列的长度小于所述输入序列，则需在5’末端和/或3’末端的相应位置补“N”。
步骤(4-4-2)中，如果步骤(1)获取的参考序列的长度小于新的输入序列，则需在5’末端和/或3’末端的相应位置补“N”。
步骤(4-3-2)依次可包括如下步骤：
(s1)将输入序列和步骤(4-3-1)获得的各个比对序列转化为二进制字符；
(s2)进行位运算的异或操作，获取输入序列与各个比对序列的匹配信息。
步骤(4-4-3)依次可包括如下步骤：
(x1)将新的输入序列和步骤(4-4-2)获得的各个比对序列转化为二进制字符；
(x2)进行位运算的异或操作，获取新的输入序列与各个比对序列的匹配信息。
所述步骤(s1)和所述步骤(x1)中，二进制字符转化方式可为：A可转化为110，T可转化为101，C可转化为010，G可转化为001，N可转化为000，一个空格可转化为111。
所述步骤(s2)和所述步骤(x2)中，进行位运算的异或操作时，相同的位为0，不同的位为1。
所述步骤(s2)和所述步骤(x2)中，匹配信息中每三位代表一个碱基的匹配结果，如果三位全为0，表示该碱基匹配，否则表示该碱基不匹配。
所述罚分规则1和所述罚分规则2可以相同，也可以不同。
所述罚分规则1具体可为：根据输入序列上各个碱基的匹配结果，计算罚分值。罚分值＝输入序列5’/3’末端罚分值+输入序列非5’/3’末端罚分值。输入序列的5’/3’末端由输入序列的“5’末端的第一个碱基和第二个碱基”和“3’末端的第一个碱基和第二个碱基”组成。输入序列非5’/3’末端罚分值：如果输入序列非5’/3’末端的碱基不匹配，每个不匹配的碱基罚分为1.6分；如果输入序列非5’/3’末端的碱基均匹配，罚分为0分。输入序列5’/3’末端罚分值：如果输入序列5’末端或3’末端的第一个碱基和第二个碱基均匹配，罚分为0分；如果输入序列5’末端或3’末端的第一个碱基匹配、第二个碱基不匹配，罚分为1.6分；如果输入序列5’末端或3’末端的第一个碱基不匹配、第二个碱基不匹配或匹配，罚分为1.1分。
与罚分规则1相比，罚分规则2的不同之处在于：输入序列非5’/3’末端罚分值中，如果输入序列非5’/3’末端的碱基不匹配，当碱基与空格不匹配时、罚分为1.8分，其它每个不匹配的碱基罚分为1.6分；输入序列5’/3’末端罚分值中，如果输入序列5’末端或3’末端的第一个碱基匹配、第二个碱基不匹配，罚分为1.1分。
所述方法还可包括步骤(5):按照比对得分的高低，对步骤(4)中的比对序列进行排序，输出比对结果。即分别按上述(4-1)、(4-3)和(4-4)获得每个输入序列的比对得分，然后再将比对得分按高低排队，比对得分高的比对序列输出。
所述步骤(5)中，输入比对结果默认可为前5条。
上述任一所述的方法具体是一个输入序列进行核酸序列比对的方法。
本发明还保护一种若干个输入序列进行核酸序列比对的方法，具体可为将每个输入序列按照上述任一所述的方法进行比对，然后输出比对结果。具体的，进行步骤(1-2)时，计算机可以依次读取若干个输入序列，经过比对，输入相应的比对结果。步骤(1-1)中，各个输入序列对应的参考序列集可以相同，也可以不同。
上述任一所述的方法在核酸序列比对中的应用也属于本发明的保护范围。
分别采用blast软件(网址为：https://blast.ncbi.nlm.nih.gov/Blast.cgi)、bowtie2软件(网址为：http://bowtie-bio.sourceforge.net/)和本发明提供的核酸比对的方法对10778条输入序列(部分见表2)和表1中的参考序列集进行比对。结果表明，blast软件的比对率为85.69％，比对位置一致的比对率为41.94％；bowtie2软件的比对率为45.28％，比对位置一致的比对率为41.28％；本发明提供的方法的比对率为99.88％，比对位置一致的比对率为99.88％。例如hsa_piR_017019_piRNA，23，1，7I的真实的比对结果是第7位插入了一个碱基，即1-6匹配，8-23匹配；在blast软件的比对结果中比对位置为8-23匹配(见表3)，在bowtie2软件的比对结果中是比对不上的(见表4)，本发明提供的方法的比对结果是真实的匹配情况(见表5)。由此可见，本发明提供的核酸比对方法可以进行序列的比对，尤其是可以进行小序列的比对和一些存在突变(如插入突变、缺失突变)的小序列的比对。本发明具有重大的应用价值。
具体实施方式
以下的实施例便于更好地理解本发明，但并不限定本发明。下述实施例中的实验方法，如无特殊说明，均为常规方法。下述实施例中所用的实验材料，如无特殊说明，均为自常规生化试剂商店购买得到的。以下实施例中的定量实验，均设置三次重复实验，结果取平均值。
二进制字符转化方式为A转化为110，T转化为101，C转化为010，G转化为001，N转化为000，一个空格转化为111。
罚分规则1为：根据输入序列上各个碱基的匹配结果，计算罚分值；
罚分值＝输入序列5’/3’末端罚分值+输入序列非5’/3’末端罚分值；
输入序列的5’/3’末端由输入序列的“5’末端的第一个碱基和第二个碱基”和“3’末端的第一个碱基和第二个碱基”组成；
输入序列非5’/3’末端罚分值：如果输入序列非5’/3’末端的碱基不匹配，每个不匹配的碱基罚分为1.6分；如果输入序列非5’/3’末端的碱基均匹配，罚分为0分；
输入序列5’/3’末端罚分值：如果输入序列5’末端或3’末端的第一个碱基和第二个碱基均匹配，罚分为0分；如果输入序列5’末端或3’末端的第一个碱基匹配、第二个碱基不匹配，罚分为1.6分；如果输入序列5’末端或3’末端的第一个碱基不匹配、第二个碱基不匹配或匹配，罚分为1.1分。
与罚分规则1相比，罚分规则2的不同之处在于：输入序列非5’/3’末端罚分值中，如果输入序列非5’/3’末端的碱基不匹配，当碱基与空格不匹配时、罚分为1.8分，其它每个不匹配的碱基罚分为1.6分；输入序列5’/3’末端罚分值中，如果输入序列5’末端或3’末端的第一个碱基匹配、第二个碱基不匹配，罚分为1.1分。
实施例1、建立核酸序列比对的方法
本发明的发明人经过大量实验，建立了建立核酸序列比对的方法。具体步骤如下：
1、将参考序列集中的各个参考序列均转化为大写字母(特殊嘌呤或嘧啶(即非ATCG)转化为N)，然后进行连续不重叠k-mer切分(k-mer切分长度需根据待比对序列的长度进行设定)，得到若干k-mer。
2、完成步骤1后，将各个k-mer和对应的位置的数据存储入计算机，获得种子库。
其中，对应的位置按k-mer对应的参考序列名称和在参考序列上的起始位置存储(k-mer对应的所有参考序列名称及在相应参考序列上的起始位置均需存储；如某k-mer的核苷酸序列为ATCG，其对应参考序列1(核苷酸序列为：5’-ATCGACGCG-3’)和参考序列15(核苷酸序列为：5’-ATCGTCCCG-3’)，则数据记录如下：ATCG，g1_1，g15_1)。
3、将输入序列(即与参考序列集进行比对的核酸序列)均转化为大写字母(特殊嘌呤或嘧啶(即非ATCG)转化为N)，然后进行连续重叠k-mer切分(k-mer切分长度和步骤1中的k-mer切分长度相同)，得到N(N为自然数)个k-mer。
4、用步骤3得到的N个k-mer和其反向互补序列分别在步骤2获得的种子库中检索，记录每个k-mer和其反向互补序列对应的参考序列名称，然后按照计分原则计分，如果N个k-mer对应的某参考序列计分最高，则记录该参考序列名称及N个k-mer在该参考序列上的起始位置和终止位置；如果N个k-mer的反向互补序列对应的某参考序列计分最高，则记录该参考序列名称及N个k-mer的反向互补序列在该参考序列上的起始位置和终止位置。
计分原则如下(片段即为k-mer)：
k1)在同一条参考序列上，如果只存在一个片段或者不连续的多个片段，计1分；
k2)在同一条参考序列上，如果存在连续的多个片段，则以连续的多个片段中片段数目最多的进行计分，每个片段计1分。例如在同一条参考序列上，分别存在N个连续的片段和M个连续的片段，N和M均为大于1的自然数且M>N，则计M分。
5、从步骤4中计分最高的参考序列中截取与输入序列等长的若干个比对序列。
如果步骤4中计分最高的参考序列的长度小于输入序列，则需在5’末端和/或3’末端的相应位置补“N”。
6、完成步骤5后，将输入序列和各个比对序列转化为二进制字符，然后进行位运算的异或操作(相同的位为0，不同的位为1)，得到输入序列的匹配信息，进而得到输入序列上各个碱基的匹配结果(匹配信息中每三位代表一个碱基的匹配结果，如果三位全为0，表示该碱基匹配，否则表示该碱基不匹配)。
7、完成步骤6后，根据罚分规则1获得各个比对序列的罚分值。选择罚分值为3以下的比对序列进行如下步骤，获得比对得分：
(1)根据输入序列和各个罚分值为3以下的比对序列上碱基的匹配结果，获得相应的匹配分值(每个匹配碱基计1分)；
(2)完成步骤(1)后，获得各个罚分值为3以下的比对序列的比对得分(比对得分＝匹配分值-罚分值)。
8、根据步骤6得到的输入序列的匹配信息，将输入序列分为四个区段，分别为5’末端的第1-2个碱基组成的输入区段S1、A个(A为1以上的自然数)碱基组成的输入区段A、B个(B为1以上的自然数)碱基组成的输入区段B和3’末端的第1-2个碱基组成的输入区段S2。输入区段A与输入区段B相邻的一端有C个碱基与步骤4中计分最高的参考序列完全匹配。输入区段B与输入区段A相邻的一端有D个碱基与步骤4中计分最高的参考序列完全不匹配。C为自然数且1＜C≤A；D为自然数且1＜D≤B。
如果步骤4中N个k-mer对应的某参考序列计分最高，则输入序列各个区段从上游至下游依次为：输入区段S1、输入区段A、输入区段B、输入区段S2。
如果步骤4中N个k-mer的反向互补序列对应的某参考序列计分最高，则输入序列各个区段从上游至下游依次为：输入区段S2、输入区段A、输入区段B、输入区段S1。
将输入序列分别进行插入比对和缺失比对：
(1)插入比对
(a1)从步骤4中计分最高的参考序列中截取并组装与输入序列等长的若干个比对序列；如果步骤4中计分最高的参考序列的长度小于输入序列，则需在5’末端和/或3’末端的相应位置补“N”；
每个比对序列由第1-2个碱基组成的比对区段V1、A个(A为1以上的自然数)碱基组成的比对区段A、一个空格、B-1个(B为1以上的自然数)碱基组成的比对区段C和第1-2个碱基组成的比对区段V2组成；其中空格在比对区段A和比对区段C之间。比对区段A与空格相邻的一端的C个碱基与输入区段A中与输入区段B相邻的一端的C个碱基完全匹配。
(a2)将输入序列和步骤(a1)获得的各个比对序列转化为二进制字符，然后进行位运算的异或操作(相同的位为0，不同的位为1)，得到输入序列的匹配信息，进而得到输入序列上各个碱基的匹配结果(匹配信息中每三位代表一个碱基的匹配结果，如果三位全为0，表示该碱基匹配，否则表示该碱基不匹配)。
(a3)完成步骤(a2)后，根据罚分规则2获得各个比对序列的罚分值；选择罚分值为3以下的比对序列再进行如下步骤，获得比对得分：
(k1)根据输入序列和各个罚分值为3以下的比对序列上碱基的匹配结果，获得相应的匹配分值(每个匹配碱基计1分)；
(k2)完成步骤(k1)后，获得各个罚分值为3以下的比对序列的比对得分(比对得分＝匹配分值-罚分值)。
(2)缺失比对
(b1)向输入序列的输入区段A和输入区段B的中间添加一个空格，作为新的输入序列；
(b2)从步骤4中计分最高的参考序列中截取与新的输入序列等长的若干个比对序列。如果步骤4中计分最高的参考序列的长度小于新的输入序列，则需在5’末端和/或3’末端的相应位置补“N”。
每个比对序列由第1-2个碱基组成的比对区段V3、A个(A为1以上的自然数)碱基组成的比对区段A、B+1个(B为1以上的自然数)碱基组成的比对区段D和第1-2个碱基组成的比对区段V4组成。比对区段A与比对区段D相邻的一端的C个碱基与所述输入区段A中与输入区段B相邻的一端的C个碱基完全匹配。
(b3)将新的输入序列和步骤(b2)获得的各个比对序列转化为二进制字符，然后进行位运算的异或操作(相同的位为0，不同的位为1)，得到新的输入序列的匹配信息，进而得到新的输入序列上各个碱基的匹配结果(匹配信息中每三位代表一个碱基的匹配结果，如果三位全为0，表示该碱基匹配，否则表示该碱基不匹配)。
(b4)完成步骤(b3)后，根据罚分规则2获得各个比对序列的罚分值；选择罚分值为3以下的比对序列再进行如下步骤，获得比对得分：
(h1)根据新的输入序列和各个罚分值为3以下的比对序列上碱基的匹配结果，获得相应的匹配分值(每个匹配碱基计1分)；
(h2)完成步骤(h1)后，获得各个罚分值为3以下的比对序列的比对得分(比对得分＝匹配分值-罚分值)。
9、将步骤7和步骤8中比对得分最高的比对序列(默认为前5条)作为比对结果输出。比对结果中，比对得分越高，则输入序列和该比对序列的相似度越高。
实施例2、实施例1建立的核酸比对的方法与现有比对软件的比较
一、参考序列集的获得
本发明的发明人从人的小RNA数据库中随机选取了表1所示的30条RNA序列作为参考序列。30个参考序列组成参考序列集。
表1.参考序列
二、输入序列的获得
从表1的30条参考序列中产生了10778条输入序列，包括随机位置的错配、首末端的滑动、随机位置的插入、随机位置的缺失等等。hsa_piR_017019_piRNA的部分输入序列见表2。
表2.输入序列
| 输入序列名称 | 核苷酸序列(5’-3’) |
| hsa_piR_017019_piRNA,19,1,1L | cCTAGGGAGCTGAGATGGC |
| hsa_piR_017019_piRNA,20,1,2L | caCTAGGGAGCTGAGATGGC |
| hsa_piR_017019_piRNA,19,1,1R | CTAGGGAGCTGAGATGGCc |
| hsa_piR_017019_piRNA,20,1,2R | CTAGGGAGCTGAGATGGCtc |
| hsa_piR_017019_piRNA,19,1,18I | CTAGGGAGCTGAGATGGgC |
| hsa_piR_017019_piRNA,17,1,18D | CTAGGGAGCTGAGATGG |
| hsa_piR_017019_piRNA,18,1,18M | CTAGGGAGCTGAGATGGg |
| hsa_piR_017019_piRNA,19,1,3I | CTtAGGGAGCTGAGATGGC |
| hsa_piR_017019_piRNA,17,1,3D | CTGGGAGCTGAGATGGC |
| hsa_piR_017019_piRNA,18,1,3M | CTtGGGAGCTGAGATGGC |
| hsa_piR_017019_piRNA,19,1,2I | CaTAGGGAGCTGAGATGGC |
| hsa_piR_017019_piRNA,23,1,7I | CTAGGGtAGCTGAGATGGCTTAG |
注：输入序列名称的格式为参考序列名称，比对长度，比对不上的碱基数，比对不上的类型。小写字母代表和参考序列比对不上：在首末端第1位和第2位需要滑动才能匹配，其它位置为错配。L为左端滑动，R为右端滑动，I为插入，D为缺失，M为错配。
三、采用blast软件比对
采用blast软件(网址为：https://blast.ncbi.nlm.nih.gov/Blast.cgi)分别对步骤二的10778条输入序列和步骤一的参考序列集进行比对。比对参数为：-p blastn-e1e-5-m 8。
部分比对结果见表3。结果表明，比对率为85.69％，比对位置一致的比对率为41.94％。
表3.blast的比对结果
四、采用bowtie2软件比对
采用bowtie2软件(网址为：http://bowtie-bio.sourceforge.net/)分别对步骤二的10778条输入序列和步骤一的参考序列集进行比对。比对参数为：-f-L 13--sensitive–quiet。
部分比对结果见表4。结果表明，比对率为45.28％，比对位置一致的比对率为41.28％。
表4.bowtie2的比对结果
五、采用实施例1建立的核酸比对的方法比对
采用实施例1建立的核酸比对的方法，分别对步骤二的10778条输入序列和步骤一的参考序列集进行比对，k-mer长度设置为5。
比对结果见表5。结果表明，比对率为99.88％，比对位置一致的比对率为99.88％；其中未比对上的13条，属于参考序列中存在部分序列是完全相似的片段，导致定位到相似片段而比对不上的。例如hsa_piR_017019_piRNA，23，1，7I的真实的比对结果是第7位插入了一个碱基，即1-6匹配，8-23匹配；在表3的比对结果中比对位置为8-23匹配，在表4的比对结果中是比对不上的，在表5的比对结果是真实的匹配情况。
由此可见，采用实施例1建立的核酸比对的方法可以进行小序列的比对，尤其是可以进行一些存在突变(如插入突变、缺失突变)的小序列的比对。
表5.实施例1建立的方法的比对结果
| 输入序列名称 | 参考序列名称 | 正负链 | 起始位置 | 得分值 | 比对细节 |
| hsa_piR_017019_piRNA,19,1,1L | hsa_piR_017019_piRNA | + | 第1位 | 16.9 | 1S18M |
| hsa_piR_017019_piRNA,20,1,2L | hsa_piR_017019_piRNA | + | 第1位 | 15.8 | 2S18M |
| hsa_piR_017019_piRNA,19,1,1R | hsa_piR_017019_piRNA | + | 第1位 | 16.9 | 18M1S |
| hsa_piR_017019_piRNA,23,1,7I | hsa_piR_017019_piRNA | + | 第2位 | 19.2 | 6M1I16M |
| hsa_piR_017019_piRNA,19,1,18I | hsa_piR_017019_piRNA | + | 第1位 | 14.8 | 17M2S |
| hsa_piR_017019_piRNA,17,1,18D | hsa_piR_017019_piRNA | + | 第1位 | 17 | 17M |
| hsa_piR_017019_piRNA,18,1,18M | hsa_piR_017019_piRNA | + | 第1位 | 15.9 | 17M1S |
| hsa_piR_017019_piRNA,19,1,3I | hsa_piR_017019_piRNA | + | 第1位 | 14.8 | 2S17M |
| hsa_piR_017019_piRNA,17,1,3D | hsa_piR_017019_piRNA | + | 第2位 | 12.8 | 2S15M |
| hsa_piR_017019_piRNA,18,1,3M | hsa_piR_017019_piRNA | + | 第1位 | 15.4 | 2M1X15M |
注：比对细节的格式为从5’至3’N个碱基，比对不上的类型。S表示和参考序列比对不上，需要滑动才能匹配，I为插入，D为缺失，M为匹配，X为错配。
Claims (10)
1.一种核酸序列比对的方法，包括如下步骤：
(1)获取与输入序列类似性积分最高的参考序列；
(2)从所述参考序列中截取与输入序列等长的若干个比对序列；
(3)将输入序列分别和比对序列比对，获取输入序列与各个比对序列的匹配信息；
(4)根据步骤(3)获取的匹配信息，获取比对序列的比对得分；
比对得分越高，则输入序列和该比对序列的相似度越高。
2.如权利要求1所述的方法，其特征在于：
所述步骤(1)包括如下步骤：
(1-1)连续不重叠k-mer切分参考序列集中的各个参考序列，得到k-mer组合甲；将k-mer组合甲中的各个k-mer及在参考序列集中相应的参考序列信息进行存储，获得种子库；
(1-2)读取输入序列，连续重叠k-mer切分输入序列，得到N个k-mer组成的k-mer组合乙；
(1-3)将k-mer组合乙中的N个k-mer和其反向互补序列分别在所述种子库中检索，获得每个k-mer和其反向互补序列对应的参考序列，然后进行类似性积分：
如果N个k-mer对应的某参考序列类似性积分最高，则记录该参考序列及N个k-mer在该参考序列上的起始位置和终止位置；
如果N个k-mer的反向互补序列对应的某参考序列类似性积分最高，则记录该参考序列及N个k-mer的反向互补序列在该参考序列上的起始位置和终止位置；
步骤(1-3)记录的参考序列即为获取与输入序列类似性积分最高的参考序列。
3.如权利要求2所述的方法，其特征在于：
所述步骤(1-1)中，参考序列信息包括参考序列的名称和在参考序列上的起始位置。
4.如权利要求2所述的方法，其特征在于：
所述步骤(1-3)中，类似性积分原则为：在同一条参考序列上，如果只存在一个片段或者不连续的多个片段，计1分；在同一条参考序列上，如果存在连续的多个片段，则以连续的多个片段中片段数目最多的进行计分，每个片段计1分。
5.如权利要求1所述的方法，其特征在于：所述步骤(3)依次包括如下步骤：
(3-1)将输入序列和各个比对序列转化为二进制字符；
(3-2)进行位运算的异或操作；
(3-3)获取输入序列与各个比对序列的匹配信息。
6.如权利要求1所述的方法，其特征在于：所述步骤(4)包括如下步骤：
(4-1)根据步骤(3)获取的匹配信息，按照罚分规则1获得各个比对序列的罚分值；然后选择罚分值为3以下的比对序列进行如下步骤，获得比对得分：
(4-1-1)根据输入序列和各个罚分值为3以下的比对序列的匹配结果，获得相应的匹配分值；每个匹配碱基计1分；
(4-1-2)获得各个罚分值为3以下的比对序列的比对得分；比对得分＝匹配分值-罚分值；
(4-2)根据步骤(3)获取的匹配信息，将输入序列分为四个区段，分别为5’末端的第1-2个碱基组成的输入区段S1、A个碱基组成的输入区段A、B个碱基组成的输入区段B和3’末端的第1-2个碱基组成的输入区段S2；A和B均为1以上的自然数；
输入区段A与输入区段B相邻的一端有C个碱基与步骤(1)获取的参考序列完全匹配；输入区段B与输入区段A相邻的一端有D个碱基与步骤(1)获取的参考序列完全不匹配；C为自然数且1＜C≤A；D为自然数且1＜D≤B；
(4-3)按照如下步骤进行插入比对；
(4-3-1)从步骤(1)获取的参考序列中截取并组装与输入序列等长的若干个比对序列；每个比对序列由第1-2个碱基组成的比对区段V1、A个碱基组成的比对区段A、一个空格、B-1个碱基组成的比对区段C和第1-2个碱基组成的比对区段V2组成；其中空格在比对区段A和比对区段C之间；A和B均为1以上的自然数；
比对区段A与空格相邻的一端的C个碱基与所述输入区段A中与输入区段B相邻的一端的C个碱基完全匹配；
(4-3-2)将输入序列分别和步骤(4-3-1)获得的比对序列比对，获取输入序列与各个比对序列的匹配信息；
(4-3-3)根据步骤(4-3-2)获取的匹配信息，按照罚分规则2获得各个比对序列的罚分值；选择罚分值为3以下的比对序列进行如下步骤，获得比对得分：
①根据输入序列和各个罚分值为3以下的比对序列的匹配结果，获得相应的匹配分值；每个匹配碱基计1分；
②获得各个罚分值为3以下的比对序列的比对得分；比对得分＝匹配分值-罚分值；
(4-4)按照如下步骤进行缺失比对；
(4-4-1)向输入序列的输入区段A和输入区段B的中间添加一个空格，作为新的输入序列；
(4-4-2)步骤(1)获取的参考序列中截取与新的输入序列等长的若干个比对序列；每个比对序列由第1-2个碱基组成的比对区段V3、A个碱基组成的比对区段A、B+1个碱基组成的比对区段D和第1-2个碱基组成的比对区段V4组成；A和B均为1以上的自然数；
比对区段A与比对区段D相邻的一端的C个碱基与所述输入区段A中与输入区段B相邻的一端的C个碱基完全匹配；
(4-4-3)将新的输入序列分别和步骤(4-4-2)获得的比对序列比对，获取新的输入序列与各个比对序列的匹配信息；
(4-4-4)根据步骤(4-4-3)获取的匹配信息，按照罚分规则2获得各个比对序列的罚分值；选择罚分值为3以下的比对序列进行如下步骤，获得比对得分：
③根据新的输入序列和各个罚分值为3以下的比对序列的匹配结果，获得相应的匹配分值；每个匹配碱基计1分；
④获得各个罚分值为3以下的比对序列的比对得分；比对得分＝匹配分值-罚分值。
7.如权利要求6所述的方法，其特征在于：
所述罚分规则1为：根据输入序列上各个碱基的匹配结果，计算罚分值；罚分值＝输入序列5’/3’末端罚分值+输入序列非5’/3’末端罚分值；输入序列的5’/3’末端由输入序列的“5’末端的第一个碱基和第二个碱基”和“3’末端的第一个碱基和第二个碱基”组成；输入序列非5’/3’末端罚分值：如果输入序列非5’/3’末端的碱基不匹配，每个不匹配的碱基罚分为1.6分；如果输入序列非5’/3’末端的碱基均匹配，罚分为0分；输入序列5’/3’末端罚分值：如果输入序列5’末端或3’末端的第一个碱基和第二个碱基均匹配，罚分为0分；如果输入序列5’末端或3’末端的第一个碱基匹配、第二个碱基不匹配，罚分为1.6分；如果输入序列5’末端或3’末端的第一个碱基不匹配、第二个碱基不匹配或匹配，罚分为1.1分；
与罚分规则1相比，罚分规则2的不同之处在于：输入序列非5’/3’末端罚分值中，如果输入序列非5’/3’末端的碱基不匹配，当碱基与空格不匹配时、罚分为1.8分，其它每个不匹配的碱基罚分为1.6分；输入序列5’/3’末端罚分值中，如果输入序列5’末端或3’末端的第一个碱基匹配、第二个碱基不匹配，罚分为1.1分。
8.如权利要求1至7任一所述的方法，其特征在于：所述方法还包括步骤(5):按照比对得分的高低，对步骤(4)中的比对序列进行排序，输出比对结果。
9.一种若干个输入序列进行核酸序列比对的方法，为将每个输入序列按照权利要求1至8任一所述的方法进行比对，然后输出比对结果。
10.权利要求1至9任一所述的方法在核酸序列比对中的应用。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810914779.6A CN110875084B (zh) | 2018-08-13 | 2018-08-13 | 一种核酸序列比对的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810914779.6A CN110875084B (zh) | 2018-08-13 | 2018-08-13 | 一种核酸序列比对的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110875084A true CN110875084A (zh) | 2020-03-10 |
| CN110875084B CN110875084B (zh) | 2022-06-21 |
Family
ID=69714103
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810914779.6A Active CN110875084B (zh) | 2018-08-13 | 2018-08-13 | 一种核酸序列比对的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110875084B (zh) |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101430741A (zh) * | 2008-12-12 | 2009-05-13 | 深圳华大基因研究院 | 一种短序列映射方法及系统 |
| CN102682226A (zh) * | 2012-04-18 | 2012-09-19 | 盛司潼 | 一种核酸测序信息处理系统及方法 |
| EP2677449A1 (en) * | 2012-06-19 | 2013-12-25 | Sjöblom, Tobias | Method and device for efficient calculation of allele ratio confidence intervals and uses thereof |
| CN105069325A (zh) * | 2012-07-28 | 2015-11-18 | 盛司潼 | 一种对核酸序列信息进行匹配的方法 |
| CN105243297A (zh) * | 2015-10-09 | 2016-01-13 | 人和未来生物科技（长沙）有限公司 | 一种参考基因组上基因序列片段的快速比对定位方法 |
| CN106202991A (zh) * | 2016-06-30 | 2016-12-07 | 厦门艾德生物医药科技股份有限公司 | 一种基因组多重扩增测序产物中突变信息的检测方法 |
| CN106676182A (zh) * | 2017-02-07 | 2017-05-17 | 北京诺禾致源科技股份有限公司 | 一种低频率基因融合的检测方法及装置 |
| CN106874709A (zh) * | 2015-12-12 | 2017-06-20 | 北京大学 | 测序结果中序列数据错误的检测和校正方法 |
| CN106951731A (zh) * | 2017-03-28 | 2017-07-14 | 上海至本生物科技有限公司 | 一种大片段插入或缺失的预测方法及系统 |
| CN107208152A (zh) * | 2015-03-06 | 2017-09-26 | 深圳华大基因股份有限公司 | 检测突变簇的方法和装置 |
| CN107488660A (zh) * | 2017-09-07 | 2017-12-19 | 中国农业科学院兰州兽医研究所 | 回文与互补回文小rna及其应用 |
-
2018
- 2018-08-13 CN CN201810914779.6A patent/CN110875084B/zh active Active
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101430741A (zh) * | 2008-12-12 | 2009-05-13 | 深圳华大基因研究院 | 一种短序列映射方法及系统 |
| CN102682226A (zh) * | 2012-04-18 | 2012-09-19 | 盛司潼 | 一种核酸测序信息处理系统及方法 |
| EP2677449A1 (en) * | 2012-06-19 | 2013-12-25 | Sjöblom, Tobias | Method and device for efficient calculation of allele ratio confidence intervals and uses thereof |
| CN105069325A (zh) * | 2012-07-28 | 2015-11-18 | 盛司潼 | 一种对核酸序列信息进行匹配的方法 |
| CN107208152A (zh) * | 2015-03-06 | 2017-09-26 | 深圳华大基因股份有限公司 | 检测突变簇的方法和装置 |
| CN105243297A (zh) * | 2015-10-09 | 2016-01-13 | 人和未来生物科技（长沙）有限公司 | 一种参考基因组上基因序列片段的快速比对定位方法 |
| CN106874709A (zh) * | 2015-12-12 | 2017-06-20 | 北京大学 | 测序结果中序列数据错误的检测和校正方法 |
| CN106202991A (zh) * | 2016-06-30 | 2016-12-07 | 厦门艾德生物医药科技股份有限公司 | 一种基因组多重扩增测序产物中突变信息的检测方法 |
| CN106676182A (zh) * | 2017-02-07 | 2017-05-17 | 北京诺禾致源科技股份有限公司 | 一种低频率基因融合的检测方法及装置 |
| CN106951731A (zh) * | 2017-03-28 | 2017-07-14 | 上海至本生物科技有限公司 | 一种大片段插入或缺失的预测方法及系统 |
| CN107488660A (zh) * | 2017-09-07 | 2017-12-19 | 中国农业科学院兰州兽医研究所 | 回文与互补回文小rna及其应用 |
Non-Patent Citations (5)
| Title |
|---|
| APANGHUANG153: "DNA序列比较与比对", 《道客巴巴：HTTPS://WWW.DOC88.COM/P-3708091261212.HTML?R=1》 * |
| BARRY G. HALL: "Comparison of the Accuracies of Several Phylogenetic Methods Using Protein and DNA Sequences", 《MOLECULAR BIOLOGY AND EVOLUTION》 * |
| SHUYINGL2266: "NCBI序列比对方法与实例操作", 《道客巴巴：HTTP://WWW.DOC88.COM/P-9703197609868.HTML》 * |
| YANG Z等: "The improvement and implementation on the algorithm for local alignment of pairs of DNA sequences", 《2016 IEEE ADVANCED INFORMATION MANAGEMENT, COMMUNICATES, ELECTRONIC AND AUTOMATION CONTROL CONFERENCE (IMCEC)》 * |
| 孟凡荣等: "在基因序列中高效引入多位点突变的方法", 《中国肺癌杂志》 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110875084B (zh) | 2022-06-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Song et al. | New developments of alignment-free sequence comparison: measures, statistics and next-generation sequencing | |
| Nilsson et al. | Five simple guidelines for establishing basic authenticity and reliability of newly generated fungal ITS sequences | |
| US10453559B2 (en) | Method and system for rapid searching of genomic data and uses thereof | |
| CN107403075B (zh) | 比对方法、装置及系统 | |
| CN107798216B (zh) | 采用分治法进行高相似性序列的比对方法 | |
| CN105243297A (zh) | 一种参考基因组上基因序列片段的快速比对定位方法 | |
| WO2018218788A1 (zh) | 一种基于全局种子打分优选的三代测序序列比对方法 | |
| Li et al. | Analysis of transposable elements in the genome of Asparagus officinalis from high coverage sequence data | |
| CN115631789A (zh) | 一种基于泛基因组的群体联合变异检测方法 | |
| He et al. | De novo assembly methods for next generation sequencing data | |
| CN115240770A (zh) | 一种检测短串联重复扩张和基因分型的方法、电子设备及存储介质 | |
| CN102841988B (zh) | 一种对核酸序列信息进行匹配的系统和方法 | |
| CN110875084B (zh) | 一种核酸序列比对的方法 | |
| CN108388772A (zh) | 一种利用文本比对分析高通量测序基因表达水平的方法 | |
| CN114424288A (zh) | 用于确定与两个突变的序列读段衍生自包含突变的相同序列的概率相关的量度的方法 | |
| CN114999572B (zh) | 一种设计引物的方法、设备、可读介质及装置 | |
| CN108733974A (zh) | 一种基于高通量测序的线粒体序列拼接及拷贝数测定的方法 | |
| CN114564306B (zh) | 一种基于GPU并行计算的第三代测序RNA-seq比对方法 | |
| CN116130001A (zh) | 一种基于k-mer定位的三代序列比对算法 | |
| CN116050348A (zh) | 一种fastq文件的拆分方法、系统、电子设备及存储介质 | |
| CN107688727B (zh) | 生物序列聚类和全长转录组中转录本亚型识别方法和装置 | |
| KR101482010B1 (ko) | 전체 유전체 서열분석을 위한 초고속 범용 검색장치 및 방법 | |
| Fu et al. | A parsimony approach to genome-wide ortholog assignment | |
| WO2019023978A1 (zh) | 比对方法、装置及系统 | |
| KR20130101711A (ko) | 시드의 길이를 고려한 염기 서열 처리 시스템 및 방법 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |