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 区切り)
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
今回はここまで