Computational Notes

Eli Niktab

Recent posts

Pruning the Pruners: Miniasm vs Hifiasm, the Graph Theory Behind Them, and How an Edge-Conditioned GIN Learned Both

2026-05-27

miniasm and hifiasm are both string-graph assemblers, and hifiasm’s graph-cleaning code is directly forked from miniasm’s — they share the same asg_t data structure and the same core function names. The differences live in what hifiasm adds on top: a parallel hierarchy of diploid-aware and trio-aware variants that consult per-read haplotype tags, multi-round error correction wrapped around the cleaning, and a richer multi-output graph (r_utg, p_utg, p_ctg, a_ctg) instead of a single GFA.

Here I walk through the graph theory underneath both tools, then describe what we learned trying to replace those heuristics with a graph neural network. Our current architecture — GINEConv (5 layers, 512 hidden dim) with focal loss — sits at 97.44 % validation edge accuracy on the real HG002 HiFi graph produced by hifiasm (~321 k edges, ~250 connected components, ~93 % / 7 % class imbalance).

1. Graph theory: the string graph and what “pruning” means

1.1 The bidirected string graph

A long-read assembly graph is a bidirected multigraph \(G = (V, E)\) where:

  • Each read \(r\) contributes two nodes, \(r^{+}\) and \(r^{-}\) (forward and reverse complement). Concretely, both miniasm and hifiasm use \(n_\text{vtx} = 2N\) for \(N\) reads, with the encoding v = read_id * 2 + dir.
  • Each pairwise overlap induces a pair of directed arcs: an arc \(u \rightarrow v\) has a twin arc \(v^{c} \rightarrow u^{c}\) where \(w^{c}\) flips the strand bit (v ^ 1). If one arc is deleted, the cleaner enforces \(v^{c} \rightarrow u^{c}\) is deleted too, as a post-cleanup pass in every del routine.

Adjacency is naturally written as a sparse matrix \(M \in \{0,1\}^{2N \times 2N}\) with \(M_{uv} = 1\) iff the arc \(u \rightarrow v\) exists.

1.2 The four topological motifs that “cleaning” targets

Motif What it looks like Why it’s spurious
Transitive edge \(u \to v\), \(v \to w\), and \(u \to w\) all present; lengths satisfy \(\ell(u{\to}v) + \ell(v{\to}w) \approx \ell(u{\to}w)\) Encodes the same overlap chain twice; the direct arc is redundant under Myers’s string-graph definition
Tip A short path of length \(\le k\) that dead-ends (no outgoing arcs at one end) Most commonly caused by sequencing errors near read ends or by chimeric reads
Bubble Two distinct paths \(u \to \dots \to w\) that diverge and reconverge within distance \(d\) Either real heterozygosity (preserve!) or local error (collapse) — disambiguating is the diploid problem
Short / weak overlap An arc whose overlap length is much smaller than the longest competing arc out of the same node Symptom of repeat-mediated false adjacency

1.3 Transitive reduction, formalised

Definition (Myers 2005). Given a directed graph \(G\) with edge lengths \(\ell\) and a tolerance \(\epsilon\), the arc \(u \to w\) is transitive at tolerance \(\epsilon\) if there is a \(v \ne u, w\) with \(u \to v\) and \(v \to w\) present and

\[ \bigl|\, \ell(u\!\to\!v) + \ell(v\!\to\!w) \;-\; \ell(u\!\to\!w) \,\bigr| \;\le\; \epsilon. \tag{1}\]

The transitive reduction \(G_r\) is obtained by removing every transitive arc. The implementation in miniasm (and its byte-identical fork in hifiasm) computes this in expected near-linear time per node by sorting outgoing arcs by length and using a length budget

\[ L = \ell\bigl(av_{n_v-1}\bigr) + \texttt{fuzz} \tag{2}\]

, then sweeping: an arc to \(w\) gets marked transitive iff some \(w\)-out-arc reaches a sibling within the budget. After every pass the assembler calls asg_cleanup() and asg_symm() to symmetrise the deletion across the twin arcs.

In sparse-matrix form the boolean version is

\[ M_r \;=\; M \;\setminus\; \Bigl( M \cap \bigcup_{k=2}^{\infty} M^{(k)} \Bigr), \qquad M^{(k)} = M^{(k-1)} \,\odot\, M, \tag{3}\]

where \(\odot\) is boolean matrix multiplication. Without a length budget this drops every topologically transitive arc — which, as we discovered the hard way in Phase 1 below, shatters repetitive regions in real genomes: shortcuts through repeat copies close a triangle topologically but the arc lengths do not satisfy Myers’s \(\epsilon\)-tolerance, so deleting the longer direct arc is biologically meaningless.

1.4 The tropical (min-plus) semiring

To enforce Myers’s fuzz tolerance as algebra rather than as a control-flow predicate, replace \((\times, +)\) with the min-plus tropical semiring \(\mathbb{T} = (\mathbb{R} \cup \{+\infty\}, \min, +)\):

\[ (M \otimes M)_{uw} \;=\; \min_{v}\!\bigl(\, M_{uv} + M_{vw} \,\bigr), \tag{4}\]

with \(M_{uv} = \ell(u\!\to\!v)\) if the arc exists and \(+\infty\) otherwise. An arc \(u \to w\) is transitive at tolerance \(\epsilon\) iff

\[ (M \otimes M)_{uw} \;\le\; \ell(u\!\to\!w) + \epsilon. \tag{5}\]

This is the exact algebraic form of the inner loop in asg_arc_del_trans. We implemented it with scatter_reduce_(..., reduce='amin') over the COO representation.

1.5 The Weisfeiler–Lehman test, message passing, and GIN

All modern GNNs fit Gilmer et al.’s (2017) message-passing template — at each step \(t\) a node updates its embedding from its previous state and an aggregation over its neighbours:

\[ \mathbf{h}_u^{(t+1)} \;=\; \phi\Bigl(\mathbf{h}_u^{(t)},\;\; \bigoplus_{v \in \mathcal{N}(u)} \psi\bigl(\mathbf{h}_v^{(t)}, \mathbf{e}_{uv}\bigr)\Bigr), \tag{6}\]

with \(\oplus\) a permutation-invariant pool. Different architectures instantiate \(\phi\), \(\psi\), \(\oplus\) differently: GCN (Kipf & Welling 2017) uses a symmetrically-normalised mean and ignores edge features; GraphSAGE (Hamilton et al. 2017) concatenates a separate node vector with a pooled neighbour vector; GAT (Veličković et al. 2018) and its corrected variant GATv2 (Brody et al. 2022) replace \(\oplus\) with attention.

Xu et al. (2019) proved that the discriminative power of any message-passing GNN is upper-bounded by the 1-dimensional Weisfeiler–Lehman (1-WL) graph-isomorphism test, and that GIN (Graph Isomorphism Network) attains that bound when \(\oplus\) is sum and \(\phi\) is a multi-layer MLP:

