原創(chuàng)：hxj7
前言
序列比對是生信領(lǐng)域的一個(gè)古老課題湾揽，在這一波NGS的浪潮中重新引起大家的廣泛關(guān)注。由于生物序列的特殊性售貌，在比對的時(shí)候允許插入缺失育特，所以往往是一種不精確匹配生逸。從本文開始，我們陸續(xù)學(xué)習(xí)前人開發(fā)出的流行算法。
全局比對算法
所謂全局比對算法槽袄，就是根據(jù)一個(gè)打分矩陣（替換矩陣）計(jì)算出兩個(gè)序列比對最高得分的算法。關(guān)于它的介紹網(wǎng)上已經(jīng)非常多了锋谐，我們只需看看其中的關(guān)鍵點(diǎn)及實(shí)現(xiàn)代碼遍尺。
關(guān)鍵點(diǎn)
1. 打分矩陣：
選用不同的打分矩陣或者罰分分值會(huì)導(dǎo)致比對結(jié)果不同，常用BLAST打分矩陣涮拗。
2. 計(jì)算比對最高得分的算法：
常用動(dòng)態(tài)規(guī)劃算法（Needleman-Wunsch算法）乾戏。
圖片引自http://www.reibang.com/p/2b99d0d224a2
3. 打印出最高得分相應(yīng)的序列比對結(jié)果：
根據(jù)得分矩陣回溯，如果最優(yōu)比對結(jié)果有多個(gè)三热，全部打印出來鼓择。
4. 理解打分系統(tǒng)背后的概率論模型：
比對分值可以理解為匹配模型和隨機(jī)模型的對數(shù)幾率比（log-odds ratio）。
實(shí)現(xiàn)代碼（C代碼及示例）
網(wǎng)上的教程給出的代碼大多不全就漾，所以我決定還是自己造輪子了呐能。
需要說明的是：
1. 沒有對輸入進(jìn)行檢查，如果有非AGCT的字符程序可能會(huì)出錯(cuò)抑堡。
2. 對空位的罰分是線性的摆出，即penalty＝gd，其中g(shù)為空位長度首妖，d為一個(gè)空位的罰分偎漫。
在以后的文章中，我們會(huì)給出仿射型罰分的代碼有缆，即
penalty=d+(g-1)e象踊，其中g(shù)為連續(xù)空位的長度，d為連續(xù)空位中第一個(gè)空位的罰分棚壁，e為連續(xù)空位中第2個(gè)及以后空位的罰分杯矩。
3. 其他未提及的錯(cuò)誤或者不足。
示例
代碼
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSEQ 1000
#define GAP_CHAR '-'
// 對空位的罰分是線性的
struct Unit {
int W1; // 是否往上回溯一格
int W2; // 是否往左上回溯一格
int W3; // 是否往左回溯一格
float M; // 得分矩陣第(i, j)這個(gè)單元的分值灌曙，即序列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
// 數(shù)組下標(biāo)是兩個(gè)字符的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)才算結(jié)束，這代表初始的兩個(gè)字符串的每個(gè)字符都被比對到了
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);
// 動(dòng)態(tài)規(guī)劃算法計(jì)算得分矩陣每個(gè)單元的分值
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);
// 打印最優(yōu)比對結(jié)果在刺，如果有多個(gè)逆害，全部打印
// 遞歸法
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);
// 釋放內(nèi)存
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);
}
（公眾號：生信了）
