Blast使用

简介
Blast,全称Basic Local Alignment Search Tool,即"基于局部比对算法的搜索工具",由Altschul等人于1990年发布。Blast能够实现比较两段核酸或者蛋白序列之间的同源性的功能,它能够快速的找到两段序列之间的同源序列并对比对区域进行打分以确定同源性的高低。
Blast的运行方式是先用目标序列建数据库(这种数据库称为database,里面的每一条序列称为subject),然后用待查的序列(称为query)在database中搜索,每一条query与database中的每一条subject都要进行双序列比对,从而得出全部比对结果。
Blast是一个集成的程序包,通过调用不同的比对模块,blast实现了五种可能的序列比对方式:
blastp:蛋白序列与蛋白库做比对,直接比对蛋白序列的同源性。
blastx:核酸序列对蛋白库的比对,先将核酸序列翻译成蛋白序列(根据相位可以翻译为6种可能的蛋白序列),然后再与蛋白库做比对。
blastn:核酸序列对核酸库的比对,直接比较核酸序列的同源性。
tblastn:蛋白序列对核酸库的比对,将库中的核酸翻译成蛋白序列,然后进行比对。
tblastx:核酸序列对核酸库在蛋白级别的比对,将库和待查序列都翻译成蛋白序列,然后对蛋白序列进行比对。
Blast提供了核酸和蛋白序列之间所有可能的比对方式,同时具有较快的比对速度和较高的比对精度,因此在常规双序列比对分析中应用最为广泛。可以毫不夸张的说,blast是做比较基因组学乃至整个生物信息学研究所必须掌握的一种比对工具。
下载
NCBI提供免费下载,网址:ftp://ftp.ncbi.nih.gov/blast/executables/release/,可根据自己得机器选择相应操作系统的版本。
安装
直接解压缩包即可。解压缩命令:
zcat *.tar.gz | tar xvf -
使用
Blast的运行分为两个步骤:第一,建立目标序列的数据库;第二,做blast比对。
1.运行建库程序formatdb:
建库的过程是建立目标序列的索引文件,所用程序是formatdb。程序允许的输入格式FASTA或者ASN.1格式,通常我们使用FASTA格式的序列作为输入。用于建库的FASTA序列是db.seq,formatdb的基本命令是:
formatdb -i db.seq [-options]
常用的参数有以下几个:
-p (T/F):-p参数的意义是选择建库的类型,"T"表示蛋白库,"F"表示核酸库。缺省值为"T"。
-o (T/F):-o参数的意义是判断是否分析序列名并建立序列名索引。"T"表示建立序列名索引,"F"表示不建立序列名索引。缺省值为"F"。
程序输出:
如果建立的是核酸库,输出为db.seq.nhr、db.seq.nin、db.seq.nsq,如果选择了参数"-o T",还会同时输出db.seq.nsd、db.seq.nsi、db.seq.nni、db.seq.nnd。
蛋白库和核酸库的输出类似,相应的输出文件为:db.seq.phr、db.seq.pin、db.seq.psq和db.seq.psd、db.seq.psi、db.seq.pni、db.seq.pnd。
除了这些结果,程序还会输出LOG文件(默认为formatdb.log),里面记录了运行时间、版本号、序列数量等信息。
几点需要注意的问题:
1、建库以后,做blast比对的输入文件就是建库所得的文件db.seq.n**或者db.seq.p**,而不是原始的FASTA序列。也就是说,建库以后,原始的序列文件是可以删除的。
2、如果命令行中选择了"-o T",并且目标序列中含有gi号重复的的序列名时,程序会停止建库并报错。例如,下列序列文件中出现了重复的序列名:
>gi|112385745|gb|DQ859020.1| Oryza sativa (japonica cultivar-group) glutathione S-transferase 2 mRNA, complete cds
ATGGCGGAGGCGGCGGGGGCGGCGGTGGCGCCGGCGAAGCTGGGTCTGTACTCGTACTGGCGGAGCTCGT
GCTCGCACCGCGTCCGCATCGCCCTCAACCTCAAAGGATTGGAGTACGAGTACAAGGCGGTGAACCTGCT
CAAGGGGGAGCACTCTGATCCAGAATTCATGAAGGTTAATCCTATGAAGTTCGTCCCGGCATTGGTCGAT
......
CAAGCAGCACTCCCAGACAGACAACCAGATGCCCCTTCCTCTACCTAG
>gi|112385745|gb|DQ859020.1| Oryza sativa (japonica cultivar-group) glutathione S-transferase 2 mRNA, complete cds
ATGGCGGAGGCGGCGGGGGCGGCGGTGGCGCCGGCGAAGCTGGGTCTGTACTCGTACTGGCGGAGCTCGT
GCTCGCACCGCGTCCGCATCGCCCTCAACCTCAAAGGATTGGAGTACGAGTACAAGGCGGTGAACCTGCT
CAAGGGGGAGCACTCTGATCCAGAATTCATGAAGGTTAATCCTATGAAGTTCGTCCCGGCATTGGTCGAT
......
运行时就会报如下错误:
[formatdb] ERROR: Failed to create index.  Possibly a gi included more than once in the database.
3、如果输入序列不符合FASTA格式或者ASN.1格式,程序会自动退出,并报错:
[formatdb] ERROR: Could not open db
4、核酸序列可以用于建核酸库和蛋白库,但是蛋白序列不能用于建核酸库。
其他参数简介:
-l:"-l 文件名"用来改变LOG文件的命名
-n:"-n 文件名"可以自定义生成的库文件命名
-a:输入文件为ASN.1格式
2.运行比对程序blastall:
Blast的主程序是blastall。程序的输入文件是query序列(-i 参数)和库文件(-d 参数),比对类型的选择(-p 参数)和输出文件(-o 参数)由用户指定。其中“-p”参数有5种取值:
-p blastp:蛋白序列与蛋白库做比对。
-p blastx:核酸序列对蛋白库的比对。
-p blastn:核酸序列对核酸库的比对。
-p tblastn:蛋白序列对核酸库的比对。
-p tblastx:核酸序列对核酸库在蛋白级别的比对。
这些元素就构成了blast的基本运行命令(以blastn为例):
blastall -i query.fasta -d database_prefix -o blast.out -p blastn
其中如果"-o"参数缺省,则结果输出方式为屏幕输出。下面以一个blastn比对为例,来说明比对全过程:
Query序列(query.fasta):
>gi|45593933|gb|AY551259.1| Oryza sativa precursor microRNA 319c gene
AGGAAGAGGAGCTCCTTTCGATCCAATTCAGGAGAGGAAGTGGTAGGATGCAGCTGCCGATTCATGGATA
CCTCTGGAGTGCATGGCAGCAATGCTGTAGGCCTGCACTTGCATGGGTTTGCATGACCCGGGAGATGAAC
CCACCATTGTCTTCCTCTATTGATTGGATTGAAGGGAGCTCCACATCTCT
>gi|45593932|gb|AY551258.1| Oryza sativa precursor microRNA 319b gene
CATATTCTTTTAATTTGATGGAAGAAGCGATCGATGGATGGAAGAGAGCGTCCTTCAGTCCACTCATGGG
CGGTGCTAGGGTCGAATTAGCTGCCGACTCATTCACCCACATGCCAAGCAAGAAACGCTTGAGATAGCGA
AGCTTAGCAGATGAGTGAATGAAGCGGGAGGTAACGTTCCGATCTCGCGCCGTCTTTGCTTGGACTGAAG
GGTGCTCCCTCCTCCTCGATCTCTTCGATCTAATTAAGCTACCTTGACAT
库文件Database(db.seq,已经运行formatdb -i db.seq -p F -o T建库):
>fake_seq
AGGAAGAGGAGCTCCTTTCGTTCCAATTCAGGAGAGGAAGTGGTAGGATGCAGCTGCCGATTCATGGATA
CCTCTGGAGTGCATGCAGCAATGCTGTAGGCCTGCACTTGCATGGGTTTGCATGACCCGGCGAGATGAAC
CCACCATTGTCTTCCTCTATTGATTGGATTGAAGGGAGCTCCACATCTCT
运行命令:
blastall -i query.fasta -d db.seq -o blast.out -p blastn
运行结果:
BLASTN 2.2.8 [Jan-05-2004]
Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.
Query= gi|45593933|gb|AY551259.1| Oryza sativa precursor microRNA 319c
gene, complete sequence
         (190 letters)
