SamtoolsMergeMany, slurm_pipe/import upate

pandora-config.example.toml

@@ -26,7 +26,7 @@ docker_max_memory_go = 400
 db_cases_path = "/data/cases.sqlite"
 
 # Path to the conda activation script.
-conda_sh = "/data/miniconda3/etc/profile.d/conda.sh"
+conda_sh = "/mnt/beegfs02/software/recherche/miniconda/25.1.1/etc/profile.d/conda.sh"
 
 
 #######################################
@@ -253,7 +253,7 @@ germline_phased_vcf = "{result_dir}/{id}/diag/{id}_variants_constit_phased.vcf.g
 #######################################
 
 # Path to Severus script.
-severus_bin = "/data/tools/MySeverus/severus.py"
+severus_bin = " /home/t_steimle/somatic_pipe_tools/Severus/severus.py"
 
 # Force Severus recomputation.
 severus_force = false
@@ -262,10 +262,10 @@ severus_force = false
 severus_threads = 32
 
 # VNTRs BED for Severus.
-vntrs_bed = "/data/ref/hs1/vntrs_chm13.bed"
+vntrs_bed = "/home/t_steimle/ref/hs1/vntrs_chm13.bed"
 
 # Path of the Severus panel of normals.
-severus_pon = "/data/ref/hs1/PoN_1000G_chm13.tsv.gz"
+severus_pon = "/home/t_steimle/ref/hs1/PoN_1000G_chm13.tsv.gz"
 
 # Paired Severus output directory.
 # {result_dir}, {id}

src/callers/clairs.rs

