十分钟吃透blast比对底层逻辑并实操






十分钟吃透blast比对底层逻辑并实操

小果  生信果  2023-10-03 19:00:14

相信做生物的小伙伴都接触过序列比对,不论是想知道测序得到的序列还是你手中的由ATCG组成的一条序列想是什么,都需要进行进行序列的比对,才能确定它的身份。

序列比对是一种用于比较和分析序列之间的相似性和差异性的方法。它在生物学和计算机科学领域被广泛应用于DNA、RNA和蛋白质序列的比较分析。

序列比对简介

序列比对的原理基于两个基本假设:保守性和演化性。保守性假设认为在不同物种或个体之间的相似序列可能在进化过程中被保留下来,而相似的序列可能具有相似的功能。演化性假设认为序列之间会发生变异和突变,通过比较序列之间的差异可以推断它们之间的共同祖先和演化关系。根据比对的目的和序列的长度,选择适当的比对算法。常用的比对算法包括全局比对算法(如Needleman-Wunsch算法)和局部比对算法(如Smith-Waterman算法)。

需要注意的是,序列比对是一个计算复杂度较高的任务,特别是当处理较长的序列或大规模的多序列比对时。因此,研究人员通常利用高性能计算和优化算法来加快比对的速度和提高比对的准确性。



序列比对定义

简单的进行匹配,可以看到仅有两个碱基匹配

如果引入一个空格或空格(gap),匹配上的状况就可以明显改善

因此,小果根据查阅文献认为以下说法是序列比对的最好解释。

根据特定的计分规则,通过一定的算法对两条或者多条DNA或者蛋白序列进行比较,找出他们之间的最优匹配或者最大相似度匹配。序列匹配的意义在于,通过序列比对可以获得一个序列相似性度比值,通过这个值可以获得两条序列相似度或者同源性的统计学判定结果。

动态规划

想要了解blastn或者序列比对的底层逻辑,我们需要把动态规划弄明白,讲清楚。


打分矩阵


首先,需要了解什么是打分矩阵。

对碱基匹配情况进行一个标准化的矩阵,可以根据

1、替代矩阵:

DNA替代矩阵:用于评估DNA序列中不同碱基之间的替代得分。常见的DNA替代矩阵包括简单的等分矩阵(例如,A和T之间得分为1,A和G之间得分为-1)以及更复杂的模型,如PAM(Point Accepted Mutation)和BLOSUM(Blocks Substitution Matrix)矩阵。

蛋白质替代矩阵:用于评估蛋白质序列中不同氨基酸之间的替代得分。常见的蛋白质替代矩阵包括BLOSUM系列和PAM系列矩阵。这些矩阵根据蛋白质家族的多序列比对结果建立,考虑了氨基酸的进化关系和替代频率。

2、缺失和插入罚分:

Gap罚分矩阵:用于惩罚序列比对中的缺失(Gap)和插入。引入Gap会导致序列不匹配和不对齐,为了惩罚Gap的引入,通常在打分矩阵中引入了开放(Gap Opening Penalty)和扩展(Gap Extension Penalty)处罚得分。这些得分会随着Gap的长度逐渐增加。

3、直接得分矩阵:

匹配-不匹配得分:在某些情况下,打分矩阵中的得分可以根据碱基或氨基酸之间的匹配或不匹配来定义。例如在DNA序列比对中,匹配得分可以设为正值1。

这里对各种矩阵不理解可以直接跳过,直接了解下边例如中的两个矩阵。

例如:下边两个矩阵,第一个只考虑匹配上的得分为1,第二个只考虑了考虑颠换和置换的情况(序列中嘌呤被嘧啶替换或反之,称为颠换(transversion),如A→C、G→T;嘌呤或嘧啶自己互换称为置换(transition),如A→G、C→T。可对核苷酸的颠换与置换给予不同的评分。)在序列比对中打分矩阵可以根据需求进行自定义~




比对算法

全局比对:给定两个序列A=a1a2…an 和B=b1b2…bm,S(i,j)表示两个序列所有比对的最好分数。在设定初值后可以用以下递归关系计算:



局部比对的具体算法与全局比对几乎一样,给定两个序列A=a1a2…an 和B=b1b2…bm,S(i,j)表示两个序列任何比对的最好分数。在设定初值后可以用以下递归关系计算:



上述的算法中,不论是全局比对还是局部比对都是来寻找最大值的过程,我们需要得到Mij的分数,是取了三个方向中的最大值。


例如:现在还是用开头时候我们用的序列,我们现在用动态规划的方式来得到下图中的匹配方式:



规定的打分矩阵为:匹配得分5分;空位罚分-1分;错配 -3分

小果建议可以在excel尝试自己算出以下小果填入的数字。

这里注意,是mismatch还是gap,如果你是从对角线移动的为错配,从上边或者右边的格子移动的即为gap,这两种情况的罚分是不一样的。


得到上述格子后,查阅到得分最高的格子,追溯到是怎么来的序列即可。


blast比对原理

1、编译一个由查询序列生成的长度固定的字段编译列表“words”

2、在数据库中扫描获得与编译列表中的字段匹配的序列记录;

