专栏首页生信了序列比对(一)全局比对Needleman-Wunsch算法

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

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

全局比对算法

所谓全局比对算法,就是根据一个打分矩阵(替换矩阵)计算出两个序列比对最高得分的算法。关于它的介绍网上已经非常多了,我们只需看看其中的关键点及实现代码。

关键点

  1. 打分矩阵: 选用不同的打分矩阵或者罚分分值会导致比对结果不同,常用BLAST打分矩阵。
  2. 计算比对最高得分的算法: 常用动态规划算法(Needleman-Wunsch算法)。

图片引自https://www.jianshu.com/p/2b99d0d224a2

  1. 打印出最高得分相应的序列比对结果: 根据得分矩阵回溯,如果最优比对结果有多个,全部打印出来。
  2. 理解打分系统背后的概率论模型: 比对分值可以理解为匹配模型和随机模型的对数几率比(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 max3(float a, float b, float c) {
    float f = a > b ? a : b;
    return f > c ? f : c;
}

// 替换矩阵:match分值为5,mismatch分值为-4
// 数组下标是两个字符的ascii码减去65之后的和
float FMatrix[] = {
    5, 0, -4, 0, 5, 0, -4, 0, -4, 0,
    0, 0, 5, 0, 0, 0, 0, 0, 0, -4,
    0, -4, 0, 0, 0, -4, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 5
};

float getFScore(char a, char b) {
    return FMatrix[a + b - 'A' - 'A'];
}

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);
    float gap = -2.5;     // 对空位的罚分
    float m1, m2, m3, maxm;
    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;
        }
    }
    aUnit[0][0]->M = 0;
    for (i = 1; i <= m; i++) {
        aUnit[i][0]->M = gap * i;
        aUnit[i][0]->W1 = 1;
    }
    for (j = 1; j <= n; j++) {
         aUnit[0][j]->M = gap * 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 + gap;
            m2 = aUnit[i - 1][j - 1]->M + getFScore(s[i - 1], r[j - 1]);
            m3 = aUnit[i][j - 1]->M + gap;
            maxm = max3(m1, m2, m3);
            aUnit[i][j]->M = maxm;
            if (m1 == maxm) aUnit[i][j]->W1 = 1;
            if (m2 == maxm) aUnit[i][j]->W2 = 1;
            if (m3 == maxm) aUnit[i][j]->W3 = 1;
        }
    }
/*
    // 打印得分矩阵
    for (i = 0; i <= m; i++) {
        for (j = 0; j <= n; j++)
            printf("%f ", aUnit[i][j]->M);
        printf("\n");
    }
*/
    printf("max score: %f\n", aUnit[m][n]->M);
    // 打印最优比对结果,如果有多个,全部打印
    // 递归法
    if ((salign = (char*) malloc(sizeof(char) * (m + n + 1))) == NULL) {
        fputs("Error: Out of space!\n", stderr);
        exit(1);
    }
    if ((ralign = (char*) malloc(sizeof(char) * (m + n + 1))) == 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);
}

本文分享自微信公众号 - 生信了(gh_ed36a29a9a9d),作者:hxj7

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-04-08

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 序列匹配(五)重复匹配问题的动态规划算法

    问题及算法描述 更具体地描述上面的问题:有序列x和y,其中y是包含结构域的序列,x是要从中找到多重匹配的序列。将x分割成一段一段的不交叠的子序列,这些子序列要...

    一只羊
  • 序列比对(四)Smith-Waterman算法之仿射罚分

    关于全局联配,局部联配以及仿射罚分模型的介绍可参见前文: 序列比对(一)全局比对Needleman-Wunsch算法 序列比对(二)Needleman-Wuns...

    一只羊
  • 序列比对(25)编辑距离

    编辑距离的求解过程和全局比对是十分相似的(关于全局比对,可以参见前文《序列比对(一)全局比对Needleman-Wunsch算法》),都需要全部符号参与比对,都...

    一只羊
  • 序列匹配(五)重复匹配问题的动态规划算法

    问题及算法描述 更具体地描述上面的问题:有序列x和y,其中y是包含结构域的序列,x是要从中找到多重匹配的序列。将x分割成一段一段的不交叠的子序列,这些子序列要...

    一只羊
  • 序列比对(25)编辑距离

    编辑距离的求解过程和全局比对是十分相似的(关于全局比对,可以参见前文《序列比对(一)全局比对Needleman-Wunsch算法》),都需要全部符号参与比对,都...

    一只羊
  • 序列比对(三)局部联配Smith-Waterman算法

    关于全局联配的介绍可参见前文: 序列比对(一)全局比对Needleman-Wunsch算法 序列比对(二)Needleman-Wunsch算法之仿射罚分

    一只羊
  • 见缝插针游戏--实现转圈

    用户2965768
  • [Setting]用VS2008将类封装为静态库library

    原文链接:http://blog.csdn.net/humanking7/article/details/50726271

    祥知道
  • 什么是渗透测试?

    渗透测试告诉系统上采用的现有防御措施是否足够强大,可以防止任何安全漏洞。渗透测试报告还建议了可以采取的对策,以减少系统被黑客入侵的风险。

    用户7466307
  • 推荐三款机器学习论文搜索利器

    简单介绍:arXiv是个提交论文预印本(preprint)的平台,里面的论文都没有经过同行评审(peer review),所以文章质量参差不齐,但却会比较新颖,...

    石晓文

扫码关注云+社区

领取腾讯云代金券