Order Assembly — Stitching Together Short Fragments Using a k-mer Graph

Upon Completion of This Topic

By combining the graph adjacency list and DFS/BFS traversal concepts learned in the textbook, you will be able to create your own assembly tool that stitches together hundreds of short sequencing reads into a single long contig. You will gain a coded understanding of the underlying principle of assembly – the De Bruijn graph – which real-world assemblers like SPAdes or Velvet seem to perform magically.

This article is an educational, general example. Real sequence assembly involves much more complex practical pipelines, but it accurately covers the core algorithmic concepts.

"Why can't we get the 3B bp genome from 300 bp reads?" — The Assembly Trap

Let's say you are sequencing a human genome. Modern sequencers (Illumina) produce hundreds of millions of short reads, each about 300 bp long. You need to reassemble these fragments into a single sequence of the 3 billion bp human genome.

The most naive approach is to find overlapping regions between the reads and stitch them together.

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

There are two fatal problems with this approach.

Problem 1: The order of the reads doesn't matter. The reads are randomly sampled from the genome, so the probability that reads[0] and reads[1] are adjacent in the actual genome is very low.

Problem 2: Overlap calculation explodes. To check which two of the 100 million reads actually overlap, in the worst case, you have to make (100 million)² = 10^16 comparisons. It won't finish before the end of the universe.

The real approach uses a graph. The reads are cut into short fragments (k-mers), and a De Bruijn graph is created, where each k-mer is a node and the relationships between them are edges. Then, finding a path in this graph gives you the assembled sequence.

From Black Box to Components — Unraveling the De Bruijn Graph

To fully understand this framework, we need to break it down into three components.

Component 1: k-mers — Slicing the Sequence

A k-mer is a short sequence of length k. If we slice ATGCAT into k=3 (3-mers), we get ATG, TGC, GCA, and CAT.

def get_kmers(sequence: str, k: int) -> list[str]:
    return [sequence[i:i+k] for i in range(len(sequence) - k + 1)]

The key property is that consecutive k-mers share a (k-1)-mer prefix/suffix. The suffix of ATG, which is TG, is equal to the prefix of TGC. This shared portion becomes the edge of our graph.

Component 2: Graph Adjacency List

In a De Bruijn graph, the nodes are (k-1)-mers, and the edges are k-mers. An edge from one node to another indicates that the two (k-1)-mers are connected as a prefix-suffix within some k-mer.

We will store this graph as an adjacency list (a 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 graph

Querying graph["ATG"] will return a list of destination nodes, representing the edges going out from the ATG node.

Component 3: DFS/BFS Traversal

Once the graph is built, we need to find a path that traverses each edge exactly once (an Eulerian path). By concatenating this path, we can reconstruct the original sequence.

Let's use DFS. Starting from a node, we recursively visit neighboring nodes, removing the edge from the graph as we traverse it.

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]

This is the Hierholzer algorithm, with a time complexity of O(number of edges).

Combining the Three Components: The Assembly Pipeline

Now, let's connect the three components into a single pipeline.

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)

Let's test it with a simple example.

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 (or a similar reconstruction)

Fading — Three Blank Spaces for You to Fill

Now it's your turn. The code above has three blank spaces that you must fill in to make it work properly in a real-world scenario. Here are the hints for each:

Blank Space 1: Filtering Sequencing Errors

Real reads contain error characters. Erroneous k-mers create extremely low-frequency, spurious nodes in the graph. These nodes contaminate the assembly.

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: Filter reads that contain k-mers with frequency below the threshold
    filtered = []
    for read in reads:
        # Your turn: Check if all k-mers in this read have a frequency greater than or equal to the threshold
        pass
    return filtered

Hint: all(kmer_counts[kmer] >= threshold for kmer in kmerize(read, k))

Blank Space 2: Handling Multiple Starting Points

In real sequencing data, the genome may not be a single, perfect Euler path. There are often cases where multiple contigs need to be generated.

def assemble_multi_contigs(reads: list[str], k: int) -> list[str]:
    graph = build_de_bruijn_graph(reads, k)
    contigs = []
    
    while graph:
        # TODO: Choose a starting node from the remaining graph and assemble one contig
        # Remove used edges after assembly
        # Repeat if the remaining graph is not empty
        pass
    
    return contigs

Hint: Repeatedly call find_start_node, and assemble one contig per iteration using eulerian_path. After each assembly, the used edges are popped, so the graph naturally decreases.

Blank Space 3: Dealing with Repeats

Genomes often contain repeats. When the same short sequence appears in multiple places in the genome, the corresponding node in the De Bruijn graph becomes an intersection of multiple edges, making the assembly ambiguous.

A full solution to this problem involves real assemblers using several sophisticated heuristics (paired-end reads, long reads, etc.). An approach you can try:

def detect_repeat_nodes(graph: dict[str, list[str]]) -> set[str]:
    """Nodes with edges going in multiple directions are candidate repeat nodes"""
    # TODO: Return nodes with out-degree > 1 or in-degree > 1
    pass

Even just marking where these repeat nodes are, you can give the user a warning: "The assembly is ambiguous here."

Reflection — How does this code differ from a production assembler?

The assembler you've created shares conceptual roots with tools like SPAdes or Velvet, but a production assembler is significantly more sophisticated. Here's a summary of the key differences.

Scale: Production assemblers handle billions of reads. Your Python implementation cannot handle this scale due to dictionary overhead. Production assemblers are written in C++ and encode k-mers as integers (2 bits/base) to compress memory.

Paired-end reads: Modern sequencers generate paired-end reads, meaning they extract two reads from each end of a fragment, knowing the approximate distance between them. This distance information is crucial for resolving repeats, and is not yet implemented in your assembler.

Long reads (PacBio/Nanopore): Recent sequencers generate much longer reads (10kb+). In this case, the approach changes completely to OLC (Overlap-Layout-Consensus), which uses a string overlap graph instead of a De Bruijn graph.

Quality scores: Real reads have quality scores associated with each base. This information allows for much more sophisticated error filtering.

Extension Projects

This assembler can be extended in several directions.

1. FASTQ File Input: Real sequencing data comes in FASTQ format. Add a function that parses this file and reads read + quality score together.

2. Statistics Report: Calculate and display statistics about the assembly result, such as N50 (a contig length distribution metric), total assembly length, and the number of contigs.

3. Visualization: Use networkx + matplotlib to visualize a small De Bruijn graph. Observe how repetitive nodes complicate the assembly.

4. Try with Real Data: Download publicly available sequencing data for a bacterial genome (E. coli is approximately 4.6M bp), run your assembler on it, and compare the results with those from SPAdes.

Component Map for This Exercise

This section outlines the individual CS concepts combined in this exercise. Detailed explanations for each concept can be found in DryBench.

• [F] Graph Adjacency List: Dictionary representation: graph[node] = [neighbor1, neighbor2, ...]. Used to store the De Bruijn graph.
• [F] DFS/BFS: Hierholzer's stack-based DFS is used to find the Eulerian path, which represents the actual assembly traversal.
• [W] Hash Table / Dictionary: defaultdict(list) is used for O(1) node lookups.
• [W] k-mer Slicing: String slicing and list comprehension.

[F] = Concepts that you implement yourself / [W] = Tool concepts provided as complete code.