상위 K개 read 뽑기 — 힙으로 100만 개 중 100개만 남기기
이 토픽을 마치면
교과서에서 배운 힙(우선순위 큐) 과 정렬 을 엮어서, 대량 sequencing read에서 quality 상위 K개만 스트리밍 방식으로 골라내는 도구를 직접 만들 수 있습니다. 왜 힙이 전체 정렬보다 극적으로 빠른지, 왜 스트리밍 데이터에 완벽한지 코드로 이해합니다.
이 글은 교육용 일반 예제 입니다. 실전 read filtering은 quality trimming, adapter 제거 등 훨씬 복잡한 파이프라인이지만, 그 심장부의 top-K 개념을 정확히 다룹니다.
"100만 개 중 상위 1,000개" — 순진한 접근의 함정
여러분에게 100만 개의 sequencing read가 있고, 각각의 평균 quality score를 계산했습니다. 이 중 상위 1,000개 만 남기고 싶습니다.
순진한 접근 A: 전체 정렬
def topk_by_sort(reads: list[dict], k: int) -> list[dict]: return sorted(reads, key=lambda r: -r["avg_quality"])[:k]이 접근은 정확하지만 시간복잡도 O(n log n). 100만 개면 약 2 × 10⁷ 비교. 파이썬에서 몇 초. 그리고 결정적으로 모든 데이터가 메모리에 있어야 합니다.
순진한 접근 B: 매번 최소값 갱신
def topk_by_scan(reads: list[dict], k: int) -> list[dict]: result = [] for read in reads: result.append(read) if len(result) > k: worst = min(range(len(result)), key=lambda i: result[i]["avg_quality"]) result.pop(worst) return result이건 O(n × k). 100만 × 1,000 = 10⁹. 훨씬 느립니다.
진짜 접근은 힙입니다. 크기 K짜리 최소 힙(min-heap) 을 유지합니다. 새 데이터가 오면 힙의 최소값과 비교해 더 좋으면 교체합니다. 시간복잡도 O(n log K). 100만 × log(1,000) ≈ 10⁷. 정렬보다 몇 배 빠르고, 스트리밍이 가능합니다.
블랙박스에서 부품으로
부품 1: 최소 힙 (min-heap)
힙은 부모가 항상 자식보다 작거나 같은(min-heap) 이진 트리입니다. 파이썬 heapq 는 리스트로 min-heap을 시뮬레이션합니다.
import heapq
nums = []heapq.heappush(nums, 5)heapq.heappush(nums, 3)heapq.heappush(nums, 8)heapq.heappush(nums, 1)
print(nums) # [1, 3, 8, 5] — 리스트 형태지만 힙 성질 만족
smallest = heapq.heappop(nums) # 1핵심 성질: heappush 와 heappop 이 O(log n). 최소값 조회는 O(1) (nums[0]).
부품 2: 상위 K개 아이디어
이제 크기 K짜리 min-heap으로 top-K를 구합니다.
def topk_streaming(scores: list[float], k: int) -> list[float]: heap: list[float] = [] for score in scores: if len(heap) < k: heapq.heappush(heap, score) else: if score > heap[0]: heapq.heapreplace(heap, score) return sorted(heap, reverse=True)왜 작동하는가? 힙의 최소값(heap[0])은 현재 상위 K개 중 K번째 값입니다. 새 값이 이것보다 크면 K번째 값을 밀어내야 합니다 — 정확히 heapreplace 가 하는 일입니다.
시간 복잡도 분석: n번의 순회 각각 O(log K). 총 O(n log K). k=1000, n=1,000,000이면 약 10⁷ 연산. 정렬 O(n log n) ≈ 2 × 10⁷보다 약간 빠르지만, 메모리는 K만큼만 씁니다.
부품 3: 튜플로 아이덴티티 유지
read는 단일 스칼라가 아닌 dict입니다. 힙에 저장할 때는 (score, read) 튜플 형태로 저장합니다.
def topk_reads(reads_iter, k: int) -> list[dict]: heap: list[tuple[float, int, dict]] = [] counter = 0 for read in reads_iter: score = read["avg_quality"] counter += 1 entry = (score, counter, read) if len(heap) < k: heapq.heappush(heap, entry) else: if score > heap[0][0]: heapq.heapreplace(heap, entry) return [entry[2] for entry in sorted(heap, reverse=True)]왜 counter가 필요한가? 두 read의 quality가 같을 때 read dict가 서로 비교되어 예외가 발생합니다. 튜플은 첫 요소부터 순차 비교하므로, counter 를 두 번째로 넣어 quality tie-break로 사용합니다.
스트리밍 처리 — 힙의 진짜 힘
100만 개 read가 한꺼번에 메모리에 없다면? 파일에서 한 줄씩 읽어야 한다면? 힙 접근은 그대로 작동합니다 — 정렬 접근과 결정적으로 다른 지점입니다.
def topk_from_fastq(fastq_path: str, k: int) -> list[dict]: heap: list[tuple[float, int, dict]] = [] counter = 0 with open(fastq_path) as f: while True: header = f.readline().strip() if not header: break seq = f.readline().strip() plus = f.readline().strip() qual = f.readline().strip() avg_quality = sum(ord(c) - 33 for c in qual) / len(qual) counter += 1 read = {"header": header, "seq": seq, "qual": qual, "avg_quality": avg_quality} entry = (avg_quality, counter, read) if len(heap) < k: heapq.heappush(heap, entry) elif avg_quality > heap[0][0]: heapq.heapreplace(heap, entry) return [entry[2] for entry in sorted(heap, reverse=True)]이 함수는 파일 크기가 100GB든 1TB든 K개 read만큼의 메모리로 top-K를 구합니다. 데이터가 실시간으로 흐르는 스트리밍 환경에서도 그대로 씁니다. 이것이 힙의 실전 킬러 응용입니다.
heapq.nlargest 로 짧게 쓰기
파이썬 표준 라이브러리는 이미 이 패턴을 제공합니다.
from heapq import nlargest
top_reads = nlargest(1000, reads_iter, key=lambda r: r["avg_quality"])내부적으로 정확히 위와 같은 알고리즘을 씁니다. 실전에서는 이걸 씁니다. 다만 이 문법이 왜 이렇게 빠른지 이해하고 있어야 합니다 — 다른 언어의 top-K 필요 시 직접 구현할 수 있고, K가 크거나 특별한 조건이 있을 때 커스터마이즈할 수 있습니다.
페이딩 — 여러분이 채워야 할 두 빈칸
빈칸 1: 조건부 top-K
quality가 특정 임계 이상이면서 길이도 특정 이상인 read 중에서 top-K를 뽑습니다.
def topk_conditional( reads_iter, k: int, min_quality: float = 20.0, min_length: int = 100) -> list[dict]: """ 조건을 만족하는 read 중에서만 top-K. 조건 미달은 힙에 넣지 않음. """ heap: list = [] counter = 0 for read in reads_iter: # TODO: quality/length 조건 확인 후 조건 미달이면 continue # 조건 만족이면 힙에 push 또는 heapreplace pass return [entry[2] for entry in sorted(heap, reverse=True)]힌트: if read["avg_quality"] < min_quality or len(read["seq"]) < min_length: continue 를 먼저.
빈칸 2: 여러 기준 정렬
quality가 같으면 짧은 read를 우선 (또는 반대). 이걸 튜플 두 번째 요소로 넣습니다.
def topk_by_multiple( reads_iter, k: int) -> list[dict]: """ 1차 정렬: avg_quality 높은 순 2차 정렬: 길이 짧은 순 (같은 quality라면 짧은 것 선호) """ heap: list = [] counter = 0 for read in reads_iter: score = read["avg_quality"] length_key = -len(read["seq"]) counter += 1 # TODO: 힙에 저장할 튜플 구성 # 우선순위: (score 큰 순, length 작은 순) # min-heap이므로 score를 그대로, length는 부호 반대로 pass return [entry[3] for entry in sorted(heap, reverse=True)]힌트: entry = (score, length_key, counter, read). min-heap에서 pop되는 순서는 score 낮은 순, 같으면 length_key 낮은 순(= length 큰 순). 그래서 top이 되려면 반대로 정렬.
성찰 — 실전 read filter와의 차이
Quality 계산의 정교화: 여러분의 avg_quality 는 단순 평균. 실전은 각 위치의 quality를 개별적으로 봅니다 — 앞·뒤 부분의 quality 저하가 흔하므로, 위치별 필터가 별개로 필요.
Adapter trimming: 실전 필터의 첫 단계는 sequencing adapter 서열 제거. fastp, trimmomatic 등이 표준. Top-K 필터는 그 이후 단계.
Paired-end 처리: 실전은 두 read가 짝(mate)으로 함께 처리되어야 합니다. 한 쪽이 필터에서 탈락하면 짝도 함께 관리해야 합니다.
메모리 매핑 파일: 초대형 FASTQ는 mmap으로 접근하거나 압축된 상태(gzip)를 스트리밍 파싱합니다. 여러분의 open 을 gzip.open 으로 바꾸면 gzip 지원.
GPU 가속: 초대형 read 처리(수십억 개)는 NVIDIA Parabricks 같은 GPU 툴체인을 씁니다. 힙 자체가 GPU 친화적이진 않지만, quality 계산은 병렬화됩니다.
확장 프로젝트
1. paired-end 지원: 두 FASTQ 파일을 동시에 스트리밍하며 짝을 유지한 채 top-K 필터.
2. gzip 자동 지원: 파일 확장자를 확인해 .gz 면 gzip.open 자동 사용.
3. quality 분포 대시보드: 필터링 전/후의 quality 분포를 matplotlib으로 비교.
4. bottom-K 추가: 낮은 품질 read를 뽑아 별도 저장 (문제 진단용).
이 편의 부품 지도
- [F] 힙 · 우선순위 큐: 크기 K의 min-heap으로 top-K 유지. O(log K) push/replace.
- [F] 정렬 대비: sort O(n log n) vs heap O(n log K). 언제 어느 것이 유리한지.
- [W] 파일 I/O: FASTQ 파싱 (완성 스크립트로 제공).
[F] = 여러분이 직접 구현 / [W] = 완성 코드로 제공.