[bedtools 사용방법] bedtools를 이용하여 vcf 의 특정 범위 꺼내기

안녕하세요

한주현 입니다

오늘은 bedtools를 이용하여 vcf에서 필요한 영역을 꺼내는 작업을 해보겠습니다.

물론 스크립트를 작성하여 작업하는 방법도 있겠습니다만,

빌게이츠가 좋아하는 명언 중 하나인,

Do not reinvent the wheel.

(바퀴를 다시 발명하지 말라.)

라는 말 처럼 이번 시간에는 있는 tool을 사용해보도록 하겠습니다

bedtools 설치에 대해서는 다음 링크를 참조해 주세요

[bedtools] bedtools  설치 방법 - How to install bedtools

http://korbillgates.tistory.com/66

여담으로 도전정신이 있으신 분들께서는 아래의 예시를 프로그래밍하여 해결해보시기를 바랍니다 :)

VCF 의 특정 범위 꺼내는 방법

먼저 이번 실습에 사용할 샘플 vcf 입니다

https://raw.githubusercontent.com/KennethJHan/Bioinformatics_smalltalk_Python50/master/sampleVCF.vcf

vcf를 받도록 합시다

sampleVCF.vcf

##fileformat=VCFv4.2
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta ##contig=<ID=20,length=62435964,assembly=B36,md5=f1

... <생략>

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  NA00001
chr20   14370   rs6054257       G       A       29      PASS    NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ     0/1:48:1:51,
chr20   17330   .       T       A       3       q10     NS=3;DP=11;AF=0.017     GT:GQ:DP:HQ     0/1:40:1:43,43
chr20   23145   rs6040355       A       T       40      PASS    NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ     1/1:48:1:51,
chr20   61423   .       T       A       30      PASS    NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ     0/1:48:1:51,51
chr20   63475   .       C       G       25      PASS    NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ     1/1:48:1:51,51
chr20   82739   .       GG      G       27      PASS    NS=3;DP=21;AF=0.5;DB;H2 GT:GQ:DP:HQ     0/1:48:1:51,51
chr20   99384   rs6054299       G       A       29      PASS    NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ     0|0:48:1:51,
chr20   102423  .       ATC     A       37      PASS    NS=3;DP=22;AF=0.5;DB;H2 GT:GQ:DP:HQ     0/1:48:1:51,51

... <생략>

chr20   271029  .       C       T       29      PASS    NS=3;DP=19;AF=0.5;DB;H2 GT:GQ:DP:HQ     1/1:48:1:51,51

vcf 내부는 이렇게 생겼는데요,

저희는 여기서 chr20의 60000 에서 90000 까지의 범위 내의 vcf라인을 가져와보겠습니다

여러가지 방법이 있겠지만,

bedtools의 intersect를 사용하겠습니다

1. 필요한 영역을 bed range로 만든다

vi 와 같은 text editor로

chr20    60000    90000

을 입력하고 저장합니다.

전 여기서 range.bed로 저장하였습니다.

이 파일의 구분자는 탭으로 합니다.

2. bedtools intersect

$ bedtools intersect -a sampleVCF.vcf -b range.bed

chr20   61423   .       T       A       30      PASS    NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ     0/1:48:1:51,51
chr20   63475   .       C       G       25      PASS    NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ     1/1:48:1:51,51
chr20   82739   .       GG      G       27      PASS    NS=3;DP=21;AF=0.5;DB;H2 GT:GQ:DP:HQ     0/1:48:1:51,51

너무나도 쉽게 바로 결과가 뚝딱 나와버렸습니다 ㅎㅎ...

이상으로 vcf에서 원하는 영역 꺼내기에 대한 포스팅을 마치겠습니다.

그럼 다음 시간에 만나요~