blast几乎是无人不知的序列比对的工具,我说的是对学生物的同学,特别是搞bioinformatics。该工具一种是在线的(http://blast.ncbi.nlm.nih.gov/Blast.cgi),一种是本地的(Local BLAST,本地版blast)

        假设你有几条序列,想要看看这条序列是不是在某个物种中有类似的序列,或者你想在整个蛋白或者核酸数据库中搜索对应的同源序列,一般采用在线BLAST,使用ncbi的服务器,很快捷,界面也很美观。对于大量序列或者是自定义的数据库,一般采用本地BLAST,效率很更高。想想成千上万次手动提交序列,简直是不可能的。而且短时间内不停对服务器发起链接说不定会被ncbi ban掉(我就被ban了好几次,不过不是因为blast)。或者你想使用自己独特的数据库,不想使用什么NR,NT等等数据库的话,你就需要构建自己本地的数据库了

一、Blast简介

相似的序列(similarity)可以预测同源(homology),同源的序列来预测功能。这就是blast的根本目的。

但这里,我们需要先认清两个基本的概念:

1.相似性和同源性

如上所述,数据库搜索的基础是序列的相似性比对,而寻找同源序列则是数据库搜索的主要目的之一。所谓同源序列,简单地说,是指从某一共同祖先经趋异进化而形成的不同序列。必须指出,相似性(similarity)和同源性(homology)是两个完全不同的概念。相似性是指序列比对过程中用来描述检测序列和目标序列之间相同DNA碱基或氨基酸残基顺序所占比例的高低。当相似程度高于50%时,比较容易推测检测序列和目标序列可能是同源序列;而当相似性程度低于20%时,就难以确定或者根本无法确定其是否具有同源性。总之,不能把相似性和同源性混为一谈。所谓“具有50%同源性”,或“这些序列高度同源”等说法,都是不确切的,应该避免使用。

相似性概念的含义比较广泛,除了上面提到的两个序列之间相同碱基或残基所占比例外,在蛋白质序列比对中,有时也指两个残基是否具有相似的特性,如侧链基团的大小、电荷性、亲疏水性等。在序列比对中经常需要使用的氨基酸残基相似性分数矩阵,也使用了相似性这一概念。此外,相似性概念还常常用于蛋白质空间结构和折叠方式的比较。

2.局部相似性和整体相似性

序列比对的基本思想,是找出检测序列和目标序列的相似性。比对过程中需要在检测序列或目标序列中引入空位,以表示插入或删除(图3.1)。序列比对的最终实现,必须依赖于某个数学模型。不同的模型,可以从不同角度反映序列的特性,如结构、功能、进化关系等。很难断定,一个模型一定比另一个模型好,也不能说某个比对结果一定正确或一定错误,而只能说它们从某个角度反映了序列的生物学特性。此外,模型参数的不同,也可能导致比对结果的不同。

序列比对的数学模型大体可以分为两类,一类从全长序列出发,考虑序列的整体相似性,即整体比对;第二类考虑序列部分区域的相似性,即局部比对。局部相似性比对的生物学基础是蛋白质功能位点往往是由较短的序列片段组成的,这些部位的序列具有相当大的保守性,尽管在序列的其它部位可能有插入、删除或突变。此时,局部相似性比对往往比整体比对具有更高的灵敏度,其结果更具生物学意义。

区分这两类相似性和这两种不同的比对方法,对于正确选择比对方法是十分重要的。应该指出,在实际应用中,用整体比对方法企图找出只有局部相似性的两个序列之间的关系,显然是徒劳的;而用局部比对得到的结果也不能说明这两个序列的三维结构或折叠方式一定相同。BLAST和FastA等常用的数据库搜索程序均采用局部相似性比对的方法,具有较快的运行速度,而基于整体相似性比对的数据库搜索程序则需要超级计算机或专用计算机才能实现。

 

 

二、程序的使用

1  安装程序:

在ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/下载
ncbi-blast-2.2.28+-x64-linux.tar.gz
解压到用户的主目录(/sam)下(该目录位置你可以自己设定,跟后面统一即可),把解压后的文件夹重新命名为blast,则BLAST+的所有程序在目录/sam/blast/bin下。

添加环境变量:

打开终端(Terminal),切换为root用户,执行sudo vi /etc/profile
 在最末尾添加export PATH=”/sam/blast/bin:$PATH”,保存退出。
 或直接找到/etc/profile这个文件,在最末尾添加export PATH=”/sam/blast/bin:$PATH”
 此处若成功,在终端则执行blastn -version会出现版本信息。

 在目录/sam/blast下新建一个文件夹,命名为db
 在/sam下新建一个文件,命名为.ncbirc
 在文件中添加内容
 [BLAST]
 BLASTDB=/sam/blast/db
 (我后来找不到路径的时候我在/home/sam也同样建立了这样的一个文件,如果你不是在/home/sam下建立数据库,你得同时在这个地方建立./nvbirc的文件)
 主要是blast识别~/.ncbirc文件,从该文件处得知数据库的位置,这个文件很重要!!!

2 下载FASTA格式的数据库:

2.1自己的数据库

A.  如果是要做自己的数据库,要有自己的fasta文件,使用blast软件中所带的程序,将fasta文件 转变为blast的数据库文件:
B.  /db下的数据库更小,下载更快;可以生成FASTA格式;在/sam/blast/bin/update_blastdb.pl可以用来升级数据库;不用formatdb数据库;各种ID可以进入这个数据库

2.2公共数据库

下载地址ftp://ftp.ncbi.nlm.nih.gov/blast/db/
里面可以看到,有很多个分开的文件,但是其中FASTA也有很多个fasta格式的数据库,是一个完整的包,这个文件和上一节文件夹里面的数据有啥区别呢

Contents of the /blast/db/ directory
 The pre-formatted BLAST databases are archived in this directory. The
 name of these databases and their contents are listed below.
 +----------------------+-----------------------------------------------+
 |File Name | Content Description |
 +----------------------+-----------------------------------------------+
 /FASTA | subdirectory for FASTA formatted sequences
 README | README for this subdirectory (this file)
 env_nr.*tar.gz | Environmental protein sequences
 env_nt.*tar.gz | Environmental nucleotide sequences
 est.*tar.gz | volumes of the formatted est database| from the EST division of GenBank, EMBL, | and DDBJ
 est_human.tar.gz | alias and mask files for human subset of the est
 est_mouse.tar.gz | alias and mask files for mouse subset of the est
 est_others.tar.gz | alias and mask files for non-human and non-mouse | subset of the est database | These alias and mask files need all volumes of | est to function properly.
 gss.*tar.gz | volumes of the formatted gss database| from the GSS division of GenBank, EMBL, and | DDBJ
 htgs.*tar.gz | volumes of htgs database with entries | from HTG division of GenBank, EMBL, and DDBJ
 human_genomic.*tar.gz | human RefSeq (NC_######) chromosome records| with gap adjusted concatenated NT_ contigs
 nr.*tar.gz | non-redundant protein sequence database with | entries from GenPept, Swissprot, PIR, PDF, PDB, | and NCBI RefSeq
 nt.*tar.gz | nucleotide sequence database, with entries| from all traditional divisions of GenBank, | EMBL, and DDBJ excluding bulk divisions (gss, | sts, pat, est, and htg divisions. wgs entries | are also excluded. Not non-redundant.
 other_genomic.*tar.gz | RefSeq chromosome records (NC_######) for | organisms other than human
 pataa.*tar.gz | patent protein sequence database
 patnt.*tar.gz | patent nucleotide sequence database| The above two databases are directly from | USPTO or from EU/Japan Patent Agencies via | EMBL/DDBJ pdbaa.*tar.gz | protein sequences from pdb protein structures, | its parent database is nr.
 pdbnt.*tar.gz | nucleotide sequences from pdb nucleic acid | structures, its parent database it nt. They are | NOT the protein coding sequences for the | corresponding pdbaa entries.
 refseq_genomic.*tar.gz | NCBI genomic reference sequences
 refseq_protein.*tar.gz | NCBI protein reference sequences
 refseq_rna.*tar.gz | NCBI Transcript reference sequences
 sts.*tar.gz | Sequences from the STS division of GenBank, EMBL, | and DDBJ
 swissprot.tar.gz | swiss-prot sequence databases (last major update), | its parent database is nr.
 taxdb.tar.gz | Additional taxonomy information for the formatted | database (contains common and scientific names)
 wgs.*tar.gz | volumes for whole genome shotgun sequence assemblies | for different organisms
 +----------------------+-----------------------------------------------+
Contents of the /blast/db/FASTA directory
 This directory contains FASTA formatted sequence files. The file names
 and database contents are listed below. These files are now archived
 in .gz format and must be processed through formatdb before they can be
 used by the BLAST programs.
 +-----------------------+-----------------------------------------------+
 |File Name | Content Description |
 +-----------------------+-----------------------------------------------+
 alu.a.gz | translation of alu.n repeats
 alu.n.gz | alu repeat elements
 drosoph.aa.gz | CDS translations from drosophila.nt
 drosoph.nt.gz | genomic sequences for drosophila
 env_nr.gz* | Environmental protein sequences
 env_nt.gz* | Environmental nucleotide sequences
 est_human.gz* | human subset of the est database (see Note 1)
 est_mouse.gz* | mouse subset of the est database
 est_others.gz* | non-human and non-mouse subset of the est
 database
 gss.gz* | sequences from the GSS division of GenBank, | EMBL, and DDBJ
 htg.gz* | htgs database with high throughput genomic | entries from the htg division of GenBank, | EMBL, and DDBJ
 human_genomic.gz* | human RefSeq (NC_######) chromosome records | with gap adjusted concatenated NT_ contigs
 igSeqNt.gz | human and mouse immunoglobulin variable region nucleotide | sequences
 igSeqProt.gz | human and mouse immunoglobulin variable region protein | sequences
 mito.aa.gz | CDS translations of complete mitochondrial | genomes
 mito.nt.gz | complete mitochondrial genomes
 month.aa.gz | newly released/updated protein sequences
 (See Note 2)
 month.est_human.gz | newly released/updated human est sequences
 month.est_mouse.gz | newly released/updated mouse est sequences
 month.est_others.gz | newly released/updated est other than | human/mouse
 month.gss.gz | newly released/updated gss sequences
 month.htgs.gz | newly released/updated htgs sequences
 month.nt.gz | newly released/updated sequences for the nt
 database
 nr.gz* | non-redundant protein sequence database with| entries from GenPept, Swissprot, PIR, PDF, | PDB, and RefSeq
 nt.gz* | nucleotide sequence database, with entries | from all traditional divisions of GenBank, | EMBL, and DDBJ excluding bulk divisions| (gss, sts, pat, est, htg divisions) and wgs | entries. Not non-redundant.
 other_genomic.gz* | RefSeq chromosome records (NC_######) for | organisms other than human
 pataa.gz* | patent protein sequence database
 patnt.gz* | patent nucleotide sequence database | The above two dbs are directly from USPTO | of from EU/Japan Patent Agency via EMBL/DDBJ
 pdbaa.gz* | protein sequences from pdb protein structures
 pdbnt.gz* | nucleotide sequences from pdb nucleic acid | structures. They are NOT the protein coding | sequences for the corresponding pdbaa entries.
 sts.gz* | database for sequence tag site entries
 swissprot.gz* | swiss-prot database (last major release)
 vector.gz | vector sequence database (See Note 3)
 wgs.gz* | whole genome shotgun genome assemblies
 yeast.aa.gz | protein translations from yeast genome
 yeast.nt.gz | yeast genomes.
 +-----------------------+-----------------------------------------------+
FASTA格式的数据库,需要写入,便于搜索
 cd /sam/balst/db
 formatdb -i input_db -p F -o T for nucleotide
 formatdb -i input_db -p T -o T for protein
 -i 输入需要格式化的源数据库名称 Optional
 -p 文件类型,是核苷酸序列数据库,还是蛋白质序列数据库
 T – protein F - nucleotide [T/F] Optional default = T
 -a 输入数据库的格式是ASN.1(否 则是FASTA)
 T - True, F - False. [T/F] Optional default = F
 -o 解析选项
 T - True: 解析序列标识并且建立目录 F - False: 与上相反
 [T/F] Optional default = F
 注:1 因为我选择的是/db下有许多个文件夹组成的某一个数据库,所以就没有写入这个过程。
 2 在/db下下载的很多个数据库,比如nr00,nr01,nr02……解压缩后要放到一个文件夹里面,用数据库的时候,blast -db /sam/blast/db/nr/nr.我在这个问题上费了很久,刚开始是选择什么样的数据库,后来是下载的nr分开的数据库如何整合到一起。其实运行的是那个唯一不带编号的Nr文件。
 3 以下的命令是网上找到的,针对FASTA格式的,但是我没用过,把它放在这上面,供大家参考吧
 建立BLAST+可用的数据库:
 打开终端(Terminal),切换到/sam/blast/db目录下,执行:
 makeblastdb –in nr -parse_seqids -hash_index -dbtype prot
 或者是对某个fasta文件进行处理:其中的makeblastdb在bin文件夹里面。
 makeblastdb -in db.fasta -dbtype prot -parse_seqids -out dbname
 参数说明:
 -in:待格式化的序列文件
 -dbtype:数据库类型,prot或nucl
 -out:数据库名
 -parse_seqids 根据序列ID来分裂
 具体参数的设置可以参见makeblastdb -help

3选择程序:

(1) blastp:用蛋白质序列搜索蛋白质序列库

(2)?blastn:用核酸序列搜索核酸库

(3)blastx:核酸序列对蛋白质库的比对,核酸序列在比对之前自动按照六个读码框翻译成蛋白质序列

(4)tblastn:蛋白质序列对核酸库的比对,核酸库中的序列按照六个读码框翻译后与蛋白质序列进行比对搜索

(5)tblastx:核酸序列对核酸库在蛋白质质级别的比对,两者都在搜索之前翻译成为蛋白质质进行比对

(为什么会是6呢,因为DNA在翻译为蛋白质的时候移位,一个aa由三个碱基构成,所以可以移动3个,同时又是因为正反链,所以会有6种情况。)

氨基酸替代矩阵
PAM=Point Accepted Mutations 相近的序列比对得出的结论
BLOSUM = Blocks Amino Acid Substitution Matrices,来自BLOCKS database (http://blocks.fhcrc.org/)

The numbers mean different things depending on the particular matrix. For the

PAM matrix, the high number matrices are better for more divergent alignments, while the

opposite is true for the BLOSUM matrices.

假如我们找距离比较远的序列,我们可以用PAM更大数目的举证或者是BLOSUM更小的矩阵

Position-Specific Iterated (PSI)-BLAST:

很灵敏的BLAST工具,可以用来找距离很远的蛋白或者找蛋白家族新的成员。就是更加灵敏呗,第一步的psiblast其实就是简单的blast,构建新的矩阵,然后在对于大于阈值重新blast,比对的结果中有用黄色标记的就是之前没有找序列。

4.使用程序:

/sam/blast/bin/blastp -query assembly.orfs.hmm.faa -db /sam/blast/db/refseq_protein/refseq_protein -evalue 1e-5 -num_threads 10 -max_target_seqs 5 -outfmt 5 -out assembly.orfs.hmm.blast.xml
 -query 输入序列
 -db 选择的数据库
 -evalue 这个就不说了
 -num_threads 服务器使用的cpu个数
 -max_target_seqs设置最多的目标序列匹配数
 -out输出的文件名
 -remote 如果加上这个参数,就是使用远程比对,再有互联网的情况下,连接到NCBI的服务器来进行远程比对,然后返回到本地,从而不消耗本地的资源,设置了该参数,则-num_threads无效
 -perc_identity 相似度,默认是大于0就显示,如果要显示大于95%,这后面接95
 -outfmt
 5代表xml格式,这个格式的信息最全,6代表table格式,该格式的结果以table的形式给出,清晰易懂,7代表有注释行的table格式,比6多了一些#开头的注释行。
 注意文件或程序的路径,如果没有cd到相应的位置的话,就需要在命令中注明路径

如果你对比对后的结果有自己的要求的话,建议xml用5这个格式的输出,然后用脚本来提取,我之前用python写过一个这样的脚本:blast结果的提取

“7 qacc sacc evalue length pident” :这个是新BLAST+中最拉风的功能了,直接控制输出格式,不用再用parser啦, 7表示带注释行的tab格式的输出,可以自定义要输出哪些内容,用空格分格跟在7的后面,并把所有的输出控制用双引号括起来,其中qacc查询序列的acc,sacc表示目标序列的acc,evalue即是e值,length即是匹配的长度,pident即是序列相同的百分比,其他可用的特征(红色字体)如下:

*** Formatting options
-outfmt
 alignment view options:
 0 = pairwise,
 1 = query-anchored showing identities,
 2 = query-anchored no identities,
 3 = flat query-anchored, show identities,
 4 = flat query-anchored, no identities,
 5 = XML Blast output,
 6 = tabular,
 7 = tabular with comment lines,
 8 = Text ASN.1,
 9 = Binary ASN.1
 10 = Comma-separated values
 Options 6, 7, and 10 can be additionally configured to produce
 a custom format specified by space delimited format specifiers.
 The supported format specifiers are:
 qseqid means Query Seq-id
 qgi means Query GI
 qacc means Query accesion
 sseqid means Subject Seq-id
 sallseqid means All subject Seq-id(s), separated by a ';'
 sgi means Subject GI
 sallgi means All subject GIs
 sacc means Subject accession
 sallacc means All subject accessions
 qstart means Start of alignment in query
 qend means End of alignment in query
 sstart means Start of alignment in subject
 send means End of alignment in subject
 qseq means Aligned part of query sequence
 sseq means Aligned part of subject sequence
 evalue means Expect value
 bitscore means Bit score
 score means Raw score
 length means Alignment length
 pident means Percentage of identical matches
 nident means Number of identical matches
 mismatch means Number of mismatches
 positive means Number of positive-scoring matches
 gapopen means Number of gap openings
 gaps means Total number of gaps
 ppos means Percentage of positive-scoring matches
 frames means Query and subject frames separated by a '/'
 qframe means Query frame
 sframe means Subject frame
 When not provided, the default value is:
 'qseqid sseqid pident length mismatch gapopen qstart qend sstart send
 evalue bitscore', which is equivalent to the keyword 'std'
 Default = `0'
调用blastn合作加-help参数可以打印出下面详细的帮助信息
blastn -help。

4 手动注释–tblastn

我们从宏基因组中聚类出我们感兴趣的微生物物种的基因组后,想要预测其潜在的新的代谢途径的时候,一般会有两种策略:
1,对该基因组进行注释,获得所有的注释结果,然后对注释结果的解读;
2,通过看文献,搜集一些感兴趣的相关代谢途径的关键基因,把这些关键基因的蛋白序列给下下来,然后在你的基因组里面找看有没有相关的基因。第一种策略的注释之前的博文零星的有一些介绍,下面主要讲一下第二种策略所涉及到的一个方法—-tblastn.

$ cd /sam/uncltured/contig7/c3000/syn/genome2/tblastn
 $ makeblastdb -in genome2.fasta -dbtype nucl -title syn1 -parse_seqids -out syn1 -logfile syn1.log #这一步是将的基因组作为数据库
 $ /sam/blast/bin/tblastn -query /sam/syn/alkylsuccinateSynthase_α-subunit_assA.fasta -db syn1 -evalue 1e-5 -num_threads 60 -outfmt 5 -out assA.xml #比对,得到的是一个xml文件
 $ /sam/blast/bin/tblastn -query /sam/syn/alkylsuccinateSynthase_α-subunit_assA.fasta -db syn1 -evalue 1e-5 -num_threads 60 -outfmt 7 -out assA #比对,得到的是一个table文件
 其实-outfmt 5得到的结果就是一个比较详尽的结果,-outfmt 7得到的是一个table文件。
 # TBLASTN 2.2.28+
 # Query: gi|299800799|gb|ADJ51097.1| alkylsuccinate synthase [Desulfoglaeba alkanexedens ALDC]
 # Database: syn1
 # Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
 # 9 hits found
 gi|299800799|gb|ADJ51097.1| 109 69.50 341 104 0 365 705 4 1026 5e-169 494
 gi|299800799|gb|ADJ51097.1| 226 68.48 257 81 0 449 705 773 3 4e-118 395
 gi|299800799|gb|ADJ51097.1| 226 73.17 82 22 0 290 371 1007 762 2e-32 132
 gi|299800799|gb|ADJ51097.1| 839 31.61 794 504 19 77 858 13150 10850 7e-91 316
 gi|299800799|gb|ADJ51097.1| 292 87.23 47 6 0 812 858 1 141 2e-20 90.1
 gi|299800799|gb|ADJ51097.1| 900 87.23 47 6 0 812 858 111138 110998 2e-19 90.5
 gi|299800799|gb|ADJ51097.1| 201 48.44 64 33 0 36 99 192 1 9e-13 68.9
 gi|299800799|gb|ADJ51097.1| 322 46.88 64 34 0 36 99 192 1 2e-12 67.4
 gi|299800799|gb|ADJ51097.1| 1072 60.00 35 14 0 824 858 7642 7746 2e-06 48.1
 # BLAST processed 1 queries
 #当然输入序列也可以使好多个蛋白序列合并的fasta文件,例如,我对alkylsuccinateSynthase感兴趣,我可以到ncbi中把相关的基因的蛋白序列都给下下来,放到一个文件里面。

这里面有两个问题:
1,我在下载蛋白序列的时候,发现有的序列长度为800个aa,有的仅为200个aa的parital,为什么partial也能放到数据库中呢?
2,比对的结果,如何设置阈值,判断结果中跟输入序列有很高的同源性呢?

参考资料:

ywpeng的个人博客 http://blog.sciencenet.cn/blog-455820-609063.html
nci官网: ftp://ftp.ncbi.nlm.nih.gov/blast/db/README
http://www.ncbi.nlm.nih.gov/books/NBK52640/
lidaof的空间

8 thoughts on “blast简介及参数设置

  1. 它是按哪个关键字索引的,为何在计算pssm时,nr.00可以计算出结果,nr.01也可以计算出结果,并且还不一样。

      1. 我下载的nr解压后近40G,但是解压后只有一个文件,名字是nr,也没有后缀名什么的,程序显示错误为:No alias or index file found,而如果是用nr.00,解压后会有12个nr.00的文件,这样运行程序就没报错。您的nr还有吗,解压后也是只有一个nr文件吗,,跪求错误原因和nr数据包,不胜感激。

        1. 是的,你把nr.01, nr.02……也都下下来,解压缩后放到一起,唯一的一个Nr就是你的索引。nr解压缩后只有一个文件,说明你下载的是fasta文件,还需要建立索引,用makeblastdb建。个人建议你还是下载nr.*系列的

          1. 下载nr.*系列的,那样就会有nr.00,nr.01,nr.02….nr.32,现在更新到30几了,要下载30多个,解压后这30多个文件,每个里面都有nr.pal,这样合并他们就会有30多个nr.pal,,这怎么办?您的意思解压后放一起只有一个Nr索引是什么意思?

  2. 而且如果合并在一起了,程序中执行命令行-db D:\\blast-2.2.30+\\bin\\nr\\nr.00,就改为\\bin\\rn\\nr就可以吗?谢谢。

    1. 可以的,你可以打开Nr.pal看一下,30个文件的Nr.pal都是一样的,作用就是指定nr.**的位置。

  3. 嗯,是一样的,,33个我已经下好了,不需要进行格式化什么的吗?,,直接把他们放在一个文件夹里,就可以用了吗?

发表评论

电子邮件地址不会被公开。 必填项已用*标注