Database: db.seq
           1 sequences; 190 total letters
Searching.done
                                                                    Score    E
Sequences producing significant alignments:                         (bits) Value
 
fake_seq                                                             339   2e-98
>fake_seq
          Length = 190
Score =  339 bits (171), Expect = 2e-98
 Identities = 188/191 (98%), Gaps = 2/191 (1%)
 Strand = Plus / Plus
Query: 1   aggaagaggagctcctttcgatccaattcaggagaggaagtggtaggatgcagctgccga 60
           |||||||||||||||||||| |||||||||||||||||||||||||||||||||||||||
Sbjct: 1   aggaagaggagctcctttcgttccaattcaggagaggaagtggtaggatgcagctgccga 60
Query: 61  ttcatggatacctctggagtgcatggcagcaatgctgtaggcctgcacttgcatgggttt 120
           |||||||||||||||||||||||| |||||||||||||||||||||||||||||||||||
Sbjct: 61  ttcatggatacctctggagtgcat-gcagcaatgctgtaggcctgcacttgcatgggttt 119
Query: 121 gcatgacccgg-gagatgaacccaccattgtcttcctctattgattggattgaagggagc 179
           ||||||||||| ||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct: 120 gcatgacccggcgagatgaacccaccattgtcttcctctattgattggattgaagggagc 179
