qnames_at_base

src/lib.rs

@@ -403,6 +403,7 @@ pub fn scan_sa(
     }
     Ok((results_start, results_end))
 }
+
 pub fn records_at_base(
     bam: &mut rust_htslib::bam::IndexedReader,
     chr: &str,
@@ -443,6 +444,48 @@ pub fn records_at_base(
     Ok(bases)
 }
 
+pub fn qnames_at_base(
+    bam: &mut rust_htslib::bam::IndexedReader,
+    chr: &str,
+    start: i32,
+    with_next_ins: bool,
+) -> Result<Vec<(String, u8)>> {
+    let stop = start + 1;
+    let mut bases = Vec::new();
+    bam.fetch((chr, start, stop))?;
+    let mut bam_pileup = Vec::new();
+    for p in bam.pileup() {
+        let pileup = p.context(format!(
+            "Can't pilup bam at position {}:{}-{}",
+            chr, start, stop
+        ))?;
+        let position = pileup.pos() as i32;
+        if position == start {
+            for alignment in pileup.alignments() {
+                match alignment.indel() {
+                    bam::pileup::Indel::Ins(_len) => bam_pileup.push(b'I'),
+                    bam::pileup::Indel::Del(_len) => bam_pileup.push(b'D'),
+                    _ => {
+                        let record = alignment.record();
+                        let qname = String::from_utf8(record.qname().to_vec())?;
+                        if record.seq_len() > 0 {
+                            if let Some(b) = hts_base_at(&record, start as u32, with_next_ins)? {
+                                bases.push((qname, b));
+                            }
+                        } else {
+                            if alignment.is_del() {
+                                bases.push((qname, b'D'));
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+    Ok(bases)
+}
+
+
 #[cfg(test)]
 mod tests {
     // Note this useful idiom: importing names from outer (for mod tests) scope.