Genome Threader文章学习

1.生物学知识

DNA由A、C、G、T四碱基组成，我们将每个DNA分子视为一个字符串。它可以很短-100，也可以很长-几百万。长的字符串代表物种的染色体，所有字符串组成该物种的基因组。这里只考虑蛋白质编码基因。遗传密码指的是核苷酸翻译成氨基酸的规则。pre-mRNA会经历一个splicing 过程。在这个过程中，内含子被切除，只有外显子留下来，形成mRNA。mRNA可以反转录形成cDNA，也可以形成ESTs。
接下来要介绍的gene finding包括把cDNA和ESTs比对到基因组DNA上，从而识别出基因的外显子和内含子。这个非常重要，因为你找到的比对并不一定会精准比对。因为自然变异和测序误差，match过程存在少许的单碱基错配和indel造成的差异。

2. Genome Threader解决的计算问题

最简单的情况下，Genome Threader的输入包括一个基因组DNA序列（通常比较长），一组cDNA/EST（可以是一个，也可以是几百万条序列组成的结合）。基因组DNA长达几百万个碱基，而cDNA/EST通常长500个碱基，最长大概2万个碱基。我们不知道会有多少个cDNA/EST能匹配到基因组的某个位置。所以，比对问题可以分成两个小部分。1.识别能与基因组DNA形成高质量match pair对的cDNA/EST 2.推导出最优比对-识别出exon和intron。
针对第一个问题，Genome Threader用了基于增强后缀矩阵的字符串快速匹配算法。
针对第二个问题，Genome Threader采用了经典动态规划算法。拿cDNA/EST或蛋白质倒推。目的是通过产物倒推基因结构。其实就是在容忍内含子存在的情况下，将产物比对到基因组上，也就是平常说的剪接比对-spliced alignment。基因组同一区域的不同剪接比对经常不同，也就是生理上的选择性剪接。因此，Genome Threader解决的最后一个问题就是推导出某一特定DNA区域的所有可能的转录本

2.1.基本概念

我们将序列堪称字母集合。定义序列的长度为。是序列的第a个碱基。如果，那么就是序列从到的那一段。如果，那么就是空序列。序列和的edit distance就是从到转换所需的最小插入、缺失或单碱基替换数。

2.2. 剪接比对问题

我们想把一条 cDNA/EST序列比对到基因组DNA。这两条序列都是A、C、G、T、N的集合。其中，N表示不确定碱基。上位置的碱基都有两种状态和。针对每个都会产生。其中。-表示缺失。针对每个都会产生。其中。举个例子，即表示基因组DNA缺失碱基，而表示cDNA缺失碱基，表示从上剪下来的碱基。
假如我们现在有一个序列，我们也不知道哪个碱基是内含子、哪个是外显子。接着假设是对应的输出。其中，是对应的输出。是和比对的结果。其实删去-和后，就是基因组DNA，就是cDNA

上面是基因组DNA，下面是cDNA。
形式的列是内含子（表示基因组里有，cDNA里没有），其他列是内含子。如果两个碱基匹配上了，用连接。插入和缺失都用表示。

因为想使比对最优化，作者设计了状态转移权重。
针对，有以下状态转移公式:
针对，有以下状态转移公式：

指基因组中缺失一个碱基的概率，指位置碱基处于基因组donor site的第一个位置的概率，指位置碱基处于基因组ac ceptor site的最后一个位置的概率。donor site是内含子的起始位点，acceptor site是内含子的终止位点。我们可以通过donor site和acceptor site对来定义内含子。因为和遵循贝叶斯剪切位点模型（BSSM），所以他们两个又被称为参数.
虽然状态转移矩阵的权重取决于，但是输出列的权重却是独立的。它的权重定义为：

是两碱基一致时的权重；是错配权重；是列中出现的权重，是缺失权重。对应列的权重就是0。把所有的状态转移权重和所有的输出列权重相加，得到权重总和。最优化问题就是寻找最大权重的，比对，即。当时，达到最优化。

3.计算最优剪接比对

跟很多生物学序列比较一样，剪接比对也可以通过动态规划DP解决。我们计算了两个的矩阵和。其中，是和的所有比对结果中以外显子结束的最大权重，是和的所有比对结果中以内含子结束的最大权重。显然，是

整体就是动态规划。不过为了使其尽可能合理，作者设置了很多参数。

4.内含子切除技术

