SNP calling using Naive Bayesian classifier

Shotgun sequencing을 이용한 Single Nucleotide Polymorphism (SNP) calling 시 Bayes' Theorem이 자주 이용되는데, 우리 연구실의 Thomas Bleazard가 정리해 놓은 문서를 기초로 간략하게 방법론을 정리해 보고자 한다.

염색체의 기본단위는 A/C/G/T로 이루어진 nucleotide이고 세포는 상동염색체로 구성이 되어 있으므로 어떤 위치(locus)에서 가능한 genotype은 총 10가지가 가능하다. 수식으로 하면,

For given locus i, true genotype,

그리고 그 locus i에 align된 read를 k개라고 가정했을 때

그리고 q[ij]는 locus i에 align되는 j-th read의 quality score를 의미한다.

여기서 우리는 locus i에서 D[i]를 관찰했을 때, 각각의 T[i]에 대한 사후확률(posterior probability)를 계산하려 한다. 이는 Bayes' Theorem에 의해 다음을 만족한다.

여기서 P(T[i])는 사전확률(prior probability)이다.

K개의 read가 각각 독립이라고 가정하면 위의 식에서 P(D[i]|T[i])는

이다. 그리고 위 수식에서 마지막 probability는 training 용으로 이용할 data에서 직접 계산한 frequency matrix를 이용한다. 즉 o[ik], q[ik], T[i] triplet으로 이루어진 frequency matrix를 training data에서 산출하고 이를 이용한다. 그래서 각각의 T[i]에 대한 사후확률을 계산한 후 최대 확률을 보이는 genotype을 선택(classification) 한다.

Principal Component Analysis (PCA)

컴퓨터 전공이면서도 대학교 3, 4학년 때 미국으로 교환학생 가면서 수강하지 못했던 Linear Algebra (선형대수), 나중에 생물정보학 공부에도 발목을 잡을 줄을 몰랐다. 좋은 책 하나 골라서 차근차근 공부를 좀 해야할 듯.

Eigenvector, eigenvalue와 관련이 되어있는 Principal Component Analysis (PCA)에 관한 아주 쉬우면서도 잘 정리되어 있는 문서 발견 ! ("A tutorial on Principal Component Analysis")

튜토리얼에 따르면 실제 data의 covariance matrix으로부터 구해지는 eigenvector 들은 그 data를 잘 나타내는 (data의 pattern을 잘 드러내는) 축이 된다. 이 축으로 기존 data를 transformation 하는 과정이 PCA analysis의 정의.

실제 data의 covariance matrix가 n x n 일 때 n개의 eigenvector가 존재하는데 이를 eigenvalue의 크기로 정렬하면 제일 큰 eigenvalue에 해당하는 eigenvector 축으로 data를 transformation 하면 그 값이 principla component 1 (PC1)이 된다. 그 다음 큰 eigenvalue에 해당하는 eigenvector를 이용하면 PC2... 이하는 같다. 여기서 n개의 PC를 모두 쓰지 않고 일부분만 사용했을 때 dimension reduction 효과를 얻을 수 있다.

여기서 내가 들었던 의문은 그럼 sample을 clustering 하기 전에 원래 data를 사용하지 않고 PCA를 통해 data를 변형하고 일부 PC만을 사용해서 clustering를 하면 더 효과적인가? 였다. 이에 대해서는 그렇지 않다고 이미 논문이 존재 ("An empirical study of Principal Component Analysis for clustering gene expression data, Bioinformatics 17(9):763-74, 2001")

lines

lines(c(x1, x2), c(y1, y2), lwd="", col="")

점 (x1, y1)에서 시작하여 점 (x2, y2)에서 끝나는 선을 그린다. lwd(line width), col(color)

물론 이 함수 이용 전 plot을 이용하여 그래프가 존재하는 상태여야 한다는 것

예제)

plot(NULL, NULL, xlim=c(0, 1000), ylim=c(0,1000))

lines(c(100, 200), c(300, 400), lwd=3, col="red")

 



Illumina sequencing library preparation, amplification and sequencing

예전부터 개념만 대충 알았지 sequencing에 이용되는 oligonucleotide sequences에 대한 정확한 원리는 몰랐었는데 이번 기회에 공부를 좀 해 보았음.
 

