BLAST的主要理念
blast是区段比对，对于给定的两个序列，blast会把具有相似性的片段（hit）找出来，显示的是hit的信息。
- 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.
一、 linux 安装BLAST
linux下的blast是一系列工具集合，包括：blastn, blastp, blastx, tblastn, tblastx等。
采用conda安装，可能出现网络连接超时等问题，重复尝试：
# 安装blast
$>>conda install -c bioconda blast
# blast安装perl模块的方法
$>>conda install perl-digest-md5
二、 BLAST使用方法
第（1）步：makeblastdb命令建立检索所需“数据库”
blast数据库有两种类型：核酸数据库（关键字：nucl）和氨基酸数据库（关键字：prot）。
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]
[-blastdb_version version] [-max_file_sz number_of_bytes]
[-logfile File_Name] [-taxid TaxID] [-taxid_map TaxIDMapFile] [-version]
DESCRIPTION
Application to create BLAST databases, version 2.9.0+
REQUIRED ARGUMENTS
-dbtype <String, `nucl', `prot'>
Molecule type of target db
……
示例：
以拟南芥氨基酸数据为例（拟南芥基因组数据方法类似），下载地址：ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
下载后解压得到氨基酸序列文件（就是一个普通的fasta格式文件）：Arabidopsis_thaliana.TAIR10.pep.all.fa
1）构建拟南芥氨基酸的blast数据库命令：
$>>makeblastdb -in Arabidopsis_thaliana.TAIR10.pep.all.fa -input_type fasta -dbtype prot -out TAIR10 -parse_seqids
- -in：输入文件名/数据库名。本例中就是上述拟南芥氨基酸fasta格式数据。
-
-input_type：输入文件中序列类型（
asn1_bin, asn1_txt, blastdb, fasta，默认是fasta）。 -
-dbtype：目标数据库的类型，包括两种：
nucl（核酸）和prot（氨基酸）。本例中是prot。 - -out：待生成的数据库的名称（默认使用输入文件名）。本例中命名为TAIR10。
-
-parse_seqids：如果输入文件是fasta格式的，则该参数会自动解析序列id，如果输入是其他格式，则序列id自动解析。（该参数非必须）
（注意：-in和-out命令后面的文件可以带上路径，比如：/data/home/chengaoxiang/blast_data_cgx/Arabidopsis_thaliana.TAIR10.pep.all.fa。没有路径则表示都在同一个目录下，包括终端。）
（说明：数据库生成后，原始文件Arabidopsis_thaliana.TAIR10.pep.all.fa可删除。）
执行结果：
Building a new DB, current time: 01/13/2023 09:34:33
New DB name: /data/home/chengaoxiang/blast_data_cgx/TAIR10
New DB title: Arabidopsis_thaliana.TAIR10.pep.all.fa
Sequence type: Protein
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 48321 sequences in 1.11516 seconds.
数据库建成后，生成了下图中的几个文件：
2）查看blast数据库信息：
- 查看信息
blastdbcmd -db TAIR10 -dbtype prot -info - 查看所有数据
blastdbcmd -db TAIR10 -dbtype prot -entry all | head - 根据具体关键字查看，如GI号（NCBI特有，自建数据库没有）
blastdbcmd -db TAIR10 -dbtype prot -entry 3 | head
还有其他有意思的参数，可以查看帮助信息：blastdbcmd -help
第（2）步：选择blast中不同的工具进行计算
根据不同的需求，比如说你的待比对序列是氨基酸还是核苷酸，搜索数据库是核甘酸还是氨基酸，选择合适的blast工具。不同需求的对应关系可以见下图(来自biostars handbook)。
（不同工具的应用范围虽然不同，但是基本参数都是一致的，学会blastn，基本上其他诸如blastp，blastx也都会了。）
- blastn：将给定的核酸序列与核酸数据库中的序列进行比对。
- blastp：将给定氨基酸序列与氨基酸数据库中的序列进行比对。作用：可以寻找较远源的序列；
- blastx：将给定的核酸序列按照六种阅读框架将其翻译成氨基酸序列，并与氨基酸数据库中的序列进行比对。作用：对分析新序列和EST（Expressed Sequence Tag的缩写，即表达序列标签，指从一个随机选择的cDNA 克隆，进行5’端和3’端单一次测序挑选出来获得的短的cDNA 部分序列。）很有用。
- tblastn：将给定的氨基酸序列与核酸数据库中的序列（双链）按不同的阅读框进行比对。作用：对于寻找数据库中序列没有标注的新编码区很有用；
- tblastx：只在特殊情况下使用，它将DNA被检索的序列和核酸序列数据库中的序列按不同的阅读框全部翻译成蛋白质序列，然后进行蛋白质序列比对。
（延伸知识：一个DNA顺序可能有3种阅读框，但只有一种具有编码作用，称为开放阅读框（open reading frame or ORF）。有的阅读框因终止密码出现频繁，故不能生成蛋白，这种阅读框称为封闭阅读框（block reading frame）。若一个顺序所有的三个阅读框都是封闭的，则它无编码蛋白的功能。一个翻译成蛋白的顺序有一个阅读框，开始于AUG起始密码子，通过一系列有义密码子，直到终止密码子结束。通常3个阅读框中总有封闭阅读框的存在。
六种阅读框架：一段DNA或RNA序列有多种不同读取方式，因此可能同时存在许多不同的开放阅读框架。开放阅读框包含一段可以编码蛋白的碱基序列，不能被终止子打断。 当一个新基因被识别，其DNA序列被解读，人们仍旧无法搞清相应的蛋白序列是什么。这是因为在没有其它信息的前提下，DNA序列可以按六种框架阅读和翻译（核酸有两条链，每条链三种，对应三种不同的起始密码子）。）
示例：
以人血红蛋白alpha亚基（HBA_HUMAN.FASTA，内含1个氨基酸序列）作为查询序列（query sequence）（所谓查询序列就是我们要用于在数据库中比对的序列，目的在于从数据库中找到与查询序列相似的序列），如下所示，：
$>>cat HBA_HUMAN.FASTA
>sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
AVHASLDKFLASVSTVLTSKYR
然后在前面建的TAIR10数据库中进行比对，根据上表，应采用blastp方法进行对比，比对时有很多参数可以选择，下面挑重点的几项说明：
1）最简洁比对命令
$>>blastp -query /data/home/chengaoxiang/blast_test/HBA_HUMAN.FASTA -db /data/home/chengaoxiang/blast_data_cgx/TAIR10 | less
- −query：查询序列（query sequence），这里是xx/HBA_HUMAN.FASTA；
- −db：检索的目标数据库，这里是xx/TAIR10；
比对结果：
列出了所有被筛出来的氨基酸序列（6个），同时后面分别列出了每个氨基酸的详细信息。
BLASTP 2.9.0+
……
Score E
Sequences producing significant alignments: (Bits) Value
AT4G37710.1 pep chromosome:TAIR10:4:17718118:17718742:1 gene:AT4G... 31.2 0.17
AT4G37710.2 pep chromosome:TAIR10:4:17718128:17718771:1 gene:AT4G... 31.2 0.17
AT2G03880.1 pep chromosome:TAIR10:2:1181512:1183703:1 gene:AT2G03... 29.6 0.99
AT4G16420.2 pep chromosome:TAIR10:4:9262621:9266157:-1 gene:AT4G1... 26.9 9.7
AT4G16420.3 pep chromosome:TAIR10:4:9262621:9266086:-1 gene:AT4G1... 26.9 9.7
AT4G16420.1 pep chromosome:TAIR10:4:9262621:9266157:-1 gene:AT4G1... 26.9 9.7
>AT4G37710.1 pep chromosome:TAIR10:4:17718118:17718742:1 gene:AT4G37710 transcript:AT4G37710.1
gene_biotype:protein_coding transcript_biotype:protein_coding
gene_symbol:VQ29 description:VQ motif-containing
protein 29 [Source:UniProtKB/Swiss-Prot;Acc:Q9SZG3]
……
2）将比对结果保存在文件中命令（关键参数-out）：
$>>blastp -query /data/home/chengaoxiang/blast_test/HBA_HUMAN.FASTA -db /data/home/chengaoxiang/blast_data_cgx/TAIR10 -out HBA_SW.TXT
-
-out：用于保存结果信息的文件名（可以包括路径），这里是HBA_SW.TXT。
示例（略）。
3）选择输出结果格式命令（关键参数-outfmt）：
$>>blastp -query /data/home/chengaoxiang/blast_test/HBA_HUMAN.FASTA -db /data/home/chengaoxiang/blast_data_cgx/TAIR10 -outfmt 7 | less
-
-outfmt：指定输出信息格式（可选格式
0、1、2……、18，共19种，可以通过blastp -help找到-outfmt说明部分了解每种格式的意义；默认是0；本例选择格式7。其中6、7、10、17四种格式可以自定义输出字段（可选字段见-outfmt说明部分），默认输出字段为：“qaccver saccver pident length mismatch gapopen qstart qend sstart send evalue bitscore”，指定输出字段的命令示例为：blastp -query xx/HBA_HUMAN.FASTA -db xx/TAIR10 -outfmt ' 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send qseq sseq evalue bitscore' | less）。
执行结果：
# BLASTP 2.9.0+
# Query: sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
# Database: /data/home/chengaoxiang/blast_data_cgx/TAIR10
# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 6 hits found
sp|P69905|HBA_HUMAN AT4G37710.1 48.276 29 15 0 78 106 45 73 0.17 31.2
sp|P69905|HBA_HUMAN AT4G37710.2 48.276 29 15 0 78 106 45 73 0.17 31.2
sp|P69905|HBA_HUMAN AT2G03880.1 31.579 38 26 0 76 113 41 78 0.99 29.6
sp|P69905|HBA_HUMAN AT4G16420.2 27.419 62 39 2 15 76 125 180 9.7 26.9
sp|P69905|HBA_HUMAN AT4G16420.3 27.419 62 39 2 15 76 125 180 9.7 26.9
sp|P69905|HBA_HUMAN AT4G16420.1 27.419 62 39 2 15 76 125 180 9.7 26.9
# BLAST processed 1 queries
上述输出格式每列的意义：
query acc.ver（查询序列）、subject acc.ver（结果序列）、% identity（相似度）、alignment length（比对长度）、mismatches（错配数）、gap opens（起始空位数）、q. start（查询序列起始位点）、q. end（查询序列终止位点）、s. start（结果序列起始位点）、s. end（结果序列终止位点）、evalue（期望值）、bit score（得分）。
4）★设置evalue值（非常重要，前面几个示例没有设置该参数，实际中得设置）
将3）中的命令增加“-evalue 0.2”，表示只输出evalue小于0.2的匹配结果。从输出结果可以看出，只输出了evalue=0.17的两条序列。
$>>blastp -query /data/home/chengaoxiang/blast_test/HBA_HUMAN.FASTA -db /data/home/chengaoxiang/blast_data_cgx/TAIR10 -evalue 0.2 -outfmt 7 | less
输出结果：
# BLASTP 2.9.0+
# Query: sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
# Database: /data/home/chengaoxiang/blast_data_cgx/TAIR10
# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 2 hits found
sp|P69905|HBA_HUMAN AT4G37710.1 48.276 29 15 0 78 106 45 73 0.17 31.2
sp|P69905|HBA_HUMAN AT4G37710.2 48.276 29 15 0 78 106 45 73 0.17 31.2
# BLAST processed 1 queries
blast对比对结果的评价主要有两个参数：evalue（E值）和bit score（S值）
- evalue（E值）：是一个对bit score值进行校正之后的值，其意义是给出比对结果（bit socre）的可靠性评价，其值越小越好，越小表示比对结果可靠性越大。（默认evalue是10）
- bit score（比特打分）（S值）：评价比对结果（同源性或相似相），其值越大越好。（bit score通过对score（原始打分）进行归一化处理，使得使用不同打分矩阵得到的结果依然可以进行对比。）
- score（原始打分）：blast的时候，根据打分矩阵（如BLOSUM62、PAM矩阵、PSSM 位置特异性矩阵）打的分，称为原始打分，其值越大越好。但不同的打分矩阵标准不同，因此无法直接比较结果，于是有了bit score。
E值的意义：
在一次blast分析过程中，假设query序列与subject序列（数据库序列）命中了N个hit（命中片段），这N个hit一定可信吗？不一定，因为即使随机产生一组序列，也有可能和query序列存在hit。基于这点，生成一个大小与原数据库相同的随机数据库，然后计算query与该随机数据库之间的hit（假设得到了M个hit）。最后判断M中有多少个hit与N中的hit相同或更好，假设有H个。H应该越小越好，H越小说明N中的hit都是通过blast方法得到的，而不是随机得到的。E值就是这种随机性的度量（公式详见后面），与这里的H等价，因此E也是越小越好，它越小，说明被随机匹配到的可能性越小。这就是E的意义。
E经验值：
E<1e-5：匹配结果基本可以确保是同源关系，这是一个非常高的置信度；
1e-5<E<0.01：匹配结果可以视为同源性关系；
0.01<E<10：匹配结果基本没有同源性，但可能存在暂时的远程同源关系，需要额外的证据来确认；
10<EE：没有同源性，或极远的同源性，但已经低于当前方法所能达到的检测极限了。
E值和S值的计算方式：
S=？？？？？？？？？？？？？（待增加）：比特分数是序列数据库的所需大小，其中当前匹配可能只是偶然发现。比特分数是log 2缩放和标准化的原始分数。每增加一倍，所需的数据库大小（2的"位分数"次方）就会增加一倍。
E=K*m*n*(e-λ*S)：（K和λ与数据库和算法有关，是常量；m是目标序列长度，n是数据库大小，S就是前面提到的S值。）【根据公式，对于某个给定的query序列，其E值随数据库增大而增大（因为参数K和λ），但该query序列和原数据库序列中某个序列的进化关系其实并未改变。这意味着，随着数据库的增大，相同query和数据库序列对之间的匹配可信度会降低，即随着数据库的增大，原来被检测出的有效结果可能会被淹没。因此，需要一种替代E值计算的方法。】
E值的局限性：
- 当目标序列过短时，由于无法得到较大的S值，导致E值偏大。
- 当两序列同源性虽然高，但有较大的gap（空隙）时，S值会下降。这个时候gap scores就非常有用。
- 有些序列的非功能区有较低的随机性时，可能会造成两序列较高的同源性。
5）设置字长“-word_size”
$>>blastp -query xx/HBA_HUMAN.FASTA -db xx/blast_data_cgx/TAIR10 -evalue 0.2 -outfmt 7 -word_size 2 | less
- -word_size：该值越大，匹配成功率越低，匹配到的结果同源性越近；
6）设置打分矩阵（-matrix）
$>>blastp -query xx/HBA_HUMAN.FASTA -db xx/blast_data_cgx/TAIR10 -evalue 0.2 -outfmt 7 -matrix PAM250 | less
- -matrix：一般默认计分矩阵为BLOSUM62（可设置PAM250）
7）psiblast（循环搜索命令） : 敏感度更高的蛋白序列与蛋白序列之间的比对
Position-Specific Iterated (PSI)-BLAST，是一种更加高灵敏的Blastp程序，对于发现远亲物种的相似蛋白或某个蛋白家族的新成员非常有效。当你使用标准的Blastp比对失败时，或比对的结果仅仅是一些假基因或推测的基因序列时（”hypothetical protein” or “similar to…”），你可以选择PSI-BLAST重新试试。
$>>psiblast -query xx/HBA_HUMAN.FASTA -db xx/blast_data_cgx/TAIR10 -evalue 0.2 -outfmt 7 -num_iterations 2 | less
- -num_iterations：迭代次数。
8）双序列比对，指定起始位点
以癌胚抗原搜索结构域为例（网上例子）：
$>>blastp -query CEA21_HUMAN.FASTA -query_loc 147-231 -subject CEAM5_HUMAN.FASTA -subject_loc 145-675 -outfmt 7 | less
- −query：主序列 ；
- −query_loc xx−xx：主序列起止位点（如“−query_loc 23-87”）；
- −subject：副序列；
- −subject_loc xx-xx：副序列起止位点（如“−subject_loc 321-456”）；
9）PHI-BLAST : 模式发现迭代BLAST
PHI-BLAST, 模式发现迭代BLAST, 用蛋白查询来搜索蛋白数据库的一个程序。仅仅找出那些查询序列中含有的特殊模式的对齐。
PHI的语法详细介绍看这里：http://www.ncbi.nlm.nih.gov/blast/html/PHIsyntax.html （不过，不是太明白）
10）NCBI-blast（在NCBI上blast）
在百度查询；
https://wenku.baidu.com/view/7a716acca1c7aa00b52acb14?aggId=138471d707a1b0717fd5360cba1aa81144318fc3&fr=catalogMain&wkts=1673851517291&bdQuery=blast%E6%AF%94%E5%AF%B9%E5%AD%97%E9%95%BF
https://blast.ncbi.nlm.nih.gov/Blast.cgi
参考资料：
http://www.360doc.com/content/21/0714/12/76149697_986499307.shtml
https://www.jianshu.com/p/e5527735f163?utm_campaign=maleskine&ivk_sa=1024320u
https://www.metagenomics.wiki/tools/blast/evalue
https://blog.csdn.net/weixin_43364556/article/details/107099313
https://blog.csdn.net/withbeginner/article/details/124125746
http://blog.shenwei.me/local-blast-installation/
https://www.ncbi.nlm.nih.gov/books/NBK62051/
