NCBI的BLAST+工具本地安装，本地建库和BLAST比对

接触过生物学的小伙伴对NCBI在线BLAST网页一定不陌生，简单介绍一下这个网页的5种比对工具：blastn、blastp、blastx、tblastn和tblastx，以及如何进行本地建库和blast比对。

blastn 用核苷酸序列检索核苷酸数据库
blastp 用蛋白质序列检索蛋白质数据库
blastx 核苷酸序列通过6种阅读框翻译不同蛋白序列后，检索蛋白数据库
tblastn 蛋白序列比对核酸库，核酸数据库通过6种开放阅读框翻译不同蛋白质
tblastx 核酸序列和核酸数据库都通过6种开放阅读框翻译后比对

平常我们用的最多的就是blastn和blastp，进入网页，选择blast方式，然后贴上自己的quary序列，选择数据库，选择比对的物种，设置参数如E值，wordlength长度等等。但是NCBI网站的BLAST在线工具有个让人特别无语的缺点：国内访问速度巨慢！不仅仅是比对过程慢，一条序列还好，大批量数据比对就不要想了，有的时候网页都打不开一直转圈圈。因此本地化blast工具还是很有必要的。

好在NCBI很贴心地提供了blast+工具，我们安装好blast+工具和下载好数据库以后，就可以不依赖网页和NCBI地服务器，在本地服务器上运行了。

1 安装blast+

最新版blast+工具可以通过ftp方式获得，点击这里

我是在集群账户下安装，集群机器都是linux操作系统的，因此我选择的最新linux版本

别忘了要连md5校验文件一起下载，谁都不知道下载过程中是否会出错，因此大的文件下载完以后都是需要验证文件完整性的！

将上面的ftp地址拼凑一下，网速好的话可以直接用wget下载，但是我这边服务器连NCBI网速实在太慢了，wget只有10Kb/s的速度，甚至还会断开重连。看的我高血压都要犯了，无奈之下挂了个梯子，在自己电脑上下载好这两份文件，通过xftp传到了服务器上。

在服务器上首先校验文件完整性:

md5sum -c ncbi-blast-2.13.0+-x64-linux.tar.gz.md5

显示结果OK后，解压：

tar -zxvf ncbi-blast-2.13.0+-x64-linux.tar.gz

名字太长了，不方便以后找，顺便改个名就叫blast：

mv ncbi-blast-2.13.0+-x64-linux.tar.gz blast

然后是配置环境变量：

1
2
3

vim ~/.bashrc       # 编辑环境变量文件
# 在.bashrc文件最后一行加入如下内容（根据自己路径修改）
export PATH="/public/home/wlxie/blast/bin:$PATH"

保存退出后重新source一下.bashrc文件，blast+工具就安装好了。

2 下载nr/nt数据库

我们比对一般用的是NCBI的非冗余蛋白/核酸数据库，有两种方法下载nr/nt数据库：

1.通过blast+工具自带的更新程序下载
2.通过aspera工具下载

同样是网速的问题，如果用第一种方法下载，我们可以在~/blast/bin目录下找到如下的perl程序

直接运行命令 perl update_blastdb.pl nt

但是10几kb/s的速度真的让人抓狂，所以我推荐第二种方法：用IBM公司开发的快速下载神器——aspera

安装aspera

在我之前写的一篇博客里推荐过sra-tools工具中的prefetch，用来下载SRA数据中存放的高通量测序原始数据。prefetch软件就是默认通过aspera工具进行下载的。

如果之前没有安装过aspera，可以用conda直接安装，命令如下：

conda install -c hcc aspera-cli

这里注意下aspera-cli是aspera的命令行版本，各种不同版本的本质上下载都是调用ascp程序，并且需要openssh公钥认证，不同版本的aspera公钥文件存放的位置不同。因为我们是通过conda安装的aspera，aspera-cli公钥文件的位置在你的conda环境目录下的etc文件夹中，比如我的aspera-cli公钥文件在/public/home/wlxie/miniconda3/envs/biosoft/etc

而且因为是conda安装的，我们不需要修改什么配置文件和依赖关系，还是挺省事的。

用aspera下载数据库nr/nt

nr/nt数据库也可以通过ftp方式获得，点击这里查看ftp网址

为了方便找到下载到本地的数据库，先在家目录新建db/blast文件夹，进入这个文件夹后，在当前目录下运行如下命令：

ascp -QT -i /public/home/wlxie/miniconda3/envs/biosoft/etc/asperaweb_id_dsa.openssh -k1 -l 500m anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/FASTA/nt.gz ./

这里是其中一个nt数据库，nr数据库只要改一个字母就行了，两个数据库都要下载。

稍微解释下各参数的含义：