Genomic DNA oligonucleotide sequences Adapters 1
5' P-GATCGGAAGAGCTCGTATGCCGTCTTCTGCTTG
5' ACACTCTTTCCCTACACGACGCTCTTCCGATCT

PCR Primers 1
5' AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
5' CAAGCAGAAGACGGCATACGAGCTCTTCCGATCT

Genomic DNA Sequencing Primer
5' ACACTCTTTCCCTACACGACGCTCTTCCGATCT



다음과 같이 주어졌을 때 도대체 library는 어떻게 만들어지는건가...



Adapter는 fragmentation 된 DNA의 양쪽 flanking한 지역에 붙인다. 당연히 DNA가 double strand 상태이니 double strand 형태의 adapter를 붙여야 할 터. 위에 언급된 두 개의 adapter가 double strand를 이루는 것이 아니라 각각의 adapter의 complement 한 sequence를 이용해 double strand를 만든다. 여기서 첫번 째 adapter 5'의 'P'는 무엇인가? 이것은 gap으로서 나중에 'A'가 overhang 된 DNA와 연결시키기 위해 존재한다. 우리가 관심있는 DNA에 adapter를 붙인 형태는 아래와 같다.

 

5' -------------------- -----ACACTCTTTCCCTAC ACGACGCTCTTCCGATCT (DNA+) AGATCGGAAGAGCTCGTATG CCGTCTTCTGCTTG 3'

3' -------------------- -----TGTGAGAAAGGGATG TGCTGCGAGAAGGCTAGA (DNA-) TCTAGCCTTCTCGAGCATAC GGCAGAAGACGAAC 5'

여기서 DNA+는 관심있는 DNA의 + strand, DNA-는 관심있는 DNA의 - strand을 의미. 그 다음 위의 adapter-ligated double strand DNA를 denaturation을 하여 single strand 로 나눈다. 그리고 나서 각각의 single strand DNA가 어떤식으로 amplification 되고 sequencing 되는지 따라가 보기로 하자. 여기서는 (DNA+)의 single strand DNA만 따라가 보기로 한다(DNA- 도 해보면 유사).

