bcftools機能メモ-01

1. Contig情報の付加

VCFファイルを扱っているときに以下のようなwarningが出ることがあります。

[W::vcf_parse] Contig 'chr1' is not defined in the header. (Quick workaround: index the file with tabix.)

これはVCFファイルのヘッダーに以下のようなContig情報がない場合に出ます。

##contig=<ID=chr1,length=123456789>

メッセージにあるように、 tabix、または bcftools index でインデックスを作成すれば解決します。 indexを作成すると、VCFファイル自体は内容に変更はないものの、 bcftools view などでVCFを読み込み出力すると、 以下のようにContig情報が付加されていることが確認できます。

$ ls
input.vcf.gz  input.vcf.gz.csi
$ bcftools view input.vcf.gz | grep 'contig'
##contig=<ID=chr1>

また、明示的にcontig情報をVCFファイルに挿入したい場合は、以下のようにリファレンスのfaiファイル( samtools faidx で作成可)を指定することで可能です。ただし、既存のContig情報と重複がある場合は上書きされてしまうようなので注意が必要です。また、出力形式は変更できないので、入力と同じ形式で出力されます。結果をパイプに流して less などで確認しようとするとテキスト形式でなくて正しく表示されないといったこともありますので、その場合は適宜 bcftools view を間に挟むなどの必要があります。

bcftools reheader でContig情報を付加した場合はindexによる方法と異なり、length情報も付加されます。

$ bcftools reheader --fai <reference.fai> <input.vcf> -o <output.vcf>

2. InDel付近のSNPを除外する

アンプリコンシーケンスなどで、SNPをターゲットとしたプライマーを設計するときなど、SNPの近傍にInDelがあると不都合なことがあります。そういった場合に、近くにInDelが存在するSNPを除外する方法です。以下の例では近傍100bp以内にInDelが存在するSNPを除外しています。また、|でつなげている下流のコマンドによって最終的にSNPのみのVCFを生成しています。

$ bcftools filter --SnpGap 100 <input.vcf> | bcftools view -v snps

3. gtcheckによるサンプル間の遺伝子型の不一致チェック

bcftools gtcheckを使うと、VCFファイルに含まれるサンプル間の遺伝子型の不一致をチェックすることができます。以下のように、VCFファイルを指定するだけで、遺伝子型の不一致がある場合はその情報を出力します。

bcftools gtcheck <input.vcf>

出力には以下のような情報が含まれます。

行頭の識別子 内容
INFO いくつのサイトを比較に使ったか、どのFORMATタグのデータを比較したかなどの情報
DCv2 discordance情報。クエリサンプル名、比較対象サンプル名、Discordance(不一致の指標)、平均-log P(HWE)、比較したサイト数、遺伝子型の一致したサイト数

VCFをひとつだけ指定した場合はそのVCFファイル内のサンプルの総当たりで比較が行われます。総当たりではなく一部のペアに対してだけ比較を行いたい場合は、-pオプションを使って、比較したいサンプルのペアを指定することができます。以下のように、 -p オプションの後に比較したいサンプル名をカンマ区切りで指定します。

bcftools gtcheck -p sample1,sample2,sample3,sample4 <input.vcf>

上記の場合はsanple1 vs sample2, sample3 vs sample4の2つの比較だけが行われます。

また、比較対象となるサンプルが別のVCFファイルにある場合は以下のように、 -g オプションを使って指定することができます。

bcftools gtcheck -g <input2.vcf> <input1.vcf>

参考