GATK FastaAlternateReferanceMakerを使って代替参照配列を作成する
GATK FastaAlternateReferanceMakerの使い方
使用難易度★★★★☆
本記事は、GATK解説シリーズのPart 11です。
GATK解説シリーズのリンクまとめは↓こちら
GATKの導入方法から、変異情報の取得までをハプロイドの病原体を使った実例とともに紹介しています。
eupatho-bioinfomatics.hatenablog.com
今回は何をする？
GATK FastaAlternateReferanceMakerを使って、元となった参照配列に変異情報を当てはめる形で代替リファレンス配列(Fasta形式）を取得する。
感覚的にはVCF→FASTAに変換しているが、実際に変換したわけでは無く、VCFファイルに含まれる代替対立遺伝子(alternate allele）の配列をリファレンス配列と交換することで、変異を取り込んだ新しいFASTAを得ている。 ゲノムを新たにアセンブリしたわけではないので、利用には少し注意が必要。
自分で使用して失敗してしまった注意点もお伝えしますので、初めて使う方はぜひ読んでから開始することをお勧めします！
GATK Forumの説明を読んだだけで実行すると、十中八九同じ失敗をすると思います。
fastaファイルをコマンドラインで操作するseqkitも少しだけ紹介する。
vcftoolsもしくはbcftoolsでも似たような機能が利用できる。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。 フォローする @harrykun_blog
GATK公式HPのFastaAlternateReferanceMakerのリンク
1. FastaAlternateReferanceMakerの概要
参照配列の塩基を、対応するVCFファイルで提供された代替対立塩基に置き換えます。
また、-L引数で代替参照配列を生成する区間を指定することができます。
出力フォーマットは、fasta配列として出力されます。
注意点
1つのサイトから始まるバリアントが複数ある場合、そのうちの1つをランダムに選択します。
インデルが重なっている場合（開始位置が異なる場合）、最初のインデルのみが選択されます。
このツールはSNPと単純なインデルに対してのみ動作します（複雑な置換のようなものには動作しません）。
入力として
リファレンスのゲノム配列、VCFファイル
出力として
選択されたサンプルの変異を取り込んだ含む新しい FASTA ファイルを得る。
2. FastaAlternateReferanceMakerの使い方
gatk FastaAlternateReferenceMaker \ -R reference.fasta \ -O output.fasta \ -L input.intervals \ -V input.vcf \ --snp-mask
各オプションについて
-O/--output; 出力されたFASTAファイルを書き込むためのパス
-R/--reference; 参照配列ファイル
-V/--variant; 参照配列に適用するためのVCFファイル
-L/--intervals; 解析対象となる1つ以上のゲノム区間
--snp-mask; このオプションを指定すると、代替リファレンスを構築する際にオーバーラップする塩基を「N」に設定する（配列にNを挿入する）。
今回使用するVCFファイルについて
解析に使用する自分のデータがない方は、本ブログのGATK解説シリーズで作成した25株のPlasmodium falciparum（マラリアを起こす病原体）の全ゲノムデータを使って検出したSNPファイルをダウンロードしてください。
こちらのリンクからダウンロードできます（200MBほど）
本ブログでいつも使用しているPlasmodium falciparum 3D7のリファレンスゲノムが必要となるので、GWASと名付けたフォルダの中にvcfと3D7_genomesというフォルダを用意して、merged.vcfファイルとリファレンスゲノムをそれぞれ格納した状態で開始する。（リファレンスゲノムの入手方法は別の記事で紹介しています。）
3. GATK FastaAlternateReferanceMakerを実行する。
FastaAlternateReferanceMakerの罠？について
自身の解析でFastaAlternateReferanceMakerを使用した際に、失敗した点です。
FastaAlternateReferanceMakerに使用するVCFファイルは、JointGenotyping後のマージVCFファイルであることが必要。
GATK HaplotypeCallerで生成した直後の個別のg.vcfファイルは使用できません。
さらに、SelectVariantsで変異を持たないサイトを取り除いておかないと、本来無関係なバリアントも新しいFASTAに含まれてしまう。
よく考えてみると、g.vcfファイルはALTフィールドに塩基が書いてないので、ALTを参照するだけのFastaAlternateReferanceMakerに注意がいるのは納得でしたが、GTフィールドの結果を反映してくれないことは盲点でした・・・。最近、この失敗で3週間ほど棒にふりました😭。
FastaAlternateReferanceMakerを使用する前に、GATK SelectVariantsで欲しいサンプル（今回はERS347575とする）のサブセットを抽出したファイルERS347575_Genome.vcfを用意します。
GATK SelectVariantsについての記事は↓こちらから
eupatho-bioinfomatics.hatenablog.com
次に、ERS347575_Genome.vcfを使って、FastaAlternateReferenceMakerを実行します。
cd /path to GWAS
index=/path to GWAS/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta
gatk SelectVariants \
-R ${index} \
--variant merged.vcf \
--output ERS347575_Genome.vcf \
--exclude-non-variants \
-sn ERS347575
gatk FastaAlternateReferenceMaker \
-R ${index} \
-O ERS347575_Genome.fasta \
-V ERS347575_Genome.vcf \
-L Pf3D7_01_v3
-Lオプションを指定しない場合、全ゲノムが出力される。
完了すると、ERS347575_Genome.fastaが得られる。
3. 生成された新しいFastaに変異情報が取り込まれたかどうか確認する。
まず、この二つの画像を見比べてみてください。
上側が--exclude-non-variants無しで出力したVCFファイル。
下側が--exclude-non-variants有りで出力したVCFファイル。
--exclude-non-variants無しで出力したファイルには、ERS34757のGTフィールドに0や.と記載されたサイトが含まれていますが、
--exclude-non-variants有りで出力したファイルには、GTフィールドに1（つまり対立遺伝子アレルを持つ）バリアントサイトだけが保存されています。
FastaAlternateReferenceMakerでは、Sampleフィールドの情報を参照する事無くALTとREFを交換するため、不必要なバリアントサイトはあらかじめ除外することが重要です。
- 最後に、出力されたFASATAを参照配列と比較してみましょう。
新しいFASTAには１番染色体だけが含まれるので、アライメントするために、参照配列のゲノムから不要な染色体を取り除きます。
ここでは、FASTAファイルの操作を簡単・迅速に実行できるseqkitを使用します。
1. seqkitをインストールする
conda install -c bioconda seqkit
- 目的の染色体だけを取り出す。
seqkit faidx PlasmoDB-52_Pfalciparum3D7_Genome.fasta -r Pf3D7_01_v3 > Reference.fasta
下図は、新しいFASTAと抽出した参照配列を比較したものです。
上にある下側の図に記載された上から２つのSNPsが反映されていることが分かりました。
お疲れ様でした。今回はこれで終わりです。
よければ他の記事のも見ていってください。
系統解析について解説した記事は↓こちら
eupatho-bioinfomatics.hatenablog.com
eupatho-bioinfomatics.hatenablog.com
VCFファイルについて解説した記事は↓こちら
eupatho-bioinfomatics.hatenablog.com
バイオインフォマティクス関連の書籍紹介は↓こちら
