原创：hxj7

前言

序列比对是生信领域的一个古老课题，在这一波NGS的浪潮中重新引起大家的广泛关注。由于生物序列的特殊性，在比对的时候允许插入缺失，所以往往是一种不精确匹配。从本文开始，我们陆续学习前人开发出的流行算法。

全局比对算法

所谓全局比对算法，就是根据一个打分矩阵(替换矩阵)计算出两个序列比对最高得分的算法。关于它的介绍网上已经非常多了，我们只需看看其中的关键点及实现代码。

关键点

1. 打分矩阵：

选用不同的打分矩阵或者罚分分值会导致比对结果不同，常用BLAST打分矩阵。

2. 计算比对最高得分的算法：

常用动态规划算法(Needleman-Wunsch算法)。

image

图片引自https://www.jianshu.com/p/2b99d0d224a2

3. 打印出最高得分相应的序列比对结果：

根据得分矩阵回溯，如果最优比对结果有多个，全部打印出来。

4. 理解打分系统背后的概率论模型：

比对分值可以理解为匹配模型和随机模型的对数几率比(log-odds ratio)。

实现代码(C代码及示例)

网上的教程给出的代码大多不全，所以我决定还是自己造轮子了。

需要说明的是：

1. 没有对输入进行检查，如果有非AGCT的字符程序可能会出错。

2. 对空位的罚分是线性的，即penalty＝gd，其中g为空位长度，d为一个空位的罚分。

在以后的文章中，我们会给出仿射型罚分的代码，即

penalty=d+(g-1)e，其中g为连续空位的长度，d为连续空位中第一个空位的罚分，e为连续空位中第2个及以后空位的罚分。

3. 其他未提及的错误或者不足。

示例

image

代码

#include

#include

#include

#define MAXSEQ 1000

#define GAP_CHAR '-'

// 对空位的罚分是线性的

struct Unit {

int W1; // 是否往上回溯一格

int W2; // 是否往左上回溯一格

int W3; // 是否往左回溯一格

float M; // 得分矩阵第(i, j)这个单元的分值，即序列s(1,...,i)与序列r(1,...,j)比对的最高得分

};

typedef struct Unit *pUnit;

void strUpper(char *s);

float max3(float a, float b, float c);

float getFScore(char a, char b);

void printAlign(pUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n);

void align(char *s, char *r);

int main() {

char s[MAXSEQ];

char r[MAXSEQ];

printf("The 1st seq: ");

scanf("%s", s);

printf("The 2nd seq: ");

scanf("%s", r);

align(s, r);

return 0;

}

void strUpper(char *s) {

while (*s != '\0') {

if (*s >= 'a' && *s <= 'z') {

*s -= 32;

}

s++;

}

}

float max3(float a, float b, float c) {

float f = a > b ? a : b;

return f > c ? f : c;

}

// 替换矩阵：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 printAlign(pUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n) {

int k;

pUnit p = a[i][j];

if (! (i || j)) { // 到矩阵单元(0, 0)才算结束，这代表初始的两个字符串的每个字符都被比对到了

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 (p->W1) { // 向上回溯一格

saln[n] = s[i - 1];

raln[n] = GAP_CHAR;

printAlign(a, i - 1, j, s, r, saln, raln, n + 1);

}

if (p->W2) { // 向左上回溯一格

saln[n] = s[i - 1];

raln[n] = r[j - 1];

printAlign(a, i - 1, j - 1, s, r, saln, raln, n + 1);

}

if (p->W3) { // 向左回溯一格

saln[n] = GAP_CHAR;

raln[n] = r[j - 1];

printAlign(a, i, j - 1, s, r, saln, raln, n + 1);

}

}

void align(char *s, char *r) {

int i, j;

int m = strlen(s);

int n = strlen(r);

float gap = -2.5; // 对空位的罚分

float m1, m2, m3, maxm;

pUnit **aUnit;

char* salign;

char* ralign;

// 初始化

if ((aUnit = (pUnit **) malloc(sizeof(pUnit*) * (m + 1))) == NULL) {

fputs("Error: Out of space!\n", stderr);

exit(1);

}

for (i = 0; i <= m; i++) {

if ((aUnit[i] = (pUnit *) malloc(sizeof(pUnit) * (n + 1))) == NULL) {

fputs("Error: Out of space!\n", stderr);

exit(1);

}

for (j = 0; j <= n; j++) {

if ((aUnit[i][j] = (pUnit) malloc(sizeof(struct Unit))) == NULL) {

fputs("Error: Out of space!\n", stderr);

exit(1);

}

aUnit[i][j]->W1 = 0;

aUnit[i][j]->W2 = 0;

aUnit[i][j]->W3 = 0;

}

}

aUnit[0][0]->M = 0;

for (i = 1; i <= m; i++) {

aUnit[i][0]->M = gap * i;

aUnit[i][0]->W1 = 1;

}

for (j = 1; j <= n; j++) {

aUnit[0][j]->M = gap * j;

aUnit[0][j]->W3 = 1;

}

// 将字符串都变成大写

strUpper(s);

strUpper(r);

// 动态规划算法计算得分矩阵每个单元的分值

for (i = 1; i <= m; i++) {

for (j = 1; j <= n; j++) {

m1 = aUnit[i - 1][j]->M + gap;

m2 = aUnit[i - 1][j - 1]->M + getFScore(s[i - 1], r[j - 1]);

m3 = aUnit[i][j - 1]->M + gap;

maxm = max3(m1, m2, m3);

aUnit[i][j]->M = maxm;

if (m1 == maxm) aUnit[i][j]->W1 = 1;

if (m2 == maxm) aUnit[i][j]->W2 = 1;

if (m3 == maxm) aUnit[i][j]->W3 = 1;

}

}

/*

// 打印得分矩阵

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", aUnit[m][n]->M);

// 打印最优比对结果，如果有多个，全部打印

// 递归法

if ((salign = (char*) malloc(sizeof(char) * (m + n + 1))) == NULL) {

fputs("Error: Out of space!\n", stderr);

exit(1);

}

if ((ralign = (char*) malloc(sizeof(char) * (m + n + 1))) == NULL) {

fputs("Error: Out of space!\n", stderr);

exit(1);

}

printAlign(aUnit, m, n, 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);

}

(公众号：生信了)