Speed up sequence retrieval by avoiding bounds checking for decoded bases (impossible to fail). Offer entirely unchecked variant of base decoding.

Author: johanneskoester

src/bam/buffer.rs

@@ -4,6 +4,7 @@
 // except according to those terms.
 
 use std::collections::{vec_deque, VecDeque};
+use std::mem;
 use std::rc::Rc;
 use std::str;
 
@@ -23,6 +24,7 @@ pub struct RecordBuffer {
     overflow: Option<Rc<bam::Record>>,
     cache_cigar: bool,
     min_refetch_distance: u64,
+    buffer_record: Rc<bam::Record>,
 }
 
 unsafe impl Sync for RecordBuffer {}
@@ -42,6 +44,7 @@ impl RecordBuffer {
             overflow: None,
             cache_cigar,
             min_refetch_distance: 1,
+            buffer_record: Rc::new(bam::Record::new()),
         }
     }
 
@@ -107,26 +110,29 @@ impl RecordBuffer {
 
             // extend to the right
             loop {
-                let mut record = Rc::new(bam::Record::new());
-                if self
+                match self
                     .reader
-                    .read(Rc::get_mut(&mut record).unwrap())
-                    .is_none()
+                    .read(Rc::get_mut(&mut self.buffer_record).unwrap())
                 {
-                    break;
+                    None => break,
+                    Some(res) => res?,
                 }
 
-                if record.is_unmapped() {
+                if self.buffer_record.is_unmapped() {
                     continue;
                 }
 
-                let pos = record.pos();
+                let pos = self.buffer_record.pos();
 
                 // skip records before the start
                 if pos < start as i64 {
                     continue;
                 }
 
+                // Record is kept, do not reuse it for next iteration
+                // and thus create a new one.
+                let mut record = mem::replace(&mut self.buffer_record, Rc::new(bam::Record::new()));
+
                 if self.cache_cigar {
                     Rc::get_mut(&mut record).unwrap().cache_cigar();
                 }

src/bam/record.rs

@@ -717,7 +717,7 @@ impl SequenceRead for Record {
     }
 
     fn base(&self, i: usize) -> u8 {
-        decode_base(encoded_base(self.seq_data(), i))
+        *decode_base_unchecked(encoded_base(self.seq_data(), i))
     }
 
     fn base_qual(&self, i: usize) -> u8 {
@@ -820,12 +820,19 @@ static ENCODE_BASE: [u8; 256] = [
     15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
 ];
 
+#[inline]
 fn encoded_base(encoded_seq: &[u8], i: usize) -> u8 {
     (encoded_seq[i / 2] >> ((!i & 1) << 2)) & 0b1111
 }
 
-fn decode_base(base: u8) -> u8 {
-    DECODE_BASE[base as usize]
+#[inline]
+unsafe fn encoded_base_unchecked(encoded_seq: &[u8], i: usize) -> u8 {
+    (encoded_seq.get_unchecked(i / 2) >> ((!i & 1) << 2)) & 0b1111
+}
+
+#[inline]
+fn decode_base_unchecked(base: u8) -> &'static u8 {
+    unsafe { DECODE_BASE.get_unchecked(base as usize) }
 }
 
 /// The sequence of a record.
@@ -842,6 +849,20 @@ impl<'a> Seq<'a> {
         encoded_base(self.encoded, i)
     }
 
+    /// Return encoded base. Complexity: O(1).
+    #[inline]
+    pub unsafe fn encoded_base_unchecked(&self, i: usize) -> u8 {
+        encoded_base_unchecked(self.encoded, i)
+    }
+
+    /// Obtain decoded base without performing bounds checking.
+    /// Use index based access seq()[i], for checked, safe access.
+    /// Complexity: O(1).
+    #[inline]
+    pub unsafe fn decoded_base_unchecked(&self, i: usize) -> u8 {
+        *decode_base_unchecked(self.encoded_base_unchecked(i))
+    }
+
     /// Return decoded sequence. Complexity: O(m) with m being the read length.
     pub fn as_bytes(&self) -> Vec<u8> {
         (0..self.len()).map(|i| self[i]).collect()
@@ -862,7 +883,7 @@ impl<'a> ops::Index<usize> for Seq<'a> {
 
     /// Return decoded base at given position within read. Complexity: O(1).
     fn index(&self, index: usize) -> &u8 {
-        &DECODE_BASE[self.encoded_base(index) as usize]
+        decode_base_unchecked(self.encoded_base(index))
     }
 }
 

src/bcf/mod.rs

@@ -1327,7 +1327,7 @@ mod tests {
         let mut reader = Reader::from_path("test/test_svlen.vcf").unwrap();
 
         let mut record = reader.empty_record();
-        reader.read(&mut record).unwrap();
+        reader.read(&mut record).unwrap().unwrap();
 
         assert_eq!(record.info(b"SVLEN").integer().unwrap(), Some(&[-127][..]));
     }
@@ -1366,7 +1366,7 @@ mod tests {
     fn test_alt_allele_dosage() {
         let path = &"test/test_string.vcf";
         let mut bcf = Reader::from_path(path).expect("Error opening file.");
-        let header = bcf.header();
+        let _header = bcf.header();
         // FORMAT fields of first record of the vcf should look like:
         // GT:FS1:FN1	./1:LongString1:1	1/1:ss1:2
         let mut first_record = bcf.records().next().unwrap().expect("Fail to read record");