Illumina sequencer flowcell에는 PCR primer들이 covalent linking 되어 있다(flow cell surface에 primer의 5' 연결). 자 우리는 아래의 sequence로 시작해 본다.

5' -------------------- -----ACACTCTTTCCCTAC ACGACGCTCTTCCGATCT (DNA+) AGATCGGAAGAGCTCGTATG CCGTCTTCTGCTTG 3'

이 sequence는 두번 째 PCR primer와 binding 한다 (빨간색: 두번 째 PCR primer, 5' 쪽이 flowcell에 고정).
 

5' -------------------- -----ACACTCTTTCCCTAC ACGACGCTCTTCCGATCT (DNA+) AGATCGGAAGAGCTCGTATG CCGTCTTCTGCTTG 3'

3' -------------------- -------------------- ------------------ (----) TCTAGCCTTCTCGAGCATAC GGCAGAAGACGAAC 5'

그리고 3' 쪽으로 extension 되어 다음과 같은 double strand가 된다.
 

5' -------------------- -----ACACTCTTTCCCTAC ACGACGCTCTTCCGATCT (DNA+) AGATCGGAAGAGCTCGTATG CCGTCTTCTGCTTG 3'

3' -------------------- -----TGTGAGAAAGGGATG TGCTGCGAGAAGGCTAGA (DNA-) TCTAGCCTTCTCGAGCATAC GGCAGAAGACGAAC 5'

 

여기서 2nd PCR primer에서 만들어진 strand는 flow cell에 고정되어 있으므로 wash 과정을 거치면 flow cell에 covalent linking 된 다음과 같은 single strand만 남게 된다.
 

3' -------------------- -----TGTGAGAAAGGGATG TGCTGCGAGAAGGCTAGA (DNA-) TCTAGCCTTCTCGAGCATAC GGCAGAAGACGAAC 5'

 

자, 이제 소위 얘기하는 bridge amplification 과정을 살펴 보자.

지금까지의 single strand는 구부러져(bridging) 3' 쪽 sequence와 상보적인 PCR primer를 찾아 결합한다(파란색 sequence는 1st PCR primer, 5' 쪽이 flow cell과 고정).
 

5' AATGATACGGCGACCACCGA GATCTACACTCTTTCCCTAC ACGACGCTCTTCCGATCT (----) -------------------- -------------- 3'

3' -------------------- -----TGTGAGAAAGGGATG TGCTGCGAGAAGGCTAGA (DNA-) TCTAGCCTTCTCGAGCATAC GGCAGAAGACGAAC 5'

 

그리고 3' 쪽으로 extension.

5' AATGATACGGCGACCACCGA GATCTACACTCTTTCCCTAC ACGACGCTCTTCCGATCT (DNA+) AGATCGGAAGAGCTCGTATG CCGTCTTCTGCTTG 3'

3' -------------------- -----TGTGAGAAAGGGATG TGCTGCGAGAAGGCTAGA (DNA-) TCTAGCCTTCTCGAGCATAC GGCAGAAGACGAAC 5'

 

그리고 denaturation을 하면 이제 5' 쪽이 flow cell에 고정된 두 개의 single strand가 된다. 그리고 각각의 strand가 같은 방법으로 반복하여 amplification 되어 cluster를 만든다.

그럼 sequencing primer는 어떻게 쓰여지는가? Cluster 중 일관되게 DNA+ sequence를 알기 위해 1st PCR primer가 있는 strand는 cleavage된다. 그리고 sequencing primer는 다음과 같이 binding 한다(녹색 sequence가 sequencing primer).
 

5' -------------------- -----ACACTCTTTCCCTAC ACGACGCTCTTCCGATCT (----) -------------------- -------------- 3'

3' TTACTATGCCGCTGGTGGCT GTCCATGTGAGAAAGGGATG TGCTGCGAGAAGGCTAGA (DNA-) TCTAGCCTTCTCGAGCATAC GGCAGAAGACGAAC 5'

 

그리고 sequencing primer를 시작으로 3' (=>) 방향으로 extension 되면서 DNA+ 를 읽게된다.

 

지금까지 제가 이해한 것을 바탕으로 끄적끄적 대 봤습니다. 혹시나 이 원리가 잘못되어 있다면 누구든지 덧글 달아주세요~ :)

 

Aspera 고속 데이터 전송 프로토콜

얼마 전 논문을 내면서 시퀀싱 데이터를 ENA(European Nucleotide Archive)의 SRA(Sequence Read Archive)에 보낼 때 썼던 프로토콜이다. 처음에는 파일용량이 하나당 수십 GByte인지라 FTP는 불가능 할 것 같고 하드디스크를 구입해서 보내려고 했으나 ENA 쪽에서 Aspera라는 commercial 고속 파일 전송 프로토콜을 추천, 사용해봤는데 정말 빠르더라~

 

ENA-SRA에서 파일 전송에 대한 설명 페이지:

http://www.ebi.ac.uk/ena/about/sra_data_upload

맨 아래에 보면 Aspera로 전송하는 방법에 대해 나와있다.

여기서 제공하는 링크에서 Aspera ascp command line client를 받아 사용했다.

(CLIENT SOFTWARE 중 aspera connect를 받아야 함)

 

Command line의 형태는 아래와 같다.

 

ascp -QT -l300M -L- <file to upload><era-drop-N>@fasp.sra.ebi.ac.uk:/.

 

-l300M: 업로드 속도를 30MB/s로 한정. 안정적 전송을 위해서 이 수치를 줄일 필요가 있을지도 모른다고 하지만 무시하고 그대로 씀.

-L-: 전송 중 logging. 이 옵션을 사용하면 전송 중에 수 많은 log line들이 나온다. 안 써도 상관은 없겠지만 그래도 전송되고 있는 것을 확인하고 싶으면 사용하는게 좋을 듯

<file to upload>: 말 그대로 전송할 파일. "*.txt.gz" 와 같은 형태로 써 줘도 된다.

<era-drop-N>: 이건 ENA-SRA에서 제공되는 Aspera login ID

 

이렇게 한 줄 써주면 파일이 쓩쓩쓩 전송되요~ 그것도 아주 빨리!