快速序列比对之BLAST

BLAST的大名对于做生物的同学可以说是如雷贯耳，哪怕非生信的同学也或多或少接触过这个东西。我们通常会使用ncbi的blast在线工具进行比对，但具有一个缺点就是很慢，因此我们常常会搭建一个本地blast系统进行比对分析。

NCBI提供了一套用于运行BLAST的命令行工具，称为BLAST+。这允许用户在自己的服务器上执行BLAST搜索，而不受大小、容量和数据库的限制。BLAST+可以使用命令行运行，对于linux用户可以说是相当友好了。

blast+的软件安装

wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.10.1+-x64-linux.tar.gz
tar -zxvf ncbi-blast-2.10.1+-x64-linux.tar.gz

# 加入环境变量
export PATH=$PATH:$PWD/ncbi-blast-2.10.1+/bin/

nr/nt数据库下载和使用

nr指的是protein sequence，而nt指的是nucleotide。所以nr/nt数据库指的是两个数据库，前者是蛋白质数据库，后者是核酸序列数据库。
计算资源足够而带宽不足时，建议下载fasta序列，然后使用blast的本地命令建立index。命令如下

wget -c ftp://ftp.ncbi.nih.gov/blast/db/FASTA/nr.gz # nr数据库 
wget -c ftp://ftp.ncbi.nih.gov/blast/db/FASTA/nt.gz # nt数据库

# 解压缩
gunzip nr.gz
gunzip nt.gz

# 构建本地blast nr/nt数据库
mkdir nr_db; mkdir nt_db
makeblastdb -in nr -dbtype prot -title my_nr -parse_seqids -out ./nr_db/nr -logfile make_nr.log
makeblastdb -in nt -dbtype nucl -title my_nt -parse_seqids -out ./nt_db/nt -logfile make_nt.log

makeblastdb常用参数的解释：

-in: 后面接输入文件，即我们要格式的fasta序列
-dbtype: 后接序列类型，nucl为核酸，prot为蛋白
-title: 给生成的blast数据库起一个别名
-parse_seqids: 帮助我们解析fa文件中，“>”后面的id信息
-out: 后接数据库名称，起一个有意义的名称，后续进行blast比对时，-db参数后跟这个名称
-logfile: 日志文件，如果没有则输出到标准输出（屏幕）上

带宽足够而计算资源不足，建议下载已经建立好索引的nr/nt库，命令如下：

mkdir nr_db; cd nr_db
wget -c ftp://ftp.ncbi.nih.gov/blast/db/nr*
for i in *gz
do
tar -zxvf $i
done

mkdir ../nt_db; cd ../nt_db
wget -c ftp://ftp.ncbi.nih.gov/blast/db/nt*
for i in *gz
do
tar -zxvf $i
done

除了可以使用wget下载已经建立好索引的nr/nt库还可以使用blast+自带的update_blastdb.pl脚本下载。

update_blastdb.pl --decompress nr [*]
update_blastdb.pl --decompress nt [*]
# 下载压缩后的nr/nt索引数据库，并进行解压

更推荐wget下载方法，可以断点继续下载。

序列比对

以核酸序列为例进行序列比对

blastn -query test.fa -out test.result -db nr -outfmt 6 -evalue 1e-5 -num_threads 4

常用参数介绍：

-query: 需要比对的数据的文件路径以及文件名
-out: 输出比对结果
-db: 后面跟着建立好的索引数据库
-outfmt: 指定输出结果的格式，此处指定为6，即常见的m8格式
-evalue: 设置输出结果的最小e-value值
-num_threads: 指定比对时使用的线程数

此外，如果我们的目的是看测序数据的污染物来源，推荐加上下面两个参数：

-subject_besthit：数据库中中的所有sequence，只输出blast最优的那个结果
-max_target_seqs：设置每条 query reads所能比对上的最多的个数，这里推荐设置为1。

输出结果的格式解析

输出的比对结果，一共有12列，分别代表：

Query id：查询序列ID标识
Subject id：比对上的目标序列ID标识
% identity：序列比对的一致性百分比
alignment length：符合比对的比对区域的长度
mismatches：比对区域的错配数
gap openings：比对区域的gap数目
q. start：比对区域在查询序列(Query id)上的起始位点
q. end：比对区域在查询序列(Query id)上的终止位点
s. start：比对区域在目标序列(Subject id)上的起始位点
s. end：比对区域在目标序列(Subject id)上的终止位点
e-value：比对结果的期望值
bit score：比对结果的bit score值

我们一般看第1，3，11，12列。其中第一列告诉我们比对到了哪个物种，剩余三列告诉我们比对的可信程度。

导师爱说，垃圾进垃圾出……但哪怕是垃圾，你也要知道它是什么垃圾。湿垃圾or干垃圾，或者像我一样是个学术垃圾~