Genome Threader文章学习

作者: 扇子和杯子 | 来源:发表于2020-05-05 18:46 被阅读0次

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.基本概念

我们将序列堪称字母集合。定义序列 $s$ 的长度为 $|s|$ 。 $s[a]$ 是序列 $s$ 的第a个碱基。如果 $a \leq b$ ，那么 $s[a...b]$ 就是序列 $a$ 从 $a$ 到 $b$ 的那一段。如果 $a \geq b$ ，那么 $s[a...b]$ 就是空序列。序列 $s$ 和 $s^{'}$ 的edit distance就是从 $s$ 到 $s^{'}$ 转换所需的最小插入、缺失或单碱基替换数。

2.2. 剪接比对问题

我们想把一条 cDNA/EST序列 $c=c[1...m]$ 比对到基因组DNA $g=g[1...n]$ 。这两条序列都是A、C、G、T、N的集合。其中，N表示不确定碱基。 $g$ 上 $t$ 位置的碱基都有两种状态 $ex_t$ 和 $in_t$ 。针对每个 $ex_t$ 都会产生 $\begin{bmatrix}\alpha\\\beta \end{bmatrix}$ 。其中 $\alpha,\beta\in \{A,C,G,T,N,-\}$ 。-表示缺失。针对每个 $in_t$ 都会产生 $\begin{bmatrix}\alpha\\\cdot \end{bmatrix}$ 。其中 $\alpha\in \{A,C,G,T,N\}$ 。举个例子，即 $\begin{bmatrix}-\\\beta \end{bmatrix}$ 表示基因组DNA $g$ 缺失碱基 $\beta$ ，而 $\begin{bmatrix}\alpha\\- \end{bmatrix}$ 表示cDNA $c$ 缺失碱基 $\alpha$ ， $\cdot$ 表示从 $gDNA$ 上剪下来的碱基。
假如我们现在有一个序列 $Q=q_1,q_2,...,q_k$ ，我们也不知道哪个碱基是内含子、哪个是外显子。接着假设 $A=\begin{bmatrix}\alpha_1 \alpha_2 ...\alpha_k\\\beta_1 \beta_2 ...\beta_k\\\end{bmatrix}$ 是对应的输出。其中， $\begin{bmatrix}\alpha_i\\\beta_i\end{bmatrix}$ 是 $q_i$ 对应的输出。 $Q和A$ 是 $c$ 和 $g$ 比对的结果。其实删去-和 $\cdot$ 后， $\alpha_1,\alpha_2,...\alpha_k$ 就是基因组DNA $g$ ， $\beta_1,\beta_2...\beta_k$ 就是cDNA

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

因为想使比对最优化，作者设计了状态转移权重。
针对 $t\in[1,n-1]$ ，有以下状态转移公式:
针对 $t\in[1,n]$ ，有以下状态转移公式：

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

$\sigma$ 是两碱基一致时的权重； $\mu$ 是错配权重； $v$ 是列中出现 $N$ 的权重， $\delta$ 是缺失权重。 $in_t$ 对应列的权重就是0。把所有的状态转移权重和所有的输出列权重相加，得到权重总和 $w(Q,A)$ 。最优化问题就是寻找最大权重的 $g$ ， $c$ 比对，即 $w(g,c)$ 。当 $w(Q,A)=w(g,c)$ 时，达到最优化。

3.计算最优剪接比对

跟很多生物学序列比较一样，剪接比对也可以通过动态规划DP解决。我们计算了两个 $(m+1)\times(n+1)$ 的矩阵 $E$ 和 $I$ 。其中， $E_i^j$ 是 $g[1..j]$ 和 $c[1...i]$ 的所有比对结果中以外显子结束的最大权重， $I_i^j$ 是 $g[1..j]$ 和 $c[1...i]$ 的所有比对结果中以内含子结束的最大权重。显然， $w(g,c)$ 是 $max(E_m^n,I_m^n)$

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

4.内含子切除技术

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

4.1. 计算match

计算match时，考虑 $g$ 和 $c$ 的最大近似匹配。其实最大近似匹配就是 $g$ 和 $c$ 的字串 $g[j...r]$ 和 $c[i...h]$ 的配对。这个配对满足最左和最右原则：最左即
$j=1$ 或 $i=1$ 或 $g[j-1]\neq c[i-1]$ ；最右即 $r=n$ 或 $h=m$ 或 $g[r+1]\neq c[h+1]$ 。我们只对满足最小长度 $l_{min}$ 、最大距离 $d_{max}$ 的最大近似匹配感兴趣。最小长度指两个子串都不得短于这个长度；最大距离指两个子串的edit distance不得大于这个距离。关于edit distance，前面已经提到过了。
通常，采用seed-and-extend方法计算近似最大匹配。因为每个最大近似匹配里肯定会至少包括一个长 $\lfloor {l_{min}\over{d_{max}+1}}\rfloor$ 的最大精准匹配。可以将这个精准匹配作为种子，将其向 $g$ 和 $c$ 两边延伸得到最大近似匹配。

4.2. 将match串起来

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

4.3.切割

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

5.计算保守剪接比对

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

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

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