Query: 180 tccacatctct 190
           |||||||||||
Sbjct: 180 tccacatctct 190
  Database: db.seq
    Posted date:  Aug 28, 2006  8:14 PM
  Number of letters in database: 190
  Number of sequences in database:  1
Lambda     K      H
    1.37    0.711     1.31
Gapped
Lambda     K      H
    1.37    0.711     1.31
Matrix: blastn matrix:1 -3
Gap Penalties: Existence: 5, Extension: 2
Number of Hits to DB: 3
Number of Sequences: 1
Number of extensions: 3
Number of successful extensions: 3
Number of sequences better than 10.0: 1
Number of HSP's better than 10.0 without gapping: 1
Number of HSP's successfully gapped in prelim test: 0
Number of HSP's that attempted gapping in prelim test: 0
Number of HSP's gapped (non-prelim): 1
length of query: 190
length of database: 190
effective HSP length: 8
effective length of query: 182
effective length of database: 182
effective search space:    33124
effective search space used:    33124
T: 0
A: 0
X1: 6 (11.9 bits)
X2: 15 (29.7 bits)
S1: 12 (24.3 bits)
S2: 6 (12.4 bits)
BLASTN 2.2.8 [Jan-05-2004]
 
Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.
Query= gi|45593932|gb|AY551258.1| Oryza sativa precursor microRNA 319b
gene, complete sequence
         (260 letters)
Database: db.seq
           1 sequences; 190 total letters
Searching.done
***** No hits found ******
  Database: db.seq
    Posted date:  Aug 28, 2006  8:14 PM
  Number of letters in database: 190
  Number of sequences in database:  1
Lambda     K      H
    1.37    0.711     1.31
Gapped
Lambda     K      H
    1.37    0.711     1.31
Matrix: blastn matrix:1 -3
Gap Penalties: Existence: 5, Extension: 2
Number of Hits to DB: 0
Number of Sequences: 1
Number of extensions: 0
Number of successful extensions: 0
Number of sequences better than 10.0: 0
Number of HSP's better than 10.0 without gapping: 0
Number of HSP's successfully gapped in prelim test: 0
Number of HSP's that attempted gapping in prelim test: 0
Number of HSP's gapped (non-prelim): 0
length of query: 260
length of database: 190
effective HSP length: 8
effective length of query: 252
effective length of database: 182
effective search space:    45864
effective search space used:    45864
T: 0
A: 0
X1: 6 (11.9 bits)
X2: 15 (29.7 bits)
S1: 12 (24.3 bits)
S2: 6 (12.4 bits)
 
Blast的结果包含的信息很丰富。每一个query的比对结果从"BLASTN"开始,记录了版本和作者信息,"Query= "之后记录了query名和序列长度。如果两条序列没有找到相关性信息,那么在"Searching.done"下方显示"***** No hits found ******";反之,则在"Searching.done"下方记录了该query序列和库中每一条subject序列的比对概况列表,包括比对得分(Score)和期望值(E value)。期望值是一个大于0的正实数,代表两条序列不相关的可能性。期望值是在整体上综合评定两条序列的相似性的参数,期望值数值越小,序列相似性就越高,反之期望值数值越大,相似性越低。比对的输出结果会按照期望值从低到高的顺序来排列。
Query序列和每一条subject序列比对结果的详细信息以">"开始。需要注意的是同一个query和同一个subject可能会有多个比对结果,每一个具体的结果从"Score ="开始,记录了比对得分、期望值、相似度百分比(identities)、比对的空位和两条序列的比对方向,之后是比对条形图,显示了比对区域内每个碱基的比对情况。列出两条序列的所有比对结果后,罗列比对的参数设置和统计信息,至此两条序列间的比对结果输出完毕。
如上述结果所示, AY551259和fake_seq同为正链相比得到一个结果,期望值为2e-98,identities为98%,中间有2个空位,两条序列相似度很高。而AY551258和fake_seq没有找到同源性。
对于蛋白相关的比对,需要在blastall的运行目录下放置取代矩阵,并在运行时指定此替代矩阵,程序才能正常运行,否则blastall会报错退出。一般来讲,蛋白比对时最常用的取代矩阵是BLOSUM62矩阵。
参数
仅仅运行blast的基本运行命令,得到的结果往往不能清晰准确的表示出有用的信息。最大的问题就是有太多的冗余,很多很短的比对都会出现在输出结果中,导致结果杂乱无章。例如:
                                                                    Score    E
Sequences producing significant alignments:                         (bits) Value
Contig3421   out.ace.2                                              2367   0.0
Contig3424   out.ace.2                                               165   1e-40
Contig3423   out.ace.2                                                30   4.9
Contig3314   out.ace.2                                                30   4.9
……
>Contig3423   out.ace.2
          Length = 148505
