Data Structures for Biology

Woolf Software
bioinformatics data structures pangenome single-cell data engineering
A lab instrument holding a glowing branching strand above a tray of thousands of individually lit tiny wells.

This essay was prompted by an episode of Boundary’s “AI That Works” series, published as the video “The Data Structure Problem AI Coding Agents Can’t See”. The first half is a talk by Avery, an engineer at Boundary. He opens with the old Fred Brooks and Rob Pike position that the tables explain the flowcharts, compressed into three words: “representation shapes implementation”. The rest of the talk tests whether that still holds when the programmer is a model. We think it does, and we think it matters more in biology than in ticketing, because in bioinformatics the data structure is usually the file format.

What the video argues

The example is a ticket seller. A seat can be open, sold, or held while a customer is in checkout. The naive representation is two booleans, held and sold, which permits a seat that is both held and sold at the same time and forces every transition to check for that. The alternative is an enumeration with three variants (a tagged union), which makes the impossible state unrepresentable and the transition code short. Adding a fourth state to the booleans means touching several transitions by hand, while adding it to the enumeration lets the type checker find every place that needs to change.

Avery had a frontier coding agent build the ticket system one feature at a time, with later requirements hidden, and ran each step 200 times. On the first feature nearly every run chose the obvious map from seat label to state. The second feature asked for a stable ordering of seats for a printable chart. About 40 percent of runs added a separate vector of seat labels beside the map, a second source of truth that can drift from the first. About 35 percent switched to an ordered map (a B-tree map in Rust), which he considers right because the ordering is then a property of the structure. A few percent added a constant-time counter of open seats that nobody asked for. On the final feature, timed holds, 87.5 percent of runs arrived at the ideal representation and 12.5 percent ended with three sources of truth and permitted invalid states. Every run had written tests, and the tests were green the whole way through.

Two conclusions follow. When all requirements are given up front the models choose well, so it is incremental delivery that produces bad structures. Since a rewrite per feature is unaffordable, the practical advice is to review the data structures rather than every downstream line, because the structures are what the downstream code inherits. Dex, one of the hosts, adds the arithmetic: a mistake that compounds at 12 percent per feature is not a small rate once a hundred features have shipped. His summary line is “all slop starts as bad data structures”. The second half of the episode covers HumanLayer’s sync engine, in which a client subscribes to a “shape” (a thin slice of a large dataset) and the architecture exists to keep those slices current.

The same argument in bioinformatics

In bioinformatics the representation was chosen early and everyone inherits it. A genome is a string and a variant set is sorted rows. An expression experiment is a matrix and a population is a tree. The sections below take these in roughly the order they appeared and ask what each makes cheap.

Strings and their indexes

The reference genome is a string of about three billion characters, and the central operation of short-read sequencing is to find where a short read of 100 to 150 characters occurs in it, allowing a few mismatches. Scanning the string per read is hopeless at hundreds of millions of reads, and the structure that made this tractable is the Burrows-Wheeler transform (BWT), a reversible permutation of the text that groups characters by the context that follows them. The FM-index of Ferragina and Manzini adds small auxiliary tables to the BWT so that a pattern can be located by a procedure called backward search 1.

BWA applied this to short reads. Li and Durbin describe backward search on the BWT as a way to mimic a top-down walk of the prefix trie of the genome without ever building the trie. Exact matches are counted in time proportional to the read length and independent of the genome size. The whole index for a 3 gigabase genome takes about 2.3 GB of memory, and the program ran 10 to 20 times faster than MAQ at similar accuracy 2. An FM-index answers “where does this substring occur” for any substring, because exact repeats collapse onto a single path of the trie. It does not naturally answer “what lies between positions 10,000 and 20,000”, which is a different question and a different structure.

Intervals

Once reads are aligned, nearly everything else in genomics is a set of intervals on a coordinate system. GTF genes, BED regions, and VCF variants are each a row with a chromosome and a start and usually an end. VCF was designed for the 1000 Genomes Project to hold polymorphisms and genotypes across thousands of samples at millions of sites, and the paper already treats bgzip compression and tabix indexing as part of the format 3. The question these files exist to answer is “which records overlap this region”.

Tabix answers it for any position-sorted, tab-delimited file. The index has two parts. A binning index assigns records to a hierarchy of bins from 512 Mb down to 16 kb and stores the file offsets of the records in each bin. A linear index records for every 16 kb window the offset of the leftmost record that could overlap it, so a query does not have to scan earlier bins. Li reports about one seek per query with the linear index against 6 to 13 without it. Since bgzip produces gzip-compatible blocks, the same index works over HTTP against a remote file 4. A coordinate-sorted file is already a data structure, and an unsorted VCF is the boolean-flag version of one: every consumer must either sort it or scan it.

Graphs

