records_at_base

src/lib.rs

@@ -4,7 +4,6 @@ use rust_htslib::{
     bam::{ext::BamRecordExtensions, record, Read, Record},
 };
 
-use average::Mean;
 pub fn get_hts_nt_pileup(
     bam: &mut rust_htslib::bam::IndexedReader,
     chr: &str,
@@ -366,13 +365,91 @@ pub fn swap_by_primary(bam_path: &str, records: Vec<Record>) -> Vec<Record> {
     res_records
 }
 
+pub fn scan_sa(
+    bam: &mut rust_htslib::bam::IndexedReader,
+    chr: &str,
+    start: i32,
+    stop: i32,
+    mapq: u8,
+) -> Result<(Vec<Vec<Record>>, Vec<Vec<Record>>)> {
+    bam.fetch((chr, start - 1, stop))?;
+    let length = stop - start;
+    let mut results_start: Vec<Vec<Record>> = Vec::new();
+    results_start.resize(length as usize, vec![]);
+    let mut results_end: Vec<Vec<Record>> = Vec::new();
+    results_end.resize(length as usize, vec![]);
+
+    for read in bam.records() {
+        let record = read.context(format!("Error while parsing record"))?;
+        let rstart = record.pos() as i32;
+        let rend = record.reference_end() as i32;
+
+        if rstart >= start && rstart < stop {
+            if record.mapq() >= mapq && record.aux(b"SA").is_ok() {
+                
+                let index = rstart - start;
+                let u = results_start.get_mut(index as usize).unwrap();
+                u.push(record.clone());
+            }
+        }
+
+        if rend >= start && rend < stop && record.aux(b"SA").is_ok() {
+            if record.mapq() >= mapq {
+                let index = rend - start;
+                let u = results_end.get_mut(index as usize).unwrap();
+                u.push(record.clone());
+            }
+        }
+    }
+    Ok((results_start, results_end))
+}
+pub fn records_at_base(
+    bam: &mut rust_htslib::bam::IndexedReader,
+    chr: &str,
+    start: i32,
+    with_next_ins: bool,
+) -> Result<Vec<(Record, 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();
+                        if record.seq_len() > 0 {
+                            if let Some(b) = hts_base_at(&record, start as u32, with_next_ins)? {
+                                bases.push((record, b));
+                            }
+                        } else {
+                            if alignment.is_del() {
+                                bases.push((record, b'D'));
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+    Ok(bases)
+}
+
 #[cfg(test)]
 mod tests {
     // Note this useful idiom: importing names from outer (for mod tests) scope.
     use super::*;
 
     #[test]
-    fn test() {
+    fn test_se() {
         let bam_path = "/data/longreads_basic_pipe/CAMARA/diag/CAMARA_diag_hs1.bam";
         let chr = "chr9";
         let start = 21839796;
@@ -384,8 +461,30 @@ mod tests {
         // println!("{a:?}");
 
         let res = get_start_end_qual_rec(&mut bam, chr, start, start + 1, mapq).unwrap();
-        let (s, e) = res;
+        let (_, e) = res;
         let res_records = swap_by_primary(bam_path, e.get(0).unwrap().clone());
         println!("{res_records:?}");
+        println!("{}", res_records.len());
+        assert!(res_records.len() == 12);
+    }
+
+    #[test]
+    fn test_sa() {
+        let bam_path = "/data/longreads_basic_pipe/CAMARA/diag/CAMARA_diag_hs1.bam";
+        let chr = "chr9";
+        let start = 21839796;
+        let mapq = 50;
+
+        let mut bam = rust_htslib::bam::IndexedReader::from_path(bam_path).unwrap();
+
+        let a = get_start_end_qual(&mut bam, chr, start, start + 1, mapq).unwrap().pop().unwrap();
+
+        let (_, e) = get_start_end_qual_rec(&mut bam, chr, start, start + 1, mapq).unwrap();
+        let res_records = swap_by_primary(bam_path, e.get(0).unwrap().clone());
+
+        let (_, e) = scan_sa(&mut bam, chr, start, start + 1, mapq).unwrap();
+        let res_records_sa = swap_by_primary(bam_path, e.get(0).unwrap().clone());
+        assert_eq!(a as usize, res_records.len());
+        assert_eq!(a as usize, res_records_sa.len());
     }
 }