今天介绍DNA或蛋白质序列的两两比对(pairwise comparisons),除了具体操作,还加入了一些理论知识。
通过比对可以确定两序列是否同源,找到序列中的结构域,比较基因及其产物等,从而获得序列中的功能、结构、进化等信息。这次内容和之前一个推送内容有重叠,前期推送里简单介绍了全局、局部比对以及打分矩阵的概念,可以一看——【现学现卖】序列比对之算法。
上一章介绍的是序列的BLAST操作——【陪你学·生信】七、在数据库中检索相似的序列,结果返回很多序列。这时,如果想要详细对比列表中的序列与我们的查询序列,就需要两两比对。
一、准备对比的序列+选择适合的方法
1. 选择序列
(1)通过NCBI-BLAST,我们可以找到很多与查询序列比对结果比较好的序列,然后基于工具给出的identity,E-value值(DNA序列(至少100bp)identity>70%,E-value<10^-4;蛋白质序列(至少100aa)identity>25%,E-value<10^-4),选择可以继续进行两两比对的序列。因为这个指标一定程度上可以说明这两条序列是同源的,同源序列的结构和功能一般相似,这样接下来的比对才更有意义。
(2)当然两两比对的序列也不总是通过NCBI-BLAST结果,从数据库中进行选择。比如,实验过程中发现两条功能相似的片段,可以进行序列比对;或者实验中发现自己的序列有某些确定的功能,那么可以查找这个功能的关键词,下载相关功能的序列进行比对。
(3)还有可能只利用一条序列进行“两两”比对,这个操作可以发现序列中重复的区域,低复杂度的重复模体,回文序列以及RNA中潜在的二级结构。
2. 选择合适的方法
序列比对有不同的方法,适用于不同情况,主要分为两类,点阵法和动态规划算法,其中动态规划算法包括局部比对和全局比对,简介见下表。
方法
适用情况
Dot plot
查找序列重复区域;查找长插入/删除片段;提取部分序列进行多序列比对(看对角线方向,很直观)
Local alignments
通过高质量的,残基-残基的比对,分析两条序列的同源区域
Global alignments
在序列全长水平进行比对,鉴别长插入/删除片段;检查数据的质量;鉴别序列中的每个变异
二、点阵法(Dot Plot)
Dotplot 点阵法是非常简单、直观的一种两条序列比对的方法,它能展示出两条序列所有可能配对的区域,由研究者选择最有意义或自己感兴趣的区域进行后续研究。通过点阵法还可以对序列自身进行比对,寻找序列内部的正向或反向重复序列。
如果序列很短的话,点阵法是仅用纸笔就可以完成的两两比对。主要方法是将两条序列一条放在X轴从左至右,另一个放在Y轴从上至下。当对应位置的字符匹配时用“点”标记,最终形成矩阵。
如果两条序列完全相同,矩阵中会有一条长对角线;如果不完全相同,那么矩阵中就会出现不连续的对角线,如上图。
上图的是一个简单的例子,两条序列很短,如果是两条长序列比对,可想而知点阵图可能会凌乱和复杂。这时候,用滑动窗口(几个位置为一组)代替单个位点,可以更加显著的观察到两条序列的相似区域。
有Dotlet,Dnadot,Dotter,Dottup等专门做dot-plot的工具,可以设置的参数很多,比如窗口大小和相似度阈值这种的(窗口大小为10,相似度阈值为8,指的是每次比较时取10个连续的字符,如相同的字符超过8个,则标记)。我用的是geneious prime,一个很综合又方便的生信软件,进行两两比对结果里会出一个简单的,窗口大小为1的dot-plot图(付费软件,但是功能很强大,而且是实验室付费,嘿嘿,那就用好了)。
三、动态规划算法(dynamic programming algorithm)
从数学角度讲,我们可以把核苷酸/氨基酸残基看作字符,不同的匹配结果用给定的打分矩阵给分,因为序列长度有限,理论上比对方式有限,我们可以计算出所有的比对结果的得分,再找出得分最高的就是最优比对。这个思路叫做枚举法。不过工作量巨大,下面我们看看动态规划法。
动态规划用于在一个复杂空间中寻找一条最优路径。接着上面枚举法的思路,我们知道最终比对分数是各个残基比对分数之和,那么最好的比对结果,就是之前最好的比对+当前位置最好的比对。即局部最优解的组合就是全局最优解。
1. 全局比对(Global Alignments)
Needleman-Wunsch算法迭代公式
适用于序列整体水平相似程度较高的2个序列,Needleman-Wunsch算法是经典的全局比对算法。回溯时是从一端到另一端,整体回溯。以两条短序列AAG和AGC为例,看看两种比对的不同。举例中使用的打分矩阵如下:
打分矩阵
以AAG和AGC为例的全局比对
2. 局部比对(Local Alignments)
Smith-Waterman算法迭代公式
局部比对适用情况:两序列亲缘关系可能较远,但在局部具有相似性(比如有相似功能域,但是亲缘关系远的蛋白质序列);分析一条序列中的重复片段;内含子存在导致序列之间出现大片段差异。
Smith-Waterman算法是局部比对算法的基础。与全局比对的区别是公式里加了0,即在迭代时给分数加了下限(没有负值)。
所以回溯时也不一定是从一端到另一端,而是局部回溯。所产生的结果就是局部比对结果,与全局不同的是,除了最优比对结果外,还会产生次优比对结果。
以AAG和AGC为例的局部比对
四、使用Lalign进行比对(以局部比对为例,全局比对操作一样的)
1. 网址
https://embnet.vital-it.ch/software/LALIGN_form.html
2. 操作
放大操作页面解释各部分:
(1)选择比对类型,这里我选了局部比对。
(2)如果是选择局部比对才会有返回比对结果个数的选项,我写了10个,所以结果这里会从最优到次优给我列出10个不同局部比对结果。
(3)E-value阈值默认10,即随机产生10次这样的结果。E的阈值越小,比对越严格。
(4)打分矩阵,可选的有DNA,RNA打分矩阵和蛋白质打分矩阵PAM,BLOSUM。一般选BLOSUM62,这里默认50,也差不多。BLOSUM矩阵后数字越大,比对越严格,越适合相似程度高的两条序列;而PAM趋势正好相反。
(5)产生空位罚分,对应BLOSUM50的默认值为-12。没有最优的值,如果我们设定空位产生时罚更高的分,那么局部比对结果将更加局部和分散,产生很多小的比对。
(6)延长空位罚分,对应BLOSUM50的默认值为-2,绝对值一般为产生空位罚分的十分之一左右。当我们比对两条亲缘关系较远的序列时,设定“产生空位罚分”相对高,而“延长空位罚分”相对更低会产生好结果,这表示空位罚分更看重空位的产生而非空位长度。
(7)(8)(9)就是输入序列啦,可以复制粘贴,也可以(8)选择数据库,然后在(9)中输入序列在该数据库中的序列号。
3. 结果
结果将返回10个不同的比对,会给出打分和E-value,看E-value更能说明比对质量,当然也要同时考虑比对的长度。
五、网上免费比对工具
1. 下面有一些比对工具链接:(括号里面是比对类型)
(1)Blast2seqs(Local BLAST)
https://blast.ncbi.nlm.nih.gov/Blast.cgi?BLAST_SPEC=blast2seq&LINK_LOC=align2seq&PAGE_TYPE=BlastSearch
(2)Lalign(Global/Local)
https://embnet.vital-it.ch/software/LALIGN_form.html
(3)xenAliTwo(Local for DNA)
https://users.soe.ucsc.edu/~kent/xenoAli/xenAliTwo.html
(4)Pal2nal(Protein against DNA)
http://www.bork.embl.de/pal2nal/
2. 有时候比对结果可能需要进一步分析,下面是推荐:
(1)LalnView——比对结果可视化软件
下载链接
http://doua.prabi.fr/software/lalnview
(2)PRSS(embnet)——评估比对的显著性
https://embnet.vital-it.ch/software/PRSS_form.html
(3)PRSS(virginia)——评估比对的显著性
https://fasta.bioch.virginia.edu/fasta_www2/fasta_www.cgi?rm=shuffle
往期相关内容:
【陪你学·生信】五、当你有一段待分析的DNA序列(基础操作介绍)
【陪你学·生信】六、当你有一段待分析的氨基酸序列(基础操作介绍)