바이오인포매틱스

NGS 데이터에서의 Genotype and SNP calling

hongiiv 2014. 3. 22. 18:41
반응형
지금까지 수천 샘플에 대한 genotype/snp calling을 수행했음에도 2011년도에 발표된 "Genotype and SNP calling from next-generation sequencing data"라는 리뷰 논문을 이제서야 꺼내어 읽어본다. 이 논문에 대한 내용은 이미 "ideas should be in papers" 블로그에 소개되었으나 나름 다시 정리하는 차원에서 여기저기 살을 붙여서 작성했다. 간혹 이해가 가지 않는 부분은 참고하여 작성했다.

전체적인  genotype/snp calling 분석

Base calling

genotype/SNP calling에 있어서 가장 기초가 되는 것은 per-base quality score로 이는 일반적으로 NGS 장비의 기본 base calling software를 통해 제공된다. Illumina의 경우 RTA (Real Time Analysis software)를 사용하여 base calling을 수행하며, miss call error rate가 1%정도로 합성되는 시점에서 sync가 맞지 않아 signal이 맞지 않아서 발생하게 된다. 금방 언급했듯이 가장 기초적인 데이터로 Illumina의 base caller인 Bustard 대신 Ibis, BayesCall을 사용하게 되면 error rate가 5~30% 정도 향상된다고 한다. 

Summary of the available applications used for base-calling on the Illumina platofrm
Brief Bioinformatics(2011), Base-calling for next-generation sequencing platforms

 
Comparison of Base-Callers for the Illumina Platform 

Recalibration of per-base quality score

앞서서 base quality score는 genotype/snp call에서 가장 기본이라 했다. 여기서는 base calling 과정에서 생산된 raw Phread-scale quality score는 정확하지 않기 때문에 여러가지 covariate를 이용하여 재조정되어야 한다. 이러한 역활이 GATK의  BQSR (Base Quality Score Recalibration)이 담당하고 있다. 
  • Reproted quality score
  • Position within the read (machine cycle)
  • Preceding and current nucleotide (sequencing chemistry effect, dinucleotide context)
이러한 base quality를 조정하는 또 다른 방법은 Heng Li의 BAQ (Base alignment quality) 간단하게 INDEL 주변의 SNP들을 재조정하는 것이다. 여기서 잠깐 뒤에서 알아보겠지만, GATK나 Samtools에서의 genotype/snp calling은 기본적으로 bayesian model을 기반으로 하고 있지만, 그 결과는 slightly different하다는 것이다. 이러한 차이는 직접  Heng Li (아이디 lh3)가 언급한 내용을 정리하는 다음과 같다.

1. Alignemnt preprocess 부분

  • GATK는 low mapping quality를 가지는 read를 사용하지 않는 반면 samtools는 기본적으로 모든 read들을 사용한다.
  • Samtools는 base quality에 alignment quality를 덧씌워 사용하며, GATK가 현재 그렇게 하는지는 모르겠다.
  • GATK는 read pair에서 끝(ends)이 중복되는 경우 축소(collapses)하게 된다. 이게 맞는 것 같다. 하지만, 분석 정확성이 크게 향상되지 않을 것 같아 이것을 samtools에 구현하지 않았다.
  • samtools는 BAQ를 기본으로 사용하는 반면 GATK는 그렇지 않다. (앞서 언급했듯이 GATK에서도 구현되었음)

2. SNP genotype likelihood model 부분

  • GATK는 sequencing error가 독립적이라고 가정하는 반면 samtools는 두번째 에러 (second error)는 더 높은 기회를 가진다고 본다. 이러한 차이는 depth가 수백x가 되지 않는 이상 거의 영향을 미치지 않는다. 

3. Indel genotype likelihood model  부분

  • GATK와 samtool의 가장 큰 차이로 samtool은 BAQ를 기반으로 heuristic한 방법을 allele를 구분하는데 사용한다. GATK는 Dindel 모델을 기반으로 시간이 지나면서 많은 변화가 있으며 나(헨리)는 정확히 어떻게 동작하는지 모르겠다. 

4. Variant calling model 부분

  • samtools과 gatk 모델이 거의 동일하며, 최초 Rasmum Nielse과 Richard Durbin이 제안한 것으로 samtools 논문에 언급했으며 gatk 동일한 공식을 사용하여 구현했다. 추후 gatk는 multi-allelic 모델을 추가했으나 samtools는 포함하지 않았다.

