序列相似度（identity）是兩個序列之間相似度的路徑成本。在測序資料中一般是測序錯誤率的反面。當我們說“兩個物種的差異是。。”或者“測序錯誤是。。”，事實上有不同的相似度的計算方法，這裡我們讨論了集中算法，每次在說相似度時，要搞清楚用的是那種算法。
首先舉個例子說明，下面兩條序列的相似度是多少呢？
CCAGTGTGGCCGATACCCCAGGTTGGCACGCATCGTTGCCTTGGTAAGC
CCAGTGTGGCCGATGCCCGTGCTACGCATCGTTGCCTTGGTAAGC
計算相似度，通常先要進行比對，我們使用match=1,mismatch=-2, gapOpen=-2 and gapExt=-1的參數進行比對，得到如下結果
有43個matching的堿基，1個mismatch，5個deletion，1個insertion，CIGAR是18M3D2M2D2M1I22M.
Gap-excluded identity
這種計算法，我們首先排除比對中gapped列，也就是隻保留match，和mismatch，相似度就=matches / (matches + mismatches)。上面例子的相似度就是43/(43+1)=97.7%.
這種計算方法明顯的問題就是沒有考慮gaps。然而這還是一種經常使用的定義。經常有人說“人和大猩猩基因組隻有百分之幾的差異”，這種說法就是使用的這種計算方法。論文（first chimpanzee genomepaper）中的原句是“Single-nucleotide substitutions occur at a mean rate of 1.23% between copies of the human and chimpanzee genome”，意思是“任何大猩猩基因組的單堿基替換率平均為1.23%”
BLAST identity
BLAST相似度定義為matching的堿基占總的比對列數的比例。上例中，比對的列數一共50列，相似度為43/50=86%。在Sam格式中，比對的列數可以用CIGAR字元串計算（M/I/D的綜述），matching堿基等于比對總列數減去NM tag，可以用下面Perl一行指令計算BLAST特異性
perl -ane 'if(/NM:i:(\d+)/){$n=$1;$l=0;$l+=$1 while/(\d+)[MID]/g;print(($l-$n)/$l,"\n")}'
$n是Mismath/I/D的總和，$l是比對的長度。PAF格式中，第10列初一第11列就是BLAST相似度
BLAST相似度也許是最普遍的定義，但是我們也必須謹慎的用它來過濾比對資料。假設我們用一條帶有~300bp 插入偏大的1000bp序列進行比對。相似度大約70%，我們可能就把他過濾掉了。但是在進化過程中這個插入可能隻是一次插入事件，而不是300次，是以把他當成300bp的差異就會有一些問題。
Gap-compressed identity
對于過濾這也許是比較好的定義：我們把連續的gaps認為是一個差異，把連續的gap壓縮成一個，例如把握們的例子壓縮成如下的比對。
相似度為43/(50-2-1)=91.5%。我已經在幾個任務中使用這種定義。最新版的minimap2在de:f tag裡輸出這種特異性。也可以用下面的perl指令計算
perl -ane 'if(/NM:i:(\d+)/){$n=$1;$m=$g=$o=0;$m+=$1 while/(\d+)M/g;$g+=$1,++$o while/(\d+)[ID]/g;print(1-($n-$g+$o)/($m+$o),"\n")}'
$m是M的綜合，$g是I和D的和，$o是gap open的數字。
罰分的影響
罰分會影響比對的結果和相似度。如果我們使用gapOpen=-4比對上面的兩條序列，會得到如下的比對結果
BLAST 相似度 83.7% ，gap-compressed相似度89.1%。即使我們使用同一定義，不同的罰分也會得到不同的特異性。
結論
不同的定義和不同的罰分都會得到不同的相似度或錯誤率。所有當和别人讨論“測序錯誤”或“相似度”等的時候，先問清楚采用的那種定義，罰分怎麼設定的，這樣才有可比性。如果自己評估可以使用下面的指令
minimap2 -c ref.fa query.fa | perl -ane 'if(/tp:A:P/&&/NM:i:(\d+)/){$n+=$1;$m+=$1 while/(\d+)M/g;$g+=$1,++$o while/(\d+)[ID]/g}END{print(($n-$g+$o)/($m+$o),"\n")}'