서열 정렬 — 동적 계획법으로 두 DNA를 정확히 겹치기
이 토픽을 마치면
교과서에서 배운 동적 계획법, 재귀-반복 대비, 2D 매트릭스 를 엮어서, 두 DNA 서열을 정확히 겹쳐 어디에 삽입·결실·치환이 있는지 찾아내는 정렬 도구를 직접 만들 수 있습니다. BLAST가 마법처럼 해주는 그 정렬의 내부 원리를 코드로 이해하게 됩니다.
이 글은 교육용 일반 예제 입니다. 변이 분석은 바이오인포매틱스에서 매일 하는 일이라 소재로 삼았습니다.
"왜 돌연변이 위치가 안 맞아?" — 순진한 비교의 함정
여러분이 BRCA1 유전자의 야생형 부분 서열과 환자에서 시퀀싱한 돌연변이 서열을 비교한다고 합시다.
야생형: ATGCTAGCATGCA
돌연변이: ATGCAGCATGCA가장 간단한 비교는 위치별로 하나씩 맞춰보는 것입니다.
def compare_naive(seq1: str, seq2: str) -> list[tuple[int, str, str]]: diffs = [] for i in range(min(len(seq1), len(seq2))): if seq1[i] != seq2[i]: diffs.append((i, seq1[i], seq2[i])) return diffs이걸 위 두 서열에 돌려보면 이렇게 나옵니다.
위치 4: T → A
위치 5: A → G
위치 6: G → C
위치 7: C → A
위치 8: A → T
위치 9: T → G
위치 10: G → C
위치 11: C → A이상하죠. 위치 4 이후로 거의 모든 위치가 다르게 나옵니다. 진짜 돌연변이가 이렇게 많이 있는 걸까요?
실제로 두 서열을 눈으로 보면 답이 보입니다. 야생형에서 위치 4의 T가 삭제(deletion) 된 것뿐입니다.
야생형: ATGC[T]AGCATGCA
돌연변이: ATGC AGCATGCA
↑ 여기서부터 오프셋이 하나씩 밀림한 글자 삭제가 이후 위치를 전부 밀어서 마치 여러 치환이 있는 것처럼 보인 것입니다. 순진한 비교는 이 삭제·삽입을 알아채지 못합니다.
이 문제를 해결하는 게 서열 정렬(sequence alignment) 입니다. 두 서열 사이에 필요한 만큼 갭(gap, - 문자) 을 넣어서 어긋난 부분을 보정하는 것입니다.
야생형: ATGCTAGCATGCA
돌연변이: ATGC-AGCATGCA ← 4번 위치에 gap이제 위치 5 이후가 완전히 맞아떨어집니다. 진짜 변이가 하나의 결실(deletion) 임이 명확해집니다.
문제는 이 gap을 어디에 넣을지 를 결정하는 것입니다. 두 서열이 100 글자쯤 되면 gap을 넣을 수 있는 위치 조합이 폭발적으로 많아집니다. 어떤 조합이 "가장 자연스러운" 정렬인지 어떻게 찾을까요? 이 글은 그 답을 밑바닥부터 만듭니다.
먼저 완성품을 봅시다 (블랙박스 먼저 돌리기)
우리가 만들 도구는 두 서열을 받아 최적 정렬 을 반환합니다.
aligned1, aligned2, score = align(seq1, seq2)
print(aligned1) # 'ATGCTAGCATGCA'print(aligned2) # 'ATGC-AGCATGCA'print(score) # 8 (예시)=== 정렬 결과 ===
Wild: ATGCTAGCATGCA
Mut: ATGC-AGCATGCA
Score: 8
=== 정렬 없이 순진하게 비교 ===
Wild: ATGCTAGCATGCA
Mut: ATGCAGCATGCA
불일치 위치: 8개 (거의 전부)
=== 정렬 후 실제 변이 ===
위치 4: T 결실 (deletion) — 진짜 변이 1개같은 두 서열인데 정렬 전후 해석이 완전히 달라집니다. 정렬은 변이 해석의 첫 단계이자 결정적 단계 입니다. 지금부터 이 정렬을 만드는 원리를 짚어갑니다.
이 도구는 어떤 부품으로 조립되는가 (부품 분해도)
서열 정렬기 (Aligner)
┌──────────────────────────────────────────────────┐
│ [입력] 두 서열 불러오기 ─── 부품: 문자열/파일 입출력 │ ← 완성 제공 (도구)
│ │ │
│ ▼ │
│ [1단계] 점수 매트릭스 채우기 │
│ 부품: 2D 배열 + 동적 계획법 │ ← 직접 만든다 ★
│ │ │
│ ▼ │
│ [2단계] 역추적으로 최적 경로 복원 │
│ 부품: 재귀·반복 선택 │ ← 직접 만든다 ★
│ │ │
│ ▼ │
│ [출력] 정렬된 두 서열 + 점수 │
└──────────────────────────────────────────────────┘| 부품 | 어디서 배웠나 | 이 도구에서 하는 일 |
|---|---|---|
| 문자열 입출력 | string-and-file-io | 서열 읽고 다루기 |
| 2D 배열 | matrix-2d-array | 하위 문제 답을 저장하는 격자 |
| 동적 계획법 | dynamic-programming | 하위 문제 답을 재사용해서 폭발 회피 |
| 재귀 vs 반복 | recursion-vs-iteration | 격자를 채울 때, 결과를 되짚을 때 방식 선택 |
📌 이 개념들 처음 본다면 (상단 진입 링크)
직접 만드는 새 개념은 동적 계획법 · 2D 매트릭스 · 재귀/반복 선택 3개. 서열 다루기는 도구라 완성본으로 드립니다. 딱 3개죠 — 인지 한계선 안입니다.
만들기 1단계 — 서열 준비 (완성 제공)
먼저 정렬할 서열 두 개를 준비합니다. 실제로는 FASTA 파일에서 읽지만, 브라우저 실습을 위해 문자열로 정의합니다.
# 야생형 BRCA1 부분 서열 (교육용 예시)wild = "ATGCTAGCATGCA"
# 돌연변이 서열 (4번 위치 T 결실)mut = "ATGCAGCATGCA"
assert len(wild) == 13assert len(mut) == 12assert all(c in "ACGT" for c in wild)assert all(c in "ACGT" for c in mut)이 두 서열을 갭 - 을 적절히 넣어서 겹치도록 만드는 게 우리의 목표입니다.
만들기 2단계 — 점수 규칙 정하기 (완성 제공)
정렬의 "품질"을 어떻게 잴 것인가부터 정해야 합니다. 표준적인 방법은 점수 시스템 입니다.
MATCH = 1 # 같은 염기끼리 정렬되면 +1MISMATCH = -1 # 다른 염기끼리 정렬되면 -1GAP = -2 # 갭 하나가 들어가면 -2
def score_pair(a: str, b: str) -> int: """두 염기(또는 갭)의 정렬 점수.""" if a == b: return MATCH return MISMATCH
assert score_pair("A", "A") == 1assert score_pair("A", "T") == -1이 점수 규칙이 정렬 알고리즘의 취향 을 결정합니다.
- gap penalty가 낮으면 (예: -0.5) → 갭이 많이 들어간 정렬 선호
- gap penalty가 높으면 (예: -5) → 갭보다 미스매치 감수
바이오 실무에서는 BLOSUM/PAM 같은 정교한 점수 행렬 을 씁니다. 여기서는 학습을 위해 단순 규칙을 씁니다.
만들기 3단계 — 재귀부터 시작해 폭발을 겪어보기 ★ (재귀 vs 반복)
✍️ 직접 채우는 구간. 부품 = 재귀. 먼저 재귀로 정렬 점수를 구해보고, 왜 실용적이지 않은지 실감한다.
정렬 점수 문제를 재귀로 표현하면 놀랍도록 단순합니다. 두 서열 끝에서 시작해서, 각 위치에서 세 가지 선택 을 재귀로 시도합니다.
- 두 서열의 마지막 글자끼리 정렬 →
score_pair(s1[-1], s2[-1])+align(s1[:-1], s2[:-1]) - s1의 마지막 글자와 gap 정렬 →
GAP+align(s1[:-1], s2) - s2의 마지막 글자와 gap 정렬 →
GAP+align(s1, s2[:-1])
세 옵션 중 가장 높은 점수 를 택합니다. 이걸 재귀로 쓰면:
def align_recursive(s1: str, s2: str) -> int: """최적 정렬의 점수만 재귀로 계산 (역추적 X).""" if not s1: return len(s2) * GAP if not s2: return len(s1) * GAP return max( align_recursive(s1[:-1], s2[:-1]) + score_pair(s1[-1], s2[-1]), align_recursive(s1[:-1], s2) + GAP, align_recursive(s1, s2[:-1]) + GAP, )
# 검증: 짧은 예제로만assert align_recursive("AA", "AA") == 2 # 두 매치assert align_recursive("AT", "A") == 1 + GAP # 매치 + 갭 = 1 - 2 = -1이게 정답을 주긴 합니다. 그런데 긴 서열에는 절대 못 씁니다. 다음 코드를 돌려보면 이유가 실감납니다.
import time
# 길이 15 정도만 되어도 재귀 지옥short1 = "ATGCTAGCATGCAAG"short2 = "ATGCAGCATGCAAG"
t0 = time.time()score = align_recursive(short1, short2)elapsed = time.time() - t0print(f"길이 15 재귀: {elapsed:.2f}초")# 몇 초 걸립니다. 길이 20이면 몇 분, 25면 몇 시간이유는 재귀가 똑같은 하위 문제를 수억 번 반복 계산 하기 때문입니다. align_recursive("ATGCT", "ATGC") 값이 재귀 트리 여러 가지에서 반복해서 계산됩니다. 매번 처음부터.
이게 O(3^n) 시간 복잡도의 위력입니다. 이 시점에서 동적 계획법 이 등장합니다.
🤔 자기설명 프롬프트
align_recursive("AT", "A")를 재귀 트리로 그려보세요. 몇 개의 서브 문제가 반복해서 등장하나요? (길이 3, 4 정도로만 확장해도 겹침이 눈에 띕니다.)
만들기 4단계 — 격자에 답을 저장하기 ★ (2D 매트릭스 + 동적 계획법)
✍️ 직접 채우는 구간. 부품 = 2D 배열 + 동적 계획법. 하위 문제 답을 격자에 저장해서 반복 계산 회피.
발상은 이렇습니다. align(s1[:i], s2[:j])의 답을 격자 dp[i][j] 에 저장합니다. 그러면 같은 하위 문제를 두 번 계산할 필요가 없습니다.
🔎 동적 계획법이란 (드로어 — dynamic-programming) 큰 문제의 답이 작은 문제들의 답으로 만들어질 때, 작은 문제의 답을 저장해두면 두 번 계산 안 해도 됩니다. 재귀의 지수 폭발이 저장(memoization) 한 번으로 다항 시간이 됩니다. 격자에 채우는 방식은 특히 "Bottom-Up 동적 계획법" 이라 부릅니다.
격자의 크기는 (len(s1)+1) × (len(s2)+1). 여분의 한 행/열은 빈 문자열과의 정렬 을 위한 기저 조건입니다.
def build_score_matrix(s1: str, s2: str) -> list[list[int]]: """dp[i][j] = s1[:i]와 s2[:j]의 최적 정렬 점수.""" n, m = len(s1), len(s2) dp = [[0] * (m + 1) for _ in range(n + 1)]
# 기저: 빈 문자열과 정렬하면 갭만 있음 for i in range(1, n + 1): dp[i][0] = dp[i-1][0] + GAP for j in range(1, m + 1): dp[0][j] = dp[0][j-1] + GAP
# 격자 채우기 (Bottom-Up) for i in range(1, n + 1): for j in range(1, m + 1): match = dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1]) delete = dp[i-1][j] + GAP # s1 글자와 gap insert = dp[i][j-1] + GAP # gap과 s2 글자 dp[i][j] = max(match, delete, insert) return dp
dp = build_score_matrix(wild, mut)
# 검증: 최종 셀이 재귀 답과 일치해야 함assert dp[len(wild)][len(mut)] == align_recursive(wild, mut)# 크기 검증assert len(dp) == len(wild) + 1assert len(dp[0]) == len(mut) + 1# 왼쪽 열은 순수 gap 누적assert dp[3][0] == 3 * GAP같은 답을 얻는데 속도는 완전히 다릅니다.
import time
# 앞서 재귀가 몇 초 걸렸던 길이 15t0 = time.time()score = build_score_matrix(short1, short2)[len(short1)][len(short2)]elapsed = time.time() - t0print(f"길이 15 동적계획법: {elapsed*1000:.2f}ms")# 밀리초 수준재귀에서 몇 초 걸리던 게 밀리초 안에 끝납니다. 이게 O(n×m) 의 힘입니다. 격자 셀 하나하나가 단 한 번씩만 계산됩니다.
🔎 왜 격자가 답인가 (드로어 — matrix-2d-array) 두 인덱스(i, j)에 의존하는 하위 문제는 자연스럽게 2D 격자 로 표현됩니다. 각 셀이 하나의 하위 문제, 셀 사이의 화살표가 의존 관계입니다. 이 구조가 눈에 보이면 dp 문제 접근이 쉬워집니다.
🤔 자기설명 프롬프트
dp[i][j]를 계산할 때 오직dp[i-1][j-1],dp[i-1][j],dp[i][j-1]만 필요합니다. 격자 전체를 저장할 필요가 있을까요? 실은 직전 한 줄만 있으면 다음 줄을 채울 수 있어서 O(min(n,m)) 공간으로도 점수만은 구할 수 있습니다. 언제 그 최적화가 안 될까요? (힌트: 역추적)
만들기 5단계 — 격자를 되짚어 정렬 복원하기 ★
✍️ 직접 채우는 구간. 부품 = 반복적 역추적. 격자의 마지막 셀에서 시작해 어떻게 그 값이 만들어졌는지 거꾸로 따라간다.
점수만 아는 것과 어떻게 정렬됐는지 정확히 아는 것 은 다릅니다. 실제 정렬(어디에 갭이 들어가는지)을 복원하려면 격자를 끝에서 시작으로 역추적 해야 합니다.
각 셀 dp[i][j]에서 우리가 세 후보 중 어느 걸 골랐는지를 되짚습니다.
dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1])를 골랐다면 → 대각선 이동 (매치/미스매치)dp[i-1][j] + GAP를 골랐다면 → 위쪽 이동 (s1 글자와 gap)dp[i][j-1] + GAP를 골랐다면 → 왼쪽 이동 (gap과 s2 글자)
def traceback(dp: list[list[int]], s1: str, s2: str) -> tuple[str, str]: aligned1: list[str] = [] aligned2: list[str] = [] i, j = len(s1), len(s2)
while i > 0 or j > 0: current = dp[i][j] if i > 0 and j > 0 and current == dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1]): aligned1.append(s1[i-1]) aligned2.append(s2[j-1]) i -= 1 j -= 1 elif i > 0 and current == dp[i-1][j] + GAP: aligned1.append(s1[i-1]) aligned2.append("-") i -= 1 else: aligned1.append("-") aligned2.append(s2[j-1]) j -= 1
return "".join(reversed(aligned1)), "".join(reversed(aligned2))
a1, a2 = traceback(dp, wild, mut)
# 검증: 길이 동일 · 갭 제거 시 원본 서열 복원assert len(a1) == len(a2)assert a1.replace("-", "") == wildassert a2.replace("-", "") == mut# 우리가 넣은 결실 지점이 실제로 발견되어야 함assert "-" in a2 # 돌연변이 쪽에 gap이 있어야 함이 함수는 반복(iteration) 으로 구현했지만, 같은 논리를 재귀로도 쓸 수 있습니다. 격자 크기가 매우 클 때는 반복이 스택 오버플로우 방지에 유리하고, 코드 가독성은 재귀가 더 나을 때도 있습니다. 선택은 문제 크기에 달렸습니다.
🔎 재귀와 반복은 언제 어느 것 (드로어 — recursion-vs-iteration) 재귀는 문제를 자연스럽게 표현할 때 우아합니다. 반복은 스택 안 쌓아서 안전하고 종종 빠릅니다. 작은 문제 · 우아함 필요 → 재귀 / 큰 문제 · 안정성 필요 → 반복. 이 도구의 3단계(격자 채우기)는 반복으로, 4단계(역추적)도 반복으로 짰습니다. 두 곳 다 재귀로 짤 수 있지만 서열이 길어지면 스택 위험이 있어서 반복을 택합니다.
🤔 자기설명 프롬프트
traceback에서 우리는 세 후보 중 하나 를 선택합니다. 만약 두 후보의 점수가 같다면(동점) 어느 걸 골라야 할까요? 실제 정렬 결과가 달라질 수 있을까요? (힌트: 여러 최적 정렬이 존재할 수 있습니다.)
부품을 하나로 — 완성된 정렬기 클래스
두 함수를 하나의 도구로 묶습니다.
class Aligner: def __init__(self, match=1, mismatch=-1, gap=-2): self.match = match self.mismatch = mismatch self.gap = gap
def _score(self, a: str, b: str) -> int: return self.match if a == b else self.mismatch
def align(self, s1: str, s2: str) -> tuple[str, str, int]: n, m = len(s1), len(s2) dp = [[0] * (m + 1) for _ in range(n + 1)] for i in range(1, n + 1): dp[i][0] = dp[i-1][0] + self.gap for j in range(1, m + 1): dp[0][j] = dp[0][j-1] + self.gap for i in range(1, n + 1): for j in range(1, m + 1): dp[i][j] = max( dp[i-1][j-1] + self._score(s1[i-1], s2[j-1]), dp[i-1][j] + self.gap, dp[i][j-1] + self.gap, ) # 역추적 a1: list[str] = [] a2: list[str] = [] i, j = n, m while i > 0 or j > 0: if i > 0 and j > 0 and dp[i][j] == dp[i-1][j-1] + self._score(s1[i-1], s2[j-1]): a1.append(s1[i-1]); a2.append(s2[j-1]); i -= 1; j -= 1 elif i > 0 and dp[i][j] == dp[i-1][j] + self.gap: a1.append(s1[i-1]); a2.append("-"); i -= 1 else: a1.append("-"); a2.append(s2[j-1]); j -= 1 return "".join(reversed(a1)), "".join(reversed(a2)), dp[n][m]
aligner = Aligner()a1, a2, score = aligner.align(wild, mut)print(a1)print(a2)print(f"score = {score}")
# 검증: 클래스 결과가 함수 결과와 일치dp = build_score_matrix(wild, mut)assert score == dp[len(wild)][len(mut)]이 도구가 오늘 우리가 이야기하는 Needleman-Wunsch 알고리즘 의 축소판입니다. 실무 도구(BLAST, EMBOSS의 needle 등)는 여기에 정교한 점수 행렬, 갭 오프닝/확장 별도 점수, 로컬 정렬 옵션 등을 추가한 것뿐입니다. 뼈대는 여러분이 방금 만든 것과 같습니다.
성능 딥다이브 — 왜 이렇게 빨라지나
import time
def bench(n: int, m: int): s1 = "AT" * (n // 2) s2 = "AG" * (m // 2) t0 = time.time() build_score_matrix(s1, s2) return time.time() - t0
for size in [10, 50, 100, 200]: t = bench(size, size) print(f"길이 {size:4d}: {t*1000:.2f}ms ({size*size} 셀)")assert bench(100, 100) < bench(200, 200)숫자는 기계마다 다르지만 방향은 언제나 같습니다. 셀 수(n×m)에 정확히 비례 합니다. 재귀 버전이 3^n으로 폭발할 때, dp 버전은 n×m으로 자란다는 결정적 차이가 있습니다.
🔎 Big-O로 정리 (드로어 — big-o-notation)
- 재귀 순진 버전: O(3^n) — 지수. 길이 30에서 이미 실용 불가
- 동적 계획법 버전: O(n·m) — 다항. 길이 10,000끼리도 몇 초
- 저장(격자) 하나가 지수를 다항으로 만들었습니다. 이것이 동적 계획법의 마법입니다.
다른 길도 있다 (멀티패스 성찰)
- 글로벌 vs 로컬 정렬: 우리가 만든 것은 전체 서열을 처음부터 끝까지 맞추는 Needleman-Wunsch (글로벌). 실무에서는 한 서열 안에서 다른 서열의 조각을 찾는 Smith-Waterman (로컬)을 자주 씁니다. 알고리즘 뼈대는 거의 같고, 점수 시스템과 초기 조건만 다릅니다. 언제 무엇: 진화 계통 서열끼리 = 글로벌 / 유전자 안에서 도메인 찾기 = 로컬.
- BLAST 계열의 heuristic: 우리 dp는 두 서열 길이 n·m에 비례하는데, 유전체(수십억 bp) 대상으로는 여전히 큽니다. BLAST는 정확 해를 포기하고 근사 해를 빠르게 얻는 heuristic (seed → extend). 트레이드오프: 100% 최적 vs 실용 속도.
- 점수 행렬 정교화: 우리는 match/mismatch 이진 규칙이지만 실무는 BLOSUM62/PAM250 같이 잔기 쌍마다 다른 점수를 매기는 행렬을 씁니다. 단백질 서열의 화학적 유사성을 반영합니다.
- Affine gap penalty: 우리는 갭 하나에 -2 고정이지만, 실제로는 갭을 처음 여는 비용이 크고 이어지는 갭은 작은 이중 규칙(예: 오프닝 -10, 확장 -1)을 씁니다. 짧은 갭 여러 개보다 긴 갭 하나가 더 자연스럽다는 진화적 관찰의 반영입니다.
핵심: "가능한 하위 문제를 격자에 저장해 반복 계산 회피." 이 원칙은 서열 정렬을 넘어 경로 찾기, 편집 거리, 지식 검색 등 무수한 곳에서 재사용됩니다. 여러분이 방금 만든 뼈대가 열어주는 문이 매우 넓습니다.
다음 단계로 (하단 출구 링크)
- 격자를 채울 때 왜 이 순서 로 셀을 방문하는지 → 의존 관계와 방문 순서
- 서열이 매우 길 때 메모리 최적화 → 1D dp 압축 기법
- 앞서 만든 인덱스와 조합 → 응용편 서열 DB 인덱싱
직접 해보기 (독립 문제)
- 삽입 시나리오: 야생형
ATGCATGCA와 돌연변이ATGCXATGCA(X = 추가 염기 삽입)를 정렬하고 삽입 위치가 정확히 잡히는지 assert로 확인하세요. - 점수 튜닝:
Aligner(gap=-5)로 하면 결과가 어떻게 달라지나요? gap=-0.5 라면? gap 페널티가 정렬의 "성격"을 어떻게 결정하는지 정리하세요. - 다중 최적 정렬: 두 개의 최적 정렬이 존재하는 서열 쌍을 만들고,
traceback을 확장해 모든 최적 정렬을 반환 하는 버전을 짜세요. - 도전 — 메모리 O(min(n,m)): 점수만 필요하다면 격자 전체 대신 직전 한 줄만 유지해도 됩니다. 이 최적화 버전을 구현하세요. (역추적을 포기해야 하는 대신, 매우 긴 서열에 씁니다.)
정리
우리는 "두 서열의 최적 정렬을 찾는" 문제를, 세 개의 부품으로 정복했습니다.
- 동적 계획법 이 지수 폭발(O(3^n))을 다항 시간(O(n·m))으로 바꿨습니다.
- 2D 매트릭스 가 하위 문제 답을 저장하는 격자가 됐습니다.
- 재귀와 반복의 선택 이 격자 채우기와 역추적 각각에서 안정성을 확보했습니다.
BLAST가 마법처럼 정렬 결과를 던져줄 때, 이제 그 안이 보일 겁니다. 마법이 아니라, 교과서에서 배운 dp 격자에 서열을 얹은 것뿐입니다.
이 글은 일반 교육 예제 입니다. 실무 서열 정렬 도구(BLAST/EMBOSS/Bowtie 등)는 여기에 정교한 점수 행렬, 갭 이중 페널티, seed-and-extend heuristic, 병렬화 등이 추가된 것입니다. 그 상세한 버전은 이 뼈대 위에 여러분이 붙이거나, 검증된 도구에 맡기면 됩니다.