vcf 파일로부터 chromosome 별로 variants 의 밀도와 분포도 표시하기

bowtie2 와 picard, GATK 를 이용하여 BAM 파일을 만들고 이것을 다시 재정리하여 vcf 파일을 만들었다. 이것을 snpEffector 를 이용하여 분석을 하였지만 window size 가 너무 큰 관계로 chromosome 별 밀도와 분포를 자세히 표시할 수가 없었다.
window size 를 10kb 로 만드는 방법을 알아내었다.
해답은 무척 쉬운데 엄청 헤멨다.

방법은 bedtools 를 이용하는 것이다. 먼저 reference 의 각 chromosome size 를 정리하여야 한다.

[chromosome name]<TAB>[size] 로 정리 되어 있어야 한다.

1       43270923
2       35937250
3       36413819
4       35502694
5       29958434
6       31248787
7       29697621
8       28443022
9       23012720
10      23207287
11      29021106
12      27531856

정리해보면 위와 같다. 내용을 파일로 저장하자. 그 다음..

bedtools makewindows -g filename -w 10000 --> 10000 으로 windows 작성하는 것이다.

그러면 결과가 다음과 같이 나온다.
1 0 10000
1 10000 20000
1 20000 30000
1 30000 40000
1 40000 50000
1 50000 60000
1 60000 70000
1 70000 80000
1 80000 90000
.................

위의 결과를 다시 저장해둔다.

bedtools intersect -a filename -b vcf_filename -c

위를 실행하면 다음과 같이 뜬다.
1       0       10000   0
1       10000   20000   0
1       20000   30000   0
1       30000   40000   0
1       40000   50000   0
1       50000   60000   0
1       60000   70000   0
1       70000   80000   0
1       80000   90000   0
1       90000   100000  0
..............

chromosome name, start position, end position, count 개수 순서이다. '-c' 옵션은 -a 의 범위내에 -b 가 위치할 경우 그 개수를 세어주는 것이다. 여기에 -f 혹은 -F 옵션을 사용하여 최소 counts 개수를 지정할 수 있지만 필요가 없는 관계로 자세한 설명은 넘어간다.