blast 算法初探
Blast 算法自从1990以来被引用了上万次，是basic local alignment search tool的缩写
blast
注：先文字+后图片的行文方式
提高计算速度的关键是：围绕最优路径计算
image.png
image.png
具体来说，blast就是先找到两条序列之间高度相似的小片段，也就是所谓的种子seed；而后向两端延伸并比对；最后为了避免假阳性，blast还会计算统计显著性。
首先，blast将输入片段切割成若干小段，称之为seed words；
image.png
然后，快速在database中搜索和seed word相似的序列；以及在候选序列中的具体位置。
image.png
通过对所有seed words重复上述操作，就可以得到查询序列和候选序列之间的hit map
image.png
最优比对路径应当平行于主对角线，
image.png
因此我们可以去掉那些零散的hits，
image.png
只留下>=2个与主对角线平行的hit cluster，以便后续进一步缩小比对范围
image.png
接着我们可以以这几个cluster为中心向两端延伸扩展，
image.png
image.png
直到总分数的下降超过一个给定的值x之后，（没懂，简单理解成到一定范围停止向两端扩展）
image.png
在扩展区域，应用动态规划算法比对，已确定最终的比对区域（图中应该是smith-waterman算法），从而显著地降低计算量
image.png
为了提高速度和灵敏度，blast还采取了其他的一些技巧，
为了提高速度
首先，blast会屏蔽重复性的低复杂度区域，以免产生太多的假阳性hit，
image.png
这里的复杂度通常是根据序列的信息量来判断，
image.png
不屏蔽的话，seed是不可用的，因为会产生太多的假阳性种子seed，不可能从众多的比对hits中去找假阳性的hit，应当在比对之前屏蔽这些low-complexity的序列，根据下面这个k值屏蔽。
image.png
例如，对于一个典型的微卫星序列,
微卫星DNA：重复单位序列最短，只有2～6bp，串联成簇，长度50～100bp，又称为短串联重复序列（Short Tandem Repeat STR)。广泛分布于基因组中。
例如微卫星序列CACACACACACACACA，我们可以简单地计算出他们的K值大小，K=0.36，可以屏蔽掉，以免影响后续的搜索。
image.png
为了提高灵敏度
为了提高灵敏度，blast在seeding的时候，blast除了考虑由查询序列分解的种子单字 seed words以外，还会考虑那些和种子单字相似的邻居单字neighbourhood words；
具体来说，就是根据seed words所有可能的变形计算替代矩阵的分数，
image.png
例如这里以DKT作为种子单字，那么DRT的分数就是，6+2+5=13
image.png
类似的，我们可以计算出其他和DKT相似的序列的分数；为了降低假阳性，我们只考虑那些和seed words高度相似的，neighbourhood words，
image.png
具体来说，在当前版本中，分数大于等于11才予以考虑，
image.png
实例
对于一个长度为L的蛋白序列，有（1/20）^L的概率match完全一样的序列；假定你的蛋白序列长度，为6个氨基酸，则概率为：
image.png
那么用这个蛋白序列在swiss-prot序列库中检索，随机匹配的情况下，百分之百匹配的期望是3，
image.png
我们需要一个客观的方法来统计比对显著性，以确保这个比对不是由随机因素引起的，
image.png
在blast中使用e-value来对显著性进行度量，
e-value：指的是在随机情况下，获得比当前比对分数，相等或者更高分数的可能比对条数。
例如，一个比对的e-value=10，就意味着会有10个随机的匹配获得比当前比对相等或者更高的分数。
image.png
需要注意的e-value不是一个概率，它是一个匹配数目的期望（期望值E=随机匹配获得比当前比对相等或者更高的分数的条数），E是>=1的，是根据下面这个公式算出来的，
image.png
m是query sequence的长度；
n是你数据库的大小；
e是自然对数；
S是分数（没懂是啥分数）
K和λ是跟打分矩阵相关的，相当于一个normalization的factor
从这个公式可以看出，E的大小和数据库的大小n是呈正比的；
也就说数据库越大，随机匹配的可能性也就越大，这与我们刚刚看到的例子是相符合的；
另一方面，E值得大小和查询序列的长度m也呈正比（m越大，seed越多，所以正比），这是因为blast是局部比对，不需要全长的匹配；
E值与比对的分数S负相关，也就是说分数越高随机遇上的可能性也就越小；
同时公式中的λ和K是与打分系统和搜索矩阵相关的修正值，用来平衡不同打分矩阵和搜索空间对于结果的影响；
image.png
为了便于解释呢，我们把P值（概率值）和E值可以进行相互转换；
在P值和E值分别小于0.1的时候，二者几乎相等；
特别地，当P取0.05的时候，对应的E value是0.0513，因此有人也常常有人将E value 0.0513作为E value 的 cut off
image.png
blast和simith waterman、niddleman wunsch算法的区别
blast它并不能确保找到最优解，但是能保证在最短时间内找到足够好的解;
只在有限的区域运用动态规划算法，从而有效地降低计算量，提高计算速率，然而速率的提高是以灵敏度下降为代价的，这也是其他的算法固有的tradeoff。
image.png
image.png
思考
How BLAST works
image.png
image.png
image.png
image.png
image.png
高于T的word hits被保留下来
image.png
image.png
HSP
HSP图：如果序列打分低于S就停止extending，高于S就保留下来，保留下来的序列就是high scoring segment pairs，这就是blast里面最重要的这个HSP的概念。
评估blast结果，用E-value
image.png
HSP的打分服从泊松分布，意义在于，我们数据库里面其他的序列和我们的query sequence要比对的序列相似性程度要高于我们现在要找到序列的相似性程度的概率（没懂）；
E-value 越小越好，如果E-value大于1的话说明我们找到的序列是随机发生的事件，是不可靠的结果；
E-value < 0.1 or 0.05 的话代表是在统计学上有意义的一个结果；
E-value < 10^(-5) 说明我们查询到的序列和我们用来比对的序列是高度一致性的，有非常高的相似度，所以在得到blast结果的时候要特别关注E-value。
image.png
image.png
image.png
Gapped BLAST 允许gap的存在，更好的比对相似程度，中间没有中断的产生（没懂）
image.png
位置特异性的迭代的blast，创新在于第一轮blast以后，以第一轮高度blast高度相似的序列去做一个打分矩阵，以重新设定的打分矩阵进行第二轮的blast，然后以此类推进行重新的搜索，重新blast，然后重新得到一个打分矩阵。
PSI-Blast在蛋白质比对过程当中用到的比较多，如果是一个高度保守的蛋白的话，它的得分是比较高的，然后如果相似度偏低的话，它的得分就越趋近于0。
blast只是比较序列相似度的一个基于算法的一个统计学工具，使用的时候需要注意：
1.统计的结果不能完全代表生物学上的意义，要考虑到结构上的一些问题；
2.如果有蛋白质序列的话最好先比对蛋白质的相似性，再去比对序列的同源性；
3.特别注意的是，相似性的概念不能代表同源性的概念；通过blast比对的，相似性很低的两个序列不能得出同源性很低的结论，不能类推；没有百分之几十的同源性这种说法，这种表达是不科学的。
image.png
image.png
image.png
网友评论