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

今回はここまで

BEDファイルの領域に含まれるVariantを抽出する