Blast打分矩阵-Motif检索的福音

发表于 | 更新于: | | 阅读次数: | 浏览量 次
字数统计: 1.9k 字 | 阅读时长 ≈ 10 分钟

前言

众所周知，对于检索精确匹配，模糊匹配，以及特殊的正则匹配。

对于序列检索而言，同样也会序列匹配的要求。至于正则匹配，那就是本文的话题，专业名词叫做Motif检索。序列检索中的正则匹配，是本文的主要话题，专业名词叫做motif。详见https://en.wikipedia.org/wiki/Sequence_motif

正常情况而言，如果公司有自己的序列库，都会采用大名鼎鼎的NCBI blast来作为自己的序列搜索引擎。就功能而言本身是强大的，但是并不能完美支持motif检索，尤其是氨基酸简码X并不能匹配到任意氨基酸简码。用其他软件做自己的motif检索又得构建另一套索引文件，未免增加了服务器的存储。本文提供一种针对NCBI blast的源码修改使之支持motif检索。

重点内容

针对蛋白质序列(氨基酸序列)检索一共实现了8种打分矩阵，具体实现定义在src/util/tables目录下：

打分矩阵

PAM 与BLOSUM，关于两者算法的定义以及区别自行谷歌，本文核心在于如何修改代码支持motif打分而已。

以PAM30算法为例，对应的源码文件是sm_pam30.c，该文件定义了序列中的氨基酸简码两两之间突变的概率，数值越大，可能性越高，可以理解为相似程度。blastp命令在匹配序列中的每个氨基酸的时候，会基于打分矩阵，算出每个序列的相似分数。从源码文件中得知，原来X匹配其他氨基酸的比分都是固定的-1，导致blastp匹配不了其他氨基酸，因此只需要修改打分矩阵即可。

