【原】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./wiki/BLAST#Output

06

练习题

首先需要自行产生一个 blast.test 文件，公众号教程推文不方便展示。

1) 查看 blast.test 文件的前6行

1. head -6  blast.test

2) 统计blast.test文件中共有多少行，多少列

1. wc -l blast.test

2. head -1 blast.test |tr '\t' '\n' |wc -l

3）新建一个文件miniblast.text

1. touch miniblast.text

4）将blast.test文件前500行内容写入到前一步新建的miniblast.text 文件

1. head -500  blast.test > miniblast.text

思考一下，有没有可能不需要新建miniblast.text 文件

5) 删除原文件 blast.test

1. rm blast.test

6) 查看QueryId为evm.model.scaffold_1.187在miniblast.text文件中的第几行

1. grep -n evm.model.scaffold_1.187 miniblast.text |less -S      

2. # less  -SN miniblast.text

7） 查看ribosomal在miniblast.text文件中出现的次数

1. grep -c ribosomal miniblast.text

8）查看miniblast.text文件中不含有protein的行

1. grep -v protein miniblast.text |less -S

9）单独查看miniblast.text文件中第12列BitScore中所有的数值

1. cut -f 12 miniblast.text |less -S

10) 将miniblast.text文件中第12列BitScore中所有的数值去冗余(重复的得分数值只保留一个)，并从小到大排列

1. cut -f 12 miniblast.text |sort -nu |wc -l

11) 按照BitScore分值(第12列)的大小对整个文件进行排序(从大到小)

1. sort -k12nr miniblast.text > miniblast.sort

12）找出比对最长的基因的ID (即QueryLen值最大)

1. cut -f 1,2 miniblast.text |sort -k2nr |head -1

13）将文件evm.model.scaffold_1.全部替换为gene

1. sed -i "s/evm.model.scaffold_1./gene/g" miniblast.text |less -S

14）将第9列比对相似率(Identities%)大于90的行所对应的ID(第一列)输出来

1. awk '{if ($9 >90) {print $1}}' miniblast.text  |less -S

15） 将比对长度大于200(QueryLen)且比对相似率(Identities%)大于90所对应的ID(第一列)输出来

1. awk '{if ($2>200 && $9 >90) {print $1}}' miniblast.text  |less -S

16) 在 miniblast.text 第100行加入注释信息#no import genes

1. sed -i '100a #no import genes'  miniblast.text

事实上，这样的练习题可以是无数个，感兴趣的朋友可以在我们出题的基础上继续丰富，欢迎留言或者邮件跟我讨论，我的邮箱是 [e-mail]