CN112825268B - 测序结果比对方法及其应用 - Google Patents
测序结果比对方法及其应用 Download PDFInfo
- Publication number
- CN112825268B CN112825268B CN201911148507.0A CN201911148507A CN112825268B CN 112825268 B CN112825268 B CN 112825268B CN 201911148507 A CN201911148507 A CN 201911148507A CN 112825268 B CN112825268 B CN 112825268B
- Authority
- CN
- China
- Prior art keywords
- nucleic acid
- acid sequence
- basic unit
- sequence
- seed
- 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.)
- Active
Links
- 238000012163 sequencing technique Methods 0.000 title claims abstract description 63
- 238000000034 method Methods 0.000 title claims abstract description 42
- 230000011218 segmentation Effects 0.000 claims abstract description 6
- 150000007523 nucleic acids Chemical group 0.000 claims description 97
- 108091028043 Nucleic acid sequence Proteins 0.000 claims description 92
- 239000012634 fragment Substances 0.000 claims description 41
- 239000011159 matrix material Substances 0.000 claims description 32
- 238000011144 upstream manufacturing Methods 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 3
- 238000004590 computer program Methods 0.000 claims description 2
- 108091032955 Bacterial small RNA Proteins 0.000 description 7
- 230000035772 mutation Effects 0.000 description 5
- 239000012297 crystallization seed Substances 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 238000002864 sequence alignment Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 108020004707 nucleic acids Proteins 0.000 description 3
- 102000039446 nucleic acids Human genes 0.000 description 3
- 108091007412 Piwi-interacting RNA Proteins 0.000 description 2
- 238000012167 Small RNA sequencing Methods 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 238000012217 deletion Methods 0.000 description 2
- 230000037430 deletion Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000003780 insertion Methods 0.000 description 2
- 230000037431 insertion Effects 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 239000004055 small Interfering RNA Substances 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- -1 Bowtie2 Proteins 0.000 description 1
- 108091030146 MiRBase Proteins 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000036438 mutation frequency Effects 0.000 description 1
- 238000001921 nucleic acid quantification Methods 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 230000009467 reduction Effects 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
- G16B30/10—Sequence alignment; Homology search
-
- 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
- G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
- G16B20/30—Detection of binding sites or motifs
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Engineering & Computer Science (AREA)
- Biophysics (AREA)
- General Health & Medical Sciences (AREA)
- Analytical Chemistry (AREA)
- Chemical & Material Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Medical Informatics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Genetics & Genomics (AREA)
- Molecular Biology (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明提出了一种测序结果比对方法。该方法包括：将参考序列按照第一预定长度进行k‑mer切分，构建索引库；基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所述种子序列集合中所述种子序列的数目大于比对所允许的错配数，并且所述种子序列的长度不超过所述第一预定长度；将所述种子序列集合的至少一部分与所述索引库进行匹配；将测序读段与匹配上的参考序列进行全局比对。
Description
技术领域
本发明涉及核酸信息领域，具体地，本发明涉及测序结果比对方法、计算机可读存储介质、电子设备以及测序结果比对系统。
背景技术
目前主流的短序列比对工具对于片段比较短(譬如18个碱基)的情况，比对的准确性较低，尤其是这类短片段中存在小的插入、缺失时，基本很难比对上。而对于小RNA分析，由于部分物种参考序列不一定完全准确，或者测序过程中引入的误差(譬如illumina机器在高GC区时测序不准确)等原因，导致测序得到的片段存在一些小的突变，却难以准确定位到相应的参考序列，使得最终的定量结果存在偏差。
目前主要应用于小RNA的比对工具，如Bowtie2、BWA等，都是基于种子序列延伸的策略进行比对；对于片段较短且可能存在的突变，导致这些比对工具在进行种子搜索时，容易检索不到完全匹配的片段，从而导致序列不能被比对上，进而影响了下游的小RNA定量。
种子序列定位及延伸算法，通过扫描参考基因组序列，对参考基因组序列建立哈希表，将序列分成一定长度的小片段，这种小片段也被称之为种子。然后，在目标序列中查找和种子序列相同的片段并标记，以这些标记点为锚点向左右按一定规律延伸比对，将不合条件的舍弃，符合条件的结果将输出保存。如果种子中存在错配，整条read的比对就不会被进行，所以就会被认定为非比对上，导致丢失部分比对信息。
因此，对于小片段核酸序列的比对方法仍需要科研工作者的进一步开发和改进。
发明内容
本发明旨在至少在一定程度上解决相关技术中的技术问题之一。
在本发明的第一方面，本发明提出了一种测序结果比对方法。根据本发明的实施例，所述方法包括：将参考序列按照第一预定长度进行k-mer(序列中包含k个碱基的字符串)切分，构建索引库；基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所述种子序列集合中所述种子序列的数目大于比对所允许的错配数，并且所述种子序列的长度不超过所述第一预定长度；将所述种子序列集合的至少一部分与所述索引库进行匹配，将测序读段与匹配上的参考序列进行全局比对。根据本发明实施例的上述测序结果比对方法通过采取比测序的平均错配数多的多个种子序列检索模式，保证了每条测序读段至少有一条种子序列的匹配，避免了因为种子序列错配而导致匹配不上的情况；同时，通过种子序列与索引库先匹配，将参考读段与匹配上的参考序列进行全局比对，而非每条测序读段与所有参考序列进行全局比对，减少了比对的运算量，极大提高了比对速度；并且通过全局比对，极大地满足了各种类型的小核酸序列比对，如目的片段比参考序列长、目的片段中存在突变等；提高了小核酸数据的利用率，同时提升了小核酸定量的准确性。
根据本发明的实施例，上述方法还可以进一步包括如下附加技术特征至少之一：
根据本发明的实施例，所述全局比对是通过如下方式进行的：获取待比对的来源于所述测序读段的第一核酸序列与来源于所述匹配上的参考序列的第二核酸序列各位置上的基本单元信息；基于所述基本单元信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列的基本单元数目，n为所述第二核酸序列的基本单元数目，其中所述得分矩阵中的元素Mij表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j个基本单元的比对得分；基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得经过所述第一核酸序列与所述第二核酸序列的比对结果，其中，所述元素Mij是基于下列公式确定的：
其中，Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；Mi-1,j表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j个基本单元的比对得分；g表示小于零的第一预定数值；S(Ri,Sj)是基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相同时，S(Ri,Sj)为第三预定数值，所述第三预定数值小于所述第二预定数值。
需要说明的是，本申请所述的“比对所允许的错配数”是指在具体比对时，所允许的容错减基数。例如，在小RNA测序中，小RNA比较短，比对所允许的错配数一般为1，如果条件更严格一些，就是不允许错配；如果所允许的错配数大于1，比对结果可能不太可靠，导致后续定量等分析结果相对不那么准确。所以在小RNA测序中，比对所允许的错配数一般是1个错配。
根据本发明的实施例，所述第一预定数值为不小于-10的整数，优选-5。
根据本发明的实施例，所述第二预定数值为1。
根据本发明的实施例，所述第三预定数值为-2。
模拟数据测试时发明人发现，所述第一预定数值为-5、所述第二预定数值为1、所述第三预定数值为-2的预定数值组合下，数据的比对效果最佳。
根据本发明的实施例，所述基本单元为碱基。
根据本发明的实施例，所述第一预定长度为不超过20的整数，优选不超过15，更优选不超过10，最优选9。发明人发现，种子序列的长度过长，可能会因为碱基错配导致整条比对不上，从而丢失部分比对结果；种子长度过短会增加运算时间，比对速度较慢，但结果更准确。本方法通过比较不同长度后优选设置长度为9，平衡了比对速度和结果准确性。
根据本发明的实施例，所述种子序列集合是通过下列步骤确定的：(1)将所述测序读段的5’末端和3’末端各去除第二预定长度的碱基；(2)将经过步骤(1)处理的所述测序读段划分为多个长度相同的子片段，所述子片段的数目超过所述测序的平均错配数；(3)在每个子片段中，由5’末端起始基于所述第一预定长度，确定所述种子序列，其中，如果所述子片段的长度小于所述第一预定长度，则将所述子片段作为所述种子序列，如果所述子片段的长度大于所述第一预定长度，则由所述子片段作的5’末端起始延伸所述第一预定长度作为所述种子序列。
根据本发明的实施例，通过数据测试，发明人发现种子序列的长度为9比较合适。发明人发现，种子序列的长度决定了比对的速度和准确性，长度越短比对越准但更慢，长度越长则快但可能会丢失部分比对信息；
根据本发明的实施例，所述第二预定长度为不超过2的整数。
根据本发明的实施例，所述子片段的数目＝所述比对所允许的错配数+1。进而保证了一条测序读段至少有一个种子序列能够比对到参考序列。
在本发明的第二方面，本发明提出了一种计算机可读存储介质，其上存储有计算机程序。根据本发明的实施例，该程序被处理器执行时实现前面所述的测序结果比对方法。
在本发明的第三方面，本发明提出了一种电子设备。根据本发明的实施例，所述电子设备包括存储器、处理器；其中，所述处理器通过读取所述存储器中存储的可执行程序代码来运行与所述可执行程序代码对应的程序，以用于实现前面所述的测序结果比对方法。
在本发明的第四方面，本发明提出了一种测序结果比对系统。根据本发明的实施例，所述系统包括：构建索引库装置，所述构建索引库装置用于将参考序列按照第一预定长度进行k-mer切分，构建索引库；种子序列集合确定装置，所述种子序列集合确定装置用于基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所述种子序列集合中所述种子序列的数目大于比对所允许的错配数，并且所述种子序列的长度不超过所述第一预定长度；匹配装置，所述匹配装置用于将所述种子序列集合的至少一部分与所述索引库进行匹配；比对装置，所述比对装置用于将测序读段与匹配上的参考序列进行全局比对。
根据本发明的实施例，所述比对装置包括：
获取基本单元信息单元，所述获取基本单元信息单元用于获取待比对的来源于所述测序读段的第一核酸序列与来源于所述匹配上的参考序列的第二核酸序列各位置上的基本单元信息；
构建得分矩阵单元，所述构建得分矩阵单元与所述获取基本单元信息单元相连，用于基于所述基本单元信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列的基本单元数目，n为所述第二核酸序列的基本单元数目，其中所述得分矩阵中的元素Mij表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j个基本单元的比对得分；
回溯单元，所述回溯单元与所述构建得分矩阵单元相连，用于基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得经过所述第一核酸序列与所述第二核酸序列的比对结果，其中，
所述元素Mij是基于下列公式确定的：
其中，
Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；
Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；
Mi-1,j表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j个基本单元的比对得分；
g表示小于零的第一预定数值；
S(Ri,Sj)是基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相同时，S(Ri,Sj)为第三预定数值，所述第三预定数值小于所述第二预定数值。
根据本发明的实施例，所述种子序列集合确定装置进一步包括：末端去除单元，所述末端去除单元用于将所述测序读段的5’末端和3’末端各去除第二预定长度的碱基；划分子片段单元，所述划分子片段单元用于将经过末端去除单元处理的所述测序读段划分为多个长度相同的子片段，所述子片段的数目超过所述测序的平均错配数；种子序列确定单元，所述种子序列确定单元用于在每个子片段中，由5’末端起始基于所述第一预定长度，确定所述种子序列，其中，如果所述子片段的长度小于所述第一预定长度，则将所述子片段作为所述种子序列，如果所述子片段的长度大于所述第一预定长度，则由所述子片段作的5’末端起始延伸所述第一预定长度作为所述种子序列。
根据本发明的实施例，所述回溯单元进一步包括；确定回溯起始位置模块，所述确定回溯起始位置模块用于确定所述矩阵Mmn中的最大值所对应的回溯起始位置；确定下一回溯位置模块，所述确定下一回溯位置模块用于基于所述回溯起始位置上游相邻三个位置的数值，确定下一回溯位置，其中，所述上游相邻三个位置包括行相邻位置、对角线相邻位置和列相邻位置，其中，选择数值最大的位置作为所述下一回溯位置，并且优先选择所述对角线相邻位置；重复确定下一回溯位置，直到所确定的所述下一回溯位置的行号和列号的至少之一为0；比对结果输出模块，所述比对结果输出模块用于基于所确定的回溯路线，确定所述第一核酸序列与所述第二核酸序列的比对结果。
根据本发明实施例的上述测序结果比对系统适于执行上述测序结果比对方法，其附加技术特征以及附加技术特征所带来的优势效果与前面类似，在此不再赘述。
根据本发明实施例的测序结果比对方法和系统以及计算机可读存储介质和电子设备具有以下有益效果：
1)通过采取比平均错配数至少多一个的多个种子序列检索模式，保证了种子序列的匹配，避免了因为种子序列错配而导致匹配不上的情况；
2)通过全局比对，极大地满足了各种类型的小核酸序列的比对，如目的片段比参考序列长、目的片段中存在突变等；
3)提高了小核酸序列数据的利用率，同时提升了其定量的准确性。
附图说明
图1是根据本发明实施例的比对装置的结构示意图；
图2是根据本发明实施例的回溯单元的结构示意图；
图3是根据本发明实施例的测序结果比对系统的结构示意图；
图4是根据本发明实施例的种子序列集合确定装置的结构示意图；以及
图5是根据本发明实施例的比对结果图。
具体实施方式
下面详细描述本发明的实施例，所述实施例的示例在附图中示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
一方面，本发明提出了测序结果的比对方法。根据本发明的实施例，所述方法可描述为：
1)参考序列按设定种子长度进行序列切分(连续位置截取相应长度的片段)，构建索引库；
2)对需要比对的测序读段，记录其5‘端第三个位置到3’端第三个位置的片段(首末端分别允许至多两个滑动)为S，比对时允许的错配数为m；分别从S的m+1等分的起点，截取种子序列长度的片段，如果到S末端的长度小于种子序列长度时，则截取S末端种子序列长度的片段，确定种子序列集合；
3)将种子序列集合的至少一部分与索引库进行匹配；
4)将测序读段与匹配上的参考序列进行Needleman-Wunsh全局比对，具体如下所述：
a)假定参考序列R的长度为n，测序片段S的长度为m，计算得分矩阵Mmn；
其中g为indel罚分-5；
在计算得分矩阵时，记录最大得分值Ss，和相应的位置Sr和Sc；
Mij是记录参考序列和测序片段每个碱基比对时对应的得分矩阵，其中i和j分别是参考序列和测序片段中第i和j位；Ri指参考序列中第i个碱基，Sj指测序片段中第j个碱基，S(Ri，Sj)是根据Ri和Sj的碱基的匹配情况给不同的分值。
打分值与常规的全局比对不一样(常规一般是匹配1，错配0，indel 0)，因为在核酸体内indel(小的插入、缺失)发生的概率是很低的，错配可以理解为发生了突变，而突变频率也是很低的，所以在进行测序读段比对时，遇到错配或者indel应该是小概率事件，所以这样进行罚分。
b)回溯得分矩阵Mmn，从最大得分值Sr和Sc往前回溯，位置Prc的计算规则如下：
Prc是回溯得分矩阵时，下一个回溯位置的的行和列，其中r和c分别代表对应的行和列；第一个判断条件是指如果对角线的值最大，则往对角线回溯；其他两种情况是分别对参考序列和测序片段开gap。
在回溯过程中，记录累计错配数，如果超过设定阈值，则退出比对；
选择最大得分值回溯，可以减少矩阵回溯运算量，同时确保保留最佳比对信息；
回溯过程统计累积错配数，可以减少非最佳参考的比对时间，从而提升软件的运行效率。
c)计算比对得分值，假定匹配数为i，错配数为j，空位数为k，参考序列长度为n，比对得分为：Score＝i-2*j-5*k+i/n；
4)按照比对得分进行排序，输出最高得分的比对结果。
另一方面，本发明提出了一种核酸序列比对装置。根据本发明的实施例，参考图1，所述比对装置包括：获取基本单元信息单元410，所述获取基本单元信息单元410用于获取待比对的第一核酸序列与第二核酸序列各位置上的基本单元信息；构建得分矩阵单元420，所述构建得分矩阵单元420与所述获取基本单元信息单元410相连，用于基于所述基本单元信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列的基本单元数目，n为所述第二核酸序列的基本单元数目，其中所述得分矩阵中的元素Mij表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j个基本单元的比对得分；回溯单元430，所述回溯单元430与所述构建得分矩阵单元420相连，用于基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得经过所述第一核酸序列与所述第二核酸序列的比对结果，其中，
所述元素Mij是基于下列公式确定的：
其中，
Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；
Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；
Mi-1,j表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j个基本单元的比对得分；
g表示小于零的第一预定数值；
S(Ri,Sj)是基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相同时，S(Ri,Sj)为第三预定数值，所述第三预定数值小于所述第二预定数值。
根据本发明的实施例，参考图2，所述回溯单元430进一步包括；确定回溯起始位置模块431，所述确定回溯起始位置模块431用于确定所述矩阵Mmn中的最大值所对应的回溯起始位置；确定下一回溯位置模块432，所述确定下一回溯位置模块432用于基于所述回溯起始位置上游相邻三个位置的数值，确定下一回溯位置，其中，所述上游相邻三个位置包括行相邻位置、对角线相邻位置和列相邻位置，其中，选择数值最大的位置作为所述下一回溯位置，并且优先选择所述对角线相邻位置；重复确定下一回溯位置，直到所确定的所述下一回溯位置的行号和列号的至少之一为0；比对结果输出模块433，所述比对结果输出模块433用于基于所确定的回溯路线，确定所述第一核酸序列与所述第二核酸序列的比对结果。
再一方面，本发明提出了一种测序结果比对系统。根据本发明的实施例，参考图3，所述系统包括：构建索引库装置100，所述构建索引库装置100用于将参考序列按照第一预定长度进行k-mer切分，构建索引库；；种子序列集合确定装置200，所述种子序列集合确定装置200用于基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所述种子序列集合中所述种子序列的数目大于测序的平均错配数，并且所述种子序列的长度不超过所述第一预定长度；匹配装置300，所述匹配装置300用于将所述种子序列集合的至少一部分与所述索引库进行匹配；比对装置400，所述比对装置400用于将所述种子序列集合的至少一部分与所述索引库进行比对，其中所述核酸序列比对装置如前面所限定的。
根据本发明的实施例，参考图4，所述种子序列集合确定装置200进一步包括：末端去除单元210，所述末端去除单元210用于将所述测序读段的5’末端和3’末端各去除第二预定长度的碱基；划分子片段单元220，所述划分子片段单元220用于将经过末端去除单元处理的所述测序读段划分为多个长度相同的子片段，所述子片段的数目超过所述测序的平均错配数；种子序列确定单元230，所述种子序列确定单元230用于在每个子片段中，由5’末端起始基于所述第一预定长度，确定所述种子序列，其中，如果所述子片段的长度小于所述第一预定长度，则将所述子片段作为所述种子序列，如果所述子片段的长度大于所述第一预定长度，则由所述子片段作的5’末端起始延伸所述第一预定长度作为所述种子序列。
下面详细描述本发明的实施例，所述实施例的示例在附图中示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
实施例
1)模拟数据
a)对Human的mature、hairpin和piRNA共28615条进行模拟(其中，Mature、hairpin从miRBase数据库下载，piRNA从piRBase中下载的)，具体规则如下：
b)从参考序列上随机位置截取长度为18～50的片段，如果片段小于50，则最大为片段长度；输出完全匹配片段；
c)截取片段随机一个位置错配，输出错配的片段；截取片段的首末端进行随机(1～2)个碱基的Clip片段，输出clip片段；
d)b)～c)的组合，产生错配和clip的片段；
e)每一条参考序列，重复1～4步骤随机深度次(同一条参考序列，随机生产模拟数据的次数；例如深度2，第一次执行1-4步骤生产了位置为8～20几个片段，第二次执行1-4步骤生成了1-19的几个片段，深度的意思就是对同一个参考序列，进行多少次1-4步骤。10以内)。
序列名称为该片段来源以及相应的位置和匹配模式，如hsa_piR_011099_piRNA_1_27M代表来源于hsa_piR_011099_piRNA的第一个位置，匹配模式为27M，其中M代表匹配，X代表错配，S代表首末端滑动。
2)比对结果
对目前小RNA比对中常用的比对软件，按照不同的参数对1)中模拟数据进行了比对，具体的软件运行参数见表2：
错配数1，种子长度9，线程数1(本方法支持多线程运行)；
默认的种子长度是9，发明人之前测试了7～13，发现9的性价比最高(比对率高和运行时间短)，所以默认选的9。
表2：软件运行参数
其中，sa为本发明对应的软件(https://github.com/zhuqianhua/sa)。
根据各个软件的比对结果，对比其性能，具体结果见图5。
其中，perfectly mapped是要求比对的序列和位置都与真实的一致，而partperfectly mapped是相应软件能比对上就算；不管是总体比对率，还是精确位置的比对，本发明的性能都是优于现有的比对方法的。
需要说明的是，这里的perfectly mapped要求比对起始位置和序列名称正确；在真实数据比较时，Star比对上而本方法比对不上的序列，基本上都是首末端有2个以上碱基滑动,(这个在本申请的方法里面是不认为比对上的，本申请的方法至多允许首末端2个滑动)。这是因为在真实的小RNA分析时，也不太可能说首末端多出3个甚至更多的碱基，还认为是同一个小RNA的片段。
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
在本发明中，除非另有明确的规定和限定，术语“相连”、“连接”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接或彼此可通讯；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
Claims (17)
1.一种测序结果比对方法，其特征在于，包括：
将参考序列按照第一预定长度进行k-mer切分，构建索引库；
基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所述种子序列集合中所述种子序列的数目大于比对所允许的错配数，并且所述种子序列的长度不超过所述第一预定长度；
将所述种子序列集合的至少一部分与所述索引库进行匹配；
将测序读段与匹配上的参考序列进行全局比对，按照比对得分进行排序，输出最高得分的比对结果，
其中，所述全局比对是通过如下方式进行的：
获取待比对的来源于所述测序读段的第一核酸序列与来源于所述匹配上的参考序列的第二核酸序列各位置上的基本单元信息；
基于所述基本单元信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列的基本单元数目，n为所述第二核酸序列的基本单元数目，其中所述得分矩阵中的元素Mij表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j个基本单元的比对得分；
基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得经过所述第一核酸序列与所述第二核酸序列的比对结果，
其中，
所述元素Mij是基于下列公式确定的：
其中，
Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；
Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；
Mi-1,j表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j个基本单元的比对得分；
g表示小于零的第一预定数值；
S(Ri,Sj)是基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相同时，S(Ri,Sj)为第三预定数值，所述第三预定数值小于所述第二预定数值；
所述回溯处理是根据下列步骤确定的；
(a)确定所述矩阵Mmn中的最大值所对应的回溯起始位置；
(b)基于所述回溯起始位置上游相邻三个位置的数值，确定下一回溯位置，其中，所述上游相邻三个位置包括行相邻位置、对角线相邻位置和列相邻位置，其中，选择数值最大的位置作为所述下一回溯位置，并且优先选择所述对角线相邻位置；
(c)重复步骤(b)，直到步骤(b)中所确定的所述下一回溯位置的行号和列号的至少之一为0；
(d)基于步骤(a)-(c)中所确定的回溯路线，确定所述第一核酸序列与所述第二核酸序列的比对结果。
2.根据权利要求1所述的方法，其特征在于，所述第一预定数值为不小于-10的整数。
3.根据权利要求1所述的方法，其特征在于，所述第一预定数值为不小于-5的整数。
4.根据权利要求1所述的方法，其特征在于，所述第二预定数值为1。
5.根据权利要求1所述的方法，其特征在于，所述第三预定数值为-2。
6.根据权利要求1所述的方法，其特征在于，所述基本单元为碱基。
7.根据权利要求1所述的方法，其特征在于，所述第一预定长度为不超过20的整数。
8.根据权利要求1所述的方法，其特征在于，所述第一预定长度为不超过15的整数。
9.根据权利要求1所述的方法，其特征在于，所述第一预定长度为不超过10的整数。
10.根据权利要求1所述的方法，其特征在于，所述第一预定长度为9。
11.根据权利要求1所述的方法，其特征在于，所述种子序列集合是通过下列步骤确定的：
(1)将所述测序读段的5’末端和3’末端各去除第二预定长度的碱基；
(2)将经过步骤(1)处理的所述测序读段划分为多个长度相同的子片段，所述子片段的数目超过所述比对所允许的错配数；
(3)在每个子片段中，由5’末端起始基于所述第一预定长度，确定所述种子序列，其中，如果所述子片段的长度小于所述第一预定长度，则将所述子片段作为所述种子序列，如果所述子片段的长度大于所述第一预定长度，则由所述子片段作的5’末端起始延伸所述第一预定长度作为所述种子序列。
12.根据权利要求11所述的方法，其特征在于，所述第二预定长度为不超过2的整数。
13.根据权利要求11所述的方法，其特征在于，所述子片段的数目＝所述比对所允许的错配数+1。
14.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现如权利要求1-13中任一所述的测序结果比对方法。
15.一种电子设备，其特征在于，包括存储器、处理器；
其中，所述处理器通过读取所述存储器中存储的可执行程序代码来运行与所述可执行程序代码对应的程序，以用于实现如权利要求1-13中任一所述的测序结果比对方法。
16.一种测序结果比对系统，其特征在于，包括：
构建索引库装置，所述构建索引库装置用于将参考序列按照第一预定长度进行k-mer切分，构建索引库；
种子序列集合确定装置，所述种子序列集合确定装置用于基于每条测序读段的序列，确定种子序列集合，所述种子序列集合由多个种子序列构成，所述种子序列集合中所述种子序列的数目大于比对所允许的错配数，并且所述种子序列的长度不超过所述第一预定长度；
匹配装置，所述匹配装置用于将所述种子序列集合的至少一部分与所述索引库进行匹配；
比对装置，所述比对装置用于将测序读段与匹配上的参考序列进行全局比对，按照比对得分进行排序，输出最高得分的比对结果，
其中，
所述比对装置包括：
获取基本单元信息单元，所述获取基本单元信息单元用于获取待比对的来源于所述测序读段的第一核酸序列与来源于所述匹配上的参考序列的第二核酸序列各位置上的基本单元信息；
构建得分矩阵单元，所述构建得分矩阵单元与所述获取基本单元信息单元相连，用于基于所述基本单元信息，构建得分矩阵Mmn，其中，m为所述第一核酸序列的基本单元数目，n为所述第二核酸序列的基本单元数目，其中所述得分矩阵中的元素Mij表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j个基本单元的比对得分；
回溯单元，所述回溯单元与所述构建得分矩阵单元相连，用于基于所述得分矩阵Mmn的数值，进行回溯处理，以便获得经过所述第一核酸序列与所述第二核酸序列的比对结果，
其中，
所述元素Mij是基于下列公式确定的：
其中，
Mi-1,j-1表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；
Mi,j-1表示所述第一核酸序列中第i个基本单元与所述第二核酸序列中第j-1个基本单元的比对得分；
Mi-1,j表示所述第一核酸序列中第i-1个基本单元与所述第二核酸序列中第j个基本单元的比对得分；
g表示小于零的第一预定数值；
S(Ri,Sj)是基于所述第一核酸序列中第i个基本单元Ri与所述第二核酸序列中第j个基本单元Sj确定的数值，其中，当Ri与Sj相同时，S(Ri,Sj)为第二预定数值，当Ri与Sj不相同时，S(Ri,Sj)为第三预定数值，所述第三预定数值小于所述第二预定数值；
所述回溯单元进一步包括；
确定回溯起始位置模块，所述确定回溯起始位置模块用于确定所述矩阵Mmn中的最大值所对应的回溯起始位置；
确定下一回溯位置模块，所述确定下一回溯位置模块用于基于所述回溯起始位置上游相邻三个位置的数值，确定下一回溯位置，其中，所述上游相邻三个位置包括行相邻位置、对角线相邻位置和列相邻位置，其中，选择数值最大的位置作为所述下一回溯位置，并且优先选择所述对角线相邻位置；重复确定下一回溯位置，直到所确定的所述下一回溯位置的行号和列号的至少之一为0；
比对结果输出模块，所述比对结果输出模块用于基于所确定的回溯路线，确定所述第一核酸序列与所述第二核酸序列的比对结果。
17.根据权利要求16所述的系统，其特征在于，所述种子序列集合确定装置进一步包括：
末端去除单元，所述末端去除单元用于将所述测序读段的5’末端和3’末端各去除第二预定长度的碱基；
划分子片段单元，所述划分子片段单元用于将经过末端去除单元处理的所述测序读段划分为多个长度相同的子片段，所述子片段的数目超过所述测序的平均错配数；
种子序列确定单元，所述种子序列确定单元用于在每个子片段中，由5’末端起始基于所述第一预定长度，确定所述种子序列，其中，如果所述子片段的长度小于所述第一预定长度，则将所述子片段作为所述种子序列，如果所述子片段的长度大于所述第一预定长度，则由所述子片段作的5’末端起始延伸所述第一预定长度作为所述种子序列。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911148507.0A CN112825268B (zh) | 2019-11-21 | 2019-11-21 | 测序结果比对方法及其应用 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911148507.0A CN112825268B (zh) | 2019-11-21 | 2019-11-21 | 测序结果比对方法及其应用 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112825268A CN112825268A (zh) | 2021-05-21 |
| CN112825268B true CN112825268B (zh) | 2024-05-14 |
Family
ID=75907369
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911148507.0A Active CN112825268B (zh) | 2019-11-21 | 2019-11-21 | 测序结果比对方法及其应用 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112825268B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114520024B (zh) * | 2022-01-17 | 2024-03-22 | 浙江天科高新技术发展有限公司 | 一种基于k-mer的序列联配方法 |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103793628A (zh) * | 2012-10-29 | 2014-05-14 | 三星Sds株式会社 | 考虑整个短片段的碱基序列比对系统及方法 |
| CN103793625A (zh) * | 2012-10-29 | 2014-05-14 | 三星Sds株式会社 | 碱基序列比对系统及方法 |
| WO2015000284A1 (zh) * | 2013-07-05 | 2015-01-08 | 中国科学院数学与系统科学研究院 | 一种测序序列映射方法及系统 |
| WO2015058120A1 (en) * | 2013-10-18 | 2015-04-23 | Seven Bridges Genomics Inc. | Methods and systems for aligning sequences in the presence of repeating elements |
| CN104794371A (zh) * | 2015-04-29 | 2015-07-22 | 深圳华大基因研究院 | 检测逆转座子插入多态性的方法和装置 |
| CN104881592A (zh) * | 2015-02-11 | 2015-09-02 | 哈尔滨工业大学深圳研究生院 | 一种dna序列比对中的打分方法 |
| CN107403075A (zh) * | 2017-08-02 | 2017-11-28 | 深圳市瀚海基因生物科技有限公司 | 比对方法、装置及系统 |
| CN107844684A (zh) * | 2016-09-18 | 2018-03-27 | 深圳华大基因研究院 | 基因序列比对方法和装置 |
| CN108604260A (zh) * | 2016-01-11 | 2018-09-28 | 艾迪科基因组公司 | 用于现场或基于云的dna和rna处理和分析的基因组学基础架构 |
| CN108866173A (zh) * | 2017-05-16 | 2018-11-23 | 深圳华大基因科技服务有限公司 | 一种标准序列的验证方法、装置及其应用 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9679104B2 (en) * | 2013-01-17 | 2017-06-13 | Edico Genome, Corp. | Bioinformatics systems, apparatuses, and methods executed on an integrated circuit processing platform |
| US10319465B2 (en) * | 2016-11-16 | 2019-06-11 | Seven Bridges Genomics Inc. | Systems and methods for aligning sequences to graph references |
-
2019
- 2019-11-21 CN CN201911148507.0A patent/CN112825268B/zh active Active
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103793628A (zh) * | 2012-10-29 | 2014-05-14 | 三星Sds株式会社 | 考虑整个短片段的碱基序列比对系统及方法 |
| CN103793625A (zh) * | 2012-10-29 | 2014-05-14 | 三星Sds株式会社 | 碱基序列比对系统及方法 |
| WO2015000284A1 (zh) * | 2013-07-05 | 2015-01-08 | 中国科学院数学与系统科学研究院 | 一种测序序列映射方法及系统 |
| WO2015058120A1 (en) * | 2013-10-18 | 2015-04-23 | Seven Bridges Genomics Inc. | Methods and systems for aligning sequences in the presence of repeating elements |
| CN104881592A (zh) * | 2015-02-11 | 2015-09-02 | 哈尔滨工业大学深圳研究生院 | 一种dna序列比对中的打分方法 |
| CN104794371A (zh) * | 2015-04-29 | 2015-07-22 | 深圳华大基因研究院 | 检测逆转座子插入多态性的方法和装置 |
| CN108604260A (zh) * | 2016-01-11 | 2018-09-28 | 艾迪科基因组公司 | 用于现场或基于云的dna和rna处理和分析的基因组学基础架构 |
| CN107844684A (zh) * | 2016-09-18 | 2018-03-27 | 深圳华大基因研究院 | 基因序列比对方法和装置 |
| CN108866173A (zh) * | 2017-05-16 | 2018-11-23 | 深圳华大基因科技服务有限公司 | 一种标准序列的验证方法、装置及其应用 |
| CN107403075A (zh) * | 2017-08-02 | 2017-11-28 | 深圳市瀚海基因生物科技有限公司 | 比对方法、装置及系统 |
Non-Patent Citations (1)
| Title |
|---|
| Zhiyin Yang,Ruibin Zhu,Lina Zhang.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)》.2017,第1316-1320页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112825268A (zh) | 2021-05-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Myers | Efficient local alignment discovery amongst noisy long reads | |
| CN105886616B (zh) | 一种用于猪基因编辑的高效特异性sgRNA识别位点引导序列及其筛选方法 | |
| CA2424031C (en) | System and process for validating, aligning and reordering genetic sequence maps using ordered restriction map | |
| WO2018218788A1 (zh) | 一种基于全局种子打分优选的三代测序序列比对方法 | |
| CN106022002B (zh) | 一种基于三代PacBio测序数据的补洞方法 | |
| JP5279897B2 (ja) | ペア文字列検索システム | |
| CN112825268B (zh) | 测序结果比对方法及其应用 | |
| US20140121983A1 (en) | System and method for aligning genome sequence | |
| US20200243162A1 (en) | Method, system, and computing device for optimizing computing operations of gene sequencing system | |
| KR20160039386A (ko) | Itd 검출 장치 및 방법 | |
| Marić | Long read RNA-seq mapper | |
| EP3663890B1 (en) | Alignment method, device and system | |
| CN114564306B (zh) | 一种基于GPU并行计算的第三代测序RNA-seq比对方法 | |
| US20210202038A1 (en) | Memory Allocation to Optimize Computer Operations of Seeding for Burrows Wheeler Alignment | |
| WO2020182175A1 (en) | Method and system for merging alignment and sorting to optimize | |
| CN115917527A (zh) | 文档检索装置、文档检索系统、文档检索程序、以及文档检索方法 | |
| Šrámek et al. | On-line Viterbi algorithm for analysis of long biological sequences | |
| CN112632343A (zh) | 一种字符串匹配方法、装置、设备及可读存储介质 | |
| CN104239749A (zh) | 碱基序列对准系统及方法 | |
| US20150066384A1 (en) | System and method for aligning genome sequence | |
| CN112825267B (zh) | 确定小核酸序列集合的方法及其应用 | |
| WO2016003283A1 (en) | A method for finding associated positions of bases of a read on a reference genome | |
| WO2020151767A1 (en) | Methods, devices, and systems for performing character matching for short read alignment | |
| CN111326216B (zh) | 一种针对大数据基因测序文件的快速划分方法 | |
| CN115762633B (zh) | 一种基于三代测序的基因组结构变异基因型校正方法 |
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 |