monoall correction

src/rna_seq.rs

@@ -1,16 +1,9 @@
-use std::{
-    collections::{BinaryHeap, HashMap},
-    fs::File,
-    io::BufReader,
-};
+use std::{collections::HashMap, fs::File, io::BufReader};
 
 use log::info;
 use noodles_gff as gff;
 
-use crate::{
-    utils::chi_square_test_for_proportions,
-    variants::{ReferenceAlternative, Variants},
-};
+use crate::variants::{ReferenceAlternative, Variants};
 
 #[derive(Debug, Default, Clone)]
 struct Gene {
@@ -29,16 +22,17 @@ pub fn find_monoallelics(const_variants: &Variants, rna_bam: &str) -> anyhow::Re
 
     let mut genes_by_contig: HashMap<String, Vec<Gene>> = HashMap::new();
     let mut curr_gene = Gene::default();
+    let mut contig = String::new();
 
     for result in reader.records() {
         let record = result?;
+        contig = record.reference_sequence_name().to_string();
 
         if record.ty() == "exon" {
             if let Some(gene) = record.attributes().get("gene") {
                 let gene_name = gene.to_string();
-                let contig = record.reference_sequence_name().to_string();
                 if curr_gene.name != gene_name {
-                    genes_by_contig.entry(contig).or_default().push(curr_gene);
+                    genes_by_contig.entry(contig.clone()).or_default().push(curr_gene);
                     curr_gene = Gene::default();
                     curr_gene.name = gene_name;
                 }
@@ -50,6 +44,8 @@ pub fn find_monoallelics(const_variants: &Variants, rna_bam: &str) -> anyhow::Re
         }
     }
 
+    genes_by_contig.entry(contig).or_default().push(curr_gene);
+
     let mut curr_contig = String::default();
     let mut genes_iter = std::slice::Iter::default();
     let mut results = Vec::new();
@@ -97,11 +93,8 @@ 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 (
             ReferenceAlternative::Nucleotide(ref_base),
             ReferenceAlternative::Nucleotide(alt_base),
@@ -131,11 +124,6 @@ pub fn find_monoallelics(const_variants: &Variants, rna_bam: &str) -> anyhow::Re
                 variant.contig,
                 variant.position,
             ));
-
-            // println!(
-            //     "{gene_name} exon_{exon}: {}:{}",
-            //     variant.contig, variant.position
-            // );
         }
     }