上个BBQ介绍了基于动态规划算法的全局序列比对和局部序列比对，今天这篇笔记简单的介绍一下什么是BLAST和BLAT~
(写的比较基础，其实还有很多可以写写，比如说打分矩阵PAM矩阵，BLOSOME矩阵等等，这些矩阵怎么来的，如何选用这些矩阵等等，欢迎一起沟通交流！)
--------------------------------------------我是分割线-----------------------------------------------------
一：我们来看一个有意思的现象：
孟叔的BBQ里有这么一个例子：
假设我有1条序列 SeqA = 100bp（这个不是很长哟~），我想找到这条序列的相似序列。
现在服务器最快的CPU，单核心1秒可以大约运行3×10^9次，假设我们的服务器有20个核心，那么 完成任务大约需要的时间 = 471930 × 10^9 ÷ 3 × 10^9 ÷ 20 = 7860秒 大约是130分钟；
我们就为了找个相似性，就要花2个小时？这个肯定不划算吧？ 所以，有了这个需求以后，一帮大佬就搞出个优化算法叫BLAST（Basic Local Alignment Search Tool）把上面这个任务的运行时间缩短到了1分钟以内。
（ps:目前已知的非冗余核酸序列库，序列有47193206条，按平均长度在0.1Mbp左右。如果我们想要找到SeqA的相似序列，使用局部比对算法，那么至少就需要100bp × 0.1Mbp × 47193206条 次比较运算等于471930 × 10^9。）
1.1减少序列比对时间的核心是什么?
如何减少比对计算的次数！
1.2.思考一下：如果是你的话你会怎么去解决这个问题呢？如果是一段很长的序列的话，那么如何进行比较呢？采用全局的算法从头比到尾吗？
（1.0版本）如果是我的话我估计会采取先把长的序列分成小的序列，然后对每个小的序列进行序列比对。但是问题来了，如果是切成小的序列比对的话，会增加比对的次数，肯定会增大计算量。
（2.0版本）如果是取一小段先去比对呢？？？估计有戏呀
-----------------------------------------------分割线君上线--------------------------------------------------------
对于BLAST算法来说，其核心实现就是先寻找到一个“种子序列”，然后在种子序列的基础上再进行两端的延伸。这样子可以减少计算的次数，然后又可以提高效率，如果比的上就延伸，比不上就放弃。非常的高效了呀！！！！
BLAST算法思想
1.2采取seed的方法，我们来看下面一个小栗子：
举个例子，对于一个具有6个氨基酸长度的seed来说，如果匹配到完全一样的序列的概率是1.56*10^-8次方。看上去这个概率是非常的小了。（因为氨基酸有20种，所以是1/20的6次方），但是往庞大的数据库里一比较，就会发现能匹配出来的有3个（SWISS-PROT数据库）。所以对于一个比对的结果我们如何进行比对结果的判断呢？
1.3 需要解决的问题：随机匹配的出来的结果如何排除假阳性？通过什么指标呢？
在BLAST里面，采取E-value去评判随机性，简单来说，E-value是指在当前比对情况下获得相对比相等或更高的分数的条数。就拿上面的例子来说，E-value就是3，E-value是大于等于1的。当数据库越大时，随机匹配的可能性也会越大。而且当seed的数量越多时，E值也会大！
而且通过以下的表格我们发现，当P=0.05时，E=0.0513,所以可以将P=0.05当成是E的cutoff值
P值和E值的关系
但是！
这种根据seed的方式，只在有限的区域进行BLAST算法，这种可以提高速度，但是这种速度的提高是以灵敏度为代价的！
二：Question & Answer
2.1:BLAST的算法核心是什么？
启发式算法
2.2 BLAST结果里的E-value是什么意思？
它表示了在随机的情况下出现相同结果的条数，是大于等于1的
（S-score是指序列比对的结果打分值，S的得分越大越好，说明相似度高）
2.3BLAST结果里的P-value是什么意思？
出现假阳性的概率，一般要求p<0.05
2.4：如果要降低假阳性，可以从哪些方面入手呢？
1）seed的长度不能太小，如果seed长度太小的话，假阳性会高
2）seed的数量不能太多
3）去掉低复杂的重复序列
2.5. 假设给你一条序列，运行结果中序列相似度最高的来自于哪个物种？
Protein Sequence
MVRAPCCEKMGLKKGPWTPEEDQILISYIQSNGHG
NWRALPKLAGLLRCGKSCRLRWTNYLRPDIKRGNF
TREEEDSIIQLHEMLGNRWSAIAARLPGRTDNEIK
NVWHTHLKKRLKNYQPPQSSKRHSKNKDSKAPCTS
QIALKSSNNFSNIKEDGPGLGSGPNSPQLSSSEMS
TVTADSLAVTMDISNSNDQIDSSENFIPEIDESFW
TDGLSTSGGGEELQVQFPFHDMKQENVEKDVGAKL
EDDMDFWYSVFIKSGDLLELPEF
使用网站
BLAST：http://blast.ncbi.nlm.nih.gov
参数设置：
Database: Non-redundant protein sequences (nr)
Algorithm: blastp
Word size: 3
Matrix: BLOSUM62
Gap Costs: Existence: 11 Extension: 1
其他参数默认
BLAST结果
三：什么是BLAT？它和BALT有什么区别？
BLAST = Basic Local Alignment Search Tool；
BLAT = BLAST-like alignment tool；
为什么要进行BLAT?进行BLAT的目的是什么？
BLAT的功能，简单来说就是我有1条序列SeqA，我想知道SeqA在某一种基因组中的定位（比如human genome）的工具。
UCSC -BLAT
UCSC-BLAT2
谈谈个人的看法：一般在利用BLAT的时候，我们都是在设计引物的时候检测一下引物的特异性的。如果设计的引物如果匹配到了基因组上的很多地方，那么说明引物的特异性不好。需要重新设计。
具体的实战问题
我们对这两个序列进行了BLAT定位，发现确实可以在UCSC上比对到较多的结果，然后我发现这两条序列都比对到了5SRNA上面。
比对结果一
比对结果二
对这两个序列进行了比对，发现这两个序列比较类似
比对结果三
解释：这两条序列有高度的相似性，而且都匹配到了核糖体RNA上面。有可能是在建库的过程中，核糖体RNA去除的不干净导致。
思考：如何去除呢？这种建库带来的核糖体RNA-seq的结果呢？欢迎一起讨论~
Reference:
生物信息学100个基础问题 —— 第14题 你应该知道的BLAST
生物信息学：导论与方法(北京大学) - BLAST
网友评论