SAM/BAMからマルチリードを除く方法

マルチリードについて

ゲノム上の複数箇所にマップされたリードをmultiple mapped reads, multi-reads, reads with multiple hits などと呼びます．もしくは，一つのリードに対して複数のアラインメントが出力されるもの，とも言えます．

このマルチリードは，しばしば解析の前に除かれます．では，どうやって除いたらいいでしょうか．

MAPQを使う

MAPQ（SAM/BAMの５列目）が適切にセットされているなら，これを使うことができます．適切というのは，ユニークなアラインメントでMAPQが50，アラインメントが２つだとMAPQが3，アラインメントが3以上だとMAPQが0に設定されている状態です．

このとき，samtools view -qを使って，以下のように書けます．

-q INT Skip alignments with MAPQ smaller than INT [0]

samtools view -q 4 hoge.bam > hoge.uniq.bam

NHを使う

NHタグは，SAM/BAMファイルに，あるリードに対するアラインメントがいくつ入力されたかを NH:i:[0-9]+という形式で記録しています．

NH i Number of reported alignments that contains the query in the current record

（The SAM Format Specification (PDF)）

そのため，スクリプト言語で（私はNHの分布が知りたかったので，こちらの方法を使ってました）．以下は，Perlで書いた例です．

open(I, "hoge.bam");
open(O, "> hoge.uniq.bam");

while(<I>){
     if(/^@/){ # header
          print O $_;
          next;
     }
     if(/NH:i:([0-9]+)/){  # NH: Number of reported alignments that contains the query in the current record
          if($1 == 1){
               print O $_;
          }
     }
}
close(I);
close(O);

FLAGは使えなさそう

掲示板などを読んでいるとFLAG（SAM/BAMの２列目）のビットを使って，すなわちsamtools view -F ***を使う方法が提案されることがしばしばあります．これはおそらく 0x100(secondary alignment)のことを指しているのだと思うのですが，以下の例で示す通り，マルチヒットのリードを除くには不適切です．

=====

read1      339     chr1    181938  3       2M107N99M       =       181891  -255    ACCGGAGGATGTTGGGATGGCAGAGGCGGTTCATGAGCTGCACCTCCCGCAGCATGTTGGCCTTGTTGCTCGCCAGCGTGTTCATCTTCAGCGCCATCACC   BDDDCACCC<<3DDDDCDDDDDD@@9BBCCDCC@49CCCCA<DB?CDFFFFHHHHJJJJJIIJIJJJJJJJGJJJJIJJJJJJIJJJJHHHHHFFEDD?@B   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:101        YT:Z:UU XS:A:-  NH:i:2  CC:Z:=  CP:i:182045     HI:i:0:0
read1      83      chr1    182045  3       101M    =       181891  -255    ACCGGAGGATGTTGGGATGGCAGAGGCGGTTCATGAGCTGCACCTCCCGCAGCATGTTGGCCTTGTTGCTCGCCAGCGTGTTCATCTTCAGCGCCATCACC   BDDDCACCC<<3DDDDCDDDDDD@@9BBCCDCC@49CCCCA<DB?CDFFFFHHHHJJJJJIIJIJJJJJJJGJJJJIJJJJJJIJJJJHHHHHFFEDD?@B   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:101        YT:Z:UU NH:i:2  HI:i:11

これは，ペアエンドのリードをTopHatでアラインメントして得られた，あるリード read1に対する2つのアラインメントです．どちらもNH:i:2であることから，各リードに対して，ファイル中に2つのアラインメントが出力されていることがわかります．これは，MAPQ（５列目）が3になっていることからも分かります．一方で，FLAG（２列目）が，read1では339（0x1, 0x2, 0x10, 0x40, 0x100）で，read2では83（0x1, 0x2, 0x10, 0x20）となっています．２つのアラインメントにEdit distanceなどで違いがないにも関わらず，片方だけに0x100が付いています．

つまり，マルチヒットのリードを除くという目的には不適切だと言えます．

参考

Remove reads which are not uniquely mapped - SEQanswers

（参考）FLAGのビットの意味

The SAM Format Specification (PDF)

Bit	Description
0x1	template having multiple segments in sequencing
0x2	each segment properly aligned according to the aligner
0x4	segment unmapped
0x8	next segment in the template unmapped
0x10	SEQ being reverse complemented
0x20	SEQ of the next segment in the template being reversed
0x40	the first segment in the template
0x80	the last segment in the template
0x100	secondary alignment
0x200	not passing quality controls
0x400	PCR or optical duplicate

Written on July 19, 2013