BLAST的一些坑

Overview

因为在使用PSIPRED和SCRATCH都使用到了BLAST，安装和使用的过程中遇到了一些问题，这里主要记录下BLAST安装和使用过程中遇到的一些问题，并没有详细记录完整的BLAST安装和使用说明，毕竟还没有直接使用BLAST，等用到的时候再详细记录。

BLAST，全称Basic Local Alignment Search Tool，即“基于局部比对算法的搜索工具”，能够实现比较两段核酸或者蛋白质序列之间的同源性的功能，具有较快的比对速度和较高的比对精度，适用于多种序列比对的情况，在常规双序列比对分析中应用最广泛。

1.`BLAST`的下载和安装

其实在PSIPRED的安装和使用中，我们已经详细介绍了BLAST的安装过程，可以先看看这篇文章了解完整的安装过程，这里更多的是剖析PSIPRED和SCRATCH对BLAST的使用。
我们依然假定从ftp://ftp.ncbi.nih.gov/blast/executables/release/2.2.26/下载了2.2.26版本的BLAST，因为BLAST提供下载的不是源码包，而是编译之后的可执行程序，因此是平台相关的，要注意选择与你当前系统匹配的版本（操作系统类型&操作系统位数）。

2.关于`BLAST`可执行程序的配置

因为下载的BLAST已经是可执行程序，因此不需要编译，可以直接运行，将下载下来的程序解压之后就会得到blast-2.2.26文件夹，打开blast-2.2.26/bin文件夹，就会看到BLAST所有的可执行程序：

bl2seq      blastclust  copymat     formatdb    impala      megablast   seedtop    blastall    blastpgp    fastacmd    formatrpsdb makemat     rpsblast

其中blastpgp就是最常被其他软件使用的BLAST程序。
将BLAST完整的程序执行路径添加到$PATH路径之后，在任何地方都可以直接使用程序名运行程序，而不必使用完整的路径名，以blastpgp为例，可以直接只用blastpgp + 参数而不是使用完整路径/blastpgp + 参数来运行blastpgp。

2.1 `PSIPRED`中对`BLAST`程序的调用

事实上，上面虽然配置好了BLAST可执行程序路径，使得我们可以在任何地方方便地运行程序，但是当我们使用runpsipred脚本运行PSIPRED时依然会报blastpgp not found的错误，查看runpsipred脚本，发现以下代码：

# Where the NCBI programs have been installed
set ncbidir = /usr/local/bin
...
$ncbidir/blastpgp -b 0 -j 3 -h 0.001 -d $dbname -i     $tmproot.fasta -C $tmproot.chk >& $tmproot.blast
...
$ncbidir/makemat -P $tmproot

runpsipred脚本使用的/usr/local/bin/blastpgp访问blastpgp，因此找不到blastpgp，改为下面的代码，直接调用blastpgp和makemat。

# Where the NCBI programs have been installed
set ncbidir = /usr/local/bin
...
# modified
blastpgp -b 0 -j 3 -h 0.001 -d $dbname -i $tmproot.fasta -C $tmproot.chk >& $tmproot.blast
...
# modified
makemat -P $tmproot

2.2 `SCRATCH`中对`BLAST`程序的调用

如果不是在执行SCRATCH程序时报了blastpgp的错误，你很难感觉到SCRATCH中BLAST的存在，因为SCRATCH的README中没要求你安装和配置BLAST. 这是因为SCRATCH自带了一个32位的linux安装包，不能在Mac或者64位的linux系统上执行。
按照 SCRATCH的安装和使用中的方法下一个正确的版本解压之后覆盖原来cd Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg中的blast-2.2.26文件夹就行了。
如果你下的版本不是blast-2.2.26或者替换之后文件夹不叫blast-2.2.26，比如你的文件夹叫blast,你需要修改SCRATCH-1D_1.1/pkg/PROFILpro_1.1/env文件夹中的PROFILpro.sh文件，将下面这一行：

export PROFILpro_BLAST_INSTALL_DIR="/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/blast-2.2.26"

改为

export PROFILpro_BLAST_INSTALL_DIR="/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/blast

为什么是在这里修改，原因是PROFILpro_1.1程序负责调用与BLAST有关的程序，生成供后面程序使用的中间文件。因此与BLAST有关的环境变量都存放在了SCRATCH-1D_1.1/pkg/PROFILpro_1.1/env文件夹中的PROFILpro.sh文件中。后面的我们要说的BLAST的数据库配置，也是需要从这里入手才能清楚整个流程。

3.关于BLAST本地数据库的配置

BLAST作为一个基于局部比对算法的搜索工具，核心功能自然是从数据库中搜索蛋白或者核苷酸序列。因此BLAST需要正确配置数据库，这样BLAST程序才可以正确访问到数据库。
仍然以BLAST中的blastpgp为例，下面这句：

