서열 조립 — k-mer 그래프로 짧은 조각들을 이어붙이기
이 토픽을 마치면
교과서에서 배운 그래프 인접 리스트 와 DFS/BFS 순회 를 엮어서, 짧은 sequencing read 수백 개를 하나의 긴 contig로 이어붙이는 조립 도구를 직접 만들 수 있습니다. SPAdes나 Velvet 같은 실제 assembler가 마법처럼 해주는 그 조립의 내부 원리 — De Bruijn 그래프 — 를 코드로 이해하게 됩니다.
이 글은 교육용 일반 예제 입니다. 실제 서열 조립은 훨씬 복잡한 실무 파이프라인이지만, 그 심장부에 있는 알고리즘 개념을 정확히 다룹니다.
"왜 3B bp 유전체가 300bp read에서 나오지?" — 조립의 함정
여러분이 사람 유전체를 시퀀싱한다고 합시다. 현대 시퀀서(Illumina)는 한 번에 약 300bp짜리 짧은 read를 수억 개 뽑아냅니다. 30억 bp짜리 사람 유전체를 이 조각들로부터 다시 한 줄로 이어붙여야 합니다.
가장 순진한 접근은 read들을 뒤로 갈수록 겹치는 부분을 찾아 이어붙이는 것 입니다.
def naive_assemble(reads: list[str]) -> str: result = reads[0] for read in reads[1:]: overlap = find_overlap(result, read) result += read[overlap:] return result이 접근에 두 가지 치명적 문제가 있습니다.
문제 1: read 순서가 의미 없습니다. 시퀀서에서 나온 read들은 유전체상 무작위 위치에서 뽑힌 것이라, reads[0] 과 reads[1] 이 실제 유전체에서 이웃일 확률이 매우 낮습니다.
문제 2: 겹침 계산이 폭발합니다. 1억 개 read 중 어떤 두 개가 실제로 겹치는지 확인하려면 최악의 경우 (1억)² = 1경 번 비교해야 합니다. 우주가 끝나기 전에 끝나지 않습니다.
진짜 접근은 그래프를 씁니다. read들을 짧은 조각(k-mer)으로 잘게 잘라, 각 k-mer를 노드로 하고 서로 이어지는 관계를 엣지로 하는 De Bruijn 그래프를 만듭니다. 그다음 이 그래프의 경로를 찾으면 그것이 조립된 서열입니다.
블랙박스에서 부품으로 — De Bruijn 그래프 열어보기
이 프레임을 정확히 이해하려면 세 부품이 필요합니다.
부품 1: k-mer — 문자열을 잘게 자르기
k-mer 는 길이 k의 짧은 서열입니다. ATGCAT 을 k=3 (3-mer)로 자르면 ATG, TGC, GCA, CAT 네 개가 나옵니다.
def get_kmers(sequence: str, k: int) -> list[str]: return [sequence[i:i+k] for i in range(len(sequence) - k + 1)]핵심 특성은 연속된 k-mer가 (k-1)-mer prefix/suffix를 공유 한다는 것입니다. ATG 의 suffix TG = TGC 의 prefix. 이 공유가 그래프의 엣지가 됩니다.
부품 2: 그래프 인접 리스트
De Bruijn 그래프의 노드는 (k-1)-mer입니다. 엣지는 k-mer입니다. 한 노드에서 다른 노드로 가는 엣지가 있다는 것은, 그 두 (k-1)-mer가 어떤 k-mer의 prefix·suffix 관계로 이어진다는 뜻입니다.
이 그래프를 인접 리스트(dictionary)로 저장합니다.
from collections import defaultdict
def build_de_bruijn(reads: list[str], k: int) -> dict[str, list[str]]: graph = defaultdict(list) for read in reads: for kmer in get_kmers(read, k): prefix = kmer[:-1] suffix = kmer[1:] graph[prefix].append(suffix) return graphgraph["ATG"] 를 조회하면 ATG 노드에서 나가는 엣지들의 목적지 노드 리스트가 나옵니다.
부품 3: DFS/BFS 순회
그래프가 만들어지면, 모든 엣지를 정확히 한 번씩 지나는 경로(오일러 경로)를 찾아야 합니다. 이 경로를 이어붙이면 원래 서열이 재구성됩니다.
DFS로 접근합시다. 시작 노드에서 재귀적으로 이웃 노드를 방문하되, 엣지를 소비할 때마다 그래프에서 제거합니다.
def find_eulerian_path(graph: dict[str, list[str]], start: str) -> list[str]: stack = [start] path = [] graph = {k: list(v) for k, v in graph.items()} while stack: node = stack[-1] if graph.get(node): next_node = graph[node].pop() stack.append(next_node) else: path.append(stack.pop()) return path[::-1]이것이 Hierholzer 알고리즘입니다. 시간복잡도 O(엣지 수).
세 부품을 합쳐 — 조립 파이프라인
이제 세 부품을 하나의 파이프라인으로 엮어봅시다.
from collections import defaultdict
def kmerize(sequence: str, k: int) -> list[str]: return [sequence[i:i+k] for i in range(len(sequence) - k + 1)]
def build_de_bruijn_graph(reads: list[str], k: int) -> dict[str, list[str]]: graph = defaultdict(list) for read in reads: for kmer in kmerize(read, k): prefix, suffix = kmer[:-1], kmer[1:] graph[prefix].append(suffix) return dict(graph)
def find_start_node(graph: dict[str, list[str]]) -> str: out_degree = {node: len(edges) for node, edges in graph.items()} in_degree: dict[str, int] = defaultdict(int) for edges in graph.values(): for target in edges: in_degree[target] += 1 for node in graph: if out_degree[node] > in_degree[node]: return node return next(iter(graph))
def eulerian_path(graph: dict[str, list[str]], start: str) -> list[str]: graph = {k: list(v) for k, v in graph.items()} stack, path = [start], [] while stack: node = stack[-1] if node in graph and graph[node]: stack.append(graph[node].pop()) else: path.append(stack.pop()) return path[::-1]
def path_to_sequence(path: list[str]) -> str: if not path: return "" return path[0] + "".join(node[-1] for node in path[1:])
def assemble(reads: list[str], k: int = 5) -> str: graph = build_de_bruijn_graph(reads, k) start = find_start_node(graph) path = eulerian_path(graph, start) return path_to_sequence(path)간단한 예제로 확인해 봅시다.
original = "ATGCATGCATGATG"reads = [original[i:i+8] for i in range(0, len(original) - 7)]# ['ATGCATGC', 'TGCATGCA', 'GCATGCAT', 'CATGCATG', 'ATGCATGA', 'TGCATGAT', 'GCATGATG']
result = assemble(reads, k=5)print(result) # ATGCATGCATGATG (또는 유사 재구성)페이딩 — 여러분이 채워야 할 세 빈칸
이제 여러분 차례입니다. 위 코드에는 실전에서 반드시 다뤄야 하는 세 빈칸이 있습니다. 각각의 힌트만 드립니다.
빈칸 1: 시퀀싱 오류 필터링
실제 read에는 오류 문자가 섞여 있습니다. 오류가 있는 k-mer는 그래프에 극소수 발생 빈도의 잘못된 노드를 만듭니다. 이 노드가 조립을 오염시킵니다.
def filter_low_frequency_kmers(reads: list[str], k: int, threshold: int = 2) -> list[str]: from collections import Counter kmer_counts: Counter[str] = Counter() for read in reads: for kmer in kmerize(read, k): kmer_counts[kmer] += 1 # TODO: threshold 미만 발생한 k-mer를 포함하는 read를 필터링 filtered = [] for read in reads: # 여러분이 채울 곳: 이 read의 모든 k-mer의 빈도가 threshold 이상인지 확인 pass return filtered힌트: all(kmer_counts[kmer] >= threshold for kmer in kmerize(read, k))
빈칸 2: 여러 시작점의 처리
실제 시퀀싱 데이터에서는 유전체가 하나의 오일러 경로로 딱 이어지지 않을 수 있습니다. 여러 contig가 생성되어야 하는 경우가 대부분입니다.
def assemble_multi_contigs(reads: list[str], k: int) -> list[str]: graph = build_de_bruijn_graph(reads, k) contigs = [] while graph: # TODO: 남은 그래프에서 시작 노드를 골라 하나의 contig 조립 # 조립 후 사용된 엣지 제거 # 남은 그래프가 비어있지 않으면 반복 pass return contigs힌트: find_start_node 를 반복 호출하고, 각 반복마다 eulerian_path 로 하나씩 조립. 각 조립 후 사용된 엣지가 이미 pop 됐으므로 자연히 그래프가 줄어듭니다.
빈칸 3: 반복 서열(repeat) 대응
유전체에는 반복 서열이 흔합니다. 같은 짧은 서열이 유전체 여러 곳에 나타나면, De Bruijn 그래프에서 그 서열에 해당하는 노드가 여러 엣지의 교차점이 되어 조립이 모호해집니다.
이 문제의 완전한 해결은 실제 assembler가 여러 정교한 휴리스틱을 씁니다(paired-end read, long read 등). 여러분이 시도할 접근:
def detect_repeat_nodes(graph: dict[str, list[str]]) -> set[str]: """엣지가 여러 방향으로 나가는 노드는 반복 서열 후보""" # TODO: out-degree > 1 또는 in-degree > 1인 노드를 반환 pass이 반복 노드들이 어디에 있는지만 표시해도, 사용자에게 "여기 조립이 모호합니다" 라는 경고를 줄 수 있습니다.
성찰 — 이 코드가 실전 assembler와 어떻게 다른가
여러분이 만든 assembler는 개념적으로 SPAdes나 Velvet과 같은 뿌리를 공유하지만, 실제 프로덕션 assembler는 훨씬 정교합니다. 주요 차이를 정리합니다.
규모 대응: 실제 assembler는 수억 read를 처리합니다. 여러분의 파이썬 구현은 dictionary 오버헤드로 이 규모를 감당할 수 없습니다. 실제 assembler는 C++로 짜여 있고, k-mer를 정수로 인코딩(2비트/염기)해 메모리를 압축합니다.
paired-end read: 실제 시퀀서는 한 조각의 양 끝에서 두 read를 뽑아 그 사이 간격을 대략 알고 있습니다. 이 간격 정보가 반복 서열 문제 해결에 결정적입니다. 여러분의 구현엔 아직 없습니다.
long read (PacBio/Nanopore): 최근 시퀀서는 훨씬 긴 read(10kb+)를 만듭니다. 이 경우 접근이 완전히 달라져서(OLC — Overlap-Layout-Consensus), De Bruijn이 아닌 문자열 겹침 그래프를 씁니다.
품질 점수: 실제 read는 각 염기에 quality score가 붙어 있습니다. 이 정보로 오류 필터링을 훨씬 정교하게 할 수 있습니다.
확장 프로젝트
이 조립기를 몇 개 방향으로 확장해 볼 수 있습니다.
1. FASTQ 파일 입력: 실제 시퀀싱 데이터는 FASTQ 형식으로 옵니다. 이 파일을 파싱해 read + quality score를 함께 읽어오는 함수를 추가.
2. 통계 리포트: 조립 결과의 N50 (contig 길이 분포 지표), 총 조립 길이, contig 수 등을 계산해 표시.
3. 시각화: networkx + matplotlib 으로 작은 De Bruijn 그래프를 그려보기. 반복 노드가 어떻게 조립을 복잡하게 만드는지 눈으로 확인.
4. 실제 데이터로 시도: 세균 유전체(대장균은 약 4.6M bp)의 공개 시퀀싱 데이터를 다운받아 여러분의 assembler를 돌려보고, SPAdes 결과와 비교.
이 편의 부품 지도
이 응용편이 어떤 개별 CS 개념들을 조합했는지 정리합니다. 각 개념의 자세한 설명은 DryBench에 있습니다.
- [F] 그래프 인접 리스트:
graph[node] = [neighbor1, neighbor2, ...]dictionary 표현. De Bruijn 그래프 저장. - [F] DFS/BFS: Hierholzer의 스택 기반 DFS로 오일러 경로. 실제 조립 순회.
- [W] 해시 테이블·dictionary:
defaultdict(list)로 O(1) 노드 조회. - [W] k-mer 슬라이스: 문자열 슬라이싱과 리스트 컴프리헨션.
[F] = 여러분이 직접 구현 / [W] = 완성 코드로 제공한 도구 개념.