\[ \mathbf{h}_u^{(t+1)} \;=\; \mathrm{MLP}\Bigl((1+\eta)\,\mathbf{h}_u^{(t)} \;+\; \sum_{v \in \mathcal{N}(u)} \mathbf{h}_v^{(t)}\Bigr). \tag{7}\]

The GINE variant (Hu et al. 2020, “Strategies for pre-training graph neural networks”) generalises GIN to incorporate edge features into every aggregation step:

\[ \mathbf{h}_u^{(t+1)} \;=\; \mathrm{MLP}\Bigl((1+\eta)\,\mathbf{h}_u^{(t)} \;+\; \sum_{v \in \mathcal{N}(u)} \mathrm{ReLU}\bigl(\mathbf{h}_v^{(t)} + \mathbf{e}_{uv}\bigr)\Bigr). \tag{8}\]

For string-graph cleaning the edge features carry the overlap geometry (length, identity, span, etc.); GINE’s \(\mathrm{ReLU}(\mathbf{h}_v + \mathbf{e}_{uv})\) injection at every layer — rather than only at the final readout, as in vanilla SAGE — is what lets the network distinguish transitive arcs from tips from bubbles in a single forward pass. GAT-family attention uses edge features only as inputs to attention coefficients, which in our setting (Wave 2b in §4.2 below) was strictly less stable than direct injection into the aggregated signal. Strictly more expressive frameworks exist — Corso et al.’s PNA (Principal Neighbourhood Aggregation, 2020) replaces the sum with a learned combination of multiple aggregators, and graph transformers (Dwivedi & Bresson 2021) extend the WL bound at higher cost — but we did not need them.

2. Miniasm’s cleaning pipeline

Step 0 — Containment prefilter (optional, -R flag). ma_hit_no_cont(...). Contained reads are dropped on a coordinate-overlap criterion (max_hang, int_frac) before any graph is built.

Steps 1–3 — Overlap loading and two-pass read selection. ma_hit_read, then a 1-pass crude selection (ma_hit_sub → ma_hit_cut → ma_hit_flt), then a 2-pass fine selection (ma_hit_sub → ma_hit_cut → ma_hit_contained). ma_hit_contained() removes contained reads automatically (not gated by -R).

Step 4 — Graph cleaning. The string graph sg is built (ma_sg_gen), then cleaned in five ordered sub-stages:

Sub-stage Routine(s) Parameter
4.1 Transitive reduction asg_arc_del_trans(sg, opt.gap_fuzz) -g (fuzz, in bp)
4.2 Initial tip + bubble asg_cut_tip(sg, max_ext) then asg_pop_bubble(sg, bub_dist) -e, -d
4.3 Iterative short-overlap removal for i in 0..n_rounds: asg_arc_del_short(sg, r_i) and tip+bubble after each -n rounds, -r min,max ramped drop ratio
4.4 Internal-sequence and bi-loop cleanup asg_cut_internal, asg_cut_biloop, then tip + bubble again -e
4.5 Aggressive final short-overlap pass asg_arc_del_short(sg, final_ovlp_drop_ratio) + tip + bubble -F

So miniasm’s short-overlap removal is iterative — Stage 4.3 is a for loop over n_rounds+1 rounds, with the drop ratio ramped from min to max and tip + bubble passes after each round. The widely-repeated claim that “miniasm is single-pass” is therefore wrong in spirit. What miniasm does not do is wrap the entire cleaning cascade around an outer base-level error-correction loop — that outer loop is what hifiasm adds.

Algorithmically, the five graph-operation primitives:

  • asg_arc_del_trans — Myers TR with explicit fuzz, ops on the 2N node space, see §1.3 above.
  • asg_arc_del_short — for each node, drop arcs whose overlap length is below drop_ratio × max_arc_length. Default r ramps from 0.5 to 0.8 across rounds; final -F defaults to 0.8.
  • asg_cut_tip — uses a forward walk of at most max_ext hops from a candidate tip. If the walk doesn’t reach a mergeable interior, the whole sub-path is deleted.
  • asg_cut_internal and asg_cut_biloop — handle small internal stubs and adjacent two-edge loops that the standard tip routine misses.
  • asg_pop_bubble — Tarjan-style BFS from each branch point; if it converges to a single sink within max_dist, the non-best alternative is deleted.

Crucially, none of these consult haplotype information. Bubbles are popped to the single best path regardless of whether the alternative encodes a real heterozygous allele.

3. Hifiasm’s additions to the same skeleton

hifiasm keeps the miniasm primitives (sometimes byte-identically, the only differences are added timing/log statements) and surrounds them with three new categories of code:

3.1 Diploid-aware del-short variants

These consult per-read haplotype evidence rather than geometry alone:

  • asg_arc_del_short_diploid_unclean — the same drop-ratio loop as miniasm, but with one rescue: when the second-longest arc out of a node is about to be deleted, it is not deleted. This preserves the alternate haplotype branch out of every divergent node.
  • asg_arc_del_short_diploid_unclean_exact — stricter version.
  • asg_arc_del_short_diploid_by_length — length-conditioned.
  • asg_arc_del_short_diploid_by_length_trio — combines length conditioning with trio data.
  • asg_arc_del_short_diploid_by_exact — exact-match conditioned.
  • asg_arc_del_short_diploid_by_exact_trio — combines exact match with trio data.
  • asg_arc_del_short_diploi_by_suspect_edge (sic) — flags edges as “suspect” by phasing evidence.

3.3 Trio-aware bubble popping and trio-edge filtering

  • asg_pop_bubble_primary_trio. Distinct from miniasm’s asg_pop_bubble: the primary path is chosen by trio assignment, and the alternate is preserved into the a_ctg.gfa output rather than discarded.

  • drop_edges_by_trio. With trio data from parental short-read \(k\)-mers, each read gets a tag trio_flag[r] ∈ {paternal, maternal, ambiguous}. Any arc whose endpoints are tagged with opposite parents is deleted unconditionally:

    if (trio_flag[qn] != AMBIGU && trio_flag[tn] != AMBIGU)
    if (trio_flag[qn] != trio_flag[tn])
        arc.del = 1;

3.4 Wrapped around an iterative error-correction loop

asg_arc_del_trans is invoked at multiple call sites. Between calls, hifiasm re-runs base-level error correction and overlap recomputation. Miniasm does graph cleaning once over the read set it was given.

3.5 Multi-output graph

Miniasm emits one GFA. Hifiasm emits at minimum four: raw unitigs (r_utg.gfa), processed unitigs (p_utg.gfa), primary contigs (p_ctg.gfa), alternate contigs (a_ctg.gfa).

The operational summary: every step in the miniasm cascade has, in hifiasm, either a diploid- or trio-aware sibling, and the whole cascade is iterated inside an error-correction outer loop. The clean-and-collapse-everything heuristic of miniasm becomes a clean-while-preserving-heterozygosity discipline in hifiasm.

