Co-lexicographic orders of pangenome automata, in Zig. For now it answers one question, which is how wide the co-lex order of a real variation graph is.
An FM-index works because the prefixes of a text can be sorted. Cotumaccio and Prezza extended the idea to any finite automaton: sort its states by the strings that reach them, in co-lexicographic order (compared from the last character). However, most automata admit only a partial order, and the index built on it costs time and space that grow with the order's width p, the largest set of states that cannot be put in sequence. A Wheeler graph is the case p = 1. The coarsest forward-stable (CFS) order of Becker et al. can be computed in polynomial time, and its width is never larger than the co-lex width of the automaton.
A pangenome built from a reference and a VCF is such an automaton. If its width were small, one index would cover every haplotype combination with no pruning, which is the gap the plan for this library started from.
colex builds the variation graph of a reference and a list of variants (nfa.fromVariants, with VCF alleles trimmed of their anchor base) and, for an acyclic graph, the least and greatest string reaching each state (order.intervals). The infima form a forest, since each state takes its best predecessor as parent, and so do the suprema. Parents are chosen in topological order by comparing two candidate chains with polynomial hashing and Myers' jump pointers, and all 2n chain strings are then ranked together by prefix doubling with counting sorts.
The intervals give two numbers.
order.widthis the width of the interval order (u before v when u's greatest string is at most v's least one). It is a co-lex relation, so this is an upper bound on the co-lex width and on the width of the CFS order.order.certifiedWidthis a lower bound for every co-lex order and for the CFS order. It keeps, among the states at the widest point, those whose string-length ranges are pairwise disjoint. Their string sets are then disjoint, and overlapping intervals make them incomparable in any co-lex order. The proof is short and sits in the source.
The tests compare the ranks with explicit string sets on small graphs, and both bounds with the exact maximum co-lex relation, computed by brute force.
On GRCh38 chr22 with the 1000 Genomes high-coverage variant sites, the width is in the thousands per megabase (median 8 180 over 1 Mb windows). With SNVs alone, the lower and upper bounds coincide, so these are exact widths: median 7 869 per megabase and 805 per 100 kb, about one state in the antichain for every three SNVs. The cause is local. A state right after an A/T SNP is reached by strings ending in both A and T, so its interval straddles every state whose strings end in C or G before the same base, and states after different SNPs straddle each other.
A co-lex index on the variation graph itself therefore pays for a width in the thousands. Indexing these graphs through co-lex orders needs another automaton for the same language (a Wheeler DFA, for instance), not a better order of this one. Full tables are in bench/RESULTS.md.
- Acyclic graphs only. Cyclic pangenome graphs (from graph builders rather than a VCF) need the infima of infinite strings.
- Alternate alleles join the reference on both sides, so combinations of overlapping variants are not spelled.
- It does not compute the CFS order itself, only bounds on its width.
Requires Zig 0.16.0.
const colex = @import("colex");
var n = try colex.nfa.fromVariants(gpa, reference, variants);
defer n.deinit(gpa);
var iv = try colex.order.intervals(gpa, n);
defer iv.deinit(gpa);
const w = try colex.order.width(gpa, iv);
const low = try colex.order.certifiedWidth(gpa, n, iv, w.at);zig build test
zig build census -Doptimize=ReleaseFast
./zig-out/bin/colex-census chr22.fa sites.vcf 1000000 # one line per 1 Mb windowSee REFERENCES.md.
GNU AGPL-3.0 (see LICENSE), copyright endlessconflict. A commercial license without the AGPL's network-copyleft terms is available on request; open an issue to ask.