uniqueness update with RC

src/callers/deep_somatic.rs

@@ -95,15 +95,28 @@ use rust_htslib::bam::{self, Read};
 use uuid::Uuid;
 
 use crate::{
-    annotation::{Annotation, Annotations, Caller, CallerCat, Sample}, collection::vcf::Vcf, commands::{
-        AnyJob, Command as JobCommand, LocalBatchRunner, LocalRunner, SbatchRunner, SlurmParams, SlurmRunner, bcftools::{BcftoolsConcat, BcftoolsIndex, BcftoolsKeepPass}, exec_jobs
-    }, config::Config, helpers::{
+    annotation::{Annotation, Annotations, Caller, CallerCat, Sample},
+    collection::vcf::Vcf,
+    commands::{
+        bcftools::{BcftoolsConcat, BcftoolsIndex, BcftoolsKeepPass},
+        exec_jobs, AnyJob, Command as JobCommand, LocalBatchRunner, LocalRunner, SbatchRunner,
+        SlurmParams, SlurmRunner,
+    },
+    config::Config,
+    helpers::{
         get_genome_sizes_in_header_order, is_file_older, remove_dir_if_exists,
         singularity_bind_flags, split_ordered_genome_into_n_regions_exact,
-    }, io::vcf::read_vcf, jobs_seq, locker::SampleLock, pipes::{Initialize, ShouldRun, Version}, run, runners::Run, variant::{
+    },
+    io::vcf::read_vcf,
+    jobs_seq,
+    locker::SampleLock,
+    pipes::{Initialize, ShouldRun, Version},
+    run,
+    runners::Run,
+    variant::{
         variant_collection::VariantCollection,
         vcf_variant::{Label, Variants},
-    }
+    },
 };
 
 /// DeepSomatic paired tumor-normal somatic variant caller.
@@ -545,7 +558,8 @@ fn merge_deepsomatic_parts(base: &DeepSomatic, n_parts: usize) -> anyhow::Result
         vec![jobs_seq![concat]],
         base.config.slurm_runner,
         base.config.slurm_max_par.into(),
-    ).context("Failed to run bcftools concat for DeepSomatic parts")?;
+    )
+    .context("Failed to run bcftools concat for DeepSomatic parts")?;
 
     fs::rename(&final_tmp, &final_passed_vcf)
         .context("Failed to rename merged DeepSomatic PASS VCF")?;
@@ -627,7 +641,11 @@ pub fn run_deepsomatic_chunked(id: &str, config: &Config, n_parts: usize) -> any
             job.log_dir = format!("{}/part{}", base.log_dir, i + 1);
             info!("Planned DeepSomatic job:\n{job}");
 
-            let pass = BcftoolsKeepPass::from_config(&job.config, job.output_vcf_path(), job.passed_vcf_path());
+            let pass = BcftoolsKeepPass::from_config(
+                &job.config,
+                job.output_vcf_path(),
+                job.passed_vcf_path(),
+            );
 
             Ok(jobs_seq![job, pass])
         })

src/uniqueness/genome_index.rs

@@ -8,8 +8,8 @@ use anyhow::Context;
 
 use super::{AnchorBias, CacheKey, UniqueInterval, UniquenessIndex};
 use crate::{
-    io::fasta::{read_fai, read_single_contig_fasta, split_fasta, FaiEntry},
-    uniqueness::cache::load_from_cache,
+    io::fasta::{FaiEntry, read_fai, read_single_contig_fasta, split_fasta},
+    uniqueness::{UniqueContext, cache::load_from_cache},
 }; // ton module existant
 
 #[derive(Debug, Clone, PartialEq, Eq)]
@@ -55,19 +55,28 @@ impl GenomeIndex {
             }
 
             let key = CacheKey::from_fai(&entry.name, entry.len, fasta_path)?;
