利用NCBI的Blast+进行本地化序列相似性检索

(一) 软件安装

推荐使用bioconda安装

1	conda install -c bioconda blast=2.9.0

在NCBI官方网站选择最新版本,进行下载安装。

(二) 创建本地数据库

1	makeblastdb -in demo_db.fasta -dbtype nucl -out dbname

参数说明:

-in 输入待建库的序列数据,如demo.fasta
-input_type 输入待建库的序列数据格式,可选asn1_bin, asn1_txt, blastdb, fasta,默认为fasta。
-dbtype 序列数据类型,即核酸还是氨基酸。可选为nucl(核酸),prot(氨基酸)
-out 生成的数据名称
其他参数,可以使用makeblastdb -help查看。

(三) 选择检索程序

根据不同的需求,比如所使用的序列是氨基酸还是核酸,要查找的数据是核酸还是氨基酸,选择合适的blast工具。不同的工具使用方法大致相同,在本篇 Demo 中,仅选择更常用的blastn和blastp工具进行介绍。

Query Sequence Type	Databse Sequence Type	Alignment level Type	Program
nucleotide	nucleotide	nucleotide	blastn
peptide	peptide	peptide	blastp
nucleotide	peptide	peptide	blastx
peptide	nucleotide	peptide	tblastn
nucleotide	nucleotide	peptide	tblastx

(四) 序列相似性搜索

# 核酸序列检索
blastn -query demo_query.fa -out demo.blast -db dbname -outfmt 6 -evalue 1e-5 -num_descriptions 10 -num_threads 8

# 氨基酸检索
blastp -query demo_query.fa -out demo.blast -db dbname -outfmt 6 -evalue 1e-5 -num_descriptions 10 -num_threads 8

参数说明:

-query 输入待搜索的序列
-out 输出结果的文件名
-db 上一步骤建立的本地数据库名,也可以是NCBI的在线数据库。
-outfmt 输出结果文件的格式,总共有18种。一般选择6,即Tab分隔的文本文件。
-evalue 设置输出结果的evalue值
-num_descriptions Tab格式输出的条数
-num_threads 程序线程数
-remote: 可以用NCBI的远程数据库, 一般与 -db nr

输出文件格式:

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
11 = BLAST archive (ASN.1)
12 = Seqalign (JSON)
13 = Multiple-file BLAST JSON
14 = Multiple-file BLAST XML2
15 = Single-file BLAST JSON
16 = Single-file BLAST XML2
17 = Sequence Alignment/Map (SAM)
18 = Organism Report

(五) 输出结果

以输出格式 6 的输出结果,是一个Tab分隔的文档文件,一般情况,我们选择提取序列的名称,然后使用脚本从构建数据库的FASTA序列文件中提取序列。

序列提取参考作者之前编写的脚本。

(六) 其他特殊应用:基因注释

对已有序列进行注释时常见的best hit only模式命令行:

blastn -query gene.fa -out gene.blast.txt -task megablast -db nt -num_threads 12 -evalue 1e-10 -best_hit_score_edge 0.05 -best_hit_overhang 0.25 -outfmt "7 std stitle" -perc_identity 50 -max_target_seqs 1

参数详解:

-task megablast : 任务执行模式,可选有’blastn’ ‘blastn-short’ ‘dc-megablast’ ‘megablast’ ‘rmblastn’
-best_hit_score_edge 0.05 : Best Hit 算法的边界值,取值范围为0到0.5,系统推荐0.1
-best_hit_overhang 0.25 : Best Hit 算法的阈值,取值范围为0到0.5,系统推荐0.1
-perc_identity 50 : 相似度大于50
-max_target_seqs 1 : 最多保留多少个联配

(七) 使用blastall进行相似性搜索

除了直接使用blastn、blastn等程序之外,还可以使用blastall来调用程序。

blastall软件程序的参数如下:

-p: 执行的程序名称
-d: 搜索的数据库名称
-i : 要查询的序列文件名(Query File)
-e:(数学)期望值(Expectation value),E值是个统计阈值,缺省值10, 意指比对结果中由于随机偶然性产生的匹配结果不大于10,E值越小结果越可靠。
-o :查询结果输出文件名
-m: 比对结果显示格式选项,缺省值为0 ,即pairwise格式。另外还可以根据不同的需要选择1~6等不同的格式。
-I :在描述行中显示gi号[T/F],缺省值F
-v :单行描述(one-line description)的最大数目,缺省值500
-b :显示的比对结果的最大数目,缺省值250
-F :对于要查询的序列做低复杂度区域(low complexity regions, LCR)的过滤[T/F],缺省值T。对blastn用的是DUST程序,其他比对用的是SEG程序。
所谓“低复杂度区域”是指某些或一些残基过多表现,短周期重复等。对于高等哺乳动物的基因组序列,可以先用RepeatMask程序遮蔽重复元件。在输出结果中,对LCR区的序列核酸用“N”代替,蛋白质序列用“X”代替。
-a:运行BLAST程序所使用的处理器的数目,缺省值1
-S:在数据库中搜索时所使用的核酸链(strand),只对blastn、blastx和tblastx有效;1表示top,2表示bottom,3表示both;缺省值3
-T: 产生HTML格式的输出[T/F],缺省值F
-n: 使用MegaBlast搜索[T/F],缺省值F
-G: 打开一个gap的罚分(0表示使用缺省设置值),默认0
-E: 扩展一个gap的罚分(0表示使用缺省设置值),默认0
-q: 一个核酸碱基的错配(mismatch)的罚分(只对blastn有效),缺省值-3
-r : 一个核酸碱基的正确匹配(match)的奖分(只对blastn有效),缺省值1
-M: 所使用的打分矩阵,缺省值BLOSUM62