继上周后，介绍完Heng Li大神的bioawk之后，继续在他的博客闲逛，读到一篇简单也很有趣的文章。下面我会总结一些该文章的要点，和大家分享分享。希望大家可以从该文章中学到一些知识。
概念介绍
首先，本文的关键词，什么是Sequence identity？简单来说Sequence identity就是两个序列之前的相似度，简单易懂。通常来说序列之间的相似度就是序列错误的反面。
引入
目前，有好几种不同的方式去计算序列的相似度。下面会和大家一步一步来探索：
首先给大家一个例子，以下两段序列的相似度是多少？
CCAGTGTGGCCGATACCCCAGGTTGGCACGCATCGTTGCCTTGGTAAGC
CCAGTGTGGCCGATGCCCGTGCTACGCATCGTTGCCTTGGTAAGC
要是按照最经典的方式去计算相似度，第一步当然是做比对，然后如果我们定义：能比对上的序列(match=1)给一分
，比对不上扣两分（mismatch=-2），存在gap的(gapOpen=-2）扣两分还有存在gap延长的(gapExt=-1)：
Ref+: 1 CCAGTGTGGCCGATaCCCcagGTtgGC-ACGCATCGTTGCCTTGGTAAGC 49
|||||||||||||| ||| || || ||||||||||||||||||||||
Qry+: 1 CCAGTGTGGCCGATgCCC---GT--GCtACGCATCGTTGCCTTGGTAAGC 45
对你Ref序列和比对序列，这里我们可以看到有43个是比对上的碱基，1个是比对错误的，然后5个缺失还有一个插入。这个比对结果对应的CIGAR是18M3D2M2D2M1I22M。这里这两个序列的相似度是45/49=91.8%
如果你不知道什么是CIGAR，这里简单介绍一下：CIGAR就是一行数字和字母，用来浓缩介绍你比对的结果。例如上面的18M3D2M2D2M1I22M，就一一对应了每一段比对的结果，前18个碱基是比对上ref序列的，所以是18M，接着是有三个缺失，对应3D，以此类推。你可以通过CIGAR快速的了解到比对的结果。
好继续回归正题,一般来说序列的相似度可以分为以下三代类：
- Gap-excluded identity
- BLAST identity
- Gap-compressed identity
Gap-excluded identity
Gap-excluded identity顾名思义就是在计算序列相似度时候不考虑gap的存在。具体让我们看看上面的例子，如果我们不考虑gap上面的碱基的话，上面的相似度就会改变为43/(43+1)=97.7%，比上面提升了不少。在一篇猩猩的基因组文章中有这样的一句话
Single-nucleotide substitutions occur at a mean rate of 1.23% between copies of the human and chimpanzee genome
这里的SNP substitution rate就是基于gap-excludeed的相似度得出的结果（差别会小很多）。
BLAST identity
BLAST identity是我们平常使用得比较多的一个。BLAST identity是根据你比对到的碱基除去比对所涉及到的columns数目，换句话来说就是比对涉及到所有的碱基数目。
用上面的例子来讲解一下，其BLAST identity会是43/50=86%。我们可以通过sam file的信息进行BLAST identity的计算，在sam file中，我们可以通过叠加CIGAR中对应的M/I/D的数目来计算。然后比对到的碱基数目等于column的总长度，减去NMtag（比对不上的碱基位置的标记）。这一步可以通过单行perl command来实现：
perl -ane 'if(/NM:i:(\d+)/){$n=$1;$l=0;$l+=$1 while/(\d+)[MID]/g;print(($l-$n)/$l,"\n")}'
虽然BLAST identity可能是最常见最常用的一种定义，但是在不同情况我们一定要小心对其的使用。下面例举一些我们需要注意的情况：
假设我们正在比对一个1000bp的序列，这个序列中间含有一个~300bp的ALU插入。BLAST identity的值大概会在70%左右，这段序列很可能会被过滤掉。但是问题来了，在进化中，一个ALT的插入往往是在一个进化事件中产生的，所以这300bp不应该被看作300个碱基独立的缺失。用BLAST indentity在这种情况下就会不太合理，有些不该被过滤的序列会被过滤掉。
Gap-compressed identity
那上述的问题该怎么解决，这里大神提出，这种Gap-compressed identity的方法会更加合理。这种方法是将连续的gap看成一个base的不同。
现在比对就变成这样了：
Ref+: CCAGTGTGGCCGATaCCCcGTtGCTACGCATC-TTGCCTTGGTAAGC
|||||||||||||| ||| || |||||||||| ||||||||||||||
Qry+: CCAGTGTGGCCGATgCCC-GT-GCTACGCATCgTTGCCTTGGTAAGC
其相似度就变成43/(50-2-1)=91.5%。将连续的gap变成单一的base，总长度分别减少了2和1。然后大神将这种计算方式也运用到了他最新的工具minimap2中。
Effect of scoring
大神在最后还说了关于分数的问题，他提到就算你只有一个定义的方法，但是不同的扣分，加分制度会非常影响你比对相似的大小。
比如，将gapOpen=-4，出现gapOpen减4分。
Ref+: 1 CCAGTGTGGCCGATaCCCcaggtTGgcACGCATCGTTGCCTTGGTAAGC 49
|||||||||||||| ||| || ||||||||||||||||||||||
Qry+: 1 CCAGTGTGGCCGATgCCC----gTGctACGCATCGTTGCCTTGGTAAGC 45
这个比对的BLAST identity是83.7% 其gap-compressed 的identity会是89.1%，对比上面相似度不同程度的降低了。
小结
当你估算序列间的相似度的时候，不同的定义计算方式还有不同的比对打分都会影响你的相似度的结果。所以在比较序列间相似度时，要选取切当的方法，并保持方法和分数一致。
原文链接：
http://lh3.github.io/2018/11/25/on-the-definition-of-sequence-identity
