yokaのblog

湖で微生物の研究してます

コンティグ・ゲノムへのリードマッピング・カバレッジの計算方法

 メタゲノムやメタトランスクリプトームの解析をやっていると、コンティグやゲノムのカバレッジを決めたい場面が多々出てくる。アセンブル前のリードを、アセンブルしたコンティグやゲノムにアライメントして、マッピングされた(貼り付いた)リードの量(正確にはアライメントされた長さで割った厚み)を評価することで、アセンブルされたゲノムや遺伝子が元のリード中にどれだけ含まれていたかを評価する、というのがカバレッジの考え方だ。
 アセンブルされた各ゲノムや遺伝子の環境中の現存量を把握したり、メタゲノムのビニングに使ったりと、使いどころの多いカバレッジだけど、これをちゃんと計算するのは結構めんどくさい。BWAやBowtie(2)といったマッピングソフトを使って、インデックスを作って、samやbam形式でアウトプットを出して、それをさらにsamtoolsでソートして、そこからさらに別のツール(bedtoolsとか)でカバレッジに変換するというのが一般的なやり方だと思う。加えて大変なのが、中間生成物であるsamやbamのファイルサイズがでかくて、普通に数GBを超えてくる点だ。巨大なファイルが次々と生成されてくるので、読み込みやハンドリングに時間がかかるし、ディスクスペースも圧迫される。カバレッジの情報が欲しいだけなのに、こんな大げさな事をやらなければならならず、解析作業の中でもカバレッジの計算は面倒くさい作業の一つになっていた。
で、

手っ取り早く、生リードとリファレンス配列を投げたらカバレッジを返してくれるソフトは無いものか・・・

とあれこれ探していて、最近見つけたCoverMというツールが素晴らしかったので紹介。
 名前を見てわかるように、CheckMなどで有名なHugenholtzとTyson のラボのM-Toolsの一つとして開発されている。ツール自体の論文はまだpublishされてないようだけど、GitHubのリンクを引用する形でこのツールを使ってすでに論文になっているものもある。
github.com
 dependencyマッピングソフトとして使われているのがminimap2だ。2018年に論文になったばかりだけどすでに広く使われていて、速いだけでなく、ロングリードのマッピングにも強いのが特長だ。BWAやBowtieの上位互換として考えてよいのかなと思い、個人的にも乗り換えたところだけど、minimap2・BWAの作者のブログによるとminimap2がすべての点においてBWAを上回っているわけではなく、まだBWAの出番もあるだろうとのこと。なお素晴らしいことに、CoverMはオプションでBWAをマッピングソフトに使うこともできる(デフォルトではminimap2を使う)。
 インストールは超簡単。dependencyはminimap2とsamtoolsだけで、Linuxならダウンロードしたcovermというバイナリにパスを通すだけで動く。基本的な使い方はとてもシンプルで、マッピング先にしたいコンティグが入ったfastaと、マッピングしたいペアエンドのリード(.fq.gz形式でOK)を指定して、

coverm contig -r [reference_contigs.fasta] -1 [read_R1.fq.gz] -2 [read_R2.fq.gz] > [output_name.tsv]

という形だ(オプションでシングルエンドも指定可)。これでfasta内の各コンティグのmean coverageの計算結果がoutput_name.tsvに出力される。indexをつくる必要もないし、bamやsamでディスクが埋まる心配をすることもない。内部ではもちろんindexを作ったりbamができたりするステップがあるのだけど、最終的なアウトプットは一つのtsvファイルだけだ。「そうそう、ただカバレッジが欲しいだけなんだよ!」という欲望に真っすぐに答えてくれる素晴らしいツールだ。
 --sharded オプションを付けることで、複数のfastaをreferenceにした結果を一つのtsvに出力することもできる。スレッド数を指定する-tと、マッピング時の相同性の閾値を指定する--min-read-percent-identityと組み合わせると、以下のような感じ。

coverm contig -r [reference1_contigs.fasta] [reference2_contigs.fasta] [reference3_contigs.fasta] ... --sharded -t [nCPUs] --min-read-percent-identity 0.99 -1 [read_R1.fq.gz] -2 [read_R2.fq.gz] > [output_name.tsv]

 デフォルトではコンティグ名とmean coverageだけの2列のtsvが出てくるけど、その他にも色々な情報を一緒に出力できる。例えばカバレッジだけでなく、マッピングされたリードの数や、それを長さで割った値も出してしてくれるので、メタトランスクリプトーム解析でも使いやすい(ちなみにカバレッジの計算ではデフォルトでリファレンスの両端の75 bpを無視する設定になっているので、短い配列にマッピングするときは注意。--contig-end-exclusionで設定可)。さらにRPKMまで計算してくれる(ただし、マッピングする相手をきちんと揃える必要があるのと、標準化された相対値なのでサンプル間で値を比較するときは注意が必要)。詳細は本家サイトのCalculation methodsを参照。僕の好きなセッティングはこんな感じだ。

 coverm contig -r [reference_contigs.fasta] -m mean covered_fraction variance length count reads_per_base rpkm -1 [read_R1.fq.gz] -2 [read_R2.fq.gz] > [output_name.tsv]

  「bam要らないって言ったけどごめんやっぱ欲しい!」というときは、--bam-file-cache-directory オプションで作業過程のbamを消さずに出力させることもできる。なので欲しいのがbamだろうがカバレッジだろうが、マッピングするならとりあえずCoverMでやってしまうのが楽なのではないかと思う。ちなみにすでにbamを持っている場合、bamをインプットにしてカバレッジを計算することもできる。またマッピングの条件を具体的に設定したい場合は、minimap2やBWAに渡すパラメータも細かく設定できる。他にも色々できることがあるので、詳しくはcoverm contig --full-helpの出力を参照。
 さらに今回は紹介しなかったけど、contig サブコマンドの代わりにgenomeサブコマンドを使うことで、 複数のコンティグをひとまとめにした「ゲノム」単位でカバレッジを計算することもできる。詳細はcoverm genome --full-helpの出力を参照。
 個人的にはseqkitやfastpと並んで感動した解析ツールにランクイン。有り難く使わせてもらいまくっている。