BioPlayground

🧬
목록으로

DP 전역 정렬: Needleman-Wunsch의 수식 유도와 파이썬 20줄 구현

1970년 프린스턴에서 태어난 서열 정렬의 원조. 점화식은 어떻게 유도되며, 왜 O(mn) 시간과 공간을 쓰는지, 파이썬 20줄로 어떻게 완전 구현하는지 손 계산부터 코드까지.

입문
|
15분
|
검증 완료 (2026-07-19)
sequence alignmentdynamic programmingRosalind
진행률0/5 (0%)

왜 이 편을 먼저 다뤄야 하는가

서열 정렬은 바이오인포매틱스의 알파벳이다. BLAST도, BWA도, 심지어 AlphaFold의 MSA(Multiple Sequence Alignment)도 결국 이 편에서 다룰 두 서열 정렬 문제의 확장에 지나지 않는다. 그런데 국내 한국어 자료를 뒤져보면 대부분 "그냥 BLAST 쓰면 된다"고 말한다. 왜 BLAST가 그렇게 빠른지, BLOSUM은 왜 그런 숫자를 쓰는지, 답할 수 있는 자료가 사실상 없다.

이 편에서는 도구를 쓰는 게 아니라 점화식을 유도하고, 격자를 손으로 채우고, 파이썬으로 밑바닥부터 구현한다. 우리가 오늘 손으로 채운 표가 내일 실무에서 다루게 될 BWA의 시드 확장이 되고, 모레 배울 HMM Viterbi의 등뼈가 된다. 서두르지 말자.

1970년 프린스턴의 오후

Saul Needleman과 Christian Wunsch가 J Mol Biol에 논문을 발표한 것이 1970년 3월이다. 그들이 던진 질문은 지금 우리한테는 너무 익숙하다.

"두 단백질 서열이 있다. 이들이 얼마나 비슷한가?"

지금은 당연해 보이지만, 그때는 "얼마나 비슷한가"의 정의부터 없었다. 어떻게 정렬해야 최선인가? 예를 들어 두 서열 GATTACA와 GCATGCU를 두고, 어떤 사람은 앞에서부터 하나씩 붙이고, 어떤 사람은 뒤에서부터 붙이고, 어떤 사람은 가운데를 벌린다. 정답이 뭘까?

Needleman과 Wunsch가 낸 답은 놀랍도록 단순하다. 모든 가능한 정렬 중에서 점수가 가장 높은 것. 다만 "모든 가능한 정렬"의 수가 서열 길이 m, n에 대해 지수적으로 폭발한다는 것이 문제였다. 그들은 이걸 동적 계획법으로 해결했다. Bellman이 DP를 이름 지은 것이 1953년이니까 겨우 17년 지난 시점에 바이오에 들어왔다.

점수 함수 · 격자 · 점화식

점수는 어떻게 매기나

두 서열을 붙일 때 세 가지 이벤트가 일어난다.

  1. 일치 (match): 같은 문자끼리 대응 → +s_match (예: +1)
  2. 불일치 (mismatch): 다른 문자끼리 대응 → -s_mismatch (예: -1)
  3. 갭 (gap): 한쪽에 빈 자리 → -d (예: -2)

예를 들어 아래와 같은 정렬이 있다면:

text
G A T T A C A
G - A T T A C

일치 3번 (G, A, T), 불일치 3번 (A-A, T-T, A-C), 갭 1번 (--). 매치 +1, 불일치 -1, 갭 -2를 대입하면 총점 3 × 1 - 3 × 1 - 1 × 2 = -2.

격자를 채우자

핵심은 이거다. 서열 X = x₁x₂…xₘ과 Y = y₁y₂…yₙ의 최적 정렬 점수 F(i, j)를 다음처럼 정의한다.

text
F(i, j) = max {
    F(i-1, j-1) + s(xᵢ, yⱼ)   ← 대각선: match / mismatch
    F(i-1, j  ) - d            ← 위에서: X에 문자, Y에 갭
    F(i  , j-1) - d            ← 왼쪽에서: Y에 문자, X에 갭
}

