Thomas 1 year ago
parent
commit
d477b4ed8d
1 changed files with 31 additions and 7 deletions
  1. 31 7
      src/rna_seq.rs

+ 31 - 7
src/rna_seq.rs

@@ -97,6 +97,9 @@ pub fn find_monoallelics(const_variants: &Variants, rna_bam: &str) -> anyhow::Re
     }
 
     let mut bam = rust_htslib::bam::IndexedReader::from_path(rna_bam)?;
+    // let mut curr_gene = String::default();
+    // let mut curr_gene_checked = 0;
+    let mut monoall: HashMap<String, Vec<(bool, String, String, u32)>> = HashMap::new();
     for (variant, gene_name, exon) in results {
         // let mut variant = variant;
         if let (
@@ -114,21 +117,42 @@ pub fn find_monoallelics(const_variants: &Variants, rna_bam: &str) -> anyhow::Re
             if rna_bases.len() < min_rna_depth {
                 continue;
             }
-            
+
             let ref_base = ref_base.clone().into_u8();
             let alt_base = alt_base.clone().into_u8();
 
             let n_rna_ref = rna_bases.iter().filter(|(_, b)| *b == ref_base).count();
             let n_rna_alt = rna_bases.iter().filter(|(_, b)| *b == alt_base).count();
 
-            if n_rna_ref == 0 || n_rna_alt == 0 {
-                println!(
-                    "{gene_name} exon_{exon}: {}:{}",
-                    variant.contig, variant.position
-                );
-            }
+            let is_monoal = n_rna_ref == 0 || n_rna_alt == 0;
+            monoall.entry(gene_name).or_default().push((
+                is_monoal,
+                exon,
+                variant.contig,
+                variant.position,
+            ));
+
+            // println!(
+            //     "{gene_name} exon_{exon}: {}:{}",
+            //     variant.contig, variant.position
+            // );
         }
     }
 
+    monoall
+        .iter()
+        .filter(|(_k, v)| v.iter().any(|v| v.0))
+        .for_each(|(k, v)| {
+            let n = v.iter().filter(|v| v.0).count();
+            println!(
+                "{k}\t{n}/{}\t{}",
+                v.len(),
+                v.iter()
+                    .map(|(_, _, c, p)| format!("{c}:{p}"))
+                    .collect::<Vec<String>>()
+                    .join("|")
+            )
+        });
+
     Ok(())
 }