前文介绍了重复匹配问题的动态规划算法,但是遗留了重复结果输出的问题。本文对该问题进行了补充说明。
(公众号:生信了)
前文《序列匹配(五)——重复匹配问题的动态规划算法》介绍了重复匹配问题的动态规划算法。重复匹配问题就是从序列
x
x
x中找到序列
y
y
y的多个拷贝,并且使这多个拷贝的“标准比对分值”之和最大。《生物序列分析》一书中给出的迭代公式是:
F
(
i
,
0
)
=
max
{
F
(
i
−
1
,
0
)
,
F
(
i
−
1
,
j
)
−
T
,
j
=
1
,
2
,
.
.
.
,
n
F
(
i
,
j
)
=
max
{
F
(
i
,
0
)
F
(
i
−
1
,
j
−
1
)
+
s
(
i
,
j
)
F
(
i
−
1
,
j
)
+
d
F
(
i
,
j
−
1
)
+
d
\begin{aligned} & F(i,0) = \max \begin{cases} F(i-1,0), \\ F(i-1,j)-T, \ \ \ \ \ j=1,2,...,n \end{cases} \\ & F(i,j) = \max \begin{cases} F(i,0) \\ F(i-1,j-1) + s(i,j) \\ F(i-1,j) + d \\ F(i, j-1) + d \end{cases} \end{aligned}
F(i,0)=max{F(i−1,0),F(i−1,j)−T, j=1,2,...,nF(i,j)=max⎩⎪⎪⎪⎨⎪⎪⎪⎧F(i,0)F(i−1,j−1)+s(i,j)F(i−1,j)+dF(i,j−1)+d
经过笔者分析,在前文中给出的迭代公式是:
F
(
i
,
0
)
=
max
{
F
(
i
−
1
,
0
)
,
F
(
i
−
1
,
j
)
−
T
,
j
=
1
,
2
,
.
.
.
,
n
F
(
i
,
j
)
=
max
{
F
(
i
,
0
)
F
(
i
−
1
,
j
−
1
)
+
s
(
i
,
j
)
F
(
i
−
1
,
j
)
+
d
F
(
i
,
j
−
1
)
+
d
F
(
i
,
0
)
+
s
(
i
,
j
)
,
when
j
>
1
F
(
i
,
j
)
=
max
{
F
(
i
,
0
)
F
(
i
−
1
,
j
)
+
d
F
(
i
,
0
)
+
s
(
i
,
j
)
,
when
j
=
1
\begin{aligned} & F(i,0) = \max \begin{cases} F(i-1,0), \\ F(i-1,j)-T, \ \ \ \ \ j=1,2,...,n \end{cases} \\ & F(i,j) = \max \begin{cases} F(i,0) \\ F(i-1,j-1) + s(i,j) \\ F(i-1,j) + d \\ F(i, j-1) + d \\ F(i,0) + s(i,j) \end{cases}\ \ \ , \ \ \ \text{when $j > 1$} \\ & F(i,j) = \max \begin{cases} F(i,0) \\ F(i-1,j) + d \\ F(i,0) + s(i,j) \end{cases}\ \ \ , \ \ \ \text{when $j = 1$} \end{aligned}
F(i,0)=max{F(i−1,0),F(i−1,j)−T, j=1,2,...,nF(i,j)=max⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧F(i,0)F(i−1,j−1)+s(i,j)F(i−1,j)+dF(i,j−1)+dF(i,0)+s(i,j) , when j>1F(i,j)=max⎩⎪⎨⎪⎧F(i,0)F(i−1,j)+dF(i,0)+s(i,j) , when j=1
但是这个公式在回溯时会出现重复结果输出的问题,比如:
经过笔者近期回顾,发现重复结果输出的可能途径有二(当然可能还有其他途径):
一是:
F
(
i
,
0
)
→
F
(
i
−
1
,
j
)
−
T
→
F
(
i
−
1
,
0
)
F
(
i
,
0
)
⟶
F
(
i
−
1
,
0
)
\begin{aligned} & F(i,0) \rightarrow F(i-1,j)-T \rightarrow F(i-1,0) \\ & F(i,0) \ \ \ \ \ \ \ \ \ \ \ \ \ \longrightarrow \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \, F(i-1,0) \end{aligned}
F(i,0)→F(i−1,j)−T→F(i−1,0)F(i,0) ⟶ F(i−1,0)
这两种回溯步骤结果一致,但意义互相矛盾。二者只能取其一。
二是:
F
(
i
,
j
)
→
F
(
i
−
1
,
j
−
1
)
+
s
(
i
,
j
)
→
F
(
i
−
1
,
0
)
F(i,j) \rightarrow F(i-1,j-1)+s(i,j) \rightarrow F(i-1,0)
F(i,j)→F(i−1,j−1)+s(i,j)→F(i−1,0)
第一个箭头表示的回溯步骤要求 x i − 1 x_{i-1} xi−1参与联配,而第二个箭头表示的回溯步骤要求 x i − 1 x_{i-1} xi−1不参与联配。二者互相矛盾,所以该回溯不成立。
校正公式和代码
笔者认为,上述重复结果输出的问题主要可能是因为
F
(
i
,
j
)
F(i,j)
F(i,j)的定义不清晰,导致
F
(
i
,
j
)
F(i,j)
F(i,j)既能表示
x
i
x_i
xi参与联配,也能表示
x
i
x_i
xi不参与联配。
为此,笔者给出新的更清晰的符号定义以及公式。
F
(
i
)
F(i)
F(i)表示以
x
i
x_i
xi结尾且
x
i
x_i
xi不参与联配的最高得分。
A
(
i
,
j
)
A(i,j)
A(i,j)表示以
x
i
x_i
xi和
y
j
y_j
yj结尾且
x
i
x_i
xi参与联配的最高得分。
那么迭代公式是:
F
(
0
)
=
0
,
F
(
i
)
=
max
{
F
(
i
−
1
)
,
A
(
i
−
1
,
j
)
−
T
,
j
=
1
,
2
,
.
.
.
,
n
A
(
0
,
0
)
=
d
,
A
(
0
,
j
)
=
d
,
j
=
1
,
2
,
.
.
.
,
n
A
(
i
,
0
)
=
F
(
i
)
+
d
,
A
(
i
,
j
)
=
max
{
F
(
i
)
+
s
(
i
,
j
)
,
A
(
i
−
1
,
j
−
1
)
+
s
(
i
,
j
)
,
A
(
i
−
1
,
j
)
+
d
,
A
(
i
,
j
−
1
)
+
d
,
\begin{aligned} & F(0) = 0, \\ & F(i) = \max \begin{cases} F(i-1), \\ A(i-1,j)-T, \ \ \ \ \ j=1,2,...,n \end{cases} \\ & A(0,0) = d, \\ & A(0,j) = d, \ \ \ \ \ j =1,2,...,n \\ & A(i,0) = F(i) + d, \\ & A(i,j) = \max \begin{cases} F(i) + s(i,j), \\ A(i-1,j-1) + s(i,j), \\ A(i-1,j) + d, \\ A(i,j-1) + d, \end{cases} \end{aligned}
F(0)=0,F(i)=max{F(i−1),A(i−1,j)−T, j=1,2,...,nA(0,0)=d,A(0,j)=d, j=1,2,...,nA(i,0)=F(i)+d,A(i,j)=max⎩⎪⎪⎪⎨⎪⎪⎪⎧F(i)+s(i,j),A(i−1,j−1)+s(i,j),A(i−1,j)+d,A(i,j−1)+d,
需要说明的是:
- 上述公式中令 A ( 0 , 0 ) = d A(0,0)=d A(0,0)=d是为了让 A ( 0 , 0 ) ≠ F ( 0 ) A(0,0) \neq F(0) A(0,0)=F(0),以免产生重复结果输出。
- A ( i , 0 ) A(i,0) A(i,0)不光意味着 x i x_i xi和一个空位联配,且代表了一个“旧联配段”的结束和一个“新联配段”的开始,只不过这个新联配是以空位开始的。
- A ( i , j ) A(i,j) A(i,j)表示 x i x_i xi参与联配。其迭代公式中第一项表示 x i x_i xi开始了一个新的局部联配;而后三项表示 x i x_i xi加入(或者说延续)前一个局部联配。
- 最终比对的最高得分就是 F ( m + 1 ) F(m+1) F(m+1)。
这样的公式目前还没有出现重复结果输出的问题:
相应的代码放在了文末。
对比对总长度的估计
在写代码的过程中,遇到一个小问题,就是要对最优比对结果的总长度
(以下简称比对总长度)进行估计。即
x
x
x序列长度为
m
m
m,
y
y
y序列长度为
n
n
n。那么比对总长度
l
l
l是在什么范围内呢?这里说的总长度包括了
x
x
x序列中未参与联配的碱基。
由于
x
x
x序列每个碱基在迭代的过程中都要被遍历到,所以
l
l
l的下界
是
m
m
m(比如考虑
x
x
x序列是
G
A
C
G
A
C
T
GACGACT
GACGACT,而
y
y
y序列是
G
A
C
T
GACT
GACT,那么
l
=
m
=
7
l=m=7
l=m=7):
l
≥
m
l \ge m
l≥m
假设比对时
m
a
t
c
h
match
match的得分是
M
M
M,而
g
a
p
gap
gap的得分是
G
G
G。那么,总长度的一个宽松上界
是:
l
≤
m
+
[
m
∣
M
G
∣
]
l \le m + \left[ m \left| \frac{M}{G} \right| \right]
l≤m+[m∣∣∣∣GM∣∣∣∣]
其中, ∣ a ∣ \left| a \right| ∣a∣表示 a a a的绝对值, [ a ] \left[ a \right] [a]表示不大于 a a a的最大整数。上述不等式证明如下:
假设将长度为 m m m的 x x x序列分割成一系列子序列 S i S_i Si,这些子序列构成一个
完整分割集
S S S,即
⋃ i S i = S and ⋂ i S i = ∅ \bigcup_{i} S_i = S \ \ \ \text{and} \ \ \ \bigcap_{i} S_i = \emptyset i⋃Si=S and i⋂Si=∅
如果将 S i S_i Si的长度记为 m i m_i mi,则有:
0 < m i ≤ m and ∑ i m i = m 0 < m_i \le m \ \ \ \text{and} \ \ \ \sum_i m_i = m 0<mi≤m and i∑mi=m
并且, S i S_i Si要么是“非联配段”,要么是“联配段”。所谓“非联配段”,即该段序列中的所有碱基都不参与最终的联配;所谓“联配段”,即该段序列中的所有碱基都参与最终的联配。
接着,将每个子序列 S i S_i Si与长度为 n n n的序列 y y y进行局部比对,比对的长度记为 L ( S i , y ) L(S_i,y) L(Si,y)。
如果 S i S_i Si是“联配段”,由于是局部比对,假设比对时 m a t c h match match的得分是 M M M,而 g a p gap gap的得分是 G G G(这里假设都是简单情形
,即所有符号的 m a t c h match match得分都是一样的,所有符号的 g a p gap gap得分也都是一样的,且空位罚分是线性的),
则有:
L ( S i , y ) ≤ m i + [ ∣ m i M G ∣ ] L(S_i,y) \le m_i + \left[ \left| \frac{m_{i}M}{G} \right| \right] L(Si,y)≤mi+[∣∣∣∣GmiM∣∣∣∣]
其中, ∣ a ∣ \left| a \right| ∣a∣表示 a a a的绝对值, [ a ] \left[ a \right] [a]表示不大于 a a a的最大整数。上述不等式是一个宽松上界,证明也很简单,其关键就是局部比对的总分值不能小于0,以至于会限制空位的数量不得高于 [ ∣ m i M G ∣ ] \left[ \left| \frac{m_{i}M}{G} \right| \right] [∣∣GmiM∣∣]。
如果 S i S_i Si是“非联配段”,那么:
L ( S i , y ) = m i ≤ m i + [ ∣ m i M G ∣ ] \begin{aligned} L(S_i,y) & = m_i \\ & \le m_i + \left[ \left| \frac{m_{i}M}{G} \right| \right] \end{aligned} L(Si,y)=mi≤mi+[∣∣∣∣GmiM∣∣∣∣]
所以,不管 S i S_i Si是“非联配段”,还是“联配段”,都有:
L ( S i , y ) ≤ m i + [ ∣ m i M G ∣ ] L(S_i,y) \le m_i + \left[ \left| \frac{m_{i}M}{G} \right| \right] L(Si,y)≤mi+[∣∣∣∣GmiM∣∣∣∣]
那么,这些子序列的比对结果拼接起来的长度为:
∑ i L ( S i , y ) ≤ ∑ i ( m i + [ ∣ m i M G ∣ ] ) ≤ ∑ i m i + [ ∑ i ∣ m i M G ∣ ] = m + [ ∑ i m i ∣ M G ∣ ] = m + [ m ∣ M G ∣ ] \begin{aligned} \sum_i L(S_i, y) & \le \sum_i \left( m_i + \left[ \left| \frac{m_{i}M}{G} \right| \right] \right) \\ & \le \sum_i m_i + \left[ \sum_i \left| \frac{m_{i}M}{G} \right| \right] \\ & = m + \left[ \sum_i m_{i} \left| \frac{M}{G} \right| \right] \\ & = m + \left[ m \left| \frac{M}{G} \right| \right] \end{aligned} i∑L(Si,y)≤i∑(mi+[∣∣∣∣GmiM∣∣∣∣])≤i∑mi+[i∑∣∣∣∣GmiM∣∣∣∣]=m+[i∑mi∣∣∣∣GM∣∣∣∣]=m+[m∣∣∣∣GM∣∣∣∣]
我们知道,重复匹配问题其实就是将 x x x序列分割成一段段互不交叠的子序列,然后每一段子序列都和 y y y序列进行局部比对。
假设重复匹配问题最终的比对结果是将 x x x序列按上面所述完整分割集分割成一系列子序列 S ^ i \hat{S}_i S^i,那么重复匹配问题的比对总长度当然也服从上述不等式,即
l = ∑ i L ( S ^ i , y ) ≤ m + [ m ∣ M G ∣ ] \begin{aligned} l & = \sum_i L(\hat{S}_i, y) \\ & \le m + \left[ m \left| \frac{M}{G} \right| \right] \end{aligned} l=i∑L(S^i,y)≤m+[m∣∣∣∣GM∣∣∣∣]
总的来说,重复匹配问题的比对总长度(在上述简单情形下)在以下范围内:
m
≤
l
≤
m
+
[
m
∣
M
G
∣
]
m \le l \le m + \left[ m \left| \frac{M}{G} \right| \right]
m≤l≤m+[m∣∣∣∣GM∣∣∣∣]
实现代码
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define MAXSEQ 1000
#define GAP_CHAR '-'
#define UNALIGN_CHAR '.'
#define max2(a, b) ((a) > (b) ? (a) : (b))
// 对空位的罚分是线性的
struct FUnit {
int W0; // X{i-1}不参与联配
int* Wj; // 跳转到A(i - 1, j)
int nj; // Wj数组的大小
float M; // F(i,0)的值
};
typedef struct FUnit* pFUnit;
// 00000001 F(i, 0) + s(i, j)
// 00000010 是否往左回溯一格
// 00000100 是否往左上回溯一格
// 00001000 是否往上回溯一格
struct AUnit {
int Wi; // 不同的bit代表不同的回溯方式
float M;
};
typedef struct AUnit* pAUnit;
float gap = -2.5; // 对空位的罚分
float match = 5;
float unmatch = -4;
void strUpper(char *s);
float maxArray(float *a, int n);
float getFScore(char a, char b);
void printFAlign(pFUnit* f, pAUnit** a, const int i, char* s, char* r, char* saln, char* raln, int n, int flag);
void printAAlign(pFUnit* f, pAUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n);
void align(char *s, char *r, float t);
int main() {
char s[MAXSEQ];
char r[MAXSEQ];
float t;
printf("The 1st seq: ");
scanf("%s", s);
printf("The 2nd seq: ");
scanf("%s", r);
printf("T (threshold): ");
scanf("%f", &t);
align(s, r, t);
return 0;
}
void strUpper(char *s) {
while (*s != '\0') {
if (*s >= 'a' && *s <= 'z') {
*s -= 32;
}
s++;
}
}
float maxArray(float *a, int n) {
float max = a[0];
int i;
for (i = 1; i < n; i++) {
if (a[i] > max)
max = a[i];
}
return max;
}
// 替换矩阵:match分值为5,mismatch分值为-4
// 数组下标是两个字符的ascii码减去65之后的和
float FMatrix[] = {
5, 0, -4, 0, 5, 0, -4, 0, -4, 0,
0, 0, 5, 0, 0, 0, 0, 0, 0, -4,
0, -4, 0, 0, 0, -4, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 5
};
float getFScore(char a, char b) {
return FMatrix[a + b - 'A' - 'A'];
}
void printFAlign(pFUnit* f, pAUnit** a, const int i, char* s, char* r, char* saln, char* raln, int n, int flag) {
// flag: 是否是从F(i+1,0)跳转过来,而不是从F(i, 0) + s(i, j)跳转过来
int k;
pFUnit p = f[i];
//printf("i=%d, n=%d, flag=%d\n", i, n, flag);
if (! i) { // 保证序列s的每个字符都比对上
for (k = n - 1; k >= 0; k--)
printf("%c", saln[k]);
printf("\n");
for (k = n - 1; k >= 0; k--)
printf("%c", raln[k]);
printf("\n\n");
return;
}
if (flag) {
saln[n] = s[i - 1];
raln[n] = UNALIGN_CHAR;
}
if (p->W0)
printFAlign(f, a, i - 1, s, r, saln, raln, n + flag, 1);
for (k = 0; k < p->nj; k++)
if (p->Wj[k])
printAAlign(f, a, i - 1, k + 1, s, r, saln, raln, n + flag);
}
void printAAlign(pFUnit* f, pAUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n) {
int k;
pAUnit p = a[i][j];
//printf("i=%d, j=%d, n=%d\n");
if (! i) { // 保证序列s的每个字符都比对上
for (k = n - 1; k >= 0; k--)
printf("%c", saln[k]);
printf("\n");
for (k = n - 1; k >= 0; k--)
printf("%c", raln[k]);
printf("\n\n");
return;
}
if (! j) { // F(i, 0) + d
saln[n] = s[i - 1];
raln[n] = GAP_CHAR;
printFAlign(f, a, i, s, r, saln, raln, n + 1, 0);
} else {
if (p->Wi & 1) { // F(i, 0) + s(i, j)
saln[n] = s[i - 1];
raln[n] = r[j - 1];
printFAlign(f, a, i, s, r, saln, raln, n + 1, 0);
}
if (p->Wi & 2) { // 向上回溯一格
saln[n] = s[i - 1];
raln[n] = GAP_CHAR;
printAAlign(f, a, i - 1, j, s, r, saln, raln, n + 1);
}
if (p->Wi & 4) { // 向左上回溯一格
saln[n] = s[i - 1];
raln[n] = r[j - 1];
printAAlign(f, a, i - 1, j - 1, s, r, saln, raln, n + 1);
}
if (p->Wi & 8) { // 向左回溯一格
saln[n] = GAP_CHAR;
raln[n] = r[j - 1];
printAAlign(f, a, i, j - 1, s, r, saln, raln, n + 1);
}
}
}
void align(char *s, char *r, float t) {
int i, j, k;
int m = strlen(s);
int n = strlen(r);
float tm[4];
float em; // F(m + 1, 0)
pFUnit* fUnit;
pAUnit** aUnit;
char* salign;
char* ralign;
int alnSize = m + floor(m * abs(match / gap));
// 初始化
if ((fUnit = (pFUnit*) malloc(sizeof(pFUnit) * (m + 1))) == NULL || \
(aUnit = (pAUnit**) malloc(sizeof(pAUnit*) * (m + 1))) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
for (i = 0; i <= m; i++) {
if ((fUnit[i] = (pFUnit) malloc(sizeof(struct FUnit))) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
fUnit[i]->W0 = 0;
fUnit[i]->nj = n;
// 创建F(i, 0)的跳转数组
if ((fUnit[i]->Wj = (int*) malloc(sizeof(int) * fUnit[i]->nj)) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
for (k = 0; k < fUnit[i]->nj; k++) {
fUnit[i]->Wj[k] = 0;
}
}
for (i = 0; i <= m; i++) {
if ((aUnit[i] = (pAUnit *) malloc(sizeof(pAUnit) * (n + 1))) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
for (j = 0; j <= n; j++) {
if ((aUnit[i][j] = (pAUnit) malloc(sizeof(struct AUnit))) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
aUnit[i][j]->Wi = 0;
}
}
fUnit[0]->M = 0;
aUnit[0][0]->M = gap; // 注意这里设置A(0,0) != 0 是很有必要的,否则A(0,0)=F(0,0)会导致重复结果输出
for (j = 1; j <= n; j++)
aUnit[0][j]->M = gap;
// 将字符串都变成大写
strUpper(s);
strUpper(r);
// 动态规划算法计算得分矩阵每个单元的分值
for (i = 1; i <= m; i++) {
// 计算F(i, 0) i >= 1
fUnit[i]->M = fUnit[i - 1]->M;
for (j = 1; j <= n; j++)
if (fUnit[i]->M < aUnit[i - 1][j]->M - t)
fUnit[i]->M = aUnit[i - 1][j]->M - t;
if (fUnit[i]->M == fUnit[i - 1]->M)
fUnit[i]->W0 = 1;
for (j = 1; j <= n; j++)
if (fUnit[i]->M == aUnit[i - 1][j]->M - t)
fUnit[i]->Wj[j - 1] = 1;
// 计算A(i, 0) i >= 1
aUnit[i][0]->M = fUnit[i]->M + gap;
aUnit[i][0]->Wi |= 1;
// 计算A(i, j) i >= 1 and j >= 1
for (j = 1; j <= n; j++) {
tm[0] = fUnit[i]->M + getFScore(s[i - 1], r[j - 1]);
tm[1] = aUnit[i - 1][j]->M + gap;
tm[2] = aUnit[i - 1][j - 1]->M + getFScore(s[i - 1], r[j - 1]);
tm[3] = aUnit[i][j - 1]->M + gap;
aUnit[i][j]->M = maxArray(tm, 4);
for (k = 0; k < 4; k++)
if (tm[k] == aUnit[i][j]->M)
aUnit[i][j]->Wi |= 1 << k;
}
}
// 计算F(m + 1, 0)
em = fUnit[m]->M;
for (j = 1; j <= n; j++)
if (em < aUnit[m][j]->M - t)
em = aUnit[m][j]->M - t;
/*
// 打印得分矩阵
for (i = 0; i <= m; i++)
printf("%f ", fUnit[i]->M);
printf("\n");
for (i = 0; i <= m; i++) {
for (j = 0; j <= n; j++)
printf("%f ", aUnit[i][j]->M);
printf("\n");
}
*/
printf("max score: %f\n", em);
// 打印最优比对结果,如果有多个,全部打印
// 递归法
if (em == 0) {
fputs("No seq aligned.\n", stdout);
} else {
if ((salign = (char*) malloc(sizeof(char) * alnSize)) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
if ((ralign = (char*) malloc(sizeof(char) * alnSize)) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
if (em == fUnit[m]->M)
printFAlign(fUnit, aUnit, m, s, r, salign, ralign, 0, 1);
for (j = 1; j <= n; j++)
if (em == aUnit[m][j]->M - t)
printAAlign(fUnit, aUnit, m, j, s, r, salign, ralign, 0);
// 释放内存
free(salign);
free(ralign);
}
for (i = 0; i <= m; i++) {
for (j = 0; j <= n; j++)
free(aUnit[i][j]);
free(aUnit[i]);
}
free(aUnit);
for (i = 0; i <= m; i++) {
free(fUnit[i]->Wj);
free(fUnit[i]);
}
free(fUnit);
}