5. 필터링 부분

  • 또다른 큰 차이로 samtools는 사람이 손으로 (hand-tuned)  필터를 사용하는 반면 gatk는 데이터로부터 filter를 학습한다.
  • 물론 model을 훈련하기 우히나 충분한 데이터와 non-polymorphic sites가 필요하지 않은 human variant calling에 있어서 gatk가 좀 더 편리하고 강력하다.
  • 게다가 효과적인 필터링을 위해 gatk는 samtools 보다  haplotype score와 같은 좀 더 많은 통계정보를 수집한다.
  • gatk가 진단의 목적 (diagnostic purposes)에 매우 유용하다.

6. GATK만의 특징

  • haplotype caller나 indel realignment 등
  • haplotype caller는 gatk v2의 가장 중요한 특징으로 longer indels에 대한 높은 민감도를 보여준다.

7 samtools만의 특징

  • genotype-free analysis와 physical phasing 등

8. technical 차이

  • gatk가 더 정교하고 더 복잡한 것은 사실이며 그럼에도 불구하고 samtools는 훨씬 간단한 프레임워크로 거대 데이터 셋 (huge data sets)에 동일하게 확장이 가능하다.

BAQ

초기 genotype calling 방법

초창기 시절의 genotype calling 방법을 알아 보기전에 우선 SNP calling과 genotype calling은 서로 다르다는 것을 염두에 두어야 할 것이다. snp calling은 polymorphisms이 있는 base를 찾는 것이며, genotype calling은 앞서 snp 부분에 대해서 개인의 genotype을 결정하는 것이다. 초창기 방법은 누가 들어도 충분히 구현이 가능한 방법으로 high-confidence bases (Q20) 남기는 fitering을 한 후 관찰된 allele를  고정된  cutoff (ex, heterozyous가 20-80%의 proportion을 갖는지?)를 사용하여 판단하는 것이다. 이러한 방법은 sequencing depth가 높은 (>20x) 경우 비교적 잘 작동하며 이를 software들로는 Roche's GSMapper, CLC Genomic WOrkbench, DNASTAR lasergene software들이 있다. 중요한 것이 바로 cutoff값인데 이는 경험적인 cutoffs를 사용하여 성능을 향상할 수가 있다.

확률론적 방법

앞서 살펴본 초장기 방법은 moderate 또는 low sequencing depth인 경우 heterozyouse gentotype에 대해 잘 작동하지 않으며 게다가 quality score에 기반한 filtering은 정보 손실로 이루어질 수 밖에 없다. 추가적으로 genotype을  추론하는데에 있어서 불확실성에 대한 어떠한 측정값을 제공하지 않는다는 것이다.
이러한 이유로 각각의 genotype의 quality score를 사용하여 사후확률 (posterior probability)를 제공하는 확률론적 방법이 개발되어 널리 사용되고 있다.

베이지안

경험을 통해 얻은 확률과 특정 선입견 (class)하에 들어온 관측치가 그 선입견 (class)에 들어맞을 확률을 통해 들어온 관측치가 특정 class일 확률을 계산하는 것이다. 예를 들어 아이가 침대에서 떨어질 수 있다는 특정 사건 (class)가 있다고 한다면, 침대 끝에 있다는 현재 상황에서 떨어질 수 있다는 class는 꽤 높은 사전확률 p(뒤로 떨어짐)을 갖는다. 이 class의 likelihood는 p(목이 뒤로 넘어감 | 뒤로 떨어짐)을 통해 목이 뒤로 넘어가는게 맞다라고 판단할 수 있다.  

바로 genotyp calling에 이 베이지안 방법을 적용하여 우리는 해당 position에 대해서 기존의 경험을 통해 해당 부분의 genotype에 대해 확률론적인 의미를 부여할 수 있게 된다. 

사전 확률 (Prior probabilty of genotype)

사건(class)에 대한 어떠한 정보도 없는 상태에서 뒤로 떨어지는 확률을 경험을 통해 알고 있다. 마찬가지로 genopype에서도 우리는 다양한 경험을 통해 정보를 알고 있다. 
  • genome wide SNP rate
  • SNP substitution type not equally likely
  • allele frequency
