implement sequence-read trait

Cargo.toml

@@ -27,7 +27,7 @@ bitflags = "0.9"
 regex = "1.0"
 linear-map = "1.2.0"
 serde = { version = "^1", optional = true }
-bio-types = ">=0.1.1"
+bio-types = ">=0.5.1"
 
 [features]
 default = ["bzip2", "lzma"]

src/bam/mod.rs

@@ -1499,7 +1499,7 @@ CCCCCCCCCCCCCCCCCCC"[..],
 
         // Load CRAM file, records
         let mut cram_reader = Reader::from_path(cram_path).unwrap();
-        cram_reader.set_reference(ref_path);
+        cram_reader.set_reference(ref_path).unwrap();
         let cram_records: Vec<Record> = cram_reader.records().map(|v| v.unwrap()).collect();
 
         // Load BAM file, records
@@ -1559,7 +1559,7 @@ CCCCCCCCCCCCCCCCCCC"[..],
             let mut cram_writer = Writer::from_cram_path(&cram_path, &header)
                 .ok()
                 .expect("Error opening CRAM file.");
-            cram_writer.set_reference(ref_path);
+            cram_writer.set_reference(ref_path).unwrap();
 
             // Write BAM records to CRAM file
             for rec in bam_records.iter() {
@@ -1574,7 +1574,7 @@ CCCCCCCCCCCCCCCCCCC"[..],
         {
             // Load written CRAM file
             let mut cram_reader = Reader::from_path(cram_path).unwrap();
-            cram_reader.set_reference(ref_path);
+            cram_reader.set_reference(ref_path).unwrap();
             let cram_records: Vec<Record> = cram_reader.records().map(|v| v.unwrap()).collect();
 
             // Compare CRAM records to BAM records

src/bam/record.rs

@@ -12,14 +12,14 @@ use std::str;
 use std::str::FromStr;
 use std::u32;
 
-use itertools::Itertools;
 use regex::Regex;
 
 use bam::{AuxWriteError, HeaderView, ReadError};
 use htslib;
 use utils;
 
 use bio_types::alignment::{Alignment, AlignmentMode, AlignmentOperation};
+use bio_types::sequence::SequenceRead;
 
 /// A macro creating methods for flag access.
 macro_rules! flag {
@@ -452,11 +452,15 @@ impl Record {
         self.inner().core.l_qseq as usize
     }
 
+    fn seq_data(&self) -> &[u8] {
+        let offset = self.qname_capacity() + self.cigar_len() * 4;
+        &self.data()[offset..][..(self.seq_len() + 1) / 2]
+    }
+
     /// Get read sequence. Complexity: O(1).
     pub fn seq(&self) -> Seq {
         Seq {
-            encoded: &self.data()[self.qname_capacity() + self.cigar_len() * 4..]
-                [..(self.seq_len() + 1) / 2],
+            encoded: self.seq_data(),
             len: self.seq_len(),
         }
     }
@@ -607,6 +611,24 @@ impl Drop for Record {
     }
 }
 
+impl SequenceRead for Record {
+    fn name(&self) -> &[u8] {
+        self.qname()
+    }
+
+    fn base(&self, i: usize) -> u8 {
+        decode_base(encoded_base(self.seq_data(), i))
+    }
+
+    fn base_qual(&self, i: usize) -> u8 {
+        self.qual()[i]
+    }
+
+    fn len(&self) -> usize {
+        self.seq_len()
+    }
+}
+
 /// Auxiliary record data.
 #[derive(Debug, Clone, Copy, PartialEq)]
 pub enum Aux<'a> {
@@ -665,6 +687,14 @@ static ENCODE_BASE: [u8; 256] = [
     15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
 ];
 
+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]
+}
+
 /// The sequence of a record.
 #[derive(Debug, Copy, Clone)]
 pub struct Seq<'a> {
@@ -676,7 +706,7 @@ impl<'a> Seq<'a> {
     /// Return encoded base. Complexity: O(1).
     #[inline]
     pub fn encoded_base(&self, i: usize) -> u8 {
-        (self.encoded[i / 2] >> ((!i & 1) << 2)) & 0b1111
+        encoded_base(self.encoded, i)
     }
 
     /// Return decoded sequence. Complexity: O(m) with m being the read length.

src/bcf/record.rs

@@ -8,7 +8,6 @@ use std::f32;
 use std::ffi;
 use std::fmt;
 use std::i32;
-use std::iter::repeat;
 use std::ptr;
 use std::rc::Rc;
 use std::slice;