登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

云之南

风声,雨声,读书声,声声入耳;家事,国事,天下事,事事关心

 
 
 

日志

 
 
关于我

专业背景:计算机科学 研究方向与兴趣: JavaEE-Web软件开发, 生物信息学, 数据挖掘与机器学习, 智能信息系统 目前工作: 基因组, 转录组, NGS高通量数据分析, 生物数据挖掘, 植物系统发育和比较进化基因组学

Blast使用方法攻略(一)  

2010-01-01 12:06:51|  分类: 生物信息学 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

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         H

    1.37    0.711     1.31

Gapped

Lambda         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         H

    1.37    0.711     1.31

Gapped

Lambda         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矩阵

  评论这张
 
阅读(1471)| 评论(0)

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018