Reference sequence reconstruction method on Record

src/bam/md_align.rs

@@ -1207,16 +1207,6 @@ mod tests {
             assert_eq!(rec.aux_md().unwrap(), TEST_MD[i]);
             assert_eq!(rec.cigar().to_string(), TEST_CIGAR[i]);
 
-            let ap_iter = AlignPosIter::new_from_record(&rec)
-                .ok()
-                .expect("Creating AlignPosIter");
-            let ap_res: Result<Vec<AlignPos>, MDAlignError> = ap_iter.collect();
-            let aps = ap_res.ok().expect("Unable to create Vec<AlignPos>");
-            let ap_read_line: String = aps.iter().map(|ap| ap.read_line_char()).collect();
-            assert_eq!(ap_read_line, TEST_READ_LINE[i]);
-            let ap_ref_line: String = aps.iter().map(|ap| ap.ref_line_char()).collect();
-            assert_eq!(ap_ref_line, TEST_REF_LINE[i]);
-
             let rap_iter = MDAlignPosIter::new(&rec)
                 .ok()
                 .expect("Creating MDAlignPosIter");

src/bam/record.rs

@@ -16,6 +16,7 @@ use itertools::Itertools;
 use regex::Regex;
 
 use bam::{AuxWriteError, HeaderView, ReadError};
+use bam::md_align::{MDAlignError, MDAlignPos, MDAlignPosIter};
 use htslib;
 use utils;
 
@@ -600,6 +601,20 @@ impl Record {
             _ => None,
         })
     }
+
+    /// Reconstruct reference sequence from Cigar and MD information
+    /// along with the read sequence.
+    ///
+    /// # Errors
+    ///
+    /// An error variant is returned when the record has no MD aux
+    /// field, when the MD aux field is malformed, or when there are
+    /// inconsistencies between the Cigar string, the MD field, and
+    /// the read sequence.
+    pub fn reference_md_seq(&self) -> Result<Vec<u8>, MDAlignError> {
+        let md_align: Result<Vec<MDAlignPos>, MDAlignError> = MDAlignPosIter::new(&self)?.collect();
+        Ok( md_align?.iter().filter_map(|pos| pos.ref_nt()).collect() )
+    }
 }
 
 impl Drop for Record {