Score = 30.2 bits (15), Expect = 4.9
 Identities = 15/15 (100%)
 Strand = Plus / Plus
Query: 571    aaagaataaaattat 585
              |||||||||||||||
Sbjct: 103697 aaagaataaaattat 103711
可以很明显的看出,query序列和Contig3423的比对结果不能表示两条序列的相关性。事实上这个比对结果只是一个偶然出现的重复。这样的结果不但会浪费大量的运算和存储资源,更给结果分析带来了沉重的负担。为了处理杂乱无章的比对结果,满足各种比对需求,blast设置了很多参数来限制比对的范围和输出的形式。以下多数结果以blastn举例,如不做特殊说明,这些参数适用于所有比对方式。
1.-e参数:
-e(value)参数是用来过滤比对较差的结果的,用"-e"参数指定一个实数,blast会过滤掉期望值大于这个数的比对结果。这样不但简化了结果,还缩短了运行时间和结果占用的空间。比如在上一个例子中,在命令行中加上限制期望值:
blastall -i query.fasta -d db.seq -o blast.out -p blastn -e 1e-10
那么结果中就会只剩下比对较好的结果:
                                                                   Score    E
Sequences producing significant alignments:                        (bits) Value
Contig3421   out.ace.2                                              2367   0.0
Contig3424   out.ace.2                                               165   1e-40
……
通常,对于不同物种间的比对,期望值设在1e-5左右即可;而对于同源性较高的物种或者同种的比对,可以适当将期望值调得更小来过滤垃圾结果。比如同一物种cDNA和染色体的比对,参数可用1e-10或更高。
2.-F参数:
-F(T/F)参数是用来屏蔽简单重复和低复杂度序列的。如果选"T",程序在比对过程中会屏蔽掉query中的简单重复和低复杂度序列;选"F"则不会屏蔽。缺省值为"T"。例如,我们将如下含有两段简单重复的序列自己和自己进行比对(重复区用小写字母表示):
>test1
TACAATAAATAAAAAAGAGCTGTCTACAGTCTTTTcgcgcgcgcgcgTTCAGAAGTAAAG
CACTATACAtttttttGTTTGTTCTTCTCAATTTAGGAAACTCAATGAACAATGAATACG
AACTATTATTACCAGTAAATACAAGTAATAC
第一次比对采用缺省参数:
blastall -i test.seq -d test.seq -o test.blast -p blastn -e 1e-5
得到的结果:
>test1
          Length = 151
Score =  186 bits (94), Expect = 1e-52
 Identities = 132/151 (87%)
 Strand = Plus / Plus
Query: 1   tacaataaataaaaaagagctgtctacagtcttttnnnnnnnnnnnnttcagaagtaaag 60
           |||||||||||||||||||||||||||||||||||            |||||||||||||
Sbjct: 1   tacaataaataaaaaagagctgtctacagtcttttcgcgcgcgcgcgttcagaagtaaag 60
Query: 61  cactatacannnnnnngtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120
           |||||||||       ||||||||||||||||||||||||||||||||||||||||||||
Sbjct: 61  cactatacatttttttgtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120
Query: 121 aactattattaccagtaaatacaagtaatac 151
           |||||||||||||||||||||||||||||||
Sbjct: 121 aactattattaccagtaaatacaagtaatac 151
……
第二次运行采用参数“-F F”:
blastall -i test.seq -d test.seq -o test.blast -p blastn -e 1e-5 -F F
得到的结果:
>test1
          Length = 151
Score =  299 bits (151), Expect = 1e-86
 Identities = 151/151 (100%)
 Strand = Plus / Plus
Query: 1   tacaataaataaaaaagagctgtctacagtcttttcgcgcgcgcgcgttcagaagtaaag 60
           ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct: 1   tacaataaataaaaaagagctgtctacagtcttttcgcgcgcgcgcgttcagaagtaaag 60
Query: 61  cactatacatttttttgtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120
           ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct: 61  cactatacatttttttgtttgttcttctcaatttaggaaactcaatgaacaatgaatacg 120
Query: 121 aactattattaccagtaaatacaagtaatac 151
           |||||||||||||||||||||||||||||||
