생명과학

Sequence Motif 등장 확률

전체 n ntd로 이루어진 target RNA 한 가닥이 있다고 하자. 이 RNA는 완전히 랜덤하게 만들어진 가닥이라고 가정한다. 즉, 각 위치에 A, U, G, C가 같은 확률로 존재할 수 있다.

이제 k-mer sequence motif가 하나 있다고 하자. 이 sequence motif는 특정되어 있다.

이 때 Target RNA 가닥에서 이 motif가 한 번이라도 등장할 확률은 얼마일까?

우연한 기회에 이 문제에 대해 생각해보게 되어 여사건의 확률 등의 방법을 생각해보다 이런 방법으로는 쉽게 풀리지 않는다는걸 알게 되었다.

검색을 해봐도 정확한 해답을 제시하는 사람은 많지 않았다 (너무 쉬운 문제여서 인지?). 여기서는 확률의 inclusion-exclusion principle을 이용한 해답을 소개하려 한다.

마주서기 문제

본 문제와 비슷하지만 비교적 단순한 문제를 먼저 살펴보자.

n쌍의 부부가 있다. 한 쪽에는 남편끼리, 다른 쪽에는 부인끼리 늘어 설 때, 적어도 한 쌍의 부부가 서로 마주보고 설 확률은 얼마일까?

이 문제 역시 여사건의 확률을 이용해서 풀기는 쉽지 않다 (재귀함수를 만들게 된다). 대신 inclusion-exclusion principle을 사용해서 풀어보자.

사건 A_i를 “i번째 부부가 서로 마주보고 서는 사건”이라고 하자. 그러면

\begin{array}{lcl} &&P(A_i) = \dfrac{(n-1)!}{n!} \text{, } i=1, ..., n \\ &&P(A_i \cap A_j) = \dfrac{(n-2)!}{n!} \text{, } i < j \\ && ... \\ &&P(A_1 \cap A_2 \cap ... \cap A_n) = \dfrac{1}{n!} \end{array}

이므로

\begin{array}{lcl} P(\cup_{i=1}^{n} A_i) &=& \sum_{i=1}^{n}P(A_i) - \sum_{i<j}P(A_i \cap A_j) \\ & & + ... + (-1)^{n+1}\sum_{i<j<...<l}P(A_i \cap A_j \cap ... \cap A_l) \\ &=& \sum_{i=1}^{n} (-1)^{i+1} \dfrac{1}{i!} \end{array}

이다.

본 문제

사건 A_i \text{, } i=1, ..., m를 “i번째 위치에 motif가 등장하는 사건”이라고 하자. 단 mnk로 나눈 몫이다. 그러면

\begin{array}{lcl} &&P(A_i) = \dfrac{4^{n-k}}{4^n} \text{, } i \leq n-(k-1) \\ &&P(A_i \cap A_j) = \dfrac{4^{n-2k}}{4^n} \text{, } i<j \text{ and } i, j \leq n-2(k-1) \\ &&P(A_i \cap A_j \cap A_k) = \dfrac{4^{n-3k}}{4^n} \text{, } i<j<k\text{ and } i, j, k \leq n-3(k-1) \\ &&... \\ &&P(A_1 \cap A_2 \cap ... \cap A_m) = \dfrac{4^{n-mk}}{4^n} \end{array}

이므로 마주서기 문제에서와 같은 방법으로

P(\cup_{i=1}^{m} A_i) = \sum_{i=1}^{m} (-1)^{i+1} \binom{n-i(k-1)}{i} 4^{-ki}

임을 알 수 있다.

이 방법으로 1~500 ntd의 random target sequence에 3, 4, 7-mer fixed motif가 존재할 확률을 그래프로 나타내보면 다음과 같다.

Thu Oct 25 00:15:46 2018

Update:

위 방법이 일부 상황에서 정확한 값을 반환하지 않는다고 한다. 더 일반적인 경우에도 사용 가능한 Markov chain 방법을 쓴 Udaqueness 님의 글을 참고.


참고

FASTA 압축기 만들기

DNA 시퀀스를 다루는 사람들에겐 FASTA/FASTQ라는 이름만큼 익숙한 포맷이 또 없을 것이다. 시퀀싱된 염기서열을 저장하는, 가히 표준이라 할 수 있을만큼 널리 사용되는 포맷이다. FASTA 포맷은 샘플 이름과 샘플 시퀀스로 이루어져있다. 예를 들어 “sample1”, “sample2″이라는 이름의 샘플이 각각 “ATGCATGC”, “TTTTTTTTT”라는 시퀀스로 시퀀싱되었다면 아래와 같은 FASTA 포맷의 파일로 표현할 수 있다.

>sample1
ATGCATGC
>sample2
TTTTTTTTT

FASTA 포맷의 가장 큰 장점은 굳이 포맷이라고 이름붙이는 것이 어색해 보일 정도의 간결성과 직관성에 있다. Parsing 작업이 굳이 필요하지 않으며, 필요하더라도 아주 쉽게 코딩할 수 있다 (물론 Biopython등 상용화된 패키지를 이용하는 것이 제일 편하다).

굳이 단점을 꼽자면 용량 측면에서의 비효율성이 있을 것이다. 각 베이스를 베이스 그대로 표기하여 직관성을 얻은 대신 하나의 베이스마다 8비트(1바이트)의 메모리를 할당해야하는 것이다.

비록 그렇게 큰 용량은 아니지만, 빠르고 원활한 시퀀스 파일 공유를 위해 FASTA 압축 알고리즘을 짜보면 어떨까 싶어 재미삼아 몇 가지 코딩을 해봤다.

이 글은 “새로운”(?) FASTA용 무손실 압축 알고리즘을 짜보고 싶다는 생각으로 시작했지만 결국 상용화된 프로그램이 짱이라는 결론에 이르는, 컴퓨터 비전공자의 분투를 담은 글이다.

1. Look-and-Say 수열을 이용한 압축 알고리즘 (.say format)

Look-and-say 수열은 현재 항 숫자의 종류와 갯수를 읽고, 읽은 내용으로 다음 항을 결정하는 수열이다. 수열의 첫 항이 1일 때, “1이 한 개 있다”고 읽을 수 있다. “1이 한 개(1)”이므로 이 수열의 둘째 항은 11이된다. 둘째 항은 11이므로 “1이 두 개 있다”고 읽을 수 있다. “1이 두 개(2)”이므로 셋째 항은 12가 된다. 이 과정을 계속해서 반복하면 다음과 같은 수열이 만들어진다.

[1]: 1
[2]: 11
[3]: 12
[4]: 1121
[5]: 122111
...

베르베르의 소설 ‘개미’에도 등장하여 우리나라에선 ‘개미 수열’이라는 이름으로도 알려져 있다. 바로 이 수열을 시퀀스 데이터 압축에 사용하면 어떨까 싶었다.

  • 1바이트(8비트) 중 3비트는 베이스 종류를 나타내는 데에 사용한다.
    • 000은 A, 001은 T, 010은 G, 011은 G, 100은 N, 101은 U
    • 110과 111은 비워두었다.
  • 나머지 5바이트는 반복수를 나타내는 데에 사용한다.
    • 00000 ~ 11111, 즉 1부터 32까지의 반복수를 표현할 수 있다.
    • 32번 이상 반복되는 베이스는 1바이트로 표현할 수 없다. 두 개 이상의 바이트를 써야한다.

예를 들어 AAAAAAATGGGGGGGGG라는 시퀀스는 이 방법으로 3 바이트로 표현할 수 있다.

[-1 byte-] [-2 byte-] [-3 byte-]
[00000110] [00100000] [01001000]

같은 시퀀스를 FASTA로 저장하려면 17바이트가 필요하다.

FASTA 포맷에 비교한 .said 포맷은 장점은 다음 두 가지로 요약할 수 있다.

  1. Homopolymer가 많은 시퀀스일수록 압축 효율이 높아진다. poly A region과 같은 시퀀스에서는 비교도 안될 정도로 작은 용량의 파일을 만드는 것을 확인할 수 있다.
  2. 반복이 거의 없더라도 같은 시퀀스가 저장된 FASTA보다는 작은(<=) 파일을 만든다. 이는 한 베이스가 반복되지 않는 경우에도 FASTA에서와 같이 1바이트 용량만을 차지하기 때문이다.
    반복수를 (대체로) 늘려주는 Burrows-Wheeler transformation을 사용한 후 `.said` 압축을 하면 좀 더 좋은 압축률을 얻을 수 있지 않을까? 했는데 생각만큼 크게 개선되진 않았다.