4. Replacing the heuristics with deep learning — what we tried

The motivating question: every routine in hifiasm’s diploid cascade encodes a hand-tuned decision boundary (drop ratios, fuzz tolerances, trio_flag inequalities, exact-match thresholds). Can a graph neural network learn those boundaries directly from the assembler’s own gold-standard output — well enough to take its place on real human data?

We attacked this in two phases. Phase 1 establishes a tensor-only pipeline on a haploid genome (yeast) where the gold standard is miniasm — a tractable warm-up because there is no heterozygosity to preserve, so any heuristic that collapses too aggressively gets caught immediately by N50. Phase 2 is the actual replacement target: hifiasm’s diploid r_utg graph on real HG002 HiFi reads.

4.1 Phase 1 — haploid groundwork (yeast / miniasm)

All yeast numbers below are evaluated against the SacCer3 reference with QUAST defaults. Reads: PBSIM3 (Ono et al. 2021) at --strategy wgs --method errhmm --depth 30 --length-mean 15000 --accuracy-mean 1.00 — a perfect-read stand-in for the output of a future learned base corrector. Overlapper: minimap2 -x ava-pb for the early GPU steps, hifiasm’s internal overlapper for the final.

Step 1 — boolean transitive reduction on GPU. Synthetic smoke test on a dense DAG (\(N = 1{,}000\), \(E = 25{,}570\), chains of 4–5 hops): torch.sparse_coo_tensor with \(M_r = M \setminus (M \cdot M)\), \(\cdot\) boolean. Result: bit-identical to networkx.transitive_reduction (both kept 2,325 arcs), GPU 0.149 s vs CPU 0.259 s. On real yeast at 21,591 raw overlaps the GPU was correct but over-aggressive — it dropped 21,102 arcs and left 489 vs miniasm’s 178, because it ran on the unfiltered hairball without the prior containment removal that miniasm does in its read-selection stages.

Lesson 1. A pure topological transitive reduction is correct in the Myers (2005) sense but useless on its own. The CPU cascade only works because every stage upstream of it (containment removal, two-pass read selection) has already done substantial filtering.

Step 2 — error-correction detour. Pure-tensor TR on uncorrected reads (PBSIM3 default ~99 %): the 1 % residual error rate baked massive error bubbles into the sparse matrix and the downstream layout produced a 117 MB hairball. The workaround was the --accuracy-mean 1.00 setting — a stand-in for the perfect-corrector limit. The real fix is a learned base corrector (HERRO-line architectures; see Stanojević et al. 2024 and DeepConsensus, Baid et al. 2023), separately developed.

Lesson 2. Graph cleaning and base correction are not separable. Hifiasm wraps cleaning inside an error-correction outer loop for a reason: if the reads contain even single-percent error, no amount of downstream topological surgery can recover.

Step 3 — min-plus tropical semiring. Boolean TR was deleting metrically-impossible “transitive” arcs through repeat regions. Replacing \((\times, +)\) with \((\min, +)\) as in §1.4 lifts the fuzz tolerance into the algebra. Implementation: torch.scatter_reduce_(..., reduce='amin') over the COO representation. Yeast N50: 97,312 → 100,730 bp.

Lesson 3. Myers’s fuzz parameter is not an implementation detail; it is the difference between Boolean reachability and metric reachability. The tropical semiring (Maclagan & Sturmfels 2015 textbook) is the right algebra for shortest-path-style operations on weighted graphs and lifts naturally onto GPU sparse kernels.

Step 4 — bidirected 2N × 2N tensors and the full miniasm cascade. A walk through the assembler internals showed that our single-strand sparse matrix was hallucinating transitive paths through what should have been reverse-complement nodes. Switching to the 2N representation with twin-arc enforcement (every \(u \to v\) deletion triggers the symmetric \(v^{c} \to u^{c}\) deletion) and adding the full miniasm-style cascade — fuzz TR + tip clip + bubble pop + iterated short-overlap removal — produced the second-to-last row below. Pushing it one step further (extra tip + bubble + short-overlap rounds beyond miniasm’s stop point) collapsed the assembly.

Configuration # contigs Total (bp) N50 (bp) Largest (bp)
miniasm (default 30× HiFi) 28 12,182,237 784,824 1,529,136
miniasm (100 % accurate reads) 19 11,925,444 784,847 1,529,136
GPU boolean TR (Step 1) 266 18,911,523 97,312 321,651
GPU + min-plus (Step 3) 121 10,447,510 100,730 225,564
GPU + bidirected (Step 4) 102 11,612,891 142,041 314,935
GPU + strict miniasm cascade 120 12,573,862 149,960 457,834
GPU + extra tip+bubble (over-pruned) 75 2,494,027 31,962 70,957
Dual-Headed GNN inference (Step 6 below) 28 12,180,284 784,824 1,529,136

Lesson 4. The DNA double helix is not optional. Any sparse-tensor representation of a string graph has to be 2N × 2N with twin-arc enforcement; otherwise the algebra hallucinates transitive paths through reverse-complement copies of the same read. Equivalently: \(M \in \{0,1\}^{2N \times 2N}\) must satisfy the involution \(M_{uv} = M_{v^c u^c}\) for every \(u, v\).

Lesson 5. Miniasm’s default parameters are a precariously-tuned local optimum. The “over-pruned” row shows that one extra tip-and-bubble round on a small repetitive genome collapses total length by 4.7× and N50 by 4.7×. The CPU cascade works because it happens to stop just before this catastrophe, not because the cascade itself is robust.

Step 5 — edge-only GNN. A 3-layer message-passing network classifying each arc:

\[ P_\text{edge}(u, v) = \sigma\!\bigl(\,\mathrm{MLP}\bigl([\mathbf{h}_u, \mathbf{h}_v, \mathbf{e}_{uv}]\bigr)\bigr). \tag{9}\]

Edge precision was excellent; assembly contiguity was destroyed. Reads with no surviving incident arcs stayed in the graph as zero-degree islands, and the layout engine terminated contigs at every dead-end.

Lesson 6. Edge classification alone is the wrong objective. The downstream layout depends on node survival; classifying edges without simultaneously deciding which reads to keep produces a graph that looks clean by per-arc precision but is structurally broken.

Step 6 — Dual-Headed GNN. Topology is jointly node-and-edge. Add a node head and mask edges by the boolean conjunction:

\[ \mathrm{Mask}(u, v) \;=\; \bigl[P_\text{node}(u) \ge \tau\bigr] \,\land\, \bigl[P_\text{node}(v) \ge \tau\bigr] \,\land\, \bigl[P_\text{edge}(u, v) \ge \tau\bigr]. \tag{10}\]

Trained on labels extracted from miniasm’s golden GFA (every arc on the unitig path → 1, else 0), this model reduces the 29,560-edge yeast graph to exactly 4,920 arcs at 100 % precision / 100 % recall — the bottom row above, N50 784,824 bp, bit-identical to miniasm.