Sbjct: 121 aactattattaccagtaaatacaagtaatac 151
比较两个结果,我们看出使用缺省参数的比对结果损失了一部分信息,得到的统计结果也出现失真,期望值和identity都没有反映出真实情况。有时较长的重复序列甚至会导致比对终止。加了"-F F"就保证了比对结果的完整性。 通常在大规模、低精度的比对中,往往用缺省参数,这样能避免程序把过多的时间浪费在无意义的简单重复上,提高运行速度;而在小规模、高精度的比对中,需要加上参数"-F F",保证比对的精确度和完整性。
3.-m参数:
“-e”参数能够做到筛选适当的比对结果,但是即使如此,blast的输出结果仍然非常庞大并且难以处理。为了精简输出、节省存储空间、实现更多功能并使结果易于处理,blast提供了参数“-m (integer)”来设定输出格式,可供选择的值为0~11之间的整数,缺省为0。下面就通过实例逐个解析“-m”参数能够实现的输出功能。
输入文件的内容(针对-m 0到-m 7),其中:加粗的区域是三条序列的重合位置,注意subject1多一个碱基。
query.fasta:
>query1
TACAATAAATAAAATAGAGCTGTCTACAGTACTTTTTCAGGAACTCCTTCAGAAGTAAAG
CACTATACAtttttttGTTTGTTCTTTTCAATTTAGGAAACTCAATGAACAATGAATACG
AACTATTATTACCAGTAAATACAAGTAATAC
database.fasta:
>subject1
TCCTTCAGAAGTAAAGCACTATACAtttttttGTTTGTTCTTTTCAATTTAGGAAACTCA
AATGAACAATGAATAC
>subject2
AATTTAGGAAACTCAATGAACAATGAATACGAACTATTATTACCAGTAAATACAAGTAAT
输出:
-m 0:缺省参数,显示一个query和一个subject两两比对的信息。
>subject1
          Length = 76
Score = 93.7 bits (47), Expect = 1e-24
 Identities = 68/76 (89%), Gaps = 1/76 (1%)
 Strand = Plus / Plus
Query: 45  tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 103
           |||||||||||||||||||||||||        ||||||||||||||||||||||||||
Sbjct: 1   tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 60
……
-m 1:显示query在所有subjects上的定位信息,并显示一致性比对信息,subject之间不同的碱基会被标出。
Sequences producing significant alignments:                        (bits) Value
subject2                                                             119   2e-32
subject1                                                              94   1e-24
QUERY 45  tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactca 104
1     1                                                ............... 15
0     1   .........................ttttttt............................ 61
                                                                     \
                                                                     |
                                                                     a
……
-m 2:显示query在所有subjects上的定位信息但是不显示一致性比对信息,subject之间不同的碱基会被标出。
Sequences producing significant alignments:                        (bits) Value
subject2                                                             119   2e-32
subject1                                                              94   1e-24
QUERY 45  tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactca 104
1     1                                                aatttaggaaactca 15
0     1   tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 61
                                                                     \
                                                                     |
                                                                     a
……
-m 3:显示query在所有subjects的定位和一致性比对信息,不显示subjects之间的差异。
Sequences producing significant alignments:                        (bits) Value
subject2                                                             119   2e-32
subject1                                                              94   1e-24
QUERY 45  tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 103
1     1                                                ..............- 14
0     1   .........................ttttttt...........................a 60
……
-m 4:显示query在所有subjects上的定位信息但是不显示一致性比对信息,不显示subjects之间的差异。
Sequences producing significant alignments:                        (bits) Value
subject2                                                             119   2e-32
subject1                                                              94   1e-24
QUERY 45  tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 103
1     1                                                aatttaggaaactc- 14
0     1   tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 60
……
-m 5:显示query在所有subjects上的定位信息但是不显示每个碱基的比对信息,补充"-"对齐比对区域,subjects之间不同的碱基会被标出。
Sequences producing significant alignments:                        (bits) Value
subject2                                                             119   2e-32
subject1                                                              94   1e-24
 
QUERY 45  tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactca 104
1     1   ---------------------------------------------aatttaggaaactca 15
0     1   tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 61
                                                                     \
                                                                     |
                                                                     a
……
-m 6:显示query在所有subjects上的定位信息但是不显示每个碱基的比对信息,补充"-"对齐比对区域,不显示subjects之间的差异。
Sequences producing significant alignments:                        (bits) Value
subject2                                                             119   2e-32
subject1                                                              94   1e-24
QUERY 45  tccttcagaagtaaagcactatacannnnnnngtttgttcttttcaatttaggaaactc- 103
1      1  ---------------------------------------------aatttaggaaactc- 14
0      1  tccttcagaagtaaagcactatacatttttttgtttgttcttttcaatttaggaaactca 60
……
-m 7:输出XML格式的blast结果。
        <Hsp_num>1</Hsp_num>
        <Hsp_bit-score>119.434</Hsp_bit-score>
        <Hsp_score>60</Hsp_score>
        <Hsp_ue>1.95644e-32</Hsp_ue>
        <Hsp_query-from>90</Hsp_query-from>
        <Hsp_query-to>149</Hsp_query-to>
        <Hsp_hit-from>1</Hsp_hit-from>
        <Hsp_hit-to>60</Hsp_hit-to>
        <Hsp_query-frame>1</Hsp_query-frame>
        <Hsp_hit-frame>1</Hsp_hit-frame>
        <Hsp_identity>60</Hsp_identity>
        <Hsp_positive>60</Hsp_positive>
        <Hsp_align-len>60</Hsp_align-len>
        <Hsp_qseq>AATTTAGGAAACTCAATGAACAATGAATACGAACTATTATTA</Hsp_qseq>
        <Hsp_hseq>AATTTAGGAAACTCAATGAACAATGAATACGAACTATTATTA</Hsp_hseq>
        <Hsp_midline>||||||||||||||||||||||||||||||||||||||||||</Hsp_midline>
