CN105243297A - 一种参考基因组上基因序列片段的快速比对定位方法 - Google Patents
一种参考基因组上基因序列片段的快速比对定位方法 Download PDFInfo
- Publication number
- CN105243297A CN105243297A CN201510648108.6A CN201510648108A CN105243297A CN 105243297 A CN105243297 A CN 105243297A CN 201510648108 A CN201510648108 A CN 201510648108A CN 105243297 A CN105243297 A CN 105243297A
- Authority
- CN
- China
- Prior art keywords
- gene order
- order fragment
- fragment
- key
- database
- 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
Landscapes
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种参考基因组上基因序列片段的快速比对定位方法，步骤包括：从参考基因组上提取基因序列片段；针对每一个基因序列片段，以基因序列片段作为key部分及其目标信息作为Value部分建立键值对，将基因序列片段采用哈希函数映射确定数据库中的目标存储位置，并将键值对写入目标存储位置，完成参考基因组的建库；当需要对待匹配基因序列片段进行快速比对定位时，首先将所述待匹配基因序列片段采用哈希函数映射确定数据库中的目标存储位置，然后从目标存储位置读取匹配基因序列片段对应键值对的基因序列片段的目标信息。本发明具有时间复杂度低、比对定位速度快、定位效率高、应用范围广、能够适用于跨物种的混合快速分析的优点。
Description
技术领域
本发明涉及基因测序数据的生物信息学分析技术，具体涉及一种参考基因组上基因序列片段的快速比对定位方法。
背景技术
随着基因测序技术的发展，测序价格呈指数下降，速度甚至超过摩尔定律。伴随而来的大量测序数据为快速准确的计算分析提出了巨大挑战。分析测序数据的第一步是序列比对，即把序列片段比对定位到参考基因组上，往往需要消耗大量的计算资源与时间。比对定位序列片段正成为一个加速数据分析流程的瓶颈。
为了解决序列片段在参考基因组上的比对定位问题，人们已经发展了很多算法并开发了应用广泛的具体实现，其中比较知名的有BWA，Bowtie，SOAP等软件。现有的比对算法和对软件都是针对在参考基因组上寻找和定位序列片段的最佳比对位置而设计的。因为在比对时考虑了位点突变的容错情况，往往序列片段上的突变越多，在参考基因组上的可能比对上的位置越多，需要的计算量也越大。人类基因组上有大量组成唯一的序列片段。以36个碱基(36bp)的片段为例，30亿个碱基组成的基因组上至少有2/3的区域上的36bp的片段是唯一的。目前的比对算法和软件并不是针对这些唯一片段区域而设计的，它们解决的是一个更普遍的序列比对和定位问题，因此其算法和实现方法针对某类具体问题都存在优化空间。而对于那些不关注突变分析，只关注唯一比对片段的测序流程来说，比对方法的设计有很大的优化空间，通过重新设计，可以节省比对时间，进一步提高效率。
发明内容
本发明要解决的技术问题：针对现有技术的上述问题，提供一种时间复杂度低、比对定位速度快、定位效率高、应用范围广、能够适用于跨物种的混合快速分析的参考基因组上基因序列片段的快速比对定位方法。
为了解决上述技术问题，本发明采用的技术方案为：
一种参考基因组上基因序列片段的快速比对定位方法，步骤包括：
1)从参考基因组上提取基因序列片段；
2)针对参考基因组上的每一个基因序列片段，以基因序列片段作为key部分、基因序列片段的目标信息作为Value部分建立键值对，将基因序列片段采用哈希函数映射确定数据库中的目标存储位置，并将键值对写入目标存储位置，最终完成参考基因组的建库；
3)当需要对待匹配基因序列片段进行快速比对定位时，首先将所述待匹配基因序列片段采用哈希函数映射查找数据库中的目标存储位置，如果查找成功，则从目标存储位置读取匹配基因序列片段对应键值对的基因序列片段的目标信息；否则返回查找失败信息。
优选地，所述步骤1)的详细步骤包括：
1.1)设定基因序列片段的长度L；
1.2)计算参考基因组上唯一的基因序列片段的位置和目标信息；
1.3)根据基因序列片段的位置分别提取基因序列片段及其目标信息。
优选地，所述步骤1.2)中唯一的基因序列片段具体是指任意两个基因序列片段之间的编辑距离大于等于设定的阈值n。
优选地，所述目标信息包括染色体、染色***置、GC含量、物种分类中的至少一种。
优选地，所述步骤2)的详细步骤包括：
2.1)从提取得到的所有基因序列片段中取出一个基因序列片段作为当前基因序列片段；
2.2)以当前基因序列片段作为key部分、当前基因序列片段的目标信息作为Value部分建立键值对(Key,Value)以描述当前基因序列片段及其目标信息之间的映射关系；
2.3)将当前基因序列片段的键值对(Key,Value)中的Key部分和Value部分进行编码，选用指定的哈希函数将当前基因序列片段映射到d个数据库中的数据库i中；
2.4)将当前基因序列片段采用哈希函数映射确定其在数据库i中的目标存储位置，并将编码后的键值对(Key,Value)写入数据库i的目标存储位置；
2.5)判断提取得到的所有基因序列片段是否已经处理完毕，如果已经处理完毕，则跳转执行步骤2.1)；否则，判定已经完成参考基因组的建库。
优选地，所述步骤2.3)中将当前基因序列片段映射到d个数据库中的数据库i中的详细步骤包括：
2.3.1)设定数据库个数d；
2.3.2)取当前基因序列片段Key部分长度为m的前缀子串，选用指定的哈希函数并采用如式(1)所示函数式计算当前基因序列片段在d个数据库中对应的数据库编号i，从而将当前基因序列片段映射到d个数据库中的数据库i中；
i＝f(Key[1:m])％d(1)
式(1)中，i为当前基因序列片段在d个数据库中对应的数据库编号，Key[1:m]为当前基因序列片段Key部分长度为m的前缀子串，d为预先设定的数据库个数，f为指定的哈希函数。
优选地，所述步骤3)的详细步骤包括：
3.1)当需要对待匹配基因序列片段进行快速比对定位时，首先计算待匹配基因序列片段的反向互补序列，将待匹配基因序列片段及其反向互补序列分别进行编码，然后将编码后长度为m的前缀子串选用指定的哈希函数并采用式(1)所示函数式分别计算待匹配基因序列片段及其反向互补序列在d个数据库中对应的数据库编号i；
3.2)针对待匹配基因序列片段及其反向互补序列，分别在计算得到的数据库编号i对应目标数据库中并行查找待匹配基因序列片段及其反向互补序列，如果待匹配基因序列片段及其反向互补序列两者中的任意一个在目标数据库中找到匹配记录，则将匹配键值对中的Value部分的目标信息作为比对定位的结果返回；否则待匹配基因序列片段及其反向互补序列两者在目标数据库中均未找到匹配记录，则返回查找失败信息。
优选地，所述步骤3.1)中计算待匹配基因序列片段的反向互补序列具体是指首先将待匹配基因序列片段反转，然后分别将反转后的基因序列片段中碱基A和碱基T互换、碱基C和碱基G互换，得到待匹配基因序列片段的反向互补序列。
优选地，所述步骤2.3)将当前基因序列片段的键值对(Key,Value)中的Key部分和Value部分进行编码具体是指：针对Key部分依次以00替换碱基A、以01替换碱基C、以10替换碱基G、以11替换碱基T进行编码，针对Value部分采用可变长编码；所述步骤3.1)将所述待匹配基因序列片段进行编码具体是指：依次以00替换碱基A、以01替换碱基C、以10替换碱基G、以11替换碱基T进行编码。
优选地，所述指定的哈希函数为MurmurHash函数。
本发明参考基因组上基因序列片段的快速比对定位方法具有下述优点：
1、本发明针对参考基因组上的每一个基因序列片段，以基因序列片段作为key部分、基因序列片段的目标信息作为Value部分建立键值对，将基因序列片段采用哈希函数映射确定数据库中的目标存储位置，用键值对的搜索代替比对算法，快速获得序列片段的比对结果(即是否比对上，如果比对上则可以获取其目标信息)，键值对搜索具有最低的时间复杂度，与现有序列比对方法相比，具有最快的比对定位速度，具有时间复杂度低、比对定位速度快、定位效率高的优点。
2、本发明针对参考基因组上的每一个基因序列片段，以基因序列片段作为key部分、基因序列片段的目标信息作为Value部分建立键值对，将基因序列片段采用哈希函数映射确定数据库中的目标存储位置，因此可以通过将多个基因组的唯一独特序列片段放在一起的方式，同时搜索比对多个基因组，可以将不同物种基因组上的唯一比对区域放到一起，同时在多物种基因组上比对定位序列片段，适合跨物种的混合快速分析。
3、对于大多数物种的参考基因组来说，其上的大部分区域里的片段都是独一无二的，正因为如此，本发明能够显著加速那些不关注非唯一比对序列片段的分析流程。
4、本发明以基因序列片段作为key部分、基因序列片段的目标信息作为Value部分建立键值对，基于目标信息的具体内容不同，能够应用于各类基因序列片段快速比对定位的基因分析应用中，例如CNV分析、菌种检定等，具有应用范围广的优点。
附图说明
图1为本发明实施例一中快速比对定位方法的基本流程示意图。
图2为本发明实施例一CNV分析得到的CNV分析结果示意图。
具体实施方式
实施例一：
下文以通过快速序列比对进行CNV分析(CopyNumberVariation分析/拷贝数变异分析)为例，对本发明参考基因组上基因序列片段的快速比对定位方法进行进一步的说明。
如图1所示，本实施例参考基因组上基因序列片段的快速比对定位方法的步骤包括：
1)从参考基因组上提取基因序列片段；
2)针对参考基因组上的每一个基因序列片段，以基因序列片段作为key部分、基因序列片段的目标信息作为Value部分建立键值对，将基因序列片段采用哈希函数映射确定数据库中的目标存储位置，并将键值对写入目标存储位置，最终完成参考基因组的建库；
3)当需要对待匹配基因序列片段进行快速比对定位时，将所述待匹配基因序列片段采用哈希函数映射查找数据库中的目标存储位置，如果查找成功，则从目标存储位置读取匹配基因序列片段对应键值对的基因序列片段的目标信息；否则返回查找失败信息。
本实施例中，参考基因组来自在UCSC网站的人类基因组(版本hg19)36bp唯一基因序列片段的位置信息来自网址：
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeCrgMapabilityAlign36mer.bigWig
该文件为人类基因组(版本hg19)36mermappability文件，将上述文件下载后使用工具bigWigToBedGraph将bigWig文件转换为BED格式文件，BED格式文件包括四列，第1列为染色体编号，第2列为染色体起始位置(以0开始计算)，第3列为染色体终止位置(以1开始计算)，第4列为染色体起始位置至染色体终止位置(不包括终止位置)区间中每个位置开始计算的36bp序列的mappability值,mappabilit为1表示基因序列片段在基因组上唯一，本实施例中在BED格式文件中选取mappability为1的记录，并根据BED文件提供的位置信息提取出hg19参考基因组上的唯一基因序列片段，提取出的目标基因序列片段集合将作为构建数据库所使用Key的集合，各个染色体上的唯一覆盖度如表1所示。
表1：各个染色体上的唯一覆盖度数据表。
| 染色体 | 碱基数目 | 唯一比对碱基数目 | 唯一覆盖度 |
| 1 | 249,250,621 | 171,109,348 | 0.69 |
| 2 | 243,199,373 | 187,599,789 | 0.77 |
| 3 | 198,022,430 | 155,568,147 | 0.79 |
| 4 | 191,154,276 | 148,235,325 | 0.78 |
| 5 | 180,915,260 | 139,351,070 | 0.77 |
| 6 | 171,115,067 | 132,615,985 | 0.78 |
| 7 | 159,138,663 | 115,076,932 | 0.72 |
| 8 | 146,364,022 | 113,561,421 | 0.78 |
| 9 | 141,213,431 | 85,976,834 | 0.61 |
| 10 | 135,534,747 | 100,580,972 | 0.74 |
| 11 | 135,006,516 | 101,426,392 | 0.75 |
| 12 | 133,851,895 | 101,134,431 | 0.76 |
| 13 | 115,169,878 | 77,517,018 | 0.67 |
| 14 | 107,349,540 | 68,732,800 | 0.64 |
| 15 | 102,531,392 | 59,987,912 | 0.59 |
| 16 | 90,354,753 | 56,576,358 | 0.63 |
| 17 | 81,195,210 | 55,582,856 | 0.68 |
| 18 | 78,077,248 | 60,876,802 | 0.78 |
| 19 | 59,128,983 | 35,479,931 | 0.6 |
| 20 | 63,025,520 | 47,555,303 | 0.75 |
| 21 | 48,129,895 | 26,757,998 | 0.56 |
| 22 | 51,304,566 | 24,675,505 | 0.48 |
| X | 155,270,560 | 103,355,975 | 0.67 |
| Y | 59,373,566 | 7,016,301 | 0.12 |
参见表1，第一列为染色体编号，第二列为每个染色体包含的碱基数目，第三列为唯一基因序列片段所包含的碱基数目，第四列为唯一基因序列片段所包含的碱基数目占每个染色体碱基数目的比例，本实施例中，36bp的唯一基因序列片段总计2,176,351,405，碱基总数占参考基因组的70.3％，染色体的平均覆盖度为67.1％。
本实施例中，步骤1)的详细步骤包括：
1.1)设定基因序列片段的长度L；本实施例中，基因序列片段的长度L为36；
1.2)计算参考基因组上唯一的基因序列片段的位置和目标信息；
1.3)根据基因序列片段的位置分别提取基因序列片段及其目标信息。
本实施例中，步骤1.2)中唯一的基因序列片段具体是指任意两个基因序列片段之间的编辑距离大于等于设定的阈值n，阈值n具体取值为2，此外也可以根据需要设置为其他值。
需要说明的是，目标信息具体和基因序列片段的快速比对定位应用相关，目标信息可以包括染色体、染色***置、GC含量、物种分类中的至少一种。本实施例中，目标信息具体包括：(1)基因序列片段的染色体编号{1,2,...,22,X,Y}；(2)基因序列片段在染色体上的起始位置(从0开始计算)；(3)基因序列片段的GC含量(即基因序列片段中碱基G、碱基C的总数目)。
本实施例中，步骤2)的详细步骤包括：
2.1)从提取得到的所有基因序列片段中取出一个基因序列片段作为当前基因序列片段；
2.2)以当前基因序列片段作为key部分、当前基因序列片段的目标信息作为Value部分建立键值对(Key,Value)以描述当前基因序列片段及其目标信息之间的映射关系；
2.3)将当前基因序列片段的键值对(Key,Value)中的Key部分和Value部分进行编码，选用指定的哈希函数将当前基因序列片段映射到d个数据库中的数据库i(1≤i≤d)中；
2.4)将当前基因序列片段采用哈希函数映射确定其在数据库i中的目标存储位置，并将编码后的键值对(Key,Value)写入数据库i的目标存储位置；
2.5)判断提取得到的所有基因序列片段是否已经处理完毕，如果已经处理完毕，则跳转执行步骤2.1)；否则，判定已经完成参考基因组的建库。
本实施例中，步骤2.3)中将当前基因序列片段映射到d个数据库中的数据库i中的详细步骤包括：
2.3.1)设定数据库个数d；本实施例中数据库个数d具体取值为12；
2.3.2)取当前基因序列片段Key部分长度为m的前缀子串，选用指定的哈希函数并采用如式(1)所示函数式计算当前基因序列片段在d个数据库中对应的数据库编号i，从而将当前基因序列片段映射到d个数据库中的数据库i中；
i＝f(Key[1:m])％d(1)
式(1)中，i为当前基因序列片段在d个数据库中对应的数据库编号，Key[1:m]为当前基因序列片段Key部分长度为m的前缀子串，d为预先设定的数据库个数，f为指定的哈希函数。
指定的哈希函数f用于实现将所有基因序列片段映射到d个数据库中，本实施例中哈希函数f具体采用MurmurHash函数，能够将所有基因序列片段均匀地映射到d个数据库中。
本实施例中，步骤3)的详细步骤包括：
3.1)当需要对待匹配基因序列片段进行快速比对定位时，首先计算待匹配基因序列片段的反向互补序列，将待匹配基因序列片段及其反向互补序列分别进行编码，然后将编码后长度为m的前缀子串选用指定的哈希函数并采用式(1)所示函数式分别计算待匹配基因序列片段及其反向互补序列在d个数据库中对应的数据库编号i(参考基因组的建库时确定的编号)；本实施例中，前缀子串m具体取值为3；
3.2)针对待匹配基因序列片段及其反向互补序列，分别在计算得到的数据库编号i对应目标数据库中并行查找待匹配基因序列片段及其反向互补序列，如果待匹配基因序列片段及其反向互补序列两者中的任意一个在目标数据库中找到匹配记录，则将匹配键值对中的Value部分的目标信息作为比对定位的结果返回；否则待匹配基因序列片段及其反向互补序列两者在目标数据库中均未找到匹配记录，则返回查找失败信息。
本实施例中，步骤3.1)中计算待匹配基因序列片段的反向互补序列具体是指首先将待匹配基因序列片段反转，然后分别将反转后的基因序列片段中碱基A和碱基T互换、碱基C和碱基G互换，得到待匹配基因序列片段的反向互补序列。例如，对于36位的待匹配基因序列片段AC…ACGT而言，反转后得到TGCA…CA，再进行碱基A和碱基T互换、碱基C和碱基G互换后，得到ACGT…GT。
本实施例中，步骤2.3)将当前基因序列片段的键值对中的Key部分和Value部分进行编码时，针对当前基因序列片段的键值对中的Key部分依次以00替换碱基A、以01替换碱基C、以10替换碱基G、以11替换碱基T进行编码，针对当前基因序列片段的键值对中的Value部分采用可变长编码；步骤3.1)将待匹配基因序列片段进行编码时，具体是指依次以00替换碱基A、以01替换碱基C、以10替换碱基G、以11替换碱基T进行编码；本实施例中，步骤2.3)将当前基因序列片段的键值对中的Key部分和Value部分进行编码时，编码字符采用ASCII码，编码后Key部分的长度为9，Value部分的长度介于3～6个字符之间；步骤3.1)将待匹配基因序列片段进行编码时，编码后Key部分的长度为9。步骤2.3)针对Value部分采用可变长编码时，编码的顺序依次为染色体编号、基因序列片段的GC含量、基因序列片段在染色体上的起始位置，染色体编号、基因序列片段的GC含量各占一个字符，基因序列片段在染色体上的起始位置部分采用可变长编码，即当位置位于区间[0,28-1)时用一个字符编码，当位置位于区间[28,216-1)时用两个字符编码，当位置位于区间[216,224-1)时用三个字符编码，当位置位于区间[224,232-1)时用四个字符编码，从而能够节约存储空间的占用。
本实施例中，待匹配基因序列片段数据集为真实人类血液样本的测序结果，测序片段长度为36bp，序列片段的数目为4,929,709，用本发明方法进行快速序列比对定位并统计比对结果如下表2所示。
表2：各个染色体上的序列数据表。
| 染色体 | 比对序列数目 | 唯一比对序列数目 | 唯一比对序列GC含量 |
| 1 | 279880 | 277196 | 0.40525 |
| 2 | 299817 | 296786 | 0.390581 |
| 3 | 248322 | 245823 | 0.385119 |
| 4 | 234336 | 231968 | 0.369009 |
| 5 | 220284 | 218086 | 0.383703 |
| 6 | 213813 | 211737 | 0.384034 |
| 7 | 185199 | 183429 | 0.392871 |
| 8 | 180296 | 178557 | 0.389799 |
| 9 | 137702 | 136404 | 0.403912 |
| 10 | 164105 | 162531 | 0.404724 |
| 11 | 161240 | 159706 | 0.407018 |
| 12 | 163420 | 161837 | 0.394989 |
| 13 | 122829 | 121550 | 0.372917 |
| 14 | 110550 | 109441 | 0.397038 |
| 15 | 99286 | 98357 | 0.409776 |
| 16 | 91390 | 90550 | 0.437317 |
| 17 | 91902 | 91116 | 0.444336 |
| 18 | 95782 | 94852 | 0.387118 |
| 19 | 57913 | 57502 | 0.477258 |
| 20 | 76343 | 75659 | 0.432834 |
| 21 | 44167 | 43733 | 0.402295 |
| 22 | 40754 | 40440 | 0.47358 |
| X | 164658 | 163060 | 0.382272 |
| Y | 481 | 476 | 0.394771 |
参见表2，第一列为染色体编号，第二列为匹配数据库基因序列片段的数目，第三列为匹配数据库中基因序列片段中唯一基因序列片段的数目，第四列为匹配数据库基因序列片段中唯一基因序列片段的GC含量，具体是指匹配数据库基因序列片段中碱基C和G所占的比例。本实施例中匹配数据库的基因序列片段有3,484,469条，唯一匹配的基因序列片段有3,450,796条，唯一匹配基因序列片段占待匹配基因序列片段总数目的比例为70％，与表1中所述的唯一基因序列片段碱基数目占参考基因组总碱基数目的比例一致。
在得到表2的基础上，标准化比对到各染色体的序列并进行GC校正和mappability校正，选取正常的样本做参照分析染色体拷贝数变异，最终在2号染色体q(Chr2)端准确发现已知CNV突变，CNV分析结果如附图2所示。参见图2，x坐标为2号染色体上以5K距离间隔滑动的60K窗口的位置；y坐标为log2ratio，log2ratio显著偏离0，表示本实施例样本在对应染色***置存在拷贝数变异。y坐标log2ratio的具体计算方法如式(2)所示；
式(2)中，x60K为本实施例样本数据中对唯一比对到2号染色体每一个60K窗口内的基因序列片段数目作标准化后的值，为正常样本数据中对唯一比对到2号染色体每一个60K窗口内的基因序列片段数目作标准化后的均值。
以本实施例36bp序列为例，快速比对算法的比对效率为1Million序列/s，与传统的基于BWT的通用型比对软件相比，省掉了耗时的SmithWaterman比对过程，比对效率大幅提高。
实施例二：
除了进行CNV分析以外，本发明参考基因组上基因序列片段的快速比对定位方法还可以用于进行菌种鉴定。本实施例中，参考基因组为跨物种参考基因组，步骤2)进行参考基因组的建库时，参考基因组为在NCBIRefSeq数据库中(http://www.ncbi.nlm.nih.gov/refseq/)下载的人参考基因组和所有细菌、病毒基因组序列。用gem-mappability软件计算(人、细菌、病毒)参考序列31mer的mappability图谱，从中选取mappability为1的唯一序列片段；在NCBI物种分类数据库(http://www.ncbi.nlm.nih.gov/taxonomy)下载物种的详细分类信息。针对参考基因组上的每一个基因序列片段，以基因序列片段作为key部分、基因序列片段的目标信息(物种分类)作为Value部分建立键值对，将基因序列片段采用哈希函数映射确定数据库中的目标存储位置，并将键值对写入目标存储位置，最终完成参考基因组的建库。在HMP网站(http://hmpdacc.org/)下载正常人粪便样本的测序数据来分析人类肠道微生物的组成，针对其中的待匹配基因序列片段使用步骤3)进行快速序列比对分析在属水平上的微生物组成，最终得到的结果如表3所示。
表3：微生物组成分析结果数据表。
| 属 | 分类序列数目 | 标准化分类序列数目 | 百分比 |
| Bacteroides | 3345216 | 153556.7 | 15.36 |
| Eubacterium | 1597681 | 73338.94 | 7.33 |
| Roseburia | 224359 | 10298.83 | 1.03 |
| Alistipes | 197067 | 9046.04 | 0.9 |
| Odoribacter | 174594 | 8014.45 | 0.8 |
| Others | 505842 | 23219.85 | 2.32 |
参见表3，第一列为属水平的微生物名称，第二列为比对到相应微生物属水平的基因序列片段的数目，第三列为每1,000,000序列片段中比对到相应微生物属水平的序列片段的数目，第四列为标准化后的分类序列数目占1,000,000序列片段的百分比，用来衡量样本中微生物的丰度，该值越高表明样本中对应微生物的组成比例越高。最终，根据定位到不同菌种的唯一序列片段来判断所分析样本中的菌种丰度，该分析结果与HMP网站给出的分析结果一致。本实施例参考基因组上基因序列片段的快速比对定位方法可针对各种具体的应用对参考基因组、目标信息进行高度化定制，灵活多变；对参考基因组的定制和简化极大提高了序列比对的效率；本实施例参考基因组上基因序列片段的快速比对定位方法在加速比对分析的同时，并未牺牲分析准确性，极其适用于大规模海量数据分析和对实时性要求较高的应用场景，应用前景广阔。
以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
Claims (10)
1.一种参考基因组上基因序列片段的快速比对定位方法，其特征在于步骤包括：
1)从参考基因组上提取基因序列片段；
2)针对参考基因组上的每一个基因序列片段，以基因序列片段作为key部分、基因序列片段的目标信息作为Value部分建立键值对，将基因序列片段采用哈希函数映射确定数据库中的目标存储位置，并将键值对写入目标存储位置，最终完成参考基因组的建库；
3)当需要对待匹配基因序列片段进行快速比对定位时，首先将所述待匹配基因序列片段采用哈希函数映射查找数据库中的目标存储位置，如果查找成功，则从目标存储位置读取匹配基因序列片段对应键值对的基因序列片段的目标信息；否则返回查找失败信息。
2.根据权利要求1所述的参考基因组上基因序列片段的快速比对定位方法，其特征在于，所述步骤1)的详细步骤包括：
1.1)设定基因序列片段的长度L；
1.2)计算参考基因组上唯一的基因序列片段的位置和目标信息；
1.3)根据基因序列片段的位置分别提取基因序列片段及其目标信息。
3.根据权利要求2所述的参考基因组上基因序列片段的快速比对定位方法，其特征在于，所述步骤1.2)中唯一的基因序列片段具体是指任意两个基因序列片段之间的编辑距离大于等于设定的阈值n。
4.根据权利要求3所述的参考基因组上基因序列片段的快速比对定位方法，其特征在于，所述目标信息包括染色体、染色***置、GC含量、物种分类中的至少一种。
5.根据权利要求1～4中任意一项所述的参考基因组上基因序列片段的快速比对定位方法，其特征在于，所述步骤2)的详细步骤包括：
2.1)从提取得到的所有基因序列片段中取出一个基因序列片段作为当前基因序列片段；
2.2)以当前基因序列片段作为key部分、当前基因序列片段的目标信息作为Value部分建立键值对(Key,Value)以描述当前基因序列片段及其目标信息之间的映射关系；
2.3)将当前基因序列片段的键值对(Key,Value)中的Key部分和Value部分进行编码，选用指定的哈希函数将当前基因序列片段映射到d个数据库中的数据库i中；
2.4)将当前基因序列片段采用哈希函数映射确定其在数据库i中的目标存储位置，并将编码后的键值对(Key,Value)写入数据库i的目标存储位置；
2.5)判断提取得到的所有基因序列片段是否已经处理完毕，如果已经处理完毕，则跳转执行步骤2.1)；否则，判定已经完成参考基因组的建库。
6.根据权利要求5所述的参考基因组上基因序列片段的快速比对定位方法，其特征在于，所述步骤2.3)中将当前基因序列片段映射到d个数据库中的数据库i中的详细步骤包括：
2.3.1)设定数据库个数d；
2.3.2)取当前基因序列片段Key部分长度为m的前缀子串，选用指定的哈希函数并采用如式(1)所示函数式计算当前基因序列片段在d个数据库中对应的数据库编号i，从而将当前基因序列片段映射到d个数据库中的数据库i中；
i＝f(Key[1:m])％d(1)
式(1)中，i为当前基因序列片段在d个数据库中对应的数据库编号，Key[1:m]为当前基因序列片段Key部分长度为m的前缀子串，d为预先设定的数据库个数，f为指定的哈希函数。
7.根据权利要求6所述的参考基因组上基因序列片段的快速比对定位方法，其特征在于，所述步骤3)的详细步骤包括：
3.1)当需要对待匹配基因序列片段进行快速比对定位时，首先计算待匹配基因序列片段的反向互补序列，将待匹配基因序列片段及其反向互补序列分别进行编码，然后将编码后长度为m的前缀子串选用指定的哈希函数并采用式(1)所示函数式分别计算待匹配基因序列片段及其反向互补序列在d个数据库中对应的数据库编号i；
3.2)针对待匹配基因序列片段及其反向互补序列，分别在计算得到的数据库编号i对应目标数据库中并行查找待匹配基因序列片段及其反向互补序列，如果待匹配基因序列片段及其反向互补序列两者中的任意一个在目标数据库中找到匹配记录，则将匹配键值对中的Value部分的目标信息作为比对定位的结果返回；否则待匹配基因序列片段及其反向互补序列两者在目标数据库中均未找到匹配记录，则返回查找失败信息。
8.根据权利要求7所述的参考基因组上基因序列片段的快速比对定位方法，其特征在于，所述步骤3.1)中计算待匹配基因序列片段的反向互补序列具体是指首先将待匹配基因序列片段反转，然后分别将反转后的基因序列片段中碱基A和碱基T互换、碱基C和碱基G互换，得到待匹配基因序列片段的反向互补序列。
9.根据权利要求8所述的参考基因组上基因序列片段的快速比对定位方法，其特征在于，所述步骤2.3)将当前基因序列片段的键值对(Key,Value)中的Key部分和Value部分进行编码具体是指：针对Key部分依次以00替换碱基A、以01替换碱基C、以10替换碱基G、以11替换碱基T进行编码，针对Value部分采用可变长编码；所述步骤3.1)将所述待匹配基因序列片段进行编码具体是指：依次以00替换碱基A、以01替换碱基C、以10替换碱基G、以11替换碱基T进行编码。
10.根据权利要求9所述的参考基因组上基因序列片段的快速比对定位方法，其特征在于，所述指定的哈希函数为MurmurHash函数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510648108.6A CN105243297A (zh) | 2015-10-09 | 2015-10-09 | 一种参考基因组上基因序列片段的快速比对定位方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510648108.6A CN105243297A (zh) | 2015-10-09 | 2015-10-09 | 一种参考基因组上基因序列片段的快速比对定位方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN105243297A true CN105243297A (zh) | 2016-01-13 |
Family
ID=55040943
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510648108.6A Pending CN105243297A (zh) | 2015-10-09 | 2015-10-09 | 一种参考基因组上基因序列片段的快速比对定位方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105243297A (zh) |
Cited By (16)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106295250A (zh) * | 2016-07-28 | 2017-01-04 | 北京百迈客医学检验所有限公司 | 二代测序短序列快速比对分析方法及装置 |
| CN107180166A (zh) * | 2017-04-21 | 2017-09-19 | 北京希望组生物科技有限公司 | 一种基于三代测序的全基因组结构变异分析方法和*** |
| WO2019023978A1 (zh) * | 2017-08-02 | 2019-02-07 | 深圳市瀚海基因生物科技有限公司 | 比对方法、装置及*** |
| CN109584965A (zh) * | 2018-12-21 | 2019-04-05 | 龙口味美思环保科技有限公司 | 一种转基因食品的交叉式基因组重组排序方法 |
| CN109698011A (zh) * | 2018-12-25 | 2019-04-30 | 人和未来生物科技（长沙）有限公司 | 基于短序列比对的Indel区域校正方法及*** |
| CN110021368A (zh) * | 2017-10-20 | 2019-07-16 | 人和未来生物科技（长沙）有限公司 | 比对型基因测序数据压缩方法、***及计算机可读介质 |
| CN110060731A (zh) * | 2019-04-12 | 2019-07-26 | 福建师范大学 | 基于分布式计算确定基因组间重叠基因对数量的方法 |
| CN110517728A (zh) * | 2019-08-29 | 2019-11-29 | 苏州浪潮智能科技有限公司 | 一种基因序列比对方法及装置 |
| CN110782946A (zh) * | 2019-10-17 | 2020-02-11 | 南京医基云医疗数据研究院有限公司 | 识别重复序列的方法及装置、存储介质、电子设备 |
| CN110875084A (zh) * | 2018-08-13 | 2020-03-10 | 深圳华大基因科技服务有限公司 | 一种核酸序列比对的方法 |
| CN111063394A (zh) * | 2019-12-13 | 2020-04-24 | 人和未来生物科技（长沙）有限公司 | 基于基因序列的物种快速查找及建库方法、***和介质 |
| CN112700819A (zh) * | 2020-12-31 | 2021-04-23 | 云舟生物科技（广州）有限公司 | 基因序列的处理方法、计算机存储介质及电子设备 |
| CN113012755A (zh) * | 2021-04-12 | 2021-06-22 | 聊城大学 | 基因组atcg的检索方法 |
| WO2022048284A1 (zh) * | 2020-09-04 | 2022-03-10 | 苏州浪潮智能科技有限公司 | 一种基因对比的哈希查表方法、装置、设备及存储介质 |
| CN115862740A (zh) * | 2022-12-06 | 2023-03-28 | 中国人民解放军军事科学院军事医学研究院 | 一种面向大规模病毒基因组数据的快速分布式多序列比对方法 |
| CN115862735A (zh) * | 2022-12-28 | 2023-03-28 | 郑州思昆生物工程有限公司 | 一种核酸序列检测方法、装置、计算机设备及存储介质 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103971031A (zh) * | 2014-05-04 | 2014-08-06 | 南京师范大学 | 一种面向大规模基因数据的读段定位方法 |
| US20150234842A1 (en) * | 2014-02-14 | 2015-08-20 | Srinivasan Kumar | Mapping of Extensible Datasets to Relational Database Schemas |
-
2015
- 2015-10-09 CN CN201510648108.6A patent/CN105243297A/zh active Pending
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150234842A1 (en) * | 2014-02-14 | 2015-08-20 | Srinivasan Kumar | Mapping of Extensible Datasets to Relational Database Schemas |
| CN103971031A (zh) * | 2014-05-04 | 2014-08-06 | 南京师范大学 | 一种面向大规模基因数据的读段定位方法 |
Non-Patent Citations (3)
| Title |
|---|
| 危扬: "基于光模型的高产作物株型定量化设计研究", 《中国优秀硕士学位论文全文数据库 农业科技辑》 * |
| 李航: "基于MPI的并行DNA序列比对***的设计与实现", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
| 王立坤: "RNA-seq数据的处理与应用", 《中国博士学位论文全文数据库 基础科学辑》 * |
Cited By (27)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106295250A (zh) * | 2016-07-28 | 2017-01-04 | 北京百迈客医学检验所有限公司 | 二代测序短序列快速比对分析方法及装置 |
| CN106295250B (zh) * | 2016-07-28 | 2019-03-29 | 北京百迈客医学检验所有限公司 | 二代测序短序列快速比对分析方法及装置 |
| CN107180166A (zh) * | 2017-04-21 | 2017-09-19 | 北京希望组生物科技有限公司 | 一种基于三代测序的全基因组结构变异分析方法和*** |
| WO2019023978A1 (zh) * | 2017-08-02 | 2019-02-07 | 深圳市瀚海基因生物科技有限公司 | 比对方法、装置及*** |
| US11482304B2 (en) | 2017-08-02 | 2022-10-25 | Genemind Biosciences Company Limited | Alignment methods, devices and systems |
| CN110021368B (zh) * | 2017-10-20 | 2020-07-17 | 人和未来生物科技（长沙）有限公司 | 比对型基因测序数据压缩方法、***及计算机可读介质 |
| CN110021368A (zh) * | 2017-10-20 | 2019-07-16 | 人和未来生物科技（长沙）有限公司 | 比对型基因测序数据压缩方法、***及计算机可读介质 |
| CN110875084A (zh) * | 2018-08-13 | 2020-03-10 | 深圳华大基因科技服务有限公司 | 一种核酸序列比对的方法 |
| CN110875084B (zh) * | 2018-08-13 | 2022-06-21 | 深圳华大基因科技服务有限公司 | 一种核酸序列比对的方法 |
| CN109584965A (zh) * | 2018-12-21 | 2019-04-05 | 龙口味美思环保科技有限公司 | 一种转基因食品的交叉式基因组重组排序方法 |
| CN109698011A (zh) * | 2018-12-25 | 2019-04-30 | 人和未来生物科技（长沙）有限公司 | 基于短序列比对的Indel区域校正方法及*** |
| CN109698011B (zh) * | 2018-12-25 | 2020-10-23 | 人和未来生物科技（长沙）有限公司 | 基于短序列比对的Indel区域校正方法及*** |
| CN110060731A (zh) * | 2019-04-12 | 2019-07-26 | 福建师范大学 | 基于分布式计算确定基因组间重叠基因对数量的方法 |
| CN110517728A (zh) * | 2019-08-29 | 2019-11-29 | 苏州浪潮智能科技有限公司 | 一种基因序列比对方法及装置 |
| CN110517728B (zh) * | 2019-08-29 | 2022-04-29 | 苏州浪潮智能科技有限公司 | 一种基因序列比对方法及装置 |
| CN110782946A (zh) * | 2019-10-17 | 2020-02-11 | 南京医基云医疗数据研究院有限公司 | 识别重复序列的方法及装置、存储介质、电子设备 |
| CN111063394A (zh) * | 2019-12-13 | 2020-04-24 | 人和未来生物科技（长沙）有限公司 | 基于基因序列的物种快速查找及建库方法、***和介质 |
| CN111063394B (zh) * | 2019-12-13 | 2023-07-11 | 人和未来生物科技（长沙）有限公司 | 基于基因序列的物种快速查找及建库方法、***和介质 |
| WO2022048284A1 (zh) * | 2020-09-04 | 2022-03-10 | 苏州浪潮智能科技有限公司 | 一种基因对比的哈希查表方法、装置、设备及存储介质 |
| CN112700819B (zh) * | 2020-12-31 | 2021-11-30 | 云舟生物科技（广州）有限公司 | 基因序列的处理方法、计算机存储介质及电子设备 |
| CN112700819A (zh) * | 2020-12-31 | 2021-04-23 | 云舟生物科技（广州）有限公司 | 基因序列的处理方法、计算机存储介质及电子设备 |
| CN113012755A (zh) * | 2021-04-12 | 2021-06-22 | 聊城大学 | 基因组atcg的检索方法 |
| CN113012755B (zh) * | 2021-04-12 | 2023-10-27 | 聊城大学 | 基因组atcg的检索方法 |
| CN115862740A (zh) * | 2022-12-06 | 2023-03-28 | 中国人民解放军军事科学院军事医学研究院 | 一种面向大规模病毒基因组数据的快速分布式多序列比对方法 |
| CN115862740B (zh) * | 2022-12-06 | 2023-09-12 | 中国人民解放军军事科学院军事医学研究院 | 一种面向大规模病毒基因组数据的快速分布式多序列比对方法 |
| CN115862735A (zh) * | 2022-12-28 | 2023-03-28 | 郑州思昆生物工程有限公司 | 一种核酸序列检测方法、装置、计算机设备及存储介质 |
| CN115862735B (zh) * | 2022-12-28 | 2024-02-27 | 郑州思昆生物工程有限公司 | 一种核酸序列检测方法、装置、计算机设备及存储介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105243297A (zh) | 一种参考基因组上基因序列片段的快速比对定位方法 | |
| McKnight et al. | microDecon: A highly accurate read‐subtraction tool for the post‐sequencing removal of contamination in metabarcoding studies | |
| Alkhnbashi et al. | Characterizing leader sequences of CRISPR loci | |
| Auch et al. | Digital DNA-DNA hybridization for microbial species delineation by means of genome-to-genome sequence comparison | |
| CN104504304B (zh) | 一种成簇的规律间隔的短回文重复序列识别方法及装置 | |
| Ju et al. | 16S rRNA gene high-throughput sequencing data mining of microbial diversity and interactions | |
| Williams et al. | A robust species tree for the alphaproteobacteria | |
| Wolf et al. | Genome trees and the tree of life | |
| Liao et al. | A new unsupervised binning approach for metagenomic sequences based on n-grams and automatic feature weighting | |
| US11347810B2 (en) | Methods of automatically and self-consistently correcting genome databases | |
| Saheb Kashaf et al. | Recovering prokaryotic genomes from host-associated, short-read shotgun metagenomic sequencing data | |
| Bicego et al. | A bioinformatics approach to 2D shape classification | |
| Rasheed et al. | Metagenomic taxonomic classification using extreme learning machines | |
| CN102521528A (zh) | 一种基因序列数据的筛选方法 | |
| Goussarov et al. | Introduction to the principles and methods underlying the recovery of metagenome‐assembled genomes from metagenomic data | |
| Wu et al. | DeepRetention: a deep learning approach for intron retention detection | |
| Chen et al. | Tree2GD: a phylogenomic method to detect large-scale gene duplication events | |
| Beiko et al. | Detecting lateral genetic transfer: a phylogenetic approach | |
| Yuan et al. | RNA-CODE: a noncoding RNA classification tool for short reads in NGS data lacking reference genomes | |
| Ji et al. | HOTSPOT: hierarchical host prediction for assembled plasmid contigs with transformer | |
| Machado et al. | Evidence of absence treated as absence of evidence: The effects of variation in the number and distribution of gaps treated as missing data on the results of standard maximum likelihood analysis | |
| Solovyev et al. | Automatic annotation of bacterial community sequences and application to infections diagnostic | |
| Zinger et al. | Assessment of microbial communities by graph partitioning in a study of soil fungi in two alpine meadows | |
| CN102663287B (zh) | 利用编码转化实现基于序列比对的攻击特征提取的方法 | |
| Gudodagi et al. | Investigations and Compression of Genomic Data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WD01 | Invention patent application deemed withdrawn after publication | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20160113 |