BioPlayground

🧬
목록으로

FM-index와 초고속 패턴 매칭: BWA · Bowtie2의 심장

왜 BWA가 인간 게놈 3GB에 3GB 인덱스로 초당 수십만 리드를 정렬하는가. FM-index의 Rank/Select · Backward search 원리를 손 계산으로 유도합니다.

심화
|
18분
|
검증 완료 (2026-07-19)
FM-indexBWABowtiebackward search
진행률0/34 (0%)

왜 이 편이 필요한가

M14에서 BWT를 다뤘습니다. 압축과 검색의 이중 성질을 가진 자료구조라고 했습니다. 그런데 정확히 어떻게 검색이 빠른가는 아직 정리하지 않았습니다.

Paolo Ferragina와 Giovanni Manzini가 2000년에 발표한 FM-index가 그 답입니다. BWT에 Rank 배열과 Count 배열을 붙이면 임의 패턴 P의 매치 위치를 O(|P|) 시간에 찾을 수 있습니다. 인간 게놈 30억 bp에서 100 bp 패턴을 찾는 데 100번의 배열 조회면 됩니다.

이 편이 정착되면 BWA가 왜 그렇게 빠른지, 왜 인간 게놈 인덱스가 원본과 비슷한 크기로 저장되는지 명확해집니다.

FM-index의 세 가지 부품

FM-index는 세 가지로 구성됩니다.

  1. BWT: 원본 문자열의 Burrows-Wheeler 변환.
  2. C 배열: 각 문자에 대해, 정렬된 문자열에서 그 문자보다 사전순으로 작은 문자의 총 개수. 즉 첫 열에서 그 문자가 시작하는 위치.
  3. Occ (Rank) 배열: Occ(c, i) = BWT의 앞 i개 문자 중 c가 나타난 횟수.

T = BANANA$의 예로 봅시다.

  • BWT = ANNB$AA
  • 정렬된 첫 열 = $AAABNN

C 배열:

text
C['`'] = 0  # `는 사전순 첫 번째
C['A'] = 1  # A 이전: $ 하나
C['B'] = 4  # B 이전: $, A, A, A (4개)
C['N'] = 5  # N 이전: $, A, A, A, B (5개)

Occ 배열 (BWT = ANNB$AA):

text
i:      0 1 2 3 4 5 6 7
BWT:    A N N B $ A A
Occ(A): 1 1 1 1 1 2 3 3
Occ(N): 0 1 2 2 2 2 2 2
Occ(B): 0 0 0 1 1 1 1 1
Occ($): 0 0 0 0 1 1 1 1

Backward Search — 패턴을 뒤에서부터 검색

FM-index의 검색은 패턴의 뒤에서부터 진행합니다. 이걸 backward search라고 부릅니다.

패턴 P = ANA를 BANANA에서 찾아봅시다.

초기: 마지막 문자 A

A가 BWT에서 나타나는 첫 열 범위를 찾습니다. C 배열에서 A는 1부터 시작하고, A의 총 개수는 3이므로 범위는 [1, 4).

  • start = C['A'] = 1
  • end = C['A'] + count('A') = 1 + 3 = 4

이 범위는 정렬된 회전들 중 A로 시작하는 것들을 가리킵니다.

text
1: A$BANAN
2: ANA$BAN
3: ANANA$B

다음: NA (뒤에서 두 번째 문자 N)

이제 범위 안의 회전들 중 바로 앞에 N이 오는 것들만 남깁니다. LF-mapping을 사용합니다.

new_start=C[N]+Occ(N,start)=5+Occ(N,1)=5+1=6\text{new\_start} = C[N] + \text{Occ}(N, \text{start}) = 5 + \text{Occ}(N, 1) = 5 + 1 = 6 new_end=C[N]+Occ(N,end)=5+Occ(N,4)=5+2=7\text{new\_end} = C[N] + \text{Occ}(N, \text{end}) = 5 + \text{Occ}(N, 4) = 5 + 2 = 7

새 범위 [6, 7). 즉 정렬된 회전 6번 하나만 남았습니다.

text
6: NA$BANA

마지막: ANA (뒤에서 세 번째 문자 A)

new_start=C[A]+Occ(A,6)=1+2=3\text{new\_start} = C[A] + \text{Occ}(A, 6) = 1 + 2 = 3 new_end=C[A]+Occ(A,7)=1+3=4\text{new\_end} = C[A] + \text{Occ}(A, 7) = 1 + 3 = 4

새 범위 [3, 4). 정렬된 회전 3번 하나.

text
3: ANANA$B

여기서 SA[3] = 1이므로, 원본에서 ANA는 위치 1에서 시작합니다.

파이썬으로 완전 구현

python
from collections import Counter


def build_fm_index(text: str):
    text = text + "$"
    n = len(text)
    sa = sorted(range(n), key=lambda i: text[i:])
    bwt = "".join(text[(sa[i] - 1) % n] for i in range(n))


    # C 배열
    sorted_chars = sorted(text)
    C = {}
    for i, c in enumerate(sorted_chars):
        if c not in C:
            C[c] = i


    # Occ 배열
    alphabet = sorted(set(text))
    Occ = {c: [0] * (n + 1) for c in alphabet}
    for i, c in enumerate(bwt):
        for ch in alphabet:
            Occ[ch][i+1] = Occ[ch][i]
        Occ[c][i+1] += 1


    return {"bwt": bwt, "C": C, "Occ": Occ, "sa": sa, "n": n}