……
-m 8:列表格式的比对结果。从左到右各列的意义依次是:query名、subject名、identity、比对长度、错配数、空位数、query比对起始坐标、query比对终止坐标、subject比对起始坐标、subject比对终止坐标、期望值、比对得分。
query1 sub24   91.11   45    3   1    198   241   502208  502252  2.7e-06 50.05
query1 sub21   98.68   151   2   0    532   682   1360665 1360515 1.0e-76 284.0
query1 sub21   86.17   94    12  1    198   290   479232  479139  4.8e-14 75.82
query1 sub21   87.04   54    7   0    238   291   1297867 1297920 6.9e-07 52.03
query2 sub21   99.44   892   3   2    28    918   1351055 1350165 0.0     1713.2
query2 sub21   87.58   153   17  1    343   495   1358110 1357960 2.1e-35 147.2
query2 sub21   84.11   107   16  1    699   805   1305723 1305618 4.0e-12 69.88
query2 sub21   89.58   48    5   0    519   566   1305968 1305921 6.0e-08 56.00
query2 sub14   88.24   153   16  1    343   495   145402  145252  8.7e-38 155.1
query2 sub24   88.08   151   16  1    345   495   567561  567709  1.4e-36 151.2
query2 sub24   87.80   123   14  1    686   808   563341  563220  1.9e-26 117.5
 
在m8格式中通过subject的比对起止位置可以判断出序列的比对方向。比如上述结果中第1行,subject的起始坐标小于终止坐标,则两条序列是同方向比对上的;第2行中subject起始坐标大于终止坐标,则query序列是和subject的互补链比上的。
-m 9:带注释行的列表格式。格式和-m 8一样,只是在每个query的比对结果前面加了注释行用以说明列表中各列的意义。
# BLASTN 2.2.8 [Jan-05-2004]
# Query: query1   out.ace.1
# Database: database.seq
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
query1 sub24   91.11   45    3   1    198   241   502208  502252  2.7e-06 50.05
query1 sub21   98.68   151   2   0    532   682   1360665 1360515 1.0e-76 284.0
query1 sub21   86.17   94    12  1    198   290   479232  479139  4.8e-14 75.82
query1 sub21   87.04   54    7   0    238   291   1297867 1297920 6.9e-07 52.03
# BLASTN 2.2.8 [Jan-05-2004]
# Query: query1   out.ace.1
# Database: database.seq
# Fields: Query id, Subject id, % identity, alignment length, mismatches, gap openings, q. start, q. end, s. start, s. end, e-value, bit score
query2 sub21   99.44   892   3   2    28    918   1351055 1350165 0.0     1713.2
query2 sub21   87.58   153   17  1    343   495   1358110 1357960 2.1e-35 147.2
query2 sub21   84.11   107   16  1    699   805   1305723 1305618 4.0e-12 69.88
query2 sub21   89.58   48    5   0    519   566   1305968 1305921 6.0e-08 56.00
query2 sub14   88.24   153   16  1    343   495   145402  145252  8.7e-38 155.1
query2 sub24   88.08   151   16  1    345   495   567561  567709  1.4e-36 151.2
query2 sub24   87.80   123   14  1    686   808   563341  563220  1.9e-26 117.5
-m 10和11:分别是ASN格式的文本文件和二进制文件,这里就不做介绍了。
“-m”参数的值从1到6都是为了便于在subjects之间做比较而设立的功能;8和9保留了所有比对结果的原貌,只是统计成了列表的格式,从而大幅度降低了存储空间的消耗,并使结果更加清晰易读。但是m8/m9格式也有相应的缺点,就是损失了一部分比对信息,除了序列长度信息和比对条形图以外,还会在blastx、tblastn和tblastx的比对中损失关键的相位信息,这是要尽量避免的。因此在大规模的blastn比对任务中,往往要采用m8格式的输出结果来节省空间;而在小规模高精度比对中,通常用默认的输出格式,再用其他程序来提取结果中的有用信息。
4.-v参数和-b参数:
这两个参数都是限制输出结果的数量的。
-v (integer):规定输出中每一个query的比对列表最多显示subject个数(即"Sequences producing significant alignments:"后面列出的subjects数目),缺省为500条。
-b (integer):规定输出中每个query最多显示与多少条subject的比对条形图(即每条query的结果中">"的个数),缺省为250条。
如果同时使用"-m 8"参数,则输出结果中的subjects数量和"-b"参数规定的数量一致。
在database数据中能和query比上的subjects过多的时候,这两个参数就能够帮助我们把其中比对结果最好的一部分挑出来,屏蔽掉相对差的结果。当然有些时候我们是不希望屏蔽掉这些结果的,比如在某个大基因组的Contig数据集中统计一条转座子的重复次数,就需要把"-v"和"-b"参数定的足够大以显示所有结果。
5.-T参数:
-T (T/F)参数用于决定是否输出html格式的比对结果,缺省值为"F"。选择"-T T"就会输出html格式的比对结果。如果在建库过程中选择了"-o T",并且database数据中的序列是以gi号命名的,那么在html结果中以gi号命名的相应序列会自动链接到NCBI的数据库上。如图3-14:
 
