counts at with tabix on counts files

src/helpers.rs

@@ -381,6 +381,21 @@ where
     sum.into() / count as f64
 }
 
+pub fn std_dev<T>(values: &[T]) -> f64
+where
+    T: Copy + Into<f64>,
+{
+    let n = values.len();
+    if n == 0 {
+        return 0.0;
+    }
+
+    let floats: Vec<f64> = values.iter().map(|&x| x.into()).collect();
+    let mean = floats.iter().sum::<f64>() / n as f64;
+    let variance = floats.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / n as f64;
+    variance.sqrt()
+}
+
 pub fn bin_data<V>(data: Vec<V>, bin_size: V) -> Vec<(V, usize)>
 where
     V: Copy + Default + PartialOrd + std::ops::AddAssign + std::ops::Add<Output = V>,

src/scan/scan.rs

@@ -1166,26 +1166,45 @@ impl Label for SomaticScan {
 #[cfg(test)]
 mod tests {
     use crate::{
-        helpers::test_init, io::readers::TabixReader, positions::GenomeRange,
-        scan::scan::CountsReader,
+        config::Config,
+        helpers::{mean, std_dev, test_init},
+        io::readers::TabixReader,
+        pipes::Initialize,
+        positions::GenomeRange,
+        runners::Run,
+        scan::scan::{CountsReader, SomaticScan},
     };
 
     #[test]
     fn tabix_counts_at() -> anyhow::Result<()> {
         test_init();
 
-        let contig = "chr16";
-        let start0 = 146763;
-        let end0 = 160013;
-        let case = "DUMCO";
-        let reader = TabixReader::open(
-            &format!("/mnt/beegfs02/scratch/t_steimle/data/wgs/{case}/diag/counts/{contig}_count.tsv.gz"),
-        )?;
+        let contig = "chr17";
+        let start0 = 32_040_661;
+        let end0 = 32_323_039;
+        let start0 = 32_177_773 - 10;
+        let end0 = start0 + 1000;
+        let case = "CHALO";
+        let time = "diag";
+
+        let config = Config::default();
+        SomaticScan::initialize(case, &config)?.run()?;
+
+        let reader = TabixReader::open(&format!(
+            "/mnt/beegfs02/scratch/t_steimle/data/wgs/{case}/{time}/counts/{contig}_count.tsv.gz"
+        ))?;
+        // TabixReader::open_smb
         let mut creader = CountsReader::from_tabix_reader(reader);
 
         let region = GenomeRange::from_1_inclusive(contig, start0, end0);
         let u = creader.depths(&region)?;
+        let cov_mean = mean(&u);
+        let cov_min = u.iter().min().cloned().unwrap_or(0);
+        let cov_max = u.iter().max().cloned().unwrap_or(u32::MAX);
+        let cov_sd = std_dev(&u);
+
         println!("{u:?}");
+        println!("mean={cov_mean:.2} min={cov_min} max={cov_max} std_dev={cov_sd:.2}");
 
         Ok(())
     }