바이오인포매틱스

GATK의 incremental joint discovery를 위한 reference model pipeline

hongiiv 2014. 3. 11. 16:58
반응형

좀 더 많은 genomes이 필요한 시대

$1000 게놈 시대가 진짜 도래했다. 이제까지 NGS 연구의 대부분이 하나의 genome 데이터를 가지고 연구(rare variant를 찾던)하던 것이 GWAS처럼 대규모의 cohort의 샘플을 수용하기 시작하면서 "Common Variant Association Study (CVAS)"에 눈을 돌리기 시작했다. 이는 가격뿐만 아니라 대량의 NGS 데이터를 다루기 위한 툴 또한 발전하면서 가능케 되었다.

이러한 CVAS 데이터는 cohort의 샘플들(individucal callsets)을 개별적으로 variant call을 하는 것이 아니라, joint callset을 만들어 joint variant discovery를 수행하여 흔히 말하는 power를 부여할 수 있다. 각각의 샘플 Sample #1을 분석해서 사용한다면 "A" site에 대해서 진짜 변이인지에 신뢰할 수 없다 하지만, Sample #1~Sample #N까지 모든 샘플을 보면 cohort내에서의 "A" site에 대한 real variation에 대한 좀 더 신뢰할 수 있게 된다.


또한 join discovry를 통해 bias 이슈 또한 어느정도 해결이 가능하다. 하나의 샘플만을 본다면, strand에 따라 FWD reads는 모두 REF와 같은 A(4개)를 REV reads는 REF allele A가 5개  VAR allele인 C가 3개 즉, allele ratio = 9: 3 = 3:1이 된다. 


자, 이제 joint analysis에서는 어떻게 systematic biases를 줄이는지를 보면, 여러샘플 데이터를 통해 allele bias와 strand bias를 filter out 할 수 있게 된다. 이것은 위와 같이 single sample에서는 해결 할 수 없는 부분이었다. 


근데 문제는 있다.

첫째 이렇게 joint analysis를 하기 위해서는 모든 site에 정보가 있어야 한다. 하나의 샘플만을 분석한다면 단지 variant에 대한 정보만 존재하면 되지만 joint calling을 위해서는 모든 샘플의 모든 site에 대한 정보를 필요로 하게 된다. 둘째, 현재까지의 접근으로 multisample을 분석할때 매우 comutationally intensive한 작업으로 computing cost와 "N+1 문제"에 봉착(모든 샘플을 동시에 분석하기 때문에 샘플이 하나 더 증가하면 +된 1개의 샘플을 포함해서 다시 N+1개의 샘플을 다시 분석해야하는...)하게 된다. 



기존의 여러개의 BAM(Analysis-Reday Reads)을 joint variant calling을 하는 경우 N+1 문제가 발생합니다. 물론 Reduced Reads를 이용해서 어느정도 computing resource에 대한 해결을 보긴 했지만,,, 새로운 방법은 각 sample별로 variant calling을 합니다. 이때 HC(Haplotyper Caller)의 gVCF mode를 통해 각 sample별 variant calling을 수행한 후 개별 gVCF를 통해 Joint Genotyping을 수행하면 기존의 Joint Varaint Calling을 동일한 결과물을 얻을 수 있게 된다.

gVCF

