depths

res_150.tsv

@@ -1,276 +0,0 @@
-chr19:0-1000000	242807
-chr19:1000000-2000000	30264
-chr19:2000000-3000000	30893
-chr19:3000000-4000000	29720
-chr19:4000000-5000000	33420
-chr19:5000000-6000000	24684
-chr19:6000000-7000000	26470
-chr19:7000000-8000000	41607
-chr19:8000000-9000000	190288
-chr19:9000000-10000000	28049
-chr19:10000000-11000000	22145
-chr19:11000000-12000000	30992
-chr19:12000000-13000000	54074
-chr19:13000000-14000000	21388
-chr19:14000000-15000000	28306
-chr19:15000000-16000000	33537
-chr19:16000000-17000000	22089
-chr19:17000000-18000000	22006
-chr19:18000000-19000000	19936
-chr19:19000000-20000000	21959
-chr19:20000000-21000000	54625
-chr19:21000000-22000000	31606
-chr19:22000000-23000000	53769
-chr19:23000000-24000000	50640
-chr19:24000000-25000000	451093
-chr19:25000000-26000000	821050
-chr19:26000000-27000000	996399
-chr19:27000000-28000000	998739
-chr19:28000000-29000000	992286
-chr19:29000000-30000000	771181
-chr19:30000000-31000000	53419
-chr19:31000000-32000000	30759
-chr19:32000000-33000000	17330
-chr19:33000000-34000000	6470
-chr19:34000000-35000000	7465
-chr19:35000000-36000000	32951
-chr19:36000000-37000000	15000
-chr19:37000000-38000000	35233
-chr19:38000000-39000000	189271
-chr19:39000000-40000000	93191
-chr19:40000000-41000000	136843
-chr19:41000000-42000000	26123
-chr19:42000000-43000000	71584
-chr19:43000000-44000000	49779
-chr19:44000000-45000000	24306
-chr19:45000000-46000000	157241
-chr19:46000000-47000000	57600
-chr19:47000000-48000000	34246
-chr19:48000000-49000000	24620
-chr19:49000000-50000000	26481
-chr19:50000000-51000000	239162
-chr19:51000000-52000000	27511
-chr19:52000000-53000000	36043
-chr19:53000000-54000000	159753
-chr19:54000000-55000000	21201
-chr19:55000000-56000000	35163
-chr19:56000000-57000000	52336
-chr19:57000000-58000000	94880
-chr19:58000000-59000000	49762
-chr19:59000000-60000000	35744
-chr19:60000000-61000000	31797
-chr19:61000000-62000000	46945
-chr8:0-1000000	68533
-chr8:1000000-2000000	49291
-chr8:2000000-3000000	27463
-chr8:3000000-4000000	23922
-chr8:4000000-5000000	14892
-chr8:5000000-6000000	20900
-chr8:6000000-7000000	123619
-chr8:7000000-8000000	241166
-chr8:8000000-9000000	17935
-chr8:9000000-10000000	14873
-chr8:10000000-11000000	28569
-chr8:11000000-12000000	472830
-chr8:12000000-13000000	816157
-chr8:13000000-14000000	18511
-chr8:14000000-15000000	20348
-chr8:15000000-16000000	26961
-chr8:16000000-17000000	25887
-chr8:17000000-18000000	19004
-chr8:18000000-19000000	13391
-chr8:19000000-20000000	11021
-chr8:20000000-21000000	19115
-chr8:21000000-22000000	13760
-chr8:22000000-23000000	15108
-chr8:23000000-24000000	25297
-chr8:24000000-25000000	20387
-chr8:25000000-26000000	30644
-chr8:26000000-27000000	20773
-chr8:27000000-28000000	26018
-chr8:28000000-29000000	12223
-chr8:29000000-30000000	27797
-chr8:30000000-31000000	15637
-chr8:31000000-32000000	28684
-chr8:32000000-33000000	19944
-chr8:33000000-34000000	26331
-chr8:34000000-35000000	30672
-chr8:35000000-36000000	27036
-chr8:36000000-37000000	33291
-chr8:37000000-38000000	13190
-chr8:38000000-39000000	18615
-chr8:39000000-40000000	20351
-chr8:40000000-41000000	33001
-chr8:41000000-42000000	21425
-chr8:42000000-43000000	19358
-chr8:43000000-44000000	69550
-chr8:44000000-45000000	782708
-chr8:45000000-46000000	999539
-chr8:46000000-47000000	348700
-chr8:47000000-48000000	81035
-chr8:48000000-49000000	32569
-chr8:49000000-50000000	30055
-chr8:50000000-51000000	30421
-chr8:51000000-52000000	32614
-chr8:52000000-53000000	26036
-chr8:53000000-54000000	29016
-chr8:54000000-55000000	34791
-chr8:55000000-56000000	29380
-chr8:56000000-57000000	39563
-chr8:57000000-58000000	77406
-chr8:58000000-59000000	21889
-chr8:59000000-60000000	38213
-chr8:60000000-61000000	34406
-chr8:61000000-62000000	31596
-chr8:62000000-63000000	27351
-chr8:63000000-64000000	31328
-chr8:64000000-65000000	38265
-chr8:65000000-66000000	31987
-chr8:66000000-67000000	22378
-chr8:67000000-68000000	30152
-chr8:68000000-69000000	38946
-chr8:69000000-70000000	41975
-chr8:70000000-71000000	26788
-chr8:71000000-72000000	49096
-chr8:72000000-73000000	40323
-chr8:73000000-74000000	24162
-chr8:74000000-75000000	32208
-chr8:75000000-76000000	45718
-chr8:76000000-77000000	30392
-chr8:77000000-78000000	25887
-chr8:78000000-79000000	25293
-chr8:79000000-80000000	28518
-chr8:80000000-81000000	30245
-chr8:81000000-82000000	30004
-chr8:82000000-83000000	27401
-chr8:83000000-84000000	28666
-chr8:84000000-85000000	43135
-chr8:85000000-86000000	38383
-chr8:86000000-87000000	863674
-chr8:87000000-88000000	21664
-chr8:88000000-89000000	53232
-chr8:89000000-90000000	40596
-chr8:90000000-91000000	57164
-chr8:91000000-92000000	29208
-chr8:92000000-93000000	50302
-chr8:93000000-94000000	31246
-chr8:94000000-95000000	25028
-chr8:95000000-96000000	27114
-chr8:96000000-97000000	37175
-chr8:97000000-98000000	28165
-chr8:98000000-99000000	41410
-chr8:99000000-100000000	48531
-chr8:100000000-101000000	46196
-chr8:101000000-102000000	30069
-chr8:102000000-103000000	26116
-chr8:103000000-104000000	29067
-chr8:104000000-105000000	28685
-chr8:105000000-106000000	26942
-chr8:106000000-107000000	21335
-chr8:107000000-108000000	22272
-chr8:108000000-109000000	23850
-chr8:109000000-110000000	37156
-chr8:110000000-111000000	32476
-chr8:111000000-112000000	17296
-chr8:112000000-113000000	23359
-chr8:113000000-114000000	22644
-chr8:114000000-115000000	25971
-chr8:115000000-116000000	57131
-chr8:116000000-117000000	22939
-chr8:117000000-118000000	18094
-chr8:118000000-119000000	22998
-chr8:119000000-120000000	25292
-chr8:120000000-121000000	45518
-chr8:121000000-122000000	26975
-chr8:122000000-123000000	38132
-chr8:123000000-124000000	26444
-chr8:124000000-125000000	19303
-chr8:125000000-126000000	13024
-chr8:126000000-127000000	18447
-chr8:127000000-128000000	31443
-chr8:128000000-129000000	35427
-chr8:129000000-130000000	37307
-chr8:130000000-131000000	29259
-chr8:131000000-132000000	12967
-chr8:132000000-133000000	43873
-chr8:133000000-134000000	22518
-chr8:134000000-135000000	12325
-chr8:135000000-136000000	19013
-chr8:136000000-137000000	29663
-chr8:137000000-138000000	28612
-chr8:138000000-139000000	32380
-chr8:139000000-140000000	34589
-chr8:140000000-141000000	11857
-chr8:141000000-142000000	22650
-chr8:142000000-143000000	19581
-chr8:143000000-144000000	19862
-chr8:144000000-145000000	37079
-chr8:145000000-146000000	46578
-chr8:146000000-147000000	18564
-chr15:0-1000000	963519
-chr15:1000000-2000000	951973
-chr15:2000000-3000000	957629
-chr15:3000000-4000000	1000000
-chr15:4000000-5000000	999069
-chr15:5000000-6000000	964149
-chr15:6000000-7000000	978855
-chr15:7000000-8000000	999218
-chr15:8000000-9000000	997631
-chr15:9000000-10000000	996622
-chr15:10000000-11000000	999366
-chr15:11000000-12000000	999668
-chr15:12000000-13000000	999696
-chr15:13000000-14000000	932346
-chr15:14000000-15000000	474988
-chr15:15000000-16000000	440571
-chr15:16000000-17000000	941262
-chr15:17000000-18000000	772487
-chr15:18000000-19000000	252891
-chr15:19000000-20000000	210773
-chr15:20000000-21000000	463315
-chr15:21000000-22000000	165392
-chr15:22000000-23000000	104408
-chr15:23000000-24000000	18671
-chr15:24000000-25000000	18380
-chr15:25000000-26000000	34570
-chr15:26000000-27000000	476150
-chr15:27000000-28000000	99471
-chr15:28000000-29000000	208830
-chr15:29000000-30000000	33212
-chr15:30000000-31000000	141351
-chr15:31000000-32000000	22255
-chr15:32000000-33000000	109464
-chr15:33000000-34000000	24037
-chr15:34000000-35000000	22320
-chr15:35000000-36000000	47449
-chr15:36000000-37000000	25336
-chr15:37000000-38000000	28685
-chr15:38000000-39000000	21243
-chr15:39000000-40000000	17935
-chr15:40000000-41000000	28988
-chr15:41000000-42000000	152518
-chr15:42000000-43000000	43394
-chr15:43000000-44000000	27120
-chr15:44000000-45000000	43748
-chr15:45000000-46000000	10121
-chr15:46000000-47000000	10470
-chr15:47000000-48000000	48879
-chr15:48000000-49000000	30549
-chr15:49000000-50000000	29068
-chr15:50000000-51000000	31980
-chr15:51000000-52000000	30361
-chr15:52000000-53000000	35331
-chr15:53000000-54000000	39643
-chr15:54000000-55000000	30611
-chr15:55000000-56000000	17320
-chr15:56000000-57000000	10778
-chr15:57000000-58000000	15933
-chr15:58000000-59000000	16828
-chr15:59000000-60000000	28083
-chr15:60000000-61000000	16235
-chr15:61000000-62000000	23519
-chr15:62000000-63000000	27947
-chr15:63000000-64000000	16295
-chr15:64000000-65000000	15349
-chr15:65000000-66000000	18196
-chr15:66000000-