@@ -2,7 +2,9 @@ use crate::{
     annotation::{Annotation, Annotations, Caller, CallerCat, Sample},
     collection::vcf::Vcf,
     commands::{
-        CapturedOutput, Command as JobCommand, Runner as LocalRunner, SbatchRunner, SlurmParams, SlurmRunner, bcftools::{BcftoolsConcat, BcftoolsKeepPass}, run_many_sbatch
+        bcftools::{BcftoolsConcat, BcftoolsKeepPass},
+        run_many_sbatch, CapturedOutput, Command as JobCommand, Runner as LocalRunner,
+        SbatchRunner, SlurmParams, SlurmRunner,
     },
     config::Config,
     helpers::{
@@ -25,7 +27,7 @@ use regex::Regex;
 use rust_htslib::bam::{self, Read};
 use std::{
     fmt, fs,
-    path::Path,
+    path::{Path, PathBuf},
     process::{Command as ProcessCommand, Stdio},
 };
 
@@ -42,9 +44,9 @@ pub struct ClairS {
     pub id: String,
     pub config: Config,
     pub log_dir: String,
-    /// Optional list of regions passed as repeated `-r REGION` args.
-    /// When empty, ClairS runs genome-wide.
-    pub regions: Vec<String>,
+    /// Optional region passed as `-r REGION` (`ctg:start-end`).
+    /// When `None`, ClairS runs genome-wide.
+    pub region: Option<String>,
 
     /// Optional part index for chunked parallel runs (1-indexed).
     ///
@@ -58,12 +60,8 @@ impl fmt::Display for ClairS {
         writeln!(f, "🧬 ClairS {}", self.id)?;
         writeln!(
             f,
-            "  Regions   : {}",
-            if self.regions.is_empty() {
-                "genome-wide".to_string()
-            } else {
-                format!("{} regions", self.regions.len())
-            }
+            "  Region   : {}",
+            self.region.as_deref().unwrap_or("genome-wide")
         )?;
         writeln!(
             f,
@@ -86,7 +84,7 @@ impl Initialize for ClairS {
             id,
             log_dir,
             config: config.clone(),
-            regions: Vec::new(),
+            region: None,
             part_index: None,
         };
 
@@ -131,15 +129,11 @@ impl JobCommand for ClairS {
         let output_dir = self.part_output_dir();
 
         // Build repeated -r REGION args if any regions were set (for batched runs)
-        let region_args = if self.regions.is_empty() {
-            String::new()
-        } else {
-            self.regions
-                .iter()
-                .map(|r| format!("-r {r}"))
-                .collect::<Vec<_>>()
-                .join(" ")
-        };
+        let region_arg = self
+            .region
+            .as_ref()
+            .map(|r| format!("-r {r}"))
+            .unwrap_or_default();
 
         let sample_name = format!("{}_diag", self.id);
 
@@ -161,7 +155,7 @@ impl JobCommand for ClairS {
              --use_longphase_for_intermediate_haplotagging true \
              --output_dir {output_dir} \
              -s {sample_name} \
-             {region_args}",
+             {region_arg}",
             singularity_bin = self.config.singularity_bin,
             image = self.config.clairs_image,
             tumor_bam = self.config.tumoral_bam(&self.id),
@@ -171,7 +165,7 @@ impl JobCommand for ClairS {
             platform = self.config.clairs_platform,
             output_dir = output_dir,
             sample_name = sample_name,
-            region_args = region_args,
+            region_arg = region_arg,
         )
     }
 }
@@ -585,8 +579,6 @@ impl Version for ClairS {
 
 /// Merge N chunked ClairS PASS VCFs into the final clairs_passed_vcf().
 fn merge_clairs_parts(base: &ClairS, n_parts: usize) -> anyhow::Result<()> {
-    use std::path::PathBuf;
-
     let mut part_pass_paths: Vec<PathBuf> = Vec::with_capacity(n_parts);
 
     for i in 1..=n_parts {
@@ -646,15 +638,15 @@ pub fn run_clairs_chunked_sbatch_with_merge(
         bam::Reader::from_path(&normal_bam).with_context(|| format!("Opening BAM {normal_bam}"))?;
     let header = bam::Header::from_template(reader.header());
     let genome_sizes = get_genome_sizes(&header)?;
-    let region_chunks = split_genome_into_n_regions(&genome_sizes, n_parts);
-    let n_parts = region_chunks.len();
 
-    // Build jobs
-    let mut jobs = Vec::with_capacity(n_parts);
-    for (i, regions) in region_chunks.into_iter().enumerate() {
+    let regions = split_genome_into_n_regions(&genome_sizes, n_parts);
+    let n_jobs = regions.len();
+
+    let mut jobs = Vec::with_capacity(n_jobs);
+    for (i, region) in regions.into_iter().enumerate() {
         let mut job = base.clone();
         job.part_index = Some(i + 1);
-        job.regions = regions;
+        job.region = Some(region);
         job.log_dir = format!("{}/part{}", base.log_dir, i + 1);
         info!("Planned ClairS job:\n{job}");
         jobs.push(job);
@@ -662,7 +654,7 @@ pub fn run_clairs_chunked_sbatch_with_merge(
 
     // Run all parts via Slurm
     let outputs = run_many_sbatch(jobs)?;
-
+    
     // Merge somatic PASS VCFs into final clairs_passed_vcf()
     merge_clairs_parts(&base, n_parts)?;
 
@@ -692,6 +684,8 @@ mod tests {
         let u = ClairS::initialize("34528", &config)?;
         info!("{u}");
 
+        let _ = run_clairs_chunked_sbatch_with_merge("34528", &config, 5)?;
+
         Ok(())
     }
 }

src/commands/samtools.rs

@@ -147,16 +147,15 @@ impl super::Command for SamtoolsMerge {
         let (intersection, frac_into) = into_qset.intersect(&bam_qset);
         let n_shared = intersection.len();
 
-        // relaxed for recovery reads
-        // if n_shared > 100_000 {
-        //     anyhow::bail!(
-        //         "Cannot merge: {} and {} share {} read name(s) ({:.4} fraction of `into`)",
-        //         self.into.display(),
-        //         self.bam.display(),
-        //         n_shared,
-        //         frac_into,
-        //     );
-        // }
+        if n_shared > 0 {
+            anyhow::bail!(
+                "Cannot merge: {} and {} share {} read name(s) ({:.4} fraction of `into`)",
+                self.into.display(),
+                self.bam.display(),
+                n_shared,
+                frac_into,
+            );
+        }
 
         let dir = self.into.parent().context(format!(
             "Failed to find parent dir of: {}",
@@ -259,6 +258,194 @@ impl super::SlurmRunner for SamtoolsMerge {
     }
 }
 
+/// Merges multiple BAM files into a single output BAM.
+///
+/// Uses the first input BAM as the header source while including
+/// reads from all input BAMs.
+pub struct SamtoolsMergeMany {
+    /// Path to samtools binary
+    pub bin: String,
+    /// Number of threads for samtools merge
+    pub threads: usize,
+    /// Output merged BAM path
+    pub into: PathBuf,
+    /// Input BAM files (first BAM provides the header)
+    pub bams: Vec<PathBuf>,
+    /// Pipeline configuration
+    pub config: Config,
+}
+
+impl super::Command for SamtoolsMergeMany {
+    /// Prepares the filesystem and validates preconditions for the multi-merge.
+    ///
+    /// This method:
+    /// 1. Ensures at least one input BAM is provided
+    /// 2. Ensures all input BAMs exist
+    /// 3. Creates the output directory if needed
+    /// 4. Ensures no QNAME overlap across any pair of BAMs
+    /// 5. Saves the union QNAME set next to the output BAM
+    fn init(&mut self) -> anyhow::Result<()> {
+        if self.bams.is_empty() {
+            anyhow::bail!("At least one input BAM is required");
+        }
+
+        if self.into.exists() {
+            anyhow::bail!(
+                "Destination BAM already exists: {}",
+                self.into.display()
+            );
+        }
+
+        // Create output directory if needed
+        if let Some(parent) = self.into.parent() {
+            if !parent.as_os_str().is_empty() {
+                fs::create_dir_all(parent).with_context(|| {
+                    format!("Failed to create output directory: {}", parent.display())
+                })?;
+            }
+        }
+
+        // Check that all input BAMs exist
+        for bam in &self.bams {
+            anyhow::ensure!(
+                bam.exists(),
+                "BAM file doesn't exist: {}",
+                bam.display()
+            );
+        }
+
+        // Skip QNAME collision check if only one BAM
+        if self.bams.len() == 1 {
+            let qset = QNameSet::load_or_create(
+                &self.bams[0],
+                self.config.bam_min_mapq,
+                self.threads,
+            )?;
+            qset.save(WGSBamStats::qnames_path_from_bam(&self.into))?;
+            return Ok(());
+        }
+
+        // QNAME uniqueness check across all BAMs
+        let mut iter = self.bams.iter();
+        let first_bam = iter.next().expect("checked non-empty above");
+
+        let mut union_qset = QNameSet::load_or_create(
+            first_bam,
+            self.config.bam_min_mapq,
+            self.threads,  // Fixed: was hardcoded to 4
+        )
+        .with_context(|| {
+            format!("Failed to load or create qnames for {}", first_bam.display())
+        })?;
+
+        for bam in iter {
+            let bam_qset = QNameSet::load_or_create(
+                bam,
+                self.config.bam_min_mapq,
+                self.threads,  // Fixed: was hardcoded to 4
+            )
+            .with_context(|| {
+                format!("Failed to load or create qnames for {}", bam.display())
+            })?;
+
+            let (intersection, frac) = union_qset.intersect(&bam_qset);
+            let n_shared = intersection.len();
+
+            if n_shared > 0 {
+                anyhow::bail!(
+                    "Cannot merge: {} shares {} read name(s) ({:.2}%) with previous BAMs",
+                    bam.display(),
+                    n_shared,
+                    frac * 100.0,
+                );
+            }
+
+            union_qset.merge(&bam_qset);
+        }
+
+        // Save the union QNAME set next to the OUTPUT BAM
+        let qnames_path = WGSBamStats::qnames_path_from_bam(&self.into);
+        union_qset.save(&qnames_path).with_context(|| {
+            format!("Failed to save union qnames to {}", qnames_path.display())
+        })?;
+
+        Ok(())
+    }
+
+    /// Builds the `samtools merge` command line.
+    ///
+    /// # Command Format
+    ///
+    /// ```text
+    /// samtools merge -@ {threads} -c -p -h {header_bam} {output} {input1} {input2} ...
+    /// ```
+    ///
+    /// - `-@`: Number of compression threads
+    /// - `-c`: Combine @RG headers with colliding IDs (use first occurrence)
+    /// - `-p`: Combine @PG headers with colliding IDs (use first occurrence)
+    /// - `-h`: Use header from specified BAM (reads from that file are ignored for header)
+    ///
+    /// The first input BAM provides the header. All input BAMs (including the first)
+    /// contribute reads to the output.
+    fn cmd(&self) -> String {
+        let header_bam = self
+            .bams
+            .first()
+            .expect("bams should be non-empty in cmd()");
+
+        let inputs = self
+            .bams
+            .iter()
+            .map(|p| p.display().to_string())
+            .collect::<Vec<_>>()
+            .join(" ");
+
+        // -h header_bam: use header from first BAM
+        // inputs: all BAMs including first (for reads)
+        format!(
+            "{bin} merge -@ {threads} -c -p -h {header} {output} {inputs}",
+            bin = self.bin,
+            threads = self.threads,
+            header = header_bam.display(),
+            output = self.into.display(),
+            inputs = inputs,
+        )
+    }
+
+    fn clean_up(&self) -> anyhow::Result<()> {
+        Ok(())
+    }
+}
+
+impl SamtoolsMergeMany {
+    /// Creates a new merge command from config defaults.
+    pub fn from_config(into: PathBuf, bams: Vec<PathBuf>, config: &Config) -> Self {
+        Self {
+            bin: config.align.samtools_bin.clone(),
+            threads: config.bam_n_threads as usize,
+            into,
+            bams,
+            config: config.clone(),
+        }
+    }
+}
+
+
+impl super::Runner for SamtoolsMergeMany {}
+
+impl super::SlurmRunner for SamtoolsMergeMany {
+    fn slurm_args(&self) -> Vec<String> {
+        SlurmParams {
+            job_name: Some("samtools_merge_many".into()),
+            cpus_per_task: Some(self.threads as u32),
+            mem: Some("60G".into()),
+            partition: Some("shortq".into()),
+            gres: None,
+        }
+        .to_args()
+    }
+}
+
 #[derive(Debug)]
 pub struct SamtoolsSplit {
     /// Path to the `samtools` executable.

src/helpers.rs

@@ -795,10 +795,14 @@ pub fn get_genome_sizes(header: &rust_htslib::bam::Header) -> anyhow::Result<FxH
 
 /// Split genome into ~n_parts equal-sized chunks (by number of bases),
 /// returning for each chunk a list of regions in `ctg:start-end` form.
+/// Split genome into ~n_parts equal-sized chunks (by bases).
+/// Each chunk is a single region `ctg:start-end` (no cross-contig regions).
+///
+/// Note: the number of regions returned will be ≈ n_parts (may differ by 1–2).
 pub fn split_genome_into_n_regions(
     genome_sizes: &FxHashMap<String, u64>,
     n_parts: usize,
-) -> Vec<Vec<String>> {
+) -> Vec<String> {
     if n_parts == 0 || genome_sizes.is_empty() {
         return Vec::new();
     }
@@ -815,47 +819,23 @@ pub fn split_genome_into_n_regions(
         return Vec::new();
     }
 
-    let target_chunk_size: u64 = total_bases.div_ceil(n_parts as u64); // ceil
+    let target_chunk_size: u64 = (total_bases + n_parts as u64 - 1) / n_parts as u64; // ceil
 
-    let mut chunks: Vec<Vec<String>> = vec![Vec::new(); n_parts];
-    let mut current_part = 0usize;
-    let mut remaining_in_part = target_chunk_size;
+    let mut regions = Vec::new();
 
     for (ctg, len) in contigs {
-        let mut remaining_ctg = len;
+        let mut remaining = len;
         let mut start: u64 = 1;
 
-        while remaining_ctg > 0 && current_part < n_parts {
-            let take = remaining_in_part.min(remaining_ctg);
+        while remaining > 0 {
+            let take = remaining.min(target_chunk_size);
             let end = start + take - 1;
+            regions.push(format!("{ctg}:{start}-{end}"));
 
-            chunks[current_part].push(format!("{ctg}:{start}-{end}"));
-
-            remaining_ctg -= take;
+            remaining -= take;
             start = end + 1;
-            remaining_in_part -= take;
-
-            if remaining_in_part == 0 {
-                current_part += 1;
-                if current_part >= n_parts {
-                    break;
-                }
-                remaining_in_part = target_chunk_size;
-            }
         }
-
-        // If we ran out of parts but contig still has bases, dump the rest into the last part
-        if remaining_ctg > 0 && current_part >= n_parts {
-            let end = start + remaining_ctg - 1;
-            chunks[n_parts - 1].push(format!("{ctg}:{start}-{end}"));
-            break;
-        }
-    }
-
-    // Remove empty chunks at the end (e.g. n_parts > total_bases)
-    while chunks.last().is_some_and(|c| c.is_empty()) {
-        chunks.pop();
     }
 
-    chunks
+    regions
 }

src/pipes/somatic_slurm.rs

@@ -8,7 +8,7 @@ use crate::{
     commands::{
         dorado::{DoradoAlign, DoradoBasecall},
         run_many_sbatch,
-        samtools::{SamtoolsIndex, SamtoolsMerge, SamtoolsSort, SamtoolsSplit},
+        samtools::{SamtoolsIndex, SamtoolsMerge, SamtoolsMergeMany, SamtoolsSort, SamtoolsSplit},
         SlurmRunner,
     },
     config::Config,
@@ -457,6 +457,13 @@ pub fn import_run(run: &Pod5sRun, config: &Config) -> anyhow::Result<()> {
             // remove unsorted chunk
             fs::remove_file(&bam)?;
 
+            let mut index_cmd = SamtoolsIndex {
+                bin: config.align.samtools_bin.clone(),
+                threads: config.align.samtools_view_threads,
+                bam: sorted_bam.to_string_lossy().to_string(),
+            };
+            SlurmRunner::run(&mut index_cmd)?;
+
             // replace path in case_bam_map with sorted version
             *bam = sorted_bam;
         }
@@ -481,6 +488,11 @@ pub fn import_run(run: &Pod5sRun, config: &Config) -> anyhow::Result<()> {
             continue;
         };
 
+        if aligned_bams.is_empty() {
+            warn!("Aligned BAM list is empty for case {}", case.case_id);
+            continue;
+        }
+
         let final_bam_path = PathBuf::from(&config.result_dir)
             .join(&case.case_id)
             .join(&sample_type_dir)
@@ -491,39 +503,106 @@ pub fn import_run(run: &Pod5sRun, config: &Config) -> anyhow::Result<()> {
 
         fs::create_dir_all(final_bam_path.parent().unwrap())?;
 
-        for bam in aligned_bams {
-            if final_bam_path.exists() {
-                let mut index_cmd = SamtoolsIndex {
-                    bin: config.align.samtools_bin.clone(),
-                    threads: config.align.samtools_view_threads,
-                    bam: bam.to_string_lossy().to_string(),
-                };
-                SlurmRunner::run(&mut index_cmd)?;
-                let mut index_cmd = SamtoolsIndex {
-                    bin: config.align.samtools_bin.clone(),
-                    threads: config.align.samtools_view_threads,
-                    bam: final_bam_path.to_string_lossy().to_string(),
-                };
-                SlurmRunner::run(&mut index_cmd)?;
-
-                // Merge into existing - clean_up() removes source bam
-                let mut merge_cmd = SamtoolsMerge::from_config(config, bam, &final_bam_path);
-                SlurmRunner::run(&mut merge_cmd)?;
-            } else {
-                // First BAM becomes the base
-                fs::rename(bam, &final_bam_path)?;
-
-                // Index the final BAM
-                let mut index_cmd = SamtoolsIndex {
-                    bin: config.align.samtools_bin.clone(),
-                    threads: config.align.samtools_view_threads,
-                    bam: final_bam_path.to_string_lossy().to_string(),
-                };
-                SlurmRunner::run(&mut index_cmd)?;
+        // -----------------------------------------------------------------
+        // 5a) Merge all chunk BAMs for this case into a single per-case BAM
+        //     using SamtoolsMergeMany (header from first chunk).
+        // -----------------------------------------------------------------
+
+        let case_merged_bam: PathBuf = if aligned_bams.len() == 1 {
+            // Only one chunk: use it directly
+            aligned_bams[0].clone()
+        } else {
+            let tmp_case_merged = tmp_dir.join(format!(
+                "{}_{}_{}_merged.bam",
+                Uuid::new_v4(),
+                case.case_id,
+                sample_type_dir
+            ));
+
+            info!(
+                "  Merging {} chunk BAMs for case {} → {}",
+                aligned_bams.len(),
+                case.case_id,
+                tmp_case_merged.display()
+            );
+
+            let mut merge_many_cmd = SamtoolsMergeMany::from_config(
+                tmp_case_merged.clone(),
+                aligned_bams.clone(),
+                config,
+            );
+            SlurmRunner::run(&mut merge_many_cmd)?;
+
+            // We can safely remove the chunk BAMs now (they were all inputs).
+            for bam in aligned_bams {
+                if bam != &tmp_case_merged {
+                    if let Err(e) = fs::remove_file(bam) {
+                        warn!(
+                            "Failed to remove chunk BAM {} after merge: {e}",
+                            bam.display()
+                        );
+                    }
+                }
             }
+
+            tmp_case_merged
+        };
+
+        // -----------------------------------------------------------------
+        // 5b) Merge the per-case BAM into the final BAM for this case
+        //     using SamtoolsMerge (incremental merge).
+        // -----------------------------------------------------------------
+
+        if final_bam_path.exists() {
+            info!(
+                "  Final BAM already exists for case {}, merging into existing: {}",
+                case.case_id,
+                final_bam_path.display()
+            );
+
+            // Index both source (per-case) and destination BAMs
+            let mut index_cmd = SamtoolsIndex {
+                bin: config.align.samtools_bin.clone(),
+                threads: config.align.samtools_view_threads,
+                bam: case_merged_bam.to_string_lossy().to_string(),
+            };
+            SlurmRunner::run(&mut index_cmd)?;
+
+            let mut index_cmd = SamtoolsIndex {
+                bin: config.align.samtools_bin.clone(),
+                threads: config.align.samtools_view_threads,
+                bam: final_bam_path.to_string_lossy().to_string(),
+            };
+            SlurmRunner::run(&mut index_cmd)?;
+
+            // Merge into existing final BAM.
+            // SamtoolsMerge::clean_up() will remove the source BAM.
+            let mut merge_cmd =
+                SamtoolsMerge::from_config(config, &case_merged_bam, &final_bam_path);
+            SlurmRunner::run(&mut merge_cmd)?;
+        } else {
+            info!(
+                "  Creating new final BAM for case {} → {}",
+                case.case_id,
+                final_bam_path.display()
+            );
+
+            // First per-case BAM becomes the base final BAM.
+            fs::rename(&case_merged_bam, &final_bam_path)?;
+
+            // Index the new final BAM.
+            let mut index_cmd = SamtoolsIndex {
+                bin: config.align.samtools_bin.clone(),
+                threads: config.align.samtools_view_threads,
+                bam: final_bam_path.to_string_lossy().to_string(),
+            };
+            SlurmRunner::run(&mut index_cmd)?;
         }
 
-        // Sort the merged final BAM
+        // -----------------------------------------------------------------
+        // 5c) Sort and index the final BAM
+        // -----------------------------------------------------------------
+
         let sorted_tmp = final_bam_path.with_extension("sorted.bam");
 
         info!(
@@ -547,6 +626,63 @@ pub fn import_run(run: &Pod5sRun, config: &Config) -> anyhow::Result<()> {
         SlurmRunner::run(&mut index_cmd)?;
 
         info!("  Output: {}", final_bam_path.display());
+
+        // for bam in aligned_bams {
+        //     if final_bam_path.exists() {
+        //         let mut index_cmd = SamtoolsIndex {
+        //             bin: config.align.samtools_bin.clone(),
+        //             threads: config.align.samtools_view_threads,
+        //             bam: bam.to_string_lossy().to_string(),
+        //         };
+        //         SlurmRunner::run(&mut index_cmd)?;
+        //         let mut index_cmd = SamtoolsIndex {
+        //             bin: config.align.samtools_bin.clone(),
+        //             threads: config.align.samtools_view_threads,
+        //             bam: final_bam_path.to_string_lossy().to_string(),
+        //         };
+        //         SlurmRunner::run(&mut index_cmd)?;
+        //
+        //         // Merge into existing - clean_up() removes source bam
+        //         let mut merge_cmd = SamtoolsMerge::from_config(config, bam, &final_bam_path);
+        //         SlurmRunner::run(&mut merge_cmd)?;
+        //     } else {
+        //         // First BAM becomes the base
+        //         fs::rename(bam, &final_bam_path)?;
+        //
+        //         // Index the final BAM
+        //         let mut index_cmd = SamtoolsIndex {
+        //             bin: config.align.samtools_bin.clone(),
+        //             threads: config.align.samtools_view_threads,
+        //             bam: final_bam_path.to_string_lossy().to_string(),
+        //         };
+        //         SlurmRunner::run(&mut index_cmd)?;
+        //     }
+        // }
+
+        // Sort the merged final BAM
+        // let sorted_tmp = final_bam_path.with_extension("sorted.bam");
+        //
+        // info!(
+        //     "  Sorting final BAM for case {} → {}",
+        //     case.case_id,
+        //     final_bam_path.display()
+        // );
+        //
+        // let mut sort_cmd = SamtoolsSort::from_config(config, &final_bam_path, &sorted_tmp);
+        // SlurmRunner::run(&mut sort_cmd)?;
+        //
+        // // Replace unsorted BAM with sorted BAM
+        // fs::rename(&sorted_tmp, &final_bam_path)?;
+        //
+        // // Index the **sorted** final BAM
+        // let mut index_cmd = SamtoolsIndex {
+        //     bin: config.align.samtools_bin.clone(),
+        //     threads: config.align.samtools_view_threads,
+        //     bam: final_bam_path.to_string_lossy().to_string(),
+        // };
+        // SlurmRunner::run(&mut index_cmd)?;
+        //
+        // info!("  Output: {}", final_bam_path.display());
     }
 
     if let Some(d) = tmp_split_dir {