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

作者: 生信了 | 来源:发表于2019-11-14 11:24 被阅读0次

序列比对（十九）——基序发现和中间字符串问题
序列比对（二十）——基序发现问题的算法及实现代码
序列比对（二十一）——中间字符串问题的算法及实现代码
10.13 GATK-1局部重比对 BQSR
序列对位排列（sequence alignment）
科普笔记03：split gene 断裂基因
序列比对（二十四）——最长公共子序列
3️⃣ 多序列比对(1):简介
STAR序列比对（测试支原体污染序列）
mRNA-seq学习（二）：Bowtie2比对

原创：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$ 是常数。这样，我们就可以看出基序发现问题和中间字符串问题在求解上其实是一回事。

小结

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

（公众号：生信了）

序列比对（十九）——基序发现和中间字符串问题
原创：hxj7 本文介绍了基序发现问题和中间字符串问题。引言：DNA调控元件我们知道，DNA调控元件往往是一段...
序列比对（二十）——基序发现问题的算法及实现代码
原创：hxj7 前文介绍了基序发现问题和中间字符串问题，本文给出了基序发现问题的具体算法和实现代码。基序发现问题...
序列比对（二十一）——中间字符串问题的算法及实现代码
原创：hxj7 前文介绍了基序发现问题和中间字符串问题。本文给出了中间字符串的算法和实现代码。中间字符串问题的简...
10.13 GATK-1局部重比对 BQSR
局部重比对局部重比对的目的是将BWA比对过程中所发现有潜在序列插入或者序列删除（insertion和deleti...
序列对位排列（sequence alignment）
1. 什么是序列比对？序列比对的目的、意义和用途是什么？ 1.1：序列比对（sequence alignment）...
科普笔记03：split gene 断裂基因
真核细胞内的结构基因并非全部由编码序列组成，而是在编码序列中间插入了无编码作用的碱基序列，这类基因称为断裂基因。 ...
序列比对（二十四）——最长公共子序列
原创：hxj7 本文介绍如何求解两个字符串的最长公共子序列。最长公共子序列问题前文《序列比对（23）最长公共子...
3️⃣ 多序列比对(1):简介
序列比对和序列特征分析总目录关于多序列比对可以参考Wiki 定义：多序列比对是对3条以上（包括3条）DNA,R...
STAR序列比对（测试支原体污染序列）
转录组测序发现唯一比对率特别低，进一步将未比对到参考基因组的序列提出来比对到NCBI NT库。发现多数序列比对到了...
mRNA-seq学习（二）：Bowtie2比对
1. 比对之前需要考虑哪些问题 1. 选什么作为参考序列基因组序列既能做表达定量，还能发现新的基因和转录本转录...

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

引言：DNA调控元件

一、基序发现问题

二、中间字符串问题

三、两个问题是等价的

小结

相关文章

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

序列比对（二十）——基序发现问题的算法及实现代码

序列比对（二十一）——中间字符串问题的算法及实现代码

10.13 GATK-1局部重比对 BQSR

序列对位排列（sequence alignment）

科普笔记03：split gene 断裂基因

序列比对（二十四）——最长公共子序列

3️⃣ 多序列比对(1):简介

STAR序列比对（测试支原体污染序列）

mRNA-seq学习（二）：Bowtie2比对

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读