src/collection/bam_stats.rs

@@ -1,13 +1,15 @@
 use std::{
     collections::BTreeMap,
-    fmt, fs,
+    fmt,
+    fs::{self, File},
     hash::{Hash, Hasher},
-    io::{Read, Write},
+    io::{BufWriter, Read, Write},
     path::{Path, PathBuf},
 };
 
 use anyhow::Context;
 use log::{debug, info};
+use rayon::prelude::*;
 use rust_htslib::{
     bam::{ext::BamRecordExtensions, record::Aux, Read as BamRead},
     htslib::{BAM_FDUP, BAM_FQCFAIL, BAM_FSECONDARY, BAM_FSUPPLEMENTARY, BAM_FUNMAP},
@@ -15,7 +17,13 @@ use rust_htslib::{
 use rustc_hash::{FxHashMap, FxHashSet, FxHasher};
 use serde::{Deserialize, Serialize};
 
-use crate::{config::Config, helpers::get_genome_sizes};
+use crate::{
+    config::Config,
+    helpers::get_genome_sizes,
+    io::{dict::read_dict, readers::TabixReader},
+    positions::GenomeRange,
+    uniqueness::{GenomeIndex, GenomicPos, UniqueContext},
+};
 
 /// Flags to skip: unmapped, secondary, QC fail, supplementary
 const SKIP_FLAGS: u16 = (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FSUPPLEMENTARY) as u16;
@@ -1348,10 +1356,262 @@ impl QNameSet {
     }
 }
 
+/// Computes a per-window SV detection quality index for a single sequenced sample.
+///
+/// For each 1 Mb genomic window, counts the number of positions at risk of
+/// producing false negatives in structural variant calling, defined as positions
+/// where **at least one** of the following conditions holds:
+///
+/// - **Low coverage**: per-base depth < `min_coverage` (default threshold: 15x,
+///   below which most long-read SV callers show increased false negative rates)
+/// - **Multimapping risk**: the minimal unique sequence length on either flank
+///   exceeds `half_median` (= `median_read_length / 2`), meaning a read of
+///   median length cannot unambiguously anchor to this position
+///
+/// The two conditions are combined as a **union** (logical OR) at the position
+/// level, avoiding double-counting while capturing both independent sources of
+/// false negatives.
+///
+/// # Multimapping criterion
+///
+/// Multimapping risk is assessed using a precomputed genome-wide uniqueness
+/// index ([`GenomeIndex`]) based on suffix arrays. For each position `pos`,
+/// the index provides:
+///
+/// - `before`: minimal unique sequence length ending at `pos`  
+///   (`genome[pos - k .. pos)` is unique)
+/// - `after`: minimal unique sequence length starting at `pos`  
+///   (`genome[pos .. pos + k)` is unique)
+///
+/// A position is considered multimappable if either flank requires more
+/// sequence than the read's half-median to achieve uniqueness:
+///
+/// ```text
+/// is_multimap = (before > half_median) || (after > half_median)
+/// ```
+///
+/// where `half_median = median_read_length / 2`. This is case-specific:
+/// a sample with longer reads will have fewer multimappable positions.
+///
+/// # Read length
+///
+/// The median read length is extracted from the cached [`WGSBamStats`] for
+/// the given `(id, time)` combination. It is computed on primary, mapped,
+/// non-duplicate reads passing the configured MAPQ threshold — excluding
+/// unmapped and supplementary alignments that would otherwise deflate the
+/// median.
+///
+/// # Output
+///
+/// Writes a TSV file with two columns (tab-separated, with header):
+///
+/// ```text
+/// region              n_at_risk
+/// chr1:0-1000000      2880
+/// chr1:1000000-2000000    0
+/// ...
+/// ```
+///
+/// - `region`: genomic window in `contig:start-end` format (0-based, half-open)
+/// - `n_at_risk`: number of positions at risk in the window (0 – window size)
+///
+/// Windows at contig boundaries may be shorter than 1 Mb.
+/// `chrM` is excluded.
+///
+/// # Arguments
+///
+/// * `id` — sample identifier (e.g. `"DUMCO"`)
+/// * `time` — time point (e.g. `"diag"`, `"norm"`)
+/// * `config` — pipeline configuration (BAM paths, count file paths, dict file)
+/// * `genome_index` — precomputed genome-wide uniqueness index (hs1/T2T-CHM13v2)
+/// * `min_coverage` — minimum per-base depth threshold (positions below are at risk)
+/// * `output_path` — path to the output TSV file (created or overwritten)
+///
+/// # Errors
+///
+/// Returns an error if:
+/// - [`WGSBamStats`] cannot be loaded for the given `(id, time)`
+/// - The dict file cannot be parsed
+/// - Any per-contig count file (`{contig}_count.tsv.gz`) cannot be opened
+/// - A count file line has an unexpected number of fields or unparseable values
+/// - The output file cannot be created or written
+pub fn sv_quality_index(
+    id: &str,
+    time: &str,
+    config: &Config,
+    genome_index: &GenomeIndex,
+    min_coverage: u32,
+    output_path: impl AsRef<Path>,
+) -> anyhow::Result<()> {
+    let bam_stats = WGSBamStats::open(id, time, config)?;
+    let half_median = bam_stats.median_read_length as u32 / 2;
+
+    info!(
+        "sv_quality_index: id={} time={} median_read_length={}bp half_median={}bp min_coverage={}x",
+        id, time, bam_stats.median_read_length, half_median, min_coverage
+    );
+
+    let contigs: Vec<(String, u32)> = read_dict(&config.dict_file)?
+        .into_iter()
+        .filter(|(ctg, _)| ctg != "chrM")
+        .collect();
+
+    info!("Processing {} contigs", contigs.len());
+
+    let interval = 1_000_000u32;
+
+    let total_windows: u32 = contigs
+        .iter()
+        .map(|(_, size)| size.div_ceil(interval))
+        .sum();
+
+    // Compteur atomique pour la progression globale
+    let window_counter = std::sync::atomic::AtomicU32::new(0);
+
+    // Traitement parallèle par contig
+    let results: Vec<(String, Vec<(u32, u32, u32)>)> = contigs
+        .par_iter()
+        .enumerate()
+        .map(|(contig_idx, (contig, size))| {
+            let depth_path = format!("{}/{}_count.tsv.gz", config.dir_count(id, time), contig);
+
+            info!(
+                "[{}/{}] Opening: {}",
+                contig_idx + 1,
+                contigs.len(),
+                depth_path
+            );
+
+            let mut tabix = TabixReader::open(&depth_path)
+                .with_context(|| format!("Failed to open tabix: {}", depth_path))?;
+
+            let mut windows: Vec<(u32, u32, u32)> = Vec::new(); // (win_start, win_end, n_at_risk)
+            let mut current_interval = 0u32;
+            let mut contig_at_risk = 0u32;
+            let mut contig_positions = 0u32;
+
+            loop {
+                let win_start = current_interval;
+                let win_end = (current_interval + interval).min(*size);
+
+                let region = GenomeRange::new(contig, win_start, win_end);
+                let mut n_at_risk = 0u32;
+                let mut n_positions = 0u32;
+
+                tabix.fetch_with(&region, |line| {
+                    let fields: Vec<&str> = line.split('\t').collect();
+                    anyhow::ensure!(
+                        fields.len() == 12,
+                        "Expected 12 fields, got {} | {:?}",
+                        fields.len(),
+                        &fields[..fields.len().min(4)]
+                    );
+
+                    let row_start: u32 = fields[1]
+                        .parse()
+                        .with_context(|| format!("Invalid start: '{}'", fields[1]))?;
+
+                    for (i, d_str) in fields[9].split(',').filter(|s| !s.is_empty()).enumerate() {
+                        let pos = row_start + i as u32;
+                        if pos < win_start || pos >= win_end {
+                            continue;
+                        }
+
+                        let depth: u32 = d_str
+                            .parse()
+                            .with_context(|| format!("Invalid depth '{}' at pos {}", d_str, pos))?;
+
+                        let low_cov = depth < min_coverage;
+
+                        let is_multimap = match genome_index.unique_context_at(&GenomicPos {
+                            contig: contig.clone(),
+                            pos,
+                        }) {
+                            Some(UniqueContext { before, after }) => {
+                                let before_ambig = before.is_none_or(|k| k > half_median);
+                                let after_ambig = after.is_none_or(|k| k > half_median);
+                                before_ambig || after_ambig
+                            }
+                            None => true,
+                        };
+
+                        n_positions += 1;
+                        if low_cov || is_multimap {
+                            n_at_risk += 1;
+                        }
+                    }
+
+                    Ok(())
+                })?;
+
+                contig_at_risk += n_at_risk;
+                contig_positions += n_positions;
+
+                let idx = window_counter.fetch_add(1, std::sync::atomic::Ordering::Relaxed) + 1;
+
+                let pct_at_risk = if n_positions > 0 {
+                    100.0 * n_at_risk as f64 / n_positions as f64
+                } else {
+                    0.0
+                };
+                let global_pct = 100.0 * idx as f64 / total_windows as f64;
+
+                info!(
+                    "  [{:>4}/{:>4} windows | {:>5.1}% total] {contig}:{win_start}-{win_end} \
+                     → {n_at_risk}/{n_positions} at risk ({pct_at_risk:.1}%)",
+                    idx, total_windows, global_pct
+                );
+
+                windows.push((win_start, win_end, n_at_risk));
+
+                current_interval += interval;
+                if current_interval >= *size {
+                    break;
+                }
+            }
+
+            let contig_pct = if contig_positions > 0 {
+                100.0 * contig_at_risk as f64 / contig_positions as f64
+            } else {
+                0.0
+            };
+            info!(
+                "  └─ {contig} done: {contig_at_risk}/{contig_positions} at risk ({contig_pct:.1}%)"
+            );
+
+            Ok((contig.clone(), windows))
+        })
+        .collect::<anyhow::Result<Vec<_>>>()?;
+
+    // Écriture séquentielle dans l'ordre original des contigs
+    let file = File::create(output_path.as_ref())
+        .with_context(|| format!("Failed to create: {}", output_path.as_ref().display()))?;
+    let mut writer = BufWriter::new(file);
+
+    writeln!(writer, "region\tn_at_risk")?;
+
+    // Réordonne selon l'ordre original des contigs
+    for (contig, _) in &contigs {
+        if let Some((_, windows)) = results.iter().find(|(c, _)| c == contig) {
+            for (win_start, win_end, n_at_risk) in windows {
+                writeln!(writer, "{contig}:{win_start}-{win_end}\t{n_at_risk}")?;
+            }
+        }
+    }
+
+    writer.flush()?;
+    info!(
+        "sv_quality_index complete → {}",
+        output_path.as_ref().display()
+    );
+
+    Ok(())
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::helpers::test_init;
+    use crate::{helpers::test_init, uniqueness::GenomeIndex};
 
     #[test]
     fn bam_stats() -> anyhow::Result<()> {
@@ -1365,4 +1625,26 @@ mod tests {
         // println!("{stats}");
         Ok(())
     }
+
+    #[test]
+    fn bam_quality() -> anyhow::Result<()> {
+        test_init();
+
+        let config = Config::default();
+
+        let genome_index = GenomeIndex::build(
+            "/home/t_steimle/ref/hs1/chm13v2.0.fa",
+            "/home/t_steimle/ref/hs1/hs1_uniqueness",
+            None, // None = tout indexer
+        )?;
+        sv_quality_index(
+            "DUMCO",
+            "diag",
+            &config,
+            &genome_index,
+            15,
+            "DUMCO_sv_qual.tsv",
+        )?;
+        Ok(())
+    }
 }

src/config.rs

@@ -820,12 +820,16 @@ impl Config {
         )
     }
 
-    /// Normal count directory: `<normal_dir>/counts`.
+    pub fn dir_count(&self, id: &str, time: &str) -> String {
+        format!("{}/{}", self.solo_dir(id, time), self.count_dir_name)
+    }
+
+    /// Normal count directory: `<normal_dir>/<count_dir_name>`.
     pub fn normal_dir_count(&self, id: &str) -> String {
         format!("{}/{}", self.normal_dir(id), self.count_dir_name)
     }
 
-    /// Tumor count directory: `<tumoral_dir>/counts`.
+    /// Tumor count directory: `<tumoral_dir>/<count_dir_name>`.
     pub fn tumoral_dir_count(&self, id: &str) -> String {
         format!("{}/{}", self.tumoral_dir(id), self.count_dir_name)
     }

src/io/readers.rs

@@ -14,7 +14,8 @@
 use std::{
     fs::{self, File},
     io::{BufRead, BufReader, Read, Write},
-    path::Path, sync::Arc,
+    path::Path,
+    sync::Arc,
 };
 
 use anyhow::Context;
@@ -220,9 +221,7 @@ impl TabixReader<SmbReader> {
             .index_cache_dir()
             .context("no index cache dir configured for tabix index")?;
 
-        let local_tbi_path = cache_dir.join(
-            bgz_rel_path.replace(['/', '\\'], "_") + ".tbi",
-        );
+        let local_tbi_path = cache_dir.join(bgz_rel_path.replace(['/', '\\'], "_") + ".tbi");
 
         fs::write(&local_tbi_path, &tbi_bytes)
             .with_context(|| format!("write cached tabix index: {}", local_tbi_path.display()))?;
@@ -288,7 +287,16 @@ impl<R: std::io::Read + std::io::Seek> TabixReader<R> {
             return Ok(());
         }
 
+        // eprintln!(
+        //     "region: {}:{}-{}",
+        //     rname, region.range.start, region.range.end
+        // );
+        // eprintln!("ref_seq_id: {}", ref_seq_id);
+        // eprintln!("num chunks: {}", chunks.len());
+
         for chunk in chunks {
+            // eprintln!("  chunk: start={:?} end={:?}", chunk.start(), chunk.end());
+
             self.reader
                 .seek(chunk.start())
                 .context("BGZF seek failed")?;
@@ -326,7 +334,9 @@ impl<R: std::io::Read + std::io::Seek> TabixReader<R> {
 ///
 /// Returns an error if the index or file cannot be read.
 pub fn fetch_tabix_lines<R: std::io::Read + std::io::Seek>(
-    reader: &mut TabixReader<R>, region: &GenomeRange) -> anyhow::Result<Vec<String>> {
+    reader: &mut TabixReader<R>,
+    region: &GenomeRange,
+) -> anyhow::Result<Vec<String>> {
     let mut lines = Vec::new();
     reader.fetch_with(region, |line| {
         lines.push(line.to_owned());

src/scan/scan.rs

@@ -60,13 +60,14 @@ use rust_htslib::bam::{self, IndexedReader, Read, Record};
 
 use crate::helpers::{bam_contigs, get_genome_sizes, is_file_older};
 use crate::io::bam::fb_inv_from_record;
-use crate::io::readers::get_gz_reader;
+use crate::io::readers::{get_gz_reader, TabixReader};
 use crate::io::tsv::TsvLine;
 use crate::io::writers::{BgzTabixWriter, IndexFormat};
 use crate::math::filter_outliers_modified_z_score_with_indices;
 
 use crate::config::Config;
 use crate::pipes::{Initialize, ShouldRun};
+use crate::positions::GenomeRange;
 use crate::runners::Run;
 use crate::scan::bin::{Bin, BinStats};
 use crate::variant::vcf_variant::Label;
@@ -345,6 +346,51 @@ impl fmt::Display for BinOutlier {
     }
 }
 
+pub struct CountsReader<R> {
+    reader: TabixReader<R>,
+}
+
+impl<R: std::io::Read + std::io::Seek> CountsReader<R> {
+    pub fn from_tabix_reader(reader: TabixReader<R>) -> Self {
+        Self { reader }
+    }
+
+    pub fn depths(&mut self, region: &GenomeRange) -> anyhow::Result<Vec<u32>> {
+        let mut out: Vec<u32> = Vec::new();
+        self.reader
+            .fetch_with(region, |line| -> anyhow::Result<()> {
+                let bin_count = BinCount::from_tsv_row(line)?;
+                if bin_count.contig != region.contig()
+                    || bin_count.end <= region.range.start
+                    || bin_count.start >= region.range.end
+                {
+                    return Ok(());
+                }
+
+                // overlap window in absolute genomic coordinates
+                let overlap_start = bin_count.start.max(region.range.start);
+                let overlap_end = bin_count.end.min(region.range.end);
+
+                // offsets into bin_count.depths, assuming depths[i] corresponds to position bin_count.start + i
+                let lo = (overlap_start - bin_count.start) as usize;
+                let hi = (overlap_end - bin_count.start) as usize;
+
+                anyhow::ensure!(
+                    hi <= bin_count.depths.len(),
+                    "depths array shorter than expected: bin {}-{}, depths.len()={}, hi={}",
+                    bin_count.start,
+                    bin_count.end,
+                    bin_count.depths.len(),
+                    hi
+                );
+
+                out.extend_from_slice(&bin_count.depths[lo..hi]);
+                Ok(())
+            })?;
+        Ok(out)
+    }
+}
+
 // Definitions matching your *current* bin enumeration (bitwise identical ranges)
 #[derive(Clone, Copy)]
 struct BinDef {
@@ -1116,3 +1162,31 @@ impl Label for SomaticScan {
         "Somatic Scan".to_string()
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use crate::{
+        helpers::test_init, io::readers::TabixReader, positions::GenomeRange,
+        scan::scan::CountsReader,
+    };
+
+    #[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 mut creader = CountsReader::from_tabix_reader(reader);
+
+        let region = GenomeRange::from_1_inclusive(contig, start0, end0);
+        let u = creader.depths(&region)?;
+        println!("{u:?}");
+
+        Ok(())
+    }
+}

src/uniqueness/mod.rs

@@ -989,19 +989,6 @@ mod tests {
 
         writer.flush()?;
 
-        // for bias in [AnchorBias::Left, AnchorBias::Right, AnchorBias::Auto] {
-        //     match index.query(&pos, bias) {
-        //         Some(iv) => println!(
-        //             "{bias:?}: {}:{}-{} ({}bp)",
-        //             pos.contig,
-        //             iv.start,
-        //             iv.end,
-        //             iv.len()
-        //         ),
-        //         None => println!("{bias:?}: non-unique"),
-        //     }
-        // }
-
         Ok(())
     }
 }