blastpgp -i example.fasta -d uniref90filt

运行了blastpgp。其中-i 指定了程序的输入文件；-d指定了程序要查询的数据库，如果没有配置数据库，这里-d需要指定数据库存放的绝对路径，很不方便。

3.1 `SCRATCH`中对`BLAST`数据库的使用方式

尽管运行程序时用-d指定数据库的完整路径不是很方便，但是SCRATCH中在使用BLAST时就是采用这种方式，因为SCRATCH在SCRATCH-1D_1.1/pkg/PROFILpro_1.1/data/uniref50下，存放了完整的uniref50数据库，直接使用绝对路径指向数据库，也就减少了用户自己配置的复杂性，这也是为什么SCRATCH程序解压之后有6G+的原因。
PROFILpro_1.1程序使用SCRATCH-1D_1.1/pkg/PROFILpro_1.1/bin/generate_profiles.sh脚本作为程序开始的入口，我们看下这个脚本的部分内容：

...
# Project installation folder
   install_dir="/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/PROFILpro_1.1"
...
# Checking project profile
profile="${install_dir}/env/PROFILpro.sh"; if [ ! -f "$profile" ]; then
echo "$m project profile ($profile) not found"; exit 1; fi

# Checking PROFILpro root script
profilpro="${install_dir}/lib/generate_profiles.sh"; if [ ! -f "$profilpro" ]; then
echo "$m cannot find PROFILpro root script ($profilpro)"; exit 1; fi
...
# Running PROFILpro with database UNIREF50
database="UNIREF50"; $profilpro $profile $database $rounds $fasta $output $threads

这里重点看最后一句代码，$profilpro $profile $database $rounds $fasta $output $threads，实际上是调用了SCRATCH-1D_1.1/pkg/PROFILpro_1.1/lib/generate_profiles.sh脚本，并传入环境变量文件SCRATCH-1D_1.1/pkg/PROFILpro_1.1/env/PROFILpro.sh等参数。
继续查看*SCRATCH-1D_1.1/pkg/PROFILpro_1.1/lib/generate_profiles.sh脚本：

...
# Checking PROFILpro scripts
check_profile="${PROFILpro_LIB_DIR}/check_project_profile.pl"
profilpro="${PROFILpro_LIB_DIR}/generate_profiles.pl"
...
# Checking requested database
db_flag=0; db_dir=""; blastdb="";
if [ "$database" == "UNIREF50" ]; then     db_dir="$PROFILpro_UNIREF50_DIR";
blastdb="$PROFILpro_UNIREF50_BLASTDB"; db_flag=1; fi
if [ $db_flag -ne 1 ]; then echo "$m unknown database : $database"; exit 1; fi

# Running PROFILpro with requested database
$profilpro $db_dir $blastdb $rounds $fasta $output $threads

其中的${PROFILpro_LIB_DIR}，$PROFILpro_UNIREF50_DIR，$PROFILpro_UNIREF50_BLASTDB在SCRATCH-1D_1.1/pkg/PROFILpro_1.1/env/PROFILpro.sh中定义：

...
# Project Installation Folder
export PROFILpro_ROOT_DIR="/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/PROFILpro_1.1"

# Project Static Folders
export PROFILpro_DATA_DIR="${PROFILpro_ROOT_DIR}/data"
export PROFILpro_LIB_DIR="${PROFILpro_ROOT_DIR}/lib"
...
# Project Data Folder
export PROFILpro_UNIREF50="uniref50"
export PROFILpro_UNIREF50_DIR="${PROFILpro_DATA_DIR}/$PROFILpro_UNIREF50"
export PROFILpro_UNIREF50_BLASTDB="$PROFILpro_UNIREF50"

因此，SCRATCH-1D_1.1/pkg/PROFILpro_1.1/lib/generate_profiles.sh脚本中最后的$profilpro $db_dir $blastdb $rounds $fasta $output $threads语句，是调用了SCRATCH-1D_1.1/pkg/PROFILpro_1.1/lib/generate_profiles.pl*脚本，并将/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/PROFILpro_1.1/data/uniref50，uniref50等信息作为参数传进去。

继续查看SCRATCH-1D_1.1/pkg/PROFILpro_1.1/lib/generate_profiles.pl脚本:

...
# Script Inputs
my $usage="database_dir blastdb rounds input_fasta output_file threads";
if($#ARGV!=5){ print "Usage : ./$script_name $usage\n"; exit 1; }
my $db_dir=$ARGV[0]; my $blastdb=$ARGV[1]; my $rounds=$ARGV[2];
my $fasta=$ARGV[3]; my $output=$ARGV[4]; my $threads=$ARGV[5];
...
# Preparing Temporary Workspace
my $lib_dir=$ENV{"PROFILpro_LIB_DIR"}; my $tmp_dir=$ENV{"PROFILpro_TMP_DIR"};
...
# Starting Jobs In Background
my $process_batch="$lib_dir/process_full_batch.pl"; if(! -f "$process_batch"){
print "$m cannot find PROFILpro lib scripts.\n"; exit 1; }     for(my $i=0;$i<$num_batchs;$i++){
my $cmd="$process_batch $db_dir $blastdb $rounds $work_dir $batchs_input[$i] ";

我们发现实际上调用了/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/PROFILpro_1.1/lib/process_full_batch.pl，并继续将/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/PROFILpro_1.1/data/uniref50，uniref50等信息作为参数传进去。

查看/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/PROFILpro_1.1/lib/process_full_batch.pl：

...
# Script Inputs
my $usage="database_dir blastdb rounds work_dir fasta output flag prefix"; if($#ARGV!=7){
print "Usage : ./$script_name $usage\n"; exit 1; } my $db_dir=$ARGV[0]; my $blastdb=$ARGV[1];
my $rounds=$ARGV[2]; my $work_dir=$ARGV[3]; my $fasta=$ARGV[4]; my $output=$ARGV[5];
my $flag=$ARGV[6]; my $prefix=$ARGV[7]; my $m="[$prefix]"; my $f="$work_dir/$flag";

# Third-Party Software
my $blastpgp_exe=$ENV{"PROFILpro_BLAST_BLASTPGP_EXE"};
my $blastpgp_opt=$ENV{"PROFILpro_BLAST_BLASTPGP_OPT"};
...
# Processing Protein Sequences
open(PROFIL,">$work_dir/$output"); my $next_display=5; for(my $p=0;$p<$num_entries;$p++){
my $id=$entries[$p]; my $prov=$provided{"$id"}; my $seq=$sequence{"$id"}; my $len=$seqlen{"$id"};
my $f_in="$prefix.$id.fa"; my $f_aln="$prefix.$id.aln"; my $f_chk="$prefix.$id.chk";
my $f_mat="$prefix.$id.mat"; open(OUT,">$work_dir/$f_in"); print OUT ">$id\n$seq\n"; close(OUT);
my $cmd="$blastpgp_exe -d $db_dir/$blastdb -i $f_in -o $f_aln -C $f_chk -Q $f_mat";

注意最后一句代码，$ENV{"PROFILpro_BLAST_BLASTPGP_EXE"}和$ENV{"PROFILpro_BLAST_BLASTPGP_OPT"}都从SCRATCH-1D_1.1/pkg/PROFILpro_1.1/env/PROFILpro.sh中得到，进行一下变量替换就会发现实际上调用了/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/blast-2.2.26/bin/blastpgp*程序，其中的-d指定了数据库的绝对路径$db_dir/$blastdb，即/Users/wangjiawei/Bioinformatics/Softwares/SCRATCH-1D_1.1/pkg/PROFILpro_1.1/data/uniref50/uniref50。

最容易出错的就是`-d`数据库路径的指定。

注意这里数据库指定的路径和实际上数据库文件的存放的关系，SCRATCH-1D_1.1/pkg/PROFILpro_1.1/data/文件夹下有一个uniref50文件夹，而这个SCRATCH-1D_1.1/pkg/PROFILpro_1.1/data/uniref50文件夹下，存放的才是BLAST可以访问的数据库文件。查看一下SCRATCH-1D_1.1/pkg/PROFILpro_1.1/data/uniref50目录下的内容：

uniref50.00.phr uniref50.00.psq uniref50.01.pin uniref50.02.phr uniref50.02.psq uniref50.03.pin uniref50.04.phr uniref50.04.psq
uniref50.00.pin uniref50.01.phr uniref50.01.psq uniref50.02.pin uniref50.03.phr uniref50.03.psq uniref50.04.pin uniref50.pal

共有5G+，因此是分了4段存储的，在-d中只需要指定数据库的名字即可，不需要指定分段信息，BLAST程序会自动查找这个数据库下所有的分段。
值得注意的是，数据库的每一段都有.phr，.psq和.pin三个文件，这三个文件是下载的标准数据库经过BLAST格式化之后的文件，也就是BLAST`可以访问的数据库格式。

3.2 `PSIPRED`中对`BLAST`数据库的使用方式

PSIPRED因为没有没有自带数据库，所以需要自己配置数据库。
由于数据库太大了，正好SCRATCH中带了uniref50数据库，就可以直接使用这个数据库，为了通用，我们另起一个专门的地方放置这个数据库，详细的步骤参见PSIPRED的安装和使用。