BioinformaticsAndMe

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

BioinformaticsAndMe

을 이해해보자.

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 알고리즘의 어떤 위치에서 진행되는지 보여주는 흐름도이다.

그 이유는 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

- 아래 그림 오른쪽 표를 이용해 베이지안 모델에 적용하여, 가장 높은 evidence를 가지는 variation을 선택하여 변이 동정.

[GATK] HaplotypeCaller 알고리즘 End.

BioinformaticsAndMe

12. PCA plot 에 샘플 좌표 표시하기

-PC1과 PC2 각각에 투영된 좌표를 연결하면, PCA plot 에서 샘플들의 위치를 알 수 있다.

-위 그림은 'Sample2' 좌표가 확인되었다.

-'Sample1'의 좌표도 확인할 수 있다.

-같은 방식으로 진행했을 때, PCA plot 에서 모든 샘플들의 좌표를 확인하게 된다.

13. 주성분의 Variation

-SS 값은 origin (0,0)에서 샘플들의 d 제곱합 이었다 (주성분분석 2 칼럼 참조).

-SS 값에 자유도 (n-1)로 나누게되면, 주성분에 대한 Variation 값이 나오게 된다.

-위 설명을 분산의 공식으로 이해하면 되겠다 (Variation 값이 클수록 샘플들의 분포가 넓게 퍼져있다는 의미다).

14. 주성분의 설명력 이해하기

-위 그림의 예처럼, PC1의 Variation 값을 15, PC2의 Variation 값을 3 이라 해보자.

-두 개 주성분의 Total variation은 15+3 = 18 이다.

-PC1의 Variation은 Total variation의 83% (15/18)을 설명한다.

-같은 방식으로 PC2의 Variation은 Total variation의 17% (3/18)을 설명한다.

-왼쪽 Bar plot이 PCA에서 주성분들이 샘플들을 어느정도 설명할 수 있느냐를 보여주는 것이다.

-설명력이란 표현을 사용하며, PC1은 변량의 83%이상 설명할 수 있고 PC2는 변량의 17%를 설명할 수 있다.

-PC1이 PC2보다 샘플들의 유사성(비유사성)을 보는데 상대적으로 더 정확하다고 할 수 있겠다.

-쉽게 설명하면, 예로든 PCA plot 에서 위아래의 차이보다 좌우의 차이가 더 큰 의미를 갖는다고 말하겠다.

마무리하며...

후 이제 끝났다....

PCA 알고리즘을 설명하는데 오래 걸렸다.

https://www.youtube.com/watch?time_continue=285&v=FgakZw6K1QQ

사실, 칼럼의 출처인 위 동영상과 함께본다면 더 쉽게 이해될 수 있다.

개인적으로 유전자 발현분석에서 PCA를 엄청나게 많이 사용했었는데, 개념과 원리정도만 이해하고

만들어지는 알고리즘에 대해서는 무심했던 것 같다.

다음 칼럼에서는 Multi-Dimensional Scaling (MDS, 다차원척도법)과의 공통점 및 차이점에 대해 살펴보겠다.  

[PCA] 주성분분석 3 End.

BioinformaticsAndMe

8. Best line 을 찾기 위한 정량적 접근 방법

-샘플들의 좌표(초록색점)를 설명하는 Best line을 찾는 과정은 두가지로 설명된다.

-먼저 위에 그림처럼 임의의 빨간선을 그엇을 때, 샘플들의 최단거리를 '초록색 X'로 표시한다.

-모든 샘플의 최단거리를 표시하고, 모든 거리의 합이 가장 최소가 되는 Best line을 찾는게 첫번째 방법이다.

-'Minimization method' 라 칭하겠다. 

-다음은 Origin (0, 0) 에서 '초록색 X'로 표시까지의 거리를 최대로하는 Best line을 찾는 것이 두번째 방법이다.

-'Maximization method' 라 칭하겠다.

-사실, Minimization이나 Maximization 모두 계산해보면 같은 의미를 가지게 된다.

-Best line에서 샘플 좌표들의 거리의 합이 최단일수록, 원점에서 각 초록색 X의 거리의 합은 최대가 된다.

-수학적으로 우리가 잘 알고있는 피타고라스 정리를 적용하면 이해가 빨라진다. 아래를 살펴보자.

9. 피타고라스 정리 적용

-'a' : 원점에서 샘플 좌표까지의 거리 (일정함, 고정된 값)

-'b' : Best line 과 샘플 좌표의 최단 거리 (Minimization 값)

-'c' : 원점에서 샘플 좌표의 최단 거리에 있는 Best line까지의 거리 (Maximization 값)

-'a'의 값이 변하지 않고, 피타고라스 정리에 따라 b값이 적어질수록 c값이 커지고, b값이 커질수록 c값이 작아진다.

-정리하자면,  PCA는 'b'를 최소로하는 혹은 'c'를 최대로하는 Best line (주성분)을 찾는 과정이라 볼 수 있겠다.

10. Sum of Squared distances (SS) 찾기

-Best line을 찾는 과정은 주로 Maximization을 활용한다.

