blast简介及格式解读及练习题
01
blast产生背景
双序列比对可以采用是基于动态规划算法的Needleman-Wunsch（NW）和Smith-Waterman algorithm（SW）算法，虽然精度高，但计算消耗大。当与数据库比对的时候，该算法就显得不切实际。因此TASTA，blast采用启发式算法使得通过大幅度丢失灵敏度来减少运行时间。与FASTA软件相比，blast通过把搜索限制在狭隘的矩阵对角线条带上，来改进FASTA进行数据库搜索的速度。
02
blast的大致原理
blast 程序首先查询query序列的所有子序列，储存在哈希表中。收索数据库中所有与子序列精确匹配的序列，作为种子，向两个方向继续延伸每个精确匹配。期间不允许有空位和错配的情况。然后在限制性区域内；连接延伸的匹配序列，期间允许空位和错配，比对分值要大于设定的阈值。阈值越大，需要匹配的计算越小，软件计算速度越快。仅仅对对延伸匹配进行连接的区域(限制性区域)，而不是整个矩阵，是blast 相对于其他算法速度提高的关键，是以牺牲对角线带以外的任何匹配信息为代价，因此并不能确保query序列与数据库比对结果是最优的比对结果。
03
blast的格式解读
因为blast可以进行本地化，网上教程很多，这里不再详细介绍。根据不同的参数可以输出多种比对格式，例如HTML, plain text, XML等。因为输出的格式多样，我们以常用的M8格式进行简单的介绍。
这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：比对结果的期望值，将比对序列随机打乱重新组合，和数据库进行比对，如果功能越保守，则该值越低；该E值越高说明比对的高得分值是由GC区域，重复序列导致的。对于判断同源性是非常有意义的几个参数。
bit score：比对结果的bit score值
04
习题
4.1）blast的产生背景是什么？
4.2）blast采用的是什么算法
4.3）和动态规划算法相比，blast 的算法有什么优点？
4.4) blast和FASTA软件算法之间的区别是什么？
4.5) blast核酸比对中默认的word是多少？
4.6）word的大小是如何决定blast的运行速度？
4.7）blast比FSAT软件搜索速度快的原因是什么？
4.8) blast是对什么建立索引的？
4.9）blast建立索引的目的是什么？
4.10）blast比对输出的结果有哪些格式
4.11）在M8格式中共有多少列，每一列代表的是什么意思？
4.12）E值的含义是什么？
4.13）统计test.blast有多少条query序列
4.14）统计比对得分最低的query序列
4.15）将比对长度大于200(QueryLen)且比对相似率(Identities%)大于90的信息输出来
4.16）找出比对最长的基因的ID (即QueryLen值最大)
4.17）按照BitScore分值(第12列)的大小对整个文件进行排序(从大到小)
4.18）找出比对长度大于100且E值小于3e-17的比对
4.19）找出序列一致性大于60%，得分高于70的比对信息
4.20）输出gap最多的比对信息
05
参考资料
https://en.wikipedia.org/wiki/BLAST#Output
06
练习题
首先需要自行产生一个 blast.test 文件，公众号教程推文不方便展示。
1) 查看 blast.test 文件的前6行
head -6 blast.test2) 统计blast.test文件中共有多少行，多少列
wc -l blast.test
head -1 blast.test |tr '\t' '\n' |wc -l3）新建一个文件miniblast.text
touch miniblast.text4）将blast.test文件前500行内容写入到前一步新建的miniblast.text 文件
head -500 blast.test > miniblast.text思考一下，有没有可能不需要新建miniblast.text 文件
5) 删除原文件 blast.test
rm blast.test6) 查看QueryId为evm.model.scaffold_1.187在miniblast.text文件中的第几行
grep -n evm.model.scaffold_1.187 miniblast.text |less -S
# less -SN miniblast.text7） 查看ribosomal在miniblast.text文件中出现的次数
grep -c ribosomal miniblast.text8）查看miniblast.text文件中不含有protein的行
grep -v protein miniblast.text |less -S9）单独查看miniblast.text文件中第12列BitScore中所有的数值
cut -f 12 miniblast.text |less -S10) 将miniblast.text文件中第12列BitScore中所有的数值去冗余(重复的得分数值只保留一个)，并从小到大排列
cut -f 12 miniblast.text |sort -nu |wc -l11) 按照BitScore分值(第12列)的大小对整个文件进行排序(从大到小)
sort -k12nr miniblast.text > miniblast.sort12）找出比对最长的基因的ID (即QueryLen值最大)
cut -f 1,2 miniblast.text |sort -k2nr |head -113）将文件evm.model.scaffold_1.全部替换为gene
sed -i "s/evm.model.scaffold_1./gene/g" miniblast.text |less -S14）将第9列比对相似率(Identities%)大于90的行所对应的ID(第一列)输出来
awk '{if ($9 >90) {print $1}}' miniblast.text |less -S15） 将比对长度大于200(QueryLen)且比对相似率(Identities%)大于90所对应的ID(第一列)输出来
awk '{if ($2>200 && $9 >90) {print $1}}' miniblast.text |less -S16) 在 miniblast.text 第100行加入注释信息#no import genes
sed -i '100a #no import genes' miniblast.text事实上，这样的练习题可以是无数个，感兴趣的朋友可以在我们出题的基础上继续丰富，欢迎留言或者邮件跟我讨论，我的邮箱是 jmzeng1314@163.com