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.test

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

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

3)新建一个文件miniblast.text

touch miniblast.text

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

head -500  blast.test > miniblast.text

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

5) 删除原文件 blast.test

rm blast.test

6) 查看QueryIdevm.model.scaffold_1.187miniblast.text文件中的第几行

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

7) 查看ribosomalminiblast.text文件中出现的次数

grep -c ribosomal miniblast.text

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

grep -v protein miniblast.text |less -S

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

cut -f 12 miniblast.text |less -S

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

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

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

sort -k12nr miniblast.text > miniblast.sort

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

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

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

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

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

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

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

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

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

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

事实上,这样的练习题可以是无数个,感兴趣的朋友可以在我们出题的基础上继续丰富,欢迎留言或者邮件跟我讨论,我的邮箱是 jmzeng1314@163.com

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-09-10

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏儿童编程

Python基础语法导图汇总(补充篇)

这一段间经过对Python的学习,在之前两次汇总后,又对遗漏的部分及新接触的内容用思维导图的形式进行了汇总。

1313
来自专栏PaddlePaddle

PaddlePaddle发布新版API,简化深度学习编程

PaddlePaddle是百度于2016年9月开源的一款分布式深度学习平台,为百度内部多项产品提供深度学习算法支持。为了使PaddlePaddle更加易用,我们...

3558
来自专栏专知

【干货】主题模型如何帮助法律部门提取PDF摘要及可视化(附代码)

【导读】本文是Oguejiofor Chibueze于1月25日发布的一篇实用向博文,详细介绍了如何将主题模型应用于法律部门。文章中,作者分析了律师在浏览大量的...

3737
来自专栏机器学习实践二三事

Caffe中均值文件的问题

关于均值文件 (1) 在Caffe中作classification时经常需要使用均值文件,但是caffe自己提供的脚本只能将图像数据转换为 binarypr...

2149
来自专栏青青天空树

取随机数

  常用于去随机数的函数为rand()(在stdlib.h头文件中,不同的编译器可能有不同),但是实际在使用这个函数时却发现每次程序运行产生的数都是一样的,这是...

1482
来自专栏机器之心

教程 | 如何用百度深度学习框架PaddlePaddle做数据预处理

3426
来自专栏大数据挖掘DT机器学习

R语言进行中文分词,并对6W条微博聚类

由于时间较紧,且人手不够,不能采用分类方法,主要是没有时间人工分类一部分生成训练集……所以只能用聚类方法,聚类最简单的方法无外乎:K-means与层次聚类。 尝...

3705
来自专栏和蔼的张星的图像处理专栏

LCT代码跑起来先文章思路总结

论文才刚开始看,但是代码先跑了一下看结果,有一点小坑,记录下: 首先去论文的github上去下载代码:点这里 readme里其实写了怎么搞:

4843
来自专栏大数据智能实战

LargeVis可视化技术学习

大图可视化一直是大数据可视化领域的一个关键技术,当前有各种办法,但是今年出来了一个LargeVis的技术,因此对这个技术进行复现和学习一下。 前面有很多基础理论...

3977
来自专栏用户2442861的专栏

Caffe学习:Blobs, Layers, and Nets

-注意:网络结构是设备无关的,Blob和Layer=隐藏了模型定义的具体实现细节。定义网络结构后,可以通过Caffe::mode()或者Caffe::set_m...

640

扫码关注云+社区

领取腾讯云代金券