【练习】双序列比对

简介: 【练习】双序列比对

还剩一个TODO 需要解决。

#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;
}
相关文章
|
算法 Python
双序列比对
双序列比对
420 0
双序列比对
|
机器学习/深度学习 自然语言处理 算法
【自然语言处理】正向、逆向、双向最长匹配算法的 切分效果与速度测评
【自然语言处理】正向、逆向、双向最长匹配算法的 切分效果与速度测评
220 0
|
9月前
|
安全
基因序列比对的注意点
基因序列比对的注意点
|
算法 大数据 数据库
《大数据算法》一2.5 串相等判定算法
本节书摘来华章计算机《大数据算法》一书中的第2章 ,第2.5节,王宏志 编著, 更多章节内容可以访问云栖社区“华章计算机”公众号查看。 2.5 串相等判定算法 本节讨论一个通信亚线性算法问题,因为在很多情况下,数据传输时间和数据量大致成正比,因而将通信亚线性算法归到本章讨论。
1926 0
m 序列(最长线性反馈移位寄存器序列)详解
m 序列(最长线性反馈移位寄存器序列)详解
639 0
|
算法 数据挖掘 索引
转录组入门(5): 序列比对
欢迎来GitHub上fork,一起进步: https://github.com/xuzhougeng 比对软件很多,首先大家去收集一下,因为我们是带大家入门,请统一用hisat2,并且搞懂它的用法。
2381 0
|
存储 安全 算法
使用jotp实现双因子验证
扫盲使用totp增强身份安全性指南,原理看懂也不用自己造轮子呀,最讨厌哪些啥也不懂的搬运工,我这里给大家解惑吧
652 0
|
监控 算法
算法训练Day37|738.单调递增的数字 ● 968.监控二叉树
算法训练Day37|738.单调递增的数字 ● 968.监控二叉树
|
9月前
leetcode-187:重复的DNA序列
leetcode-187:重复的DNA序列
65 0

热门文章

最新文章