【陪你学·生信】八、序列两两比对的理论和操作
?今天介绍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序列（基础操作介绍）
【陪你学·生信】六、当你有一段待分析的氨基酸序列（基础操作介绍）
【陪你学·生信】七、在数据库中检索相似的序列