-위에 'd1' 과 같이 최대거리를 통해 주성분을 만드는게 실제 계산과정에서 용이하다.

-Minimization은 데이터의 최단거리를 이용해 주성분을 만든다는 해석의 측면으로 접근한다.

-우리가 가진 샘플 6개의 d 값 (d1, d2, d3, d4, d5, d6)를 구할 수 있고 제곱하여 합한다.

-제곱하는 이유는 d 값은 상대적이므로 음수의 값을 가질 수 있기에 제곱한다.

-제곱의 합은 'sum of squared distances', 줄여서 'SS'라 부른다.

-우리는 저 SS가 최대인 Best line (주성분)을 찾으면 끝난다.

11. Principal Component 1 (PC1)  &  Principal Component 2 (PC2)

-우리가 찾아낸 Best line이 위에서 언급했던대로 '주성분 1 (Principal Component 1)'이 된다.

-위 그림의 화살표는 Gene1과 Gene2의 분포정도와 PC1 의 관계를 설명한다.

-'우리가 PC1을 만들었을때, 샘플들의 분포는 Gene1보다 Gene2의 스케일을 따르는 경향이 있다' 라 이해하면 되겠다.

-우리의 주성분 PC1 은 Gene1의 값에 크게 영향받고 있으며,

-다시 해석하여, 샘플들은 Gene1 발현에 큰 편차를 보인다.

-'주성분 2 (Principal Component 2)' 는 PC1과 직교하는 라인이라 생각하면 되겠다.

-PC2의 자세한 설명은 동영상을 참고하시면 된다.

-앞서 과정에서 만든 PC1과 그에 수직하는 PC2을 rotation하면 위와 같은 그림을 볼 수 있다.

-우리가 그동안 봐왔던 PCA의 그림의 두 축이 생성된 순간이다.

마무리하며..

PCA 알고리즘 마지막 파트에서는 PC 축에 늘 함께 붙어있는 ('%' ,설명력) 을 설명하는 시간을 갖겠다.

[PCA] 주성분분석 2 End.

BioinformaticsAndMe

사실 PCA는 워낙 유명하기에, 구글이나 여러 블로그에서 매우 자세하게 설명하고 있다.

본인 역시 그곳에서 공부하였다.

그런데 한 3개월 전에....   PCA 원리를 아주 쉽고 자세하게 설명해주는 20분짜리 유투브 영상을 발견하였다!!

주성분분석을 한번도 접해보지 못했다면 아래 영상을 20분을 들여 보는 것을 매우 추천한다.

물론 나역시 영상을 보고 이해하고 덧붙일 내용을 아래아래에 정리하겠다.

1. 예제 데이터시트

-예제로 사용할 데이터 시트가 되겠다. 만든이가 genetics하시는 분이라 마우스와 유전자 발현 테이블로 설명한다.

-Mouse는 PCA plot 에 점으로 찍힐 Sample들이 되겠고, Gene은 축으로 사용될 Variable이 될 것이다.

2. 하나의 Variable부터 설명 시작 ! (1차원)

-먼저, 마우스 실험에서 한개의 유전자발현만 확인했다면, 6개의 마우스 각각이 하나의 발현값을 갖고 있을 것이다 (표 참고).

-이를 오른쪽 아래 수평선에 값으로 나열할 수 있겠고, 'Gene1'의 발현값이 유사한 마우스들을 볼 수 있다.

-마우스 6;5;4 가 상대적으로 낮은 발현을 보이고, 마우스 3;1;2가 상대적으로 높은 발현을 보인다.

3. 유전자 하나가 더해져 두개의 Variable ! (2차원)

-Variable이 두개가 됐으니, 축이 하나 더필요하겠다. 'Gene1'을 X축, 'Gene2'를 Y축으로 뒀다.

-마우스샘플들은 좌표 (X,Y)에 맞게 오른쪽 plot 처럼 찍혔다.

-마우스 3;1;2 가 좌표상으로 봤을 때 두 유전자의 비슷한 발현패턴을 보였다 (5;6;4도 역시 지들끼리 비슷한 발현패턴).

4. 이제 3개의 Variable ! (3차원)

-'Gene3'가 추가되어 Z축 까지 생겼다면 우리는 오른쪽 plot과 같이 3차원 형태로 놓인 샘플들을 볼 수 있다.

-예를 들어, 1번마우스가 (X, Y, Z) = (10, 6, 12) 좌표를 가지므로 저 멀리 오른쪽 구석즈음에 표시될 것이다.

5. 4차원을 그려보자..

-4차원은.. 인터스텔라로

-Array같은 마우스 실험에서 유전자 한개, 두개만 보는 일은 없다. 수천개 혹 2만개 유전자의 발현을 확인해야할텐데...

-모든 값을 표현하여 온전히 그릴수도 없거니와 샘플의 유사성도 보기 어렵다. 따라서,

-우리는 필연적으로 모든 유전자들의 발현값 분포를 설명할 수 있는 두개의 대표축을 뽑아 시각화하고 싶어한다.