이 스크립트를 사용하면 .said 포맷으로 FASTA 변환할 수 있다: [압축, 압축해제]

# COMPRESS
python fa2said.py input.fa
# DECOMPRESS
python said2fa.py input.said

2. 한 베이스를 2 비트만으로 표현하기 (.two format)

genome / CDS 시퀀스는 생각보다 반복이 그리 많지 않았다. .said format으로는 높은 압축률을 얻기 어려운 상황이다. 그래서 아예 반복수를 기록하지 않고 한 베이스가 차지하는 용량을 극도로 줄이면 어떨까 싶었다.

압축 할 시퀀스에 A, T, G, C 네 종류의 베이스밖에 없다고 가정하면 (2^2개 이므로) 2비트만 사용하면 한 베이스를 인코딩할 수 있다. 한 베이스를 표현하는 데에 8비트를 사용해야 했던 FASTA 포맷에 비해 용량을 무려 1/4로 줄일 수 있는 것이다. 이렇게 압축한 파일 포맷을 .two 포맷이라고 부르도록 하자.

ATGCAATC의 경우 .two 포맷으론 2바이트(16비트)로 인코딩할 수 있다.

A T G C    A A T C
[00011011] [00000111]

이 스크립트를 사용하면 FASTA파일을 .two 포맷으로 압축할 수 있다: [압축, 압축해제]

# COMPRESS
python fa2two.py input.fa
# DECOMPRESS
python two2fa.py input.say

맺으며

이 기회에 야매로 압축 알고리즘을 공부하면서 첫 번째 시도(.said)가 run-length encoding이라 불리는 가장 원시적인 압축 알고리즘이라는 것을 알았다.

두 번째 시도 (.two)는 Huffman coding과 상당히 닮은 구석이 있다. 차이점은 DNA/RNA의 경우 인코딩하려는 데이터의 종류(A, T, U, G 또는 C, 추가로 N)가 이미 알려져있다는 것이다. DNA의 경우 데이터의 종류가 4개(A, T, G, C) 뿐이기 때문에 사람이 직접 2비트만으로 모든 데이터를 인코딩할 수 있었다.

두 번째 시도는 생물정보학계의 영웅 Jim Kent가 과거에 이미 .2bit라는 이름으로 상용화해서 쓰이고 있는 포맷이었다 [Jim Kent’s Repo]. C로 짜인 코드여서 내 것보다 훨씬 빠르게 압축할 수 있다. .2bit 포맷은 원하는 샘플의 시퀀스만 가져올 수도 있도록 각 샘플의 해시 딕셔너리도 함께 저장한다.

DNA 시퀀스와 달리 데이터 종류가 알려져 있지 않은 경우엔 Huffman coding을 사용할 수 있다. Huffman coding은 binary tree를 사용해서 자동으로 최적의 인코딩을 찾아준다. .two는 binary tree를 사용하지 않고 사람이 직접 인코딩 딕셔너리를 정해주는 특수한 변형이라고 할 수 있겠다.

.zip 압축에 사용되는 Lempel-Ziv 알고리즘은 인코딩 딕셔너리가 Huffman coding과 같이 고정되어 있지 않고 현재까지 압축한 데이터에 따라 유동적으로 변한다. 비슷한 패턴이 반복해서 나타나는 파일에 있어서 효율적인 압축을 가능케 한다. 압축을 풀기 위해서는 데이터 인코딩에 사용한 binary tree의 정보를 어떤 형태로든 압축된 결과 파일과 함께 저장해야하는 Huffman coding과 달리, 압축된 데이터 그 자체에 인코딩 딕셔너리가 포함된 것과 마찬가지이기 때문에 별도의 정보를 저장할 필요도 없다.

코딩한 결과물 자체는 성능 면에서나 속도 면에서나 무쓸모하지만 그 과정에서 많이 배웠다.

참고

23andMe로 무엇을 알 수 있나