-            let cache_path = cache_dir.join(key.file_name());
-
-            if cache_path.exists() {
-                match load_from_cache(&cache_path, entry.len) {
-                    Ok(min_unique_len) => {
+            let fwd_path = cache_dir.join(key.file_name());
+            let rev_path = cache_dir.join(format!("rev.{}", key.file_name()));
+
+            if fwd_path.exists() && rev_path.exists() {
+                match (
+                    load_from_cache(&fwd_path, entry.len),
+                    load_from_cache(&rev_path, entry.len),
+                ) {
+                    (Ok(min_unique_len), Ok(rev_min_unique_len)) => {
                         info!("[genome_index] cache hit  : {}", entry.name);
                         cached.insert(
                             entry.name.clone(),
-                            UniquenessIndex::from_raw(min_unique_len),
+                            UniquenessIndex::from_raw(min_unique_len, rev_min_unique_len),
                         );
                     }
-                    Err(e) => {
-                        warn!("[genome_index] cache miss : {} ({e})", entry.name);
+                    (fwd, rev) => {
+                        warn!(
+                            "[genome_index] cache miss : {} (fwd={:?}, rev={:?})",
+                            entry.name,
+                            fwd.err(),
+                            rev.err()
+                        );
                         need_build.push(entry);
                     }
                 }
@@ -116,12 +125,6 @@ impl GenomeIndex {
         Ok(Self { contigs: cached })
     }
 
-    pub fn query(&self, pos: &GenomicPos, bias: AnchorBias) -> Option<UniqueInterval> {
-        self.contigs
-            .get(&pos.contig)?
-            .minimal_unique_interval_containing(pos.pos, bias)
-    }
-
     pub fn contig_names(&self) -> impl Iterator<Item = &str> {
         self.contigs.keys().map(String::as_str)
     }
@@ -129,4 +132,14 @@ impl GenomeIndex {
     pub fn contig_len(&self, contig: &str) -> Option<usize> {
         self.contigs.get(contig).map(|idx| idx.genome_len())
     }
+
+    pub fn unique_context_at(&self, pos: &GenomicPos) -> Option<UniqueContext> {
+        Some(self.contigs.get(&pos.contig)?.unique_context_at(pos.pos))
+    }
+
+    pub fn query(&self, pos: &GenomicPos, bias: AnchorBias) -> Option<UniqueInterval> {
+        self.contigs
+            .get(&pos.contig)?
+            .minimal_unique_interval_containing(pos.pos, bias)
+    }
 }

src/uniqueness/mod.rs

@@ -1,7 +1,75 @@
-use log::{info, warn};
-use std::fs;
-use std::io::{self, BufReader, BufWriter, Read, Write};
-use std::path::{Path, PathBuf};
+//! Genome uniqueness index based on suffix arrays.
+//!
+//! This module computes, stores, and queries minimal unique sequence lengths
+//! across a genome.
+//!
+//! For each genomic position `pos`, the index stores the shortest sequence
+//! starting at `pos` that occurs exactly once in the genome:
+//!
+//! ```text
+//! genome[pos .. pos + k)
+//! ```
+//!
+//! where `k` is the minimal unique length at `pos`.
+//!
+//! Internally, uniqueness is computed from a suffix array and LCP
+//! (Longest Common Prefix) array, allowing construction in `O(n log n)`
+//! time and constant-time uniqueness queries.
+//!
+//! # Query types
+//!
+//! The index supports two complementary views of uniqueness:
+//!
+//! ## Unique intervals containing a position
+//!
+//! Given a genomic position, the index can return a unique interval that
+//! contains that position.
+//!
+//! ```text
+//! start <= pos < end
+//! ```
+//!
+//! Several intervals may satisfy this condition; [`AnchorBias`] controls
+//! how a candidate interval is selected.
+//!
+//! ## Independent left/right uniqueness
+//!
+//! The index can also report how much sequence is required on either side
+//! of a position to obtain uniqueness independently.
+//!
+//! For a position `pos`:
+//!
+//! ```text
+//! before: genome[pos - k .. pos)
+//! after : genome[pos .. pos + k)
+//! ```
+//!
+//! [`UniqueContext`] returns the smallest unique sequence ending at `pos`
+//! (`before`) and the smallest unique sequence starting at `pos` (`after`).
+//!
+//! To support efficient left-side queries, the index stores uniqueness
+//! information for both the forward genome and its reverse complement-free
+//! reversal.
+//!
+//! # Caching
+//!
+//! Construction can be expensive for large genomes. The module therefore
+//! supports disk-backed caching of precomputed uniqueness vectors.
+//!
+//! Cached data are keyed by genome metadata and loaded transparently when
+//! available.
+//!
+//! # Main types
+//!
+//! - [`UniquenessIndex`] — uniqueness index for a single sequence.
+//! - [`GenomeIndex`] — multi-contig genome-wide index.
+//! - [`UniqueInterval`] — unique interval containing a queried position.
+//! - [`UniqueContext`] — independent left/right uniqueness lengths.
+//! - [`GenomicPos`] — genomic coordinate used by [`GenomeIndex`].
+
+use log::warn;
+use std::io::{self};
+use std::path::Path;
 
 mod cache;
 mod genome_index;
@@ -35,6 +103,14 @@ fn build_min_unique_len(genome: &[u8]) -> Vec<usize> {
     mul
 }
 
+pub struct UniqueContext {
+    /// Bases required before `pos`.
+    pub before: Option<usize>,
+
+    /// Bases required from `pos` onward.
+    pub after: Option<usize>,
+}
+
 #[derive(Debug, Clone, PartialEq, Eq)]
 pub struct UniqueInterval {
     pub start: usize,
@@ -56,6 +132,7 @@ pub enum AnchorBias {
 
 pub struct UniquenessIndex {
     min_unique_len: Vec<usize>,
+    rev_min_unique_len: Vec<usize>,
     genome_len: usize,
 }
 
@@ -64,33 +141,83 @@ impl UniquenessIndex {
     /// Cache lives in `cache_dir` (created if absent).
     /// If the cache file exists and is valid, loading is O(n) read — no recomputation.
     pub fn build_with_cache(genome: &[u8], key: &CacheKey, cache_dir: &Path) -> io::Result<Self> {
-        let cache_path = cache_dir.join(key.file_name());
-        if cache_path.exists() {
-            match load_from_cache(&cache_path, genome.len()) {
-                Ok(mul) => return Ok(Self::from_raw(mul)),
-                Err(e) => warn!("[cache] invalid ({e}), rebuilding..."),
+        std::fs::create_dir_all(cache_dir)?;
+
+        let fwd_path = cache_dir.join(key.file_name());
+        let rev_path = cache_dir.join(format!("rev.{}", key.file_name()));
+
+        let min_unique_len = match load_from_cache(&fwd_path, genome.len()) {
+            Ok(mul) => mul,
+            Err(e) => {
+                if fwd_path.exists() {
+                    warn!("[cache] invalid forward ({e}), rebuilding...");
+                }
+
+                let mul = build_min_unique_len(genome);
+                save_to_cache(&fwd_path, &mul)?;
+                mul
             }
-        }
-        let mul = build_min_unique_len(genome);
-        save_to_cache(&cache_path, &mul)?;
-        Ok(Self::from_raw(mul))
+        };
+
+        let rev_min_unique_len = match load_from_cache(&rev_path, genome.len()) {
+            Ok(mul) => mul,
+            Err(e) => {
+                if rev_path.exists() {
+                    warn!("[cache] invalid reverse ({e}), rebuilding...");
+                }
+
+                let mut rev = genome.to_vec();
+                rev.reverse();
+
+                let mul = build_min_unique_len(&rev);
+                save_to_cache(&rev_path, &mul)?;
+                mul
+            }
+        };
+
+        Ok(Self::from_raw(min_unique_len, rev_min_unique_len))
     }
 
     /// Construct directly from a precomputed `min_unique_len` vector.
     /// Used when loading from disk cache — skips SA/LCP computation entirely.
-    pub fn from_raw(min_unique_len: Vec<usize>) -> Self {
+    pub fn from_raw(min_unique_len: Vec<usize>, rev_min_unique_len: Vec<usize>) -> Self {
+        assert_eq!(min_unique_len.len(), rev_min_unique_len.len());
+
         let genome_len = min_unique_len.len();
+
         Self {
             min_unique_len,
+            rev_min_unique_len,
             genome_len,
         }
     }
 
     /// Build without cache (e.g. for tests or one-shot use).
     pub fn build(genome: &[u8]) -> Self {
-        Self::from_raw(build_min_unique_len(genome))
+        let min_unique_len = build_min_unique_len(genome);
+
+        let mut rev = genome.to_vec();
+        rev.reverse();
+        let rev_min_unique_len = build_min_unique_len(&rev);
+
+        Self::from_raw(min_unique_len, rev_min_unique_len)
     }
 
+    /// Returns the number of nucleotides before and after `pos` needed to form
+    /// a unique sequence containing `pos`.
+    ///
+    /// The returned flank defines the half-open interval:
+    ///
+    /// `[pos - front, pos + back)`
+    ///
+    /// where `front` is the number of nucleotides before `pos`, and `back` is the
+    /// number of nucleotides from `pos` onward, including the nucleotide at `pos`.
+    ///
+    /// Returns `None` if no unique sequence containing `pos` can be found.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `pos >= self.genome_len()`.
     pub fn minimal_unique_interval_containing(
         &self,
         pos: usize,
@@ -161,6 +288,44 @@ impl UniquenessIndex {
     pub fn genome_len(&self) -> usize {
         self.genome_len
     }
+
+    /// Returns how much sequence is needed before and after `pos`
+    /// to obtain uniqueness independently on each side.
+    ///
+    /// `before = Some(k)` means:
+    ///
+    /// `genome[pos - k .. pos)` is the shortest unique sequence ending at `pos`.
+    ///
+    /// `after = Some(k)` means:
+    ///
+    /// `genome[pos .. pos + k)` is the shortest unique sequence starting at `pos`.
+    ///
+    /// Returns `None` for a side if no unique sequence exists on that side.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `pos >= self.genome_len()`.
+    pub fn unique_context_at(&self, pos: usize) -> UniqueContext {
+        assert!(pos < self.genome_len);
+
+        let after = match self.min_unique_len[pos] {
+            usize::MAX => None,
+            k => Some(k),
+        };
+
+        let before = if pos == 0 {
+            None
+        } else {
+            let rev_pos = self.genome_len - pos;
+
+            match self.rev_min_unique_len[rev_pos] {
+                usize::MAX => None,
+                k => Some(k),
+            }
+        };
+
+        UniqueContext { before, after }
+    }
 }
 
 // tests/uniqueness.rs
@@ -558,7 +723,7 @@ mod tests {
         let index = GenomeIndex::build(
             "/home/t_steimle/ref/hs1/chm13v2.0.fa",
             "/home/t_steimle/ref/hs1/hs1_uniqueness",
-            None, // None = tout indexer
+            Some(&["chr1"]), // None = tout indexer
         )?;
 
         let pos = GenomicPos {