왜 이 편을 먼저 다뤄야 하는가
서열 정렬은 바이오인포매틱스의 알파벳입니다. 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년 지난 시점에 바이오에 들어온 셈입니다.
점수 함수 · 격자 · 점화식
점수는 어떻게 매기나
두 서열을 붙일 때 세 가지 이벤트가 일어납니다.
- 일치 (match): 같은 문자끼리 대응 → +s_match (예: +1)
- 불일치 (mismatch): 다른 문자끼리 대응 → -s_mismatch (예: -1)
- 갭 (gap): 한쪽에 빈 자리 → -d (예: -2)
예를 들어 아래와 같은 정렬이 있다고 해봅시다.
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)를 이렇게 정의합니다.
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 정도로 해봅시다.
| − | G | A | C | |
|---|---|---|---|---|
| − | 0 | −2 | −4 | −6 |
| G | −2 | +1 | −1 | −3 |
| A | −4 | −1 | +2 | 0 |
| T | −6 | −3 | 0 | +1 |
F(1,1):GvsG— 매치 →F(0,0)+1 = 1. 대각선 승.F(2,2):AvsA— 매치 →F(1,1)+1 = 2.F(3,3):TvsC— 불일치 →max(F(2,2)−1, F(2,3)−2, F(3,2)−2) = max(1, −2, −2) = 1.
우하단 값 +1이 최적 점수입니다. 여기서부터 화살표를 거꾸로 따라가면 (traceback) 실제 정렬이 튀어나옵니다.
G A T
G A C매치 2번 + 불일치 1번 = 2 − 1 = 1. 격자 값과 일치합니다.
파이썬 20줄로 완전 구현
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 함수를 추가합니다.
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⁹)과 다른 게놈 서열을 이대로 비교하겠다는 건 현실적으로 불가능합니다.
그래서 실무는 두 갈래로 갈라집니다.
- 국소 정렬: 두 서열의 전체가 아니라 유사한 부분만 뽑아냅니다. 다음 편에서 다룰 Smith-Waterman (M07)이 그 갈래입니다.
- 휴리스틱: 격자를 안 채우고 시드(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까지 풀어 오면 다음 편 진입 준비가 끝납니다. 손을 움직여야 합니다.