여기서 s(xᵢ, yⱼ)는 두 문자가 같으면 매치 점수, 다르면 불일치 점수다.

경계 조건은:

  • F(0, 0) = 0
  • F(i, 0) = -i · d
  • F(0, j) = -j · d

첫 행과 첫 열은 계속 갭만 쌓는다는 뜻이다.

왜 이 점화식이 옳은가

F(i, j)가 X의 앞 i글자와 Y의 앞 j글자의 최적 정렬이라면, 그 정렬의 마지막 열은 세 가지 중 하나다. 대각선(둘 다 문자), 위(X만 문자), 왼쪽(Y만 문자). 그 세 가지 각각에서 이전 문제 F(i-1, j-1), F(i-1, j), F(i, j-1)의 최적해를 이미 안다면, 그중 최대값이 지금의 최적해다. 이게 최적 부분 구조 (optimal substructure) — DP가 성립하는 이유다.

손으로 채우기

X = GATTACA, Y = GCATGCU, 매치 +1, 불일치 -1, 갭 -2로 격자를 채워보자. 지면 낭비를 줄이기 위해 짧게 자르면 X = GAT, Y = GAC 정도로 해보자.

-GAC
-0-2-4-6
G-2+1-1-3
A-4-1+20
T-6-30+1
  • F(1,1): G vs G — 매치 → F(0,0)+1 = 1. 대각선 승.
  • F(2,2): A vs A — 매치 → F(1,1)+1 = 2.
  • F(3,3): T vs C — 불일치 → max(F(2,2)-1, F(2,3)-2, F(3,2)-2) = max(1, -2, -2) = 1.

우하단 값 +1이 최적 점수다. 여기서부터 화살표를 거꾸로 따라가면 (traceback) 실제 정렬이 튀어나온다.

text
G A T
G A C

매치 2번 + 불일치 1번 = 2 - 1 = 1. 격자 값과 일치한다.

파이썬 20줄로 완전 구현

python
def needleman_wunsch(x, y, match=1, mismatch=-1, gap=-2):
    m, n = len(x), len(y)
    F = [[0] * (n + 1) for _ in range(m + 1)]


    # 경계 조건
    for i in range(m + 1):
        F[i][0] = i * gap
    for j in range(n + 1):
        F[0][j] = j * gap


    # 격자 채우기
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            s = match if x[i-1] == y[j-1] else mismatch
            F[i][j] = max(
                F[i-1][j-1] + s,
                F[i-1][j]   + gap,
                F[i][j-1]   + gap,
            )
    return F[m][n], F


score, table = needleman_wunsch("GATTACA", "GCATGCU")
print(score)  # 최적 점수

이게 다다. 20줄. 실무에서 쓰이는 EMBOSS needle도 결국 이 뼈대다.

Traceback 붙이기

정렬을 실제로 뽑고 싶다면 traceback 함수를 추가한다.

python
def traceback(F, x, y, match=1, mismatch=-1, gap=-2):
    i, j = len(x), len(y)
    ax, ay = "", ""
    while i > 0 or j > 0:
        s = match if (i > 0 and j > 0 and x[i-1] == y[j-1]) else mismatch
        if i > 0 and j > 0 and F[i][j] == F[i-1][j-1] + s:
            ax = x[i-1] + ax; ay = y[j-1] + ay; i -= 1; j -= 1
        elif i > 0 and F[i][j] == F[i-1][j] + gap:
            ax = x[i-1] + ax; ay = "-" + ay; i -= 1
        else:
            ax = "-" + ax; ay = y[j-1] + ay; j -= 1
    return ax, ay

우하단에서 시작해서 어느 방향에서 왔는지 되짚어가면서 매치·갭·불일치를 재구성한다.

Rosalind로 채점받기