原始文件定义如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
static const TNCBIScore s_Pam30PSM[25 * 25] = {
    /*       A,  R,  N,  D,  C,  Q,  E,  G,  H,  I,  L,  K,  M,
             F,  P,  S,  T,  W,  Y,  V,  B,  J,  Z,  X,  *        */ 
    /*A*/    6, -7, -4, -3, -6, -4, -2, -2, -7, -5, -6, -7, -5,
            -8, -2,  0, -1,-13, -8, -2, -3, -6, -3, -1,-17,
    /*R*/   -7,  8, -6,-10, -8, -2, -9, -9, -2, -5, -8,  0, -4,
            -9, -4, -3, -6, -2,-10, -8, -7, -7, -4, -1,-17,
    /*N*/   -4, -6,  8,  2,-11, -3, -2, -3,  0, -5, -7, -1, -9,
            -9, -6,  0, -2, -8, -4, -8,  6, -6, -3, -1,-17,
    /*D*/   -3,-10,  2,  8,-14, -2,  2, -3, -4, -7,-12, -4,-11,
           -15, -8, -4, -5,-15,-11, -8,  6,-10,  1, -1,-17,
    /*C*/   -6, -8,-11,-14, 10,-14,-14, -9, -7, -6,-15,-14,-13,
           -13, -8, -3, -8,-15, -4, -6,-12, -9,-14, -1,-17,
    /*Q*/   -4, -2, -3, -2,-14,  8,  1, -7,  1, -8, -5, -3, -4,
           -13, -3, -5, -5,-13,-12, -7, -3, -5,  6, -1,-17,
    /*E*/   -2, -9, -2,  2,-14,  1,  8, -4, -5, -5, -9, -4, -7,
           -14, -5, -4, -6,-17, -8, -6,  1, -7,  6, -1,-17,
    /*G*/   -2, -9, -3, -3, -9, -7, -4,  6, -9,-11,-10, -7, -8,
            -9, -6, -2, -6,-15,-14, -5, -3,-10, -5, -1,-17,
    /*H*/   -7, -2,  0, -4, -7,  1, -5, -9,  9, -9, -6, -6,-10,
            -6, -4, -6, -7, -7, -3, -6, -1, -7, -1, -1,-17,
    /*I*/   -5, -5, -5, -7, -6, -8, -5,-11, -9,  8, -1, -6, -1,
            -2, -8, -7, -2,-14, -6,  2, -6,  5, -6, -1,-17,
    /*L*/   -6, -8, -7,-12,-15, -5, -9,-10, -6, -1,  7, -8,  1,
            -3, -7, -8, -7, -6, -7, -2, -9,  6, -7, -1,-17,
    /*K*/   -7,  0, -1, -4,-14, -3, -4, -7, -6, -6, -8,  7, -2,
           -14, -6, -4, -3,-12, -9, -9, -2, -7, -4, -1,-17,
    /*M*/   -5, -4, -9,-11,-13, -4, -7, -8,-10, -1,  1, -2, 11,
            -4, -8, -5, -4,-13,-11, -1,-10,  0, -5, -1,-17,
    /*F*/   -8, -9, -9,-15,-13,-13,-14, -9, -6, -2, -3,-14, -4,
             9,-10, -6, -9, -4,  2, -8,-10, -2,-13, -1,-17,
    /*P*/   -2, -4, -6, -8, -8, -3, -5, -6, -4, -8, -7, -6, -8,
           -10,  8, -2, -4,-14,-13, -6, -7, -7, -4, -1,-17,
    /*S*/    0, -3,  0, -4, -3, -5, -4, -2, -6, -7, -8, -4, -5,
            -6, -2,  6,  0, -5, -7, -6, -1, -8, -5, -1,-17,
    /*T*/   -1, -6, -2, -5, -8, -5, -6, -6, -7, -2, -7, -3, -4,
            -9, -4,  0,  7,-13, -6, -3, -3, -5, -6, -1,-17,
    /*W*/  -13, -2, -8,-15,-15,-13,-17,-15, -7,-14, -6,-12,-13,
            -4,-14, -5,-13, 13, -5,-15,-10, -7,-14, -1,-17,
    /*Y*/   -8,-10, -4,-11, -4,-12, -8,-14, -3, -6, -7, -9,-11,
             2,-13, -7, -6, -5, 10, -7, -6, -7, -9, -1,-17,
    /*V*/   -2, -8, -8, -8, -6, -7, -6, -5, -6,  2, -2, -9, -1,
            -8, -6, -6, -3,-15, -7,  7, -8,  0, -6, -1,-17,
    /*B*/   -3, -7,  6,  6,-12, -3,  1, -3, -1, -6, -9, -2,-10,
           -10, -7, -1, -3,-10, -6, -8,  6, -8,  0, -1,-17,
    /*J*/   -6, -7, -6,-10, -9, -5, -7,-10, -7,  5,  6, -7,  0,
            -2, -7, -8, -5, -7, -7,  0, -8,  6, -6, -1,-17,
    /*Z*/   -3, -4, -3,  1,-14,  6,  6, -5, -1, -6, -7, -4, -5,
           -13, -4, -5, -6,-14, -9, -6,  0, -6,  6, -1,-17,
    /*X*/   -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
            -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,-17,
    /***/  -17,-17,-17,-17,-17,-17,-17,-17,-17,-17,-17,-17,-17,
           -17,-17,-17,-17,-17,-17,-17,-17,-17,-17,-17,  1
};

