“看到BLAST算法，我一下子想起了盲人摸象這個(gè)小故事，只不過(guò)這一次，要講一個(gè)新故事了，權(quán)且叫它‘新盲人摸象’吧?！?/h3>

國(guó)王有一天突然想和當(dāng)初玩摸象的盲人再玩?zhèn)€新游戲，就把一個(gè)自己最喜愛(ài)的盲人叫了過(guò)來(lái)。

國(guó)王對(duì)他說(shuō)：“嘿，你還記得上次我告訴你們的象的全貌嗎？”

盲人：“當(dāng)然記得啦，它的牙齒像一個(gè)大蘿卜，鼻子像一段長(zhǎng)長(zhǎng)的水管，耳朵像大蒲扇，四條腿像柱子，尾巴像條草繩，···”

國(guó)王：“好啦好啦，看來(lái)你都記住了。我下面帶你到動(dòng)物園里面去，你要靠摸準(zhǔn)確地找出大象哦?！?/p>

盲人：“謹(jǐn)遵王命，喵喵喵?！?/p>

國(guó)王：“汪汪汪?！?/p>

于是國(guó)王帶了盲人到動(dòng)物園去，剛進(jìn)去，一只雞咯咯叫著從他們身邊跑過(guò)，盲人一欠身就把它撈了起來(lái)上下其手：“喔，嘴又尖又短像個(gè)小圖釘，摸不到鼻子和耳朵，渾身都是雞毛撣子的毛毛，只有兩條小細(xì)腿，尾巴似乎也摸不到，看來(lái)完全不是哦?！?/p>

雞奮力掙脫了他的咸豬手，委屈地咯咯著。

接著他們又路遇鴨、羊、狗、雞，盲人通過(guò)撫摸牙齒、鼻子、耳朵、腿、尾巴的方式都排除了。

于是二人頭頂著雞毛、鴨毛、狗毛繼續(xù)往前走。

在驢馬館，盲人有了驚喜的發(fā)現(xiàn)：“嗨呀，尾巴像草繩！嗨呀，四條腿！”

國(guó)王站在一旁，憋著笑，要看盲人出丑。

盲人沿著尾巴開(kāi)始往別的地方摸，“背上這些毛怎么肥四？啊，說(shuō)不定是長(zhǎng)了草。這耳朵好小哦，根本不能扇風(fēng)嘛。也摸不到長(zhǎng)鼻子和大牙，看來(lái)根本不是嘛?！?/p>

盲人又依次因?yàn)橐柏i的獠牙、長(zhǎng)頸鹿的腿、大食蟻獸的長(zhǎng)鼻子而空歡喜一場(chǎng)。終于到了大象館，盲人此時(shí)已謹(jǐn)慎了許多，“鼻子對(duì)上了”，他接著往上摸，“牙也能對(duì)上”，“耳朵是蒲扇，腿是柱子，尾巴是草繩，身上沒(méi)毛”。

“是大象！”他轉(zhuǎn)身興奮地告訴國(guó)王?！懊畎。　眹?guó)王也為自己看到了盲人笨拙而有趣的比對(duì)而興奮。“走走走，我要給你升官，就······去管光墊鬃菊吧，干這個(gè)，眼睛好不好使都是次要的?！?/p>

盲人也很興奮：“國(guó)王萬(wàn)歲，王啊，我的腿有些疼，不知道是不是被什么咬了?！?/p>

國(guó)王：“啊，好像是被你摸的那條狗咬了，我剛剛似乎也中招了。不要緊，太醫(yī)剛剛從極寒之地帶來(lái)幾針異喵，一會(huì)就順便給我們打上?！?/p>

后來(lái)，再也沒(méi)人見(jiàn)過(guò)國(guó)王和盲人，這個(gè)王國(guó)，從此就沒(méi)有新的故事了。新盲人摸象，是流傳下來(lái)的最后一個(gè)故事。

故事講完了。大象，就對(duì)應(yīng)著我們提交的DNA或是氨基酸序列，動(dòng)物園，就是包含著許許多多其他序列的數(shù)據(jù)庫(kù)。而我們的目標(biāo)，就是應(yīng)用BLAST算法，給出提交序列與數(shù)據(jù)庫(kù)序列的最佳局部聯(lián)配結(jié)果。

具體過(guò)程

• 做查詢序列的k-字母單詞表，對(duì)于DNA序列通常我們?nèi)為11，這里我們就令k為3。這一步就好比把牙齒、鼻子、耳朵、腿、尾巴這些大象獨(dú)有的特征都提取出來(lái)。

