【5.1】blast原理介绍(序列数据库的搜索)

Homology is the central concept for all of biology.”——David Wake. Science, 1994

博主跟实验室交流的时候,被问得最多的就是blast这个参数怎么设,这些参数有什么意义?作为一个‘资深的生物信息工程师’,有必要好好梳理blast究竟都干了一些什么。

一、 为什么blast

BLAST,Basic Local Alignment Search Tool. 通过序列的相似性,我们能知道其在结构,功能,进化地位的相似性。通过用一个未知的查询序列,在包含已知序列及注释的数据库中检索,已推断查询序列的功能,结构以及演化信息。具体来说就是将未知序列与已知序列的数据库逐一进行比对,并根据统计上显著的hit序列所对应的注释,来推测查询序列的结构,功能以及演化。

我们可以想象到,只要在之前博客中提到的双序列比对上加一个for循环就可以了。因为数据库中庞大的数据,这个算法所用的时间会比较长。我们发现,在迭代比对和全局比对中,最优比对的最优路径在主对角线小范围的区域,而在局部比对过程中,我们发现最优比对对应的最优路径平行于主对角线。如果我们事先知道最优比对所对应的最优路径在表格中所对应的位置,我们就能直接填空,有效的降低需要填写的格子数目,从而减少计算量。这就是blast的基本原理与思路,该方法自1990年首次提出后,是目前被广泛认可和使用的算法之一。

二、如何blast

在NCBI中进行序列比对, 粘贴完querry sequence 后,选择数据库。第一个nr数据库,也就是所有蛋白的去入围结合,这个数据库一般用在数据的筛选中,以确定查询之前该序列是否被人研究过,虽然内容丰富,但并非每条序列都有详尽的注释。

我们选择的SWISSPROT数据库是UniProt数据库中的一部分,其中的每一条序列都经过专家团队的手工注释,包含了功能,修饰以及结构全方面的信息以及到其他相关数据库的链接。如果需要对新序列尽可能详尽和准确进行注释的时候,通常会用SWISSPROT数据库最为入手数据库。然后选择物种。点击blast 这个过程比较傻瓜式,多点击点击就明白怎么用了。

更多数据库的介绍,见专题

三、blast的原理

blast首先找到两条序列的高度相似小片段,也就是所谓的种子,而后以此为基础向两端延伸,并构建比对,最后为了避免假阳性,还需要计算统计显著性。

3.1 将输入序列切割成若干小段,成为seed word,也就是种子单字

 -w  种子长度
长度为n的序列,种子个数为n-w+1

3.2 然后根据已经建立好的索引表

快速定位其在数据库中候选序列的位置(据说设置好索引结构,这一步可以在线性甚至常数的时间内完成,不明白啊!每日一生信–序列数据库的搜索(blast)

3.3 通过对所有的SEED重复上述的操作

我们可以获得查询序列和比对序列的hit-map,根据前面讨论,我们知道最优比对应该平行于主对角线,因此我们可以去掉那些零散的HITS,而只容许沿对角线的两个或两个以上的hits所组成的Hit Clusters,从而减少搜索空间

3.4 然后我们以HIT cluster 向两边延伸

直到总分数的下降超过一个指定的值X后停止,在扩展的区域我们可以用我之前博文提到的动态规划从而显著的降低计算量。

3.5 为了提高blast的速度和灵敏度,还有一些其他的技巧。

首先我们需要屏蔽重复性的低复杂性的片段,从而减少假阳性。这里的低复杂度是根据序列的信息量来判断。没有这个SEED是不能用的,因为有太多的种子。

-F parameter  完成过滤
Ns for nudeotide residues
Xs for amino acid residues

例如下面的,我们可以屏蔽掉K=0.36,以免影响后面的结构

另一方面,为了提高灵敏度,在seeding的时候,除了考虑由查询序列分离而成的种子单字,我们还会考虑与种子单字相识的临近单字。根据替代矩阵来计算分数,例如以DKT作为种子单字,那DRT就是6+2+5=13。我们可以计算出其他的分数,为了避免其他的假阳性,我们只考虑其高度相似的临近单字。在当前版本中要求其分数大于11。

3.6 一个aa跟另一个aa随机匹配上的概率是1/20,而如果有一段6个aa的序列,则随机匹配上的概率是(1/20)^6,在一个有192206270个序列的数据库中,就会有3个是百分之百的随机匹配上。因此我们需要有一个方法来客观的评价统计显著性。所以我们在获得一个比对结果后,需要确定其统计显著性,确保这个比对不是由随机性引起的。我们引入了E-value这个值,E-value指在随机情况下获得跟当前比对相等或高于当前得分的比对条数。具体来说,如果一个比对结果的E-value值为10,则意味着有10个随机的比对等于或高于当前比对的得分。他不是一个概率,而是一个期望,所以可以大于1。m是querry sequence的长度,n数据库的大小,k,λ与你打分矩阵有关系,s是分数。

总之,blast是一种启发式的算法,他并不能确保一定能找到最优解,但尽力在最到的时间类找到最优解,然而速度的提高是以牺牲灵敏度为代价的。

四、说明

4.3 两种blast

Gapped BLAST:
Adopts a lower neighborhood word score threshold to maintain
the same level of sensitivity for detecting sequence similarity.
Gaps allowed(容许gap的存在,更好的表示序列相似程度,中间没有中断产生)
 
PSI-BLAST:(位置特异性迭代blast)
Position-Specific Iterative BLAST (blastpgp)
Constructs PSSMs (position specific scoring matrix) automatically
Searches protein database with PSSMs
Used to find distant relatives of a protein, and is much more
sensitive in picking up distant evolutionary relationships than a
standard protein-protein BLAST.
(高度比对之后的序列,我们重新设打分举证,然后第二轮比对,高度保守的蛋白序列得分高)

4.4总结

BLAST is based on the similarity statistics, but statisticalsignificance doesn’t mean biological significance.blast是基于统计学上的相似性,但统计显著性不代表生物学上的显著性。

参考资料

  • 北大高歌老师的课件
药企,独角兽,苏州。团队长期招人,感兴趣的都可以发邮件聊聊:tiehan@sina.cn
个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn