test ins

Cargo.toml

@@ -14,5 +14,3 @@ log = "^0.4.22"
 env_logger = "^0.11.3"
 uuid = { version = "1.10.0", features = ["v4"] }
 noodles-fasta = "0.41.0"
-
-

src/lib.rs

@@ -239,7 +239,7 @@ pub fn primary_record(bam_path: &str, record: &Record) -> Option<Record> {
                             let position: i32 = s.get(1).unwrap().parse().unwrap();
                             (chr, position)
                         })
-                        .collect();
+                       .collect();
                 }
                 _ => (),
             },
@@ -711,10 +711,37 @@ pub fn write_fasta(path: &str, records: &Vec<Record>) -> anyhow::Result<()> {
     Ok(())
 }
 
+pub fn calculate_end_position(start: u32, cigar_str: &str) -> anyhow::Result<u32> {
+    let mut end = start;
+    let mut count = String::new();
+
+    for c in cigar_str.chars() {
+        if c.is_ascii_digit() {
+            count.push(c);
+        } else {
+            let len: u32 = count.parse().map_err(|e| anyhow!(format!("Parse error: {}", e)))?;
+            count.clear();
+            match c {
+                'M' | 'D' | 'N' | '=' | 'X' => {
+                    end += len;
+                }
+                'I' | 'S' | 'H' | 'P' => {
+                    // These operations do not consume reference positions
+                }
+                _ => return Err(anyhow!(format!("Invalid CIGAR operation: {}", c))),
+            }
+        }
+    }
+
+    Ok(end)
+}
+
+
 #[cfg(test)]
 mod tests {
     use env_logger::Env;
     use log::info;
+    use rust_htslib::bam::IndexedReader;
     use uuid::Uuid;
 
     // Note this useful idiom: importing names from outer (for mod tests) scope.
@@ -833,4 +860,19 @@ mod tests {
         info!("{} group of records", records.len());
         Ok(())
     }
+
+    #[test]
+    fn ins() {
+        init();
+        let id = "SALICETTO";
+        let result_dir = "/data/longreads_basic_pipe";
+        let diag_bam_path = format!("{result_dir}/{id}/diag/{id}_diag_hs1.bam");
+        let ins_pos = ("chr7", 511_746);
+        let mut bam = IndexedReader::from_path(diag_bam_path).unwrap();
+
+        let reads = qnames_at_base(&mut bam, ins_pos.0, ins_pos.1, true).unwrap();
+
+        let n_ins = reads.iter().filter(|(_, b)| *b == b'I').count();
+        assert_eq!(n_ins, 5);
+    }
 }