vcfファイルのサンプル名を変えたい

vcfファイルのサンプル名がbamファイル名なんですけど…

bcftoolsで何も考えずにバリアントコールをしたとき、vcfのサンプル名はbamファイル名になっています。 別にこのままでもいいかもしれませんがなんかダサい感じがします。

そこでvcfの列名の部分にサンプル名が来るようにしたいと思います。バリアントコールは以下のように行う ことにします。

bcftools mpileup -f reference.fa sample1.bam sample2.bam sample3.bam | bcftools call -mv -Ob -o calls.bcf

ここで、vcfファイルのヘッダーに対し、任意のサンプル名をつけるためには bcftools mpileup コマンドにおいて追加で -G または --read-groups というオプションを用います。

-G オプションはファイル名を引数にとり、このファイルにbamファイルとサンプル名の対応関係を記録しておくこと でその情報を bcftools mpileup にわたします。

この際に指定するファイルの形式は以下のようになります(ただし、1つのbamファイルに1つのサンプル由来の リードのみ含まれているとします)。

* sample1.bam sample1
* sample2.bam sample2
* sample3.bam sample3

1列目には*を書きます。これは1つのbamファイルに含まれるリードはすべて1つのサンプル由来であることを 表しています。もしそうでない(1つのbamファイルに複数サンプルのリードが含まれる)ときは リードグループIDを1列目に指定することになります(詳しくは公式マニュアル参照)。 続いて2列目にはbamファイルのファイル名を指定し、 それぞれに対応するサンプル名を3列目に記載するといったファイル形式です。区切り文字は空白を用います。

上記のファイルを read_groups.txt としてカレントディレクトリに保存した場合、具体的な コマンドラインは以下のようになります。

bcftools mpileup -f reference.fa -G read_groups.txt sample1.bam sample2.bam sample3.bam | bcftools call -mv -Ob -o calls.bcf

これでvcfファイルの列名がちゃんとサンプル名になったと思います。

すでにあるvcfのサンプル名を変えたいんですけど…

もしかしたら、もうすでにあるvcfのサンプル名を変えたいといった状況に遭遇するかもしれません。 その場合は bcftools reheader が有用です。-s または --samples オプションを用います。

-s オプションも同様にファイル名を引数にとります。このファイルには各行ごとにサンプル名を 記載するだけで大丈夫です。このとき、vcfの列順とファイルに記載するサンプル名の順が一致していること が必要になります。

以下に例を示します。

sample1
sample2
sample3

上記のようなファイルを samples.txt というファイル名でカレントディレクトリに保存したら、 以下のコマンドを実行します。

bcftools reheader -s samples.txt -o new.bcf old.vcf

次回はこのあたりを踏まえ、マッピングからバリアントコールまでのパイプラインを snakemake を用いて 記述していきたいと思います。