implement various bio-types (interval, locus, strand)

Author: johanneskoester

Cargo.toml

@@ -28,7 +28,7 @@ regex = "1.3"
 linear-map = "1.2"
 serde_base = { version = "^1", optional = true, package = "serde" }
 serde_bytes = { version = "0.11", optional = true }
-bio-types = ">=0.5.1"
+bio-types = ">=0.6"
 snafu = ">= 0.5.0, <= 0.6.0"
 hts-sys = { version = "^1.10", path = "hts-sys", default-features = false }
 

src/bam/mod.rs

@@ -20,6 +20,7 @@ use std::ffi;
 use std::path::Path;
 use std::slice;
 use std::str;
+
 use url::Url;
 
 use crate::htslib;
@@ -268,7 +269,11 @@ impl Read for Reader {
             -1 => Ok(false),
             -2 => Err(Error::TruncatedRecord),
             -4 => Err(Error::InvalidRecord),
-            _ => Ok(true),
+            _ => {
+                record.set_header(self.header().clone());
+
+                Ok(true)
+            }
         }
     }
 
@@ -458,7 +463,11 @@ impl Read for IndexedReader {
                     -1 => Ok(false),
                     -2 => Err(Error::TruncatedRecord),
                     -4 => Err(Error::InvalidRecord),
-                    _ => Ok(true),
+                    _ => {
+                        record.set_header(self.header().clone());
+
+                        Ok(true)
+                    }
                 }
             }
             None => Ok(false),
@@ -862,6 +871,10 @@ impl HeaderView {
         }
     }
 
+    pub fn tid2name(&self, tid: u32) -> &[u8] {
+        unsafe { ffi::CStr::from_ptr(htslib::sam_hdr_tid2name(self.inner, tid as i32)).to_bytes() }
+    }
+
     pub fn target_count(&self) -> u32 {
         self.inner().n_targets as u32
     }

src/bam/record.rs

@@ -22,7 +22,9 @@ use crate::htslib;
 use crate::utils;
 
 use bio_types::alignment::{Alignment, AlignmentMode, AlignmentOperation};
+use bio_types::genome;
 use bio_types::sequence::SequenceRead;
+use bio_types::strand::ReqStrand;
 
 /// A macro creating methods for flag access.
 macro_rules! flag {
@@ -46,6 +48,7 @@ pub struct Record {
     pub inner: htslib::bam1_t,
     own: bool,
     cigar: Option<CigarStringView>,
+    header: Option<HeaderView>,
 }
 
 unsafe impl Send for Record {}
@@ -109,6 +112,7 @@ impl Record {
             inner: unsafe { MaybeUninit::zeroed().assume_init() },
             own: true,
             cigar: None,
+            header: None,
         }
     }
 
@@ -128,6 +132,7 @@ impl Record {
             },
             own: false,
             cigar: None,
+            header: None,
         }
     }
 
@@ -162,6 +167,10 @@ impl Record {
         }
     }
 
+    pub fn set_header(&mut self, header: HeaderView) {
+        self.header = Some(header);
+    }
+
     pub(super) fn data(&self) -> &[u8] {
         unsafe { slice::from_raw_parts(self.inner().data, self.inner().l_data as usize) }
     }
@@ -224,6 +233,16 @@ impl Record {
         self.inner_mut().core.qual = mapq;
     }
 
+    /// Get strand information from record flags.
+    pub fn strand(&mut self) -> ReqStrand {
+        let reverse = self.flags() & 0x10 != 0;
+        if reverse {
+            ReqStrand::Reverse
+        } else {
+            ReqStrand::Forward
+        }
+    }
+
     /// Get raw flags.
     pub fn flags(&self) -> u16 {
         self.inner().core.flag
@@ -706,6 +725,39 @@ impl SequenceRead for Record {
     }
 }
 
