序列比对（19）基序发现和中间字符串问题

本文介绍了基序发现问题和中间字符串问题。

引言：DNA调控元件

我们知道，DNA调控元件往往是一段相似的DNA序列。理想情况下这些序列完全一致，比如下面这样：

图片引自《生物信息学算法导论》

但实际上，这些序列不会完全一样，总会有若干位点发生“变异”，从而不同，比如下面这样：

图片引自《生物信息学算法导论》

如果给定一组DNA序列（暂且假定它们长度相等），那么如何找出这些相似的序列呢？由此可以引出两个问题，即基序发现问题和中间字符串问题。

一、基序发现问题

要说明基序是什么，首先介绍一下序列剖面（Profile）。

假设给定$t$条长度为$n$的$DNA$序列，我们为其中每条序列选择一个起点$s_i(i=1,2,...,t)$，截取该序列中以$s_i$为起点的长度为$l$的一个序列，称为一条$l$-元组序列。那么这$t$个$l$-元组序列就构成了一个$t\times l$的联配矩阵（alignment matrix），统计该矩阵的每一列中各个碱基出现次数，则构成了一个新的$4\times l$的矩阵，称为剖面矩阵（profile matrix）。剖面矩阵各列最大值对应的碱基放到一起构成了一条长度为$l$的序列，称为共有序列（consensus）。

比如下图中$t=7,l=8$，序列剖面（Profile）就是一个$4\times 8$的一个矩阵：

图片引自《生物信息学算法导论》

接下来我们给出一系列符号定义，以便下文的讨论：

我们将这$t$条长度为$n$的序列记为$DNA$；
定义 s ⃗ = { s 1 , s 2 , . . . , s t } ,   1 ≤ s i ≤ n − l + 1 ,   f o r   ( i = 1 , 2 , . . . , t ) \boldsymbol{\vec{s}} = \{s_1,s_2,...,s_t\}, \ 1 \leq s_i \leq n-l+1, \ for \ (i=1,2,...,t) s ={s1​,s2​,...,st​}, 1≤si​≤n−l+1, for (i=1,2,...,t)为起始位点向量；
定义$\mathbf{P}(\vec{\mathbf{s}})$为$t\times l$的剖面矩阵；
将剖面矩阵$\mathbf{P}(\vec{\mathbf{s}})$中第$j$列的最大值记为$M_{\mathbf{P}(\vec{\mathbf{s}})}(j),for(j=1,2,...,l)$；
将共有序列的得分记为$Score(\vec{\mathbf{s}},DNA)=\sum _{j=1}^l{M}_{\mathbf{P}(\vec{\mathbf{s}})}(j)$。

那么，基序发现问题用我们上面的符号表示就是要找到：$$\underset{\vec{\mathbf{s}}}{\mathrm{argmax}}Score(\vec{\mathbf{s}},DNA)\text{\hspace{1em}(1.1)}$$
也就是说我们要计算得到:
$$\underset{\vec{\mathbf{s}}}{\max }Score(\vec{\mathbf{s}},DNA)\text{\hspace{1em}(1.2)}$$

二、中间字符串问题

同样地，要讲清楚中间字符串问题，我们首先给出一些符号：

将一个$l$-元组序列$v$和一个以$s_i$为起始位点的$l$-元组序列的汉明距离记为$d_H(v,s_i)$，表示这两个序列中不相同位点的数目；
将一个$l$-元组序列$v$以及一组分别以 s ⃗ = { s 1 , s 2 , . . . , s t } \boldsymbol{\vec{s}} = \{s_1,s_2,...,s_t\} s ={s1​,s2​,...,st​}作为起始位点的$l$-元组序列的总汉明距离表示为$d_H(v,\vec{\mathbf{s}})=\sum _{i=1}^t{d}_H(v,s_i)$；
将一个$l$-元组序列与一组$DNA$序列的任意起始位点的总汉明距离的最小值记为$TotalDistance(v,DNA)=\underset{\vec{\mathbf{s}}}{\min }{d}_H(v,\vec{\mathbf{s}})$。

