在李恒的github博客 On the definition of sequence identity 當(dāng)中看到這樣一個Perl單行代碼：
$ perl -ane 'if(/NM:i:(\d+)/){$n=$1;$l=0;$l+=$1 while/(\d+)[MID]/g;print(($l-$n)/$l,"\n")}'
這行代碼的目的是為了計算SAM文件中每條記錄BLAST identity
BLAST identity是根據(jù)你比對到的堿基除去比對所涉及到的columns數(shù)目，換句話來說就是比對涉及到所有的堿基數(shù)目
例如這樣的雙序列比對：
Ref+: 1 CCAGTGTGGCCGATaCCCcagGTtgGC-ACGCATCGTTGCCTTGGTAAGC 49 |||||||||||||| ||| || || |||||||||||||||||||||| Qry+: 1 CCAGTGTGGCCGATgCCC---GT--GCtACGCATCGTTGCCTTGGTAAGC 45它們的BLAST identity就是
43/50=86%
那么要計算SAM文件中每條reads的BLAST identity进鸠，總長可以通過疊加CIGAR中對應(yīng)的M/I/D的數(shù)目得到腐缤，比對到的堿基數(shù)目等于總長減去NMtag（比對不上的堿基位置的標(biāo)記）
SAM文件中的NMtag
李恒的這行代碼中有一部分一開始沒有讀懂述吸，就是下圖紅框中的那部分：
其實這是一種簡寫方式壁却，正規(guī)完整且更容易讀懂的形式可以寫成下面這樣：
# 這里為了更好看巧婶，添加了適當(dāng)?shù)膿Q行和縮進(jìn)
$ perl -ane \
'if(/NM:i:(\d+)/){
$n=$1;
$l=0;
while(/(\d+)[MID]/g){
$l+=$1;
}
print(($l-$n)/$l,"\n");
}
'
while(/(\d+)[MID]/g)中的正則表達(dá)式/(\d+)[MID]/g缸血，引起了我極大的好奇：它是在正則表達(dá)式后面添加了一個g字符物蝙，即開啟了全局匹配，又由于是在while( )中進(jìn)行的正則匹配蚜迅，等于是開啟了循環(huán)匹配舵匾，即對于CIGAR字符串18M3D22M，正則表達(dá)式/(\d+)[MID]/g谁不，先會匹配上18M坐梯，然后會匹配上3D，最后匹配上22M
很有意思的用法
參考資料：