+impl genome::AbstractInterval for Record {
+    /// Return contig name. Panics if record does not know its header (which happens if it has not been read from a file).
+    fn contig(&self) -> &str {
+        let tid = self.tid();
+        if tid < 0 {
+            panic!("invalid tid, must be at least zero");
+        }
+        str::from_utf8(
+            self.header
+                .as_ref()
+                .expect(
+                    "header must be set (this is the case if the record has been read from a file)",
+                )
+                .tid2name(tid as u32),
+        )
+        .expect("unable to interpret contig name as UTF-8")
+    }
+
+    /// Return genomic range covered by alignment. Panics if `Record::cache_cigar()` has not been called first or `Record::pos()` is less than zero.
+    fn range(&self) -> ops::Range<genome::Position> {
+        let end_pos = self
+            .cigar_cached()
+            .expect("cigar has not been cached yet, call cache_cigar() first")
+            .end_pos() as u64;
+
+        if self.pos() < 0 {
+            panic!("invalid position, must be positive")
+        }
+
+        self.pos() as u64..end_pos
+    }
+}
+
 /// Auxiliary record data.
 #[derive(Debug, Clone, Copy, PartialEq)]
 pub enum Aux<'a> {
@@ -1091,6 +1143,28 @@ impl CigarStringView {
         pos
     }
 
+    /// Get number of bases softclipped at the beginning of the alignment.
+    pub fn leading_softclips(&self) -> i64 {
+        self.first().map_or(0, |cigar| {
+            if let Cigar::SoftClip(s) = cigar {
+                *s as i64
+            } else {
+                0
+            }
+        })
+    }
+
+    /// Get number of bases softclipped at the end of the alignment.
+    pub fn trailing_softclips(&self) -> i64 {
+        self.last().map_or(0, |cigar| {
+            if let Cigar::SoftClip(s) = cigar {
+                *s as i64
+            } else {
+                0
+            }
+        })
+    }
+
     /// For a given position in the reference, get corresponding position within read.
     /// If reference position is outside of the read alignment, return None.
     ///

src/bcf/buffer.rs

@@ -52,14 +52,14 @@ impl RecordBuffer {
         self.ringbuffer2.clear();
     }
 
