【练习】双序列比对

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/* 打分矩阵
 * 使用两个碱基的ASSIC码差的绝对值作为索引
 */
int point_matrix[20] = {
    5,  /* same */ 0, -3,          /* A/C */ 0, -3, /*  C/G*/
    0,  -3,           /* A/G */ 0, 0,           0,  0,           0, 0,
    -3, /* G/T */ 0,  0,           0,           -3, /* C/T */ 0, -3 /* A/T */
};
/**
 * @brief 打分函数
 *
 * @param base1
 * @param base2
 * @param point_matrix 打分矩阵
 * @return int
 */
int point(char base1, char base2) {
  int s = base1 - base2;
  if (s < 0) {
    s = -s;  // abs
  }
  return *(point_matrix + s);
}
// gap
/* TODO:gap extend 还不知道该怎么用 */
struct {
  int GO; /* gap open */
  int GE; /* gap extend */
} gap = {-2, -1};
/**
 * @brief 释放内存
 *
 * @param matrix
 * @param len1
 */
void destory(int** matrix, int len1) {
  for (int i = 0; i < len1; i++) {
    free(matrix[i]);
  }
  free(matrix);
}
/**
 * @brief 初始化得分矩阵
 *
 * @param len1
 * @param len2
 * @return int**
 */
int** init_matrix(int len1, int len2) {
  len1 += 1;  // 第一列全空
  len2 += 1;  // 第一行全空
  int** _matrix = NULL;
  _matrix = (int**)malloc(len1 * sizeof(int*));
  if (_matrix == NULL) return NULL;
  int pos;  // 记录当前是第几行
  int* row;
  for (int i = 0; i < len1; i++) {
    pos = i;
    row = (int*)calloc(len2, sizeof(int));
    if (row == NULL) goto clean;
    _matrix[i] = row;
  }
  for (int i = 0; i < len1; i++) _matrix[i][0] = gap.GO * (i + 1);
  for (int i = 0; i < len2; i++) _matrix[0][i] = gap.GO * (i + 1);
  return _matrix;
clean:
  destory(_matrix, pos);
  return NULL;
}
/**
 * @brief 获取三个数中最大的那一个数，同时把路径记录下来
 *
 * @param s1
 * @param s2
 * @param s3
 * @param 
 * @return int
 */
int max3(int s1, int s2, int s3) {
  return s1 > s2 ? s1 > s3 ? s1 : s3 : s2 > s3 ? s2 : s3;
}
/**
 * @brief 计算得分矩阵
 *
 * @param score_matrix
 * @param seq1
 * @param len1
 * @param seq2
 * @param len2
 */
void computer_matrix(int** score_matrix, char* seq1, int len1, char* seq2,
                     int len2) {
  for (int i = 1; i <= len1; i++) {
    for (int j = 1; j <= len2; j++) {
      int s1 = score_matrix[i - 1][j - 1] + point(seq1[i - 1], seq2[j - 1]);
      int s2 = score_matrix[i][j - 1] + gap.GO;
      int s3 = score_matrix[i - 1][j] + gap.GO;
      int max = max3(s1, s2, s3);
      score_matrix[i][j] = max;
    }
  }
}
/**
 * @brief 打印结果
 *
 * @param score_matrix
 * @param len1
 * @param len2
 */
void print_matrix(int** const score_matrix, int len1, int len2) {
  for (int x = 0; x <= len2; x++) {
    printf("%-5d", x + 1);
  }
  printf("\n");
  for (int x = 0; x <= len2; x++) {
    printf("%-5c", '|');
  }
  printf("\n");
  for (int i = 0; i <= len1; i++) {
    for (int j = 0; j <= len2; j++) {
      printf("%-5d", score_matrix[i][j]);
    }
    printf("--%-5d", i + 1);
    printf("\n");
  }
}
/**
 * @brief 开始回溯比对路径
 *
 * @param score_matrix
 * @param seq1
 * @param len1
 * @param seq2
 * @param len2
 */
void traceback(int** score_matrix, char* seq1, int len1, char* seq2, int len2) {
  int len_max = len1 + len2;
  char tseq[len_max];    // top seq
  char center[len_max];  // center char
  char bseq[len_max];    // bottom seq
  int x = len1, y = len2, count = 0;
  while (x > 0 && y > 0) {
    if (score_matrix[x][y] ==
        score_matrix[x - 1][y - 1] + point(seq1[x - 1], seq2[y - 1])) {
      tseq[count] = seq1[x - 1];
      bseq[count] = seq2[y - 1];
      center[count] = seq1[x-1] == seq2[y-1] ? '|' : '*';
      x--;
      y--;
    } else if (score_matrix[x][y] == score_matrix[x - 1][y] + gap.GO) {
      tseq[count] = seq1[--x];
      bseq[count] = center[count] = '-';
    } else {
      bseq[count] = seq2[--y];
      tseq[count] = center[count] = '-';
    }
    count++;
  }
  for (x = count - 1; x >= 0; x--) {
    printf("%c", tseq[x]);
  }
  printf("\n");
  for (x = count - 1; x >= 0; x--) {
    printf("%c", center[x]);
  }
  printf("\n");
  for (x = count - 1; x >= 0; x--) {
    printf("%c", bseq[x]);
  }
  printf("\n");
}
int main() {
  char seq1[] = "CGTTGCTTGG";
  char seq2[] = "CGTTGCTCTTTG";
  int len1 = strlen(seq1);
  int len2 = strlen(seq2);
  int** score_matrix = init_matrix(len1, len2);
  if (score_matrix == NULL) {
    printf("fail!");
    return 0;
  }
  computer_matrix(score_matrix, seq1, len1, seq2, len2);
  print_matrix(score_matrix, len1, len2);
  traceback(score_matrix, seq1, len1, seq2, len2);
  destory(score_matrix, len1);
  score_matrix = NULL;
  return 0;
}