アセンブルされた2つの配列を比較してVCFを出力する

これまで私は、VCFを出力するケースといえば、リファレンスに対してリードをマッピングしてバリアントコーラーを走らせるという場合しか考えていませんでした。 一方で、アセンブルされた配列同士をアラインメントするときは、MummerやMinimap2を用いており、出力されるのは基本的にはアラインメント結果のファイルでした。 ここで、アラインメントの結果から、配列の異なる部分をVCF形式で出力する方法が必要になったので検索してみると、GATKにその機能があるようなので使ってみました(4.3.0.0以降?)。

VCFファイルを出力するためには、gatk CompareReferencesを使用します。今回試したGATKのバージョンは4.4.0.0です。 このツールは現段階でEXPERIMENTALです。

# a.faをリファレンスとして、b.faと比較する
samtools faidx a.fa
samtools faidx b.fa
gatk CreateSequenceDictionary -R a.fa
gatk CreateSequenceDictionary -R b.fa
mkdir output
gatk CompareReferences -R a.fa -refcomp b.fa --base-comparison FULL_ALIGNMENT --base-comparison-output output

このコマンドを実行すると、outputディレクトリにa.fa_b.fa.vcfというファイルが出力されます。 このファイルには、リファレンスと比較したときの差異がREF,ALT列にVCF形式で格納されています。

最近はVCFを入力にとる便利ツールも増えてきている印象なので、この機能は意外と使える気がしています。