fix: allow leading deletions in read_pos method of CigarStringView.

src/bam/mod.rs

@@ -322,7 +322,7 @@
     /// # Arguments
     ///
     /// * `start` - Optional starting virtual offset to seek to. Throws an error if it is not
     /// a valid virtual offset.
     ///
     /// * `end` - Read until the virtual offset is less than `end`
     pub fn iter_chunk(&mut self, start: Option<i64>, end: Option<i64>) -> ChunkIterator<'_, Self> {
@@ -1122,7 +1122,7 @@
 
             //println!("{}", str::from_utf8(&header_string).unwrap());
             let rec = htslib::sam_hdr_parse(
                 ((l_text + 1) as usize).try_into().unwrap(),
                 text as *const c_char,
             );
 
@@ -1288,7 +1288,7 @@
             }
         };
 
         match self.reader.read(&mut record) {
             None => None,
             Some(Ok(_)) => Some(Ok(Rc::clone(&self.record))),
             Some(Err(err)) => Some(Err(err)),
@@ -1392,7 +1392,7 @@
                 header_string.len(),
             );
 
             let rec = htslib::sam_hdr_parse((l_text + 1), text as *const c_char);
             (*rec).text = text as *mut c_char;
             (*rec).l_text = l_text;
             rec
@@ -2997,7 +2997,7 @@
     fn test_bam_header_sync() {
         let reader = Reader::from_path("test/test_issue_156_no_text.bam").unwrap();
         let header_hashmap = Header::from_template(reader.header()).to_hashmap();
-        let header_refseqs = header_hashmap.get("SQ".into()).unwrap();
+        let header_refseqs = header_hashmap.get("SQ").unwrap();
         assert_eq!(header_refseqs[0].get("SN").unwrap(), "ref_1",);
         assert_eq!(header_refseqs[0].get("LN").unwrap(), "10000000",);
     }

src/bam/record.rs

@@ -14,7 +14,7 @@
 use std::rc::Rc;
 use std::slice;
 use std::str;
 use std::u32;
 
 use byteorder::{LittleEndian, ReadBytesExt};
 
