Gibbs sampling을 이용한 Multiple alignment implementation (C++)

/*

 * http://sosal.kr/

 * made by so_Sal

 */

 

Gibbs sampling

 

1. Multiple alignment

2. Gibbs sampling

3. Motif란?

4. Gibbs sampling의 동작

5. Gibbs sampling의 실제 구현 (C++ programming)

6. 프로그램 실행 결과

 

1. Multiple alignment

DNA alignment는 DNA들을 정렬하여 비슷한 인자들을 찾아 유전학적으로 유용한 DNA서열을 뽑아내는 작업중 하나인데, 진핵생물의 특성상 유전자들 사이에서도 반복 DNA서열이 존재하기 때문에 분자생물학에서 많은 DNA들을 정렬하는건 큰 의미가 있는 작업입니다.

 

DNA alignment는 각각의 N1, N2 갯수만큼 유전자 정보(A,C,G,T 4가지 염기)가 있다고 가정했을때 N by N matrix를 만들어 유전사 서열이 일치하는 부분을 찾아내는 방법을 기본적으로 사용하는데 (일반적인 waterman 알고리즘) 이러한 방법은 아무리 최적화를 시켜도 DNA가 3개이상 들어가게되면 O(n^3) 알고리즘을 갖게되어 엄청나게 긴 시간이 요구됩니다.

 

이렇게 3가지 이상의 DNA를 정렬하여 유전정보를 찾아내는 작업이 Multiple alignment입니다.

 

 

figure 1.1 Multiple alignment

 

2. Gibbs sampling

위에서 설명한 Multiple alignment problem은 NP이며 이를 해결하기 위한 방법중 하나가 바로 gibbs sampling입니다.

gibbs sampling은 stochastic algorithm이며 100% 정확한 해답을 찾아내진 않지만 (전체 탐색이 아닌 확률을 사용한 부분 탐색) 빠르고 쉽게 다중 DNA의 반복서열을 찾아낼 수 있습니다.

 

 

3. Motif: 단백질이나 DNA서열에서 반복적으로 나타나는 짧은 패턴

 

ex)

CTACCTTTGACTAGTCCGTTCTAGCG

AGCTAAGTCGGTTCCTAGCATGCGTG

TTCTACCGCAGTCAGCTCCGATGCGC

GAGTGAGTTCACGTAGCCGTAGCTAA

TTGCTAAGTGAGTTCACGGCTTGGCT

data 3.1 DNA sequences

 

100% 일치하진 않더라도 비슷한 반복서열은 DNA의 Multiple alignment를 하기위한 중요한 key가 됩니다.

 

이러한 motif는 유전적으로도 의미가 있습니다. 전사인자의 합성부위의 DNA 서열을 찾아내는 방법으로도 사용할 수 있고 DNA와 protein의 결합부위도 예측할 수 있는 좋은 자료입니다.

 

4. Gibbs sampling의 동작

 

1) 유저에서 요구되는 Data

  1. Multiple DNA / RNA / Protein sequence

  2. Motif width.

  3. The number of Motif

 

2) input으로 주어진 DNA / RNA / Protein sequence에서 각 단위체들이 몇% 포함되어있는지 측정해야합니다.

DNA를 예로 들자면 AAGTC / 5개의 짧은 서열에서 A는 40%, 나머지 G,T,C는 20%만큼 포함되어있습니다.

 

3). 무작위적으로 sequence에서 motif의 시작위치를 선택합니다.

예로 100개의 염기가 존재하는 DNA 서열에서 motif의 길이가 10이라고 가정한다면

1~91, 총 91개의 위치에 해당되는 경우의수가 존재합니다.

 

4). PSSM (Position specific score matirx) 를 구성합니다.

 

예로 아래의 5개의 sequence를 준비하여 gibbs sampling을 이용하여 pssm을 구성해보겠습니다.

 

Seq 1. GTAAACAATATTTATAGC

Seq 2. AAAATTTACCTTAGAAGG

Seq 3. CCGTACTGTCAAGCGTGG

Seq 4. TGAGTAAACGACGTCCCA

Seq 5. TACTTAACACCCTGTCAA

 

data 4.1. input DNA seuqnce of gibbs sampling

 

 

4. 1-3) 작업에서 The number of Motif = 8로 선택하고

4. 3) 작업에서 motif의  위치를 {s1, s2, s3, s4, s5 | 7, 11, 9, 4, 1} 라고 선택헀다고 가정했을 시

 

Seq 1. GTAAACAATATTTATAGC

Seq 2. AAAATTTACCTTAGAAGG

Seq 3. CCGTACTGTCAAGCGTGG

Seq 4. TGAGTAAACGACGTCCCA

Seq 5. TACTTAACACCCTGTCAA

각 5개의 squence에서 무작위적으로 선택된 주황색 부분이 motif이며

이들중 단 하나의 문장을 choosed sequence로 지정하여 PSSM 테이블에 포함시키지 않습니다. 

