omics began while I was writing the
chainfile crate, which translates
genomic locations between genome builds. I was building chainfile from scratch and
using rust-htslib and noodles as references. As I moved through their APIs, I kept
having to stop and answer a basic question: did this position name a nucleotide or a
boundary between nucleotides? The number looked valid in either system, and I had to
remember what it meant whenever code moved between them. I wanted that distinction to
survive function calls and refactors as part of the type, where the compiler could
check it. That requirement became the first part of omics.
omics is the foundation I wanted while building chainfile: a collection of Rust
crates that keeps biological meaning attached to coordinates, molecules, and
variation as those values move through a program. We are building it so those
conventions survive library boundaries and refactors without giving up performance.
Coordinates are where that work began, and the rest of this post uses them to show
what the principle looks like in code, how it performs, and where we plan to take it
next.
Coordinates from first principles
In genomics, a coordinate system is a convention for naming locations and spans along
a reference sequence. Systems can differ in where counting begins, whether a position
names a nucleotide or a boundary, and how interval endpoints are included. I will
leave the full treatment of these systems and their terminology to the
omics-coordinate documentation.
This shorter explanation is good enough for our purposes here (the docs call the
first system in-base, while its Rust marker is Base).
========================== seq0 =========================
• G • A • T • A • T • G • A •
║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║
║[--1--]║[--2--]║[--3--]║[--4--]║[--5--]║[--6--]║[--7--]║ In-base Positions
0 1 2 3 4 5 6 7 Interbase Positions
Base positions name nucleotides: 1 through 7 name GATATGA, while 0 is invalid.
Interbase positions name the boundaries between nucleotides: 0 falls before G,
1 between G and A, and 7 after the final A. Seven nucleotides therefore
cover eight interbase positions. One base position identifies a nucleotide; an
interbase representation needs the two adjacent boundaries.
Software expresses collections of contiguous nucleotides using intervals. A base
interval is often one-based and fully closed (GAT occupies [1, 3]). For an
interbase endpoint, inclusion determines whether the nucleotide immediately after
that boundary belongs to the interval. The same span is [0, 3) in an interbase
interval, which is often zero-based and half-open.
Both coordinate systems exist for good reason. Base positions tend to be easier for
people to reason about because they name nucleotides directly. Interbase positions
tend to be easier for software to work with: a zero-based, half-open interval has
length end - start, and adjacent intervals can share a boundary without
overlapping.
UCSC
displays one-based, fully closed coordinates in its web interface and stores
zero-based, half-open coordinates in database tables. Keeping the system with each
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 tool or API documents its local coordinate convention but the value itself has no common representation for that meaning.
bedtools provides a prime example of where
mixing coordinate systems can go wrong. In
issue #311, intersecting a BED
interval with a one-based VCF record at position 32 produced a zero-length BED
interval [32, 32) instead of the equivalent zero-based, half-open interval [31, 32). The
first fix
converted the VCF start, while related GFF/BED paths still mixed systems. The
follow-up refactor
gave record types an isZeroBased() method so output code could obtain the convention
from the model. The case shows that even widely used bioinformatics software can make
this mistake when coordinate-system meaning lives outside the value being handled.
rust-htslib
rust-htslib is the simplest case. Its current
Record::pos
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
SAM/BAM specification defines SAM
coordinates as one-based and BAM coordinates as zero-based. In a
comment quoted here, original
SAM/BAM author Heng Li explains that the human-readable formats were intended for
biologists, while BAM and BCF largely mirror the zero-based in-memory representations
used by programmers. rust-htslib preserves that deliberate choice, but the type
retains only a signed integer, so a receiving function must recover the convention
from documentation or surrounding code.
noodles
noodles-core represents positions with a
Position
that is nonzero and one-based. This makes the indexing convention part of the type.
noodles-vcf still needs record-level logic to decide how far a variant extends:
variant_end
derives a version-aware end from combinations of END, REF, SVLEN, or sample
LEN, and the
regional query implementation
turns that result into an inclusive start..=end interval for intersection testing.
A VCF breakend also needs a remote coordinate and orientation.
noodles-vcf keeps alternate alleles as String values,
so both remain inside VCF ALT rather than a typed adjacency. Together, the two cases
show where a mixed representation would help: the inclusive query span covers
reference bases, while the breakend marks an oriented boundary between them.
rust-htslib and noodles, two prominent Rust libraries for working with
next-generation sequencing data, use different coordinate conventions, and those
differences still matter downstream. Passing the one-based, inclusive VCF span that
noodles
derives
into a
coitrees::IntervalNode
works because it
also uses inclusive endpoints,
but the values lose their one-based origin; conversion to a zero-based, half-open BED
interval remains the caller's responsibility. Rust's strong type system should let
us represent both conventions cleanly and make conversion explicit at library
boundaries.
A shared foundation for coordinates
The coordinate layer in omics is meant to do that: carry coordinate meaning
through Rust code until an API boundary requires a different representation.
Position<S> is the smallest piece of that design. It attaches the coordinate system
to the integer:
use omics::coordinate::Position;
use omics::coordinate::system::Base;
use omics::coordinate::system::Interbase;
let base = Position::<Base>::try_new(8)?;
let interbase = Position::<Interbase>::new(8);
Both values store 8; Position<Base> names a nucleotide, while
Position<Interbase> names a boundary. Passing one to a function that accepts the
other is a compiler error without a deliberate conversion.
Position::<Base>::try_new(0) returns an error, while
Position::<Interbase>::new(0) creates the valid boundary before the first nucleotide.
When a function preserves S in its signature, it also preserves the coordinate
system. Base position 0 fails at construction.
A position identifies a location within a coordinate system, but not the contig or
strand where it occurs. Coordinate<S> adds that context by combining a contig and
strand with a Position<S>. The system parameter remains part of the resulting type.
An interval describes the span between two coordinates. Interval<S> accepts
endpoints from the same coordinate system, requires matching contigs and strands, and
verifies strand-aware order: start cannot exceed end on +, while end cannot exceed
start on -. One checked value retains the convention, contig, strand, and
orientation.
Unlike interval representations that always order endpoints by reference position,
omics keeps them in the molecule's 5′-to-3′ direction. Coordinates ascend on the
positive strand but descend on the negative strand, so chr1:-:20-10 is ordered
from start to end even though its numeric values decrease.
A system conversion still requires a scientific decision. For a base coordinate at
10, the caller must select an adjacent boundary relative to the strand. On +,
nudge_forward() produces interbase 10 and nudge_backward() produces 9, whereas
on - the results reverse. The methods name movement along the molecule and return
None on overflow or base position 0, avoiding unchecked + 1 or - 1
arithmetic.
Applying the foundation in chainfile
chainfile shows how the model works when a format already defines its coordinate
system. UCSC defines chain alignment endpoints as zero-based, half-open
intervals, which correspond
directly to Interval<Interbase>. The crate
parses each header span into that type,
advances interbase coordinates through every aligned block and gap,
and accepts the same type in
Machine::liftover().
Its internal representation therefore matches the format during the calculation,
leaving conversion to code that actually needs another coordinate system.
The
liftover example
illustrates how the typed interval passes directly into the liftover machine:
use std::{fs::File, io::BufReader};
use chain::liftover;
use chainfile as chain;
use flate2::read::GzDecoder;
use omics::coordinate::interval::interbase::Interval;
let interval = "chr1:+:1000000-1000100".parse::<Interval>()?;
let reader = File::open("hg38ToHg19.over.chain.gz")
.map(GzDecoder::new)
.map(BufReader::new)
.map(chain::Reader::new)?;
let machine = liftover::machine::builder::Builder.try_build_from(reader)?;
let results = machine.liftover(interval);
A conversion between noodles and omics occurs in the
chain_check example.
At the FASTA lookup, Sequence::slice() receives an inclusive
noodles_core::region::Interval whose
Position endpoints are one-based.
parse_interval() first calls into_equivalent_base(), then orders the endpoints by
reference position: start then end on +, end then start on -. Here the code
constructs the inclusive start..=end range only inside parse_interval(), and a
library conversion should remain at the function that crosses that API boundary
rather than spread through the liftover calculation. This is what a shared coordinate
foundation should allow: each library uses the representation its API expects, and
the conversion remains explicit.
Correctness without sacrificing performance
A coordinate model can represent the biology correctly and still impose enough
overhead to be impractical. I expect each abstraction to add minimal storage or
runtime cost. Coordinate-system information can exist entirely in Rust's type system
and disappear during compilation; Base and Interbase should therefore add no
storage to a Position. Constructing an owned contig or checking interval endpoints
can require runtime work. We enforce layout properties at compile time and
benchmark the operations that remain.
The position benchmarks compare typed operations with their raw u32 counterparts,
and on my machine, typed position construction, base validation, and checked addition
all remained within about 4% of the corresponding raw operations. The compiler can
distinguish the coordinate systems for little measurable cost.
Memory benchmarking immediately showed its value for intervals. By inspecting the representation, I recognized that each interval stored its contig and strand twice, once in each owned endpoint, and changed it to store each only once.
I then benchmarked construction separately because memory use and construction time
measure different costs. Interval::try_new was about 7% slower than a raw-endpoint
control that mirrors the same validation and structured error paths.
Upon examining the assembly for that remaining difference, I recognized that the
typed path keeps more intermediate state on the stack and masks the Strand value
before branching, likely accounting for the remaining overhead, which I considered
reasonable for the validation and type information the operation preserves.
Performance will remain a key design objective for omics as these representations
expand.
Moving forward
With the coordinate model in place, we are considering other representations that might naturally extend from it. Variation is one of them. A variant's coordinate system must agree with what the alteration does to the reference and alternate sequences, so it is a useful test of whether the same approach can prevent contradictory descriptions of the event.
We have started exploring that question in recent work on small variants, structural variants, and copy-number variants. I will leave the implementation details in those pull requests; the broader question is which biological relationships should remain attached to the values that move through a program.
An invitation to test the design
Coordinates and variants now give omics enough of a foundation to evaluate the
design as a system. I want community feedback now, while its APIs can still change
before downstream projects depend on them.
Next, we plan to continue work on copy-number variants and begin work on deeper variation normalization and an immutable interval index. None has a release date yet.
The quickest contribution is an API review. If a type encodes the wrong biological
distinction or leaves one implicit, open an
issue. Trying to integrate
omics into an existing project would
be even better. Real input and output boundaries can show which abstractions hold and
which need to change.
In my mind, success means seeing developers adopt these shared types across projects
and share responsibility for their direction and maintenance. I hope omics can
become infrastructure that the Rust bioinformatics community builds and owns
together.