Rosalind GLOB에 들어가서 데이터셋을 다운로드하고, 위 코드를 그대로 붙여 실행한 뒤 결과 정수 하나를 제출해보자. 5분 안에 정답을 확인할 수 있다. 자기가 짠 알고리즘이 실제로 통과하는 경험은 튜토리얼 100번 보는 것보다 강하다.

복잡도와 실무의 벽

Needleman-Wunsch는 **시간 O(mn) · 공간 O(mn)**이다. 두 서열 길이가 10,000씩만 되어도 10⁸개의 셀을 채워야 한다. 인간 게놈(3 × 10⁹)과 다른 게놈 서열을 이대로 비교하겠다는 건 미친 짓이다.

그래서 실무는 두 갈래로 갈라진다.

  1. 국소 정렬: 두 서열의 전체가 아니라 유사한 부분만 뽑아낸다. 다음 편에서 다룰 Smith-Waterman (M07).
  2. 휴리스틱: 격자를 안 채우고 시드(seed) + 확장 방식으로 근사해를 빠르게 찾는다. 세 편 뒤의 BLAST (M10).

BWA 같은 매핑 도구도 결국 Needleman-Wunsch 격자를 국소적으로만, 필요한 부분만 채우는 것이다. 즉, 앞으로 만나게 될 모든 정렬 도구는 이 편의 격자를 어떻게 잘라내고 접어 넣느냐의 문제다.

CS 매핑 — 이건 DP의 정수다

DryBench에서 다룬 DP 격자 최적 경로 문제를 기억하자. Needleman-Wunsch는 그것의 문자열 판이다.

  • 격자 좌표 = 두 서열의 부분 서열 조합
  • 경로 = 하나의 정렬
  • 점수 함수 = 각 셀의 승계 비용
  • 최적 부분 구조 = 격자를 채워가면서 이전 결과 재사용
  • Traceback = 최적 경로 복원

CS에서 편집 거리 (Levenshtein distance) 배웠다면 정확히 그거다. 편집 거리는 삽입·삭제·치환에 각각 +1을 매기고 최소를 구하는 반면, Needleman-Wunsch는 진화적 관점에서 매치를 보상하고 갭·불일치를 벌하는 방식으로 최대를 구한다. 부호만 뒤집으면 완전히 같은 알고리즘이다.

다음 편으로 이어지는 갈래

  • 다음 편 (M07): Smith-Waterman — 국소 정렬로 확장. 격자에서 음수를 0으로 자르는 순간 문제가 완전히 달라진다.
  • 두 편 뒤 (M08): 아핀 갭 페널티 — 갭을 여는 것과 이어가는 것에 다른 벌점을 줘야 실제 진화적 이벤트와 맞는다.
  • 세 편 뒤 (M09): BLOSUM/PAM — match=+1, mismatch=-1은 너무 거칠다. 아미노산의 진화 확률로 점수를 정하자.

더 깊게 파고 싶다면

본문은 BPD가 직접 재구성한 서술이다. 원리를 이미 이해했다면 아래 정통 강의로 심화해도 좋다.

  • MIT OCW 7.91J — Christopher Burge 교수의 Local Alignment (BLAST) and Statistics 강의 (자막 완비, 자동 번역 95%+). 격자 유도를 수학적으로 더 엄밀하게 다룬다.
  • Harvard STAT115 — Xiaole Shirley Liu 교수의 Early Sequence Alignment (1 with 1) (자막 완비, 자동 번역 우수). 판서 방식으로 격자를 하나씩 채워 보인다.
  • Stanford CS262 — Gill Bejerano 교수의 재생목록 (자막 완비). CS 관점에서의 최적화 경로 문제 해석이 강하다.
  • 참고 교재: Compau & Pevzner Bioinformatics Algorithms — Rosalind와 연동된 무료 문제집. 이 편의 파이썬 구현을 그대로 이어갈 수 있다.

Rosalind에서 파이썬으로 GLOB · LOCA · GAFF까지 풀어 오면 다음 편 진입 준비가 끝난다. 손을 움직여야 한다.