@@ -137,7 +137,7 @@
         Record {
             inner: {
                 #[allow(clippy::uninit_assumed_init)]
                 let mut inner = unsafe { MaybeUninit::uninit().assume_init() };
                 unsafe {
                     ::libc::memcpy(
                         &mut inner as *mut htslib::bam1_t as *mut ::libc::c_void,
@@ -811,7 +811,7 @@
                     ctag,
                     b'A' as c_char,
                     size_of::<u8>() as i32,
                     [v].as_mut_ptr() as *mut u8,
                 ),
                 Aux::I8(v) => htslib::bam_aux_append(
                     self.inner_ptr_mut(),
@@ -825,7 +825,7 @@
                     ctag,
                     b'C' as c_char,
                     size_of::<u8>() as i32,
                     [v].as_mut_ptr() as *mut u8,
                 ),
                 Aux::I16(v) => htslib::bam_aux_append(
                     self.inner_ptr_mut(),
@@ -1565,7 +1565,7 @@
                     self.aux = &self.aux[offset..];
                     (tag, aux)
                 })
                 .map_err(|e| {
                     // In the case of an error, we can not safely advance in the aux data, so we terminate the Iteration
                     self.aux = &[];
                     e
@@ -1620,7 +1620,7 @@
 
     /// Return encoded base. Complexity: O(1).
     #[inline]
     pub unsafe fn encoded_base_unchecked(&self, i: usize) -> u8 {
         encoded_base_unchecked(self.encoded, i)
     }
 
@@ -1628,7 +1628,7 @@
     /// 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))
     }
 
@@ -2097,10 +2097,15 @@
                     }
                     break;
                 },
-                Cigar::Del(_) => {
-                    return Err(Error::BamUnexpectedCigarOperation {
-                        msg: "'deletion' (D) found before any operation describing read sequence".to_owned()
-                    });
+                Cigar::Del(l) => {
+                    // METHOD: leading deletions can happen in case of trimmed reads where
+                    // a primer has been removed AFTER read mapping.
+                    // Example: 24M8I8D18M9S before trimming, 32H8D18M9S after trimming
+                    // with fgbio. While leading deletions should be impossible with
+                    // normal read mapping, they make perfect sense with primer trimming
+                    // because the mapper still had the evidence to decide in favor of
+                    // the deletion via the primer sequence.
+                    rpos += l;
                 },
                 Cigar::RefSkip(_) => {
                     return Err(Error::BamUnexpectedCigarOperation {
@@ -2245,7 +2250,7 @@
     /// This function allocates memory for the state structure
     /// and initializes the iterator to the start of the modification
     /// records.
     fn new<'a>(r: &'a Record) -> Result<BaseModificationState<'a>> {
         let mut bm = unsafe {
             BaseModificationState {
                 record: r,
@@ -2269,7 +2274,7 @@
 
         let types = bm.recorded();
         bm.buffer.reserve(types.len());
         return Ok(bm);
     }
 
     pub fn buffer_next_mods(&mut self) -> Result<usize> {
@@ -2297,7 +2302,7 @@
             // not update the length so needs to be manually set
             self.buffer.set_len(ret as usize);
 
             return Ok(ret as usize);
         }
     }
 
@@ -2322,7 +2327,7 @@
     /// and the ascii code for the canonical base.
     /// If there are multiple modifications with the same code this will return the data
     /// for the first mod.  See https://github.com/samtools/htslib/issues/1635
     pub fn query_type<'a>(&self, code: i32) -> Result<BaseModificationMetadata> {
         unsafe {
             let mut strand: i32 = 0;
             let mut implicit: i32 = 0;
@@ -2337,13 +2342,13 @@
                 &mut canonical,
             );
             if ret == -1 {
                 return Err(Error::BamBaseModificationTypeNotFound);
             } else {
                 return Ok(BaseModificationMetadata {
                     strand,
                     implicit,
                     canonical: canonical.try_into().unwrap(),
                 });
             }
         }
     }
@@ -2364,7 +2369,7 @@
 }
 
 impl BaseModificationsPositionIter<'_> {
     fn new<'a>(r: &'a Record) -> Result<BaseModificationsPositionIter<'a>> {
         let state = BaseModificationState::new(r)?;
         Ok(BaseModificationsPositionIter { mod_state: state })
     }
@@ -2373,8 +2378,8 @@
         return self.mod_state.recorded();
     }
 
     pub fn query_type<'a>(&self, code: i32) -> Result<BaseModificationMetadata> {
         return self.mod_state.query_type(code);
     }
 }
 
@@ -2391,13 +2396,13 @@
         match ret {
             Ok(num_mods) => {
                 if num_mods == 0 {
                     return None;
                 } else {
                     let data = (self.mod_state.buffer_pos, self.mod_state.buffer.clone());
                     return Some(Ok(data));
                 }
             }
             Err(e) => return Some(Err(e)),
         }
     }
 }
@@ -2410,7 +2415,7 @@
 }
 
 impl BaseModificationsIter<'_> {
     fn new<'a>(r: &'a Record) -> Result<BaseModificationsIter<'a>> {
         let state = BaseModificationState::new(r)?;
         Ok(BaseModificationsIter {
             mod_state: state,
@@ -2422,8 +2427,8 @@
         return self.mod_state.recorded();
     }
 
     pub fn query_type<'a>(&self, code: i32) -> Result<BaseModificationMetadata> {
         return self.mod_state.query_type(code);
     }
 }
 
@@ -2457,7 +2462,7 @@
             self.mod_state.buffer[self.buffer_idx],
         );
         self.buffer_idx += 1;
         return Some(Ok(data));
     }
 }