那么，中间字符串用上述符号表示就是要找到：
$$\underset{v}{\mathrm{argmin}}\underset{\vec{\mathbf{s}}}{\min }{d}_H(v,\vec{\mathbf{s}})\text{\hspace{1em}(2.1)}$$
也就是说我们要计算得到:
$$\underset{v}{\min }\underset{\vec{\mathbf{s}}}{\min }{d}_H(v,\vec{\mathbf{s}})\text{\hspace{1em}(2.2)}$$

三、两个问题是等价的

我们可以证明计算式子（1.2）和计算（2.2）是一回事。
首先，根据第一部分的定义，式（1.2）其实就是：
$$\underset{\vec{\mathbf{s}}}{\max }Score(\vec{\mathbf{s}},DNA)=\underset{\vec{\mathbf{s}}}{\max }\sum _{j=1}^l{M}_{\mathbf{P}(\vec{\mathbf{s}})}(j)\text{\hspace{1em}(3.1)}$$

而式（2.2）也可以做变换。我们再给出一些符号，假定：

$l$-元组序列 v = { v 1 , v 2 , . . . , v l } v = \{v_1,v_2,...,v_l\} v={v1​,v2​,...,vl​}
第一部分涉及到的$t\times l$阶的联配矩阵为$A_{\vec{\mathbf{s}}}^{ij}$, 其中$i=1,2,...,t;j=1,2,...,l。$
定义$S(x,y)$来判断碱基$x$和碱基$y$是否相同，即：
$$S(x,y)=\{\begin{array}{ll}1,& \text{if x=y}\\ 0,& \text{if x̸=y}\end{array}$$
定义$D(x,y)$来判断碱基$x$和碱基$y$是否不同，即：
$$D(x,y)=\{\begin{array}{ll}0,& \text{if x=y}\\ 1,& \text{if x̸=y}\end{array}$$

那么：
$$\begin{array}{rl}\underset{v}{\min }\underset{\vec{\mathbf{s}}}{\min }{d}_H(v,\vec{\mathbf{s}})& =\underset{\vec{\mathbf{s}}}{\min }\underset{v}{\min }{d}_H(v,\vec{\mathbf{s}})\\ & =\underset{\vec{\mathbf{s}}}{\min }\underset{v}{\min }\sum _{i=1}^t{d}_H(v,s_i)\\ & =\underset{\vec{\mathbf{s}}}{\min }\underset{v}{\min }\sum _{i=1}^t\sum _{j=1}^{l}D(A_{\vec{\mathbf{s}}}^{ij},v_j)\\ & =\underset{\vec{\mathbf{s}}}{\min }\underset{v}{\min }\sum _{j=1}^l\sum _{i=1}^{t}D(A_{\vec{\mathbf{s}}}^{ij},v_j)\\ & =\underset{\vec{\mathbf{s}}}{\min }\sum _{j=1}^l\underset{v_j}{\min }\sum _{i=1}^{t}D(A_{\vec{\mathbf{s}}}^{ij},v_j)\\ & =\underset{\vec{\mathbf{s}}}{\min }\sum _{j=1}^l\underset{v_j}{\min }[t-\sum _{i=1}^{t}S(A_{\vec{\mathbf{s}}}^{ij},v_j)]\\ & =\underset{\vec{\mathbf{s}}}{\min }\sum _{j=1}^l[t-\underset{v_j}{\max }\sum _{i=1}^{t}S(A_{\vec{\mathbf{s}}}^{ij},v_j)]\\ & =\underset{\vec{\mathbf{s}}}{\min }\sum _{j=1}^l[t-M_{\mathbf{P}(\vec{\mathbf{s}})}(j)]\\ & =lt-\underset{\vec{\mathbf{s}}}{\max }\sum _{j=1}^l{M}_{\mathbf{P}(\vec{\mathbf{s}})}(j)\\ & =lt-\underset{\vec{\mathbf{s}}}{\max }Score(\vec{\mathbf{s}},DNA)\end{array}$$

上式中$lt$是常数。这样，我们就可以看出基序发现问题和中间字符串问题在求解上其实是一回事。

小结

本文内容基于《生物信息学算法导论》，笔者所作的工作就是将算法推导过程补充详细。至于实现代码，我们会在后续文章中讨论。

（公众号：生信了）