A linear reference represents one version of each locus. Reads from a person whose sequence differs from that version align worse or not at all. This is reference bias, and it is worst at structural variants. Garrison and colleagues describe variation graphs as bidirected sequence graphs in which paths correspond to haplotypes, so that population variation including inversions and duplications lives in the reference itself. Their vg toolkit maps reads to such graphs using generalized compressed suffix arrays, an extension of the FM-index idea from a string to a graph, with better accuracy than linear alignment and with reference bias largely removed 5. Giraffe, the later mapper, aligns to a collection of haplotypes threaded through the graph and maps to thousands of genomes at a speed comparable to mapping against one. That allowed genotyping of 167,000 structural variants in 5,202 short-read genomes 6.

The Human Pangenome Reference Consortium built the data to fill such a graph. The 2023 draft contains 47 phased diploid assemblies and adds 119 million base pairs of polymorphic sequence and 1,115 gene duplications relative to GRCh38. Used for short-read analysis, it reduced small-variant discovery errors by 34 percent and increased the structural variants detected per haplotype by 104 percent 7. The cost is that a coordinate stops being a single integer. “Position 12,345 on chromosome 1” becomes a statement about one path through the graph, and every file keyed by linear coordinates has to be projected onto that path. The graph represents variation the string could not, and it turns the interval structure of the previous section into a view rather than the truth.

Trees

Evolutionary history is a tree, and the question a tree answers is “where does this new sample belong, and what does it share with its neighbours”. During the SARS-CoV-2 pandemic the phylogeny grew faster than placement methods could keep up with. UShER addressed this with a mutation-annotated tree, a structure in which each branch carries the mutations inferred along it. A new sample is placed by walking the tree and comparing its mutation set to the branch labels rather than by realigning it against everything. Turakhia and colleagues report placement fast enough for real-time genomic contact tracing and a continuously updated reference phylogeny in the UCSC browser 8.

Population genetics has gone further and made the tree the storage format. Along a recombining genome the genealogy changes from one segment to the next, and a tree sequence stores that series of correlated trees compactly by sharing the nodes and edges that persist across neighbouring segments. Kelleher and colleagues describe an inference method for such genealogies that handles four orders of magnitude more sequences than earlier methods. They call the result an “evolutionary encoding” of the data, because statistics that would take a pass over the genotype matrix become efficient traversals of the trees 9. This is a strong version of the video’s claim: the tree is at once a compression, an index, and a query engine. The algorithms on it are short because the structure has done the work.

Sparse matrices

A single-cell RNA experiment produces a matrix of cells by genes, with tens of thousands of columns and up to millions of rows. Most entries are zero, because most genes are not detected in most cells. Stored dense, a million cells by twenty thousand genes at four bytes each is 80 GB. Stored as compressed sparse rows (CSR), where only the nonzero values and their column indices are kept, it is a small fraction of that. SCANPY’s authors report handling more than a million cells in Python this way 10. Alongside it they introduced AnnData, a container whose layout encodes the shape of the experiment. There is a central matrix X with an obs table of per-cell annotations and a var table of per-gene annotations. Beside them sit layers for alternative versions of X such as raw counts, obsm for per-cell embeddings, and obsp for cell-by-cell graphs such as neighbour networks. It is written to disk as HDF5 (the .h5ad file) or Zarr, and it can be opened in backed mode so that a dataset larger than memory is read a slice at a time 11.

The sparse layout is itself a choice about access. CSR makes it cheap to pull all genes for one cell and expensive to pull one gene across all cells, and compressed sparse columns (CSC) is the reverse. A per-cell normalization wants CSR and a per-gene test wants CSC. No layout makes both cheap, which is why converting between them is a routine pipeline step.

Chunked and columnar stores for population data

VCF stores one variant per row with every sample’s genotype and quality fields on that row. This is efficient when the question is “give me everything at this site” and inefficient when the question is “give me this one sample across all sites”, because the second requires reading and parsing every row. Czech and colleagues make this case for biobank-scale data, where hundreds of thousands of genomes produce hundreds of terabytes of compressed VCF. They propose VCF Zarr: the same data model encoded as chunked multidimensional arrays with variants and samples as separate dimensions. Compression and single-threaded performance are competitive with specialised genotype stores. The case studies are subsets of Genomics England (78,195 samples), Our Future Health (651,050 samples), and All of Us (245,394 samples) as well as about 4.5 million SARS-CoV-2 genomes 12. Each chunk is an independent object, so a query for one sample or one field reads only the chunks it touches.

This is the closest biological analogue to the “shape” in the sync-engine half of the video. A client rarely wants the whole dataset, and a chunked columnar layout makes the slice it does want the unit of I/O. The same idea underlies the tabix linear index and the AnnData backed mode: the designer decided in advance which slices would be common and laid the bytes out so those slices are contiguous.

What this means for one person’s longitudinal data

