CN116130001A - 一种基于k-mer定位的三代序列比对算法 - Google Patents
一种基于k-mer定位的三代序列比对算法 Download PDFInfo
- Publication number
- CN116130001A CN116130001A CN202211653043.0A CN202211653043A CN116130001A CN 116130001 A CN116130001 A CN 116130001A CN 202211653043 A CN202211653043 A CN 202211653043A CN 116130001 A CN116130001 A CN 116130001A
- Authority
- CN
- China
- Prior art keywords
- sequence
- mer
- genome
- comparison
- 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.)
- Pending
Links
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
Abstract
本发明提供了一种基于k‑mer定位的三代序列比对算法。首先利用哈希表构建基因组序列的k‑mer位置文库;然后根据待比对序列的每个k‑mer,通过哈希函数转换,可方便快速查找待比对序列的k‑mer在基因组中的所有位置;对待比对序列的每个k‑mer进行打分,用来衡量作为比对起始位置的可信度,得到每个k‑mer的得分值后,选取得分值最大的k‑mer即可快速找到每条待比对序列在基因组中的比对起始位置;比对起始位置可将序列和基因组划分为上游序列对和下游序列对,采用列降维带状打分分别对上游序列对和下游序列对进行比对,可避免传统带状比对的大规模矩阵存储问题,降低比对阶段的内存消耗;最后合并上有序列对和下游序列对的比对结果,得到最终的序列比对结果。
Description
技术领域
本发明涉及一种DNA基因序列处理方法,主要是一种基于k-mer定位的三代序列比对算法。
背景技术
序列比对是序列分析的重要研究内容,也是后续生物信息挖掘的基础。第三代单分子测序(single molecule sequencing,SMS)技术产生的序列长,但错误率较高(~15%)。现有的大多数序列比对算法都是针对第二代测序技术产生的数据(序列短,错误率低),并不适用于处理三代序列数据,因此需要开发新的三代序列比对算法。目前,针对三代序列的比对方法主要包括基于哈希搜索种子的比对方法、基于BWT-FM索引的种子搜索比对方法及基于已有比对工具的种子搜索比对方法。各类方法各具优势但也存在一些局限性,普遍存在的问题是这些方法得到的比对结果都是局部序列比对,导致每个方法的比对灵敏度和比对覆盖率较低,且对测序误差较敏感。
发明内容
为了克服现有方法的不足,本发明提供基于k-mer定位和列降维带状打分的三代序列比对方法(简称smsMap比对方法)。
本发明目的在于针对三代序列长度长、错误率高,现有比对方法序列比对灵敏度较低、对测序误差鲁棒性差以及比对覆盖率较低等问题,提出一种基于k-mer定位和列降维带状打分(SMS sequence mapping,smsMap)的三代序列比对方法,该方法有较高的比对灵敏度和比对覆盖率,且对测序误差鲁棒性强,为第三代单分子测序数据分析提供有效技术保障。
为实现上述目的,本发明技术方案的基本思想是:针对基因组序列,首先提取基因组序列的所有k-mer子片段并对每个k-mer进行哈希转换得到哈希值,根据哈希值存储在基因组中的位置,构建基因组k-mer位置文库;然后,对于待比对序列,提取k-mer并通过基因组k-mer位置文库找到待比对序列每个k-mer在基因组中的位置,对待比对序列每个k-mer进行打分,衡量作为比对起始位置的可信度,根据得分最高的k-mer得到在待比对序列和基因组中的比对起始位置;最后,采用列降维带状打分得到上游序列对和下游序列对的比对结果,合并后完成整条序列的比对结果。
本发明基于k-mer定位和列降维带状打分三代序列比对方法包括如下步骤:
步骤1:构建基因组k-mer位置文库
基因组k-mer位置文库指的是存储基因组k-mer子片段位置的哈希表,首先提取基因组序列的所有k-mer子片段,然后采用哈希函数进行k-mer转换,并将k-mer在基因组中的位置存储到哈希表中,具体实现过程为:
1)根据k-mer大小(Γ)创建长度为4Γ的哈希表,即数组,用来存储相应k-mer在基因组中的位置;
2)提取基因组序列的所有k-mer,k-mer是指基因组序列中包含k个碱基的子片段,对于一条长度为L的基因组DNA序列,在k-mer长度为Γ的情况下,基因组所有k-mer个数为L-Γ+1;
3)对基因组序列第一个k-mer进行哈希转换,假设该k-mer(长度为Γ)可表示为:w=c1,c2,...,cΓ,其在哈希表中的存储地址可通过以下哈希函数计算得到:
式中4Γ-γ是k-mer中第γ个位置上碱基(cγ)的权重,I(cγ)是索引函数,定义为:
每个k-mer的哈希编码可看作是Γ位四进制数的一个转换,通过公式(1)计算基因组序列每个k-mer哈希值即为该k-mer在位置文库中的索引位置,然后将此k-mer在基因组中的位置存储到该索引下的数组中;
4)重复步骤3)计算基因组所有k-mer的哈希值并存储其在基因组中的位置,当存储完所有k-mer位置信息后,即为构建的基因组k-mer位置文库;
步骤2:定位比对起始位置
在定位比对起始位置阶段,首先提取待比对序列的所有k-mer;然后根据上一步构建的基因组k-mer位置文库找到待比对序列每个k-mer在基因组中出现的位置;进而计算待比对序列每个k-mer作为比对起始位置的可信度得分值,最后找出得分值最大的k-mer,即可得到在待比对序列和基因组中的位置信息,作为比对的起始位置;具体实现步骤如下:
1)对于待比对序列r,提取r的所有k-mer,然后通过公式(1)计算哈希值,即可找到该k-mer在基因组中的位置信息,序列r每个k-mer及其在基因组中的位置可用一个三元组表示:
式中i表示序列r中第i个k-mer,是序列r第i个k-mer在基因组中第l个匹配的位置,Li是序列r第i个k-mer在基因组中匹配的总个数,是的修正位置,即序列r第i个k-mer在基因组中的位置减去在序列r中的位置;
式中Θ是序列r匹配到基因组中的k-mer个数,Lj是序列r第j个k-mer在基因组中匹配的总个数,函数δ为指示函数,L(r)是序列r的容错长度,定义为:
L(r)=0.2len(r) (6)
式中len(r)是序列r的长度。
3)根据公式(4)计算序列r中每个(i=1,2,L,Θ,l=1,2,L,Li)的可信度得分值,然后选取得分值最大的作为比对的起始k-mer,再根据公式(3)找出该k-mer在序列r和基因组中的相应位置,作为下一步的比对起始位置;
步骤3:列降维带状打分比对
根据步骤2找出的比对起始位置一般位于序列r和基因组中间,可将基因组和待比对序列划分为上游序列对(ru和gu)和下游序列对(rd和gd),接下来采用列降维带状打分比对法对上游序列对和下游序列对进行详细比对,然后合并形成最终的序列比对结果;具体实现步骤如下:
1)根据比对起始位置将对序列和基因组划分为上游序列对(ru和gu)和下游序列对(rd和gd);
2)首先创建大小为l(rd)×2b的列降维打分矩阵M,l(rd)为序列rd的长度,b为带状比对的带宽:b=0.1×l(rd);
3)采用动态规划算法对矩阵M进行打分,打分公式为:
v′=v+sci(u)-sci(u-1) (9)
v″=v+sci(v) (10)
sci(u)=maxfloor(ldown(u)),0 (11)
ldown(u)=1.2×u-b (12)
其中F(u,v)是矩阵M第u行第v列元素的得分值,Score[rd(u),gd(v″)]是序列rd中第u个碱基和基因组gd中第v″个碱基的匹配得分函数,sci(u)是列降维矩阵第u行起始列的索引值,floor()是向下取整函数;然后根据回溯路径即可得到序列rd和序列gd的比对结果;
4)再采用步骤2)和3)对上游序列对(ru与gu)进行打分得到比对结果,然后将上游比对结果和下游比对结果进行合并,得到序列r与基因组最终完整的序列比对结果。
本发明具有如下有益效果:
1、对每条待比对序列的每个k-mer进行打分,根据得分值最高的k-mer即可得到待比对序列和基因组的比对起始位置,可将每条序列定位到基因组中,增强方法定位序列的灵敏度,提高测序误差的鲁棒性。
2、采用列降维带状打分方式实现上游序列对和下游序列对的得分矩阵,避免只对序列局部区域比对,提高比对覆盖率,同时对打分矩阵进行列降维,进一步降低比对阶段的内存消耗。
附图说明
图1是smsMap比对方法流程图。其中,图(a)是基因组k-mer位置文库构建,图(b)是定位比对起始位置,图(c)是列降维带状比对打分矩阵。
图2是基因组序列k-mer位置文库构建流程图,其中,图(a)是提取基因组的所有k-mer子片段,图(b)是k-mer哈希值计算,图(c)是构建的基因组k-mer位置文库。
图4是比对起始位置将序列和基因组划分为上游序列对(ru,gu)和下游序列对(rd,gd)。
图5是列降维带状比对转换示意图。其中,图(a)是传统带状矩阵得分值,图(b)是提取图(a)中计算得分的区域,图(c)是列降维带状比对得分矩阵。
图6是根据图5得到的比对结果,其中“|”表示匹配状态,“-”表示插入或删除状态。
图7是不同测序误差下的比对结果。其中,图(a)是比对序列率(FAR)变化曲线,图(b)是比对碱基率(FAB)变化曲线,图(c)是平均比对覆盖率(ACR)变化曲线。
具体实施方式
在Intel Xeon E5-2667V4@3.2GHz、128GB运行内存服务器上,基于Ubuntu16.04.5版本的Linux平台,选取六组不同误差率的模拟数据和四组真实PacBio测序数据集进行比对仿真实验,模拟数据集由模拟软件产生,真实数据集由三代测序平台PacBio产生。
设置k-mer长度为11,用E.coli序列数据集第一条序列r(长度为2824bp)阐述具体比对步骤。
步骤1:构建基因组k-mer位置文库
1)由于k-mer长度设定为11,创建长度为411(4194304)的数组;
2)提取基因组g所有k-mer,由于g长度为4681865bp,k-mer长度为11,所以总的k-mer个数为4681855。
3)计算基因组g第一个k-mer哈希值,第一个k-mer序列为CAAGCCAGCCA,通过公式(1)可以得到哈希值为1086612。由于此k-mer是基因组第一个k-mer,将0存储到第1086612个位置的哈希表中。
4)重复步骤3),直到计算基因组所有k-mer哈希值并将其在基因组中的位置存储到哈希表中,即为构建的基因组k-mer位置文库。
步骤2:定位比对起始位置
1)提取序列r第一个k-mer,w1=CTTGTGGTGAT,根据公式(1)得到其哈希值H(w1)=2079459,通过查找基因组k-mer位置文库第2079459个存储的元素即可得知在基因组g中出现1次,在第445058位置,根据公式(3)即可构建w1的三元组信息:
3)选取得分值最大的k-mer,得到最大的可信度得分为486,根据三元组信息得知是序列r中第548个k-mer,匹配到基因组中第785562位置。
步骤3:列降维带状打分比对
1)根据可比对起始将序列分为上游序列对(ru,gu)和下游序列对(rd,gd),其中ru长度为l(ru)=547,rd长度为l(rd)=2277。
2)对下游序列对,创建大小为2277×558的打分矩阵,采用公式(7)得到矩阵的打分值,并根据回溯路径得到比对结果。
3)对上游序列对,创建大小为547×108的打分矩阵,采用公式(7)得到矩阵的打分值,并根据回溯路径得到比对结果。
4)合并上游比对结果和下游比对结果,最终得到序列r的完整比对结果。
图5为列降维打分比对转换示意图,相对于图(a)传统带状比对,图(c)只存储打分部分,避免了构建原始大小矩阵,可以得到和图(a)中一样的打分结果,降低了内存消耗。
图6是根据图5得出的详细比对结果,可以得到序列的完整比对结果,即将序列每个碱基都比对到基因组上,提高比对灵敏度。
表1是比对阶段和传统带状比对方法比较的内存消耗,可以看出通过列降维带状比对,smsMap方法避免了大规模矩阵的存储,可将内存降低到原来的五分之一左右。
根据公式(13)计算四组真实序列数据的比对序列率(FAR)、比对碱基率(FAB)和比对覆盖率(ACR),如表2所示。
表1比对阶段内存消耗大小(GB)
表2 smsMap序列比对结果
FAR和FAB可以反映比对灵敏度,即反映每个方法比对的序列条数和碱基数,ACR可反应比对完整性。由表2可知,smsMap方法有较高的FAR、FAB和ACR值,表明smsMap有较高的比对灵敏度和比对完整性。
图7是对不同测序误差下的比对结果,可以看出随着测序误差率增大,smsMap方法均表现较高的FAR、FAB和ACR值,表明smsMap方法对测序误差具有较强的鲁棒性。
以上结果表明,smsMap三代序列比对方法可以对序列长、测序误差率高的三代测序数据进行序列比对,可以得到更多的比对序列和比对碱基,同时比对完整性好,且在序列比对阶段可有效降低内存消耗。适合各种基因组三代测序数据的比对,潜在应用价值大。
Claims (3)
1.一种基于k-mer定位的三代序列比对算法,其特征在于,包括下述步骤:
步骤1:构建基因组k-mer位置文库
基因组k-mer位置文库指的是存储基因组k-mer子片段位置的哈希表,首先提取基因组序列的所有k-mer子片段,然后采用哈希函数进行k-mer转换,并将k-mer在基因组中的位置存储到哈希表中,具体实现过程为:
1)根据k-mer大小(Γ)创建长度为4Γ的哈希表,即数组,用来存储相应k-mer在基因组中的位置;
2)提取基因组序列的所有k-mer,k-mer是指基因组序列中包含k个碱基的子片段,对于一条长度为L的基因组DNA序列,在k-mer长度为Γ的情况下,基因组所有k-mer个数为L-Γ+1;
3)对基因组序列第一个k-mer进行哈希转换,假设该k-mer(长度为Γ)可表示为:w=c1,c2,...,cΓ,其在哈希表中的存储地址可通过以下哈希函数计算得到:
式中4Γ-γ是k-mer中第γ个位置上碱基(cγ)的权重,I(cγ)是索引函数,定义为:
每个k-mer的哈希编码可看作是Γ位四进制数的一个转换,通过公式(1)计算基因组序列每个k-mer哈希值即为该k-mer在位置文库中的索引位置,然后将此k-mer在基因组中的位置存储到该索引下的数组中;
4)重复步骤3)计算基因组所有k-mer的哈希值并存储其在基因组中的位置,当存储完所有k-mer位置信息后,即为构建的基因组k-mer位置文库;
步骤2:定位比对起始位置
在定位比对起始位置阶段,首先提取待比对序列的所有k-mer;然后根据上一步构建的基因组k-mer位置文库找到待比对序列每个k-mer在基因组中出现的位置;进而计算待比对序列每个k-mer作为比对起始位置的可信度得分值,最后找出得分值最大的k-mer,即可得到在待比对序列和基因组中的位置信息,作为比对的起始位置;具体实现步骤如下:
1)对于待比对序列r,提取r的所有k-mer,然后通过公式(1)计算哈希值,即可找到该k-mer在基因组中的位置信息,序列r每个k-mer及其在基因组中的位置可用一个三元组表示:
式中i表示序列r中第i个k-mer,是序列r第i个k-mer在基因组中第l个匹配的位置,Li是序列r第i个k-mer在基因组中匹配的总个数,是的修正位置,即序列r第i个k-mer在基因组中的位置减去在序列r中的位置;
式中Θ是序列r匹配到基因组中的k-mer个数,Lj是序列r第j个k-mer在基因组中匹配的总个数,函数δ为指示函数,L(r)是序列r的容错长度,定义为:
L(r)=0.2len(r) (6)
式中len(r)是序列r的长度;
步骤3:列降维带状打分比对
根据步骤2找出的比对起始位置一般位于序列r和基因组中间,可将基因组和待比对序列划分为上游序列对(ru和gu)和下游序列对(rd和gd),接下来采用列降维带状打分比对法对上游序列对和下游序列对进行详细比对,然后合并形成最终的序列比对结果;具体实现步骤如下:
1)根据比对起始位置将对序列和基因组划分为上游序列对(ru和gu)和下游序列对(rd和gd);
2)首先创建大小为l(rd)×2b的列降维打分矩阵M,l(rd)为序列rd的长度,b为带状比对的带宽:b=0.1×l(rd);
3)采用动态规划算法对矩阵M进行打分,打分公式为:
v′=v+sci(u)-sci(u-1) (9)
v″=v+sci(v) (10)
sci(u)=max[floor(ldown(u)),0] (11)
ldown(u)=1.2×u-b (12)
其中F(u,v)是矩阵M第u行第v列元素的得分值,Score[rd(u),gd(v″)]是序列rd中第u个碱基和基因组gd中第v″个碱基的匹配得分函数,sci(u)是列降维矩阵第u行起始列的索引值,floor()是向下取整函数;然后根据回溯路径即可得到序列rd和序列gd的比对结果;
4)再采用步骤2)和3)对上游序列对(ru与gu)进行打分得到比对结果,然后将上游比对结果和下游比对结果进行合并,得到序列r与基因组最终完整的序列比对结果。
2.根据权利1所述的基于k-mer定位的三代序列比对算法,其特征在于:所述步骤2通过哈希函数转换,可方便快速查找待比对序列每个k-mer在基因组中的所有位置,然后对待比对序列的每个k-mer作为比对起始位置的可信度进行打分,得到每个k-mer的得分值,选取得分值最大的k-mer即可快速找到每条待比对序列在基因组中的比对起始位置,增强方法的比对灵敏度。
3.根据权利1所述的基于k-mer定位的三代序列比对算法,其特征在于:所述步骤3根据列降维带状打分函数分别对上游序列对和下游序列对进行比对,然后合并得到最终的序列比对结果,保证序列的每个碱基都能比对到基因组上,可以得到更多的碱基比对结果,同时降低了传统带状比对阶段的内存消耗。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211653043.0A CN116130001A (zh) | 2022-12-21 | 2022-12-21 | 一种基于k-mer定位的三代序列比对算法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211653043.0A CN116130001A (zh) | 2022-12-21 | 2022-12-21 | 一种基于k-mer定位的三代序列比对算法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN116130001A true CN116130001A (zh) | 2023-05-16 |
Family
ID=86298474
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202211653043.0A Pending CN116130001A (zh) | 2022-12-21 | 2022-12-21 | 一种基于k-mer定位的三代序列比对算法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116130001A (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116665772A (zh) * | 2023-05-30 | 2023-08-29 | 之江实验室 | 一种基于内存计算的基因组图分析方法、装置和介质 |
-
2022
- 2022-12-21 CN CN202211653043.0A patent/CN116130001A/zh active Pending
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116665772A (zh) * | 2023-05-30 | 2023-08-29 | 之江实验室 | 一种基于内存计算的基因组图分析方法、装置和介质 |
CN116665772B (zh) * | 2023-05-30 | 2024-02-13 | 之江实验室 | 一种基于内存计算的基因组图分析方法、装置和介质 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US20210193257A1 (en) | Phase-aware determination of identity-by-descent dna segments | |
CN107133493B (zh) | 基因组序列的组装方法、结构变异探测方法和相应的系统 | |
CA2424031C (en) | System and process for validating, aligning and reordering genetic sequence maps using ordered restriction map | |
WO2018218788A1 (zh) | 一种基于全局种子打分优选的三代测序序列比对方法 | |
CN110299185B (zh) | 一种基于新一代测序数据的插入变异检测方法及系统 | |
EP3084426B1 (en) | Iterative clustering of sequence reads for error correction | |
CN110621785B (zh) | 基于三代捕获测序对二倍体基因组单倍体分型的方法和装置 | |
CN116130001A (zh) | 一种基于k-mer定位的三代序列比对算法 | |
WO2009155443A2 (en) | Method and apparatus for sequencing data samples | |
CN110692101A (zh) | 用于比对靶向的核酸测序数据的方法 | |
JP2006075162A (ja) | 遺伝子の転写物マッピング方法及びシステム | |
CN115083521A (zh) | 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及系统 | |
CN115631789A (zh) | 一种基于泛基因组的群体联合变异检测方法 | |
US20140121983A1 (en) | System and method for aligning genome sequence | |
US9323889B2 (en) | System and method for processing reference sequence for analyzing genome sequence | |
US20160098517A1 (en) | Apparatus and method for detecting internal tandem duplication | |
US20150142328A1 (en) | Calculation method for interchromosomal translocation position | |
CN114564306B (zh) | 一种基于GPU并行计算的第三代测序RNA-seq比对方法 | |
US9348968B2 (en) | System and method for processing genome sequence in consideration of seed length | |
US20170132361A1 (en) | Sequence assembly method | |
WO2019023978A1 (zh) | 比对方法、装置及系统 | |
CN117292751A (zh) | 一种基于最长路径搜索的三代序列比对方法 | |
KR101584857B1 (ko) | 염기 서열 정렬 시스템 및 방법 | |
Tang et al. | Integration of hybrid and self-correction method improves the quality of long-read sequencing data | |
CN110875084B (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 |