Clay McLeod

Why we're building omics | Clay McLeodSkip to main content All postsomics began while I was writing the<br>chainfile crate, which translates<br>genomic locations between genome builds. I was building chainfile from scratch and<br>using rust-htslib and noodles as references. As I moved through their APIs, I kept<br>having to stop and answer a basic question: did this position name a nucleotide or a<br>boundary between nucleotides? The number looked valid in either system, and I had to<br>remember what it meant whenever code moved between them. I wanted that distinction to<br>survive function calls and refactors as part of the type, where the compiler could<br>check it. That requirement became the first part of omics.

omics is the foundation I wanted while building chainfile: a collection of Rust<br>crates that keeps biological meaning attached to coordinates, molecules, and<br>variation as those values move through a program. We are building it so those<br>conventions survive library boundaries and refactors without giving up performance.<br>Coordinates are where that work began, and the rest of this post uses them to show<br>what the principle looks like in code, how it performs, and where we plan to take it<br>next.

Follow omics on GitHub stjude-rust-labs/omics

Omics provides foundational abstractions for working with omics data in Rust. Follow<br>the project while its coordinate and variation APIs are still changing. Share design<br>feedback in issues or contribute code through pull requests.

View and starDiscuss in an issueContribute a PR

Coordinates from first principles

In genomics, a coordinate system is a convention for naming locations and spans along<br>a reference sequence. Systems can differ in where counting begins, whether a position<br>names a nucleotide or a boundary, and how interval endpoints are included. I will<br>leave the full treatment of these systems and their terminology to the<br>omics-coordinate documentation.<br>This shorter explanation is good enough for our purposes here (the docs call the<br>first system in-base, while its Rust marker is Base).

========================== seq0 =========================<br>• G • A • T • A • T • G • A •<br>║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║<br>║[--1--]║[--2--]║[--3--]║[--4--]║[--5--]║[--6--]║[--7--]║ In-base Positions<br>0 1 2 3 4 5 6 7 Interbase Positions

Base positions name nucleotides: 1 through 7 name GATATGA, while 0 is invalid.<br>Interbase positions name the boundaries between nucleotides: 0 falls before G,<br>1 between G and A, and 7 after the final A. Seven nucleotides therefore<br>cover eight interbase positions. One base position identifies a nucleotide; an<br>interbase representation needs the two adjacent boundaries.

Software expresses collections of contiguous nucleotides using intervals. A base<br>interval is often one-based and fully closed (GAT occupies [1, 3]). For an<br>interbase endpoint, inclusion determines whether the nucleotide immediately after<br>that boundary belongs to the interval. The same span is [0, 3) in an interbase<br>interval, which is often zero-based and half-open.

Both coordinate systems exist for good reason. Base positions tend to be easier for<br>people to reason about because they name nucleotides directly. Interbase positions<br>tend to be easier for software to work with: a zero-based, half-open interval has<br>length end - start, and adjacent intervals can share a boundary without<br>overlapping.<br>UCSC<br>displays one-based, fully closed coordinates in its web interface and stores<br>zero-based, half-open coordinates in database tables. Keeping the system with each<br>value prevents one-nucleotide endpoint errors.

How bioinformatics software represents coordinates today

Note: This is not an exhaustive survey. These examples show what can happen when a<br>tool or API documents its local coordinate convention but the value itself has no<br>common representation for that meaning.

bedtools provides a prime example of where<br>mixing coordinate systems can go wrong. In<br>issue #311, intersecting a BED<br>interval with a one-based VCF record at position 32 produced a zero-length BED<br>interval [32, 32) instead of the equivalent zero-based, half-open interval [31, 32). The<br>first fix<br>converted the VCF start, while related GFF/BED paths still mixed systems. The<br>follow-up refactor<br>gave record types an isZeroBased() method so output code could obtain the convention<br>from the model. The case shows that even widely used bioinformatics software can make<br>this mistake when coordinate-system meaning lives outside the value being handled.

rust-htslib

rust-htslib is the simplest case. Its current<br>Record::pos<br>documents the returned value as zero-based but exposes a plain i64:

pub fn pos(&self) -> i64

The zero-based value follows BAM's convention. The<br>SAM/BAM specification defines SAM<br>coordinates as one-based and BAM coordinates as zero-based. In a<br>comment quoted here, original<br>SAM/BAM author Heng Li explains that the human-readable formats were intended for<br>biologists, while BAM and BCF largely mirror the zero-based in-memory...