Posted at

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


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

今回はここまで:smile: