上位K個のreadを抽出 — 優先キューを使って100万個の中から100個だけ残す
このトピックを終えると
教科書で学んだヒープ（優先度キュー）とソートを組み合わせて、大量のシーケンスリードから品質が上位のK個だけをストリーミング方式で抽出するツールを自分で作成できます。ヒープが完全なソートよりもはるかに高速である理由、そしてなぜストリーミングデータに最適なのかをコードを通して理解します。
この記事は教育用の一般的な例です。実際のリードフィルタリングは、品質トリミング、アダプターの除去など、はるかに複雑なパイプラインですが、その中心となるトップKの概念を正確に扱います。
「100万個中上位1000個」— naiveなアプローチの落とし穴
100万個のシーケンスリードがあり、それぞれの平均クオリティスコアを計算しました。この中から上位1000個だけ残したいとします。
naiveなアプローチ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⁷回の比較が必要です。Pythonで数秒かかります。そして、決定的に、すべてのデータがメモリに収まる必要があります。
naiveなアプローチ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万 × 1000 = 10⁹。はるかに遅いです。
実際のアプローチは、ヒープを使用することです。**サイズKの最小ヒープ（min-heap）**を維持します。新しいデータが来たら、ヒープの最小値と比較して、より良い値であれば置き換えます。時間計算量はO(n log K)です。100万 × log(1000) ≈ 10⁷。ソートよりもはるかに高速で、ストリーミングが可能です。
ブラックボックスからコンポーネントへ
コンポーネント 1: 最小ヒープ (min-heap)
ヒープは、親ノードが常に子ノード以下である（最小ヒープ）二分木です。Python の heapq は、リストを使って最小ヒープをシミュレートします。
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 の最小ヒープを使って上位 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 が必要なのか？ 2 つの read の品質が同じ場合、read dict が互いに比較され、例外が発生します。タプルは最初の要素から順に比較するため、counter を 2 番目の要素として追加し、品質が同じ場合の tie-break として使用します。
ストリーミング処理 — 優先キューの真の力
100万件のリードが一度にメモリに収まらない場合、またはファイルから1行ずつ読み込む必要がある場合はどうでしょうか？ 優先キューへのアクセスは引き続き機能します。これは、ソートによるアクセスとは根本的に異なります。
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個のリード分のメモリでトップKを求めます。 また、データがリアルタイムで流れるストリーミング環境でも、そのまま使用できます。これが優先キューの、実用的な強力な応用例です。
heapq.nlargest を使って簡潔に書く
Python の標準ライブラリには、このパターンを実装する関数が既に用意されています。
from heapq import nlargest
top_reads = nlargest(1000, reads_iter, key=lambda r: r["avg_quality"])内部的には、まさに上記のようなアルゴリズムを使用しています。実用的な場面では、これを使用します。ただし、この構文がなぜ高速なのかを理解しておく必要があります。他の言語で同様の処理が必要になった場合、自分で実装したり、K が大きい場合や特別な条件がある場合にカスタマイズしたりすることができます。
フェーディング — あなたが埋めるべき2つの空白
空白1：条件付きトップK
品質が特定の閾値以上で、かつ長さも特定の閾値以上のreadの中からトップKを選び出す。
def topk_conditional( reads_iter, k: int, min_quality: float = 20.0, min_length: int = 100) -> list[dict]: """ 条件を満たすreadの中からのみトップKを選択する。 条件を満たさないreadはヒープに含めない。 """ heap: list = [] counter = 0 for read in reads_iter: # TODO: 品質/長さの条件を確認し、条件を満たさない場合は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：複数の基準によるソート
品質が同じであれば、短いreadを優先する（または逆）。これをタプルの2番目の要素として格納する。
def topk_by_multiple( reads_iter, k: int) -> list[dict]: """ 1次ソート：avg_qualityの高い順 2次ソート：長さが短い順（同じ品質の場合は短いreadを優先） """ 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の長い順）。したがって、トップにするためには、逆順にソートする。
考察 — 実際のリードフィルタとの違い
品質計算の高度化: あなたの avg_quality は単純な平均です。実際のフィルタは、各位置の品質を個別に評価します。リードの前部と後部で品質が低下することが多いため、位置ごとのフィルタが個別に必要です。
アダプタートリミング: 実際のフィルタの最初のステップは、シーケンシングアダプター配列の除去です。fastp、trimmomaticなどが標準的です。Top-Kフィルタはその後のステップです。
ペアエンド処理: 実際のフィルタでは、2つのリードをペアとして一緒に処理する必要があります。片方のリードがフィルタで除外された場合、ペアも一緒に管理する必要があります。
メモリマッピングファイル: 非常に大きなFASTQファイルは、mmapでアクセスするか、圧縮された状態で（gzip）ストリーミング解析します。あなたの open を gzip.open に置き換えると、gzipがサポートされます。
GPUによる高速化: 非常に大きなリード（数十億個）の処理には、NVIDIA ParabricksなどのGPUツールチェーンを使用します。ヒープ自体はGPUに最適化されているわけではありませんが、品質計算は並列化されます。
拡張プロジェクト
1. ペアエンドサポート: 2つのFASTQファイルを同時にストリーミングし、ペアを維持したままトップKフィルターを適用する。
2. gzip自動サポート: ファイル拡張子を確認し、.gzの場合、gzip.openを自動的に使用する。
3. 品質分布ダッシュボード: フィルタリング前後の品質分布をmatplotlibで比較する。
4. ボトムKの追加: 低品質のリードを抽出し、別途保存する（問題診断用）。
このセクションの構成要素
- [F] 優先度付きキュー (ヒープ)：サイズ K の最小ヒープとして上位 K 個を保持。O(log K) でのプッシュ/置換が可能。
- [F] ソートとの比較：ソート O(n log n) とヒープ O(n log K)。どちらが有利か。
- [W] ファイル I/O：FASTQ ファイルの解析（完成したスクリプトとして提供）。
[F] = 自分で実装 / [W] = 完成したコードとして提供。