A personal molecular dataset is small by the standards above and awkward in a different way. It holds a few hundred measurements spread over years: blood panels at irregular intervals, a genome measured once, and glucose sampled every few minutes for a few weeks at a time. The tempting representation is a spreadsheet with a row per date and a column per analyte. That is the two-boolean version. It cannot hold two measurements of the same analyte on the same day by different methods. It equates a blank cell with “not measured” and with “below detection”. It has nowhere to put the unit or the assay, so those live in column headers and drift.

We would use a long table as the single source of truth with one row per measurement. Each row carries a timestamp, an analyte identifier, and a value with its unit. Each row also carries the assay or method and the specimen and report it came from. Missing values are absent rows rather than empty cells, and a value below the limit of detection is a row with a flag and the limit. Wide matrices for plotting or modelling are derived from this table and never edited. The genome stays in its own structures, a gVCF against a named reference build with its index, and joins to the table only through variant identifiers. A genotype is not a time series, and its linear coordinates are what every annotation resource is keyed on. The glucose trace is dense within a wearing period and absent between periods, so it is stored as one array per period with a start time and a sampling interval and summarised into the long table.

The last question is what an agent handed this data needs to be told. Left alone it will choose the convenient representation, which is the wide spreadsheet, and its tests will pass. So the shape has to be stated before any code is written. State the grain of the long table (one row per measurement of one analyte at one time by one method) and the meaning of an absent row versus a flagged one. State the rule that a value is never stored without a unit. State the genome’s reference build, and state that every wide table is derived. That is a short document, and reviewing it is worth more than reviewing the code generated from it. In biology as in ticketing, the structure is the part worth reading.

Questions people also ask

Why does bioinformatics have so many file formats? Each one is a data structure tuned for a question. FASTQ holds raw reads and BAM holds alignments sorted by coordinate. VCF holds variants across samples, BED and GTF hold intervals, and .h5ad holds annotated matrices.

Is a pangenome graph replacing the linear reference? It is used alongside it. Graph mapping reduces reference bias 7, but most annotation resources are keyed to linear coordinates, so results are projected back onto a reference path.

Should I store my own lab results as a spreadsheet? As a view, yes. As the source of truth, use a long table with one row per measurement carrying its unit and method, and derive the spreadsheet from it.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. P. Ferragina, G. Manzini. Opportunistic data structures with applications. Proceedings 41st Annual Symposium on Foundations of Computer Science, 2000. https://doi.org/10.1109/SFCS.2000.892127 ↩

  2. Heng Li, Richard Durbin. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics, 2009. https://doi.org/10.1093/bioinformatics/btp324 ↩

  3. Petr Danecek, Adam Auton, Goncalo Abecasis, et al. The variant call format and VCFtools. Bioinformatics, 2011. https://doi.org/10.1093/bioinformatics/btr330 ↩

  4. Heng Li. Tabix: fast retrieval of sequence features from generic TAB-delimited files. Bioinformatics, 2011. https://doi.org/10.1093/bioinformatics/btq671 ↩

  5. Erik Garrison, Jouni Sirén, Adam M. Novak, et al. Variation graph toolkit improves read mapping by representing genetic variation in the reference. Nature Biotechnology, 2018. https://doi.org/10.1038/nbt.4227 ↩

  6. Jouni Sirén, Jean Monlong, Xian Chang, et al. Pangenomics enables genotyping of known structural variants in 5202 diverse genomes. Science, 2021. https://doi.org/10.1126/science.abg8871 ↩

  7. Wen-Wei Liao, Mobin Asri, Jana Ebler, et al. A draft human pangenome reference. Nature, 2023. https://doi.org/10.1038/s41586-023-05896-x ↩ ↩2

  8. Yatish Turakhia, Bryan Thornlow, Angie S. Hinrichs, et al. Ultrafast Sample placement on Existing tRees (UShER) enables real-time phylogenetics for the SARS-CoV-2 pandemic. Nature Genetics, 2021. https://doi.org/10.1038/s41588-021-00862-7 ↩

  9. Jerome Kelleher, Yan Wong, Anthony W. Wohns, et al. Inferring whole-genome histories in large population datasets. Nature Genetics, 2019. https://doi.org/10.1038/s41588-019-0483-y ↩

  10. F. Alexander Wolf, Philipp Angerer, Fabian J. Theis. SCANPY: large-scale single-cell gene expression data analysis. Genome Biology, 2018. https://doi.org/10.1186/s13059-017-1382-0 ↩

  11. Isaac Virshup, Sergei Rybakov, Fabian J. Theis, et al. anndata: Access and store annotated data matrices. Journal of Open Source Software, 2024. https://doi.org/10.21105/joss.04371 ↩

  12. Eric Czech, Will Tyler, Tom White, et al. Analysis-ready VCF at Biobank scale using Zarr. GigaScience, 2025. https://doi.org/10.1093/gigascience/giaf049 ↩