如何让BLAST返回最优的一个搜索结果,看看你没有有进坑

大部分时候,我们都是看着别人的教程,然后尝试处理自己的数据,结果跑完了,如果和预期相符合就不会怀疑这个工具有啥问题。如果你要学习生物信息学,那么有一个信条一定要记住,不要盲目相信软件的输出结果,多多检查。

BLAST作为最常用的一个序列比对工具,对于大多数用户而言,都是用于找找一个基因在其他物种的同源基因,或者看看一条引物在基因组上是否有多处匹配,大多不会限制输出的匹配数目。

但是对于一些分析流程,你可能希望只要输出一个最佳匹配就好,这个时候你可能会用到两个参数,一个是-max_target_seqs,另一个是-num_alignments, 它们的参数说明如下

 -max_target_seqs <Integer, >=1>                                                                                                                
   Maximum number of aligned sequences to keep                                                                                                  
   Not applicable for outfmt <= 4
  * Incompatible with:  num_descriptions, num_alignments
 -num_alignments <Integer, >=0>                                                                                                                 
   Number of database sequences to show alignments for
   * Incompatible with:  max_target_seqs

如果你认为" aligned sequences" 或者"database sequences" 指的每个序列只会输出一个匹配的话,那么你掉进坑了。

首先给定一个检索序列, 命名为"query.fa"

>query1
CTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC

选择拟南芥的前1000行序列作为 测试集1

head -n 1000 Athaliana.fa > test1.fa
makeblastdb -in  test1.fa -out test/test1  -dbtype nucl

直接搜索,找不到任何的结果

blastn -query query.fa -db test/test1 -outfmt 6

将搜索的序列加入 测试集1,得到 测试集2

cat query.fa test1.fa > test2.fa
makeblastdb -in test2.fa -out test/test2  -dbtype nucl

此时搜索,可以找到一个结果

blastn -query query.fa -db test/test2 -outfmt 6

如果将同一段序列,重复多次,中间加入一些无意义的序列(如下),之后再加入到 测试集1 中,得 测试集3

>query2
CTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGCCCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAATCTTTAAATCCTACATCCAT
GAATCCCTAAATACCTAATTCCCTAAACCCGAAACCGGTTTCTCTGGTTGAAAATCATTGTGTATATAATGATAATTTT
ATCGTTTTTATGTAATTGCTTATTGTTGTGTGTAGATTTTTTAAAAATATCATTTGAGGTCAATACAAATCCTATTTCT
TGTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTA
CATTTGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTA
TGTTTGGACATTTATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTATCAATTTATCTTTTGTGGGAAAATTATT
TAGTTGTAGGGATGAAGTCTTTCTTCGTTGTTGTTACGCTTGTCATCTCATCTCTCAATGATATGGGATGGTCCTTTAG
CATTTATTCTGAAGTTCTTCTGCTTGATGATTTTATCCTTAGCCAAAAGGATTGGTGGTTTGAAGACACATCATATCAA
AAAAGCTATCGCCTCGACGATGCTCTATTTCTATCCTTGTAGCACACATTTTGGCACTCAAAAAAGTATTTTTAGATGTCTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC
CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAATCTTTAAATCCTACATCCAT
GAATCCCTAAATACCTAATTCCCTAAACCCGAAACCGGTTTCTCTGGTTGAAAATCATTGTGTATATAATGATAATTTT
ATCGTTTTTATGTAATTGCTTATTGTTGTGTGTAGATTTTTTAAAAATATCATTTGAGGTCAATACAAATCCTATTTCT
TGTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTA
CATTTGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTA
TGTTTGGACATTTATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTATCAATTTATCTTTTGTGGGAAAATTATT
TAGTTGTAGGGATGAAGTCTTTCTTCGTTGTTGTTACGCTTGTCATCTCATCTCTCAATGATATGGGATGGTCCTTTAG
CATTTATTCTGAAGTTCTTCTGCTTGATGATTTTATCCTTAGCCAAAAGGATTGGTGGTTTGAAGACACATCATATCAA
AAAAGCTATCGCCTCGACGATGCTCTATTTCTATCCTTGTAGCACACATTTTGGCACTCAAAAAAGTATTTTTAGATGT
CTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC

此时搜索,会得到三个结果

cat query2.fa test1.fa > test3.fa
makeblastdb -in test3.fa -out test/test3 -dbtype nucl
blastn -query query.fa -db test/test3 -outfmt 6
# 结果如下
query1  query2  100.000 100 0   0   1   100 1   100 8.11e-50    185
query1  query2  100.000 100 0   0   1   100 812 911 8.11e-50    185
query1  query2  100.000 100 0   0   1   100 1623    1722    8.11e-50    185

加入我希望只要一个结果,根据我之前的经验,我加入一个参数-num_alignments 1, 帮助里写的是 "Number of database sequences to show alignments for" 我一直认为是,只返回第一个找到的序列,但是运行如下代码,结果出乎我意料

blastn -query query.fa -db test/test3 -outfmt 6 -num_alignments 1
# 结果如下
query1  query2  100.000 100 0   0   1   100 1   100 8.11e-50    185
query1  query2  100.000 100 0   0   1   100 812 911 8.11e-50    185
query1  query2  100.000 100 0   0   1   100 1623    1722    8.11e-50    185

他居然会得到三个结果. 为什么会出现这种情况?

