これまで私は、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を入力にとる便利ツールも増えてきている印象なので、この機能は意外と使える気がしています。