(아래 예제에서 2번의 sequence를 선택한다고 가정하겠습니다.)

이 motif들의 PSSM은 바로 아래와 같습니다.

 

 

Figure 4.2 - data 4.1을 사용하여 나타낸 PSSM 테이블의 예

 

무작위적으로 선택된 motif를 가져와 motif의 첫번째에 해당되는 염기중 가장 많은것을 선택하여(Consensus string) 그 결과로 데이터를 담둡니다.

 

5) PSSM 에서 주어진 결과의 score를 측정합니다.

 

선택된 DNA 2번서열은 AAAATTTACCTTAGAAGG 입니다.

 

가장 왼쪽에서부터 motif의 크기만큼 우측으로 이동하면서 각각의 확률을 구해줍니다. 

가장 위의 경우인 AAAATTTACCTTAGAAGG 에서는 아래처럼 계산해주시면 됩니다.

 

Figure 4.2 scoring - result ex1)

 

0.25 * 0.5 * 0.5 * 0.75 * 0.5 * 0.25 * 0.25 * 0.5 = 0.000732421875 가 나옵니다.

이런식으로 선택된 DNA에서 DNA_size - size_motif+1 개만큼 경우의 수에 해당되는 score를 계산해주어 가장 높은 score를 선택합니다.

 

5) score를 구한 후 또 다른 위치에서의 motif를 생성하는 3)의 경우부터 다시 시작하여 5)의 scoring까지 계속 반복합니다. 그렇게 나온 score 값을 비교하여 if ( new score >= oldsocre ) than oldscore=newscore로 갱신하여 높은 scoring 값을 선택하고, 그때의 motif들의 값들을 best_motifs로 저장하여 원하는만큼 반복한 후 가장 socre가 높은 motif가 multiple alignment의 결과가 되겠습니다.

 

반복작업 6)을 하면 할수록 더 높은 scoring값이 선택될 확률이 높아질것입니다.

 

 

5. Gibbs sampling의 실제 구현 (C++ programming)

 

 

#include <iostream>

#include <string.h>

#include <process.h>

#include <vector>

#include <iomanip>

 

using namespace std;

 

#define size_motif 8

#define count_DNA 5

#define DNA_size 18

#define n_repeat 20000

 

class DNAs{

private:

int n_DNA; //DNA의 개수

vector<string> DNA; //DNA 서열이 저장될 vector

 

public:

DNAs():n_DNA(0){}; //생성자: n_DNA 초기화

        void add_string(string input){

DNA.push_back(input);

n_DNA++;

}

 

void print_strings(){

for(int i=0;i<n_DNA;i++){

cout<<DNA[i].c_str()<<endl;

}

}

 

string get_DNA(int input){

if(input>n_DNA){

return false;

}

return DNA[input];

}

 

string get_subDNA(int i,int pos){

return DNA[i].substr(pos,size_motif);

} // i번째 sequence에서 pos부터 motif 길이만큼 출력합니다.

};

 

