[BCFtools] 유전변이(variant) 찾기 - VCF 만들기 - BCFtools를 이용한 방법

VCF의 유전변이(variant) 결과를 염기서열(DNA sequence)에 반영하는 방법은 다음과 같다.

선행되어야 할 것은 

1. 참조유전체(reference genome)

2. BCFtools

3. 맵핑되고, sorting 되어 있는 BAM과 indexing된 BAI

해당 3가지 정보가 확보되면 다음과 같이 서열을 얻게 된다.

bcftools로 변이를 찾고 vcf로 바꾸는 과정을 linux pipe ('|') 이용하여 한번에 진행하게 된다. 아쉽게도 아직은 bcftools에서 mpileup을 얻는 과정에 병렬 처리(parallel)가 되지 않는다. 

[s_ref]: 참조유전체 파일 (fasta 형식)
[bam]: binary 형태의 맵핑 파일 
[sample]: VCF 파일에서 이용할 샘플의 이름 
-Ou: 리눅스에서 pipe ('|')로 이용하기 위한, uncompressed BCF로 내보내는 output
-f: 참조유전체 input을 받기 위해 입력
-vmO: v는 variant만 출력한다는 의미이며 변이가 없는 mpileup 영역은 출력하지 않음
      m의 경우에는 여러 allele에 대한 정보를 출력한다는 의미이고, 
      O는 뒤의 vcf 형태로 출력한다는 것을 위해 넣음

mpileup이 genomic position에서 read의 수치적 특성을 나타내는 것이고
call이 mpileup input을 이용하여 유전변이를 찾아주는 명령어 이다.

bcftools mpileup -Ou -f [s_ref] [bam] |bcftools call -vmO v -o [sample].vcf

mpileup의 역사는 길다

추후 아래 링크를 설명하면서 이 코드의 탄생에 대하여 알아보도록 하겠다.

The original mpileup calling algorithm plus mathematical notes (mpileup/bcftools call -c):

• Li H, A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data, Bioinformatics (2011) 27(21) 2987-93. [21903627]
• Li H, Mathematical Notes on SAMtools Algorithms (2010) [link]