我的猜想是:这里的 Number of database sequencessequences 不是匹配上的序列,而是匹配上的序列对应所在整条序列。那么根据我的猜想, 当数据库中中有A,B,C三条序列,其中A有三个潜在匹配,B有两个潜在匹配,C没有匹配,A, B, C在建立搜索数据库的顺序会影响最终的输出结果,也就是当-num_alignments 1时,blastn搜索到A后,就会输出A里面的三个检索结果就不再继续,如果先搜索到B,就会输出一个匹配,不再继续。当-num_alignments 1会输出四个匹配。

为了验证假设,我将query.fa 和 query2.fa 以不同排列顺序构建测试数据集,发现输出结果的确和顺序有关。

#顺序1
cat query.fa query2.fa test1.fa > test4.fa
makeblastdb -in test4.fa -out test/test4 -dbtype nucl
blastn -query query.fa -db test/test4 -outfmt 6 -num_alignments 1
# 结果如下
query1  query2  100.000 100 0   0   1   100 1   100 8.12e-50    185
query1  query2  100.000 100 0   0   1   100 812 911 8.12e-50    185
query1  query2  100.000 100 0   0   1   100 1623    1722    8.12e-50    185
# 顺序2 
cat   query2.fa  query.fa   test1.fa > test5.fa
makeblastdb -in test5.fa -out test/test5 -dbtype nucl
blastn -query query.fa -db test/test5 -outfmt 6 -num_alignments 1
# 结果如下
query1  query1  100.000 100 0   0   1   100 1   100 8.12e-50    185

这个情况是我用一条短序列搜索全基因组数据时出现(如果是CDS序列库则不会出现一条序列出现多个匹配的情况), 在你设计引物的时候,如果在一条染色体上能够找到多个相似序列,那么这种引物序列就是不合格的。

但是如果你一定要保证每条序列在每个染色体都只要第一个高质量匹配,那么如何解决呢?

我们用到一个不怎么起眼的参数,-max_hsps 1, 根据说明"Set maximum number of HSPs per subject sequence to save for each query", 也就是保留query中最佳的前几个匹配

继续测试,随机删除query2的一些序列得到query3,保证最后一个是最佳匹配,用来验证返回的是一个序列中最好连配,而不是位置相关结果

>query3
CTCAAGAAAAGGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGCCCCTAAACCCTAAACCCTAAACCCT    AAACCTCTGAATCCTTAATCCCTAAATCCCTAAATCTTTAAATCCTACATCCAT]GAATCCCTAAATACCTAATTCCCTAAACCCGAAACCGGTTTCTCTGGTTGAAAATCATTGTGTATATAATGATAATTTT
ATCGTTTTTATGTAATTGCTTATTGTTGTGTGTAGATTTTTTAAAAATATCATTTGAGGTCAATACAAATCCTATTTCT
TGTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTA
CATTTGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTA
TGTTTGGACATTTATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTATCAATTTATCTTTTGTGGGAAAATTATT
TAGTTGTAGGGATGAAGTCTTTCTTCGTTGTTGTTACGCTTGTCATCTCATCTCTCAATGATATGGGATGGTCCTTTAG
CATTTATTCTGAAGTTCTTCTGCTTGATGATTTTATCCTTAGCCAAAAGGATTGGTGGTTTGAAGACACATCATATCAA
AAAAGCTATCGCCTCGACGATGCTCTATTTCTATCCTTGTAGCACACATTTTGGCACTCAAAAAAGTATTTTTAGATGTCTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTC    TCTGTGAATTTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC
CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAATCTTTAAATCCTACATCCAT
GAATCCCTAAATACCTAATTCCCTAAACCCGAAACCGGTTTCTCTGGTTGAAAATCATTGTGTATATAATGATAATTTT
ATCGTTTTTATGTAATTGCTTATTGTTGTGTGTAGATTTTTTAAAAATATCATTTGAGGTCAATACAAATCCTATTTCT
TGTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTA
CATTTGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTA
TGTTTGGACATTTATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTATCAATTTATCTTTTGTGGGAAAATTATT
TAGTTGTAGGGATGAAGTCTTTCTTCGTTGTTGTTACGCTTGTCATCTCATCTCTCAATGATATGGGATGGTCCTTTAG
CATTTATTCTGAAGTTCTTCTGCTTGATGATTTTATCCTTAGCCAAAAGGATTGGTGGTTTGAAGACACATCATATCAA
AAAAGCTATCGCCTCGACGATGCTCTATTTCTATCCTTGTAGCACACATTTTGGCACTCAAAAAAGTATTTTTAGATGT
CTCAAGAAAAGCAAGGTTTTTGTTTCTCTCTCTCTCTCTGTGAATTGACTATGTTATTCTTCCATCTCCTAAAAATCTTGGTTTCTGCAGCTTACATGGC

构建索引数据库,进行检索,会找到每个序列中的最好匹配,如果你设置了-num_alignments 1,那么最后结果就真的只有一个。

#顺序1
cat query.fa query3.fa test1.fa > test5.fa
makeblastdb -in test5.fa -out test/test5 -dbtype nucl
blastn -query query.fa -db test/test5 -outfmt 6 -max_hsps 1
# 结果如下
query1  query3  100.000 100 0   0   1   100 1606    1705    8.11e-50    185
query1  query1  100.000 100 0   0   1   100 1   100 8.11e-50    185

最后的结论如下:

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

推荐阅读更多精彩内容