-Q 用于自适应流量控制，磁盘限制所需；-T 是取消加密，否则有时候数据下载不了。两个参数是搭配一起使用的
-i 输入私钥文件，注意下载的ascp版本不一样文件位置也不一样
-k1 这里是加上了断点传续功能
-l 限制最大下载速度
后面一串是账户@ftp地址:路径。注意@和冒号。NCBI公共账号是anonftp，也就是你下载SRA数据库数据也可以用这个账号；EBI公共账号是era-fasp
最后指定下载文件的路径，我用了当前路径

可以看到下载速度杠杠的，提升了不知道多少倍…下载大数据都可以用ascp命令。

下载好之后同样别忘了校验md5文件，校验后gunzip直接解压到当前文件夹。

3 本地建库

解压完成以后我们可以看到这两个数据库总大小在980G

现在还不能用这两个数据库，需要对这两个超大的数据文件建索引，也就是本地建库。

使用如下命令：

makeblastdb -in nt -dbtype nucl -input_type fasta -out nt

makeblastdb -in nr -dbtype prot -input_type fasta -out nr

-in: 待格式化的序列文件
-dbtype: 数据类型，prot为蛋白序列，nucl为核酸序列
-input_type: 输入数据的类型，默认为fasta格式
-out: 自定义的数据库名称

这一步需要非常长时间，在目录下能看到有文件生成并且没有报错就行了，同样的操作方法可以用自己的基因组数据建库。

这里有两条核苷酸序列可能有问题，序列录到了开头第一行，不过就只有两条序列应该不影响。nt库录入了0.8亿条序列，nr库录入了4.8亿条序列。

4 创建blast全局配置文件

在家目录下创建blast全局配置文件：

$vim .ncbirc        # 家目录下创建一个新文件.ncbirc，输入如下内容

; Start the section for BLAST configuration

[BLAST]

; Specifies the path where BLAST databases are installed

BLASTDB=/public/home/wlxie/db/blast

; Specifies the data sources to use for automatic resolution

; for sequence identifiers

DATA_LOADERS=blastdb

; Specifies the BLAST database to use resolve protein sequences

BLASTDB_PROT_DATA_LOADER=/public/home/wlxie/db/blast/nr

; Specifies the BLAST database to use resolve protein sequences

BLASTDB_NUCL_DATA_LOADER=/public/home/wlxie/db/blast/nt

BATCH_SIZE=10G

; Windowmasker settings

[WINDOW_MASKER]

WINDOW_MASKER_PATH=/public/home/wlxie/db/blast/windowmasker

; end of file

以上设置中定义了blastn和blastp默认的地址，这样我们在比对数据库的时候可以直接输入数据库的名称而不用给出绝对路径，方便一点（这步不是必须的，可选）。

5 运行blast程序

以上准备工作完成后，准备一段query序列试一下，我的query序列名称是gene.fna

运行blastn程序：

blastn -query gene.fna -out gene_blastn_nr.out -db nt -outfmt 6 -evalue 1e-5 -num_threads 10

-query: 用来查询的输入序列
-db: 指定的数据库名称
-out: 自定义输出的结果文件，最好统一格式。我是基因名_比对方法_数据库.out，这样比较直观知道比对了什么，怎么比对的
-outfmt: blast结果的呈现形式，一般用6比较多，也就是m8格式，以制表符为分隔符，有部分信息会缺失。5是XML格式比较适合解析，7在6基础上加了表头。
-evalue: 限定E值
-num_threads: 指定多少个核运行blast程序

还有其他参数比如就不一一介绍了，说明一下，一个序列的blast可以用上面的命令，多个序列的blast同样适用，把多个fasta格式的序列放进去即可。

当然，批量blast的结果需要限定匹配的结果数量，毕竟我们不可能几百上千个序列一一查看，可以指定参数-max_target_seqs 5限制每个序列的最大匹配数量，这个数值推荐是在5以上，5以下会有警告信息。

blast结果m8格式如下：

一共12列，分别能获得如下信息：

1、Query id：查询序列ID标识

2、Subject id：比对上的目标序列ID标识

3、% identity：序列比对的一致性百分比

4、alignment length：符合比对的比对区域的长度

5、mismatches：比对区域的错配数

6、gap openings：比对区域的gap数目

7、q. start：比对区域在查询序列(Query id)上的起始位点

8、q. end：比对区域在查询序列(Query id)上的终止位点

9、s. start：比对区域在目标序列(Subject id)上的起始位点

10、s. end：比对区域在目标序列(Subject id)上的终止位点

11、e-value：比对结果的期望值，解释是大概多少次随即比对才能出现一次这个score，Evalue越小，表明这种情况从概率上越不可能发生，那么发生了即说明这更有可能是真实的相似序列

12、bit score：比对结果的bit score值，越高越好