기술이나 생명과학에 관심있는 이라면 스쳐지나가는 말로 한 번 쯤은 23andMe라는 것을 들어본 적이 있을 것이다. 23andMe는 개인 유전체 분석 서비스를 제공하는 생명과학 기업이다. 이전의 유전체 분석이 주로 ‘표준(reference)’이라 불리는 소수의 인간 개체에서 비롯된 유전체를 분석하는, 일반화가 필요한 작업이었다면, 개인 유전체 분석은 말 그대로 분석을 원하는 개개인의 유전체를 분석해서 개별화된 결론을 제공하는 서비스다.

kitlaydown.f45ac485ec3c

23andMe kit (23andme.com)

150~200 USD 상당의 서비스를 구입하면 23andMe에서 여러분의 집으로 키트를 보내준다. 그 키트에 침을 뱉어 샤카샤카 흔든 다음 다시 회사로 전송하면 그 침에 있는 유전정보를 바탕으로 가능한 모든 정보(심지어 여러분의 귓밥이 말랐는지 젖었는지까지!)를 리포트해준다,는 개념이다. 23andMe는 침을 기반으로 한 유전체 분석 키트를 상용화시켜 타임 지의 ‘올해의 발명품’에도 이름을 올린 바 있다. 그 외에도 구글 설립자 중 한명인 세르게이 브린의 아내의 회사로도 알려져있고 서비스 특성상 FDA와의 아슬아슬한 줄타기로도 여러차례 언론에 등장한 바 있다.

요 평범한듯 기발한 사업구상의 핵심 서비스는 크게 두 가지라고 할 수 있을게다.

첫째는 생물학적 계보 리포트다. 나의 뿌리는 어느 대륙에서 시작했는지, 어떤 민족의 유전형질이 얼마나 섞여있는지를 리포트해준다는 것이다. 더욱 거슬러 올라가 자신에게 네안데르탈인의 유전자가 얼마나 섞여있는지 또한 알 수 있다. 실용성보다는 재미를 목표로 만든 서비스가 아닐까 싶다.

둘째는 질병 리포트다. 이러이러한 대립형질을 가지고 있는 사람들은 어떠어떠한 질병을 가질 가능성이 높으니 미리 관리 좀 하시라!는 식의 개인화된 질병 리포트를 받는다는 것은 가타카를 비롯한 공상과학 영화에서도 자주 만날수 있었던 신나는 상상이다. 150~200달러라는 (굉장히 많이 저렴해졌지만) 여전히 값나가는 서비스를 제공받고자 검토하는 소비자의 입장에서는 대부분 이 질병 리포트를 목적으로 23andMe에 등록하는 것으로 보인다.

당신과 비슷한 유전체를 가지는 사람들이 잘 걸리는 질병을 바탕으로 개개인의 질병을 예측하는 것이기에 23andMe의 질병 리포트가 더욱 정확해지고 쓸모가 있어지려면 이 회사의 서비스를 이용하는 사람이 많아져야 한다. 다행(?)히도 올해 4월 서비스 이용자가 200만 명을 돌파하는 등 승승장구하는 모습을 보면 데이터의 신뢰성을 어느정도 높게 쳐줄 수 있을 듯 하다.

23andMe는 자체 데이터베이스를 기반으로 제공하는 리포트 외에도 사용자가 직접 자신의 유전체를 분석할 수 있도록 raw data를 제공하고 있다. 이 글은 23andMe에서 제공하는 서비스 외, 내가 직접 정보를 얻어오기 위해 작성한 파이썬 스크립트에 대한 이야기다. (서론이 매우 길었다)

23andMe Raw Data

본인은 아직 23andMe를 이용해보지 않았으므로(…) 다른 분이 인터넷에 공유해 둔 raw data로 몇 가지 간단한 리포트를 제공하는 파이썬 모듈을 만들어보고자 한다. 23andMe에서 제공하는 raw data는 다음과 같이 생겼다.

스크린샷 2017-05-20 오후 6.46.44

