­

序列比对（25）编辑距离

  • 2019 年 10 月 4 日
  • 笔记

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

编辑距离

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

效果如下：

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

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

解编辑距离的代码

#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 != '') {          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("nn");          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: %dn", 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);  }

VirMach 便宜 VPS

QNews

QNews