Blast
Blast
Blast使用启发式搜索来找出相关的序列，Blast算法不能够像动态规划保证搜索到的序列和要找的序列之间的相关性，BLAST的工作就是尽可能找出数据库中和所要查询的序列相关的信息而已，精确度稍微低一点。
算法
- 1 移除Query序列中之低复杂度以及有串街重复现象的区域
- 低复杂度是指由很少种类的元素（如氨基酸）所组成的一个区域；而串接重复现象是指在一个Query序列中，有两段串连的区域它们组成的方式一模一样。这两种在序列中的区域可能会让BLAST找出一些虽然分数够高，但是其实和Query序列并不相关之序列，所以在我们执行搜索之前，要先把Query序列中的这两种区域滤掉。BLAST的实际作法是，它会把这些区域用符号代替，并且在搜索的时候忽略这些符号。蛋白质序列中，就用X符号标示；而DNA序列中，则用N符号标示。
- 2 将Query序列中每k个字的组合做成一个表
- 以k=3为例（DNA序列中，我们则常以k=11为例），我们"依序"将Query序列中每3个字的组合视为一个字组，并将这些字组列在一张字组表上，直到Query序列中最后一个字也被收入进表上为止.
- 3 列出我们所关心的所有可能的字组
- 4 将这些经筛选后剩下的高分字组组织成快速搜索树的结构
- 5 对每个Query序列中的字组重复步骤1到4，并得到所有相应的高分字组
- 6 扫描数据库中的序列，看看是否有跟剩下的高分字组完全匹配的情形
- 7 将这些完全匹配的地方扩展为高分序对（high-scoring segment pair, HSP）
- 8 将所有分数够高的HSP列出来
- 9 评估这些留下来的HSP它们的分数是否具有意义
- 10 将一个数据库序列中的多个HSP区域结合成一个更长的排比
- 11 显示Query序列和所有之前找到可能相关的数据库序列之有间隙区域排比
- 12 列出上一步骤中期望分数E低于我们所要求的门槛值之数据库序列
简要
Filtering
在进行比对之前，我们需要对query的序列进行过滤（通常序列比对是将未知序列与已知的序列进行比对，我们把未知的序列叫做query序列），主要过滤的是低复杂度区域（"Low-complexity region"），也就是序列的一段区域只包含种类很少的碱基或者氨基酸，即有大量重复，因为这些区域会导致最后的计算得分很高或者混淆Blast的运算。BLAST的实际作法是，它会把这些区域用符号代替，并且在搜索的时候忽略这些符号。蛋白质序列中，就用X符号标示；而DNA序列中，则用N符号标示。低复杂度区域的部分，BLAST是用一个叫做SEG的程序来处理蛋白质序列，而用叫做DUST的程序来处理DNA序列。而且这一步直接就在seeding之前就完成了。
Seeding
首先，需要把比对的序列按照一定长度拆分成多个连续‘seed words’（通常来说对于蛋白质序列是以三个氨基酸为一个seed word, DNA碱基序列以11个长度为一个seed word，这个阈值可以自己定义），word的长度越小，最终的比对准确率就越高，所需要的时间也越长。比如对于一个蛋白质序列PQFEFG，设定seed word长度为3，那么就可以依次形成四个seed word。
Matching
Matching步骤就是与参考序列进行索引的过程，因为我们参考序列一般已经导入到数据库中，并且已经做好了索引了，所以matching这个过程是非常快的，常用的如哈希索引和后缀树都是非常高效的索引方式，哈希是典型的空间换时间的算法，而后缀树面对这种字符也有很好的查询性能，并且能节省储存空间。
这样我们就可以得到每个seed word在参考比对序列对应的位置。
Extending
这个步骤是整个Blast算法里面最耗时的，主要思想就是希望把匹配得到的seed word延伸成更长的片段，同时，计算延伸后的得分，当得分小于指定的阈值，停止延伸。如下图所示，我们知道最优的匹配结果和主对角线方向是平行的，因此我们沿着主对角线方向进行双方向的延伸，直到打分低于一定的阈值停止，我们将最后的比对结果叫做高分片段对（high-scoring segment pair, HSP），这里延伸的方法和Smith-Waterman算法基本一致，就是只计算MATCH的得分，因为Blast初始版本是不允许有GAP的。
显著性分析
Blast使用进行序列比对
download blast software from the website
解压
tar -xvf ncbi-blast-2.11.0+-x64-linux.tar
配置环境变量
cd ncbi-blast-2.11.0+/
vim ~/.bashrc
# export at the end
export PATH=/home/lilingling/hml/blast/ncbi-blast-2.11.0+/bin:$PATH
检测环境配置是否成功
source ~/.bashrc
blastn -version
构建本地blast nr/nt数据库
mkdir data
cd data
# -in 后面接输入文件，即我们要格式的fasta序列
# -dtype 后接序列类型，nucl为核酸，prot为蛋白
# -title 给生成的blast数据库起一个别名
# -parse_seq 帮助我们解析fa文件中，“>”后面的id信息
# -out 后接数据库名称，起一个有意义的名称，后续进行blast比对时，-db参数后跟这个名称
# -logfile: 日志文件，如果没有则输出到标准输出（屏幕）上
makeblastdb -in completeGenme.fasta -dbtype nucl -parse_seqids -out ./database
序列对比
# -query: 需要比对的数据的文件路径以及文件名
# -out: 输出比对结果
# -db: 后面跟着建立好的索引数据库
# -outfmt: 指定输出结果的格式，此处指定为6，即常见的m8格式
# -evalue: 设置输出结果的最小e-value值
# -num_threads: 指定比对时使用的线程数
blastn -query completeGenme.fasta -db ./database -evalue 1e-6 -outfmt 6 -num_threads 6 -out out_file
结果分析
# 一共有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值
head -10 out_file
AC_000192 AC_000192 100.000 31526 0 0 1 31526 1 31526 0.0 58218
AC_000192 FJ647226.1 99.597 31479 112 13 26 31495 1 31473 0.0 57413
AC_000192 FJ647219.1 99.539 31433 130 13 48 31471 1 31427 0.0 57228
AC_000192 FJ647221.1 95.894 31489 1238 34 29 31495 1 31456 0.0 50935
AC_000192 FJ647220.1 95.051 31480 1477 57 40 31489 1 31429 0.0 49426
AC_000192 JX169867 99.529 25285 104 13 22 25297 1 25279 0.0 46041
AC_000192 JX169867 99.809 5772 11 0 25721 31492 25280 31051 0.0 10604
AC_000192 FJ647222.1 98.674 25407 313 22 48 25441 1 25396 0.0 45029
AC_000192 FJ647222.1 96.607 5895 196 4 25591 31483 25391 31283 0.0 9775
AC_000192 FJ647227.1 98.425 25390 376 22 65 25441 1 25379 0.0 44717