예를 들어, 다음과 같은 정보를 미리 알고 있다면 각 genotype에 대한 prior값을 계산할 수 있다.
heterozygous SNP rate: 0.001
homozygous SNP rate: 0.005
Reference allele: G
Transition/transverstion ratio: 2
라고 한다면 다음과 같은 prior probability을 갖게 된다. G/G가 나올 확률이 가장 높다는 것을 사전에 우리는 알게 되었다.


Use dbSNP prior probability

이외에도 다른 정보(reference sequence, SNP data-bases, available population sample) priors에 사용할 수 있다. 아래 그림은 dbSNP의 G/T site에 대한 prior probability이며(각 genotype에 대해 GG와 TT는 0.454, GT는 0.0909, 그외에는 10-4 이하)  G/G, T/T의 확률이 같은 것, 즉 G/T hetero site라는 것을 알 수 있다. 이러한 방법은 SOAPsnp가 사용하는 방법으로 MAQ또한 유사한 방법을 사용한다.

Use differnet polymorphism rate for different genomic regions
Consider different Ti/Tv rate for exonic regions


위의 경우에는 하나의 샘플 (single sample)에 대한 prior값을 구했지만, multiple sample을 고려하는 경우에는 prior값을 좀 더 improve할 수 있다. (추가해야함)이제 genotype G에 대한  사전확률 p(G)를 알았으니, genotype에 대한 likelihood p(X|G)값을 계산하기만 하면 우리는 Bayes' formula를 이용하여 p(G|X)를 계산할 수 있게 된다.

Genotype Likelihood

앞서 이야기 했듯이 read의 quality score를 이용해서 genotype likelihood 값을 계산한다. 한 개체에서 나온 특정 genomic site의 read i의 genotype G에 대한 확률 p(Xi|G)는 Xi의 quality score의 rescaling으로 계산이 가능하다. 그리고 genotype likelihood p(X|G)는 이 p(Xi|G)들의 product로 계산되어진다.

genotype likelihood를 계산하고 genotype prior(P(G))을 이용하여 Bayes' fomula로 각 genotype G에 대한 posterior probability (P(G|X))를 계산해서 최대의 확률을 갖는 genotype,G1을 그 개체의 genotype으로 할당한다. 이때 P(G1|X) 혹은 P(G1|X)와 P(G2|X) 사이의 비율(G2는 posterior probability가 두번째로 높은 genotype을 genotype confidence로 사용한다.

LD 정보를 포함한 calling

population genetics 관점에서 imputation은 SNP data set에서 missing data를 찾는 것으로 간단히 어떤 population에서 ATA, CGC라는 두개의 haplotype만 관찰된다고 하자. 어떤 개체의 샘플에서 첫번째 site에서 A나 C가 나오고 세번째 site에서 A와 C가 나온다면 이 샘플의 두번째 site를 모른다고 한다면 손쉽게 우리는 중간의 빈곳이 T또는 G라는 것을 알 수 있다. 바로 이것을 NGS 데이터에 접목하면 되는 것이다. 아래 그림은 다음의 3가지 방법에 대한 genotype call accuracy를 보여준다.
  • 각 개체에 대해 따로 genotype calling을 수행 (파란선)
  • LD기반의 분석을 사용하지 않고 모든 개체들을 jointly analysis (붉은선)
  • LD기반의 분석을 이용하여 모든 개체들을 jointly analysis (검정색)
multiple individual을 사용하는 것이 single sample에 비교해서 각각 80%에서 87%로 향상되는 것을 볼 수 있다. 약 40% 정도의 non-calls이 있다고 할때 LD 패턴을 고려하지 않는 분석과 LD를 고려한 분석이 비슷한 accuracy를 보이는 것으로 보아 이는 약 40% 정도의 missing data가 존재한다고 볼 수 있다. 이렇게 LD 패턴을 고려하는 경우 genotype calling에 있어서 많은 향상으로 보여주며 이는 hapmap이나 SeattleSNPs처럼 고품질의 LD정보가 있는 레퍼런스 데이터를 사용하는 경우 많은 도움을 받을 수 있다. (23andMe 데이터 또한 이러한 imputation을 통해 비어 있는 정보를 채워 넣을 수 있다) 하지만, reference panel에 없는 rare mutations에 대한 calling에서는 이러한 LD 정보는 성능 향상에 도움을 주지 못한다.

Filtering

