BAMファイルのMDタグについて
概要
BAMファイルには、リファレンス配列に関連する情報を含む「MD」というタグが存在することがあります。MDタグについて調べた結果、英語のブログ「MD Tags in BAM Files - vince buffalo」で詳しい解説がされていました。
Vince Buffalo さん敬意を払いつつこの記事の内容を要約してQiita記事にまとめます。
MDタグの役割
MDタグは、リファレンスを参照しなくてもSNPやインデルを特定できるように設計されています。特定のアラインメントのリファレンス情報を保持する。
SNPの記録
例えば「10A3T0T10」というMDタグがあった場合、以下の3つのSNPが含まれる：
- 11番目のリファレンス塩基はAです。
- 15番目のリファレンス塩基はTです。
- 16番目のリファレンス塩基もTです。（0は隣接した塩基をあらわします）
欠損
MDタグは欠損情報を保持します。例：「85^A16」では、85個のマッチ後にリファレンスがAである位置に欠損があり、その後16個のマッチが続きます。
挿入の扱い
挿入はMDタグに記録されません。例えばリード長が101で、「78M1I22M」というCIGAR文字列がある場合、CIGARは1つの塩基の挿入を表現していますがMDタグは単に「100」です。これはリファレンスの情報は失われていないからです。
複雑なケース
挿入がある場合、例えば「リードの長さ：101、CIGAR: 9M1I91M、MD: 48T42G8」のような場合、MDタグは挿入を無視してマッチ部分のみを表し、リード内で挿入された塩基はMDタグには反映されません。
49番目の塩基（リファレンスはT）は実際にはリードの50番目の塩基でです。ミスマッチしたTの上流に1塩基の挿入があります。もう1つのミスマッチ（リファレンスはG）についても同様で、MDタグによると48 + 1 + 42 + 1 = 92塩基目ですが、実際には93塩基目です。
考え方とツール
MDタグの解析のアプローチとしては、MDタグとCIGAR文字列を利用してリファレンス文字列を再構築し、それをもとにリードとの比較を行います。
また、samtools calmd コマンドを -e オプション付きで使用すると、塩基が = にマスキングされ検証が容易になります。
※こんな感じの表示が得られる
=A=======G==A======================C============C=====================
IGVによるリードの観察も有用です。
この記事は以上です。
なお、Vince Buffaloさんは、名著 Bioinformatics Data Skills の作者です。買おう。