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

2019-11-14  本文已影响0人  生信了

原创：hxj7

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

引言：DNA调控元件

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

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

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

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

如果给定一组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的一个矩阵：

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

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

我们将这t条长度为n的序列记为DNA；
定义\boldsymbol{\vec{s}} = \{s_1,s_2,...,s_t\}, \ 1 \leq s_i \leq n-l+1, \ for \ (i=1,2,...,t)为起始位点向量；
定义\boldsymbol{P}(\boldsymbol{\vec{s}})为t \times l的剖面矩阵；
将剖面矩阵\boldsymbol{P}(\boldsymbol{\vec{s}})中第j列的最大值记为M_{\boldsymbol{P}(\boldsymbol{\vec{s}})}(j), \ for \ (j=1,2,...,l)；
将共有序列的得分记为Score(\boldsymbol{\vec{s}}, DNA) = \displaystyle \sum_{j=1}^l M_{\boldsymbol{P}(\boldsymbol{\vec{s}})}(j)。

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

二、中间字符串问题

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

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

那么，中间字符串用上述符号表示就是要找到：
\underset{v}{\mathrm{argmin}} \, \underset{\boldsymbol{\vec{s}}}{\min} \, d_H(v,\boldsymbol{\vec{s}}) \tag{2.1}
也就是说我们要计算得到:
\underset{v}{\min} \, \underset{\boldsymbol{\vec{s}}}{\min} \, d_H(v,\boldsymbol{\vec{s}}) \tag{2.2}

三、两个问题是等价的

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

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

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

那么：
\begin{aligned} \underset{v}{\min} \, \underset{\boldsymbol{\vec{s}}}{\min} \, d_H(v,\boldsymbol{\vec{s}}) & = \underset{\boldsymbol{\vec{s}}}{\min} \, \underset{v}{\min} \, d_H(v,\boldsymbol{\vec{s}}) \\ & = \underset{\boldsymbol{\vec{s}}}{\min} \, \underset{v}{\min} \displaystyle \sum_{i=1}^t d_H(v, s_i) \\ & = \underset{\boldsymbol{\vec{s}}}{\min} \, \underset{v}{\min} \displaystyle \sum_{i=1}^t \sum_{j=1}^l D(A_{\boldsymbol{\vec{s}}}^{ij}, v_j) \\ & = \underset{\boldsymbol{\vec{s}}}{\min} \, \underset{v}{\min} \displaystyle \sum_{j=1}^l \sum_{i=1}^t D(A_{\boldsymbol{\vec{s}}}^{ij}, v_j) \\ & = \underset{\boldsymbol{\vec{s}}}{\min} \displaystyle \sum_{j=1}^l \underset{v_j}{\min} \sum_{i=1}^t D(A_{\boldsymbol{\vec{s}}}^{ij}, v_j) \\ & = \underset{\boldsymbol{\vec{s}}}{\min} \displaystyle \sum_{j=1}^l \underset{v_j}{\min} \, \Bigg[t - \sum_{i=1}^t S(A_{\boldsymbol{\vec{s}}}^{ij}, v_j) \Bigg] \\ & = \underset{\boldsymbol{\vec{s}}}{\min} \displaystyle \sum_{j=1}^l \Bigg[t - \underset{v_j}{\max} \sum_{i=1}^t S(A_{\boldsymbol{\vec{s}}}^{ij}, v_j) \Bigg] \\ & = \underset{\boldsymbol{\vec{s}}}{\min} \displaystyle \sum_{j=1}^l \Bigg[t - M_{\boldsymbol{P}(\boldsymbol{\vec{s}})}(j) \Bigg] \\ & = lt - \underset{\boldsymbol{\vec{s}}}{\max} \displaystyle \sum_{j=1}^lM_{\boldsymbol{P}(\boldsymbol{\vec{s}})}(j) \\ & = lt - \underset{\boldsymbol{\vec{s}}}{\max} \, Score(\boldsymbol{\vec{s}},DNA) \end{aligned}

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

小结

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

（公众号：生信了）

上一篇 下一篇

猜你喜欢

热点阅读