bedファイルの領域を元にvcfのVariantを抽出

ポイント

bedtools を利用する方法
- 出力をvcfファイルとして利用するにはheaderを後から付与する必要がある
bcftools を利用する方法
- vcfファイルをbgzipで圧縮し、tabixでインデックスを作成しておく必要がある

まずはサンプルとなるvcfファイルの取得

$ curl -L -O https://raw.github.com/vcflib/vcflib/blob/master/samples/sample.vcf

bed ファイルとしては下記を利用する(tab 区切り)

region.bed

19	111	112
20	14369	14400
20	1230237	1234567

bedのスタートポジションは0-base

bedtoolsを利用して抽出

sample.vcfファイルのvariantからregion.bedの領域に含まれるvariantを抽出してみる

bedtoolsのintersectコマンドを利用する

$ bedtools intersect -a sample.vcf -b region.bed

結果は次の通り

19	112	.	A	G	10	.	.	GT:HQ	0|0:10,10	0|0:10,10	0/1:3,3
20	14370	rs6054257	G	A	29	PASS	NS=3;DP=14;AF=0.5;DB;H2	GT:GQ:DP:HQ	0|0:48:1:51,51	1|0:48:8:51,51	1/1:43:5:.,.
20	1234567	microsat1	G	GA,GAC	50	PASS	NS=3;DP=9;AA=G;AN=6;AC=3,1	GT:GQ:DP	0/1:.:4	0/2:17:2	1/1:40:3

bed のstart position の指定は 0-baseなので、vcfのchromosome:19 position:111の変異は抽出されていないこと等に注意

注意点

headerをそのまま出力する事はできないので、出力されたvariantをvcfファイルとして取り扱うためには、元のvcfファイルのheaderを利用する必要がある

コマンドシェルにbashを利用しているなら、その助けを借りて次のようにしてヘッダー付与する事ができる
行頭のvcf_file="処理したいvcfファイル"として実行

vcf_file='sample.vcf'; cat <(grep '^#' ${vcf_file}) <(bedtools intersect -a ${vcf_file} -b region.bed)

出力結果

## fileformat=VCFv4.0
## fileDate=20090805
## source=myImputationProgramV3.1
## reference=1000GenomesPilot-NCBI36
...途中省略...
# CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA00001	NA00002	NA00003
19	112	.	A	G	10	.	.	GT:HQ	0|0:10,10	0|0:10,10	0/1:3,3
20	14370	rs6054257	G	A	29	PASS	NS=3;DP=14;AF=0.5;DB;H2	GT:GQ:DP:HQ	0|0:48:1:51,51	1|0:48:8:51,51	1/1:43:5:.,.
20	1234567	microsat1	G	GA,GAC	50	PASS	NS=3;DP=9;AA=G;AN=6;AC=3,1	GT:GQ:DP	0/1:.:4	0/2:17:2	1/1:40:3

そのままファイルに出力するなら

vcf_file='sample.vcf'; cat <(grep '^#' ${vcf_file}) <(bedtools intersect -a ${vcf_file} -b region.bed) > result.vcf

result.vcfに処理結果が出力される

bcftools を利用して抽出

sample.vcf.gz (及びそのインデックスファイルのsample.vcf.gz.tbi）ファイルのvariantからregion.bedの領域に含まれるvariantを抽出してみる

bcftoolsを利用するにはbgzipで圧縮され、tabixでインデックスされたvcfファイルを利用する必要がある

ここでは念の為圧縮してインデックスを付与する方法についても書いておこう
sample.vcfを処理する手順は次の通り

bgzip sample.vcf
tabix sample.vcf.gz

bgzipコマンドでsample.vcf.gzが作成され、tabixコマンドでsample.vcf.gz.tbiが作成される

準備ができたらbcftoolsを実行

bcftools view sample.vcf.gz -R region.bed

結果は次の通り

## fileformat=VCFv4.0
## FILTER=<ID=PASS,Description="All filters passed">
...途中省略...
## bcftools_viewVersion=1.7+htslib-1.7-2
## bcftools_viewCommand=view -R region.bed sample.vcf.gz; Date=Fri Mar  1 10:34:59 2019
# CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA00001	NA00002	NA00003
19	112	.	A	G	10	.	.	GT:HQ	0|0:10,10	0|0:10,10	0/1:3,3
20	14370	rs6054257	G	A	29	PASS	NS=3;DP=14;AF=0.5;DB;H2	GT:GQ:DP:HQ	0|0:48:1:51,51	1|0:48:8:51,51	1/1:43:5:.,.
20	1234567	microsat1	G	GA,GAC	50	PASS	NS=3;DP=9;AA=G;AN=6;AC=3,1	GT:GQ:DP	0/1:.:4	0/2:17:2	1/1:40:3

出力先にファイルを指定するなら-oオプションを利用する

bcftools view sample.vcf.gz -R region.bed -o result.vcf

result.vcfに処理結果が出力される

利点

ヘッダーもちゃんと処理できる

注意点

vcfファイルは予めbgzipで圧縮し、tabixでインデックスを付与しておく必要がある

まとめ

このようにして、bedtoolsでもbcftoolsでも必要な領域のvcfファイルを抽出する事ができた

各ツールをubuntuに導入する方法

実行環境は以下の通り

$ grep 'VERSION' /etc/*release 
/etc/os-release:VERSION="18.04.1 LTS (Bionic Beaver)"
/etc/os-release:VERSION_ID="18.04"
/etc/os-release:VERSION_CODENAME=bionic

基本的にはapt-getでインストールすれば良い

bedtools

# apt-get update
# apt-get install -y bedtools

tabix(bgzipを含む), bcftools

# apt-get update
# apt-get install -y tabix
# apt-get install -y bcftools

今回はここまで