-그 두개의 축이 PC1 (Principal Component 1, 첫번째주성분)과 PC2 (Principal Component 2, 두번째주성분) 이다.

-PC1 이 전체 발현데이터 분포의 91%를 설명할 수 있고, PC2 가 4% 를 설명할 수 있다 (계산법은 뒤에서..)

-글 초반부에 정의한대로, PCA는 고차원 데이터를 차원축소하여 높은 설명력으로 샘플들의 유사성을 확인하는 기법이 되겠다.

6. 좌표 (0,0) 으로 센터맞추기

-'Gene1'의 발현값 평균을 X축에 찍어보자 (빨간 X 마크).

-같은 방법으로 'Gene2'의 발현값 평균을 Y축에 찍어보자 (빨간 X 마크).

-이제 우리는 위 그림에서 두 평균의 좌표 (Xmean, Ymean) 을 확인할 수 있다 (파란색 X 마크).

-우리는 파란색 X 마크 (평균의 좌표)를 Origin (0, 0)으로 아래 그림처럼 옮길 수 있다.

-주목할 사실은 좌표를 (0,0)으로 옮겼어도 샘플들 사이의 상대적 거리는 변하지 않았다는 것이다.

-'Gene1'에서 가장 발현이 높은 샘플은 여전히 제일 오른쪽에, 'Gene2'에서 가장 발현이 높은 샘플은 여전히 제일 위에 존재한다.

7. 샘플 좌표에 가장 적합한 라인 그리기

-샘플 좌표에 잘 맞는 선을 그렸다. 손으로 찍 그리기는 쉽겠지만,,

-우리는 이렇게 Fitness가 높은 선을 정량적인 접근으로 찾아낼 것이다.

마무리하며..

내용이 길어져 다음 포스트에 이어서 작성하겠다.

[PCA] 주성분분석 1 End.

BioinformaticsAndMe

마무리하며..

Imputation 과정에서 가장 신경써야할 부분 중 하나는 '어떤 Reference를 사용할 것인가?' 이다.

일반적으로 HAPMAP과 1000G이 사용되고 있으며, 그나마 1000G에 중국과 일본의 population이 들어가 있어,

한국인 GWAS 분석에서 적합한 reference라 볼 수 있겠다.

그러나,, 한국인 GWAS 분석에서 'Korean reference'를 사용해야하는 건 아주 당연하다.

질병관리본부에서 'KRGDB(Korean Reference Genome) project'로 한국인으로 구성된 ref가 있다.

(http://152.99.75.168/KRGDB/menuPages/intro.jsp) - raw 데이터 사용하려면 따로 신청해야 하는듯..

끝.

[GWAS] Imputation End.

BioinformaticsAndMe

BWT 알고리즘 Start.

BioinformaticsAndMe

Alignment라는 과정은 Human 분석에서 특히 빡세다..

시퀀싱을 통해서 얻은 길지 않은 리드로 30억 염기서열과 매칭하는 일은 쉽지 않아보인다.

(음. Pacbio나, 요즘핫해보이는? Nanopore는 리드 길이가 길다)

NGS sequence를 alignment하는 tool로  bowtie 와 BWA 가 널리 알려져 있다. 이 tool들이 주력으로 사용하는 알고리즘은 Burrows-Wheeler Transform이다. 30억 염기서열 중에 ACGTACGTACGT를 align하기 위해서 모든 sequence를 처음부터 검색하면 시간이 매우 오래 걸린다. 이를 해결하기 위해서 사용한 알고리즘이 BWT이다.

BWT(Burrows-Wheeler Transform)

마이클 버로우즈와 데이비드 휠러가 1994년에 발표한 블록 정렬 알고리즘으로 변환 결과에 Index 정보가 포함되어 있어, 다른 정보가 없더라도 변환된 문자열의 경우 유사한 문자열들끼리 뭉쳐진 형태로 나타나는 경우가 많아 압축을 위한 전처리 알고리즘으로 사용된다.

#BWT 과정을 이해해보자.

#그래서 BWT 특징이 뭐야

Ⅰ. 같은 single nucleotide끼리 붙어있는 경우의 수 증가

Ⅱ. 문자의 순서가 유지

#BWT를 이용한 문자열 매칭

-First column은 last column의 뒤에 있었던 문자.

-$는 원래 문자의 맨 마지막에 붙어 있었음. 따라서 A가 맨 마지막 글자임을 알 수 있음.

-Last column과 First column의 문자의 순서는 동일함.

-A0의 앞에는 N0임을 알 수 있음.

-L은 F앞의 문자임.

#Bowtie (FM-Index)

Bowtie 알고리즘인 FM-Index는 BWT 행렬의 Last First mapping 법칙을 이용하여 문자열 Index를 만든 방법으로 문자열에 대해 BWT를 수행함으로써 얻어짐.

#BWA (BWT + Suffix array)

BWA는 BWT와 접미사 배열(Suffix array)을 이용하여 정렬을 수행하는 알고리즘.

BWT 알고리즘 End.

BioinformaticsAndMe