盲人摸象与BLAST 算法

“看到BLAST算法，我一下子想起了盲人摸象这个小故事，只不过这一次，要讲一个新故事了，权且叫它‘新盲人摸象’吧。”

国王有一天突然想和当初玩摸象的盲人再玩个新游戏，就把一个自己最喜爱的盲人叫了过来。

国王对他说：“嘿，你还记得上次我告诉你们的象的全貌吗？”

盲人：“当然记得啦，它的牙齿像一个大萝卜，鼻子像一段长长的水管，耳朵像大蒲扇，四条腿像柱子，尾巴像条草绳，···”

国王：“好啦好啦，看来你都记住了。我下面带你到动物园里面去，你要靠摸准确地找出大象哦。”

盲人：“谨遵王命，喵喵喵。”

国王：“汪汪汪。”

于是国王带了盲人到动物园去，刚进去，一只鸡咯咯叫着从他们身边跑过，盲人一欠身就把它捞了起来上下其手：“喔，嘴又尖又短像个小图钉，摸不到鼻子和耳朵，浑身都是鸡毛掸子的毛毛，只有两条小细腿，尾巴似乎也摸不到，看来完全不是哦。”

鸡奋力挣脱了他的咸猪手，委屈地咯咯着。

接着他们又路遇鸭、羊、狗、鸡，盲人通过抚摸牙齿、鼻子、耳朵、腿、尾巴的方式都排除了。

于是二人头顶着鸡毛、鸭毛、狗毛继续往前走。

在驴马馆，盲人有了惊喜的发现：“嗨呀，尾巴像草绳！嗨呀，四条腿！”

国王站在一旁，憋着笑，要看盲人出丑。

盲人沿着尾巴开始往别的地方摸，“背上这些毛怎么肥四？啊，说不定是长了草。这耳朵好小哦，根本不能扇风嘛。也摸不到长鼻子和大牙，看来根本不是嘛。”

盲人又依次因为野猪的獠牙、长颈鹿的腿、大食蚁兽的长鼻子而空欢喜一场。终于到了大象馆，盲人此时已谨慎了许多，“鼻子对上了”，他接着往上摸，“牙也能对上”，“耳朵是蒲扇，腿是柱子，尾巴是草绳，身上没毛”。

“是大象！”他转身兴奋地告诉国王。“妙啊！”国王也为自己看到了盲人笨拙而有趣的比对而兴奋。“走走走，我要给你升官，就······去管光垫鬃菊吧，干这个，眼睛好不好使都是次要的。”

盲人也很兴奋：“国王万岁，王啊，我的腿有些疼，不知道是不是被什么咬了。”

国王：“啊，好像是被你摸的那条狗咬了，我刚刚似乎也中招了。不要紧，太医刚刚从极寒之地带来几针异喵，一会就顺便给我们打上。”

后来，再也没人见过国王和盲人，这个王国，从此就没有新的故事了。新盲人摸象，是流传下来的最后一个故事。

故事讲完了。大象，就对应着我们提交的DNA或是氨基酸序列，动物园，就是包含着许许多多其他序列的数据库。而我们的目标，就是应用BLAST算法，给出提交序列与数据库序列的最佳局部联配结果。

具体过程

• 做查询序列的k-字母单词表，对于DNA序列通常我们取k为11，这里我们就令k为3。这一步就好比把牙齿、鼻子、耳朵、腿、尾巴这些大象独有的特征都提取出来。

序列：TGCGGAGCGG
获取单词表：TGC、GCG、CGG、GGA、GAG、AGC（GCG、CGG）<最后两个前面已经出现过了，这里舍去>

• 以上面的单词表为基础，获取查询树。这一步就好比放宽搜索的范围，比如说，认为拥有短一些的鼻子的动物也可能是大象，认为只要是四条腿的动物也可能是大象。