-    fn drain_left(&mut self, rid: u32, window_start: u32) -> usize {
+    fn drain_left(&mut self, rid: u32, window_start: u64) -> usize {
         // remove records too far left or from wrong rid
         // rec.rid() will always yield Some(), because otherwise we won't put the rec into the
         // buffer.
         let to_remove = self
             .ringbuffer
             .iter()
-            .take_while(|rec| rec.pos() < window_start || rec.rid().unwrap() != rid)
+            .take_while(|rec| (rec.pos() as u64) < window_start || rec.rid().unwrap() != rid)
             .count();
         self.ringbuffer.drain(..to_remove);
         to_remove
@@ -69,7 +69,7 @@ impl RecordBuffer {
     /// the start coordinate of any previous `fill` operation.
     /// Coordinates are 0-based, and end is exclusive.
     /// Returns tuple with numbers of added and deleted records compared to previous fetch.
-    pub fn fetch(&mut self, chrom: &[u8], start: u32, end: u32) -> Result<(usize, usize)> {
+    pub fn fetch(&mut self, chrom: &[u8], start: u64, end: u64) -> Result<(usize, usize)> {
         // TODO panic if start is left of previous start or we have moved past the given chrom
         // before.
         let rid = self.reader.header.name2rid(chrom)?;
@@ -107,7 +107,7 @@ impl RecordBuffer {
 
         // move overflow from last fill into ringbuffer
         if self.overflow.is_some() {
-            let pos = self.overflow.as_ref().unwrap().pos();
+            let pos = self.overflow.as_ref().unwrap().pos() as u64;
             if pos >= start {
                 if pos <= end {
                     self.ringbuffer.push_back(self.overflow.take().unwrap());
@@ -129,7 +129,7 @@ impl RecordBuffer {
                 // EOF
                 break;
             }
-            let pos = rec.pos();
+            let pos = rec.pos() as u64;
             if let Some(rec_rid) = rec.rid() {
                 if rec_rid == rid {
                     if pos >= end {

src/bcf/errors.rs

@@ -16,7 +16,7 @@ pub enum Error {
     #[snafu(display("error setting threads for writing BCF/VCF file(s)"))]
     SetThreads,
     #[snafu(display("error seeking to {}:{} in indexed BCF/VCF file", contig, start))]
-    Seek { contig: String, start: u32 },
+    Seek { contig: String, start: u64 },
     #[snafu(display("error writing record to BCF/VCF file"))]
     Write,
     #[snafu(display("tag {} undefined in BCF/VCF header", tag))]

src/bcf/mod.rs

@@ -166,7 +166,7 @@ pub struct IndexedReader {
     header: Rc<HeaderView>,
 
     /// The position of the previous fetch, if any.
-    current_region: Option<(u32, u32, u32)>,
+    current_region: Option<(u32, u64, u64)>,
 }
 
 unsafe impl Send for IndexedReader {}
@@ -225,7 +225,7 @@ impl IndexedReader {
     ///           contig name to ID.
     /// * `start` - `0`-based start coordinate of region on reference.
     /// * `end` - `0`-based end coordinate of region on reference.
-    pub fn fetch(&mut self, rid: u32, start: u32, end: u32) -> Result<()> {
+    pub fn fetch(&mut self, rid: u32, start: u64, end: u64) -> Result<()> {
         let contig = self.header.rid2name(rid).unwrap();
         let contig = ffi::CString::new(contig).unwrap();
         if unsafe { htslib::bcf_sr_seek(self.inner, contig.as_ptr(), start as i64) } != 0 {
@@ -269,7 +269,7 @@ impl Read for IndexedReader {
                     Some((rid, _start, end)) => {
                         if record.rid().is_some()
                             && rid == record.rid().unwrap()
-                            && record.pos() <= end
+                            && record.pos() as u64 <= end
                         {
                             Ok(true)
                         } else {
@@ -349,7 +349,7 @@ pub mod synced {
         headers: Vec<Rc<HeaderView>>,
 
         /// The position of the previous fetch, if any.
-        current_region: Option<(u32, u32, u32)>,
+        current_region: Option<(u32, u64, u64)>,
     }
 
     // TODO: add interface for setting threads, ensure that the pool is freed properly
@@ -500,7 +500,7 @@ pub mod synced {
         ///           contig name to ID.
         /// * `start` - `0`-based start coordinate of region on reference.
         /// * `end` - `0`-based end coordinate of region on reference.
-        pub fn fetch(&mut self, rid: u32, start: u32, end: u32) -> Result<()> {
+        pub fn fetch(&mut self, rid: u32, start: u64, end: u64) -> Result<()> {
             let contig = {
                 let contig = self.header(0).rid2name(rid).unwrap(); //.clone();
                 ffi::CString::new(contig).unwrap()

src/bcf/record.rs

@@ -13,6 +13,7 @@ use std::rc::Rc;
 use std::slice;
 use std::str;
 
+use bio_types::genome;
 use ieee754::Ieee754;
 use itertools::Itertools;
 
@@ -160,8 +161,8 @@ impl Record {
     }
 
     // Return 0-based position.
-    pub fn pos(&self) -> u32 {
-        self.inner().pos as u32
+    pub fn pos(&self) -> i64 {
+        self.inner().pos
     }
 
     /// Set 0-based position.
@@ -622,6 +623,21 @@ impl Record {
     }
 }
 
+impl genome::AbstractLocus for Record {
+    fn contig(&self) -> &str {
+        str::from_utf8(
+            self.header()
+                .rid2name(self.rid().expect("rid not set"))
+                .expect("unable to find rid in header"),
+        )
+        .expect("unable to interpret contig name as UTF-8")
+    }
+
+    fn pos(&self) -> u64 {
+        self.pos() as u64
+    }
+}
+
 /// Phased or unphased alleles, represented as indices.
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
 pub enum GenotypeAllele {