本文介绍两个字符串的编辑距离并给出代码。

编辑距离

所谓编辑距离，就是给定两个字符串后，将一个字符串变为另一个字符串所需要花费的最少步骤。这个改变包括“插入一个字符”、“删除一个字符”，“替换一个字符”。比如：$\bm{v}=TGCATAT$v=TGCATAT与$\bm{w}=ATCCGAT$w=ATCCGAT这两个字符串的编辑距离为4。

编辑距离的求解过程和全局比对是十分相似的（关于全局比对，可以参见前文《序列比对（一）全局比对Needleman-Wunsch算法》），都需要全部符号参与比对，都允许插入、缺失和错配。所以，编辑距离可以用动态规划算法求解，其迭代公式是：

$\begin{aligned} & \text{$F(i,j)$ is the minimum score of alignments between $x_{1 \ldots i}$ and $y_{1 \ldots j}$.} \\ & F(i, 0) = i \quad \text{for $i = 0 \ldots m$.} \\ & F(0, j) = j \quad \text{for $j = 1 \ldots n$.} \\ & s(i,j) = \begin{cases} 0 & \text{if $x_i = y_j$,} \\ 1 & \text{otherwise.} \end{cases} \\ & F(i, j) = \min \begin{cases} F(i - 1, j) + 1, \\ F(i, j - 1) + 1, \\ F(i - 1, j - 1) + s(i, j) \end{cases} \end{aligned}$​F(i,j) is the minimum score of alignments between x1…i​ and y1…j​.F(i,0)=ifor i=0…m.F(0,j)=jfor j=1…n.s(i,j)={01​if xi​=yj​,otherwise.​F(i,j)=min⎩⎪⎨⎪⎧​F(i−1,j)+1,F(i,j−1)+1,F(i−1,j−1)+s(i,j)​​

效果如下：

（公众号：生信了）

编辑距离与最长公共子序列

在只允许插入和缺失而不允许错配的情况下，两个字符串的编辑距离可以通过最长公共子序列的长度（关于最长公共子序列，可以参看前文《序列比对（24）最长公共子序列》）间接算出来。

在只允许插入和缺失而不允许错配的情况下，假设给定两个字符串$\bm{v}$v和$\bm{w}$w，长度分别为$m$m和$n$n。将最长公共子序列的长度记为$LC(\bm{v}, \bm{w})$LC(v,w)而将编辑距离记为$DE(\bm{v}, \bm{w})$DE(v,w)的话，那么：
$DE(\bm{v}, \bm{w}) = m + n - 2 \times LC(\bm{v}, \bm{w})$DE(v,w)=m+n−2×LC(v,w)

求解编辑距离的代码

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSEQ 1000
#define GAP_CHAR '-'
#define min(a, b) ((a) < (b) ? (a) : (b))
#define diff(a, b) ((a) == (b) ? 0 : 1)

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

void strUpper(char *s);
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++;
    }
}

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);
    int m1, m2, m3, minm;
    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;
        }
    }
    for (i = 0; i <= m; i++) {
        aUnit[i][0]->M = i;
        aUnit[i][0]->W1 = 1;
    }
    for (j = 1; j <= n; j++) {
        aUnit[0][j]->M = 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 + 1;
            m2 = aUnit[i - 1][j - 1]->M + diff(s[i - 1], r[j - 1]);
            m3 = aUnit[i][j - 1]->M + 1;
            minm = min(min(m1, m2), m3);
            aUnit[i][j]->M = minm;
            if (m1 == minm) aUnit[i][j]->W1 = 1;
            if (m2 == minm) aUnit[i][j]->W2 = 1;
            if (m3 == minm) aUnit[i][j]->W3 = 1;
        }
    }
/*
    // 打印得分矩阵
    for (i = 0; i <= m; i++) {
        for (j = 0; j <= n; j++)
            printf("%d ", aUnit[i][j]->M);
        printf("\n");
    }
*/
    printf("min score: %d\n", aUnit[m][n]->M);
    // 打印最优比对结果，如果有多个，全部打印
    // 递归法
    if (aUnit[m][n]->M == 0) {
        fputs("Two seqs are totally the same.\n", stdout);
    } else {
        if ((salign = (char*) malloc(sizeof(char) * (m + n))) == NULL || \
            (ralign = (char*) malloc(sizeof(char) * (m + n))) == 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);
}