CN112687334B - 一种可应用于传染病病原体测序的读段映射延伸方法 - Google Patents
一种可应用于传染病病原体测序的读段映射延伸方法 Download PDFInfo
- Publication number
- CN112687334B CN112687334B CN202011597128.2A CN202011597128A CN112687334B CN 112687334 B CN112687334 B CN 112687334B CN 202011597128 A CN202011597128 A CN 202011597128A CN 112687334 B CN112687334 B CN 112687334B
- Authority
- CN
- China
- Prior art keywords
- moving
- scanning
- fragment
- sequence
- scan
- 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
- 238000000034 method Methods 0.000 title claims abstract description 43
- 238000013507 mapping Methods 0.000 title claims abstract description 23
- 238000012163 sequencing technique Methods 0.000 title claims abstract description 17
- 244000052769 pathogen Species 0.000 title claims abstract description 5
- 208000035473 Communicable disease Diseases 0.000 title claims abstract description 4
- 208000015181 infectious disease Diseases 0.000 title claims abstract description 4
- 230000001717 pathogenic effect Effects 0.000 title claims abstract description 4
- 230000001360 synchronised effect Effects 0.000 claims abstract description 27
- 239000012634 fragment Substances 0.000 description 37
- 238000004364 calculation method Methods 0.000 description 8
- 238000005516 engineering process Methods 0.000 description 4
- 230000004807 localization Effects 0.000 description 4
- 239000012297 crystallization seed Substances 0.000 description 3
- 238000001712 DNA sequencing Methods 0.000 description 2
- 108091028043 Nucleic acid sequence Proteins 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 238000012165 high-throughput sequencing Methods 0.000 description 2
- 238000007481 next generation sequencing Methods 0.000 description 2
- 108090000623 proteins and genes Proteins 0.000 description 2
- 208000025721 COVID-19 Diseases 0.000 description 1
- 241000700605 Viruses Species 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000012217 deletion Methods 0.000 description 1
- 230000037430 deletion Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000004907 flux Effects 0.000 description 1
- 238000003780 insertion Methods 0.000 description 1
- 230000037431 insertion Effects 0.000 description 1
- 230000003447 ipsilateral effect Effects 0.000 description 1
- 238000012177 large-scale sequencing Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 108020004707 nucleic acids Proteins 0.000 description 1
- 150000007523 nucleic acids Chemical class 0.000 description 1
- 102000039446 nucleic acids Human genes 0.000 description 1
- 230000002265 prevention Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000007480 sanger sequencing Methods 0.000 description 1
- 238000010561 standard procedure Methods 0.000 description 1
Images
Landscapes
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明涉及一种可应用于传染病病原体测序的读段映射延伸方法。本发明提出的方法通过将同步扫描与广度优先搜索结合，大大减少了计算量，提高了映射延伸的速度。
Description
技术领域
本发明涉及一种DNA测序过程中读段映射延伸方法。
背景技术
DNA测序是分子生物学等学科研究的基础工作。例如研究病原体的DNA序列。在目前全球的COVID-19疫情中，分析病毒DNA序列是防疫工作的基础。荧光标记的Sanger法是第一代测序技术的标准方法。但是该方法通量低，不适合用于大规模测序工作。高通量测序(High-Throughput Sequencing,HTS)是对传统Sanger测序的革命性变革，其解决了一代测序一次只能测定一条序列的限制，一次运行即可同时得到几十万到几百万条核酸分子的序列，因此也被称为新一代测序(Next Generation Sequencing,NGS)或第二代测序。第二代测序技术虽然测序的通量大大增加，但是其获得单条序列长度很短，想要得到准确的基因序列信息依赖于较高的测序覆盖度和准确的序列拼接技术。其中对短序列的读取和识别是非常大的工作量。
在利用第二代测序方法进行基因测序时，首先会得到许多的读段，随后，需要把读段映射到参考基因组上。而映射过程分为两步：第一步是种子序列定位，第二步是延伸扩展。
目前常用的映射延伸方法是动态规划方法。该方法的缺陷是当数据量增大时，总的计算量及存贮量急剧增大。
例如：当读段为：ACCAGTCAACTGTGCATGTCGCATGTATGCATCAATGCG
参考基因组为：
GCTAGAACGTATACCAGTAACTGTGCATGTCGCATGTATGCATGTAATGCGGCTACATGCCGATGCG时，利用现有方法进行读段比对时，需要进行如下操作。
对读段中种子序列两侧中的一侧的选取片段(以下称片段A)，从参考基因组中的种子序列的同侧也取一个片段(以下称片段B)，长度为片段A的两倍。然后使用动态规划将片段A映射到片段B前部的某个片段(以下称片段C)上。
种子序列定位后的结果为如图1所示。其中种子序列以阴影标出。若现在要对种子序列右侧的片段进行延伸，则片段A为斜体片段，长度为11bp，片段B为加粗片段，长度为22bp。
随后执行以下步骤(以表格对方方式表示)：
1.准备一张12行23列的表格(不包括行标题和列标题)，第1列的列标题为空，其他各列的列标题依次为片段B中的一个bp；第1行的行标题为空，其他各行的行标题依次为片段A中的一个bp。准备好的表格如下：
表1 12行23列的初表格
2.将表的第1行从左至右依次填入0到22，将表的第1列从上至下依次填入0到11。填入后的表格如下：
表2填数字表格
3.从第2行第2列开始，依次填充表中的每一个单元格。填充的顺序是：从左向右填充，当右侧没有单元格时，从下一行的第2列开始，继续从左向右填充。每个单元格填充的内容为以下三个值中的最小值：
1)该单元格上面的那个单元格的值加1；
2)该单元格左边的那个单元格的值加1；
3)该单元格左上的那个单元格的值加V。其中，V的值为：如果该单元格所在的行标题和列标题相同，则V的值为0；否则V的值为1。
填充后的表格如下：
表3填充后的表格
4.看最后一行中的最小值在第几列(称为第i列)，则片段B的前i-1个bp即为所求的片段C。本例中，最后一行的最小值为2，在第13列，所以片段B的前12个bp(GCATGTAATGCG)即为所求的片段C。
可以看到，传统方法需要填充一个12*23大小的表格，共276个单元格，共需计算276次。可见传统方法的计算量较大，速度较慢。
本申请中涉及到的术语定义如下。
读段：一个字符串。例如：ACCAGTCAACTGTGCA。
参考基因组：一个字符串。例如：GTAACTGTGCATGTCGCATGTATGCATGTAATGC
bp：字符串的长度单位，一个字母就是一个bp，字符串中的第i个字符就是第ibp。例如：AAACTTGGA，长度为9bp，其中字母“C”是第4bp。
编辑距离：两个字符串A和B，如果要把A变成B，最少需要i次编辑操作，则称A与B的编辑距离是i。其中，一次编辑操作是指下列三种操作之一：1.修改一个字母，2.插入一个字母，3.删除一个字母。例如：两个字符串分别为：
字符串A：ACTCTAGTATGTGCATGCGCGCCATGTGTGCATGGGCAT
字符串B：ACTCGTAGTATGAGCATGTGCGCCATGTGTGCTGGGCAT
至少需要4次编辑操作，才可将以将字符串A变成字符串B。具体的编辑操作如下：
ACTC-TAGTATGTGCATGCGCGCCATGTGTGCATGGGCAT
ACTCGTAGTATGAGCATGTGCGCCATGTGTGC-TGGGCAT
1)插入字母“G”(以红色标出)，2)修改字母“T”为“A”(以蓝色标出)，3)修改字母“C”为“T”(以黄色标出)，4)删除字母“A”(以绿色标出)。
所以这两个字符串的编辑距离为4。
读段映射：将读段映射到参考基因组的某个片段上。映射的结果就是要在参考基因组上找一个片段，使得该片段与读段之间的编辑距离尽可能地小。例如：
当读段为：
ACCAGTCAACTGTGCATGTCGCATGTATGCATGAATGCG
参考基因组为：
GCTAGAACGTATACCAGTAACTGTGCATGTCGCATGTATGCATGTAATGCGGCTACATGCCGATGCG时，映射后的结果如图2所示。
在该映射结果中，在参考基因组上找到的片段是第13bp至第51bp之间(含第13bp和第51bp)的片段(以灰底标出)，读段与该片段的编辑距离为3，仅有1处修改(以红字标出),1处插入(以绿字标出)和1处删除(以蓝字标出)。
参考基因组上除此片段之外的其他任何片段，与读段的编辑距离均大于3。
种子序列定位：在读段和参考基因组中都找一个固定长度(称为种子长度，种子长度预先设定)的片段(称为种子序列)，使得二者完全相同。例如，在示例一中，如果预先设定的种子长度为20bp，则可以在基因组上找到第20bp至第39bp之间(含第20bp和第39bp)的片段，使其与读段上第9bp至第28bp之间(含第9bp和第28bp)的片段完全相同(两者均以黄字标出)。
映射延伸：在完成种子序列定位以后，继续把读段中种子序列两侧的片段映射到参考基因组同侧的片段上。
发明内容
本发明的目的是优化现有测序过程中的读段映射延伸方法，以期降低计算量，提高测序的速度。
本发明提供的测序读段映射延伸方法包括如下步骤：
(1)对种子序列的两侧分别进行延伸，延伸方法相同；
(2)初始化扫描起始位置。在参考基因组和读段上各确定一个扫描起始位置。在参考基因组上的扫描起始位置称为位置A，在读段上的扫描起始位置称为位置B。两者的初始位置都是延伸侧与种子序列相邻的第一个bp；
(3)进行同步扫描。向延伸方向同时移动位置A和位置B，直到有一个到达序列末尾，或者两者指向的字母不相同为止；
(4)当位置A和位置B指向的字母不同时，进行广度优先搜索；
i.如果有某种位置移动方式可以继续进行同步扫描，则搜索结束；或
ii.如果都不能继续进行同步扫描，则对每一种位置移动方式，在其移动的结果上继续进行上述三种位置移动，直到可以继续进行同步扫描；或者
iii.某个位置到达片段尾为止。
根据本发明的实施例，在种子序列定位后，对右侧进行了映射延伸，共进行了两次同步扫描和一次广度优先搜索，两次同步扫描分别需要计算5次和7次，广度优先搜索需要计算8次，共需计算20次。与传统方法需要计算276次相比，本发明提出的方法通过将同步扫描与广度优先搜索结合，大大减少了计算量，提高了映射延伸的速度。
附图说明
图1为种子序列定位后的结果，其中种子序列以阴影标出。若现在要对种子序列右侧的片段进行延伸，则片段A为斜体片段，长度为11bp，片段B为加粗片段，长度为22bp；
图2为读段映射后的结果；
图3为本发明提供的测序读段映射延伸方法中单侧读段映射延伸方法流程图；
图4为实施例1初始化扫描起始位置种子序列定位后的结果；其中，种子序列以黄字标出。若现在要对种子序列右侧的片段进行延伸，则位置A为红色箭头且圈起来所指的字母“G”，位置B为绿色箭头所指的字母“G”。
图5为进行同步扫描位置A和位置B移动后将到达的位置示意图；
图6为进行广度优先搜索可能的第一种位置移动方式；
图7为进行广度优先搜索可能的第二种位置移动方式；
图8为进行广度优先搜索可能的第三种位置移动方式；
图9为搜索过程的结果1；
图10为搜索过程的结果2；
图11为搜索过程的结果3；
图12为搜索过程的结果4；
图13为搜索过程的结果5；
图14为循环执行同步扫描和广度优先搜索直到有一个位置到达序列末尾的结果；
图15为在实施例1中位置A的初始位置和最终位置示意图。
具体实施例
实施例1
本实施例以读段为：ACCAGTCAACTGTGCATGTCGCATGTATGCATGAATGCG
参考基因组为：
GCTAGAACGTATACCAGTAACTGTGCATGTCGCATGTATGCATGTAATGCGGCTACATGCCGATGCG
为例进行说明。
按现有技术方法确定种子序列片段后，对种子序列的两侧分别进行读段延伸，读段方法相同，并按照图3的方法进行。
1.初始化扫描起始位置。在参考基因组和读段上各确定一个扫描起始位置。在参考基因组上的扫描起始位置称为位置A，在读段上的扫描起始位置称为位置B。两者的初始位置都是延伸侧与种子序列相邻的第一个bp。因此，读段和参考基因组，种子序列定位后的结果为图4所示。其中，种子序列以黄字标出。若现在要对种子序列右侧的片段进行延伸，则位置A为红色箭头且圈起来所指的字母“G”，位置B为绿色箭头所指的字母“G”。
2.进行同步扫描。向延伸方向同时移动位置A和位置B，直到有一个到达序列末尾，或者两者指向的字母不相同为止。位置A和位置B移动后将到达的位置如图5所示。其中位置A达到红色箭头且圈起来所指的字母“G”，位置B达到绿色箭头所指的字母“C”。
在同步扫描过程中，每次要对位置A和位置B指向的字母进行一次比较，所需的计算量为扫描的长度。在本例中，需要计算5次。
3.进行广度优先搜索。当位置A和位置B指向的字母不同时，可能的位置移动方式有三种：
第一种：位置A不动，把位置B向后移动一个bp。例如，图二中的位置按这种方式移动后如图6所示。其中位置A仍留在红色箭头且圈起来所指的字母“G”，位置B移动到绿色箭头所指的字母“A”。
第二种：位置B不动，把位置A向后移动一个bp。例如，图二中的位置按这种方式移动后如图7所示。其中位置A移动到红色箭头且圈起来所指的字母“T”，位置B仍留在绿色箭头所指的字母“C”。
第三种：把位置A和位置B同时向后移动一个bp。例如，图二中的位置按这种方式移动后如图8所示。其中位置A移动到红色箭头且圈起来所指的字母“T”，位置B移动到绿色箭头所指的字母“A”。
在搜索过程中，尝试以上三种位置移动方式，看哪一种方式能够继续进行同步扫描。如果有某种位置移动方式可以继续进行同步扫描，则搜索结束。如果都不能继续进行同步扫描，则对每一种位置移动方式，在其移动的结果上继续进行上述三种位置移动，直到可以继续进行同步扫描，或者某个位置到达片段尾为止。例如：图二中的位置，无论按上述三种方式的哪一种移动，位置A和位置B指向的字母都不同，均不能继续进行同步扫描，所以在每种移动方式的基础上再进行上述三种方式的移动，此时会产生九种结果。实际上，这九种结果中，有一些是重复的，例如在第一种方式的基础上再按第二种方式移动，和在第二种方式的基础上再按第一种方式移动，其最终结果是一样的。所以实际上只有如下五种结果，分别如图9-13所示。
在这五个结果中再搜索，看是否可以继续进行同步扫描。可以发现，在结果四中，位置A和位置B指向的字母及其后面连续多个字母都相同，都是AATGCG，可以继续进行同步扫描，搜索结束，最终的位置A和位置B为结果4所示的位置。
在搜索过程中，每次尝试移动位置，都需要判断一次位置A和位置B指向的字母是否相同，在本例中，首先尝试了3次，接着又尝试了5次，共需要计算8次。
4.循环执行同步扫描和广度优先搜索，直到有一个位置到达序列末尾。例如，在上述结果四上继续进行同步扫描，直到位置B到达序列末尾，如图14所示。
5.结束。位置A的初始位置和最终位置之间的片段(含初始位置，不含最终位置)即为读段中种子序列一侧的片段映射到参考基因组上的片段。在本例中，位置A的初始位置和最终位置分别如图15所示。其中位置A的初始位置为绿色箭头所指的字母“G”，最终位置为红色箭头且圈起来所指的字母“G”。所以读段中种子序列右侧的片段映射到参考基因组上的片段为CATGTAATGCG。
可以看到，在本例中，对示例二中的读段和参考基因组，在种子序列定位后，对右侧进行了映射延伸，共进行了两次同步扫描和一次广度优先搜索，两次同步扫描分别需要计算5次和7次，广度优先搜索需要计算8次，共需计算20次。与传统方法需要计算276次相比，本发明提出的方法大大减少了计算量，提高了映射延伸的速度。
Claims (1)
1.一种可应用于传染病病原体测序的读段映射延伸方法，其特征在于
对种子序列的两侧分别进行延伸，延伸方法相同，所述的延伸方法如下；
（1）初始化扫描起始位置；在参考基因组和读段上各确定一个扫描起始位置；在参考基因组上的扫描起始位置称为位置A，在读段上的扫描起始位置称为位置B；两者的初始位置都是延伸侧与种子序列相邻的第一个bp；
（2）进行同步扫描；向延伸方向同时移动位置A和位置B，直到有一个到达序列末尾，或者两者指向的字母不相同为止；
（3）进行广度优先搜索；当位置 A 和位置 B 指向的字母不同时，可能的位置移动方式
有三种：
第一种：位置 A 不动，把位置 B 向后移动一个 bp；
第二种：位置 B 不动，把位置 A 向后移动一个 bp；
第三种：把位置 A 和位置 B 同时向后移动一个 bp；
在搜索过程中，尝试以上三种位置移动方式，看哪一种方式能够继续进行同步扫描；
如果有某种位置移动方式可以继续进行同步扫描，则搜索结束；如果都不能继续进行同
步扫描，则对每一种位置移动方式，在其移动的结果上继续进行上述三种位置移动，直
到可以继续进行同步扫描，或者某个位置到达片段尾为止。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011597128.2A CN112687334B (zh) | 2020-12-29 | 2020-12-29 | 一种可应用于传染病病原体测序的读段映射延伸方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011597128.2A CN112687334B (zh) | 2020-12-29 | 2020-12-29 | 一种可应用于传染病病原体测序的读段映射延伸方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112687334A CN112687334A (zh) | 2021-04-20 |
| CN112687334B true CN112687334B (zh) | 2022-09-23 |
Family
ID=75454092
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011597128.2A Active CN112687334B (zh) | 2020-12-29 | 2020-12-29 | 一种可应用于传染病病原体测序的读段映射延伸方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112687334B (zh) |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109234267B (zh) * | 2018-09-12 | 2021-07-30 | 中国科学院遗传与发育生物学研究所 | 一种基因组组装方法 |
| CN109097458A (zh) * | 2018-09-12 | 2018-12-28 | 山东省农作物种质资源中心 | 基于ngs读段搜索实现序列延伸的虚拟pcr方法 |
| US20200234795A1 (en) * | 2019-01-22 | 2020-07-23 | The Regents Of The University Of Michigan | Pruning Pair-HMM Algorithm And Hardware Architecture |
| CN111583998B (zh) * | 2020-05-06 | 2023-05-02 | 西安交通大学 | 一种考虑拷贝数变异因素的基因组结构变异分型方法 |
-
2020
- 2020-12-29 CN CN202011597128.2A patent/CN112687334B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN112687334A (zh) | 2021-04-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Hochsmann et al. | Local similarity in RNA secondary structures | |
| Benson et al. | Reconstructing the duplication history of a tandem repeat. | |
| KR20230008877A (ko) | Dna 기반 데이터 스토리지의 프로그램 및 기능 | |
| CN112735528A (zh) | 一种基因序列比对方法及系统 | |
| CN115014362B (zh) | 一种基于合成单元的牛耕式全覆盖路径规划方法和装置 | |
| CN105760706A (zh) | 一种二代测序数据的压缩方法 | |
| CN112687334B (zh) | 一种可应用于传染病病原体测序的读段映射延伸方法 | |
| CN104778687A (zh) | 一种图像匹配方法和装置 | |
| US20140121983A1 (en) | System and method for aligning genome sequence | |
| CN113724783B (zh) | 一种短串联重复序列重复数的检测和分型方法 | |
| CN102841988B (zh) | 一种对核酸序列信息进行匹配的系统和方法 | |
| WO2020052101A1 (zh) | 基于ngs读段搜索实现序列延伸的虚拟pcr方法 | |
| CN112712850A (zh) | 一种可应用于传染病病原体测序读段映射的种子序列定位方法 | |
| CN115148290A (zh) | 一种基于三代测序数据的补洞方法 | |
| CN114758721A (zh) | 一种基于深度学习的转录因子结合位点定位方法 | |
| JP5516880B2 (ja) | 配列解析装置、配列解析方法およびコンピュータプログラム | |
| CN112802553A (zh) | 一种基于后缀树算法的基因组测序序列与参考基因组比对的方法 | |
| CN111597856A (zh) | 一种基于光敏变色材料的混凝土标识提取方法 | |
| KR101953663B1 (ko) | 하나의 올리고뉴클레오티드를 이용해서 올리고뉴클레오티드 풀을 생산하는 방법 | |
| KR102380935B1 (ko) | 유전체 영역 검색 시스템 및 방법 | |
| Schonfeld et al. | Filtration and depth annotation improve non-linear projection for rna motif discovery | |
| Bucur | A stochastic de novo assembly algorithm for viral-sized genomes obtains correct genomes and builds consensus | |
| CN104239749A (zh) | 碱基序列对准系统及方法 | |
| Quan et al. | A Bidirectional Fuzzy Index and Approximate Search Algorithm for Next Generation Sequencing | |
| CN106681702A (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 | ||
| CB03 | Change of inventor or designer information |
Inventor after: Fan Liangliang Inventor after: Luo Mufeng Inventor after: Huang Hao Inventor after: Xiang Rong Inventor before: Xiang Rong Inventor before: Luo Mufeng Inventor before: Huang Hao Inventor before: Fan Liangliang |
|
| CB03 | Change of inventor or designer information |