数学建模基因重组.docx
- 文档编号:27843983
- 上传时间:2023-07-05
- 格式:DOCX
- 页数:32
- 大小:322.92KB
数学建模基因重组.docx
《数学建模基因重组.docx》由会员分享,可在线阅读,更多相关《数学建模基因重组.docx(32页珍藏版)》请在冰豆网上搜索。
数学建模基因重组
基于deBruijn图的DNA重组模型
摘要
本文在运用MATLAB编程的基础上,运用了deBruijn图模型设计出相应算法对基因小片段进行重组,以达到组装序列的连续性、完整性和准确性。
首先,由于测序中可能出现的个别碱基对识别错误,两个read间最小重叠长度,最小正确率重叠关系,每都要比对一次。
这时,要设定好如下几个参数:
最小重叠长度,最小正确率。
最小重叠长度是指两个read之间的重叠部分的长度的最小值,如果两个read之间的重叠长度小于最小重叠长度,则被认为这两个read是没有重叠的。
由于测序数据存在错误,故最小正确率是必要的,即允许一定数量的碱基不匹配,以纠正测序过程中存在的错误。
当然,这也可能导致发现错误的重叠关系,即两个read本来没有重叠,但在允许错误的条件下被误认为它们有重叠关系。
问题二:
基因组测序是生物信息学中最基本的研究方向之一,然而大多数生物的基因组都不可能一次性获得,需要利用序列拼接技术对实
验中获得的DNA片段进行拼接操作。
目前,测序过程中获得的DNA片段越来越短,基于Euler路径的拼接算法在处理这种短片段拼接时具有
优势。
在Euler路径算法中,一个关键的步骤是deBruijn图的构建,一直以来,构建deBruijn图的方式总是让后一个k-met与前一个k-mer
之间有一1个碱基的交叠,相邻的两个k-met之间相互错开一位。
但文中的研究发现,如果有边连接的两个k—mer之间有一2个或者更少的
碱基相交叠,会对deBruijn图结构复杂性产生重要影响。
针对这些影响进行详细分析,并设计实验进行验证,实验结果表明,k-mer之间的错
位数变化对deBmijn图结构复杂性有显著影响。
提出一种新的deBruijn图,并且引入了决策表的概念,通过决策表里的信息来选取后继k-mer,并在适当的时候更新决策表。
值得一提的是,本文在已有的成熟理论和模型的基础上加以应用和改进,形成可推广于研究短片段DNA重组模型。
并使用MATLAB软件使得模型更具合理性,可行性和科学性。
关键词:
重叠图模型,deBruijn图模型,决策表,重叠片段优化
一、问题重述
生物生长发育的本质是遗传信息的传递和表达,基因组包含了生物体的遗传信息。
对于生物体而言,想要探索其遗传信息的奥秘,就必须对基因组中基因间的相互作用、基因的结构及其功能等方面进行全面研究,以此获得生物体基因组的全部序列,来探索生命的本质。
因此基因组测序重组之于生物信息学领域有重要的研究意义。
基因组测序和重组的应用极其广泛,测序产生的读取片段reads经过序列拼接组装,生成基因组的碱基序列。
随着时代的进步与科技的发展,基因组测序技术从第一代发展到普遍应用的第二代,并在近年来兴起了第三代基因组测序技术。
基因组测序的发展过程在逐步中向着高通量、低成本的方向前进。
新一代测序技术正在引领生命科学研究进入一个崭新的阶段,为基因组重组的实现提供了不少进步。
与之同时,我们也面临着不小的挑战。
在这样的背景下,我们提出问题:
问题一:
建立数学模型模拟读长序列重组基因组。
设计将读长序列组装成基因组的算法并尽可能使其取得良好效果。
结合设计出的组装算法,编制程序,利用计算机将模型和算法的结果得以体现。
并在设计过程中考虑到可能出现的个别碱基对识别错误、基因组中存在重复片段等问题,对建立的模型或算法进行改进。
问题二:
利用题目所给定的Hiseq2000测序仪进行测序,尝试运用问题一中提出的算法对现有的一个全长约为120,000个碱基对的细菌人工染色体(BAC)进行程序组装,如何才能对基因组进行最好的组装?
是否需要对读长序列的组合算法进行进一步的改进?
二、问题的分析
在拼接中,当前拼接片段与read关系显然,当前拼接的碱基片段必有有如下性质:
对于碱基片段上任何一个碱基,都存在一条成功拼接read,使得该碱基出现在该read上。
即碱基片段上每一个碱基都能被一条成功拼接read所覆盖。
现在要找一个U的子集S及一个碱基序列g,满足如下条件:
(1)给定阈值quality和contentA,S中的read上碱基质量平均值不小于quality,read上不含未知碱基,且read上碱基A的含量不超过contentA;
(2)S中所有read都是成功的,U-S中所有read都不是成功的;
(3)S中的read尽可能的多;
(4)g满足性质1。
此时,g就是最终想要得到的基因组。
当然希望g与G越接近越好,最好的结果就是g与G相同。
实际上,如果没有条件(3)的限制,则满足条件的S和g有很多,若S中read
较少,就会导致g和原基因组差别较大,原因是没有充分利用测序仪所测出的read。
由于U中的read是能够完整覆盖原基因组的,故理想情况下,S中的read也能完
整覆盖原基因组,即使是这样,也不能保证g与原基因组完全一样,主要原因是
基因组中存在大量重复片段。
另外,测序过程中的错误也不容忽视。
设ra,rb是S
中的两条完全一样的read,则有如下几种可能:
rb是ra的复本(或ra是rb的复本);rb不是ra的复本,他们是来自重复片段的两条read;rb中包含测序错误,ra中无错误,导致rb碱基与ra相同;ra中包含测序错误,rb中无错误,导致rb碱基与ra相同;ra和rb中都有错误,它们的碱基相同纯属巧合。
所以,现在能做的是使碱基序列g尽可能与原始基因组一致,存在误差在所难免,本文要使该误差尽可能的小。
其中,集合S的选取很重要,实际上,S需要在拼接过程中动态生成,每当有成功的read出现,就把它并入S中。
从而生成策略的好坏会直接影响S中read数量以及g上的每一个碱基。
设R是U中所有满足条件
(1)的read,显然S是R的子集,故S的基数不会超过R的基数。
实际上,要想让R中所有read成功参与拼接是很困难的,以至于无法在有效的时间内判断是否存在一个拼接策略,使得R中所有read都成功参与拼接。
但先可以先构建集合R,采用启发式搜索策略,让R中的read逐条参与拼接。
三、基本假设
1.如果两个read之间的重叠长度小于最小重叠长度,则被认为这两个read是没有重叠的。
2.进行基因重组时,初始read肯定是正确的。
3.假设每个k-mer长度为k,对一个read划分时,先以read的任意一端为起始位置,截取k个碱基,再将起始位置向后移动一个碱基,再截取个碱基,依此类推,直至截取得到的k-mer的尾部到达read的另一端,这些k-reeF组成了deBruijn图中的顶点。
四、主要变量符号说明
为了便于描述问题,我们用一些符号来代替问题中涉及的一些基本变量,如表1所示。
表1主要变量符号说明一览表
质量值
碱基被测错的概率
多重集
基因组
【注】其余没有列出的符号,我们将在文章第一次出现时给出具体说明
五、问题模型的建立和求解
5.1read预处理
5.1.1ASCII码转换
从计算机语言上,对数字识别比对字符串识别容易,因此对于需利用ASCII码表(见表1)对每条read上的字符串进行转换。
表1ASCII码字符表
Dec
缩写字符
Dec
缩写字符
Dec
缩写字符
Dec
缩写字符
Dec
缩写字符
Dec
缩写字符
33
!
49
1
65
A
81
Q
97
a
113
q
34
"
50
2
66
B
82
R
98
b
114
r
35
#
51
3
67
C
83
S
99
c
115
s
36
$
52
4
68
D
84
T
100
d
116
t
37
%
53
5
69
E
85
U
101
e
117
u
38
&
54
6
70
F
86
V
102
f
118
v
39
'
55
7
71
G
87
W
103
g
119
w
40
(
56
8
72
H
88
X
104
h
120
x
41
)
57
9
73
I
89
Y
105
i
121
y
42
*
58
:
74
J
90
Z
106
j
122
z
43
+
59
;
75
K
91
[
107
k
123
{
44
60
<
76
L
92
\
108
l
124
|
45
-
61
=
77
M
93
]
109
m
125
}
46
.
62
>
78
N
94
^
110
n
126
~
47
/
63
?
79
O
95
_
111
o
127
DEL(delete)
48
0
64
@
80
P
96
`
112
p
5.1.2正向“完整”read
DNA分子由两条单链组成,在图中表现为两条平行直线,两条直线上相对位置的两个碱基相互结合形成碱基对(bp),并且与碱基A结合的碱基必为T,与碱基C结合的碱基必为G。
将一个含120,000个bp的完整基因组,随机打断成500bp的片段,然后对500bp的片段进行测序。
分别从500bp片段的两端,对两条单链进行测序,测得的读长记为reads1,reads2。
reads1,reads2的长度均为88bp。
基因重组需要在一条同一链上完成,因此我们需要对read1补充成一条头尾两端用88bp碱基对表示的“完整链”。
根据DNA碱基对的互补配对准则,我们利用read2对read1尾部进行互补(如图1),然后只研究一条链的重组。
图1利用read2对read1尾部进行互补过程图
5.1.1read选取
由于一些read错误率较高,不宜直接进行基因重组,故需要现将这些read过滤掉再进行基因重组。
有一下三种类型的read需要被滤掉:
类型一:
质量值低于质量阈值
数据文件中不仅保存了read的碱基,还保存了read的质量值。
每个read的每个碱基都有一个质量值,该质量值能反映该碱基的正确率。
质量值(
)越高,则碱基被测错的概率(
)越小,其计算公式为
对于每条read都是从5’端测到3’端,其中5’端的碱基正确率较高,而3’端正确率稍低。
当然,我们要事先指定一个质量阈值,选择那些质量较高的read。
碱基的质量值13,错误率为5%,20的错误率为1%,30的错误率为0.1%。
根据资料显示,Illumina官方(基因组学研究领域的技术与市场领导者)一般以Q30作为评价标准。
因此,本文也选取Q30作为标准,当该read质量值低于30,则被剔除。
类型二:
含有未知碱基
文件中有的read包含未知碱基,这是由于测序过程中有些碱基没有被测出导致的。
未知碱基用N表示,故我们要过滤掉那些含N的read。
类型三:
碱基A含量过高
有些read上碱基A含量过高(超过90%),甚至所有碱基都是A。
据测序仪性质,这样的read错误率较高,在拼接过程中应尽量避免使用该类read。
最终我们用所有符合上述条件的read导入优良read备选库S进行后续运算。
5.2基于重叠图算法
根据DNA测序过程,设G是一个字符串,称之为基因组。
取G的N个复本,然后将每个复本在随机的位置断开,形成一些小的字符串片段。
接着扔掉一些稍长的和稍短的字符串片段,将剩余字符串片段切成定长500bp的字符串保存在多重集U中,称U中的字符串为read。
字符串上的每个字符都来自字符集{A,C,G,T},每个字符就是一个碱基。
这些字符串长度都是L=500bp。
同时,每个字符串上的每个字符都有一个正确率,并且这些字符串能完整的覆盖原始基因组G。
现在要利用集合U中的read拼接成原基因组。
有时候,需要使用一条read的反向互补序列,即将read上每个碱基互补之后再将所得到的字符串反转。
实际上,当前拼接的碱基片段始终处于DNA的上链,但一条read可能来自DNA的上链或下链。
故一条read本身及其反向互补序列都要参与拼接过程,称read本身的序列为正向read,其反向互补序列为反向read。
不管是正向read还是反向read,其碱基方向都是从5’到3’的。
当然,多数情况下read上的碱基错误率很小(小于万分之一),只有每个read的3’处几个碱基出错的概率相对高些,故正向read的3’处碱基错误率较高,反向read的5’处碱基错误率较高。
这些碱基可能被丢弃,此时实际上只有该read的一个前缀或后缀在参与拼接。
若read上参与拼接的碱基长度大于事先指定的阈值minL,则称该前缀(后缀)为成功拼接前缀(后缀),称该read为成功拼接read。
所以,当前拼接的碱基片段是由成功拼接read相互重叠而成(如图)。
图2当前拼接片段与read关系
当前拼接片段与read关系显然,当前拼接的碱基片段有如下性质:
对于碱基片段上任何一个碱基,都存在一条成功拼接read,使得该碱基出现在该read上。
即碱基片段上每一个碱基都能被一条成功拼接read所覆盖。
现在要找一个U的子集S及一个碱基序列g,满足如下条件:
(1)给定阈值quality和contentA,S中的read上碱基质量平均值不小于quality,read上不含未知碱基,且read上碱基A的含量不超过contentA;
(2)S中所有read都是成功的,U-S中所有read都不是成功的;
(3)S中的read尽可能的多;
(4)g满足性质1。
此时,g就是最终想要得到的基因组。
当然希望g与G越接近越好,最好的结果就是g与G相同。
5.2.1重叠图三阶段
基于重叠图的算法主要有3个阶段:
overlap-layout-consensus。
在overlap阶段,要发现所有read之间的重叠关系,每两个read都要比对一次。
这时,要设定好如下几个参数:
最小重叠长度,最小正确率。
最小重叠长度是指两个read之间的重叠部分的长度的最小值,如果两个read之间的重叠长度小于最小重叠长度,则被认为这两个read是没有重叠的。
由于测序数据存在错误,故最小正确率是必要的,即允许一定数量的碱基不匹配,以纠正测序过程中存在的错误。
当然,这也可能导致发现错误的重叠关系,即两个read本来没有重叠,但在允许错误的条件下被误认为它们有重叠关系。
在layout阶段,建立整个重叠图。
由于重叠图本身可以不存放碱基序列,仅仅表示read之间的重叠关系,故可将整个图都读入内存。
在consensus阶段,需要遍历重叠图,图中每一条路径都有可能成为一条contig。
在得到所有contig之后,要进行contig之间的组装,这时要利用配对数据将contig组装成supercontig,从而得到整个基因组的碱基序列。
5.2.2贪婪算法
贪婪算法是通过一系列的选择来得到一个问题的解,而它所做的每一次选择都是当前状态下某种意义的最好选择,即贪婪选择。
贪婪算法的基本思路是从问题的某一个初始解出发一步一步地进行根据某个优化测度,每一步都要确保能获得局部最优解。
每一步只考虑一个数据,他的选取应该满足局部优化的条件。
若下一个数据和部分最优解连在一起不再是可行解时,就不把该数据添加到部分解中,直到把所有数据枚举完,或者不能再添加算法停止。
基本贪婪算法的实现过程:
1.从问题的某一初始解出发;
2.while朝给定总目标前进一步do;
3.求出可行解的一个解元素;
4.由所有解元素组合成问题的一个可行解。
从穷举法角度,所有的read小片段已完全覆盖整个基因组,因此对于任意一个read进行左右延伸,必然可以得到一条包含此read的真实基因组,记为G。
但是对于穷举,这是一个耗时的过程,不可行。
本文给出了一个基于贪婪算法的DNA重组算法(见图二)。
步骤一:
任取一个read,发现所有与该read右端有重叠关系的其他read;
步骤二:
选取后续read
选取与该read右端有重叠关系的其他read中最大的一个read’,使用好read’后并将其剔除read库;
步骤三:
对read’进行步骤一。
步骤四:
直至无法向右端继续延伸,重复步骤一对左端进行向左延伸。
图3贪婪算法图解
5.2.2算法改进
步骤一:
将导入的read按一定条件(该段是否有未确定的碱基N、碱基A占该段碱基序列的百分比是否大于90%、该段的的质量值是否大于等于30)进行筛选,得到筛选后的read库。
步骤二:
对筛选完的read库均取第一段碱基序列(数据库中第一行)作为连接的初始碱基序列。
步骤三:
在库中的其他碱基序列(剔除基因组序列库中已含的碱基序列(初始基因组序列库中只含第一段碱基序列))中进行检索与前一段碱基序列重合部分碱基数大于等于一定数值(暂定为10)的碱基序列。
步骤四:
在检索出的碱基序列中选取其中序列最长的作为下一段碱基序列并放入基因组序列库。
步骤五:
回到第三步直至基因组序列库中的序列无法再向两端延伸。
5.3DeBruijn模型
利用reads之间的重叠,通过公共路径的方法可以解决拼接问题。
而新一代测序产生的数据read更短、覆盖度更高、序列精度较低,为此这种“read为中心”的方法面临海量计算的困境,似乎不可能找到恰当的启发式方法来处理大量的重叠。
DeBruijn图框架为处理高覆盖、短序列提供了很好思路,该框架借鉴了Pevzner和Waterman等人针对传统的长reads提出的欧拉遍历方法,并在此基础上针对新一代测序数据的特点进行了改进。
第一个利用此技术的拼接算法是Newbler,之后陆续出现了一些直接或间接使用DeBruijn图的拼接算法,不同的方法在处理测序错误方式和使用reads信息的程度上(例如是否使用mate-pair信息)有所不同。
最新提出的拼接算法还考虑了覆盖率信息。
reads之间比对拼接算法是影响整个测序速度的和质量的关键因素。
目前新一代测序数据用于从头测序的短序列拼接组装算法普遍采用deBruijn图数据结构。
在deBruijn图上,每一个k-mer都构成图的节点,如果两个k-merread中相邻,那么这两个节点之间就有一条边。
reads集合中的每个read都对它所含的节点和边加权,这样reads集合产生一个节点和边都具有权值的deBruijn图。
在存储每一个k-mer时,往往要建一个无冲突的哈希表,以加快查找速度。
而建立哈希表可能会消耗更多的内存。
但是,由于每个k-mer在哈希表中只存储一次,不管该k-mer在read中出现了多少次,所以实际消耗的内存小于存储所有read所需要的空间。
另外,基因组中的重复片段会在deBruijn图中产生环路。
环路将在遍历deBruijn图时产生障碍。
5.3.1基本算法
基于deBruijn图数据结构的reads之间比对拼接算法可概括为以下几个步骤。
1)从read库中筛选read。
由于一些read错误率较高,不宜直接拿来建立de
Bruijn图,故需要现将这些read过滤掉。
2)确定k值,建立deBruijn图。
这时需要扫描所有read数据,将每一个长为L的read拆分成L-k+1个kmer,并用所有read的所有k-mer来累加,建立节点和边都加权的deBruijn图;
3)化简deBruijn图,连续线性延伸节点合并为单一节点,产生一些碱基序列更长的节点;
4)错误校正,删去由于测序错误产生的尖端和泡状结构;
5)通过read的配对末端(pair-end)、环化配对(mate-pair)信息伸展或者删去一些环;
6)依据环上节点和边的权值(覆盖深度信息)进一步伸展或者删去一些环;
7)遍历deBruijn图产生contig。
以上是对基于deBruijn图的算法做了简单介绍,为处理新一代测序仪所产生的海量测序数据,deBruijn图应运而生。
实际上,deBruijn图是一种特殊的加权图,不仅图的结点上有权值,而且图的边上也有权值。
化简deBruijn图是非常关键的一个步骤,通过对deBruijn图化简,可降低算法的时间复杂性以及空间复杂性,同时可以保证错误校正顺利的进行。
为得到质量较好的contig,必然会对deBruijn图进行错误校正。
使用deBruijn图只能将较短的read拼接成较长的contig,要得到全基因组,还需要contig的组装过程。
提出一种新的deBruijn图,并且引入了决策表的概念,通过决策表里的信息来选取后继k-mer,并在适当的时候更新决策表。
具体结构如下图
图4决策过程
5.3.2deBruijn图建立
我们要将测序仪得到的read合理的保存在内存中。
并且给定k-mer后,要快速的找到该k-mer出现在哪些read的哪些位置上。
由于read很多,我们不可能将read的碱基保存在内存中,故要给每个read一个编号,称为readID,我们只需保存readID即可。
如果确实需要保存碱基时,我们将碱基按照如下规则转化成整数:
A-00,C-01,G-10,T-11(其中整数都是二进制的)。
我们设计了如下的数据结构,称为readID_pos数组,每个k-mer都关联一个readID_pos数组。
其中readID是我们给每条read的一个编号,在debruijn图中,一个readID能唯一标识一条read。
后面我们会看到,在决策表中要想唯一标识一条read,需要二元组(readID,ori),即决策表中的read是有方向的,而debruijn图中的read没有方向。
readID是升序保存的,便于以后的二分查找操作。
pos是该k-mer出现编号为readID的read中的位置。
另外,我们要经常删除拼接成功的read,故在结构中保存了删除标记。
由于我们可能误删了某个read,所以使用删除标记的好处是便于我们及时恢复误删的read。
我们将所有k-mer存在一个哈希表中,取k=12,哈希表中共有16M个元素。
该哈希表实际上是一个超大数组,以每个k-mer的哈希值作为下标就可找到该k-mer的所有信息。
数组中每个元素都是指向如下所示的数据结构的指针,该数据结构就是k-mer结构。
其中,kmer_seq是转化成整数之后的k-mer碱基,实际上该整数就是该k-mer的哈希值。
。
addr是该k-mer所对应的readID_pos数组的首地址,num的是readID_pos数组中元素个数。
addr和num可唯一标识一个数组。
cur记录了readID_pos数组中下一个空闲的位置,因为填写readID_pos数组时不是一次性从头填到尾的,在遍历数据文件时,其中的的数据会随机地落到某个readID_pos数组中,故每一个数组都要对应一个cur。
deBruijn图建立完成后,cur会有其它用途
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学 建模 基因 重组