列举所有的k_字母表：如AAA、AAT、AAC、AAG等等。
依次将这些字母依据打分矩阵和上一步获得的k_字母匹配打分，获得一个标志着匹配度的分数。
设置一个匹配值T，将所有分数超过T的k_字母加入查询树当中去。T越大，越严格，当T设为最大时，种子中的k_字母就是上一步所获得的所有K字母。

• 确定数据库当中所有包含查询树中的k_字母的序列。当然，这一步并不是那种一个一个的比对，数据库里面已经提前建好了哈希索引。找到的所有序列中的那一小段k_字母，称之为种子(seed)。这一步就好比挨个对动物上下其手，看看有没有哪个动物身上恰好有与大象身上的特征完全一样或者相近的特征。

• 将种子延伸成高打分匹配片断(high-scoring segment pair, HSP）。也就是说对于每一个匹配的种子，我们以动态规划的算法向两侧延伸，在延伸过程中，联配值S也在不断地变化，当它低于提前设定好的一个临界值X时，联配结束，获得一段HSP。这就相当于摸到某个特征后继续往旁边摸，比较它与大象的其他特征是否吻合，如果不吻合之处太多了，那就不摸了。

• 当然，要注意到，我们要解决的问题的本质与一开始所说的故事还是有些区别的。我们的目标并不是找到一模一样的大象，而是找到不同的动物身上某个与大象足够接近的特征。

部分实现

这次比较新的东西主要是构建查询树这一部分，其他在前面写过，此处略去不表。

计分规则，匹配计3分，AT或CG配对计1分，其他不匹配情况计-3分。

def mark(word1, word2):
    # 给两个单词打分
    ld = ["AT", "TA", "CG", "GC"]
    score = 0
    for a, b in zip(word1, word2):
        if a == b:
            score += 3
        elif a + b in ld:
            score += 1
        else:
            score += -3
    return score

def find_all(L,K,A,word):
    # 参数分别为字符列表、单词长度、包含结果的列表、当前单词
    if len(word) == K:
        A.append(word)
        return
    for i in L:
        find_all(L,K,A,word+i)

def get_tree(seq, K, T):
    # 参数分别为DNA序列，单词长度K，分数阈值T
    F = []
    for i in range(0, len(seq)-K+1):
        F.append(seq[i:i+K])
    # 先把列表转换为集合实现去重，再转换为列表
    F = list(set(F))
    A = []
    find_all(['A','T','C','G'], K, A, "")
    TREE = []
    for i in F:
        for j in A:
            if mark(i,j) >= T:
                TREE.append(j)
    print("查询树中共包含", len(TREE), "个节点")
    for t in TREE:
        print(t)

调用函数

get_tree("ATCGTCCA",3,9)
get_tree("ATCGTCCA",3,6)
get_tree("ATCGTCCA",3,3)

结果如下

查询树中共包含 6 个节点
CGT
CCA
TCG
TCC
GTC
ATC
查询树中共包含 24 个节点
CCT
CGA
CGT
GGT
CCA
CCT
CGA
GCA
ACG
TCC
TCG
TGG
ACC
TCC
TCG
TGC
CTC
GAC
GTC
GTG
AAC
ATC
ATG
TTC
查询树中共包含 84 个节点
AGT
TGT
CAT
CTT
CCA
CCT
CGA
CGT
CGC
CGG
GCA
GCT
GGA
GGT
ACA
TCA
CAA
CTA
CCA
CCT
CCC
CCG
CGA
CGT
GCA
GCT
GGA
GGT
ACC
ACG
AGC
AGG
TAG
TTG
TCA
TCT
TCC
TCG
TGC
TGG
CCG
GCG
ACC
ACG
AGC
AGG
TAC
TTC
TCA
TCT
TCC
TCG
TGC
TGG
CCC
GCC
ATC
TTC
CAC
CAG
CTC
CTG
GAC
GAG
GTA
GTT
GTC
GTG
GCC
GGC
AAC
AAG
ATA
ATT
ATC
ATG
ACC
AGC
TAC
TAG
TTC
TTG
CTC
GTC

可以看到，当K降低时，查询树会更加庞大。