はじめに
GATKという重厚長大型のツールがあります。
このツールはJavaで書かれた非常に広範囲に及ぶツールで、主にBamファイルから変異を呼び出すときに使われます。
このツールの中で比較的よく使われる機能に、MarkDuplicates があります。Picardというツールの機能です。
MarkDuplicates
これは、bamファイルから重複を削除するツールだと言われていますが、「重複をマーク付ける」という曖昧なキーワードで説明が終わっていることが多かったので、実際には何が起きているのか見ることにしました。
忙しい人のために結論だけ説明すると、MarkDuplicatesは、Flagを変更していました。
よく「MarkDuplicateは重複をマークする」「タグを付ける」という説明があります。確かにタグをつけるオプションはありますが、これだとMarkDuplicatesが何をしているのか不明瞭で説明としてはあまりよろしくないと思います。
テストデータの準備
まずは、テストデータを準備します。ここでは、wgsimを使って簡単にでっち上げたsamファイルを使います。
@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:1000
@SQ SN:chr2 LN:900
read1 73 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222
read1 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222
はい。最小限のヘッダーに、レコードが２個しかないSAMファイルでございます。
えっと、思った人もいるかも知れませんが、上記は正常なSAMファイルであり、samtoolsに投げても認識してくれます。
$ samtools flagstat moo.sam
2 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
1 + 0 mapped (50.00% : N/A)
2 + 0 paired in sequencing
1 + 0 read1
1 + 0 read2
0 + 0 properly paired (0.00% : N/A)
0 + 0 with itself and mate mapped
1 + 0 singletons (50.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
１番目のレコードと、２番目のレコードは、qname が同じです。それぞれのflagは、73, 133であり1番目のレコードはマッピングされているが、2番目のレコードはうまくマッピングされてなかったことがわかります。
samtools flags 73 # PAIRED,MUNMAP,READ1
samtools flags 133 # PAIRED,UNMAP,READ2
このSAMファイルの２つしかないアラインメント・レコードを増やします。1番目のレコードを3本に、2番目のレコードを4本に増やしてみましょう。そのままだとqnameが同じでおかしなことになるので、名前を少し変えます。
@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:1000
@SQ SN:chr2 LN:900
read10 73 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222
read11 73 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222
read12 73 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222
read20 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222
read21 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222
read22 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222
read23 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222
最後のリード(read23)に関しては、ペアがないですが、別にいいとしましょう。単に重複させたので、sortされた状態は維持されています。
@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:1000
@SQ SN:chr2 LN:900
read10 73 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222
read11 73 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222
read12 73 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222
read10 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222
read11 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222
read12 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222
read23 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222
このファイルに対してMarkDuplicatesしてみましょう。
gatk MarkDuplicates -I moo2.sam -M metrics.txt -O out.sam
すると次のような結果が出ます。
@HD VN:1.6 SO:coordinate
@SQ SN:chr1 LN:1000
@SQ SN:chr2 LN:900
@PG ID:MarkDuplicates VN:Version:4.4.0.0 CL:MarkDuplicates --INPUT moo2.sam --OUTPUT out.sam --METRICS_FILE metrics.txt --MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP 50000 --MAX_FILE_HANDLES_FOR_READ_ENDS_MAP 8000 --SORTING_COLLECTION_SIZE_RATIO 0.25 --TAG_DUPLICATE_SET_MEMBERS false --REMOVE_SEQUENCING_DUPLICATES false --TAGGING_POLICY DontTag --CLEAR_DT true --DUPLEX_UMI false --FLOW_MODE false --FLOW_QUALITY_SUM_STRATEGY false --USE_END_IN_UNPAIRED_READS false --USE_UNPAIRED_CLIPPED_END false --UNPAIRED_END_UNCERTAINTY 0 --FLOW_SKIP_FIRST_N_FLOWS 0 --FLOW_Q_IS_KNOWN_END false --FLOW_EFFECTIVE_QUALITY_THRESHOLD 15 --ADD_PG_TAG_TO_READS true --REMOVE_DUPLICATES false --ASSUME_SORTED false --DUPLICATE_SCORING_STRATEGY SUM_OF_BASE_QUALITIES --PROGRAM_RECORD_ID MarkDuplicates --PROGRAM_GROUP_NAME MarkDuplicates --READ_NAME_REGEX <optimized capture of last three ':' separated fields as numeric values> --OPTICAL_DUPLICATE_PIXEL_DISTANCE 100 --MAX_OPTICAL_DUPLICATE_SET_SIZE 300000 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false PN:MarkDuplicates
read10 73 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222 PG:Z:MarkDuplicates
read11 1097 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222 PG:Z:MarkDuplicates
read12 1097 chr1 179 60 40M = 179 0 TTGGATACCATTCCCCACAAAGGTACATAATGATGTCCTC 2222222222222222222222222222222222222222 PG:Z:MarkDuplicates
read10 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222 PG:Z:MarkDuplicates
read11 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222 PG:Z:MarkDuplicates
read12 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222 PG:Z:MarkDuplicates
read23 133 chr1 179 0 * = 179 0 TGTGCAGGATCCCAGTCTGTTGTTACAGAAAACTATACAC 2222222222222222222222222222222222222222 PG:Z:MarkDuplicates
どこが変わったか見てみましょう。
はい。meldで差分を表示すると一目瞭然ですね。(vimdiff, diff -y, delta -s)
①ヘッダー部分に@PGタグで長大なコマンド履歴がつきました。何もオプションを指定していないのですが、使用されたデフォルトのオプションがここに記録されるようです。
②PG:Z:MarkDuplicates というストリング型のタグがAuxのPGに追加されました。
このタグはPG（Program）のタグであり、MarkDuplicateを実行したという意味ですべてのレコードに追加されます。この例ではすべてのリードが（Mapされてないだけで）重複しているのでわかりにくいですが、重複していないリードに対してもこのタグは付加されます！
③最後に、Read11, Read12の2本のリードのフラグに 1024（DUP）が追加されました。73 + 1024 = 1097 ビット演算ですね。Read10はオリジナルのリードと見なされてDUPは付加されませんでした。興味深いことにDUPが追加されたのはマップされたRead1x系のみであり、マップされてないRead2x系のflagは変更されませんでした。マッピングされてないリードの重複は無視されるようです。
「重複をマークする」の本当の意味
MarkDuplicates で重複しているリードにタグをつけるという説明はあまり適切ではありません。タグという単語はAuxを想起させますが、MarkDuplicateは Flag を立てています。コンピュータの気持ちになれば、たしかに010101...で表現されるバイナリの 1024のビットのところにマークをつけていることになりますが、人間としては「MarkDuplicatesはFlagを変更する」と理解した方がよいと思います。
--REMOVE_DUPLICATES をつけてみる
gatk MarkDuplicates -I moo2.sam -M metrics.txt -O out.sam --REMOVE_DUPLICATES
すると、単に重複したリードが消されてしまいます。何本の重複リードがあったかといった情報はオプションをつけないと記録に残らないようです。
どうすればいいんだ…
バイオインフォのファイルフォーマットとツールは複雑です。MarkDuplicatesひとつとっても、多数のオプションがあります。またGATKはバージョンによって挙動に違いがあることでも知られています。このような複雑性に対してどのように向き合っていけばよいのでしょうか？
実際のところ、こんな感じで毎回実験しながら確かめていくしかないのではないかと思います。そもそもバイオインフォの難しさの大半は、データやファイルの重厚長大さから来ているため、なるべく小さなダミーデータを用意して、しつこく実験していく。そして、そういった実験を、少ないエネルギーで簡単にできる環境を整えていくことが大事ではないかなと思いました。
この記事は以上です。