Expand description
§bitnuc
A library for efficient nucleotide sequence manipulation using 2-bit encoding.
§Summary
This a SIMD-accelerated two-bit encoding library for nucleotide sequences. It is meant to have fast encode/decode routines for small and large sequences and to provide a fairly unstructured interface for working with nucleotide sequences in memory.
It provides:
- 2-bit nucleotide encoding (A=00, C=01, G=10, T=11) with SIMD dispatched at runtime via
fearless_simd - A little-endian-pinned
u64kmer boundary (as_2bit/from_2bit) for hashing and fixed-width integer storage of sequences up to 32 bases - Ambiguous base detection (
ambiguous_bases) for tracking non-ACGTacgtbases - Extraction routines for extracting aligned and unaligned ranges from packed sequences
§Encoded Format
Sequences encode to plain bytes: byte k holds bases 4k..4k+4, base j at
bits 2*(j % 4) of its byte. A sequence of n bases occupies exactly
n.div_ceil(4) bytes and trailing pad bits are always zero.
Input Sequence: [A][C][G][T]
Output Buffer: [11 10 01 00]
[]: byte boundary
Bit pairs are MSB-first: base j occupies bits 2j..2j+1,
so the first base is the rightmost pair (A=00).Note: Encoding is lossy for bytes outside
ACGTacgt: invalid bases map to an unspecified code rather than an error. If you need to preserve ambiguous bases, detect and track them separately (seeambiguous_bases).
§Encoding and Decoding
The core functions operate on byte slices, with _resize variants that
manage the buffer length for you:
use bitnuc::{encode_resize, decode_resize};
fn main() -> Result<(), bitnuc::BitnucError> {
let seq = b"ACGTACGTAC"; // 10 bases -> 3 encoded bytes
let mut ebuf = Vec::new();
encode_resize(seq, &mut ebuf);
assert_eq!(ebuf.len(), 3);
let mut dbuf = Vec::new();
decode_resize(&ebuf, seq.len(), &mut dbuf)?;
assert_eq!(&dbuf, seq);
Ok(())
}The slice-based variants write into caller-provided buffers, which lets consumers control allocation and padding (e.g. file formats that pad encoded sequences to 8-byte words):
use bitnuc::{encode, decode};
fn main() -> Result<(), bitnuc::BitnucError> {
let seq = b"ACGTACGTAC";
// Pad the encoded buffer to an 8-byte multiple: the layout is identical
// to the legacy u64 packing serialized little-endian
let mut ebuf = vec![0u8; seq.len().div_ceil(4).next_multiple_of(8)];
encode(seq, &mut ebuf)?;
let mut dbuf = vec![0u8; seq.len()];
decode(&ebuf, seq.len(), &mut dbuf)?;
assert_eq!(&dbuf, seq);
Ok(())
}§u64 Kmer Packing
For hashing and fixed-width storage of short sequences (barcodes, UMIs,
k-mers up to 32 bases), as_2bit and from_2bit pack to and from a u64.
Note: These are pinned to little-endian internally.
use bitnuc::{as_2bit, from_2bit};
use std::collections::HashMap;
fn main() -> Result<(), bitnuc::BitnucError> {
let packed = as_2bit(b"ACGT")?;
assert_eq!(packed, 0b11100100);
// Efficient k-mer counting
let mut kmer_counts = HashMap::new();
for window in b"ACGTACGT".windows(4) {
*kmer_counts.entry(as_2bit(window)?).or_insert(0) += 1;
}
assert_eq!(kmer_counts.get(&packed), Some(&2));
// Unpacking returns a stack array of all 32 bases; slice to your length
let unpacked = from_2bit(packed);
assert_eq!(&unpacked[..4], b"ACGT");
Ok(())
}Packed kmers can be compared with hdist_scalar:
use bitnuc::{as_2bit, hdist_scalar};
fn main() -> Result<(), bitnuc::BitnucError> {
let u = as_2bit(b"ACGT")?;
let v = as_2bit(b"ACGA")?;
assert_eq!(hdist_scalar(u, v, 4)?, 1);
Ok(())
}§Identifying ambiguous bases
Ambiguous bases (non-ACGT) bases are unable to be represented with this two-bit encoding scheme.
The encoding algorithm also remaps lowercase to upper case (so acgt -> ACGT internally).
Use ambiguous_bases to track all positions of unrepresentable nucleotides.
The position buffer is generic over its element type (usize, u64, or u32):
use bitnuc::ambiguous_bases;
fn main() {
let seq = b"ACgTNACYAaTH"; // has unrepresentable bases (N/Y/H)
let mut pos: Vec<usize> = Vec::default();
ambiguous_bases(seq, &mut pos);
assert_eq!(
pos,
vec![4, 7, 11],
);
}Note:
ambiguous_basesonly tracks non-ACGTacgtbases. It does not identify lowercase letters which are also not representable but which are remapped to their uppercase variants through encoding/decoding.
§Extraction
It is oftentimes useful to extract subsequences from packed sequences, e.g. for hashing or extracting ranges.
Use extract to pull out ranges of packed bytes into an aligned buffer.
The range is specified in the basepair range (e.g. 0..10 for the first 10 bases) which is converted to a byte range internally.
The bits relevant to the range are extracted and packed into the output buffer so that the output buffer is correctly aligned.
use bitnuc::{encode_resize, extract_resize, decode_resize, BitnucError};
fn main() -> Result<(), BitnucError> {
let seq = b"ACCAAGGTTACATGAAGTTAACCAAGAGAC";
// encode the sequence
let mut packed = Vec::new();
encode_resize(seq, &mut packed);
// extract some subsequence
let range = 4..11;
let mut extracted = Vec::new();
extract_resize(&packed, range.clone(), &mut extracted)?;
// decode the extracted sequence
let mut decoded = Vec::new();
decode_resize(&extracted, range.len(), &mut decoded)?;
// validate expected output
assert_eq!(decoded, seq[range]);
Ok(())
}§Memory Usage
The 2-bit encoding provides significant memory savings:
Standard encoding: 1 byte per base
ACGT = 4 bytes = 32 bits
2-bit encoding: 2 bits per base
ACGT = 1 byte = 8 bits§Performance
Throughput by sequence length, measured on an Apple M3 Pro with
target-cpu=native (criterion mean, 1 byte per base):
| bp | encode (GB/s) | decode (GB/s) |
|---|---|---|
| 10 | 2.3 | 2.2 |
| 100 | 12.5 | 13.8 |
| 1000 | 35.4 | 29.7 |
| 10000 | 38.2 | 34.0 |
| 100000 | 38.6 | 33.5 |
To regenerate the table on your machine:
RUSTFLAGS="-C target-cpu=native" cargo bench --bench simd_comparison -- coding_2bit
uv run scripts/perf_table.py§SIMD Acceleration
The 2-bit encode and decode are SIMD accelerated via
fearless_simd, with the instruction set
(NEON, SSE, AVX2, AVX-512) selected at runtime.
§Related Work
I highly recommend checking out packed-seq. They are currently the highest performance 2-bit encoding library in Rust as far as I can tell. They follow a different bit-packing scheme than this library and they can shave off a few instructions in their SIMD routines.
If you’re interested in 2-bit encoding in general make sure to check out cute-nucleotides which has an excellent overview of different algorithms and their performance characteristics.
Enums§
Traits§
- Pos
- A position element type that a sequence index can be losslessly stored as.
Functions§
- ambiguous_
bases - Identifies all ambiguous bases (non-ACGTacgt) and appends their positions to a user-supplied buffer (generic index type).
- as_2bit
- Packs a sequence of up to 32 bases into a
u64, 2 bits per base. - decode
- Decodes a 2-bit encoded buffer back into
nASCII nucleotides. - decode_
resize - Decodes a 2-bit encoded buffer back into
nASCII nucleotides with a resizableseqbuffer. - encode
- Two-bit encodes an input sequence into an encoding buffer
- encode_
resize - Two-bit encodes an input sequence into an encoding buffer which can be resized.
- extract
- extracts the bits belonging to the given basepair
rangefrompackedintointo - extract_
resize - from_
2bit - Unpacks a 2-bit packed
u64into a stack array of 32 ASCII bases. - hdist_
scalar - Calculates the hamming distance between two 2-bit packed
u64kmers (as produced byas_2bit) of lengthlenbases.