こんにちは。h--ishiです。
さて、みなさん！
みなさんは、bamファイルのカバレッジ状況をひとまず簡便に確認したいといったことはありませんでしょうか。
今回は、samtools coverage についてご紹介いたします。
なにかと多機能なsamtools ですが、実はsamtools 1.10 から便利な機能が実装されていました。
その機能とは、
bamファイルのマッピングのstatsを計算してくれる機能です。
具体的には、マッピング結果のカバレッジ領域長、カバレッジ率、平均デプスなどの値がわかります。
実際に使ってみた方がわかりやすいかとおもいますので、早速使ってみましょう！
基本的な使い方は以下のように、statsを知りたいbamファイルを指定するだけです。
samtools coverage input.bam -o output.tsv
出力ファイルは以下のようなタブ区切りテキストになります。
less output.tsv #rname startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq contig1 1 631 260 631 100 57.6282 39.7 17.6 contig2 1 3590 782 1623 45.2089 29.8148 40 32.2 ...
上記結果について、みていきます。
contig1 についてみてみましょう。
260本リードがマップしていて（numreads列参照）、
それらのリードはリファレンス配列の631塩基にマッピングしていること（covbases 列参照）がわかります。
そのため、カバレッジが100%になってます。
contig2 のカバレッジは、45%程度でした。
各カラムの説明です。
- rname: リファレンスの染色体名（or コンティグ名）
- startpos: 集計開始位置
- endpos: 集計終了位置（対象の配列長）
- numreads: 指定した領域にアライメントされたリード数
- covbases: リードデプスが1以上の塩基数
- coverage: 集計対象の全塩基数に対するリードデプスが1以上の塩基数の割合
- meandepth: 平均デプス
- meanbaseq: マップされたリードの平均ベースクオリティ
- meanmapq: 平均マッピングクオリティ
集計に使用するリードをリード長・ベースクオリティ・マッピングクオリティなどでフィルタリングすることもできるようです。
より詳細な使用方法は以下の公式マニュアルをご参照ください。
公式マニュアルはこちら