Lesson 7. Topological masking — requiring both endpoints and the connecting edge to survive — is the minimal joint-classification scheme that respects the fact that the downstream consumer is a layout engine, not a per-edge auditor.

Phase 1 ends here. The yeast prototype proves the tensor pipeline can match a haploid heuristic exactly. Nothing in it tests the harder claim — that a GNN can replace diploid heuristics, where the right answer is not “collapse one branch” but “preserve heterozygosity while removing chimeric bridges.”

4.2 Phase 2 — diploid replacement (HG002 / hifiasm)

We switched ground truth from miniasm to hifiasm’s raw unitig graph (.r_utg.gfa) on real HG002 30× PacBio HiFi reads (Ashkenazi Jewish, NA24385) from the HPRC (Liao et al. 2023). Hifiasm’s internal cleanup produces ~321 k edges across ~250 disjoint connected components. Class imbalance: ~93 % “true biological” arcs / 7 % chimeric or repeat-collapse arcs.

Ground-truth labels. Each L-line endpoint is mapped against a diploid reference (HPRC assemblies; T2T-CHM13 v2.0 from Nurk et al. 2022 for the haploid fallback). An L-line is labeled 1 iff both endpoints share the same target sequence at mapq > 10, else 0. Two pipelines exist:

  • HG002 only — minimap2 (Li 2018) against haploid CHM13v2.0. Fast, heuristic, sufficient for initial training.
  • Tri-ancestry (HG005, HG00733) — wfmash (Guarracino et al. 2021) against the per-sample diploid <paternal + maternal> reference. Wavefront-algorithm (Marco-Sola et al. 2021) gives mathematically optimal base-level alignments, which we need for the cross-haplotype boundary calls.

The choice of wfmash for diploid labeling is deliberate: on a heterozygous SNP the two haplotypes are 99 %+ identical and minimap2’s heuristic chaining can resolve the mismap either way. WFA-based alignment is exact and the resulting labels are essentially noise-free — which we confirm empirically below by ablating label smoothing.

Architecture sweep. Three architectural waves on a single H200-class GPU. Class imbalance handled by BCEWithLogitsLoss with a dynamic pos_weight derived from each component’s edge-label distribution. Subgraph batching: connected-component batching rather than random NeighborLoader sampling, because edge dropping changes the graph Laplacian and breaks mathematical equivalence with full-batch training (verified to the 6th decimal place against a chr22 baseline). Hifiasm’s .r_utg is naturally fragmented into hundreds of disjoint components, so this is essentially Cluster-GCN (Chiang et al. 2019) with a pre-existing partition — no METIS step needed.

Wave Layers × hidden Conv Edge-feature handling Val edge acc
1 3 × 64–128 SAGEConv concat at output 87.81 %
2a 4 × 128 SAGEConv + dropout 0.2 concat at output 92.16 %
2b 4 × 128 GATv2Conv attention over edges unstable: peaks 95.5 %, collapses
3 5 × 512 GINEConv edge-conditioned messaging (§1.5) 97.08 %

The Wave 3 result is the empirical confirmation of the 1-WL expressivity bound: the network has to distinguish transitive arcs from tips from bubbles from cross-haplotype bridges in a single forward pass, and GINE’s edge-conditioned summation is the only one of the three that preserves enough overlap-geometry information through 5 layers of message passing to do that reliably.

Lesson 8. Architecture matters more than depth. SAGE 4×128 (Wave 2a) and GINE 5×512 (Wave 3) differ by 35× in parameter count, but the 5-point accuracy gap is dominated by how edge features enter message passing, not by capacity. Edge features concatenated at the readout only see the post-aggregation node embeddings; edge features injected into every layer of \(\psi\) propagate through the full receptive field.

Lesson 9. Attention is not free. GATv2 (Wave 2b) was theoretically the most expressive of the non-GINE options, but training was unstable in our setting — it spiked to 95.5 % mid-run and collapsed. We did not chase the instability; we suspect it is the standard “attention is sensitive to learning rate × batch composition” problem (Xiong et al. 2020) compounded by extreme class imbalance.

Ablation tournament. With the architecture fixed at GINE 5×512, we ran an independent ablation (each intervention solo) to identify what closes the residual ~3 % gap. Results:

Variant Val edge acc Δ
Baseline GINE (BCE) 95.50 % —
+ Focal Loss only (\(\gamma = 2\)) 97.44 % +1.94
+ Deep MLP (3-layer GINE-MLP) 97.08 % +1.58
+ GraphNorm only 96.52 % +1.02
+ Label smoothing only (\(\epsilon = 0.05\)) 93.43 % −2.07

Focal loss (Lin et al. 2017) reweights the BCE by \((1 - p_t)^\gamma\) where \(p_t\) is the model’s probability of the true class:

\[ \mathcal{L}_\text{focal}(p_t) \;=\; -\,(1 - p_t)^{\gamma} \,\log p_t,\qquad \gamma = 2. \tag{11}\]

When \(p_t\) is high (confident-and-correct) the prefactor \((1 - p_t)^\gamma \to 0\) and the loss zeros out; gradient capacity concentrates on the hard cases. For our 93 / 7 imbalance this dominates: a cumulative variant combining all four interventions scored 97.17 %, below focal alone — the regularisers actively interfere with one another once focal is in place. (Class-balanced loss, Cui et al. 2019, gives a similar reweighting via the effective-number-of-samples lens; we did not separately ablate it.)

Lesson 10. Focal loss is the right tool for the topological-imbalance problem and the gradient-imbalance problem at once. The 93 % of “easy” true biological edges would otherwise dominate the gradient and starve learning on the 7 % chimeric/bridge edges that are what we need the model for.

Lesson 11. Label smoothing is informative in reverse. A −2.07 pp regression from \(\epsilon = 0.05\) smoothing means the wfmash labels are clean enough that forcing the network to express uncertainty actively destroys it. This is unusual — most production ML labels benefit from ε ∈ [0.05, 0.1] — and it argues for keeping wfmash, not minimap2, as the source-of-truth aligner for diploid label generation.

Final locked architecture for the tri-ancestry scale-up:

\[ \text{GINEConv}\;(5\,\text{layers},\;512\,\text{hidden},\;\text{Focal}\,\gamma{=}2,\;\text{wfmash labels}) \;\rightarrow\; 97.44\%. \tag{12}\]

5. What we did not solve — open problems