각 site마다 사후확률이 정확히 계산되어졌다면야 필터링 등의 추가적이 데이터 조작은 필요없겠지만, 어디 그게 현실에 가능하단 말인가? 실제 데이터에서는 필터링을 통해 genotype/snp call이 향상된다. 1000 genomes project의 경우 이미 genotype이 알려진 HapMap 데이터와의 높은 차이로 인해 모든 sequencing batch를 제거해야 했다. 이러한 필터링은 genome-wide SNP genotyping이 이미 수행된 프로젝트에서나 가능하다. 인간 이외의 organisms이나 genotyping 데이터가 없는 경우라면 1000 genomes project 보다 높은 에러율이 나타날것으로 예상해야 할 것이다.

다른 형태의 필터링으로는 HWE로 인한 devaiations (일반적으로 low-quality scores, major/minor allele에 대한 sysmatic differnece in quality scores, LD aberrant, 너무 높은  read depth, strand bias 등)가 genotype과 snp calling의 정확도를 높이는데 활용될 수 있다.  적절한 필터는 시퀀싱 프로토콜과 upstream 분석에 의존한다. 예를 들어 strand bias (plus와 minus strand가 불균형적인 부분)가 존재하는 부분에서는 에러가 나기 쉬운 부분으로 이 부분은 필터 아웃 되어야 한다. exome 캡처를 사용한는 captured 시퀀스인 경우 이러한 bias는 capture array의 artefact로 인할 가능성이 크다.

Incorporating genotype uncertainty in analyses of NGS data

NGS 데이터 분석에서의 genotype에 대한 불확실성을 아예 안고 살아라, 즉 분석에 이러한 불확실성을 포함해서 분석을 하라는 것입니다. 다른 application은 서로 다른 genotype calling 전략을 사용하라는 것이다. 여기서는 association-mapping studies와 population-genetic studies에서 이러한 불확실성에 대한 이야기를 하고 있다. 

Recommendations

NGS 데이터 분석은 빠르게 변화하는 분야로 분석을 위한 새로운 통계 방법들이 계속 개발되고 있다. 더구나 1000 genomes project로 부터 많은 툴들이 개발되었지만 아직 publish되지 않은 것들도 많다. 이러한 시점에서 몇가지 사항에 대해 권장한다.

1. Base calling과 quality scores 계산은 테스트/벤치마크된 방법을 사용해라.
  • per-base quality scores의 recalibration은  GATK나 SOAPsnp를 추천
  • reference genome/short-read를 align하는 경우 Novoalign이나 Stampy와 같은 sensitive aligner를 사용해라.
  • 후에 bwa와 같이 효율적인 aligner와의 hybrid 사용 추천
2. SNP calling은 모든 샘플을 동시에 할 수 있는 것을 사용한다.
  • snp calling에 있어서 likelihood ratio나 베이지안을 사용하기 때문에 데이터로 부터 allele freqeuncy에 사전 분포를 사용하기에 모든 샘플을 동시에 분석하는 것이 좋다. 
3. genotype calling 또한 다수의 individual을 합쳐서 한다.

4. 가능하다면 LD-기반의 방법을 통해  genotype과 snp calling의 accuracy를 향상시킬 수 있다.

몇가지 추가적인 고정을 통해 성능을 향상 시킬 수 있는데, 가령 local realignments, 다수의 snp/genotype calling 알고리즘과 quality score기반의 post hoc 필터링을 사용할 수 있다. 추가로 아래는 1000 genomes project에서 variant/genotypes call의 성능 향상을 위해 사용한 방법으로 이미 위에서 설명한 방법들이 다수 포함되어 있다.
  • INDEL realignment
  • Per-base alignment quality (BAQ) adjustment
  • Robust consensus SNP selection strategy
    • variant quality score recalibration (VSQR)
    • support vector machine (SVM)
  • imporved genotype likelihood calculation
    • BAM-specific Binomial Mixture Model (BBMM)
    • leveraging off-target exome reads 
NGS 데이터에서 genotype/snp calling은 불확실성에 대한 확률적인 방법을 제공하는 것으로 allele의 count에 따라 간단한 방법에서 보다 정교한 많은 방법들이 존재한다. 확률적인 방법은 genotype likelihoods의 계산에 의존하므로 더 많은 연구가 필요하다. 또한 LD 기반의 방법의 개선을 포함한 다운 스트림 분석에서 genotype 확률을 통합하는 방법이 필요하다. NGS는 어찌되었든 향후 몇년동안 medical genetics 연구의 중심이 될 것이며 미래를 위해 충분히 가치가 있다.

반응형