你知道什么叫做BLAST吗

我第一次听说BLAST是在大二的结束的那个暑假。我不知道BLAST是干嘛的，但是当自己在NCBI上输入一些序列之后，点击确定，过了几分钟跳出一页常常的结果时，感觉自己突然就很厉害了。

这就是传说中的好用的工具，你不需要知道他是什么原理，你只要随便敲敲，随便摸索一下，就能学会使用了，而且结果还符合自己的需要。

关于BLAST的教程，其实在连载系列Biostar: 课程11、12 中就介绍过了，内容都其实差不多，其实还比我更全面，但是没办法，我的排版好看呀。

• Search may take place in nucleotide and/or protein space or translated spaces where nucleotides are translated into proteins.

• Searches may implement search “strategies”: optimizations to a certain task. Different search strategies will return different alignments.

• Searches use alignments that rely on scoring matrices

• Searches may be customized with many additional parameters. BLAST has many subtle functions that most users never need.

1. 用makeblastdb为BLAST提供数据库

2. 选择blast工具，blastn,blastp

3. 运行工具，有必要的还可以对输出结果进行修饰

第一步：建立检索所需数据库

BLAST数据库分为两类，核酸数据库和氨基酸数据库，可以用makeblastbd创建。可以用help参数简单看下说明。

$ makeblastdb -help
USAGE
  makeblastdb [-h] [-help] [-in input_file] [-input_type type]
    -dbtype molecule_type [-title database_title] [-parse_seqids]
    [-hash_index] [-mask_data mask_data_files] [-mask_id mask_algo_ids]
    [-mask_desc mask_algo_descriptions] [-gi_mask]
    [-gi_mask_name gi_based_mask_names] [-out database_name]
    [-max_file_sz number_of_bytes] [-logfile File_Name] [-taxid TaxID]
    [-taxid_map TaxIDMapFile] [-version]

-dbtype <String, `nucl', `prot'>

具体以拟南芥基因组作为案例，介绍使用方法：
注： 拟南芥的基因组可以在TAIR上下在，也可在ensemblgenomes下载。后者还可以下载其他植物的基因组

# 下载基因组
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gzip -d Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# 构建核酸BLAST数据库
makeblastdb -in Arabidopsis_thaliana.TAIR10.dna.toplevel.fa -dbtype nucl -out TAIR10 -parse_seqids

# 下载拟南芥protein数据
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
# 构建蛋白BLAST数据库
gzip -dArabidopsis_thaliana.TAIR10.pep.all.fa.gz
makeblastdb -in  Arabidopsis_thaliana.TAIR10.pep.all.fa -dbtype prot -out TAIR10 -parse_seqids

如果你从NCBI或者其他渠道下载了格式化过的数据库，那么可以用blastdbcmd去检索blast数据库，参数很多，常用就如下几个

• db string :  string表示数据库所在路径

• dbtype string,: string在（guess, nucl, prot)中选择一个

• 检索相关参数

• -entry all 或 555, AC147927 或 gnl|dbname|tag’

• -entry_batch 提供一个包含多个检索关键字的文件

• -info 数据库基本信息

• 输出格式 -outfmt %f %s %a %g …默认是%f

• out 输出文件

• show_blastdb_search_path: blast检索数据库路径

使用案例

# 查看信息
blastdbcmd -db TAIR10 -dbtype nucl -info
# 所有数据
blastdbcmd -db TAIR10 -dbtype nucl -entry all | head
# 具体关键字，如GI号
blastdbcmd -db TAIR10 -dbtype nucl -entry 3 | head

还有其他有意思的参数，可以看帮助文件了解

第二步：选择blast工具

根据不同的需求，比如说你用的序列是氨基酸还是核苷酸，你要查找的数据是核甘酸还是干计算，选择合适的blast工具。不同需求的对应关系可以见下图（来自biostars handbook）

不同工具的应用范围虽然不同，但是基本参数都是一致的，学会blastn，基本上其他诸如blastp,blastx也都会了。

blastn的使用参数很多blastn [-h]，但是比较常用是如下几个

• -db : 数据库在本地的位置，或者是NCBI上数据库的类型,如 -db nr
  

• -query: 检索文件

• -query_loc : 指定检索的位置

• -strand： 搜索正义链还是反义链，还是都要

• out : 输出文件

• -remote： 可以用NCBI的远程数据库， 一般与 -db nr

• -evalue 科学计数法，比如说1e3，定义期望值阈值。E值表明在随机的情况下，其它序列与目标序列相似度要大于这条显示的序列的可能性。 与S值有关，S值表示两序列的同源性，分值越高表明它们之间相似的程度越大

  E值总结： 1.E值适合于有一定长度，而且复杂度不能太低的序列。2. 当E值小于10-5时，表明两序列有较高的同源性，而不是因为计算错误。3. 当E值小于10-6时，表时两序列的同源性非常高，几乎没有必要再做确认。

• 与联配计分相关参数： -gapopen，打开gap的代价；-gapextend， gap延伸的代价；-penalty：核算错配的惩罚； -reward， 核酸正确匹配的奖励；

• 结果过滤：-perc_identity, 根据相似度

注 BLAST还提供一个task参数，感觉很有用的样子，好像会针对任务进行优化速度。

第三步：运行blast，调整输出格式。

我随机找了一段序列进行检索

echo '>test' > query.fa 
echo TGAAAGCAAGAAGAGCGTTTGGTGGTTTCTTAACAAATCATTGCAACTCCACAAGGCGCCTGTAATAGACAGCTTGTGCATGGAACTTGGTCCACAGTGCCCTACCACTGATGATGTTGATATCGGAAAGTGGGTTGCAAAAGCTGTTGATTGTTTGGTGATGACGCTAACAATCAAGCTCCTCTGGT >> query.fa

用的是blastn 检索核酸数据库。最简单的用法就是提供数据库所在位置和需要检索的序列文件

blastn -db BLAST/TAIR10 -query query.fa
# 还可以指定检索序列的位置
blastn -db BLAST/TAIR10 -query query.fa  -query_loc 20-100
# 或者使用远程数据库
blastn -db nr -remote -query query.fa 
blastn -db nt -remote -query query.fa

以上是默认输出，blast的-outfmt选项提供个性化的选择。一共有18个选择，默认是0。
    0 = Pairwise,
    1 = Query-anchored showing identities,
    2 = Query-anchored no identities,
    3 = Flat query-anchored showing identities,
    4 = Flat query-anchored no identities,
    5 = BLAST XML,
    6 = Tabular,
    7 = Tabular with comment lines,
    8 = Seqalign (Text ASN.1),
    9 = Seqalign (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，7，10，17可以自定输出格式。默认是

qaccver saccver pident length mismatch gapopen qstart qend sstart send evalue bitscore

AC： accession

以格式7为实例进行输出，并且对在线数据库进行检索

blastn -task blastn -remote -db nr -query query.fa  -outfmt 7 -out query.txt
head -n 15 query.txt

PS: 这是一篇欠了我师兄好久的教程，现在终于补上了

更多更全的内容点击原文链接，看看biostars handbook系列。