The 97.44 % is not 100 %. A few things remain:

  1. Cross-haplotype bridges in alpha-satellite arrays. The hardest 2.5 % of edges concentrate in centromeric repeats where the paternal and maternal arrays differ structurally rather than by SNPs. The current GINE features (overlap length, identity, match-length) do not include any signal about sequence-level repeat structure; the model has to infer it from topology alone. A natural extension is to add per-edge \(k\)-mer-frequency features (cf. Vrček et al. 2022’s compositional features) or to integrate a small frozen DNA-LM embedding of the overlap sequence as an edge feature.
  2. No trio data. Hifiasm’s _trio family of routines uses parental short-read \(k\)-mers to tag each read as paternal/maternal/ambiguous. We are training without trio for ancestry generality — HG002 has trio data, HG005 and HG00733 do not — but this means we are implicitly asking the network to infer haplotype assignment from graph topology alone. Adding a phasing pre-step (e.g. inheriting trio_flag labels from hifiasm’s own assignment, where available) would likely close another 1–2 percentage points.
  3. Layout, not just pruning. This work replaces hifiasm’s graph-cleaning cascade. The downstream unitig layout (the path-walking that emits p_ctg and a_ctg) is still C, and remains the next obvious GPU target.
  4. The base-corrector is upstream and absent. Everything above is conditioned on either real HiFi reads (which hifiasm has already error-corrected internally) or PBSIM3 perfect simulation. A standalone learned corrector — a permutation-invariant transformer over the read pileup, trained with wfmash CIGAR ground truth (HERRO-line; Stanojević et al. 2024; DeepConsensus, Baid et al. 2023) — is the upstream piece that closes the full pipeline.

6. Full ledger

# Step Outcome Lesson
1 Boolean GPU TR (Phase 1) ✓ bit-identical on synthetic; ❌ over-pruned on real No containment prefilter; topology ≠ cleaning
2 Assembly from uncorrected reads ❌ 117 MB hairball Cleaning and base correction inseparable
3 Simulated 100 % accurate reads ✓ proves topology math Stand-in for a learned base corrector
4 Min-plus tropical semiring ✓ +3 kb N50 \((\times,+)\) ignores fuzz; tropical algebra lifts it
5 Bidirected 2N × 2N ✓ +42 kb N50 Single-strand tensors hallucinate transitive paths through RC nodes
6 Full miniasm cascade in PyTorch ✓ matches miniasm’s N50 Faithful port of fuzz TR + tip + bubble + iterated short-ovlp
7 Push cleanup further (tip+bubble extras) ❌ 11.6 MB → 2.5 MB Miniasm defaults are a precariously-tuned local optimum
8 Edge-only GNN ❌ fragmentation Edge precision is the wrong objective when downstream is layout
9 Dual-Headed GNN (node + edge head) ✓ 100 % parity on yeast Topological mask requires both endpoints
10 Scale to HG002 / hifiasm r_utg (Phase 2) ✓ ~250 components, 93/7 imbalance Diploid replacement is the real test
11 Arch sweep SAGE → GATv2 → GINE ✓ 87.81 → 92.16 → 97.08 % Edge-conditioned messaging dominates concat-after
12 GATv2 attention ❌ unstable Attention sensitivity × extreme imbalance
13 ERA: focal loss ✓ 95.50 → 97.44 % Class-imbalance fix; dominates other interventions
14 ERA: deep MLP, GraphNorm ✓ smaller wins; clash with focal WL-expressivity gain saturates fast
15 ERA: label smoothing (\(\epsilon = 0.05\)) ❌ −2.07 % wfmash labels too clean for smoothing — empirical proof of label quality
16 NeighborLoader random subgraph batching ❌ rejected Edge dropping changes graph Laplacian
17 Connected-component batching ✓ bit-identical to full-batch Hifiasm .r_utg is naturally fragmented (Cluster-GCN with free partition)

7. References

OLC and string graph foundations

  1. Myers, E.W. (2005) The fragment assembly string graph. Bioinformatics 21 (Suppl. 2): ii79–ii85. — Transitive-reduction theorem used in §1.3.
  2. Li, H. (2016) Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences. Bioinformatics 32 (14): 2103–2110. — Reference for §2.
  3. Cheng, H., Concepcion, G.T., Feng, X., Zhang, H., Li, H. (2021) Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nature Methods 18: 170–175. — Reference for §3.
  4. Chin, C.-S. et al. (2016) Phased diploid genome assembly with single-molecule real-time sequencing (FALCON-Unzip). Nature Methods 13: 1050–1054. — Earliest diploid-aware long-read assembler; design ancestor of hifiasm’s p_ctg/a_ctg split.
  5. Koren, S., Walenz, B.P., Berlin, K., Miller, J.R., Bergman, N.H., Phillippy, A.M. (2017) Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Research 27: 722–736. — Long-read assembler with explicit correction-then-assembly stages.
  6. Koren, S. et al. (2018) De novo assembly of haplotype-resolved genomes with trio binning. Nature Biotechnology 36: 1174–1182. — Source of the parental-\(k\)-mer trio approach hifiasm’s _trio routines implement.
  7. Rautiainen, M. et al. (2023) Telomere-to-telomere assembly of diploid chromosomes with Verkko. Nature Biotechnology 41: 1474–1482. — Current state-of-the-art for diploid assembly; produced the HPRC diploid references used here as ground truth.

Alignment

  1. Li, H. (2018) Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34 (18): 3094–3100. — Heuristic aligner; used for the all-vs-all ava-pb overlapper in §4.1 and HG002 haploid labelling in §4.2.
  2. Marco-Sola, S., Moure, J.C., Moreto, M., Espinosa, A. (2021) Fast gap-affine pairwise alignment using the wavefront algorithm. Bioinformatics 37 (4): 456–463. — WFA, the optimal-alignment kernel used by wfmash.
  3. Guarracino, A., Heumos, S., Nahnsen, S., Prins, P., Garrison, E. (2021) wfmash: whole-chromosome pairwise alignment using the hierarchical wavefront algorithm. Bioinformatics. — Diploid-reference labeller in §4.2.
  4. Marçais, G. et al. (2018) MUMmer4: A fast and versatile genome alignment system. PLoS Comput Biol 14: e1005944. — Canonical genome-alignment benchmark.

Reference data

  1. Liao, W.-W. et al. (2023) A draft human pangenome reference. Nature 617: 312–324. — HPRC source for the HG002 / HG005 / HG00733 reads and diploid references.
  2. Nurk, S. et al. (2022) The complete sequence of a human genome. Science 376: 44–53. — T2T-CHM13 v2.0, the haploid reference used for HG002 labelling.