预测脊椎动物或植物的基因结构时，经常会遇到长内含子的问题。有些内含子甚至长达10000到100000个碱基。如果硬用动态规划的话，需要的时间和空间就太多了。其实，内含子不会对剪接比对的权重产生太大影响。因此，在已知内含子位置的前提下，我们就可以跳过内含子的绝大部分，减轻计算量。按道理来讲，外显子部分应该与EST高度匹配，而内含子则不应与EST有任何的匹配。按这个想法，作者设计了相似性过滤策略。首先，找到基因组DNA和EST近似匹配的位置。如果几个match是一个剪接比对的几部分，就把它们合起来形成一条链。外显子将从基因组上的这些链中产生，而那些没被链覆盖的基因组DNA片段就是候选内含子。在用动态规划之前，先把这些内含子切掉；在动态规划的回溯阶段，再把这些内含子插进去。缺点就是如果一个外显子既不太长又不那么保守，可能EST match不上去，这种方法就会导致错误的基因结构

4.1. 计算match

计算match时，考虑和的最大近似匹配。其实最大近似匹配就是和的字串和的配对。这个配对满足最左和最右原则：最左即
或或；最右即或或。我们只对满足最小长度、最大距离的最大近似匹配感兴趣。最小长度指两个子串都不得短于这个长度；最大距离指两个子串的edit distance不得大于这个距离。关于edit distance，前面已经提到过了。
通常，采用seed-and-extend方法计算近似最大匹配。因为每个最大近似匹配里肯定会至少包括一个长的最大精准匹配。可以将这个精准匹配作为种子，将其向和两边延伸得到最大近似匹配。

4.2. 将match串起来

通常会有几个近似匹配以共线性的方式出现。因此，接下来，我们就要把这些匹配连起来。
定义match对。两者之间的gap函数为。当且仅当并且时，我们认为在，即。是用户自定义的最大gap宽度。仔细观察这个定义，我们会发现这些match之间可以有重叠区域。
定义重叠程度为。现在给定一个match集合和一条链，链内的match需要满足和。叫做起始片段，叫做中止片段。为了比较每条链，给每条链分配一个分数。链分数由链内的match决定。每个match都会有一个正分数。这个分数由和的最优比对决定，该比对中不包括exon和intron状态，也没有缺失碱基。最优比对中，匹配得2分，错配扣1分，插入和缺失都扣2分。简化后的得分公式为，是match的edit distance。链得分。我们想要寻找一条全局最优链，也就是具有最高得分的链。最直接的方法就是构建一个加权有向无环图，节点代表match，边代表，边的权重为。将求最优链的问题转换为求解图中具有最高权重的路径。求解方法跟运筹学中的最优路径求解类似。不同的是，作者寻求的是所有具有生物学意义的链。根据片段所属最优链的起始片段，将片段分类。如果两个片段对应最优链共享一个起始片段，那么这两个片段属于一类。每个类中都会有一条链，这条链具有用户自定义的最低EST覆盖度（默认是50%）。保留每个类中具有最高覆盖度的链，避免输出太多相似的链。这种办法可以让我们识别到与基因组不同区域匹配的链。

4.3.切割

给定EST和基因组DNA，我们可以得到全局链。为了留下临近内含子的剪接位点（动态规划需要），我们将链覆盖的基因组区域向两端各延伸几个碱基（用户自定义）。将那些重叠或离得非常近的延伸区域合并起来。动态规划将分析这些基因组子区域。按照动态规划的方向，把这些DNA子串串起来形成一个假的，也可以称为spliced gDNA。子串边界的数字（假设是）是他们在原来基因组上的距离。对spliced gDNA做动态规划。与没有内含子的动态规划不同的是，回溯到spliced gDNA的不同DP区域的边界时，将长的内含子插到剪接比对里。

5.计算保守剪接比对

由EST推导的剪接比对不能覆盖所有的基因，因为EST通常不超过500个碱基，而基因往往更长。为了得到完整的基因结构，我们需要连接同一区域的多个剪接比对。这样做可能会产生多个基因结构，这跟选择性剪接一致。所以我们不能简单的把剪接比对合起来。

根据上面的9个剪接比对得到下面的两个保守剪接比对。圈出来的短外显子证明该基因发生了选择性剪接。

假设手上有一个DNA 、一些EST序列和剪接比对结果。用在中的位置对代表一个剪接比对。同时，也需要记载中哪里是外显子，哪里是内含子。
如果并且，则认为和重叠。像前面一样，构建一个重叠图。其中，节点是中的剪接比对；边是剪接比对间的关系，有重叠就有边。假设重叠图是全连通的。（作者是这么解释的：假设手上有一些剪接比对结果，我们很容易就可以把它们分成不连接的子集，从而使得子集的重叠图全连通）
如果和重叠，并且重叠部分关于外显子和内含子的分配也是一致的，那么我们说这两个剪接比对是compatible。注意compatible关系不能递进，下图就是一个例子

和是compatible，和也是compatible。但和不是compatible的，因为他们两个关于内含子、外显子的分配不同

保守剪接比对问题其实就是寻找满足以下条件的最小集合:

• 
• 和不能重叠也不能compatible，，
• 对每个，就compatibility而言，是最大的