序列：TGCGGAGCGG
獲取單詞表：TGC、GCG、CGG、GGA、GAG、AGC（GCG、CGG）<最后兩個(gè)前面已經(jīng)出現(xiàn)過(guò)了，這里舍去>

• 以上面的單詞表為基礎(chǔ)，獲取查詢樹(shù)。這一步就好比放寬搜索的范圍，比如說(shuō)，認(rèn)為擁有短一些的鼻子的動(dòng)物也可能是大象，認(rèn)為只要是四條腿的動(dòng)物也可能是大象。

列舉所有的k_字母表：如AAA、AAT、AAC、AAG等等。
依次將這些字母依據(jù)打分矩陣和上一步獲得的k_字母匹配打分，獲得一個(gè)標(biāo)志著匹配度的分?jǐn)?shù)。
設(shè)置一個(gè)匹配值T，將所有分?jǐn)?shù)超過(guò)T的k_字母加入查詢樹(shù)當(dāng)中去。T越大，越嚴(yán)格，當(dāng)T設(shè)為最大時(shí)，種子中的k_字母就是上一步所獲得的所有K字母。

• 確定數(shù)據(jù)庫(kù)當(dāng)中所有包含查詢樹(shù)中的k_字母的序列。當(dāng)然，這一步并不是那種一個(gè)一個(gè)的比對(duì)，數(shù)據(jù)庫(kù)里面已經(jīng)提前建好了哈希索引。找到的所有序列中的那一小段k_字母，稱(chēng)之為種子(seed)。這一步就好比挨個(gè)對(duì)動(dòng)物上下其手，看看有沒(méi)有哪個(gè)動(dòng)物身上恰好有與大象身上的特征完全一樣或者相近的特征。

• 將種子延伸成高打分匹配片斷(high-scoring segment pair, HSP）。也就是說(shuō)對(duì)于每一個(gè)匹配的種子，我們以動(dòng)態(tài)規(guī)劃的算法向兩側(cè)延伸，在延伸過(guò)程中，聯(lián)配值S也在不斷地變化，當(dāng)它低于提前設(shè)定好的一個(gè)臨界值X時(shí)，聯(lián)配結(jié)束，獲得一段HSP。這就相當(dāng)于摸到某個(gè)特征后繼續(xù)往旁邊摸，比較它與大象的其他特征是否吻合，如果不吻合之處太多了，那就不摸了。

• 當(dāng)然，要注意到，我們要解決的問(wèn)題的本質(zhì)與一開(kāi)始所說(shuō)的故事還是有些區(qū)別的。我們的目標(biāo)并不是找到一模一樣的大象，而是找到不同的動(dòng)物身上某個(gè)與大象足夠接近的特征。

部分實(shí)現(xiàn)

這次比較新的東西主要是構(gòu)建查詢樹(shù)這一部分，其他在前面寫(xiě)過(guò)，此處略去不表。

計(jì)分規(guī)則，匹配計(jì)3分，AT或CG配對(duì)計(jì)1分，其他不匹配情況計(jì)-3分。

def mark(word1, word2):
    # 給兩個(gè)單詞打分
    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):
    # 參數(shù)分別為字符列表、單詞長(zhǎng)度、包含結(jié)果的列表、當(dāng)前單詞
    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):
    # 參數(shù)分別為DNA序列，單詞長(zhǎng)度K，分?jǐn)?shù)閾值T
    F = []
    for i in range(0, len(seq)-K+1):
        F.append(seq[i:i+K])
    # 先把列表轉(zhuǎn)換為集合實(shí)現(xiàn)去重，再轉(zhuǎn)換為列表
    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("查詢樹(shù)中共包含", len(TREE), "個(gè)節(jié)點(diǎn)")
    for t in TREE:
        print(t)

調(diào)用函數(shù)

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

結(jié)果如下

查詢樹(shù)中共包含 6 個(gè)節(jié)點(diǎn)
CGT
CCA
TCG
TCC
GTC
ATC
查詢樹(shù)中共包含 24 個(gè)節(jié)點(diǎn)
CCT
CGA
CGT
GGT
CCA
CCT
CGA
GCA
ACG
TCC
TCG
TGG
ACC
TCC
TCG
TGC
CTC
GAC
GTC
GTG
AAC
ATC
ATG
TTC
查詢樹(shù)中共包含 84 個(gè)節(jié)點(diǎn)
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

可以看到，當(dāng)K降低時(shí)，查詢樹(shù)會(huì)更加龐大。