Graph neural networks

  1. Gilmer, J., Schoenholz, S.S., Riley, P.F., Vinyals, O., Dahl, G.E. (2017) Neural message passing for quantum chemistry. ICML. — The MPNN framework that subsumes all the architectures in §4.2.
  2. Kipf, T.N., Welling, M. (2017) Semi-supervised classification with graph convolutional networks. ICLR. — GCN, baseline architecture not directly tried but referenced in §1.5.
  3. Hamilton, W., Ying, R., Leskovec, J. (2017) Inductive representation learning on large graphs. NeurIPS. — GraphSAGE, our Wave-1 baseline.
  4. Veličković, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., Bengio, Y. (2018) Graph attention networks. ICLR. — Original GAT.
  5. Brody, S., Alon, U., Yahav, E. (2022) How attentive are graph attention networks? ICLR. — GATv2, our Wave 2b attempt.
  6. Xu, K., Hu, W., Leskovec, J., Jegelka, S. (2019) How powerful are graph neural networks? ICLR. — Defines GIN; proves the 1-WL bound used in §1.5.
  7. Hu, W. et al. (2020) Strategies for pre-training graph neural networks. ICLR. — Introduces GINEConv, the edge-conditioned variant in the final architecture.
  8. Corso, G., Cavalleri, L., Beaini, D., Liò, P., Veličković, P. (2020) Principal Neighbourhood Aggregation for graph nets (PNA). NeurIPS. — Strictly more expressive aggregation family than GIN; not tried in this work.
  9. Dwivedi, V.P., Bresson, X. (2021) A generalization of transformer networks to graphs. AAAI workshop. — Graph transformer; another architecture not tried.
  10. Cai, T. et al. (2021) GraphNorm: a principled approach to accelerating graph neural network training. ICML. — The normalisation ablated in §4.2.

Graph batching

  1. Chiang, W.-L., Liu, X., Si, S., Li, Y., Bengio, S., Hsieh, C.-J. (2019) Cluster-GCN: an efficient algorithm for training deep and large graph convolutional networks. KDD. — The method our connected-component batching is a free-partition special case of.
  2. Zeng, H., Zhou, H., Srivastava, A., Kannan, R., Prasanna, V. (2020) GraphSAINT: graph sampling based inductive learning method. ICLR. — Alternative sampling-based batching; rejected here for the Laplacian-distortion reason in §4.2.

Loss and regularisation

  1. Lin, T.-Y., Goyal, P., Girshick, R., He, K., Dollár, P. (2017) Focal loss for dense object detection. ICCV. — The 95.50 → 97.44 % win in §4.2.
  2. Cui, Y., Jia, M., Lin, T.-Y., Song, Y., Belongie, S. (2019) Class-balanced loss based on effective number of samples. CVPR. — Alternative imbalance treatment; not ablated.
  3. Müller, R., Kornblith, S., Hinton, G. (2019) When does label smoothing help? NeurIPS. — Theoretical companion to our −2.07 % label-smoothing result.
  4. Xiong, R. et al. (2020) On layer normalization in the transformer architecture. ICML. — Background for the GATv2 instability hypothesis in §4.2.

Learning-based assembly

  1. Vrček, L., Burdyshaw, B., Mihaljević, M., Bresson, X., Stošić, B.D., Šikić, M. (2022) Learning to untangle genome assembly with graph convolutional networks (the GNNome system). — First clean formulation of edge classification on overlap graphs; supervised-classification framing we extend with the dual-headed mask.

Long-read error correction (upstream)

  1. Stanojević, D. et al. (2024) HERRO: end-to-end accurate read correction with deep learning. bioRxiv preprint. — Reference base-correction architecture.
  2. Baid, G. et al. (2023) DeepConsensus improves the accuracy of sequences with a gap-aware sequence transformer (DeepConsensus). Nature Biotechnology 41: 232–238. — Production DL HiFi consensus, a parallel design to the HERRO line.

Simulation

  1. Ono, Y., Asai, K., Hamada, M. (2021) PBSIM3: a simulator for all types of PacBio and ONT long reads. NAR Genomics and Bioinformatics 4: lqac092. — Generator of the simulated reads used throughout §4.1.

Tropical semiring

  1. Maclagan, D., Sturmfels, B. (2015) Introduction to Tropical Geometry. Graduate Studies in Mathematics, AMS. — Reference for the \((\min, +)\) algebra used in §1.4.

Data note. HG002 (NA24385, Ashkenazi), HG005 (Han Chinese), HG00733 (Puerto Rican) 30× PacBio HiFi reads from the Human Pangenome Reference Consortium (HPRC). T2T-CHM13 v2.0 (Nurk et al. 2022) used as the haploid reference; per-sample diploid references from HPRC (Liao et al. 2023) used for the diploid labelling pipeline.

From Assembly Heuristics to Blazing-Fast Matrix Operations: Achieving 100% Genome Assembly Parity

2026-05-25

For decades, the standard paradigm of de novo genome assembly has relied on highly complex, rigid heuristics. Gold-standard tools like miniasm and hifiasm are genuine engineering marvels, but beneath the hood, their topological graph-cleaning phases are dominated by massive, serialised for-loops. These heuristics scan raw overlaps to filter out contained reads, chimeric sequences, and spurious alignments, one edge at a time.

As we enter the era of pangenomics and massive datasets, these serialised CPU operations create severe scaling bottlenecks. A single human genome produces billions of pairwise overlaps. Walking that graph sequentially on a CPU takes days, or at best about a day with heavy multithreading.

So we asked: what if we stripped away the heuristics entirely? What if we reframed topological graph cleaning not as a series of imperative rules, but as a pure mathematical operation, namely a massive, parallelised matrix classification problem that runs natively on a GPU?

Here, I document how we got there. It is also the story of every wall we hit, every hypothesis that failed, and every pivot we made along the way until we achieved 100% bitwise structural parity with the gold-standard CPU assembly on a yeast dataset using a dual-headed graph neural network.

The Hypothesis: Assembly as Matrix Operations

The core hypothesis was deceptively simple: structural heuristics are just unlearned classification thresholds.

In a traditional assembly graph, reads are nodes (\(V\)) and pairwise overlaps are edges (\(E\)). Tools like miniasm define rigid geometric rules, including span thresholds, overhang ratios, and containment checks, to decide which edges to delete. We hypothesised that if we converted the raw overlap geometry (read lengths, overlap spans, mapping qualities, and sequence identities) into a feature tensor \(X\), a neural network could learn the exact multidimensional decision boundary required to prune the graph optimally.

Because GPUs excel at large-scale matrix multiplications, processing a 50-million-edge graph becomes a single, highly parallelised millisecond-scale tensor operation rather than a multi-day CPU loop.

The Graph Theory Underpinnings: Abstracting OLC Topology

To appreciate why this matters, you need to understand what is happening inside the Overlap-Layout-Consensus (OLC) pipeline. After a sequence aligner like minimap2 dumps all pairwise overlaps into a PAF file, the resulting string graph is incredibly noisy, even after base correction, with thousands of redundant paths, spurious edges from repeat sequences, and chimeric artefacts from sequencing errors.

Historically, graph-theory algorithms are deployed in sequential cascades to clean this mess:

  1. Transitive Reduction (Myers, 2005): If directed edges \(A \rightarrow B\), \(B \rightarrow C\), and \(A \rightarrow C\) all exist, the edge \(A \rightarrow C\) is mathematically redundant (transitive) and must be pruned to reveal the true linear unitig path.
  2. Tip Clipping: Dead-end paths shorter than a read length are assumed to be sequencing errors and are recursively sheared off.
  3. Bubble Popping: Bifurcating paths that immediately reconverge represent heterozygosity or local errors and are collapsed into a single consensus path.

