[GATK] Base Quality Score Recalibration (BQSR)

[GATK] Base Quality Score Recalibration (BQSR) Start.

BioinformaticsAndMe

GATK 파이프라인의 데이터프로세싱 과정인

Base Quality Score Recalibration (염기서열점수 재보정)

을 이해해보자.

1. 왜 Base Quality Score Recalibration이 필요할까?

Base quality score는 시퀀싱 머신에서 각 base마다 발생하는 error의 추정치이다.

예를 들어, 특정 base가 Q20이라 하면, Phred-score 개념으로 99% 정확한 base라는 의미이다.

다시 말해 염기서열 100개가 있으면 1개 정도는 틀릴 수 있음을 말한다.

100개에서 1개 틀리는 거면 괜찮아 보이지 않은가?

그러나 30억 Human genome을 생각해보자. 보통 WGS는 30X 정도 되므로,

30억 * 30X = 900억개의 base call이 발생할 것이다.

1%의 error로 봤을 때, 약 9억개 base가 error를 가지고 발생한 call이라는 것이다.

많은 variant calling 알고리즘이 각 base에 할당된 quality score에 크게 의존한다.

뭐 당연할 것이다. 스코어가 높을수록 검출된 variant가 우연이 아닌 진짜일 확률에 가까워지니 말이다.

그런데 운이없게도 9억개의 error base를 가지고 variant를 뽑아낸다면, 우리는 잘못된 결론을 도출할 가능성이 높다.

그래서! 각각의 base score를 다시한번 Recalibration(재측정) 하여, 좀 더 정확한 base quality score를 부여하는 과정이

'Base Quality Score Recalibration' 이다.

아래는 base recalibration이 GATK 알고리즘의 어떤 위치에서 진행되는지 보여주는 흐름도이다.

#원래 GATK에서 Base reacalibration은 'Indel realignment' 과정 후에 진행하였는데, GATK4로 넘어가면서 'Indel realignment'가 없어졌다.

#요즘 데이터 퀄리티가 realignment할 정도로 나쁘지 않고, 파이프라인의 다른 과정에서 realignment 보정 과정을 만회할 수 있다는 것 같다.

#Indel realignment가 시간이 오래걸리는데, 굳이 할 필요가 없으니 GATK4에서 그냥 뺀 것으로 보인다.

2. BQSR(Base Quality Score Recalibration)은 어떤 과정으로 진행되나?

먼저 BQSR을 돌리기 위해 필요한 것은 Known Single Nucleic Polymorphisms (SNPs) 이다 (dbSNP).

그 이유는 BQSR algorithm에서 dbSNP에 매칭되지 않는 base는 '에러'라 가정하고 진행되기 떄문이다.

A) Finding errors

# 아래 예제를 살펴보자 (어떤 염색체 위에 있고 0~9 base까지 길이 10인 리드가 존재한다)

BQSR에서 에러라고 여겨지는 포지션은 3번과 7번이다. 3번과 7번은

1)read base가 reference와 다르고, 2)dbSNP에도 없기 때문이다.

B) Aggregating the reported phred score

위 phred score (10  11  11  20  22  22  30  20  20  10)는 아래처럼 확률로 변환시켜, 10개의 평균 확률을 구할 수 있다.

(0.1 + 0.079 + 0.079 + 0.01 + 0.006 + 0.006 + 0.001 + 0.01 + 0.01 + 0.1)/10 = 0.0401

이것을 다시 phred score로 변환하면,  phred score = -10 * log10(0.0401) ~= 14

따라서, 예제 리드의 reported phred score는 약 14이다.

C) Calculating the empirical phred score

이번에는 empirical(경험적) phred score를 구해보자.

10개의 염기 중 2개를 에러로 가정했으므로, 시퀀스는 2/10 = 0.2 정도의 에러 확률을 갖을 것이고,

(경험적이라는 의미는 실제 측정된 phred score를 이용하는 것이 아닌, 에러인지(true) 아닌지(false)를 true개수/전체개수의 빈도확률 형태로 나타낸 것)

phred score로 변환하면 -10 log10(.02) ~= 7.

따라서, empirical phred score는 약 7이다.

empirical phred score가 7 정도의 에러가 있다면, 원래 관찰된 phred score는 +7정도의 에러가 더해진 값이라 생각하자.

따라서, 각 value에 -7씩 감해주어 recalibration을 진행하자.

10-7=3, 11-7=4....

#GATK에서 위 과정은

- BaseRecalibrator로 recalibration 모델을 만들고

- ApplyBQSR로 score를 재조정하는

두 가지 command로 이루어진다.

#https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_bqsr_BaseRecalibrator.php

#https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_bqsr_ApplyBQSR.php

#아래 그림은 recalibration 후, 측정된 quality와 경험적으로 확인된 quality의 일치도가 높아진 모습이다.

#다시 말해, read quality의 오류가 적절하게 보정되었다는 의미.

마무리하며..

시퀀싱데이터에서 Recalibration은 quality 보정을 위해 반드시 필요한 작업이다.

특히, dbSNP과 같이 여러 Known 데이터베이스 갖춘 Human에서는 error model을 다른종보다 쉽게 제작할 수 있기에,

정확한 downstream variant calling을 위해선, recalibration 과정이 error model에 기반하여 확실히 수행되어야 한다.

recalibration 모델에 사용되는 feature들은 아래 4가지 정도가 되겠는데, 자세한 사항은 아래 사이트를 참고하자.

• read group the read belongs to
• quality score reported by the machine
• machine cycle producing this base (Nth cycle = Nth base from the start of the read)
• current base + previous base (dinucleotide)

#참고 사이트

1) https://software.broadinstitute.org/gatk/documentation/article?id=11081

2) http://zenfractal.com/2014/01/25/bqsr/

[GATK] Base Quality Score Recalibration (BQSR) End.

BioinformaticsAndMe