3、以编译列表中的字段对为中心向两端延伸以寻找超过阈值分数S的高分值片段HSP。

4、对于每部分HSP使用Smith-Watermans算法进行局部性比对,为每部分打分,作为最终结果。

下图,是blast的文章中的原理

E值衡量了在随机情况下,数据库存在的比当前匹配分数更好的比对的数目,因此可以用该值作为指标评价HSP比对序列的可信度,E值计算公式:E=Kmne-λS

了解了blast的原理后,下面介绍一些blast的应用和参数。


blast比对实际操作

比对

blastn -query input.fasta -db database -out output.txt -evalue 0.001 -outfmt 6blastn -query input.fasta -db database -out output.txt -outfmt " 7 qaccver saccver pident length mismatch qlen qstart qend sstart send qseq evalue bitscore" -num_threads 10 -word_size#也可以尝试7的格式。


在这个示例中,blastn是BLAST程序的名称,input.fasta是输入文件,database是BLAST数据库,output.txt是输出结果文件,-evalue参数指定过滤结果的期望值阈值,-outfmt指定输出格式为表格格式(6)


构建数据库

makeblastdb -dbtype nucl  -in ${name}.fa -out ${name}#完成后会显示以下提示


#会生成以下的几个文件



参数解读


1、E-value(期望值):

E-value是一个衡量比对结果的期望随机发生次数的指标。它表示在进行大量随机比对的情况下,期望得到与查询序列的比对相似性至少一样好(或更好)的结果的次数。

较小的E-value表示比对结果的相似性更显著和显著。通常,E-value小于某个预设的阈值(例如0.001或0.05)被认为是显著的比对结果。在BLAST输出中,E-value一般以科学计数法的形式给出,例如1.0e-10,表示10的负10次方(即0.0000000001)。

2、Score是基于比对中的匹配和不匹配进行计算的一个值,用于衡量比对的质量和相似性。具体的得分计算方式取决于所使用的打分矩阵。一般来说,较高的Score表示比对结果具有更好的相似性和匹配程度。通常情况下,Score越高,比对的质量越好。在BLAST输出中,Score是以具体数值的形式给出。

3、“Identity”(相同性)是输出结果中的一个重要指标,用于描述比对序列的相似性程度。它表示在比对过程中,两个序列中存在相同的碱基或氨基酸的百分比。相同性(Identity):在比对结果中,相同性通常以百分比的形式表示,反映了两个序列中相同位置的碱基或氨基酸的百分比。较高的相同性表示比对结果具有更高的相似性。相同性的计算:相同性的计算是通过比对过程中的匹配(Matches)和不匹配(Mismatches)数量来得出的。计算公式:Identity(相同性)= (Matches / Alignment Length) × 100%

其中,Matches表示比对中完全匹配的位置数,Alignment Length表示比对的序列长度。

例如,假设有一个比对结果,序列的长度为100,其中有80个位置是完全匹配的,那么相同性为:

Identity = (80 / 100) × 100% = 80%

请注意,相同性只考虑了完全匹配的位置,不考虑缺失(Gaps)和插入(Insertions)。缺失和插入会影响相似性和比对质量的其他方面,例如Gap Coverage和Gap Length。

Identity、E-value、Score、Coverage等。这些指标的综合分析可以提供更全面的了解和解释比对结果的相似性和可靠性。


小果附上blasn命令的详细解释:

-query <filename>:指定输入文件,其中包含待比对的查询序列。

-db <database_name>:指定比对所使用的数据库。

-out <output_file>:指定输出文件的名称或路径。

-evalue <evalue>:设置过滤结果的期望值(E-value)阈值。默认值为10。

-outfmt <format>:指定输出结果的格式。常用的格式参数包括:

0:Pairwise(默认格式)

5:XML格式

6:Tabular格式

7:Tabular格式(注释)

更多格式参数可以通过blastn –help查看。

除了上述常见的参数外,还有一些其他参数可以进一步调整BLAST比对的行为:

-num_threads <num>:设置使用的线程数(并行处理)。默认值为1。

-task <task_name>:指定BLAST任务类型,例如blastn、blastn-short、blastn-fast等。不同的任务类型针对不同的查询和目标序列类型进行了优化。

-word_size <int>:设置匹配单词的长度。默认值为11。

-gapopen <int>:设置开启一个gap(间隙)的惩罚分数。默认值为5。

-gapextend <int>:设置扩展一个gap(间隙)的惩罚分数。默认值为2。

-max_target_seqs <int>:设置从数据库中返回的最大比对序列数。默认值为500。

-perc_identity <float>:设置比对结果中的最低相似度百分比。默认值为0。


好了小果今天的分享就到这里了,小伙伴们,今天有没有学到新知识呢,想要继续了解R语言内容可以持续关注小果哦~或者也可以关注我们的官网也会持续更新的哦~

往期推荐

1.搭建生信分析流水线,如工厂一样24小时运转Snakemake——进阶命令
2.比blast还优秀的序列比对工具?HMMER来了
3.对单细胞分析毫无头绪?让popsicleR领你入门
4.小果带你绘制ROC曲线评估生存预测能力