生物信息就该这么学（5）：其实就是比对那点事儿
今天继续我们的生物信息就该这么学内容，其实生物信息分析每天做的工作就是各种比对。序列比对可以说是整个生物信息的核心，因为你会发现几乎每个生物信息分析过程都需要用到序列比对。判断两个基因或两段基因组片段是否相似是序列分析的基本工作。从序列数据库搜索，序列拼接到基因蛋白质功能注释，以及进化树构建等，都依赖于分子序列相似性的比较，也就是序列比对。测序得到一段序列，判断它是否属于新冠病毒，需要与已知病毒序列进行比对。你可以会反问，基因组拼接，变异检测也用到序列比对了吗，答案是肯定的。
同源与相似
序列比对的核心作用就是判断是否同源。所谓同源(homology)，是整个生物信息分析中最为重要的一个概念，只有同源的比对分析才是有意义的。同源也就是指来自于同一个祖先，两个物种从同一个祖先分化后，与不同的环境发生相互作用，其相应的DNA序列将各自发生一些替换或者插入缺失突变，也就是说序列不在精确相同。比如同一个大肠杆菌的祖先，经过不同的时间和空间的差异累积，最终就分化成不同类型的大肠杆菌。不同样品之间是同源的关系。
与同源概念想对应的是序列相似，相似性(similarity)和同源性(homology)是两个完全不同的概念。相似性仅仅是指字符串的相似 ，并不具有不具有生物学意义 ，因为DNA序列一共就有ATCG四种碱基，由于组合造成两段片段字符串组合比较接近。同源序列一般是相似的，但是相似的序列不一定同源。那么该如何判断序列是相似还是不相似，相似的序列是否满足同源关系呢。这些都需要序列比对来判断，并且使用一些方法和标准来进行评价。通常的解决方法是将两条序列或者两条序列中的部分进行序列比对，然后基于得到的比对结果判断相似性是由于序列之间具有进化渊源，还是纯粹的随机巧合。
比对算法
很多一看到比对算法就感觉很高深，不敢去了解，其实并没有那么复杂了。对于序列比对需要首先明确以下几个问题。
首先，比对的对象是什么，是DNA序列还是氨基酸序列？
第二，如何判断相似性？
第三，使用哪种算法找到全局或者局部的最优匹配？
最后，如何判断序列同源性高低？
两条序列比对碱基之间可能有以下几种情况。两个碱基类型完全一致，没有发生变化，另外就是碱基发生了变化，这种变化包括替换、插入和删除。插入和删除也被称为空位，我们平时在比对过程中的错配其实就是替换，gap就是插入或者删除。因为突变是随机的，但是选择是具有偏向性的，这就使得某些突变发生的可能性远大于其他类型。
我们根据比对过程中不同的情况给定分值，例如完全匹配加1分，错配减一分，空位减一分，那么最终两条序列比对我们就会计算得到一个分值。这个分值越高就说明完全匹配的碱基越多，比对的情况越好。这就是打分矩阵。比对后根据打分矩阵的值就能判断出比对的情况。打分矩阵有很多种，那么最简单的就是这种等价矩阵表。但是实际上我们知道这样明显会有很多问题，因为DNA上四种碱基替换比率是不一样的，由于化学结构上的相似性，嘌呤与嘌呤之间，嘧啶与嘧啶之间更容易发生替换，这种称为转换，而嘌呤与嘧啶之间发生替换稍微难一些，这种称为颠换。所以，我们将打分矩阵修改一下。完全匹配加2分，转换扣一分，颠换扣两分，这个打分矩阵要比前面一个更加精确一些。
不过这两种矩阵我们都还没有考虑到空位罚分，对于插入或缺失突变，目前为止，还没有适当的统计模型。通常用的是与替换矩阵相对应的经验分值，一般空位罚分有两个参数值，即起始空位罚分和空位延伸罚分。对于一个长度为k的连续空位，罚分可以表示为score=a+b*k，其中a为起始空位罚分，b为空位延伸罚分。
对于序列比对罚分矩阵有很多，但是没有一种矩阵能够适用于所有的情况，都会存在一定的误差，对于DNA序列比对，主要有mat50和mat70两种矩阵，对于氨基酸序列比对，主要采用PAM和BLOSUM两种矩阵。blast比对中默认使用的就是 BLOSUM62打分矩阵。其中62表示用来构建该矩阵的匹配数据集中精确匹配位点要占62%。
此外，还有转移矩阵，等价矩阵，遗传密码矩阵，疏水性矩阵等。
在确定了打分矩阵之后，接下来就是设计有效的算法，寻找最优的序列比对。如果不允许插入空位，序列比对的算法是非常简单的。对于全局比对，只有一种比对方法，然而由于序列突变过程的特点，不考虑空位插入的模型是不符合实际情况的，但是如果要允许存在空位，就会一下子复杂很多。所以开发出了很多比对的算法。前面我们介绍过短序列比对使用BWT算法。全局比对使用Needleman算法，而局部比对一般使用Smith-Waterman算法，另外，还有启发式联配算法等。对于算法感兴趣的可以自行搜索研究。这里面我们不做过多介绍。一般的比对过程是先用一小段序列，作为种子序列，进行精确匹配，这样一个种子序列可能与基因组上多个区域匹配上。接下来就进行逐渐延伸，延伸的过程中由于错配和空位的存在，比对分值就会降低，很多位置匹配的就比较差，最终只有一个可以延伸的很长，那么这个就是最优的匹配。
最后我们要解决的就是如何来判断序列同源性的高低，也就是如何评价序列匹配的的显著性。也就是哪种匹配是具有生物学意义的匹配，哪种是纯粹的随机巧合。一般情况下同源性越高，比对的长度越长，错配少，相对应的indentity值就越高，比对的分值也越高。在blast比对中采用E值来衡量显著性。E值综合考虑了真实匹配与随机匹配的相对概率，序列的长度，数据库大小和序列组分的偏向性等数据而计算出来的。例如，我们常用的1e-5。表示这样一组最优匹配是随机巧合而不是具有生物学意义匹配的概率为十万分之一，比对的同源性很高。
全局比对与局部比对
全局比对是用来衡量两条序列整体的相似性，满足整体相似性最大化。若两条序列长度不同，则必须插入一些空位使所有位点都能对应起来。而局部比对则不同，两条亲缘关系较远的DNA或氨基酸可能只在一些片段上相似，这就需要找到这些相似性的片段，和其相应的匹配方式。通常这样的分析就需要进行局部比对，而不是全局比对。
全局比对与局部比对有什么不同呢。全局序列比对尝试找到两个完整的序列之间的最佳比对。而局部序列比对不必对两个完整的序列进行比对；可以在每个序列中使用某些部分来获得最大得分。两种比对采取不同的比对算法和策略，因此，同样的一段序列，采用全局比对和局部比对不同的比对方法结果也会有很大的不同。
例如我们现在有两条序列S1和S2，如果采用全局比对，会得到这种比对效果，而采用局部比对，序列中间的GCG满足了最优比对。大家可以理解为，全局比对需要从全局出发，是需要全局达到最佳效果，而局部比对则不需要考虑全局，只要局部达到最佳效果即可。
全局比对主要用来比较比较两个基因组之间的同源性，绘制共线性图等，另外，全局比对也常常用于基因组结构变异的检测。因为，局部比对的话，遇到大的空位往往就断开了，例如上面的例子，采用局部比对的算法中，只追求局部的最优比对，而不会考虑整体的空位等。所以，基因组的大片段的插入或者缺失检测，可以使用全局比对软件。而局部比对软件主要搜索同源序列，例如判断那两个基因是否同源，寻找一段序列的同源序列等，就可以使用局部比对。
多序列比对
多序列比对（Multiple sequence alignment；MSA），是对三个以上的生物学序列，如蛋白质序列、DNA序列或RNA序列所作的序列比对。一般来说，是输入一组假定拥有演化关系的序列。从多序列比对的结果可推导出序列的同源性，而种系发生关系也可引导出这些序列共同的演化始祖。如图可视化展示的多序列比对，例如单碱基的错配，或是如删除突变与插入突变。都可以显示出来。多序列比对的结果常常用于系统发育分析。
短序列比对
短序列比对是伴随着高通量测序数据应运而生的。所谓短序列比对就是将这些测序的reads重新定位到基因组上，这个过程也叫做回贴或者叫做mapping。
我们可以将测序的reads重新定位到自身的基因组，也可以定位到参考基因组上，这个根据不同的研究目的而定。同样，不仅可以定位到基因组序列，也可以定位到基因集序列。或者rRNA序列等。短序列比对，至少要需要一个fasta格式的目标序列，然后是fastq格式的测序reads，可以单端也可以是双末端。注意短序列比对只能DNA与DNA进行比对，不能在氨基酸水平继续比对。
短序列比对主要是为了得到这种堆叠效果（pileup），通过堆叠可以计算每个位点的比对细节，包括覆盖度，覆盖比率，reads利用率，每个位点具体比对情况。这样就可以用于计算基因差异表达，丰度差别等。RNAseq，变异检测，宏基因组，组装结果纠错等都要利用短序列比对的结果进行后续分析，可以说短序列比对是整个高通量测序的核心工作。