图3-16 html格式的blast结果
6.-M参数:
做有关蛋白的比对时,需要用"-M"参数指定取代矩阵,比如BLOSUM45、BLOSUM62、BLOSUM80等,缺省值为BLOSUM62。这三个矩阵都可以在blast安装目录的data目录下找到。BLOSUM矩阵后面的数字代表比对结果允许的最低相似度百分比,我们可以根据不同的精度需求选择不同的取代矩阵。BLOSUM62的内容如下:
#  Matrix made by matblas from blosum62.iij
#  * column uses minimum score
#  BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
#  Blocks Database = /data/blocks_5.0/blocks.dat
#  Cluster Percentage: >= 62
#  Entropy =   0.6979, Expected =  -0.5209
   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  *
A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1  0 -4
R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1  0 -1 -4
N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  3  0 -1 -4
D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4  1 -1 -4
C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4
Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0  3 -1 -4
E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4
G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -2 -1 -4
H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0  0 -1 -4
I -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3 -3 -3 -1 -4
L -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1 -4 -3 -1 -4
K -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2  0  1 -1 -4
M -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1 -3 -1 -1 -4
F -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1 -3 -3 -1 -4
P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2 -2 -1 -2 -4
S  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2  0  0  0 -4
T  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0 -1 -1  0 -4
W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3 -4 -3 -2 -4
Y -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1 -3 -2 -1 -4
V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4 -3 -2 -1 -4
B -2 -1  3  4 -3  0  1 -1  0 -3 -4  0 -3 -3 -2  0 -1 -4 -3 -3  4  1 -1 -4
Z -1  0  0  1 -3  3  4 -2  0 -3 -3  1 -1 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4
X  0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2  0  0 -2 -1 -1 -1 -1 -1 -4
* -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4  1
 
7.-W参数:
-W(integer):指定做比对时的“字”的长度。缺省值是0(代表blastn的搜索字长为11,megablast是28,其他是3)。这个参数多数时候不用调整,但是需要做短序列的比对时,可能要适当调短字长,来增加比对的敏感度。
以上为blastall 的常用参数,对于一些不常用的参数,可以查找blast的参数表,此参数表可以通过直接运行blastall得到。
常见问题
1.运行formatdb时报下列错误:
[formatdb] ERROR: 1.seq.nhrOutput
Blast-def-line-set.E.<title>
Invalid value(s) [9] in VisibleString [seq1#ss ...]
这种情况通常是subject序列的标题行中含有非法字符,比如Tab。
 
2.运行blastall时报下列警告:
[blastall] WARNING: seq1: Could not find index files for database sub.seq
出现这种问题通常是没有做formatdb建库或者建的库不符合比对类型要求,需要重新建库。尤其注意tblastx的建库应该使用核酸库。
 
3.运行blastall时报下列警告和错误:
[blastall] WARNING:  [000.000]  seq1: SetUpBlastSearch failed.
[blastall] ERROR:  [000.000]  seq1: BLASTSetUpSearch: Unable to calculate Karlin-Altschul params, check query sequence
[blastall] ERROR:  [000.000]  seq1: BLASTSetUpSearch: Unable to calculate Karlin-Altschul params, check query sequence
 
这种情况多出在reads等短序列的比对中,某一个query序列的有效长度是0,则会导致这个错误。但是这个错误不会影响到其他query序列的正常比对,通常情况下可以忽略。
 
4.运行blastall时报下列警告并退出:
[blastall] WARNING:  [000.000]  seq1: Unable to open BLOSUM62
[blastall] WARNING:  [000.000]  seq1: BlastScoreBlkMatFill returned non-zero status
[blastall] WARNING:  [000.000]  seq1: SetUpBlastSearch failed.
这是因为做蛋白相关比对的时候目录下没有蛋白比对矩阵BLOSUM62。
实例
其他四种比对方式的使用和blastn大同小异,下面通过几个例子加以概括:
1.用BLOSUM45矩阵、“-F F”参数对两条相似的蛋白做blastp比对:
输入的query序列:query.seq
输入的database序列:db.seq
运行命令:blastall -i db.seq -d db.seq -o blastp.out -p blastp -F F -M BLOSUM45
输出结果存放在:blastp.out
2.用一条核酸序列和一条蛋白序列做交叉比对blastx和tblastn:
输入的核酸序列:cdna.seq
输入的蛋白序列:pep.seq
运行命令:
blastall -i cdna.seq -d pep.seq -o blastx.out -p blastx
blastall -i pep.seq -d cdna.seq -o tblastn.out -p tblastn
输出结果存放在:
blastx.out
tblastn.out
3.用2条核酸序列蛋白比对tblastx:
输入的query序列:query.seq
输入的database序列:db.seq
运行命令:blastall -i query.seq -d db.seq -o tblastx.out -p tblastx -e 0.5
输出结果存放在:tblastx.out
 
练习
1.用blast搜索人类基因组Chr10的重复序列(提示,为了减少比对时间,可以采取分割比对的方法)。
2.用blast查找蛋白序列ENSP00000328808在人类基因组Chr10上的位置。
3.用blast检查以下引物primer1和primer2是否可以用作在参考序列refseq上扩增。(提示,注意引物可能发生错配的条件和blast的比对参数)
>refseq
CTTAATTCGCCTCGTGAAAGAATATCATCTGCTGAACCCGGTCATTGTTG
ACTGCACCTCCAGCCAGGCAGTGGCGGATCAATATGCCGACTTCCTGCGC
GAAGGTTTCCACGTTGTCACGCCGACGGCACCTCCAGCCAGGCAGTCGAT
GGATTACTACCATCTGTTGCGTCATGCGGCTGAAAAATCGCGGCGTAAAT
TCCTCTATGACACCAACGTTGGGGCTGGATTACCGGTTATTGAGAACCTG
CAAAATCTGCTCAATGCTGGTGATGAATTGATGAAGTTCTCCGGCATTCT
TTCAGGTTCGCTTTCTTATATCTTCGGCAAGTTAGACGAAGGCATGAGTT
 
>primer1
ACTGCACCTCCAGCCAGGCAG
 
>primer2
GATATAAGAAAGCGAACCTG
 
参考文献
1.   Altschul, S.F., Gish, W., Miller, W., Myers, E.W., and Lipman, D.J. 1990. Basic local alignment search tool. J. Mol. Biol.215: 403–410.
2.   Altschul, S.F., Madden, T.L., Schaffer, A.A., Zhang, J., Zhang, Z.,Miller, W., and Lipman, D.J. 1997. Gapped BLAST and PSI-BLAST: A new generation of protein database search programs. Nucleic Acids Res. 25: 3389–3402.
 
From : BGI-生物信息学培训教材
  • 3
    点赞
  • 22
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
在Linux系统下安装和使用BLAST(Basic Local Alignment Search Tool)可以按照以下步骤进行: 1. 安装依赖项:首先,确保您的系统已经安装了必要的依赖项,包括gcc编译器、make工具和zlib库。您可以使用包管理器(如apt、yum或dnf)来安装这些依赖项。例如,在Ubuntu系统上,可以使用以下命令进行安装: ```shell sudo apt update sudo apt install build-essential zlib1g-dev ``` 2. 下载BLAST软件包:访问NCBI BLAST官方网站(https://blast.ncbi.nlm.nih.gov/Blast.cgi)并下载适用于Linux系统的BLAST软件包。您可以选择下载预编译的二进制文件或源代码。 3. 解压缩文件:如果您下载的是预编译的二进制文件,解压缩下载的文件。如果您下载的是源代码,则需要解压缩并进入解压后的目录。 4. 配置和编译:打开终端,进入解压后的BLAST目录,并运行以下命令来配置和编译BLAST: ```shell ./configure make ``` 5. 安装:运行以下命令以root权限安装BLAST: ```shell sudo make install ``` 6. 设置环境变量:为了能够在终端中随时使用BLAST命令,您需要将BLAST可执行文件所在的路径添加到系统的PATH环境变量中。您可以编辑~/.bashrc文件并在其中添加以下行(假设您安装的是blast+软件包): ```shell export PATH="/path/to/blast+/bin:$PATH" ``` 替换"/path/to/blast+"为您实际安装的BLAST软件包的路径。 7. 验证安装:重新启动终端或运行`source ~/.bashrc`使环境变量生效。然后,尝试运行以下命令来验证BLAST是否成功安装: ```shell blastn -version ``` 如果成功安装,将显示BLAST软件的版本信息。 现在您已经成功安装了BLAST。您可以使用各种BLAST命令(如blastn、blastp等)来执行不同类型的序列比对和搜索操作。请参考BLAST官方文档以了解更多关于使用BLAST的详细信息和示例。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值