bcftoolsでvcfファイルの編集(個人メモ・追記予定)

※2023年7月19日 一部修正

VCF形式

詳しい説明は省略。GATK公式による解説がおそらく一番分かりやすい。 VCFファイルを整形したり、ほしい多型情報を抽出するツールとしてはbcftoolsとVCFtoolsが定番となっている。個人的にはbcftoolの方が、

  • biocondaで配布している
  • bamファイルからの多型抽出コマンドも備えていて、多型抽出からフィルタリングをシームレスに行える

という点で扱いやすいと感じた。

bcftoolsのインストール方法

anaconda環境があればbiocondaからインストールできる。

conda install -c bioconda bcftools

もしくはgithubからインストールする。

VCFの連結

bcftools concat を使う。引数-aをつけると最後のファイルので1行目を上書き。サンプル名、フィールド等が同じであればこれをつける。

# make list
echo -e "file-1.vcf.gz\nfile-2.vcf.gz\nfile-3.vcf.gz" > vcflist

# concatenate tables
bcftools concat -a -f vcflist -o vcf_conc.vcf.gz

# add index
bcftools index -f vcf_conc.vcf.gz

VCFの情報を確認

# list samples
bcftools query -l file.bcf

# number of samples
bcftools query -l file.bcf | wc -l

# list of positions
bcftools query -f '%POS\n' file.bcf

VCFのフィルタリング

-i '条件'を引数に取ることで、条件に合う多型サイトを抽出できる。条件に合うものを除外するには-e '条件'とする。 複数サンプルのmulti-vcfファイルの場合、一つでも条件を満たすサンプルがあれば、そのサイトは抽出される。全サンプルが満たす条件で抽出する場合には、-eを使う。 パイプで繋げる時は-Ou (output uncompressed BCF)をつけると早い。

(例) 1. FILTERフィールドがPASSであり(GTAKなどで事前にフィルターをかけていた場合)、かつQAUL (quality)が100以上であり、 2. 1~4番目の個体のうちgenotypeが0/0である個体が3個体以上の多型を抽出

bcftools view -O u -i 'FILTER="PASS" & QUAL>100' file.vcf.gz |
bcftools view -i 'COUNT(FMT/GT[0-3]="0/0")>2' -o file_filtered.vcf

※フィールドのINFOとFMTの違い INFOは全サンプルを総合した数値。FMT(FORMAT)は個別サンプルにつけられる値を表す。たとえばINFO/DPは全サンプルのdepth合計だが、FMT/DPとすると1サンプルごとのdepthを示す。 ※&と&&の違い &で条件を繋げると、全条件が一つのサンプル内で満たされなければいけない。&&で繋げた場合、各条件がいずれかのサンプルで満たされていればよい。

VCFをテーブルに変換

bcftools query で必要な情報を取り出してテーブルに保存できる。-f '条件' の形式。-H:ヘッダーをつけて保存

(例) タブ区切りでCHROM, POS, REF, ALT, QUAL, AC, AFの値、その右にGT, AD, DPの値を個体ごとにタブで区切って出力、テキスト形式で保存。

bcftools query -H -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%AC\t%AF[\t%GT %AD %DP]\n' file_filered.vcf.gz > file_filtered_GTADDP.txt

ジェノタイプ変換

setGTプライグインが提供されている。ただしデータの変換はvcfにあった情報を破壊することになるので慎重に行うべき。一度保存してから外部のツールを使うべきかもしれない。 (例) GQ (Genotype Quality: 個体ごとのGTの確からしさ)が20未満の場合、そのGTを./.に変換

bcftools +setGT input.vcf.gz -o transformed.vcf.gz -- -t q -i 'GQ<20' -n "./."

参考サイト

http://samtools.github.io/bcftools/howtos/filtering.html

。これに関しては情報が少ない。