📂알고리즘

니들맨-분쉬 알고리즘: 전역 서열 정렬

개요

두 염기서열의 공통 부분이 가장 많아지는 정렬을 찾는 것을 전역 정렬^{global Alignment}이라 하는데, 그 방법으로 가장 널리 쓰이는 니들맨-분쉬 알고리즘^{needleman-Wunsch Algorithm}을 소개한다.

서열정렬에는 너무나 많은 경우의 수가 있기 때문에 다이내믹 프로그래밍 을 통해 효율적이고 빠르게 계산할 필요가 있다.

알고리즘¹

Input

두 문자열 $\textbf{v}, \textbf{w}$ 과 치환행렬 $\delta$ 가 주어져 있고 다음과 같이 나타낸다고 하자.

$$ \textbf{v} := \left( v_{1} , \cdots , v_{m} \right) \\ \textbf{w} := \left( w_{1} , \cdots , w_{n} \right) $$

$\textbf{v}$ 의 길이는 $m$, $\textbf{w}$ 의 길이는 $n$ 이다. 두 캐릭터 $x,y$ 의 패널티는 $\delta (x,y)$ 와 같이 나타내고, 특히 갭을 나타내는 캐릭터는 $-$ 로 나타내기로 하자.

---

Step 1. 초기화

빈 행렬 $\left( s_{i,j} \right)_{i=0, \cdots, m \\ j = 0, \cdots , n}$ 을 만들고 다음과 같이 첫번째 행, 첫번째 열을 초기화한다.

$$ s_{i,0} \gets i \cdot \delta \left( v_{i},- \right) \\ s_{0,j} \gets j \cdot \delta \left( -, w_{j} \right) $$

---

Step 2. 동적계획법

$s_{i-1,j}$ 과 $s_{i,j-1}$ 이 이미 계산되어 있으면 $s_{i,j}$ 를 다음과 같이 계산한다. (물론, $s_{i-1,j-1}$ 은 계산되어 있을 것이다.)

$$ s_{i,j} \gets \max \begin{cases} s_{i-1,j} + \delta \left( v_{i}, - \right) \\ s_{i,j-1} + \delta \left( -, w_{j} \right) \\ s_{i-1, j-1} + \delta \left( v_{i} , w_{j} \right) \end{cases}$$

---

Ouput

$s_{mn}$ 부터 행과 열이 오른쪽 아래에서 왼쪽 위로 따라가며 일치하는 원소를 고르면 두 문자열과 가장 유사한 전역 정렬을 얻는다. 그 점수는 $s_{mn}$ 이다.

예제

Input

julia> s1 = "ATTAGCTCAC";

julia> s2 = "ACTATAGCGCAG";

위와 같은 두 염기서열을 전역 정렬하자.

julia> EDNAFULL
SubstitutionMatrix{BioSymbols.DNA,Int64}:
     A  C  M  G  R  S  V  T  W  Y  H  K  D  B  N
  A  5 -4  1 -4  1 -4 -1 -4  1 -4 -1 -4 -1 -4 -2
  C -4  5  1 -4 -4  1 -1 -4 -4  1 -1 -4 -4 -1 -2
  M  1  1 -1 -4 -2 -2 -1 -4 -2 -2 -1 -4 -3 -3 -1
  G -4 -4 -4  5  1  1 -1 -4 -4 -4 -4  1 -1 -1 -2
  R  1 -4 -2  1 -1 -2 -1 -4 -2 -4 -3 -2 -1 -3 -1
  S -4  1 -2  1 -2 -1 -1 -4 -4 -2 -3 -2 -3 -1 -1
  V -1 -1 -1 -1 -1 -1 -1 -4 -3 -3 -2 -3 -2 -2 -1
  T -4 -4 -4 -4 -4 -4 -4  5  1  1 -1  1 -1 -1 -2
  W  1 -4 -2 -4 -2 -4 -3  1 -1 -2 -1 -2 -1 -3 -1
  Y -4  1 -2 -4 -4 -2 -3  1 -2 -1 -1 -2 -3 -1 -1
  H -1 -1 -1 -4 -3 -3 -2 -1 -1 -1 -1 -3 -2 -2 -1
  K -4 -4 -4  1 -2 -2 -3  1 -2 -2 -3 -1 -1 -1 -1
  D -1 -4 -3 -1 -1 -3 -2 -1 -1 -3 -2 -1 -1 -2 -1
  B -4 -1 -3 -1 -3 -1 -2 -1 -3 -1 -2 -1 -2 -1 -1
  N -2 -2 -1 -2 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1
(underlined values are default ones)

치환행렬에 갭 페널티는 리니어하게 $-5$ 점이라고 하자. 그러면 쓰이는 캐릭터는 $A,T,C,G$ 네 개의 주요 염기뿐이고, 캐릭터가 일치할 때 $5$ 점을 받고 틀렸을 때 $-4$ 점, 갭 패널티로 $-5$ 를 받는다.

이해를 돕기 위해 엑셀에서 간단하게 구현해보자.

Step 1. 쿼리 시퀀스를 행으로, 레퍼런스 시퀀스를 열로 적고 초기화한다.

---

Step 2. 각 셀에서 스코어를 계산한다.

C3(3행 3열)에서 다음과 같이 함수를 입력하고 모든 칸으로 복사하면 된다.

=MAX(B2+IF(C$1=$A3,5,-4), C2-5,B3-5)

---

Output. 우측 하단부터 트레이싱한다.

오른쪽 아래 모서리부터 다음의 규칙에 따라 거꾸로 타고간다.

• 적색 ↖: 행과 열이 매치되지 않았으니 둘 중 어느쪽을 골라도 상관 없다.
• 황색 ↖: 행과 열이 매치되었으니 해당 염기를 선택한다.
• 회색 ← or ↑: 갭, 인델^indel을 의미한다. 상단과 좌측 중 상대적으로 점수가 높은 곳을 고른다.

트레이싱 과정에서 알 수 있듯 최적의 서열 정렬은 딱히 유일하지 않을 수 있다. 이렇게 찾을 수 있는 최적의 전역 정렬은 A-T-TAGCCAC 혹은 A-T-TAGCGAC 심지어 ACT-TAGCGAC도 될 수 있다. 그러나 찾을 수 있는 서열정렬들의 점수는 모두 22점으로 동일하다.

코드

다이내믹 프로그래밍인만큼 위와 같이 엑셀로도 간단하게 구현할 수 있다. 예제에서 보았듯 공부한다는 목적에는 괜찮지만 이걸 계속 재사용하는 것은 분명 미련한 짓이다. 다음은 줄리아에서 쉽고 빠르게 같은 작업을 수행하는 코드다.

using BioAlignments

s1 = "ATTAGCTCAC";
s2 = "ACTATAGCGCAG";

scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=-5, gap_extend=0);
global_result = pairalign(GlobalAlignment(), s1, s2, scoremodel)

실행 결과

using BioAlignments

julia> s1 = "ATTAGCTCAC";

julia> s2 = "ACTATAGCGCAG";

julia> scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=-5, gap_extend=0);

julia> global_result = pairalign(GlobalAlignment(), s1, s2, scoremodel)
PairwiseAlignmentResult{Int64,String,String}:
  score: 22
  seq:  1 A-T-TAGCTCAC 10
          | | |||| || 
  ref:  1 ACTATAGCGCAG 12