int main(){

DNAs DNA;

string tmp_motif;

string best_motifs[count_DNA]; //random으로 선택된 motif중 가장 좋은 motif

string best_motifs_choosed_DNA; // PSSM에 사용되지 않았던 DNA의 motif

string best_consensus_motif;

double best_score=0;

 

DNA.add_string("GTAAACAATATTTATAGC");

DNA.add_string("AAAATTTACCTTAGAAGG");

DNA.add_string("CCGTACTGTCAAGCGTGG");

DNA.add_string("TGAGTAAACGACGTCCCA");

DNA.add_string("TACTTAACACCCTGTCAA");

//////////////////////////////////////// DNA input

 

for(int k=0;k<n_repeat;k++){ //n_repeat만큼 반복합니다.

double PSSM[4][size_motif]={0,};

double ratio[4]={0,};

tmp_motif.clear();

double score=1;

int pos_motif[count_DNA]={0,};

string get_sequence[count_DNA];

 

for(int i=0;i<count_DNA;i++){

for(int j=0;j<DNA_size;j++){

switch((char)DNA.get_DNA(i)[j]){

case 'A':

ratio[0]++;

break;

case 'C':

ratio[1]++;

break;

case 'G':

ratio[2]++;

break;

case 'T':

ratio[3]++;

break;

}

}

}

ratio[0]/=(count_DNA*DNA_size); //A 비율

ratio[1]/=(count_DNA*DNA_size); //C 비율

ratio[2]/=(count_DNA*DNA_size); //G 비율

ratio[3]/=(count_DNA*DNA_size); //T 비율

int choosed_sequence = rand() % count_DNA; //랜덤으로 하나의 DNA를 선택 

                // 이 DNA는 PSSM 테이블에 포함되지 않으며, PSSM에서 만들어진 각 염기의 확률을 이용하여

                // choosed_sequence의 최적의 mltif 위치를 찾아낼 것입니다.

 

for(int i=0;i<count_DNA;i++)

pos_motif[i]=rand()%(DNA_size - size_motif+1) + 1;

                //무작위적으로 motif의 위치를 선택한다.

 

for(int i=0;i<count_DNA;i++){

if( i == choosed_sequence )

continue;

get_sequence[i] = DNA.get_subDNA(i,pos_motif[i]-1);

// cout<<get_sequence[i].c_str()<<endl;

for(int j=0;j<size_motif;j++){

if(get_sequence[i][j] == 'A')

PSSM[0][j]++;

else if(get_sequence[i][j] == 'C')

PSSM[1][j]++;

else if(get_sequence[i][j] == 'G')

PSSM[2][j]++;

else if(get_sequence[i][j] == 'T')

PSSM[3][j]++;

}

}

for(int i=0;i<size_motif;i++){

for(int j=0;j<4;j++){

PSSM[j][i]/=count_DNA-1;//choosed sequence를 제외한 나머지 DNA 개수로 나눠줌

}

}

/////////////  PSSM 테이블 만들기

 

/////////////////////////////////////////

                /*cout.setf(ios::fixed,ios::floatfield);

cout.setf(ios::showpoint);

for(int i=0;i<4;i++){

for(int j=0;j<size_motif;j++){

cout.precision(3);

cout<<setw(4)<<PSSM[i][j]<<" ";

}cout<<endl;

}

                *////////////// PSSM 출력

 

for(int i=0;i<size_motif;i++){

double max=0;

int pos;//ACGT 0,1,2,3

for(int j=0;j<4;j++){

if(PSSM[j][i] > max){

max=PSSM[j][i];

pos=j;

}

}

switch(pos){

case 0:

tmp_motif+='A';

break;

case 1:

tmp_motif+='C';

        break;

case 2:

tmp_motif+='G';

break;

case 3:

tmp_motif+='T';

break;

}

}

//PSSM 테이블에서 MOTIF 추출 (가장 빈노가 높은 염기를 MOTIF로 지정): Consensus string

best_consensus_motif = tmp_motif.c_str();

 

for(int i=0;i<DNA_size-size_motif+1;i++){

score=1;

string score_motif = DNA.get_subDNA(choosed_sequence,i);

for(int j=0;j<size_motif;j++){

switch( score_motif[j] ){

case 'A':

score *= PSSM[0][j];

break;

case 'C':

score *= PSSM[1][j];

break;

case 'G':

score *= PSSM[2][j];

break;

case 'T':

score *= PSSM[3][j];

break; 

}//각 확률을 곱해준다.

}

if(score > best_score){

cout<<score_motif.c_str()<<endl;

cout<<score<<endl;

best_score = score;

best_motifs_choosed_DNA = score_motif;

for(int i=0;i<count_DNA;i++){

best_motifs[i] = get_sequence[i];

}

}//best score라면 각 motif와 score 등을 저장

}

}

cout<<"Best scores: "<<best_score<<endl;

cout<<"Consensus motif: "<<best_consensus_motif.c_str()<<endl;

for(int i=0;i<count_DNA;i++){

if(best_motifs[i].size() > 0)

cout<<i+1<<": "<<best_motifs[i].c_str()<<endl;

else

cout<<i+1<<": "<<best_motifs_choosed_DNA.c_str()<<" (choosed)"<<endl;

 

}

}

 

 

6. 프로그램 실행 결과

 

 

 

Seq 1. GTAAACAATATTTATAGC

Seq 2. AAAATTTACCTTAGAAGG

Seq 3. CCGTACTGTCAAGCGTGG

Seq 4. TGAGTAAACGACGTCCCA

Seq 5. TACTTAACACCCTGTCAA

 

2만번의 랜덤 motif를 설정하여 score를 계산하였더니 이러한 결과룰 도출해냈습니다. 

 

 

이 글의 상단부에 있는 data 3.1을 대입해볼까요?

CTACCTTTGACTAGTCCGTTCTAGCG

AGCTAAGTCGGTTCCTAGCATGCGTG

TTCTACCGCAGTCAGCTCCGATGCGC

GAGTGAGTTCACGTAGCCGTAGCTAA

TTGCTAAGTGAGTTCACGGCTTGGCT

 

DNA_size = 26

size_motif = 6

n_repeat = 50000 으로 수정하여 프로그램을 실행해보았습니다.

 

 

sequence 1,3,4,5번은 정확하게 찾아냈지만 seq.2는 제가 원했던 motif가 아니네요

하지만 n_repeat을 더 크게 조정한다면 정확한 값을 찾아낼거라 생각합니다.

 

프로그램 소스에 모르는 내용이 없도록 주석을 달려고 했지만

제 프로그래밍 버릇상 가독성이 그렇게 좋아보이진 않네요.

 

프로그램 소스에 의문이 있거나 혹은 gibbs sampling에 대한 질문이 있으시다면

언제든지 문의주세요. 답변해드리겠습니다.