In traditional tools, each of these is solved by iterative Breadth-First Searches (BFS) and sub-graph traversals. The CPU must “walk” the graph, inspecting neighbourhoods sequentially. Our bet was that a graph neural network, a multi-layer message-passing network, could aggregate the structural neighbourhood of each edge into a latent vector \(\mathbf{h}_v\). Instead of hand-coding traversal rules, the GNN would learn to recognise the topological signatures of transitive edges, tips, and bubbles implicitly, compressing all those graph-walking algorithms into a single feed-forward classification pass.

But before we got to the GNN, we went on quite a detour.

Step 1: GPU Sparse Boolean Transitive Reduction (The First “Aha”)

Our first move was to attack transitive reduction directly with pure matrix algebra. The math is elegant:

Let \(M\) be the adjacency matrix of a Directed Acyclic Graph. The transitive closure \(M^*\) can be iteratively computed by \(M_k = M_{k-1} \times M\), where \(M_k\) contains all paths of exactly length \(k\). The fully reduced graph is:

\[ M_{reduced} = M \setminus \left( M \cap \bigcup_{k=2}^{\infty} M_k \right) \tag{13}\]

We implemented this entirely in torch.sparse_coo_tensor, with no Python loops and no CPU dispatch overhead, just raw matrix multiplications on the GPU.

The Smoke Test: Bit-Identical to NetworkX

We generated a dense synthetic DAG (\(N=1{,}000\) nodes, \(E=25{,}570\) edges) with deep transitive structures intentionally designed to break naive pathing (chains of 4-5 hops with all possible shortcuts). Then we ran our GPU kernel against the CPU ground truth (networkx.transitive_reduction):

Metric CPU (networkx) GPU (PyTorch Sparse)
Edges kept 2,325 2,325
Time 0.259s 0.149s

The result was bit-identical. The math was proven correct on synthetic data.

Yeast Benchmark vs Miniasm: Reality Check

Then we ran it on real data, the yeast 30× overlaps generated by minimap2, and things got interesting.

Metric CPU (miniasm) GPU (PyTorch Sparse)
Raw Input Overlaps 21,591 21,591
Pre-Reduction Edges 1,021 21,591
Edges Dropped 843 21,102
Final Spanning Edges 178 489

The GPU was far more aggressive because it operated on the raw adjacency matrix without any prior containment filtering. miniasm first removes contained reads (coordinate-based matching on the CPU), which strips the graph from 21,591 to 1,021 edges before it even begins transitive reduction. Our GPU algorithm just saw the full dense hairball and mathematically pruned everything it could find.

Was it correct? Mathematically, yes. Was it useful? Not quite yet.

Step 2: The Error-Correction Problem (A Big Detour)

When we tried to actually assemble from the GPU-reduced edges, things went sideways.

Even at 99% read accuracy, the remaining 1% of errors bake massive error bubbles and spurious tips directly into the sparse adjacency matrix. The downstream layout engine either crashed or produced a 117 MB hairball assembly, orders of magnitude worse than what miniasm produces.

So, we needed clean reads, which will be the subject of the next post. Traditional assemblers like hifiasm spend enormous CPU time on error correction before they even touch the graph.

Since our deep learning error-correction pipeline was not ready yet, we sidestepped the problem by simulating 100% accurate reads using PBSIM3 with --accuracy-mean 1.00. This represents the exact output of a hypothetical perfect error corrector: zero sequencing errors and pure structural signal.

With these clean reads:

pbsim --strategy wgs --method errhmm --depth 30.0 --length-mean 15000 \
      --accuracy-mean 1.00 --genome yeast_clean.fa --prefix yeast_100pct
minimap2 -x ava-pb -t 20 yeast_100pct_full.fq.gz yeast_100pct_full.fq.gz > yeast_100pct_ava.paf
python3 scripts/filter_containments.py yeast_100pct_ava.paf yeast_100pct_filtered.paf
python3 scripts/gpu_transitive_reduction.py --paf yeast_100pct_filtered.paf --out yeast_100pct_reduced.paf

We were skipping that step entirely and expecting the graph math to somehow compensate. It could not.

The GPU transitive reduction achieved:

  • N50: 97,312 bp (2x better than the uncorrected run)
  • Total assembled length: 18.9 MB
  • Transitive edges pruned: 13,633 redundant paths cleanly removed

This confirmed the structural math was sound. But the N50 was still far short of the CPU engine’s 812 kb. Something else was wrong.

Step 3: The Min-Plus Tropical Semiring (Fixing the Algebra)

We traced the N50 gap to a fundamental problem with boolean algebra in repetitive genomes.

Standard boolean matrix multiplication treats all edges as equal. A 500 bp overlap and a 15,000 bp overlap are both just 1 in the adjacency matrix. But in transitive reduction, path length matters. Miniasm’s CPU algorithm uses a “fuzz” parameter: an edge \(A \rightarrow C\) is only transitive if the sum of the intermediate edges (\(A \rightarrow B + B \rightarrow C\)) is within a distance threshold.

Our boolean GPU kernel was ignoring distance entirely, which caused it to shatter repetitive genomes by pruning edges that were topologically but not metrically transitive. In other words, it saw a mathematical shortcut through a repeat region and deleted the longer direct edge, even though that “shortcut” was biologically meaningless because it jumped through a completely different genomic region.

We upgraded to a Min-Plus Tropical Semiring, where matrix “multiplication” computes \(\min\) instead of \(\sum\) and “addition” is \(+\) instead of \(\times\). In PyTorch, we implemented this using scatter_reduce_ with amin.

The Min-Plus engine correctly preserved hundreds of essential repetitive structural paths that the boolean logic had wrongfully deleted. N50 improved from 97 kb to 114 kb. But 114 kb was still far from 812 kb.

Step 4: The Bidirected Tensor Discovery (The Root Cause)

This was the breakthrough in understanding, and also the most frustrating one, because it revealed a fundamental architectural limitation we could not easily fix.

A biological overlap graph is inherently bidirected: every read has a forward (+) and reverse complement (−) node. In the 2N × 2N adjacency matrix that miniasm uses internally, each overlap creates two directed arcs:

U = q_id * 2 + u_dir    →    V = t_id * 2 + v_dir
V' = V ^ 1              →    U' = U ^ 1

Our GPU tensor had been modelling the graph using simple sequence IDs, effectively treating it as undirected. In an undirected graph, if Read A overlaps Read B and Read C overlaps Read B, the tensor detects a mathematical transitive path \(A \to B \to C\).

But biologically, A might connect to the 3’ end of \(B^+\) while C connects to the 3’ end of \(B^-\). In a strictly bidirected tensor, these edges live on different nodes (\(B^+\) vs \(B^-\)), so no transitive path exists between A and C.

Because our GPU tensor was strand-agnostic, it was hallucinating millions of biologically impossible transitive paths through nodes, causing the Min-Plus algorithm to erroneously prune true structural edges.

