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

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

 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:核酸序列对核酸库在蛋白质质级别的比对,两者都在搜索之前翻译成为蛋白质质进行比对

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个我已经下好了,不需要进行格式化什么的吗?,,直接把他们放在一个文件夹里,就可以用了吗?

发表评论

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