浏览代码

bam_compo

Thomas 1 年之前
父节点
当前提交
f1fca549e7
共有 1 个文件被更改,包括 25 次插入5 次删除
  1. 25 5
      src/lib.rs

+ 25 - 5
src/lib.rs

@@ -35,10 +35,8 @@ pub fn get_hts_nt_pileup(
                             if let Some(b) = hts_base_at(&record, start as u32, with_next_ins)? {
                                 bases.push(b);
                             }
-                        } else {
-                            if alignment.is_del() {
-                                bases.push(b'D');
-                            }
+                        } else if alignment.is_del() {
+                            bases.push(b'D');
                         }
                     }
                 }
@@ -703,6 +701,22 @@ pub fn sa_ratio(
         .collect()
 }
 
+pub fn bam_compo(file_path: &str, sample_size: usize) -> Result<Vec<(String, f64 )>> {
+    let mut bam = bam::Reader::from_path(file_path)?;
+
+    let mut rgs: HashMap<String, u64> = HashMap::new();
+    for result in bam.records().filter_map(Result::ok).take(sample_size)
+    {
+        if let std::result::Result::Ok(t) = result.aux(b"RG") {
+            if let record::Aux::String(s) = t {
+                *rgs.entry(s.to_string()).or_default() += 1;
+            }
+        }
+    }
+
+    Ok(rgs.into_iter().map(|(rg, n)| (rg.to_string(), n as f64 * 100.0 / sample_size as f64)).collect())
+}
+
 #[cfg(test)]
 mod tests {
     use env_logger::Env;
@@ -771,8 +785,14 @@ mod tests {
             start,
             stop,
             mapq,
-            1_000
+            1_000,
         );
         println!("{r:?}");
     }
+
+    #[test]
+    fn rg() {
+        let rg = bam_compo("/data/longreads_basic_pipe/KENNOUCHE/mrd/KENNOUCHE_mrd_hs1.bam", 20000);
+        println!("{rg:#?}");
+    }
 }