序列比对(二)——Needleman-Wunsch算法之仿射罚分

原创： hxj7

关于全局比对Needleman-Wunsch算法的介绍及线性罚分的实现代码请参考序列比对（一）全局比对Needleman-Wunsch算法。

概念引用如下：
线性罚分，即penalty＝g*d，其中g为空位长度，d为一个空位的罚分。
仿射型罚分，即penalty=d+(g-1)*e， 其中g为连续空位的长度，d为连续空位中第一个空位的罚分，e为连续空位中第2个及以后空位的罚分。

这样的话，对某个位置的空位的罚分的分值需要根据前一个位置是否也是空位来定。代码会稍稍复杂一点（有一篇网上文章提到了可以利用有限状态机来解决，但是没有给出代码）。得分矩阵计算如下：

图片引自https://fengcong97.cn/sequece-alignment/
需要特别说明的是X(0, j)以及Y(i, 0)分值的选择，至少需要保证不影响X(1, j)及Y(i, 1)分值的计算**。
这套计分系统与线性罚分系统最终效果的不同主要体现在仿射罚分会倾向于让临近的空位尽可能挨在一起形成连续空位，比如下面这个例子：

本文给出了对空位进行仿射罚分的代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSEQ 1000
#define GAP_CHAR '-'

// 对空位的罚分是仿射的
struct Unit {
   
    int W1;   // 是否往上回溯一格
    int W2;   // 是否往左上回溯一格
    int W3;   // 是否往左回溯一格
    float X;
    float Y;
    float M;
    float O;      // 得分矩阵第(i, j)这个单元的分值，即序列s(1,...,i)与序列r(1,...,j)比对的最高得分
};
typedef struct Unit *pUnit;

void strUpper(char *s);
float max2(float a, float b);
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 max2(float a, float b) {
   
    return a > b ? a : b;
}

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