Question

基因组序列分析中，序列一致性指的是什么？

Answer 1

基因组序列分析中，序列一致性指的是什么？

何谓序列一致性（Sequence identity）？我想大家或多或少应该知道，它指的是两个不同序列之间的相似程度，这是对序列比对结果的一个描述，也常常被用于基因组之间的比较。

测序read比对于参考基因组的时候，这个一致性也可以用来反映测序错误率的大小。不过，不知道大家是否注意到这样的情况，那就是当我们说“两个物种之间的序列差异是多少多少”或者说“测序错误率是什么”时，我们实际上是默认大家都清楚地知道这个“序列一致性”是怎么计算的。而要达到这个“默认”，就意味着这个计算是唯一的。但事实上，序列一致性的计算方法是有多种的，并非唯一，所以这就有必要引起我们的注意了。

李恒博士在他最新开发的长序列比对软件——minimap2时就谈到过这个问题。这篇文章我们借用李恒博士举过的例子来对这个问题进行具体说明。

我们用下面这两条序列为例子，其中第一条是参考序列（Ref），另一个是要去比对的Query序列（Qry），我们来计算一下它与参考序列之间的一致性情况。

Ref: CCAGTGTGGCCGATACCCCAGGTTGGCACGCATCGTTGCCTTGGTAAGC

Qry: CCAGTGTGGCCGATGCCCGTGCTACGCATCGTTGCCTTGGTAAGC

当然要计算一致性就需要先进行序列比对，但既然涉及到序列比对，那么就需要先规定一个打分矩阵，否则我们就不知道该如何比了。按照最经典的做法，比对的打分矩阵一般是这样定义的：

匹配得1分

错配得-2分

开gap是-2分

gap的延伸是-1分

按照这样的打分矩阵，我们可以得到这两个序列的最佳比对情况将如下：

结果有了之后，仔细数一下，我们就可以知道，在上面这个比对例子中，完全匹配的比对碱基是43个、错配1个、同时还有5个序列删除和1个序列插入。如果我们用BAM中的CIGAR字符串来表示这个比对，那么是：18M3D2M2D2M1I22M。

根据这个比对结果，我们目前有三个不同的方法来计算序列的一致性比例。如下：

1. 不算gap

顾名思义，这个方法在计算的时候是不把gap计算在内的，只考虑匹配和错配这两种情况。所以，这时序列一致性的计算公式是这样的：

序列一致性 = 完全匹配数/(完全匹配数+碱基错配数)

在这个例子中，按照这个方法来计算的话，其序列一致性就是：43/(43+1) = 0.977，也就是这个两个序列之间有97.7%的一致性。

这个计算方式也是我们以前比较人类基因组和黑猩猩基因组一致性程度时，得到98.77%这个一致性值所采用的方法，在那篇黑猩猩基因组文章中，原文的描述是这样的：“Single-nucleotide substitutions occur at a mean rate of 1.23% between copies of the human and chimpanzee genome”。

2. BLAST的序列一致性计算方法

BLAST有它自己的一套方法，而且很简单。它在计算序列一致性的时候，直接只计算在总比对长度中，匹配上的碱基数目所占到的比例。需要特别说一下的是，这个总比对长度除了指匹配和错配这两个情况的碱基数目之外，还包括：序列删除和序列插入的长度。比如，在上面的例子里面，虽然Ref序列的长度是49，但是总比对长度却是50（就是要多加1个插入碱基的长度）。所以，按照这个方法，序列一致性的结果就是：43/50=0.86(也就是86%)了。

BLAST的计算方法应该是我们用得最多的一个。不过在用它来对序列比对结果进行过滤的时候就要特别注意了，尤其是当遇到ALU这一类转座子元件的时候更是如此。

因为由于生物本身的机制，这类转座元件在基因组上发生序列插入或者删除时，往往是一次性整体发生的，不会分成多个独立事件分开进行。

举个例子，假设我们有一个1000bp长的序列，中间有一个300bp的ALU序列插入，那么完成比对之后，按照BLAST一致性的计算方式，我们就将只得到70%的序列一致性，但这是不对的，因为，这个300bp的ALU序列实际上是一次性发生的，而不是分成了300次。这时如果仅仅只是按照序列一致性的高低作为条件来过滤比对结果，那么很可能会挑到错误的比对。

3. gap压缩

这个计算方法就比较特别了，其实它的目的是为了解决上面这个BLAST计算一致性时所出现的问题而提出的。这也是李恒在minimap2算法中所使用的方法，它的值记录在minimap2比对结果文件的ed:f标签中，可以直接调出查看。

它的定义和BLAST类似，分母也是要求比对的总长度，但这个“总长度”并不是像BLAST那样的硬长度！而是要把那些连续的gap区域压缩起来，当做只有1个字符来计算，比如下图这样：

左边是原始总长度，右边就是计算的时候把gap压缩成一个字符来看待的长度。不过要注意的是右边这个形式，并不是真实的比对结果，而只是在计算时，从逻辑上将它处理成这样而已。实际上，就是把每一个的gap比对，不管它延伸出来的区域有多长，都只看做是一次整体事件而不是多个——实际上这应该也更符合实际的生物过程，这就可以避免在BLAST中可能碰到的类似情况。

所以，如果按照这个方法来计算，我们上面那个例子的序列一致性就变成了 43/(43+1+2+1) = 91.5%。看这个分母43+1+2+1，后面1+2+1分别代表的是：1个错配，2个序列删除和1个序列插入。

比对所用的打分矩阵也会对序列一致性的计算结果产生影响

比对结果的好坏需要通过打分矩阵来计算和衡量，既然如此，打分矩阵的变换往往将导致比对结果的变化，而比对结果变了，那么序列一致性的计算结果也必然会跟着改变。

比如，我们把上述例子一开始所用的打分矩阵略作修改，将开gap的分数改为 -4 分，其它不变。做完这个修改之后，我们会发现，这一次最佳的比对结果就变成了这样：

按照这个新的比对结果，这两个序列之间的一致性情况如何呢？

按照BLAST的计算方式是：42/50=83.7%，按照gap压缩的方法则是，41/(41+4+1)=89.1%，可以看到这个结果就和原先的都不同了。

小结

举这些例子不是为了在这篇文章里对这些计算方法比个优劣，定个高低，而是想说，有时候对看似简单的问题，也值得我们多留个心眼，不应该只是做一个“规矩”的工具党。当你碰到这类有关序列一致性或者讨论测序错误率的时候，应该先搞清楚前提条件是什么，是用什么方法计算的，打分矩阵是什么样的，于你关心的是否相同等，不要想当然，这个序列一致性并没有一个唯一的计算方式（话虽如此，我觉得gap压缩的方法应该是更优的一种）。如果你不怎么关注它的计算细节和条件，那么在做研究和分析的过程中踩到不该踩的坑也只是迟早的事了。