序列比对(一)——全局比对Needleman-Wunsch算法

原创

已于 2024-08-29 11:16:14 修改 · 2.6w 阅读

于 2019-04-08 17:51:15 首次发布

本文介绍了序列比对在生物信息学中的重要性，特别是全局比对算法——Needleman-Wunsch算法。该算法利用动态规划计算两个序列的最高得分，并探讨了打分矩阵、回溯过程以及概率论模型。文中还提供了C语言的实现代码，但未涵盖所有可能的错误处理和罚分策略。

原创：hxj7

前言
序列比对是生信领域的一个古老课题，在这一波NGS的浪潮中重新引起大家的广泛关注。由于生物序列的特殊性，在比对的时候允许插入缺失，所以往往是一种不精确匹配。从本文开始，我们陆续学习前人开发出的流行算法。

全局比对算法
所谓全局比对算法，就是根据一个打分矩阵（替换矩阵）计算出两个序列比对最高得分的算法。关于它的介绍网上已经非常多了，我们只需看看其中的关键点及实现代码。
关键点
1. 打分矩阵：
选用不同的打分矩阵或者罚分分值会导致比对结果不同，常用BLAST打分矩阵。
2. 计算比对最高得分的算法：
常用动态规划算法（Needleman-Wunsch算法）。
在这里插入图片描述
图片引自https://www.jianshu.com/p/2b99d0d224a2
3. 打印出最高得分相应的序列比对结果：
根据得分矩阵回溯，如果最优比对结果有多个，全部打印出来。
4. 理解打分系统背后的概率论模型：
比对分值可以理解为匹配模型和随机模型的对数几率比（log-odds ratio）。

实现代码（C代码及示例）
网上的教程给出的代码大多不全，所以我决定还是自己造轮子了。
需要说明的是：
1. 没有对输入进行检查，如果有非AGCT的字符程序可能会出错。
2. 对空位的罚分是线性的，即penalty＝g*d，其中g为空位长度，d为一个空位的罚分。
在以后的文章中，我们会给出仿射型罚分的代码，即
penalty=d+(g-1)*e，其中g为连续空位的长度，d为连续空位中第一个空位的罚分，e为连续空位中第2个及以后空位的罚分。
3. 其他未提及的错误或者不足。

示例
在这里插入图片描述
代码

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSEQ 1000
#define GAP_CHAR '-'

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

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