티스토리 툴바


Computer science에서 말하는 pattern matching technique에는 매우 많은 종류의 algorithm이 존재하지만, Bioinformatics 에서 nucleotide sequence 를 reference genome sequence 에 alignment 사용할 때 주로 사용되는 방법은 크게 다음의 두 가지로 나눌 수 있다.

1. Hash Table 을 이용하는 방법
2. BWT Index 를 이용하는 방법


1번 Hash Table을 이용하는 방법의 경우, (human genome 을 기준으로 했을 때) 기본적으로 필요로 하는 메모리의 양이 일반 PC 또는 서버의 Main memory 한계를 훨씬 넘어서고, Performance 측면에서도 많이 사용되지 않고 있다. (MAQ, SOAP1 등이 1번의 방식을 사용하고 있음.)

1번 Hash Table의 경우는 명칭만으로도 대략적인 방법을 유추할 수 있지만 2번의 경우는 BWT (Burrows-Wheeler Transform)을 알고있다고 하여도 이것이 Pattern Matching 과 어떤식으로 연관되어 사용되는지 유추하기가 어렵다. 이 내용을 이해하기 위해서는 먼저 BWT (Burrows-Wheeler Transform), SA (Suffix Array) 두가지 개념의 활용을 Implementation 측면에서 확실히 이해해야 한다.

1. Suffix Array


- 먼저, Suffix 라 함은, 한글로 접미사와 동의어이며 Pattern Matching 에서는 주어진 문자열 T에서 마지막문자를 포함하는 모든경우의 sub sequence 를 의미한다. 예를 들어 문자열 T를 "
BASDEGA" 라고 하였을 때 이 문자열의 Suffix 들은 아래와 같다.
              $  (start position : 7)
            A$  (start position : 6)
          GA$  (start position : 5)
        EGA$  (start position : 4)
      DEGA$  (start position : 3)
    SDEGA$  (start position : 2)
  ASDEGA$  (start position : 1)
BASDEGA$  (start position : 0)


- 문자열 T의 가장 마지막에 덧붙이는 "$" 는 알파벳에 존재하지 않는 특수문자이며 알파벳순서상 모든 알파벳 (A~Z) 보다 앞순서이라고 가정한다. 문자열 제일 마지막에 "$"를 붙이는 이유는 아래의 내용을 전부 이해하고 나면 알 수 있다.

- Suffix Array 라 함은, (Integer, long 또는 그 이상의 범위를 나타내는 type의) 정수를 원소로 갖는 배열이며, 이 배열의 각각의 원소로 들어가는 값은 주어진 문자열 T 상에서의 어떠한 suffix 의 start position (0-based 즉, 위의 예에서 B의 position 값은 0이다.) 값이다.

- 각각의 배열의 원소의 값은 아래와 같은 기준으로 정해진다.

- i 번째 원소의 값은 주어진 suffix 들을 알파벳순으로 정렬했을때 i 번째로 정렬되는 suffix 의 start position 값이다.
즉, 위의 예제에서 나온 문자열 T 의 suffix 들을 기준으로 했을 때, 알파벳순서에 의해서 정렬하게 되면

              $  (start position : 7)
            A$  (start position : 6)
  ASDEGA$  (start position : 1)
BASDEGA$  (start position : 0)
      DEGA$  (start position : 3)
        EGA$  (start position : 4)
          GA$  (start position : 5)
    SDEGA$  (start position : 2)


와 같게 되고, suffix array 의 각각의 원소의 값은 SA = {7, 6, 1, 0, 3, 4, 5, 2} 가 된다.



2. BWT (Burrows-Wheeler Transform)


- 원재 BWT는 "압축기법" 알고리즘으로서 처음 고안되었는데, 2000년경으로 접어들면서 Ferragina 와 Manzini 에 의해서 pattern matching 응용되기 시작하였다.


- BWT 그 자체의 내용에 대해서 궁금할 경우에는 http://hosang.tistory.com/entry/Burrows-Wheeler-Transform-BWT-Algorithm 를 참고하면 된다.

- 이 Post 에서는 Pattern Matching 에 있어서 사용되는 BWT 문자열 (character array) 에 대해서만 자세하게 설명한다.

- BWT index 기법에서 사용되는 BWT 라는 단어는 하나의 문자열을 의미하며 더욱 정확하게 말해서 Character Array 를 의미한다. 이 character 배열의 i 번째 원소로 들어가는 값은 suffix array 의 i 번째 원소값에서 1을 뺀 숫자를 a 라고 했을때 문자열 T 의 a 번째 문자이다 (0-based). Suffix array 의 i 번째 값이 0 일 경우는 문자열 T의 가장 마지막 문자 "$" 가 된다.

