基于WGS和CBC测序策略的DNA序列拼接算法研究（三）
第一章 基因组测序与DNA序列拼接简介
2.1 基因组测序与序列拼接算法研究
2.1.1 基因组的概念
众所周知，生物由细胞组成，而细胞中含有一种与遗传有关的高分子化合物，称为脱氧核糖核酸（DNA）。DNA的基本组成单位为腺嘌呤、鸟嘌呤、胞嘧啶和胸腺嘧啶四种核苷酸，分别由A、G、C和T四个字母（称为碱基）表示，这些核苷酸按照一定的顺序和方向排列，就形成了DNA的序列。一般的DNA分子都有两条互补的分子链，两条链间通过A-T、G-C间的配对形成了双螺旋结构（如图2.1）。在解旋酶的作用下，DNA的两条链分离开，分别作为一个模板，又在聚合酶的作用下合成一条新链。这个过程称为DNA的复制（如图2.2）。
在解旋酶的作用下，DNA的两条链分离开，分别作为一个模板，又在聚合酶的作用下合成一条新链。这个过程称为DNA的复制（如图2.2）。
核糖核酸（RNA）是与DNA相似的一种高分子化合物，它的组成中以尿嘧啶（用字母U表示）替代胸腺嘧啶，并只形成单链结构。
通过一定化学作用下的A-U、G-C间的配对，一段RNA序列可以获取DNA部分离散区域（称为编码区或基因）的序列信息，这个过程称为转录。RNA序列中连续的三个字母可以确定唯一的氨基酸，而氨基酸的序列进一步确定了蛋白质的组成进而决定蛋白质的功能，这一过程称为翻译。
生命活动始终遵循从DNA序列到RNA序列，再从RNA序列到蛋白质这一规律，考虑到蛋白质是生命活动的基本单位，因此，DNA序列在某种程度上“主宰”了生命活动，这便是分子生物学中著名的中心法则（如图2.3），
它指出，生物活动都是通过由DNA序列转录成为RNA序列再由RNA序列翻译成为蛋白质这一过程进行的，即DNA序列蕴含了包括遗传在内的生物体所有活动的秘密。
基因是DNA上具有特定功能的一个片断，负责一种特定性状的表达。一般来讲，一个基因只编码一个蛋白质。
任何一条染色体上都带有许多基因，一条高等生物的染色体上可能带有成千上万个基因。一个细胞中的全部基因序列（包括编码区和非编码区）及其间隔序列统称为基因组（genomes）, 我们可以认为基因组就是生物所含DNA分子的全部序列。
可以看出，基因组在生物的生命活动中起着基础性作用，因此基因组DNA测序工作就显得尤为重要。
2.2 基因组测序方法简介
DNA快速测序的方法有两种：
1、链终止法：这一方法的原理是，通过合成与单链DNA互补的多核苷酸链来读取待测DNA分子的顺序，合成可在不同位置随机终止反应。
2、化学降解法：双链DNA分子被化学试剂处理，可在特定的核苷酸位点产生切口。用同位素标记测序碱基，以此确定DNA分子的序列。
由于链终止法更易自动化，所以，基因组的测序一般都使用链终止法。但是，一次测序反应只能测得几百个碱基。而复杂生物的DNA序列通常都在1M 碱基以上，有的甚至几十亿碱基。因此，不可能用测序仪一次测定全基因组的序列。现代大规模基因组测序普遍使用了shotgun和BAC方法.下面将就着两种算法的流程进行说明与介绍。
2.2.1 BAC简介
BAC算法简单来说就是把基因组打碎成200－300kb的片段，并制成BAC文库，再选择一些BAC进一步打碎成3kb左右的小片段，测序并拼接。
如图2.4是一个BAC测序的流程图：
2.2.2 WGS简介
WGS全称为whole-genome shotgun，它的作法是把基因组直接打碎成3kb左右的小片段，测序并拼接。并且，WGS在现在的测序项目中使用得越来越广泛。例如水稻基因的测序，就是使用的WGS策略。如图2.5，是一个WGS的测序流程图
2.3 拼接算法研究
2.3.1 为什么要拼接
目前，对于DNA测序的技术只能一次测序大约500左右长度的DNA序列，可是人类的DNA全长大约为30亿个，大约是600万个500。如果是能够顺序的进行分析，一次500的进行，这个也不是问题。只是一个重复的过程。但是，实际上的每次测序是在整个DNA片段中随机的一小段。怎样把这些随机的小片段拼接成一个完整的DNA序列，这个就是拼接过程需要解决的问题。
2.3.2 拼接算法的数学描述
DNA序列可以看成是字母集｛a, c, g, t｝上的字符串，对于字符s和序列A，称
为s的互补字符， 为A的逆互补序列，则
， ， ， 。
若 ，则
。
2.3 拼接算法研究
2.3.1 为什么要拼接
目前，对于DNA测序的技术只能一次测序大约500左右长度的DNA序列，可是人类的DNA全长大约为30亿个，大约是600万个500。如果是能够顺序的进行分析，一次500的进行，这个也不是问题。只是一个重复的过程。但是，实际上的每次测序是在整个DNA片段中随机的一小段。怎样把这些随机的小片段拼接成一个完整的DNA序列，这个就是拼接过程需要解决的问题。
2.3.2 拼接算法的数学描述
DNA序列可以看成是字母集｛a, c, g, t｝上的字符串，对于字符s和序列A，称
为s的互补字符， 为A的逆互补序列，则
， ， ， 。
若 ，则
。
2.2.2 WGS简介
WGS全称为whole-genome shotgun，它的作法是把基因组直接打碎成3kb左右的小片段，测序并拼接。并且，WGS在现在的测序项目中使用得越来越广泛。例如水稻基因的测序，就是使用的WGS策略。如图2.5，是一个WGS的测序流程图
2.3 拼接算法研究
2.3.1 为什么要拼接
目前，对于DNA测序的技术只能一次测序大约500左右长度的DNA序列，可是人类的DNA全长大约为30亿个，大约是600万个500。如果是能够顺序的进行分析，一次500的进行，这个也不是问题。只是一个重复的过程。但是，实际上的每次测序是在整个DNA片段中随机的一小段。怎样把这些随机的小片段拼接成一个完整的DNA序列，这个就是拼接过程需要解决的问题。
2.3.2 拼接算法的数学描述
DNA序列可以看成是字母集｛a, c, g, t｝上的字符串，对于字符s和序列A，称
为s的互补字符， 为A的逆互补序列，则
， ， ， 。
若 ，则
。
2.3 拼接算法研究
2.3.1 为什么要拼接
目前，对于DNA测序的技术只能一次测序大约500左右长度的DNA序列，可是人类的DNA全长大约为30亿个，大约是600万个500。如果是能够顺序的进行分析，一次500的进行，这个也不是问题。只是一个重复的过程。但是，实际上的每次测序是在整个DNA片段中随机的一小段。怎样把这些随机的小片段拼接成一个完整的DNA序列，这个就是拼接过程需要解决的问题。
2.3.2 拼接算法的数学描述
DNA序列可以看成是字母集｛a, c, g, t｝上的字符串，对于字符s和序列A，称
为s的互补字符， 为A的逆互补序列，则
， ， ， 。
若 ，则
。
定义1：若A、B为两条DNA序列，如果我们将从A转换到B所需的插入（insertion）、删除（deletion）和置换（substitution）操作的次数的最小值称作编辑距离（edit distance），记作d(A, B)。
片段拼接问题可以描述为确定一个包含了所有片段的最小的DNA字符集上的字符串，问题的数学模型描述如下：
定义2: DNA序列拼接问题(the DNA sequence reconstruction problem) ，给定一个片段序列的集合F和一个误差
， ，寻找一条最短的序列S，对于任意一个片段序列A，A∈F，在S中都能找到一个子序列B，
DNA序列拼接问题与组合数学中的最短超串问题相似。最短超串问题即给定一个字符串的集合，找出一最短的字符串称为超串，使得它将集合中的任何一元素作为其子串。在
＝0的情况下，DNA序列拼接问题实质上就是最短超串问题。
由于序列拼接问题的复杂性以及在生物信息学中的特殊重要性，人们对这一问题进行了积极的研究，并提出了许多算法。总的来说，序列拼接算法可以大致分成两类：一类是利用Hamiltonian path 的方法，另一类是利用Eulerian path 的方法。两类算法都是采用了基于图论的近似方法，其主要思想是根据片段间的重叠信息构造某种类型的图，然后利用基于图论的一些算法得到结果，以下几节将对这个问题加以讨论。
DNA序列拼接问题与组合数学中的最短超串问题相似。最短超串问题即给定一个字符串的集合，找出一最短的字符串称为超串，使得它将集合中的任何一元素作为其子串。在
＝0的情况下，DNA序列拼接问题实质上就是最短超串问题。
由于序列拼接问题的复杂性以及在生物信息学中的特殊重要性，人们对这一问题进行了积极的研究，并提出了许多算法。总的来说，序列拼接算法可以大致分成两类：一类是利用Hamiltonian path 的方法，另一类是利用Eulerian path 的方法。两类算法都是采用了基于图论的近似方法，其主要思想是根据片段间的重叠信息构造某种类型的图，然后利用基于图论的一些算法得到结果，以下几节将对这个问题加以讨论。
2.3.3 基于Hamilton路径方法的拼接算法
基于Hamilton路径方法的拼接算法基本思想是：首先将所有的Read构成一个有向图G，每个Read看成一个结点，如果两个read之间存在有重叠Overlap，那么在相应的结点之间就存在有一条边。然后通过寻找经过每个Read一次且仅一次的一条路径，就将序列拼接问题转化成Hamilton path 问题。这种方法可以分为如下三步：
1）找出序列片段间的重叠信息；
2）将存在重叠的片段组合起来，形成一个Contig结构；
3）根据片段中每个碱基的质量值，在Contig结构中寻找一条最终序列，称作“Consensus”序列。
这一类算法主要有Phrap,TIGR,CAP3/CAP4,GigAssemble等，它们都遵循”Overlap-Layout -Consensus”三步曲，但是都有一些细微的差别和特点。
以Phrap为例，Phrap算法的三个步骤描述如下：
(1) Overlap：如果两个Read（序列片段）有Overlap（重叠区域），那么它们之间必定有一定长度的精确匹配区域，称为Match。
因此算法首先寻找存在长度至少为某个阈值的Match的所有Read Pair（两个相关的片段），对每个Read Pair来说，一个Match就在Smith－Waterman比对矩阵上定义了一条对角线，将这条对角线扩展成Band（矩阵上的一个带状区域），对这个Band运行Smith－Waterman算法，如果Smith－Waterman得分超过某个阈值，则认为Overlap。Phrap寻找Match的算法是：首先建立一个指针集合，包含指向每个Read中的每个字符位置的指针 ，然后按照字母序列对集合中指针排序 ，最后在相邻的指针中寻找Match。
(2) Layout：在这一步中，Phrap首先对Overlap进行打分，衡量两个Read之间的不同是真不同，还是由于测序错误造成的，然后将所有Overlap按照该打分值由高到低进行排序。
该步操作中，Phrap算法定义了一种称为Contig的结构来表示一组片段的组合关系。Contig 包含一个片段集合F，这些片段最终将形成一条模糊的序列S（因为某些重叠并不是完全匹配的），而每一条片段相对于S起始位置都有一个偏移量。片段在进行组合时，每个片段先形成一个单独的Contig，当发现某两个片段fi和fj存在有重叠（可以形成一个实际的拼接时），Phrap算法将Contig(fi)和Contig(fj)进行合并(称为Merge)，形成一个新的结构Contig(fi, fj)。我们用Ffi表示contig(fi)中的片段集合，用Sfi表示Ffi所形成的序列。
在将Contig(fi)和Contig(fj)进行合并形成Contig(fi, fj)时，首先需要对原结构Contig(fi)和Contig(fj)的片段集合Ffi和Ffj进行并操作形成Contig(fi, fj)中的片段集合F(fi ,fj)，对原结构所表达的序列Sfi和Sfj通过合并重叠形成新的序列S(fi ,fj)，并重新计算Contig(fi, fj)中每条片段相对于S(fi ,fj)的偏移位置，最后还需要删除原来的结构Contig(fi)和Contig(fj)。这样，Phrap算法将存在重叠的片段组合起来形成一个Contig结构。
(3) Consensus：最终的Consensus序列将由各个Read的高质量部分来构成。
在将存在重叠的片段形成一种组合关系后，Phrap利用一种称作“Mosaic”的方法将Contig结构中一些高质量的序列片段连接起来形成一条Consensus序列。使用这种方法需要先在Contig结构上构造一个加权有向图G，将求解Consensus序列问题转化为在G中寻找一条从起始节点到终止节点的权重最大的路径。加权有向图G的具体构造方法如下：
1）节点：在Contig结构的片段集合中，定义偏移量最小的片段的最左端碱基字符为图的起始节点，偏移量最大的片段的最右端碱基字符为图的终止节点，从起始节点开始每间隔若干个碱基字符定义一个图的节点。
图2.6 加权有向图G示意图
2）边：在同一条片段的两个相邻节点之间定义一条单向边，边的方向与起始节点到终止节点的方向相同，该边的权重为两个节点间所有碱基字符所对应的质量值之和。在两条片段重叠区域的相对应节点间定义一条权重为0双向边。如图2.6所示：
Phrap采用Tarjan算法来寻找权值最大的路径，虽然该算法的时间复杂度和空间复杂度均为问题规模的线性函数，但是由于加权有向图G对内存的需求较大，尤其在处理大规模数据时该步的运行时间以及对内存的需求量都是非常庞大的。因此，如何提高Phrap的运行时间以及减少程序对内存空间的需求就显得尤为重要。
2.3.4 基于Euler路径方法的拼接算法
此类拼接算法是2001年由California大学计算机系的Pavel A. Pevzner等提出的，他们认为传统的“重叠－排列－生成共有序列”的思维模式导致了将拼接问题抽象为NP-完全的Hamilton路径问题，他们从另一种从来没有在实际测序工程中应用但是直接导致了基因芯片工业的杂交测序方法SBH（Sequencing by Hybridization）中得到启示，提出了容易计算的在de Bruujin图中寻找Euler超路径的拼接方法，并用研制的软件得到了差错和Repeat匹配错误均优于Phrap的结果。
以下介绍此类拼接算法的基本思想。
1. 对每个片段fi，设它的长度为ni，则将fi转化为一个含有ni＋1－l个相互重叠l-子串的集合F(fi)。即，从fi的第一个字母起，依次取以此字母开头的l-子串，这样相邻两次取到的l-子串至少有长度为l－1的重叠区域。
2. 依次将每个F(fi)合并进de Bruujin图G，其中G的节点为每个l-子串，边集为每两个有长为l－1的重叠区域的l-子串有序对。于是，每个F(fi)可以看作是依次通过某个l-子串的路径，并且不同的F(fi)、F(fj)可通过他们所含的公共节点和节点间的边相联系。
通过2. 拼接问题转化为Euler 超路径问题，即给定一个图G及此图中一些路径的集合
，给出图G的一个Euler路径P（称为Euler超路径），使得它将F中的每一条路径都作为自己的子路径。
对系统<G, F>进行一系列“等价”的T变换
其中对T变换的最基本要求就是变换前后两个系统的Euler超路径相同。最终得到的变换结果<Gk, Fk>就提供了所需要的拼接结果。
该类拼接算法以一种新颖的观点解决片段拼接问题，在剔除重复子序列及矫正片段数据中错误方面有其独到之处，这一类拼接算法最著名的是EULER，下面详细描述之。
以EULER为例
EULER将拼接问题归结为寻找Eulerian SuperPath问题， 主要特点有：1）不尝试各个Read之间的Overlap，没有Overlap步骤；2）将Read切割成规整的小片段。
EULER算法流程：
Step1：
1）将所有的Read切割成小片段（即k-子串），以排除Read中的错误，获得准确度较高的Read 。思路是：如果知道该Read的原始序列G，那么直接使用Read和G进行比较，但是不知道G，那么使用G的k-子串集合Gk亦可，但是Gk亦不可求，所以使用Gk的近似。所有出现在至少m（阈值）个Read中的k-子串形成的集合称为Gk的近似。
2）将每个Read和Gk的近似进行比对，寻求Read的最小改变使得Read的所有k-子串在Gk的近似中。
Step2：构造de Bruujin图。
结点：每个k-子串是结点 。
边： 结点v和u有一条边，当前仅当v的尾部和u的头部相同。
图的构造结果是：每个Read表示成一条Path，每个Repeat表示成一个多入口、多出口的单一链，但是不知道出入口之间的对应关系。
Step3：通过图的等价变换.，在de Bruujin图中寻找一条Eulerian Path，使得能够覆盖所有的Path。
解决DNA序列拼接问题的关键点在于如何提高拼接速度和准确度。事实上，影响拼接速度的主要因素是片段（Read）数目巨大、片段比对操作非常耗时，影响拼接准确度的主要因素是测序错误和重复序列（Repeat）的存在。目前的两类拼接算法都使用了图这一复杂的数据结构记录拼接的中间信息，无论是存储开销还是时间开销都是非常大的。另外，两者还需要保留大量的临时和过渡数据，特别是第二类拼接算法，在将片段分割成更小的串后以及进行图变换时，存储的要求将数倍于第一类拼接算法。在时间开销上，虽然Euler路径的计算复杂性低，但其图节点数百倍于Hamilton图，抵消了一部分这方面的优势。
在基于Hamilton路径的拼接算法中，人们把每个片段看成一个图的顶点，两个片段有重叠就连一条边，所有片段根据重叠信息排列，算法的目标是寻求一条合理的Hamilton路径。在基于Euler路径方法的拼接算法中，人们把每个片段看成一个路径，许多相同的路径组成了de Bruujin图的边，算法的目标是寻求一条合理的Euler路径。（本小结内容参考《DNA序列拼接的分布式并行处理》方小永 骆志刚 国防科学技术大学）
