twobitreader

This crate provides fast DNA sequence extraction from 2bit files, a standard format in bioinformatics.

The motivation for this crate is speed. Extracting sequences is consistently faster than the best alternative. The focus is raw reading from 2bit, but fast concatenation and reverse-complement methods are also provided to make higher-level use cases easier.

Examples

Extracting sequences is straightforward:

let tbr = TwobitReader::open("hg38.2bit")?; // Human genome, build 38
let seq = tbr.get("chr1", 10000, 10005);    // -> String ("TAACC")

Concatenation works by iterating over (start, end) pairs. For example, assembling a spliced transcript:

// Exon ranges for human FKHL6 gene transcript (Gencode v43)
let exons = [(1389575, 1391118),             // Exon 1 (start, end)
             (1394695, 1395603)];            // Exon 2 (start, end)
let transcript = tbr.concat("chr6", &exons); // -> String

Parallelism is easy with crates like rayon. For example, batch extraction of sequences:

use rayon::prelude::*;
let args = [("chr1", 10000, 15000),
            ("chr1", 30000, 35000), /* ... */ ];
let seqs = args.into_par_iter()
    .map(|(chrom, start, end)| tbr.get(chrom, start, end))
    .collect::<Vec<_>>(); // -> Vec<String>

Or, assembling a batch of spliced transcripts in parallel:

use twobitreader::reverse_complement;
use rayon::prelude::*;

fn stranded(seq: String, strand: char) -> String {
    if strand == '+' { seq } else { reverse_complement(seq) }
}

let transcripts = [                   // (transcript_id, chromosome, exons)
    ("ENST00000407983.7", "chr2", '+', vec![(264899, 265007),     // Exon 1
                                            (271865, 271939),     // Exon 2
                                            (272036, 272557)]),   // Exon 3
    ("ENST00000319331.4", "chr3", '+', vec![(3799430, 3799919),   // Exon 1
                                            (3844363, 3849834)]), // Exon 2
    /* ... */
];
// Concatenate exons and then reverse-complement if necessary.
// (Correct if exons listed in genome-coordinate order.)
let seqs = transcripts.into_par_iter()
    .map(|(id, chrom, strand, exons)| (id, stranded(tbr.concat(chrom, exons), strand)))
    .collect::<HashMap<_, _>>();      // HashMap<&str, String>
let seq = &seqs["ENST00000407983.7"]; // -> &String to transcript sequence

Cold files are an order of magnitude slower to access than files already in memory ("hot"). Use prefetching to dramatically improve single-threaded speed:

let exons = [("chr1", 10000, 10200),
             ("chr1", 10500, 10700), /* ... */ ];
tbr.prefetch(&exons);             // Ask the operating system to start paging this data from disk.
let seqs = tbr.get_batch(&exons); // Access the memory as it arrives.

Speed

Two tasks were benchmarked:

• exons: extract 133,388 distinct human exon sequences;
• transcripts: concatenate 319,468 exons into 29,180 human spliced transcript sequences.

Speed depends on parallelism and page cache (hot vs cold):

• hot runs represent repeated or interactive dna extraction scenarios;
• cold runs represent a first run of a genomics pipeline, bound by disk speed;
• prefetch runs are cold but with a prefetch call preceding extraction.

The table below shows running times in milliseconds. Experimental details are BENCH.md.

Dependencies

• byteorder for handling endian-ness
• memmap2 for memory mapping the 2bit file
• seq-macro for generating 2bit decoder lookup table
• libc for prefetching file ranges on Apple targets
• windows-sys for prefetching file ranges on Windows targets

License: MIT OR Apache-2.0