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

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

## ポイント

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

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

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

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

```text: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
```

結果は次の通り

```text
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ファイル"として実行

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

出力結果

```text
##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
```

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

```bash
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を処理する手順は次の通り

```terminal
bgzip sample.vcf
tabix sample.vcf.gz
```

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


準備ができたらbcftoolsを実行

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

結果は次の通り

```text
##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`オプションを利用する

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

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

利点

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

注意点

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


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


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

実行環境は以下の通り

```terminal
$ 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
```terminal
# apt-get update
# apt-get install -y bedtools
```

### tabix(bgzipを含む), bcftools

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

今回はここまで:smile: