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是做比较基因组学乃至整个生物信息学研究所必须掌握的一种比对工具。
BLAST程序常用的两个评价指标
Score:
使用打分矩阵对匹配的片段进行打分,这是对各对氨基酸残基(或碱基)打分求和的结果,一般来说,匹配片段越长、相似性越高则Score值越大,结果越可信。
E-value:
BLAST程序在搜索空间中可随机找到获得这样高分的序列的可能性(期望值),因此E-value越高,则代表结果越有可能是随机获得的,也就越不可信。搜寻空间大小约略等于查询序列的长度乘以全部database序列长度的总和,再乘以一些系数。我们在获得一个Blast结果时需要看这两个指标。
如果Blast获得的目标序列的Score值越高并且E-value越低表明结果越可信,反之越不可信。
Blast的运行分为两个步骤:第一,建立目标序列的数据库;第二,做blast比对。
一.运行建库程序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格式
命令示例:
formatdb -i ecoli.nt -p F -o T
运行此命令就会在当前目录下产生用于BLAST搜索的7个文件,一旦如上的formatdb命令执行完毕,就不再需要ecoli.nt,可以移除。此时,blastall可以直接使用。
二.运行比对程序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"参数缺省,则结果输出方式为屏幕输出。
Blast的结果
包含的信息很丰富。每一个query的比对结果从"BLASTN"开始,记录了版本和作者信息,"Query= "之后记录了query名和序列长度。如果两条序列没有找到相关性信息,那么在"Searching.done"下方显示"* No hits found **";反之,则在"Searching.done"下方记录了该query序列和库中每一条subject序列的比对概况列表,包括比对得分(Score)和期望值(E value)。期望值是一个大于0的正实数,代表两条序列不相关的可能性。期望值是在整体上综合评定两条序列的相似性的参数,期望值数值越小,序列相似性就越高,反之期望值数值越大,相似性越低。比对的输出结果会按照期望值从低到高的顺序来排列。
Query序列和每一条subject序列比对结果的详细信息以">"开始。需要注意的是同一个query和同一个subject可能会有多个比对结果,每一个具体的结果从"Score ="开始,记录了比对得分、期望值、相似度百分比(identities)、比对的空位和两条序列的比对方向,之后是比对条形图,显示了比对区域内每个碱基的比对情况。列出两条序列的所有比对结果后,罗列比对的参数设置和统计信息,至此两条序列间的比对结果输出完毕。
对于蛋白相关的比对,需要在blastall的运行目录下放置取代矩阵,并在运行时指定此替代矩阵,程序才能正常运行,否则blastall会报错退出。一般来讲,蛋白比对时最常用的取代矩阵是BLOSUM62矩阵。
参数
仅仅运行blast的基本运行命令,得到的结果往往不能清晰准确的表示出有用的信息。最大的问题就是有太多的冗余,很多很短的比对都会出现在输出结果中,导致结果杂乱无章。
1.-e参数:
-e(value)参数是用来过滤比对较差的结果的,用"-e"参数指定一个实数,blast会过滤掉期望值大于这个数的比对结果。这样不但简化了结果,还缩短了运行时间和结果占用的空间。比如在上一个例子中,在命令行中加上限制期望值:
blastall -i query.fasta -d db.seq -o blast.out -p blastn -e 1e-10
2.-F参数:
-F(T/F)参数是用来屏蔽简单重复和低复杂度序列的。如果选"T",程序在比对过程中会屏蔽掉query中的简单重复和低复杂度序列;选"F"则不会屏蔽。缺省值为"T"。
3.-m参数:
“-e”参数能够做到筛选适当的比对结果,但是即使如此,blast的输出结果仍然非常庞大并且难以处理。为了精简输出、节省存储空间、实现更多功能并使结果易于处理,blast提供了参数“-m (integer)”来设定输出格式,可供选择的值为0~11之间的整数,缺省为0。下面就通过实例逐个解析“-m”参数能够实现的输出功能。输入文件的内容(针对-m 0到-m 7),其中:加粗的区域是三条序列的重合位置,注意subject1多一个碱基。
输出:
-m 0:缺省参数,显示一个query和一个subject两两比对的信息。
-m 1:显示query在所有subjects上的定位信息,并显示一致性比对信息,subject之间不同的碱基会被标出。
-m 2:显示query在所有subjects上的定位信息但是不显示一致性比对信息,subject之间不同的碱基会被标出。
-m 3:显示query在所有subjects的定位和一致性比对信息,不显示subjects之间的差异。
-m 4:显示query在所有subjects上的定位信息但是不显示一致性比对信息,不显示subjects之间的差异。
-m 5:显示query在所有subjects上的定位信息但是不显示每个碱基的比对信息,补充"-"对齐比对区域,subjects之间不同的碱基会被标出。
-m 6:显示query在所有subjects上的定位信息但是不显示每个碱基的比对信息,补充"-"对齐比对区域,不显示subjects之间的差异。
-m 7:输出XML格式的blast结果。
-m 8:列表格式的比对结果。从左到右各列的意义依次是:query名、subject名、identity、比对长度、错配数、空位数、query比对起始坐标、query比对终止坐标、subject比对起始坐标、subject比对终止坐标、期望值、比对得分。
-m 9:带注释行的列表格式。格式和-m 8一样,只是在每个query的比对结果前面加了注释行用以说明列表中各列的意义。
-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的数据库上。
6.-M参数:
做有关蛋白的比对时,需要用"-M"参数指定取代矩阵,比如BLOSUM45、BLOSUM62、BLOSUM80等,缺省值为BLOSUM62。这三个矩阵都可以在blast安装目录的data目录下找到。BLOSUM矩阵后面的数字代表比对结果允许的最低相似度百分比,我们可以根据不同的精度需求选择不同的取代矩阵。
7.-W参数:
-W(integer):指定做比对时的“字”的长度。缺省值是0(代表blastn的搜索字长为11,megablast是28,其他是3)。这个参数多数时候不用调整,但是需要做短序列的比对时,可能要适当调短字长,来增加比对的敏感度。
以上为blastall 的常用参数,对于一些不常用的参数,可以查找blast的参数表,此参数表可以通过直接运行blastall得到。
Blastall常用参数简析
BLAST (Basic Local Alignment Search Tool) 基本局部比对搜索工具,是一套在蛋白质数据库或DNA数据库中进行相似性比较的分析工具,它是基于Altschul等人在J.Mol.Biol上发表的方 法(J.Mol.Biol.215:403-410(1990)),在序列数据库中对查询序列进行同源性比对工作。BLAST程序 能迅速与公开数据库进行相似性序列比较,利用比较结果中的得分对序列相似性进行说明。 BLAST可以 对一条或多条序列(可以是 任何形式的序列)在一个 或多个核酸或蛋白序列库中进行比对,并且从最初的BLAST发展到现在NCBI提供的BLAST2.0,已将有缺口 的比对序列也考虑在内了。BLAST可处理任何数量的序列,包括蛋白序列和核酸序列;也可选择多个数据库但数据库必须是同一类型的,即要么都是蛋白数据库要么都是核酸数据库。所查询的序 列和调用的数据库则可以是任何形式的组合,既可以是核酸序列到蛋白库中作查询,也可以是蛋白序列到蛋白库中作查询,反之亦然。 由于Blast功能强大,检索速度快, Blast工具流行于世界上几乎所有的生物信息中心。
BLAST 提供的检索功能:
BLASTn: 核酸序列到核酸库中的一种查询。库中存在的每条已知序列都将同所查序列作一对一地核酸序列比对。 BLASTp: 蛋白序列到蛋白库中的一种查询。库中存在的每条已知序列将逐一地同每条所查序列作一对一的序列比对。
BLASTx: 核酸序列到蛋白库中的一种查询。先将核酸序列翻译成蛋白序列(一条核酸序列会被翻译成可能的6条蛋白),再对每一条作一对一的蛋白序列比对。
TBLASTn: 蛋白序列到核酸库中的一种查询。与BLASTx相反,它是将库中的核酸序列翻译成蛋白序列,再同所查序列作蛋白与蛋白的比对。
TBLASTx: 核酸序列到核酸库中的一种查询。此种查询将库中的核酸序列和所查的核酸序列都翻译成蛋白(每条核酸序列会产生6条可能的蛋白序列),这样每次比对会产生36种比对阵列。
在使用blastall对测试序列在序列数据库中进行查询之前 ,用户需要对blastall命令涉及的主要常用参数有充分的理解。
用户可以在命令行方式下运行:blastall –
下面对blastall主要常用参数进行说明:
blastall -p blastn d db.fasta -i input.fasta -o output.blast -e 1e-30 -b 2 -v 2 -m 8 -I T -a 2
-p Program Name [String]
所用程序名称[String],用户可以根据需要从blastn,blastp,blastx, tblastn,tblastx中任选一程序。
-d Database [String] default = nr
所用序列数据库的名称 [String],默认为:nr,本文例为:ecoli.nt
-i Query File [File In] default = stdin
所用查询序列文件[File In],默认为:stdin,本文例为 test.txt
-e Expectation value (E) [Real] default = 10.0
期望值[Real] 默认为10.0描述搜索某一特定数据库时,随机出现的匹配序列数目。
-m alignment view options:
比对显示选项,其具体的说明可以用以下的比对实例说明 0 = pairwise, 显示具体匹配信息(缺省)
-o BLAST report Output File [File Out] Optional default
=stdout
BLAST报告的输出文件[File Out] 默认为:stdout
-F Filter query sequence (DUST with blastn, SEG with others) [String] default = T 查询序列过滤,将那些给出影响比对结果的低 复杂度区域过滤掉。用blastn进行查询的序列用DUST程序过滤,其他的用SEG过滤。对DUST和SEG的详细情况,用户可以自己查询资料。
-G Cost to open a gap (zero invokes default behavior)
[Integer] default = 0 空位开放罚分[Integer]
(设为0则调用默认行为) 默认为0分.
-E Cost to extend a gap (zero invokes default behavior)
[Integer] default = 0
空位扩展罚分[Integer] (设为0则调用默认行为) 默认为0分
-X X dropoff value for gapped alignment (in bits) (zero invokes default behavior) blastn 30, megablast 20, tblastx 0, all others 15 [Integer],default = 0
-I Show GI's in deflines [T/F] default = F
提示行显示GI number 默认不显示
-q Penalty for a nucleotide mismatch (blastn only)
[Integer] default = -3
核酸序列基对不匹配所罚分数(blastn only) [Integer]默认罚3分
-r Reward for a nucleotide match (blastn only)
[Integer] default = 1
核苷酸序列基对匹配所加分数(blastn only) [Integer]默认加1分
-g Perfom gapped alignment (not available with tblastx) [T/F] default = T
是否执行带缺口的比对(not available with tblastx)默认为是& nbsp;
-a Number of processors to use [Integer] default = 1
使用处理器的数目[Integer] 默认为单机
-B Number of concatenated queries, for blastn and tblastn
[Integer] Optional default = 0
需要联配查询的序列数目 for blastn and tblastn
[Integer]默认为单序列
以上所列只是blastall命令部分参数的说明,用户在对自己的序列进行BLAST时可根据自己的需要选择参数,以便得到自己需要的查询报告。同时,参数选择的正确与否也是blastall程序能否顺利执行的关键。