gVCF는 GATK 2.7버전부터 HC에 들어있는 기능으로, --emitRefConfidence (-ERC) 옵션을 통해 사용이 가능하며 이 옵션을 통해 이전에 언급한 모든 site에 대한 정보가 있어야 joint analysis가 가능하다라는 부분을 해결할 수 있다.
  • NONE: don't emit anythong (default) 아무것도 넣지 않는것으로 기본
  • BP_RESOLUTION: emit detailed information for each BP 각 Base Position의 모든 정보를 포함 (ref와 같던 같지 않던 모든 bp에 대한 정보를 포함하며 마치 UG(Unified Genotyper의 --output_mode의 EMIT_ALL_SITES와 유사한)
  • GVCF:  emit a block summarized version of the BP_RESOLUTION data BP_RESOLUTION의 모든 bp에 대한 정보를 넣는 반면 block 단위로 데이터를 압축(summarize)하는 것으로 마치 Reduced Reads가 불필요한(정보성이 없는) read에 대해서 필요한 정보만 남겨 데이터를 압축하는 것처럼, 일정 block이 reference와 같으면 해당 region 정보를 bp레벨로 다 정보를 입력하지 않고 합쳐서 보여주는 형태
마지막 GVCF를 다시 보면, 첫라인의 염색체 20번의 10000435 bp는 ref가 'T'이고 <NON_REF> 즉, 변이가 없고 해당 block은 길이가 4 즉, 10000435~10000438까지가 모두 Reference와 같다는 의미이다. (BLOCK_SIZE=4;END=10000438) 이 block은 대충 GQ가 비슷한 놈들끼리 묶이게 되며 MIN_DP와 MIN_GQ값으로 해당 block의 정보를 표현하게 된다.

20 10000435 . T . . BLOCK_SIZE=4;END=10000438 GT:DP:GQ:MIN_DP:MIN_GQ 0/0:53:99:53:159
20 10000439 . T G 1553.77 . AC=2;AF=1.00;AN=2;DP=56;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=221.40;MQ0=0;QD=27.75 GT:AD:DP:GQ:PL 1/1:0,55:55:99:1582,164,0
20 10000440 . T . . BLOCK_SIZE=6;END=10000445 GT:DP:GQ:MIN_DP:MIN_GQ 0/0:56:99:55:165

그럼, 왜 BP_RESOLTION이 있는데 왜 굳이  gVCF모드를 따로 두었느냐? 바로 데이터의 크기 때문이다. 모든 bp의 정보를 표현하다 보면 데이터가 커지게 되며 이는 곧 cohort내의 모든 샘플을 분석하는데에 비용이나 시간적인 문제로 직결되기 때문이다. 실제로 exome 데이터를 NONE, BP_RESOLUTION, GVF 모드로 각각 수행했을때 결과 파일의 용량을 비교해 보면 gvcf에 비해 bp resolution이 60배 이상 나는것을 확인 할 수 있다. 즉 몇천샘플을 가지고 분석을 한다면 gvcf가 훨씬 데이터 reduction에 의한 효과를 볼 수 있게된다.


gVCF 파일의 joint genotyping

각 샘플별로 HC를 통해 calling된 gVCF (genomic VCF, record every single locus in the genome, 각 샘플은 variant가 아닌 sites에 대해서도 reference에 대한 confidence estimation 정보를 갖게 됨(reference model)) 파일은 마지막으로  GenotypeGVCFs를 통해 하나의 joint된 VCF  파일을 만들게 된다. 이후 recalibration을 수행하면 끝

정리하자면, joint discovery(gVCF와 reference model을 통해)는 추후 지속적으로 늘어나는 샘플에 대해서 incremental하게 함으로써 cost를 줄이고, 기존에는 가능하지 못했던 large cohorts에 대한 study가 가능하다.

좀 혼동되어 사용된면이 있는데, 우선 나혼자라도 용어정리를 하고 난 주욱 이걸로 커뮤니케이션 하겠다.
  • individual callsets: 샘플별로 call된  variant들 셋트 (복수 결과들)
  • joint callset: 다수의 샘플이 한번에 call된 세트 (단수 결과)
  • joint variant discovery (join discovery): joint callset을 얻기 위한 일련의 과정(gvcf,,,reference model등을 활용)
  • join variant calling: 각 개별 샘플을 한번에 묶어서 calling하는 방법(sample #1 + ... + sample #N = one bam)
  • joint genotyping:  각 개별 샘플별로 calling 한 후 이를 묶는것 (sample #1, sample #N---->joint genotyping)
  • multi-sample calling: 다수의 샘플을 joint variant calling이나 joint genotyping등의 기법을 통해 분석하는 것
  • reference model: reference와 같은 site에 대해서 reference에 대한 confidence est 정보
     
참고
GATK HaplotypeCaller document
GATK blog - The GATK reference model pipeline for incremental joint discovery, in full detail
GATK Using the reference confidence calculation model & generating a GVCF document
GATK GenotypeFVCFs document
GATK presentations in DropBox

반응형