4.2. GATK HaplotypeCaller 수행으로 genome VCF 만들기 (Make the genome VCF by using GATK HaplotpeCaller)

오랜만에 본업으로 글을 작성하게 되었다.

GATK germline variant call의 경우 올해 안에 완성하도록 노력하겠다.

주업과 취미를 고루 잘하는 블로거가 되도록 하겠다.

우선 small variant ( < 50 bp variant) 에 대하여 먼저 언급하겠다.

1.아니 Human genomics와 Next generation sequencing의 연구에 대하여 이야기 해보겠다.

현재의 Cancer Moon Shot project (https://www.cancer.gov/research/key-initiatives/moonshot-cancer-initiative) 와 Precision Medicine Initiative (https://obamawhitehouse.archives.gov/precision-medicine) 이전에

인간 유전체 지도를 만드는 큰 노력이 있었다.

2. 이 때 유전체를 깬 후 작은 조각을 읽어내고, 그 정보를 통계적으로 그리고 de bruijn graph (미국 시민권자에게, 그리고 google에 물어본 결과 "드 브론" 식으로 읽어야 한다고 했다. 브루진은 절대 아니고 브루인도 아니다!!!) 를 통해서 유전체 전체 분석이 가능하다는 아이디어를 제시한 George Church 교수가 있었다.

이러한 아이디어로 서열 분석기 (sequencer)에 대한 새로운 기술이 시작되었다. 그리고 인간 유전체 지도 중 대표성을 지닌 참조 유전체 (reference genome) 도 완성되었다. 

2.1. 첫번째 기기 기술의 발전에 있어서

이에 따라 Sanger가 제시한 서열 분석방법인 Sanger sequencing이 아닌 차세대 서열 분석기가 발명되기 시작한다 

2.2. 초기부터 2010년 까지의 역사를 가장 잘 정리한 논문은 아래의 논문이다. 

이 분야 공부할 꺼면 10번은 읽자.

Metzker, M. L. (2010). Sequencing technologies - the next generation. Nature Reviews. Genetics, 11(1), 31–46. http://doi.org/10.1038/nrg2626

 

두번째 기기 기술의 응용과 관련하여 기준점이 되는 참조유전체의 이용에 있어서

짧은 DNA 서열의 해독, DNA 서열을 기준에 대한 좌표화, 좌표에 대하여 대립형질의 유무를 분석하는 프로그램이 나오게 되었다.

3. 이와 관련하여서 대립형질의 유무를 알기 위해서 참조유전체와 차이가 있는 유전변이 (variant)를 추출하는 방법, 그리고 저장하는 방법에 대한 발전이 진행된다.

3.0.1. 초기의 단순 서열의 좌표화 --> SAM format
(Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., … Durbin, R. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078–2079. http://doi.org/10.1093/bioinformatics/btp352)

3.0.2. 좌표에 대한 집단적 유전변이의 파악 --> VCF format

(Danecek, P., Auton, A., Abecasis, G., Albers, C. a., Banks, E., DePristo, M. a., … Durbin, R. (2011). The variant call format and VCFtools. Bioinformatics, 27(15), 2156–2158. http://doi.org/10.1093/bioinformatics/btr330)

두 가지 표준 포맷이 탄생하였다.

3.1. 초기에 만들어진 유전변이 발굴 소프트웨어는 

지금도 선두그룹 중 하나인 Wellcome Trust Sanger Institute의 samtools (http://www.htslib.org/) 

아직도 빠르게 분석된다.

간단한 분석은 충분하다.

3.2. 그리고 추후에 집단 분석이 가능해진, Bayesian statistics가 도입되어 통계적으로 quality 확인이 명확해진 Broad institute의 Unified Genotyper (GATK의 시작; DePristo, M. a, Banks, E., Poplin, R., Garimella, K. V, Maguire, J. R., Hartl, C., … Daly, M. J. (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics, 43(5), 491–498. http://doi.org/10.1038/ng.806)

대부분의 유전변이 발굴은 이와 비슷한 통계적 분석이 강화된 프로그램이다.

3.3. 세 번째로 local re assembly 기술이 가미된 FreeBayes (Boston College, 원조맛집), Platypus (Wellcome Trust Sanger Institute; Rimmer, A., Phan, H., Mathieson, I., Iqbal, Z., Twigg, S. R. F., Wilkie, A. O. M., … Lunter, G. (2014). Integrating mapping-, assembly- and haplotype-based approaches for calling variants in clinical sequencing applications. Nature Genetics, 46(8), 912–918. http://doi.org/10.1038/ng.3036), HaplotypeCaller (현재 최고 맛집, Broad Intstitute)가 있다. 앞서 통계적 기법에다가 부분적인 재조립을 통하여 small variant ( < 50 bp) 를 잘 찾는 방식이다.

이 중 우리는 3번째로 나온 local re assembly 기술이 적용된 최신 유전변이 발굴 프로그램인 Haplotype의 이용법을 알아보자!

4.0. [본문 시작]

haplotype caller를 이용하면 2단계의 분석이 진행된다. 

1. gvcf를 만든다. (variant가 아닌 부분도 read가 있었다라는 정보를 지니게 된다. 추후 multi-sample analysis의 재료로 이용 가능하다.

2. vcf를 만든다. (single sample 분석에 대하여 진행하게 된다., 해당 과정은 선택이 있기에 다음 포스팅에 기록할 예정이다.)

4.1. gvcf generation

준비물
Reference genome
해당 reference로 분석된 VCF
BAM 파일 

표 만들고, BAM 보정한다. 

gatk --java-options "-XX:+UseParallelGC -XX:ParallelGCThreads=8" HaplotypeCaller -R [] -I bam -O g.vcf -ERC GVCF
bgzip -@ 8 -c gvcf > gvcfgz
tabix -p vcf gvcfgz

첫번째 명령어는 BAM을 VCF 형태인 gVCF를 만드는 과정이다. 해당 과정에서 read를 읽었지만, reference genome과 같은 결과를 갖는 부분까지 기록하게 된다.

이러한 기록을 시작하게 된 이유는 

변이가 없는 유전좌(loci)에서 

1. 실제로 reference genome과 genotype이 같아서 정보가 없는 것인지

2. 아니면 read 정보가 없기 때문에 variant 정보가 없는 것인지

해당 사항에 대한 궁금증이 발생하였기 때문이다.

이러한 정보를 BAM에서는 가지고 있으나, VCF로 변환하는 과정에서 소실된다.

비어있는 부분은 InDel의 일종으로 볼 수 있는 것이다.

그렇기 때문에 이미 생산 완료된 각각의 VCF에서 여러 샘플 단위의 group의 VCF를 만들게 되면 InDel 때문인지 아니면 Reference와 같은 genotype이어서 Variant가 없는 가를 궁금하게 되었다. 이에 BROAD에서는 multi sample variant calling을 더 잘하기 위해서 genome vcf를 도입하게 되었다. 

코드의 3 줄의 의미는 다음과 같다.

1. gvcf 만들기

2. gz 형태로 압축 (그래야 후수 과정 분석이 빠르다)

3. vcf.gz.에 대한 index 만들기 (gz를 사용하는 프로그램들이 indexing을 요구한다.)

------------------------

늦은 수정 완료였다.

앞으로 포스팅 빈도를 높이겠다.!