修改后的打分矩阵，

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
static const TNCBIScore s_Pam30PSM[25 * 25] = {
    /*       A,  R,  N,  D,  C,  Q,  E,  G,  H,  I,  L,  K,  M,
             F,  P,  S,  T,  W,  Y,  V,  B,  J,  Z,  X,  *        */ 
    /*A*/    6, -7, -4, -3, -6, -4, -2, -2, -7, -5, -6, -7, -5,
            -8, -2,  0, -1,-13, -8, -2, -3, -6, -3,  5,-17,
    /*R*/   -7,  8, -6,-10, -8, -2, -9, -9, -2, -5, -8,  0, -4,
            -9, -4, -3, -6, -2,-10, -8, -7, -7, -4,  7,-17,
    /*N*/   -4, -6,  8,  2,-11, -3, -2, -3,  0, -5, -7, -1, -9,
            -9, -6,  0, -2, -8, -4, -8,  6, -6, -3,  7,-17,
    /*D*/   -3,-10,  2,  8,-14, -2,  2, -3, -4, -7,-12, -4,-11,
           -15, -8, -4, -5,-15,-11, -8,  6,-10,  1,  7,-17,
    /*C*/   -6, -8,-11,-14, 10,-14,-14, -9, -7, -6,-15,-14,-13,
           -13, -8, -3, -8,-15, -4, -6,-12, -9,-14,  9,-17,
    /*Q*/   -4, -2, -3, -2,-14,  8,  1, -7,  1, -8, -5, -3, -4,
           -13, -3, -5, -5,-13,-12, -7, -3, -5,  6,  7,-17,
    /*E*/   -2, -9, -2,  2,-14,  1,  8, -4, -5, -5, -9, -4, -7,
           -14, -5, -4, -6,-17, -8, -6,  1, -7,  6,  7,-17,
    /*G*/   -2, -9, -3, -3, -9, -7, -4,  6, -9,-11,-10, -7, -8,
            -9, -6, -2, -6,-15,-14, -5, -3,-10, -5,  5,-17,
    /*H*/   -7, -2,  0, -4, -7,  1, -5, -9,  9, -9, -6, -6,-10,
            -6, -4, -6, -7, -7, -3, -6, -1, -7, -1,  8,-17,
    /*I*/   -5, -5, -5, -7, -6, -8, -5,-11, -9,  8, -1, -6, -1,
            -2, -8, -7, -2,-14, -6,  2, -6,  5, -6,  7,-17,
    /*L*/   -6, -8, -7,-12,-15, -5, -9,-10, -6, -1,  7, -8,  1,
            -3, -7, -8, -7, -6, -7, -2, -9,  6, -7,  6,-17,
    /*K*/   -7,  0, -1, -4,-14, -3, -4, -7, -6, -6, -8,  7, -2,
           -14, -6, -4, -3,-12, -9, -9, -2, -7, -4,  6,-17,
    /*M*/   -5, -4, -9,-11,-13, -4, -7, -8,-10, -1,  1, -2, 11,
            -4, -8, -5, -4,-13,-11, -1,-10,  0, -5, 10,-17,
    /*F*/   -8, -9, -9,-15,-13,-13,-14, -9, -6, -2, -3,-14, -4,
             9,-10, -6, -9, -4,  2, -8,-10, -2,-13,  8,-17,
    /*P*/   -2, -4, -6, -8, -8, -3, -5, -6, -4, -8, -7, -6, -8,
           -10,  8, -2, -4,-14,-13, -6, -7, -7, -4,  7,-17,
    /*S*/    0, -3,  0, -4, -3, -5, -4, -2, -6, -7, -8, -4, -5,
            -6, -2,  6,  0, -5, -7, -6, -1, -8, -5,  5,-17,
    /*T*/   -1, -6, -2, -5, -8, -5, -6, -6, -7, -2, -7, -3, -4,
            -9, -4,  0,  7,-13, -6, -3, -3, -5, -6,  6,-17,
    /*W*/  -13, -2, -8,-15,-15,-13,-17,-15, -7,-14, -6,-12,-13,
            -4,-14, -5,-13, 13, -5,-15,-10, -7,-14, 12,-17,
    /*Y*/   -8,-10, -4,-11, -4,-12, -8,-14, -3, -6, -7, -9,-11,
             2,-13, -7, -6, -5, 10, -7, -6, -7, -9,  9,-17,
    /*V*/   -2, -8, -8, -8, -6, -7, -6, -5, -6,  2, -2, -9, -1,
            -8, -6, -6, -3,-15, -7,  7, -8,  0, -6,  6,-17,
    /*B*/   -3, -7,  6,  6,-12, -3,  1, -3, -1, -6, -9, -2,-10,
           -10, -7, -1, -3,-10, -6, -8,  6, -8,  0,  5,-17,
    /*J*/   -6, -7, -6,-10, -9, -5, -7,-10, -7,  5,  6, -7,  0,
            -2, -7, -8, -5, -7, -7,  0, -8,  6, -6,  5,-17,
    /*Z*/   -3, -4, -3,  1,-14,  6,  6, -5, -1, -6, -7, -4, -5,
           -13, -4, -5, -6,-14, -9, -6,  0, -6,  6,  5,-17,
    /*X*/    5,  7,  7,  7,  9,  7,  7,  5,  8,  7,  6,  6, 10,
             8,  7,  5,  6, 12,  9,  6,  5,  5,  5, 13,-17,
    /***/  -17,-17,-17,-17,-17,-17,-17,-17,-17,-17,-17,-17,-17,
           -17,-17,-17,-17,-17,-17,-17,-17,-17,-17,-17,  1
};

为什么X与其他氨基酸的分值不设置成相同的正整数呢，大家可以思考下。

剩下的工作就是编译出来一个新的blastp命令即可。

详细参考另一篇博文： Compile and debug NCBI blast+ source code - NCBI blast源码编译调试教程