So, in the next step, we used the full bidirected 2N × 2N representation with explicit strand-aware arc computation, symmetric twin enforcement (V' → U' whenever U → V is deleted), and all four of miniasm’s graph-cleaning phases:

  1. Strict Transitive Reduction (fuzz-aware, distance-weighted)
  2. Tip Clipping (dead-end path removal up to max_ext hops)
  3. Bubble Popping (bifurcation/reconvergence detection via BFS)
  4. Short Overlap Removal (span-ratio based filtering)

This was no longer a cute matrix multiply. It was a faithful port of miniasm’s entire state machine into PyTorch tensors.

Step 5: The Edge-Only GNN (High Precision, Shattered Assembly)

With the transitive reduction journey behind us, we pivoted to a learning-based approach. Instead of trying to replicate miniasm’s exact heuristic cascade, why not let a neural network learn the whole thing end-to-end?

We trained a standard GNN to predict the survival probability of each overlap:

\[ P_{edge}(u, v) = \sigma\left( \text{MLP}_{edge}( \mathbf{h}_u, \mathbf{h}_v, \mathbf{e}_{u,v} ) \right) \tag{14}\]

where \(\mathbf{h}_u\) and \(\mathbf{h}_v\) are learned node embeddings propagated through 3 message-passing layers, and \(\mathbf{e}_{u,v}\) are the 10-dimensional geometric features of the overlap (distance, span, matching bases, block length, mapping quality, identity, and normalised coordinate ratios).

The model achieved extremely high precision, but the assembly fragmented. The N50 plummeted.

The Mathematical Flaw

By only classifying edges, the network successfully deleted spurious overlaps, but it ignored the biological reality of the reads themselves. Reads that were deeply contained or chimeric lost all their valid edges, leaving them as “dead-end” floating nodes in the graph. When the layout algorithm encountered these dead-ends, it assumed the contig had terminated, artificially fracturing the assembly path.

This was a classic case of optimising the wrong objective. The edge-level precision was great, but the graph-level contiguity was destroyed.

Step 6: The Dual-Headed GNN (The One That Worked)

The fix was conceptually simple but architecturally significant: structural topology is fundamentally joint. A valid edge cannot exist between invalid nodes.

We architected a Dual-Headed GNN (DualHeadedGNN in train_gnn.py) with 3 message-passing layers (each with its own msg_mlp and update_mlp) and two classification heads that simultaneously predict:

1. Node Head: whether a read should survive (not contained, not chimeric):

\[ P_{node}(u) = \sigma\left( \text{MLP}_{node}( \mathbf{h}_u ) \right) \tag{15}\]

2. Edge Head: whether an overlap is structurally valid:

\[ P_{edge}(u, v) = \sigma\left( \text{MLP}_{edge}( \mathbf{h}_u, \mathbf{h}_v, \mathbf{e}_{u,v} ) \right) \tag{16}\]

Ground Truth Generation

To train this network, we needed labels. We built the training set, which takes two inputs: 1. The raw PAF overlaps from minimap2 (noisy, 29,560 edges). 2. The golden GFA unitig file from miniasm (the CPU gold standard).

The script parses the GFA’s a-records (which list the reads in each unitig and their strand orientations), traces the consecutive read pairs, maps them back to bidirected arc coordinates (\(U = \text{read\_id} \times 2 + \text{dir}\)), and creates binary labels for both nodes and edges. Every edge in the raw graph gets a 1 if it appears in the golden assembly, 0 otherwise. Same for nodes.

The Topological Masking Function

During inference, instead of naively thresholding edges, we enforce a strict topological mask. An edge is only retained if both of its connecting nodes survive and the edge itself is classified as valid:

\[ \text{Mask}_{final} = \big[ P_{node}(u) \ge 0.5 \big] \land \big[ P_{node}(v) \ge 0.5 \big] \land \big[ P_{edge}(u, v) \ge 0.5 \big] \tag{17}\]

This boolean intersection is natively vectorised in a PyTorch implementation, executing across millions of edges in milliseconds. And because node validity is checked first, we never get the “dead-end contig” problem that killed the edge-only GNN.

After the GNN mask is applied, we run the surviving edges through the bidirected topological cleanup cascade (TR → tips → bubbles → short overlaps) for mathematical completeness.

Loss Formulation

Because false negatives (deleting a true edge) irrevocably fracture an assembly, and you cannot stitch a broken contig back together, we heavily weighted the positive class. We used joint Weighted Binary Cross Entropy with pos_weight set to \(2\times\) the class imbalance ratio:

\[ \mathcal{L}_{total} = \mathcal{L}_{BCE\_node}(\mathbf{w}_{node}) + \mathcal{L}_{BCE\_edge}(\mathbf{w}_{edge}) \tag{18}\]

100% Parity: The Results

We trained this Dual-Headed GNN on ground-truth decisions extracted from miniasm across a highly complex yeast dataset.

The raw input generated by minimap2 contained a massive 29,560-edge noisy graph. During inference, the GNN executed the forward pass and applied the topological masking function, perfectly stripping the noise.

The GPU outputted exactly 4,920 edges, a 100% precision and 100% recall match with the gold-standard CPU heuristic filtering.

The Future: Scaling to the Mammalian Era

The yeast benchmark proves the concept: rigid, serial C-heuristics can be successfully replaced by vectorised deep learning matrices without sacrificing a single base pair of contiguity. But yeast is 12 Mb. The human genome is 3 Gb, 250× larger.

The gap between our prototype and hifiasm is enormous, but it is entirely algorithmic. The bottleneck is graph cleaning, which is exactly what the GNN replaces. And unlike hifiasm’s CPU heuristics, our tensor operations scale linearly across thousands of CUDA cores.

Two architectural requirements must be met for a production-tier GPU assembler: 1. Bidirected Sparse Tensors: The adjacency matrix must be 2N × 2N to strictly separate forward and reverse sequence strands. 2. Min-Plus Tropical Semirings: Transitive reduction must use distance-weighted algebra, not unweighted boolean multiplication.

The era of deep-learning-driven assembly is here.

References & Literature Context

Our Dual-Headed architecture builds upon an emerging body of research that reframes genome assembly graphs as deep learning problems:

  1. GNNome (2021): One of the pioneering frameworks that formulated haploid genome assembly as a supervised binary edge classification task, proving that Graph Convolutional Networks (GCNs) could learn true adjacencies from simulated overlap graphs.
  2. COGRAM: Demonstrated the feasibility of using lightweight GCNs to score and classify edges on k-mer overlap graphs.
  3. DipGNNome (2023): Extended the GNNome architecture to diploid genome assembly, highlighting the power of GNNs in haplotype-aware graph processing.

Our Dual-Headed GNN advances this literature by explicitly solving the “dead-end contig” fragmentation problem through joint node-and-edge classification, allowing the pipeline to execute topological masking that perfectly mirrors the structural integrity of C-based heuristics at a fraction of the serialised compute time.