序列比对(四)——Smith-Waterman算法之仿射罚分

原创: hxj7

前言： 本文介绍的是采用仿射罚分模型的Smith-Waterman算法。

关于全局联配，局部联配以及仿射罚分模型的介绍可参见前文：
序列比对（一）全局比对Needleman-Wunsch算法
序列比对（二）Needleman-Wunsch算法之仿射罚分
序列比对（三）局部联配Smith-Waterman算法
其中，《序列比对（二）Needleman-Wunsch算法之仿射罚分》一文中对**X(i, j)以及Y(i, j)**的计算需要补充说明一下（以X(i, j)为例，Y(i, j)与之类似）：

在上文中X(i, j)的计算方式如下：

图片引自https://fengcong97.cn/sequece-alignment/
其前提是假设一个缺失和一个插入不会连在一起。否则，X(i, j)的计算就会变成下面这样：

局部联配线性罚分模型与仿射罚分模型结果的区别
（alignSW代表局部联配线性罚分模型，alignSW2代表局部联配仿射罚分模型）：

C代码

#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