4.1. GATK BQSR 실행 (Base Quality Score Recalibration) [optional,but recommended]

GATK 본격 시작

https://software.broadinstitute.org/gatk/best-practices/

10번은 읽어야 NGS variant calling을 안다고 할 수 있지 않겠는가?

Illumina 기기 이전부터, Sanger sequencing을 할 때에도 Phred-score 로 염기서열 읽은 것에 대한 수치를 기록해 왔다.

우리는 이러한 Phred-score 기반으로 변이가 맞는 가 아닌가를 베이지안 추론을 진행하게 된다. 그러기 때문에 Illumina 기기 특유의 bias를 줄이지 않는다면 계통적 오류(systemic error)에 의한 위양성(false positive)를 얻게 된다. 

물론 위음성(false negative)도 생긴다.

* 잠깐 프레드 스코어는 해당 결과를 얼마나 믿을 수 있는가에 대한 수치이다.

Phred quality scores ${\displaystyle Q}$ are defined as a property which is logarithmically related to the base-calling error probabilities ${\displaystyle P}$.^[2]

${\displaystyle Q=-10\ \log _{10}P}$

or

${\displaystyle P=10^{\frac {-Q}{10}}}$

For example, if Phred assigns a quality score of 30 to a base, the chances that this base is called incorrectly are 1 in 1000.

Phred quality scores are logarithmically linked to error probabilities

Phred Quality Score	Probability of incorrect base call	Base call accuracy
10	1 in 10	90%
20	1 in 100	99%
30	1 in 1000	99.9%
40	1 in 10,000	99.99%
50	1 in 100,000	99.999%
60	1 in 1,000,000	99.9999%

이러한 사유로, Phred score를 보정하기 위해서 BQSR이라는 과정을 거치게 된다.

그런데, 이것을 하려면 해당 참조유전체에 대하여 분석된 VCF가 있어야 한다.

처음 분석이라 VCF가 없으면 어떻게 하나????

이거 안하고 바로 다음 과정을 해야 한다.

그래서 optional한 과정이지만, 재대로 분석하려면 필요하기 때문에 recommended라고 하였다.

GATK 홈페이지에서는 VCF 없으면 우선 없이 VCF 만들고 그것을 input으로 BQSR 돌리고 다시 VCF 만들라고 한다.!!!!!! (이렇게 시간 보내고 컴퓨터 쓰는 걸 설득시키는 과정이 필요하다.)

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

표 만들고, BAM 보정한다. 

gatk --java-options "-XX:+UseParallelGC -XX:ParallelGCThreads=6 -Xmx64g" BaseRecalibrator -R [reference genome fasta] -I [input bam] -O [table file] --known-sites [해당 reference로 분석된 VCF]
gatk --java-options "-XX:+UseParallelGC -XX:ParallelGCThreads=6 -Xmx64g" ApplyBQSR -I [input bam] -O [output bam] -bqsr [table file]




하고 안하고의 차이가 다소 있다.

하자. 그래야 best practice다.