def fm_search(index, pattern: str) -> list[int]:
    C, Occ, sa, n = index["C"], index["Occ"], index["sa"], index["n"]
    start, end = 0, n
    for c in reversed(pattern):
        if c not in C:
            return []
        start = C[c] + Occ[c][start]
        end   = C[c] + Occ[c][end]
        if start >= end:
            return []
    return sorted(sa[start:end])


# 사용
index = build_fm_index("BANANA")
print(fm_search(index, "ANA"))  # [1, 3]

패턴의 길이만큼(3번) 반복하고, 각 반복에서 배열 조회 몇 번만 하면 매치 위치를 얻습니다. 총 시간 O(|P|) — 참조 크기와 무관.

실제 BWA의 최적화

위 구현은 교육용이라 Occ 배열이 원본 크기의 수 배가 됩니다. 실제 BWA는 다음 최적화들을 씁니다.

  • Wavelet Tree: Occ를 압축된 형태로 저장 (알파벳 4글자 DNA에 특화).
  • Sample SA: 접미사 배열을 전부 저장하지 않고, 일정 간격으로만 저장. 매치 위치는 LF-mapping으로 인접 샘플까지 이동해 복원.
  • Checkpoint: Occ 배열도 일정 간격 체크포인트만 저장.

결과: 인간 게놈 3GB의 BWA 인덱스가 대략 3GB (원본과 비슷한 크기)입니다. 매우 인상적인 압축비입니다.

실무 응용 — 짧은 리드 정렬

BWA (2009), Bowtie (2009), Bowtie2 (2012) 모두 FM-index를 씁니다. 100bp Illumina 리드 하나의 정렬 시간이 초당 수십만 개 규모입니다. 인간 게놈 30× 커버리지의 시퀀싱 결과(약 3억 리드)를 몇 시간 만에 정렬합니다.

  • BWA-backtrack (BWA 1세대): 짧은 리드(70bp 미만)용.
  • BWA-SW: 중간 길이(70bp～10kb)용. Smith-Waterman 확장 결합.
  • BWA-MEM (실무 표준): 짧은 리드부터 중간 롱리드까지. FM-index로 seed를 찾고 아핀 갭 확장.

S03 편에서 BWA-MEM의 실무 사용법을 다룹니다.

자주 만나는 실무 함정

  • 인덱스 파일: BWA는 여러 파일(.amb, .ann, .bwt, .pac, .sa)로 인덱스를 저장합니다. 모두 있어야 하며, 참조 게놈 파일과 같은 디렉토리에 두는 것이 관례.
  • 인덱스 재사용: 참조 게놈이 갱신되지 않는 한 인덱스는 계속 재사용합니다. bwa index 명령은 한 번만.
  • 버전 호환: BWA 인덱스는 대개 호환되지만, BWT 압축 방식이 바뀌면 재인덱싱이 필요할 수 있습니다.
  • 경계 케이스: 반복 서열이 많은 게놈은 seed 매치가 매우 많아 검색이 느려질 수 있습니다. BWA-MEM은 -c 옵션으로 최대 후보를 제한합니다.

CS 매핑

  • Rank/Select 자료구조: 정보 이론과 자료구조 이론의 정통. 압축된 표현에서 임의 위치의 문자 · 개수를 빠르게 접근.
  • Wavelet Tree: 알파벳이 큰 문자열에 대한 Rank/Select 확장. Grossi · Gupta · Vitter가 2003년에 발표.
  • 압축된 인덱스: 원본과 인덱스의 크기 합이 원본 두 배 이하가 되는 것이 정보 이론적으로 흥미로운 성취.
  • 백워드 서치: LF-mapping의 반복 응용. 자료구조의 우아함이 알고리즘 설계의 핵심.

다음 편으로 이어지는 갈래

  • 다음 편 (M16): minimap2 · 미니마이저 — 롱리드용 인덱스 사고. FM-index와 다른 접근.
  • 두 편 뒤 (M17): DIAMOND — 단백질 검색용 인덱스 최적화.
  • 여섯 편 뒤 (S03): BWA-MEM 실무 — 지금 배운 자료구조의 실제 사용.
  • 열 편 뒤 (S07～S12): GATK4 — BWA로 정렬된 결과 위에서 변이 검출.

더 깊게 파고 싶다면

본문은 BPD가 자체 재구성한 서술입니다. 심화는 아래로.

  • CMU 02-510 — Ben Langmead 교수(Bowtie 저자)의 FM-index and Read Alignment (자막 완비, 자동 번역 우수). Bowtie 저자가 직접 설명하는 정통 강의.
  • 원 논문: Ferragina, P. & Manzini, G. (2000), Opportunistic data structures with applications, FOCS. FM-index 시조.
  • 원 논문: Li, H. & Durbin, R. (2009), Fast and accurate short read alignment with Burrows-Wheeler transform, Bioinformatics 25, 1754–1760. BWA 시조.
  • 원 논문: Langmead, B. et al. (2009), Ultrafast and memory-efficient alignment of short DNA sequences to the human genome, Genome Biology 10, R25. Bowtie 시조.
  • 참고 교재: Mäkinen et al. Genome-Scale Algorithms Chapter 8. 정통 교과서.

Rosalind에서 관련 문자열 검색 문제를 풀어봅시다. 다음 편 M16의 minimap2가 왜 FM-index를 안 쓰는지 자연스럽게 이해됩니다.