데이터에 대한 간단한 설명과 함께 탭으로 분리된 텍스트 형식의 데이터가 나온다. 내가 사용한 다른 분의 데이터는 2011년, 아마도 일루미나 HumanHap550 BeadChip을 이용해서 만들어진 것으로 보이며 레퍼런스 지놈으로는 NCBI36(hg18) 버전을 사용했다.

  • chromosome, position은 그 말 그대로 각각 염색체와 염색체상의 위치이다.
  • rsid(Reference SNP cluster ID)는 특정 SNP(single nucletide polymorphism; 연속하지 않은 한 염기에 변이가 일어난 것)의 구별자이다. 23andMe 뿐만 아니라 국제적으로 널리 사용되는 명명법이라고 한다. 주로 ‘rs’로 시작한다.
  • rsid가 등록되어있지 않은 SNP의 경우 23andMe 내부에서 자체적으로 할당한 id를 사용한다. 이건 ‘i’로 시작한다.
  • genotype은 해당 SNP에서의 genotype을 말한다. 두 글자는 각각 엄빠로부터 받은 SNP을 의미한다. 23andMe에서는 모든 genotype을 (+) strand를 기준으로 작성한다고 한다. dbSNP 등 다른 데이터베이스에는 (-) strand의 genotype으로 기재되어있는 경우도 있어 주의를 요한다.

MyGeno: 23andMe 분석 툴

이 데이터로 세 가지 종류의 리포트를 제공하는 python module (GenoMe)을 작성해 보았다. 튜토리얼 따라 명령어 몇 줄로 쉽게 이용할 수 있다. 23andMe 서비스를 이용하고 있는 사람이라면 raw data를 다운로드받아 사용할 수 있다.

1. Wellness Report

알코올 분해를 잘 하는지, 카페인 섭취는 어떻게 되는지, 유당 분해는 가능한지 등 약 8개 항목에 대해 리포트를 출력한다. 23andMe API를 이용했다. 출력 예시는 아래와 같다.

============================== WELLNESS REPORT ============================== 

< Alcohol Flush Reaction >
------------------------------------------------------------------------------------
The marker we tested comes in two different forms, the G variant and the A variant. The A variant results in an enzyme that is less efficient at breaking down acetaldehyde. The A variant is also known as c.1510G>A, Glu487Lys, and Glu504Lys. This marker has been studied the most in people of East Asian descent.
  rsID: rs671
  Yours: G;G - homozygous
  Detail: 2 | Alcohol Flush: Normal, doesn't flush. Normal hangovers. Normal risk of Alcoholism. Normal risk of Esophageal Cancer. Disulfiram is effective for alcoholism. | 60.2% of JPT

...

2. rsid로 검색

원하는 rsid의 SNP에 대해 나의 genotype과 함께 연관 형질을 출력한다.

rs28897696 | G;G | 0 | normal | 0.0% of JPT

3. 형질로 검색

궁금한 형질을 입력하면 이와 관련된 모든 SNP와 연관 정보를 출력한다. 다음은 “baldness”를 입력한 출력 결과.

rs6152 | G;G | 0.5 | able to go bald | 100.0% of JPT
rs2223841 | T;T | 1.2 | more likely to go bald before age 40 | 100.0% of JPT
rs6152 | G;G | 0.5 | able to go bald | 100.0% of JPT
rs2223841 | T;T | 1.2 | more likely to go bald before age 40 | 100.0% of JPT
rs2180439 | C;T | 2.5 | Increased risk of Male Pattern Baldness. | 44.2% of JPT
rs11683401 | C;T |  | --% of JPT
rs2073963 | G;G | 2.5 | increased risk of baldness | 18.6% of JPT
rs1511061 | T;T |  | --% of JPT
 * more details at https://www.snpedia.com/index.php/baldness

“food allergy” 등도 좋은 예시가 될 듯.

SNP 정보가 주로 서구권에 대해서 축적된 것을 감안하여 1000 genomes project의 데이터를 포함시켰다. 기본값으로 JPT(일본 도쿄의 사람들)를 사용해서, 군집 전체에서 나와 같은 genotype을 갖는 사람들의 비율을 보여준다. 23andMe 또한 데이터가 서구권에만 집중되었다는 비판을 받곤 하는데, 이에 아시아권 군집에서의 genotype 비율을 제공하여 참고치로 사용할 수 있도록 했다. JPT 외의 다른 군집을 사용하려면 여기여기를 참고.

이미 200만 명이 넘는 사람들의 유전체 데이터를 축적한 23andMe의 리포트에 비할 바는 아니지만, 새롭게 표현형과의 연관성이 밝혀진, 아직 23andMe에는 등록되지 않은 SNP, 혹은 23andMe에서 제공하지 않는 궁금한 형질에 대해서는 쉽고 빠르게 대략의 정보를 얻을 수 있을 것이다.