Selecting the Top K Reads — Using a Heap to Keep 100 Out of 1 Million

After completing this topic

By combining the heap (priority queue) and sorting concepts learned in the textbook, you can create a tool that selects and streams the top K highest-quality reads from a large sequencing dataset. Understand through code why a heap is dramatically faster than a full sort and why it is perfect for streaming data.

This article is a general educational example. Real-world read filtering involves much more complex pipelines, including quality trimming and adapter removal, but it accurately covers the core concept of top-K.

"Top 1,000 out of 1 Million" - The Pitfalls of a Naive Approach

You have 1 million sequencing reads, and you've calculated the average quality score for each. You want to keep only the top 1,000.

Naive Approach A: Sort Everything

def topk_by_sort(reads: list[dict], k: int) -> list[dict]:
    return sorted(reads, key=lambda r: -r["avg_quality"])[:k]

This approach is accurate, but it has a time complexity of O(n log n). For 1 million reads, this is approximately 2 × 10⁷ comparisons. In Python, this might take a few seconds. Crucially, all the data must be in memory.

Naive Approach B: Update Minimum Each Time

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

This has a time complexity of O(n × k). 1 million × 1,000 = 1 billion. This is much slower.

The real approach is to use a heap. Maintain a min-heap of size K. When new data arrives, compare it to the minimum value in the heap, and replace it if the new data is better. The time complexity is O(n log K). 1 million × log(1,000) ≈ 10 million. This is several times faster than sorting, and it is streamable.

From Black Box to Components

Component 1: The Min-Heap

A heap is a binary tree where the parent is always less than or equal to its children (min-heap). Python's heapq simulates a min-heap using a list.

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] — It's a list, but it satisfies the heap property

smallest = heapq.heappop(nums)  # 1

Key property: heappush and heappop are O(log n). Finding the minimum is O(1) (nums[0]).

Component 2: Top-K Ideas

Now, let's find the top-K using a min-heap of size 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)

Why does it work? The minimum value in the heap (heap[0]) is the Kth largest among the top-K. If a new value is greater than this, it needs to replace the Kth value – which is exactly what heapreplace does.

Time complexity analysis: Each of the n iterations takes O(log K). Total: O(n log K). If k=1000 and n=1,000,000, this is about 10⁷ operations. Faster than sorting O(n log n) ≈ 2 × 10⁷, but it uses only K memory.

Component 3: Maintaining Identity with Tuples

A read is a dictionary, not a single scalar. When storing it in the heap, store it as a (score, read) tuple.

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)]

Why do we need counter? When two read dictionaries have the same quality, comparing them directly will raise an exception. Tuples compare elements sequentially, so we put counter as the second element to use it as a tie-breaker for quality.

Stream Processing — The Real Power of Heaps

What if you can’t fit a million reads into memory at once? What if you have to read them one line at a time from a file? Heap access still works — this is a crucial difference from sorted access.

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)]

This function finds the top-K reads using memory for only K reads, regardless of whether the file size is 100GB or 1TB. It works the same way in a streaming environment where data flows in real time. This is the real-world killer application of heaps.

Use heapq.nlargest for a concise solution

The Python standard library already provides this pattern.

from heapq import nlargest

top_reads = nlargest(1000, reads_iter, key=lambda r: r["avg_quality"])

Internally, it uses exactly the algorithm described above. In practice, use this. However, it's important to understand why this syntax is so fast — this allows you to implement it yourself when a top-K solution is needed in other languages, and it allows for customization when K is large or there are special conditions.

Fading — Two Blanks You Need to Fill In

Blank 1: Conditional Top-K

Select the top-K reads from those with a quality above a certain threshold and a length above a certain threshold.

def topk_conditional(
    reads_iter,
    k: int,
    min_quality: float = 20.0,
    min_length: int = 100
) -> list[dict]:
    """
    Top-K selection, considering only reads that satisfy the conditions.
    Reads that do not meet the conditions are not added to the heap.
    """
    heap: list = []
    counter = 0
    
    for read in reads_iter:
        # TODO: Check quality/length conditions. If the conditions are not met, continue.
        # If the conditions are met, push or heapreplace into the heap.
        pass
    
    return [entry[2] for entry in sorted(heap, reverse=True)]

Hint: Start with if read["avg_quality"] < min_quality or len(read["seq"]) < min_length: continue.

Blank 2: Sorting by Multiple Criteria

If the qualities are the same, prioritize shorter reads (or vice versa). Incorporate this into the second element of the tuple.

def topk_by_multiple(
    reads_iter,
    k: int
) -> list[dict]:
    """
    Primary sorting: in descending order of avg_quality.
    Secondary sorting: in ascending order of length (if the quality is the same, prefer shorter reads).
    """
    heap: list = []
    counter = 0
    
    for read in reads_iter:
        score = read["avg_quality"]
        length_key = -len(read["seq"])
        counter += 1
        
        # TODO: Construct the tuple to be stored in the heap.
        # Priority: (larger score, smaller length)
        # Since it's a min-heap, use the score as is and reverse the sign of the length.
        pass
    
    return [entry[3] for entry in sorted(heap, reverse=True)]

Hint: entry = (score, length_key, counter, read). In the min-heap, the order of popping is: (lower score, if the scores are the same, lower length_key = larger length). Therefore, to get the top, sort in reverse.

Reflection — Differences from a Production-Ready Read Filter

Refinement of Quality Calculation: Your avg_quality is a simple average. In production, the quality of each position is considered individually—as quality degradation at the beginning and end is common, a position-specific filter is needed.

Adapter Trimming: The first step in a production filter is the removal of sequencing adapter sequences. Fastp and Trimmomatic are standard tools for this. The Top-K filter is a subsequent step.

Paired-End Processing: In production, the two reads must be processed together as a pair. If one read is filtered out, its mate must also be handled accordingly.

Memory Mapping of Files: For very large FASTQ files, memory mapping (mmap) is used for access, or the compressed file (gzip) is streamed and parsed. You can add gzip support by changing your open to gzip.open.

GPU Acceleration: For processing very large numbers of reads (billions), GPU toolchains like NVIDIA Parabricks are used. Although the heap itself is not inherently GPU-friendly, quality calculation can be parallelized.

Extension Project

1. Paired-End Support: Stream two FASTQ files simultaneously and apply the top-K filter while maintaining read pairing.

2. Automatic gzip Support: Check the file extension; if it's .gz, automatically use gzip.open.

3. Quality Distribution Dashboard: Compare the quality distribution before and after filtering using matplotlib.

4. Bottom-K Addition: Extract low-quality reads and save them separately (for troubleshooting).

Component Guide for This Section

• [F] Heap/Priority Queue: Maintain the top-K elements using a min-heap of size K. Push/replace operations take O(log K) time.
• [F] Comparison with Sorting: Compare sort (O(n log n)) with heap (O(n log K)). Determine when each approach is advantageous.
• [W] File I/O: FASTQ parsing (provided as a complete script).

[F] = You will implement this yourself / [W] = Provided as complete code.