- 즉, 위의 예제 문자열을 기준으로 생각해보면 suffix array 의 값이 {7, 6, 1, 0, 3, 4, 5, 2} 이므로 BWT의 각 원소의 값은 문자열 T의 (0-based)

{6번째 문자, 5번째 문자, 0번째 문자, -1번째 문자($), 2번째 문자, 3번째 문자, 4번째 문자, 1번째 문자}

이며, 이것은 다시말해
{"A", "G", "B", "$", "S", "D", "E", "A"} = "
AGB$SDEA" 가 된다.

- 다시 정리하자면, 아래와 같다.

문자열 T     : "BASDEGA$"
Suffix Array : {7, 6, 1, 0, 3, 4, 5, 2}
BWT           :
"AGB$SDEA"


위의 내용이 BWT Index 를 Pattern matching 에 사용하는 방법을 이해하기 위한 기본개념이지만, 이 것만으로는 정확하게 방법을 이해할 수는 없다.
지금까지의 내용이 다소 이해하기 쉬웠다고 한다면, 이제부터가 조금 이해하기 난해한 부분이다.



3. Backward Search Technique


먼저 Pattern X 에대한 SA[i] 와 SA[j] 값을 아래와 같다고 가정하자.

SA[i] = The smallest suffix of T that have X as the prefix.
SA[j] = The largest suffix of T that have X as the prefix.

그러면 range [i, j] 를 우리는 SA Interval (또는 range) of X 라고 부를 수 있다.

Pattern X에 대한 SA interval [i, j] 가 주어지면 우리는 pattern zX (z: any character in T) 의 SA interval [p, q]를 "Backward Search Technique" 에 의해 constant time 내에 구할 수 있다. (Ferragina, P. and Manzini G. (2000) Opportunistic data structures with applications. FOCS. 390-398)

Lemma 1. X 를 string, z 를 character 라고하고, X의 SA interval을 [i..j], zX의 SA interval 을 [p..q] 라고 가정하면 p, q 는 다음과 같이 구해질 수 있다.
p = C(z) = Occ(z, i-1) = 1
q = C(z) = Occ(z, j)
여기서 C(z) 는 문자열 T에서 알파벳순서상 z 보다 작은 character 들의 총 수, Occ(z, i) 는 BWT[0..i] 에서 나타나는 character z 수 를 의미한다.


위의 Lemma 1에 의해면, C(z) 는 pre-compute를 통해서 미리 계산해 두고 constant time 내에 retrieve 할 수 있고, Ozz(z, i) 역시 constant time 내에 계산 가능하므로, [i, j] 를 통해 [p, q]를 구하는 것은 constant time 내에 가능하다고 말할 수 있다.


Reference : T. W. Lam et al. (2008) Compressed indexing and local alignment of DNA. Bioinformatics, 24, 791-797
Posted by  Hosang Jeon
Posted by Abraham

MAQ


Speed O(c1*NlogN + c2*L + c3*2^- k*NL).

 

N : # of reads (typically 2 million)
L : length of reference
k : used bit size for indexing (24 bits)

NlogN : time spent on sorting the indexes
L : time spent on scanning the whole reference sequence
2^-k*NL : time spent on processing the alignment when there is a seed hit.

ie.)
         sorting the index : O(NlogN)
      scanning reference : O(L)
         processing alignment : O(2^-k * NL)

Memory
                            
                 - Alignment : 1GB
                 - Consensus calling : 2GB

Accuracy :
            
                 - Mapping Rate : 97.44 %
                 - Mapping Error Rate : less than 1%

Read Length : ~ 63bp

Paired-End : YES

Advantages : It takes quality scores for each base. 

Preliminary functions to handle ABI SOLiD data (color-space data). 
It has a graphical viewer, but again for its own alignment format.

Disadvantages : For long reads (63bp~ ) MAQ can loose the limit at the cost of speed and memory usage. 
Only cover ungapped alignment.
Does not support gapped alignment for single-end reads
Speed problem
Output format is not SAM(Sequence Alignment/Map) format

Description : Hash table-based method

Posted by Abraham

1. BWT compression index 통해 reference indexing 하는 방식
2. space complexity : O(n/4) bytes
3. local alignment : Smith-Waterman
4. inexact match : "split-read strategy" - 1 mismatch 허용하기 위해 read 2 fragments split, 2 mismatch 허용하기 위해 read 3 fragments split
5. Paired-end alignment : MAQ 유사한 방식. 먼저, 각각의 reads 독립적으로 alignment , orientation relationship proper distance 통해 pair 찾음.
6. read data 5'에서 3' 으로 가며 quality 떨어지기 때문에 high-quality 5' region 에서는 2 mismatch 허용하고, 3' region 에서는 많은 mismatch 허용하는 option 추천하고 있음.
7. 1024 bp read data 까지 handling 있음.
8. read quality score 고려 .
Posted by Abraham