盲人摸象与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降低时,查询树会更加庞大。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 201,681评论 5 474
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,710评论 2 377
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 148,623评论 0 334
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,202评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,232评论 5 363
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,368评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,795评论 3 393
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,461评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,647评论 1 295
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,476评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,525评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,226评论 3 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,785评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,857评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,090评论 1 258
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,647评论 2 348
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,215评论 2 341

推荐阅读更多精彩内容

  • 工作闲暇之余,来趟大宋武侠城,参观了民俗结婚博物馆,看了演出《七侠五义》,还有投资上千万的大型马战《三打祝家庄》。
    贾家甲天下阅读 461评论 0 1
  • 以前我很喜欢写点东西,现在却不愿意写,不知该怎么写,感谢简书提醒我开始~ 提笔就想起了那个人,那些事,已经整整16...
    Amy玉阅读 186评论 1 0
  • PHP之父,毕业于加拿大滑铁卢大学,45岁。 勒道夫开发了PHP的前两个版本并参与了后续版本的开发,2013年加入...
    梁杰_numbbbbb阅读 2,126评论 0 6
  • 我不知道该不该给你写这封信,因为我不知道收到信后你会是怎样的心情。你会开心?还是难过?抑或是什么想法也没有,一如往...
    青苔依旧520阅读 267评论 0 3