pandora_lib_promethion/src/lib.rs

//! # 🧬 Long-read Somatic Variant Calling and Analysis Framework
//!
//! This Rust library provides a modular, parallelizable framework for somatic variant calling, annotation, and interpretation from long-read sequencing data. It is designed to support full pipelines for research and clinical workflows across multiple variant callers and analysis stages.
//!
//! The library also serves as an extensible platform that developers can leverage to add custom features, integrate new tools, and tailor workflows to specific use cases.
//!
//! ## 🧩 Key Features
//!
//! - **POD5 Demultiplexing and Alignment**: End-to-end support for processing ONT POD5 files:
//!     - Barcode-aware demultiplexing using metadata CSVs
//!     - POD5 subsetting and organization by case
//!     - Integration with basecallers (e.g., Dorado) for read alignment
//! - **Pipeline Management**: Full orchestration of Dockerized execution pipelines for tools such as ClairS, Nanomonsv, DeepVariant, Savana, Modkit, and Severus.
//! - **Flexible Configuration**: Centralized configuration system (`Config`, `CollectionsConfig`) for all modules and pipelines.
//! - **Input Abstraction**: Unified handling of BAM, POD5, and VCF file collections across cohorts and directories.
//! - **Variant Processing**: Modular loading, filtering, statistical analysis, and annotation of somatic and germline variants.
//! - **Haplotype Phasing and Methylation**: Support for LongPhase-based phasing and Modkit methylation pileups with support for multi-threaded pileup and aggregation.
//! - **Parallel Execution**: Uses `rayon` for efficient multicore parallelization over large cohorts and tasks.
//!
//! ## 📚 Module Highlights
//!
//! - `callers`: Interfaces to variant calling tools (ClairS, DeepVariant, Nanomonsv, Savana, etc...)
//! - `runners`: Pipeline runners (e.g. `Somatic`, `SeverusSolo`, `LongphasePhase`) that manage end-to-end execution.
//! - `collection`: Organizes input data across BAMs, VCFs, and POD5 files with auto-detection of completed runs.
//! - `annotation`: VEP line parsing and high-level annotation aggregation.
//! - `pipes`: Composition modules for executing pipelines across callers and post-processing steps.
//! - `functions`: Custom logic for genome assembly, entropy estimation, and internal tooling.
//! - `positions`, `variant`, `helpers`: Utilities for SV modeling, variant filtering, position overlap logic, and helper methods.
//!
//! ## ⚡ Workflow Overview
//!
//! ### 1. 📦 From POD5 to BAM Alignment
//!
//! - **Demultiplexing**: POD5 files are subset and demuxed using barcodes (via CSV metadata).
//! - **Flowcell Case Management**: Each sample is identified by a [`collection::pod5::FlowCellCase`] containing its ID, time point, and POD5 directory.
//! - **Alignment**: The [`commands::dorado::Dorado`] module handles alignment of POD5 reads to reference genome, producing BAMs.
//!
//! ```rust
//! let case = FlowCellCase { id: "PATIENT1", time_point: "diag", barcode: "01", pod_dir: "...".into() };
//! Dorado::init(case, Config::default())?.run_pipe()?;
//! ```
//!
//! ### 2. 🧬 Variant Calling (BAM ➝ VCF)
//!
//! Using the aligned BAMs, multiple variant callers can be run in parallel. The [`callers`] and [`runners`] modules support:
//!
//! - **ClairS** – somatic small variant calling with LongPhase haplotagging  
//! - **Nanomonsv** – structural variants (SV)  
//! - **DeepVariant** – germline small variants  
//! - **Savana** – SVs and copy number variations (CNV)  
//! - **Modkit** – methylation pileups  
//! - **LongPhase** – phasing and modcalling
//!
//! All workflows can be triggered per-case or per-cohort using `Collections` or `Somatic` runners.
//!
//! ```rust
//! ClairS::initialize("PATIENT1", Config::default())?.run()?;
//! NanomonSV::initialize("PATIENT1", Config::default())?.run()?;
//! ```
//!
//! ### 3. 📈 Aggregation & Statistics (VCF ➝ JSON / Stats)
//!
//! After variant calling:
//!
//! - Annotate with VEP ([`annotation`] module)
//! - Load and filter with [`variant::variant_collection`]
//! - Compute variant and region-level stats (e.g., mutation rates, alteration categories, coding overlaps)
//!
//! ```rust
//! let variants = Variants::load_from_json("/path/to/somatic_variants.json.gz")?;
//! let stats = VariantsStats::new(&variants, "PATIENT1", &config)?;
//! stats.save_to_json("/output/path/stats.json.gz")?;
//! ```
//!
//! ### 4. 🧠 Intelligent Task Management (`collection` module)
//!
//! - Auto-discovers available samples, POD5s, BAMs, and VCFs
//! - Detects missing outputs and creates task lists
//! - Tasks are parallelizable using Rayon and can be run on-demand
//!
//! ```rust
//! let mut collections = Collections::new(CollectionsConfig::default())?;
//! collections.todo()?;      // Identify missing steps
//! collections.run()?;       // Run them automatically
//! ```
//!
//! ## 🔬 Testing
//!
//! Integration tests demonstrate the entire pipeline. Run with logging enabled:
//!
//! ```bash
//! export RUST_LOG=debug
//! cargo test -- --nocapture
//! ```
//!
//! ## 🧪 Example Use Cases
//!
//! - Full somatic variant calling pipeline on matched tumor/normal samples
//! - POD5-based pipeline from raw signal to variants
//! - Aggregation and annotation of SVs across a clinical cohort
//! - Methylation analysis using nanopore-specific tools
//! - Variant calling and analysis in large-scale longitudinal studies
//!
//! ## 🚀 Getting Started
//!
//! All workflows are initialized from `Config` and driven by the `Collections` structure:
//!
//! ```rust
//! let collections = Collections::new(CollectionsConfig::default())?;
//! collections.todo()?;
//! collections.run()?;
//! ```
//!
//! ## 🔗 References
//!
//! **Basecalling and alignment**
//! - Dorado: <https://github.com/nanoporetech/dorado>
//!
//! **Variants Callers**
//! - ClairS: <https://github.com/HKU-BAL/ClairS>
//! - Nanomonsv: <https://github.com/friend1ws/nanomonsv>
//! - Savana: <https://github.com/cortes-ciriano-lab/savana>
//! - DeepVariant: <https://github.com/google/deepvariant>
//! - DeepSomatic: <https://github.com/google/deepsomatic>
//! - LongPhase: <https://github.com/PorubskyResearch/LongPhase>
//! - Modkit: <https://github.com/nanoporetech/modkit>
//!
//! **Variants annotation**
//! - VEP: <https://www.ensembl.org/info/docs/tools/vep/index.html>
//!
//! ---

use std::sync::{Arc, Mutex};

pub mod aligner;
pub mod annotation;
pub mod callers;
pub mod collection;
pub mod commands;
pub mod config;
pub mod de_novo;
pub mod functions;
pub mod helpers;
pub mod io;
pub mod locker;
pub mod math;
pub mod pipes;
pub mod positions;
pub mod runners;
pub mod scan;
pub mod slurm_helpers;
pub mod uniqueness;
pub mod variant;

#[macro_use]
extern crate lazy_static;

// Define DOCKER_ID lock for handling Docker kill when ctrlc is pressed
lazy_static! {
    static ref DOCKER_ID: Arc<Mutex<Vec<String>>> = Arc::new(Mutex::new(Vec::new()));
    static ref TEST_DIR: String = "/mnt/beegfs02/scratch/t_steimle/test_data".to_string();
}

#[cfg(test)]
mod tests {
    use std::{collections::HashMap, path::Path};

    use annotation::{
        vep::{VepLine, VEP},
        Annotations,
    };

    use collection::bam::{counts_at, counts_ins_at, nt_pileup, WGSBam};
    use functions::assembler::{Assembler, AssemblerConfig};
    use helpers::estimate_shannon_entropy;
    use io::bed::read_bed;
    use log::{error, info};
    use positions::{overlaps_par, GenomePosition, GenomeRange};
    use rayon::prelude::*;
    use variant::{variant_collection, vcf_variant::VcfVariant};

    use self::{/* collection::pod5::{FlowCellCase, Pod5Collection}, */ config::Config};
    use super::*;
    use crate::{
        annotation::{Annotation, colorscb::annotate_colorsdb},
        collection::{
            bam::{self},
            flowcells::{FlowCells, scan_archive},
            vcf::VcfCollection,
        },
        helpers::find_files,
        io::{bed::bedrow_overlaps_par, dict::read_dict, gff::features_ranges},
        positions::{merge_overlapping_genome_ranges, range_intersection_par, sort_ranges},
        scan::scan::somatic_scan,
        variant::{
            variant_collection::{
                ExternalAnnotation, VariantCollection, group_variants_by_bnd_desc, group_variants_by_bnd_rc
            },
            variants_stats::{self, VariantsStats, somatic_depth_quality_ranges},
            vcf_variant::{AlterationCategory, BNDDesc, BNDGraph, ToBNDGraph},
        },
    };

    // export RUST_LOG="debug"
    fn init() {
        let _ = env_logger::Builder::from_env(env_logger::Env::default().default_filter_or("info"))
            .try_init();
    }

    // #[test]
    // fn run_dorado() -> anyhow::Result<()> {
    //     init();
    //     let case = FlowCellCase {
    //         id: "CONSIGNY".to_string(),
    //         time_point: "mrd".to_string(), barcode: "07".to_string(), pod_dir: "/mnt/beegfs01/ data/run_data/20240326-CL/CONSIGNY-MRD-NB07_RICCO-DIAG-NB08/20240326_1355_1E_PAU78333_bc25da25/pod5_pass/barcode07".into()
    //     };
    //     dorado::Dorado::init(case, Config::default())?.run_pipe()
    // }

    // #[test]
    // fn pod5() -> anyhow::Result<()> {
    //     let _ = env_logger::Builder::from_env(env_logger::Env::default().default_filter_or("info"))
    //         .build();
    //
    //     let coll = Pod5Collection::new(
    //         "/data/run_data",
    //         "/data/flow_cells.tsv",
    //         "/data/longreads_basic_pipe",
    //     )?;
    //     println!("{coll:#?}");
    //     // let runs = Runs::import_dir("/home/prom/store/banana-pool/run_data", "/data/flow_cells.tsv")?;
    //     Ok(())
    // }

    #[test]
    fn bam() -> anyhow::Result<()> {
        init();
        let bam_collection = bam::load_bam_collection("/data/longreads_basic_pipe");

        bam_collection
            .bams
            .iter()
            // .filter(|b| matches!(b.bam_type, BamType::Panel(_)))
            .for_each(|b| println!("{b:#?}"));

        let u = bam_collection.get("PARACHINI", "mrd");
        println!("{u:#?}");

        Ok(())
    }

    #[test]
    fn vcf() -> anyhow::Result<()> {
        init();
        let mut vcf_collection = VcfCollection::new("/data/longreads_basic_pipe");
        vcf_collection.sort_by_id();
        vcf_collection
            .vcfs
            .iter()
            .for_each(|v| v.println().unwrap());

        Ok(())
    }

    // pod5 view -I /data/run_data/20240903-CL/ARMEM-DG-N02_ASSJU-DG-N03/20240903_1428_1B_PAW47629_fc24c3cf/pod5/PAW47629_fc24c3cf_77b07847_0.pod5 | head -5000 | awk '{if(NR==1){print "target,"$0}else{print "subset_1.pod5,"$0}}' > /tmp/subset_ids.csv
    // pod5 subset /data/run_data/20240903-CL/ARMEM-DG-N02_ASSJU-DG-N03/20240903_1428_1B_PAW47629_fc24c3cf/pod5/PAW47629_fc24c3cf_77b07847_0.pod5 --csv /tmp/subset_ids.csv -o /data/test_suite/pod5/muxed/
    // #[test]
    // fn mux() -> anyhow::Result<()> {
    //     init();
    //     let result_dir = "/data/test_suite/results".to_string();
    //     let cases = vec![
    //         FlowCellCase { id: "test_02".to_string(), time_point: "diag".to_string(), barcode: "02".to_string(), pod_dir: "/data/test_suite/pod5/muxed".into() },
    //         FlowCellCase { id: "test_03".to_string(), time_point: "diag".to_string(), barcode: "03".to_string(), pod_dir: "/data/test_suite/pod5/muxed".into() },
    //     ];
    //
    //     cases.iter().for_each(|c| {
    //         let dir = format!("{result_dir}/{}", c.id);
    //         if Path::new(&dir).exists() {
    //             fs::remove_dir_all(dir).unwrap();
    //         }
    //     });
    //     let config = Config { result_dir, ..Default::default()  };
    //     Dorado::from_mux(cases, config)
    // }

    // #[test_log::test]
    // fn clairs() -> anyhow::Result<()> {
    //     let config = ClairSConfig {
    //         result_dir: "/data/test".to_string(),
    //         ..ClairSConfig::default()
    //     };
    //     ClairS::new("test_a", "/data/test_data/subset.bam", "/data/test_data/subset_mrd.bam", config).run()
    // }

    // #[test]
    // fn nanomonsv() -> anyhow::Result<()> {
    //     init();
    //     let id = "HAMROUNE";
    //     NanomonSV::initialize(id, Config::default())?.run()
    // }
    //
    // #[test]
    // fn nanomonsv_version() -> anyhow::Result<()> {
    //     init();
    //     let v = NanomonSV::version(&Config::default())?;
    //     println!("NanomonSV version: {v}");
    //     let v = DeepVariant::version(&Config::default())?;
    //     println!("DeepVariant version: {v}");
    //     let v = Savana::version(&Config::default())?;
    //     println!("Savana version: {v}");
    //     let v = DeepSomatic::version(&Config::default())?;
    //     println!("DeepSomatic version: {v}");
    //     let v = ClairS::version(&Config::default())?;
    //     println!("ClairS version: {v}");
    //     Ok(())
    // }
    //
    // #[test]
    // fn nanomonsv_solo() -> anyhow::Result<()> {
    //     init();
    //     NanomonSVSolo::initialize("LAKHDHAR", "diag", Config::default())?.run()
    // }

    #[test]
    fn run_assemblers() -> anyhow::Result<()> {
        Assembler::new(
            "CAMEL".to_string(),
            "diag".to_string(),
            AssemblerConfig::default(),
        )
        .run()
    }

    // #[test]
    // fn run_dmr_par() -> anyhow::Result<()> {
    //     init();
    //     let collections = Collections::new(
    //         CollectionsConfig::default()
    //     )?;
    //     let tasks = collections.todo_dmr_c_diag_mrd();
    //     tasks.iter().for_each(|t| info!("{t}"));
    //     let len = tasks.len();
    //     // let pool = ThreadPoolBuilder::new().num_threads(10).build().unwrap();
    //     // pool.install(|| {
    //     //     tasks.par_iter().enumerate().for_each(|(i, t)| {
    //     //         let config = ModkitConfig {threads: 2,  ..Default::default() };
    //     //         if let collection::CollectionsTasks::DMRCDiagMrd { id, .. } = t { let _ = dmr_c_mrd_diag(id, &config); }
    //     //         println!("⚡ {i}/{len}");
    //     //     });
    //     // });
    //     Ok(())
    // }

    // #[test]
    // fn run_severus() -> anyhow::Result<()> {
    //     init();
    //     Severus::initialize("BANGA", Config::default())?.run()
    // }
    //
    // #[test]
    // fn run_severus_solo() -> anyhow::Result<()> {
    //     init();
    //     SeverusSolo::initialize("LAKHDHAR", "diag", Config::default())?.run()
    // }
    //
    // #[test]
    // fn check_versions() -> anyhow::Result<()> {
    //     init();
    //     let config = Config::default();
    //     let v = Savana::version(&config)?;
    //     info!("Savanna version {v}");
    //     let v = Severus::version(&config)?;
    //     info!("Severus version {v}");
    //     Ok(())
    // }

    // #[test]
    // fn run_deepvariant() -> anyhow::Result<()> {
    //     init();
    //     DeepVariant::initialize("HAMROUNE", "diag", Config::default())?.run()
    // }

    // #[test]
    // fn run_clairs() -> anyhow::Result<()> {
    //     init();
    //     ClairS::initialize("ADJAGBA", Config::default())?.run()
    // }

    // #[test]
    // fn run_longphase() -> anyhow::Result<()> {
    //     init();
    //     let id = "BECERRA";
    //     let diag_bam = format!("/data/longreads_basic_pipe/{id}/diag/{id}_diag_hs1.bam");
    //     let vcf = format!("/data/longreads_basic_pipe/{id}/diag/ClairS/clair3_normal_tumoral_germline_output.vcf.gz");
    //     let mrd_bam = format!("/data/longreads_basic_pipe/{id}/mrd/{id}_mrd_hs1.bam");
    //
    //     LongphaseHap::new(id, &diag_bam, &vcf, LongphaseConfig::default()).run()?;
    //     LongphaseHap::new(id, &mrd_bam, &vcf, LongphaseConfig::default()).run()
    // }

    // #[test]
    // fn run_longphase_modcall() -> anyhow::Result<()> {
    //     init();
    //     let id = "ADJAGBA";
    //     let time = "diag";
    //     LongphaseModcallSolo::initialize(id, time, Config::default())?.run()
    // }
    //
    // #[test]
    // fn run_longphase_phase() -> anyhow::Result<()> {
    //     init();
    //     let id = "ADJAGBA";
    //     LongphasePhase::initialize(id, Config::default())?.run()
    // }

    #[test]
    fn snv_parse() -> anyhow::Result<()> {
        init();
        let config = Config::default();
        // ClairS
        let row = "chr1\t10407\t.\tA\tG\t10.1\tPASS\tF\tGT:GQ:DP:AD:AF\t0/1:10:31:10,20:0.6452";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);

        let mut variant_col = VariantCollection {
            variants: vec![variant],
            vcf: collection::vcf::Vcf::new(
                "/data/longreads_basic_pipe/ACHITE/diag/ClairS/ACHITE_diag_clairs_PASSED.vcf.gz"
                    .into(),
            )?,
            caller: Annotation::Callers(annotation::Caller::ClairS, annotation::Sample::Somatic),
        };
        let annotations = Annotations::default();
        variant_col.annotate_with_constit_bam(
            &annotations,
            "/data/longreads_basic_pipe/ACHITE/mrd/ACHITE_mrd_hs1.bam",
            &config,
        )?;
        // DeepVariant
        let row =
            "chr1\t10407\t.\tA\tG\t7.4\tPASS\t.\tGT:GQ:DP:AD:VAF:PL\t1/1:5:9:1,8:0.888889:5,5,0";
        variant_col.variants.push(row.parse()?);

        let anns: Vec<Vec<Annotation>> = variant_col
            .variants
            .iter()
            .filter_map(|e| {
                annotations
                    .store
                    .get(&e.hash())
                    .map(|v| v.value().to_owned())
            })
            .collect();
        assert_eq!(anns[0], anns[1]);
        Ok(())
    }

    #[test]
    fn deletion_parse() -> anyhow::Result<()> {
        init();
        // Clairs
        let row = "chr1\t16760\t.\tAaag\tA\t8.48\tPASS\tF\tGT:GQ:DP:AD:AF\t0/1:8:39:22,17:0.4359";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        // case are not keeped
        assert_eq!(
            &var_string,
            "chr1\t16760\t.\tAAAG\tA\t8.48\tPASS\tF\tGT:GQ:DP:AD:AF\t0/1:8:39:22,17:0.4359"
        );
        assert_eq!(AlterationCategory::DEL, variant.alteration_category());

        let mut variant_col = VariantCollection {
            variants: vec![variant],
            vcf: collection::vcf::Vcf::new(
                "/data/longreads_basic_pipe/ACHITE/diag/ClairS/ACHITE_diag_clairs_PASSED.vcf.gz"
                    .into(),
            )?,
            caller: Annotation::Callers(annotation::Caller::ClairS, annotation::Sample::Somatic),
        };
        let annotations = Annotations::default();

        // DeepVariant, VAF is not an f32 at it should be, sometimes too long removed the last number for eq
        let row = "chr12\t326911\t.\tGTGTA\tG\t4.7\tPASS\t.\tGT:GQ:DP:AD:VAF:PL\t0/1:5:8:5,3:0.37500:2,0,21";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        variant_col.variants.push(row.parse()?);

        assert_eq!(AlterationCategory::DEL, variant.alteration_category());

        // Severus, 0000 added to VAF and hVAF
        let row = "chr10\t108974982\tseverus_DEL885\tN\t<DEL>\t60\tPASS\tPRECISE;SVTYPE=DEL;SVLEN=474;END=108975456;STRANDS=+-;INSIDE_VNTR=TRUE;MAPQ=60\tGT:GQ:VAF:hVAF:DR:DV\t0/1:61:0.30000:0.30000,0.00000,0.00000:7:3";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::DEL, variant.alteration_category());

        variant_col.variants.push(row.parse()?);

        // Nanomonsv dont parse last format remove: \t22:0 and nt putted in uppercase
        let row = "chr14\t5247080\td_106\tG\t<DEL>\t.\tPASS\tEND=5247255;SVTYPE=DEL;SVLEN=-175;SVINSLEN=39;SVINSSEQ=ATAACCCAGGTGATATAACACTTCTTTAGGCTCTGCCTA\tTR:VR\t18:3";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::DEL, variant.alteration_category());

        variant_col.variants.push(row.parse()?);

        let row  = "chr1\t20654667\t.\tTG\tT\t3.5\tPASS\t.\tGT:GQ:DP:AD:VAF:PL\t0/1:3:23:20,3:0.13044:0,0,16";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::DEL, variant.alteration_category());

        variant_col.variants.push(row.parse()?);

        let row = "chr1\t2160094\t.\tTCTGACAGCCTGGAACAGCACCCACAACCGCAGGTGAGCATCTGACAGCCCGCAGCAGCACCCACACGCACAGGTGAGAATCTGACAGCCCGGAGCAGCACCCACACAGGCAGGGGAGCATCTGACATCCTGGAGCAGCACCGACAACCCCAGGTGAGCAACTGAGAGCCTGGAACAGCACCCACACCCCCAGGTGAGAATCTGACAGCCTGGAAGAGCACCCCACATCCCCGGGTGAGCATCTGACAGCCTGGAACAGCATCAACACCCCCAGGTGAGCATCCGATAGCCTGGAGCAGCACCCACACCCTCAGGTGAGCATCTGACAGCCTGGAACAGCAACCACACCCCCAGGTGAACATCTGACAGCCCGGAGCAGCACCCACACCCCCAGGTGAGCATCTGACAGCCTGGAACAGCACCCACACCCCCAGGAGAGCATCCGGTAGCCTGGAGCAGAACCCACACCCACAGGCGAGCATCTGACAGCCTGGGTTGGCACCCACACCCCCAGGTGAGCATCTGATGGTCTGGAGCAGCACCCACACCTACAGGAGAGCATCTGACAACCTGCAACAGAACCCAAACCCCCAGGTGAGCATCTGACAGACTGGAACAGCACCCTGCACCCCCAGGTGAGTATCTGACGGCCTGGAACAGAACACACAAGCCCAGGTGAGCATCCGACAGCCTGGAGCAGCACCCACACCCCCAGGTGAGCATCAGACAGCCTGGAGTAGCACCCCACACCCCCAGGTGAGCATCCGACAGCCTGGAGCAGCACCCACACTCACCAGGTGAGCATGTGACATGCTGGAACAGCACCCACACCCCCAGGCGAGCATCTGACAGCCTGGAGCGGCACCCCACACCCCCAGGTGAGCATCGGACAGCCTGGAGCAGTACCCACACCCCCAGGTGAGCATCCGACAGCCTGGAACAGCACCCACACCCCCAGGTGAGCATCTGACAGACAGGAACAGCACCCACATGTCCAGGTGAGCCTCTGACAGACTGGAACAGCACGCGCACCCCCAGGTGAGCATCTGACAGGCTGGAACAGCACCCACACCCCCAGGAGAGCATCTTACAGCATGTAACAGCACCCACACACCCACGTGAGCATCTGACAGCCTGGAACAGCACCCTGCACCCCCAGGTGCGCACGTGACAGCCTGGAACAGCACCCACACCCCCAGGCGAGCATCGGACAGCCTGGAGCAGCACCCCACACCCGCAGGTGAGCATCCGACAGCCTGGAGCAGCACCCACACCCCCAGGTGAGCATGTGACAGCCTGGAACAGCACCCACACCCCCAGGCGAGCATCTGACAGCCTGGAGCAACACCCCACACCCCCAGGTGAGCATCGAACTGTCTGGAGCAGCACCCACAACCACAGGTGGGCATCGGAGAGAGTGGAGCAGCGCCCAGACCACCAGGCAAGCATCTGACAGCCTGGAGCAGTGCCCACACCCCCATGTGAGCATCTGACAGTATGGAGCAGCACCCACAGCCCAAGGTGAGCATCTGACAACCTGGAGCAGCACCCACACCCCCAGGCGAGCATCTGAACGCACAGAGCAGCACCCACACCCCCAGGCGAGCATCCGACAGCGTGGAGCAGCACCCACACTCCCAGGTGCGCGTGTGATGGTCTGGGGCAGCACCCACACACACAGGTGAGCCTGCGACAGCCTGGAGCAGCACCCACAGCCCCAGGTGAGCATCTGATGGTCTGGAGCAGCACCCACAACCACAGGTGAGCATCGGAGAGACTGGAGCAGCGCCGAAACCCCCAGGCGAACATTTGAGAGCCTGGAGCAGTGCCGACACCCCCAGGTGAGCATCTGACACCGTGGAGCAGCACCCACAGCCCAAGGTGAGCATCTGACAACCTGGAGCAGCACCCACAGCCCCAGGCGAGCATTTGAACGCACGGAGCAGGACCTACAGCCCCAGGCGAGCATCCGACAGCCTGGAGCAGCACCCACACACCCAGGTGAGCATCTGACAGCCTGGAGCAGCACCCACAACCCCAGGGGAGCATCTGACCGCATGGAATGGCATCCTCACCCGTAGGTGAACGTCCGACAGCCTGGAGCAGCACCCACACCCCCAGGTGAGCATCTGACAGCCTGGAACAGCACCTGCACCCCCGGGTGAGGATCAGATAGCCTGGAGCAGCACCCACACTCCAGGTGAGCATCTGACAGCCTGAAGCAGCACCCACACCAACAGGTGAGCATCTGACAGCCTGGAAGAGCACCCACAACCCCACGTGAGCATCTGACAGCCTGGAACAGGACCCTGCACCCCAAAGTGAGCATCTGACAACCTGGAGCAGGAACCACAACCCCAGGTGAGCATCTGATAGCCTGGAATAGCACCCACACACCCAGGTGAGCATCTGAGAGCCTGGAGCAGCACCCACACCCCCAGGTGAGCATCCGACAGCCTGGAACAGTACCCACACACCCAGGCGAGCATCTGACAGCCTGAAACAGCACACACACTCCCAGGTGAGCATCTCATATGCTGGAACAGCACCCACACCCCCAGGTGAGCATCTGACTGCCCGGAGCAGCACGCACACCCCCGGGTGAGCATCTGATAGCCTGGAACAGCACCCACACCCCCAGATGAGCATCCGACAGCCTGGAGCGGAGCCCACAGCCCCAGGCGAGCATCTGACAGCCTGGAACAGCACCCTGCATCCCCGGGTGAGGATCAGACAGCCTGGAGCAGCACCCACACTCCAGGTGAGCATCTGACAGCCTGAAGCAGCACCCACACCAACAGGTGAGCCTCTGACAGCCTGCAACAGCACCCACACCCCAAGGTGAGCATCTGACAGCCTGGAAGAGGACCCTGCAGCCCCAAGTGAGCATCTGACAACCTGAAGCAGCAACCACACTCCAAGGTGAGCATCCAACAGCCTGGAACAGCACCCACACACCCAGGTGAGCATCTGACAGCCTGGAGCAGCACCACACCCCCAGGTGAGCATCCGACAGCCTGGAACAGCACCTACAAACCCAGGAGAGCATCCGACAGCCTGGAGCAGCACCCACACCCCCAGGCGAGCATCTGACAGCCCGGAGCAGCACGCAAACCCCCAGGTGAGCATCTGATACCCTGGAACAGCACCCACACACCCAGGTGAGCATCCGACAGCCTAGAGCTGAACCCACACCAACAGGAGAGCATCTGACAGCCTGGGTCGGCACCCACACCCCCAGGTGAGCATCTGACAGCCTGGAACAGCACCCTGCACCCTCAGGTGCGCACGTGACAGCCTGGAACAGCACCCACACACCCAGGCGAGCATCTGACGGCCTGGAACGGCACCCACACCCCGAGGCGAGCATGGGACAGCCTGGAGAGCAGCCACACCCTCAGATGAGCATCTGACAGCCTGGAACAGCACCCTGCACCCCCAGGTGAGCATCTGACAGCCTGGAACAGGACCCACGTCCCCAGGCGAGCAA\tT\t3.1\tPASS\t.\tGT:GQ:DP:AD:VAF:PL\t0/1:3:15:11,4:0.26667:0,0,22";

        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::DEL, variant.alteration_category());

        variant_col.variants.push(row.parse()?);

        let row = "chr1\t27073397\t.\tTGGGG\tT\t28.9\tPASS\t.\tGT:GQ:DP:AD:VAF:PL\t./.:2:13:1,3:0.23077:24,2,26";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::DEL, variant.alteration_category());

        variant_col.variants.push(row.parse()?);

        let row = "chrY\t639648\t.\tCTTTTTT\tC\t10.4\tPASS\t.\tGT:GQ:DP:AD:VAF:PL\t0/1:1:7:0,2:0.28571:4,0,2";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::DEL, variant.alteration_category());

        variant_col.variants.push(row.parse()?);

        let row = "chr11\t31260259\tseverus_BND1077_1\tN\t]chr11:36214171]N\t60\tPASS\tPRECISE;SVTYPE=BND;SVLEN=4953912;MATE_ID=severus_BND1077_2;STRANDS=+-;MAPQ=60\tGT:GQ:VAF:hVAF:DR:DV\t0/1:68:0.67000:0.67000,0.00000,0.00000:4:8";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::DEL, variant.alteration_category());
        println!("{:#?}", variant.deletion_desc());

        variant_col.variants.push(row.parse()?);

        variant_col.annotate_with_constit_bam(
            &annotations,
            "/data/longreads_basic_pipe/ACHITE/mrd/ACHITE_mrd_hs1.bam",
            &Config::default(),
        )?;

        Ok(())
    }

    #[test]
    fn insertion_parse() -> anyhow::Result<()> {
        init();
        // Clairs
        let row = "chr1\t23005\t.\tT\tTCAC\t3.1\tPASS\tF\tGT:GQ:DP:AD:A\t0/1:3:13:9,4:0.3077";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(&var_string, row);
        assert_eq!(AlterationCategory::INS, variant.alteration_category());

        let mut variant_col = VariantCollection {
            variants: vec![variant],
            vcf: collection::vcf::Vcf::new(
                "/data/longreads_basic_pipe/ACHITE/diag/ClairS/ACHITE_diag_clairs_PASSED.vcf.gz"
                    .into(),
            )?,
            caller: Annotation::Callers(annotation::Caller::ClairS, annotation::Sample::Somatic),
        };
        let annotations = Annotations::default();

        // DeepVariant, VAF is not an f32 at it should be, sometimes too long removed the last number for eq
        let row = "chrX\t152732993\t.\tG\tGTTTTTTTTT\t18.9\tPASS\t.\tGT:GQ:DP:AD:VAF:PL\t1/1:3:6:0,5:0.50000:15,7,3";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        // variant_col.variants.push(row.parse()?);

        assert_eq!(AlterationCategory::INS, variant.alteration_category());

        // Severus, 0000 added to VAF and hVAF
        let row = "chr11\t26669932\tseverus_INS9403\tN\tT\t60\tPASS\tPRECISE;SVTYPE=INS;SVLEN=460;INSIDE_VNTR=TRUE;MAPQ=60\tGT:GQ:VAF:hVAF:DR:DV\t0/1:71:0.27000:0.27000,0.00000,0.00000:8:3";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::INS, variant.alteration_category());

        // variant_col.variants.push(row.parse()?);

        // Nanomonsv dont parse last format remove: \t22:0 and nt putted in uppercase
        let row = "chr8\t87940084\td_333\tT\t<INS>\t.\tPASS\tEND=87940201;SVTYPE=INS;SVINSLEN=172;SVINSSEQ=TA\tTR:VR\t9:5";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::INS, variant.alteration_category());

        // variant_col.variants.push(row.parse()?);

        //
        let row = "chr5\t36122736\t.\tT\tTGCTCCG\t3.9\tPASS\t.\tGT:GQ:DP:AD:VAF:PL\t0/1:4:28:11,16:0.57143:1,0,37";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);
        assert_eq!(AlterationCategory::INS, variant.alteration_category());

        variant_col.variants = Vec::new();
        variant_col.variants.push(row.parse()?);
        variant_col.annotate_with_constit_bam(
            &annotations,
            "/data/longreads_basic_pipe/PASSARD/mrd/PASSARD_mrd_hs1.bam",
            &Config::default(),
        )?;
        println!("{variant_col:?}");
        println!("{annotations:?}");

        Ok(())
    }

    #[test]
    fn trl_parse() -> anyhow::Result<()> {
        init();

        let row = "chr2\t207968575\tID_16420_1\ta\ta]chr11:41497080]\t.\tPASS\tSVTYPE=BND;MATEID=ID_16420_2;TUMOUR_READ_SUPPORT=4;TUMOUR_ALN_SUPPORT=4;NORMAL_READ_SUPPORT=0;NORMAL_ALN_SUPPORT=0;SVLEN=0;BP_NOTATION=++;SOURCE=SUPPLEMENTARY;CLUSTERED_READS_TUMOUR=4;CLUSTERED_READS_NORMAL=0;ORIGIN_STARTS_STD_DEV=0.433;ORIGIN_MAPQ_MEAN=56.25;ORIGIN_EVENT_SIZE_STD_DEV=0;ORIGIN_EVENT_SIZE_MEDIAN=0;ORIGIN_EVENT_SIZE_MEAN=0;END_STARTS_STD_DEV=17.754;END_MAPQ_MEAN=56.25;END_EVENT_SIZE_STD_DEV=0;END_EVENT_SIZE_MEDIAN=0;END_EVENT_SIZE_MEAN=0;TUMOUR_DP_BEFORE=8,11;TUMOUR_DP_AT=4,11;TUMOUR_DP_AFTER=4,11;NORMAL_DP_BEFORE=7,16;NORMAL_DP_AT=7,16;NORMAL_DP_AFTER=7,16;TUMOUR_AF=1,0.364;NORMAL_AF=0,0;TUMOUR_TOTAL_HP_AT=1,3,0;NORMAL_TOTAL_HP_AT=4,3,0;TUMOUR_ALT_HP=2,1,1;TUMOUR_PS=207946665;NORMAL_ALT_HP=0,0,0;CLASS=PREDICTED_SOMATIC\tGT\t0/1";
        let variant: VcfVariant = row.parse()?;
        // let var_string = variant.into_vcf_row();
        let u = variant.n_alt_depth();
        println!("{u:?}");

        Ok(())
    }

    #[test]
    fn dup_parse() -> anyhow::Result<()> {
        init();

        let row = "chr9\t91352364\tID_14667_1\tC\t]chr9:98584394]C\t.\tPASS\tSVTYPE=BND;MATEID=ID_14667_2;TUMOUR_READ_SUPPORT=3;TUMOUR_ALN_SUPPORT=3;NORMAL_READ_SUPPORT=0;NORMAL_ALN_SUPPORT=0;SVLEN=7232030;BP_NOTATION=-+;SOURCE=SUPPLEMENTARY;CLUSTERED_READS_TUMOUR=3;CLUSTERED_READS_NORMAL=0;ORIGIN_STARTS_STD_DEV=0.471;ORIGIN_MAPQ_MEAN=60;ORIGIN_EVENT_SIZE_STD_DEV=8.524;ORIGIN_EVENT_SIZE_MEDIAN=7232030;ORIGIN_EVENT_SIZE_MEAN=7232020;END_STARTS_STD_DEV=8.731;END_MAPQ_MEAN=60;END_EVENT_SIZE_STD_DEV=8.524;END_EVENT_SIZE_MEDIAN=7232030;END_EVENT_SIZE_MEAN=7232020;TUMOUR_DP_BEFORE=16,38;TUMOUR_DP_AT=20,33;TUMOUR_DP_AFTER=20,33;NORMAL_DP_BEFORE=26,22;NORMAL_DP_AT=26,22;NORMAL_DP_AFTER=26,22;TUMOUR_AF=0.15,0.091;NORMAL_AF=0,0;TUMOUR_TOTAL_HP_AT=7,12,1;NORMAL_TOTAL_HP_AT=13,13,0;TUMOUR_ALT_HP=0,2,1;TUMOUR_PS=90796185;NORMAL_ALT_HP=0,0,0;CLASS=PREDICTED_SOMATIC\tGT\t0/1";

        let variant: VcfVariant = row.parse()?;
        assert_eq!(AlterationCategory::DUP, variant.alteration_category());

        let row = "chr9\t98584394\tID_14667_2\tC\tC[chr9:91352364[\t.\tPASS\tSVTYPE=BND;MATEID=ID_14667_1;TUMOUR_READ_SUPPORT=3;TUMOUR_ALN_SUPPORT=3;NORMAL_READ_SUPPORT=0;NORMAL_ALN_SUPPORT=0;SVLEN=7232030;BP_NOTATION=-+;SOURCE=SUPPLEMENTARY;CLUSTERED_READS_TUMOUR=3;CLUSTERED_READS_NORMAL=0;ORIGIN_STARTS_STD_DEV=0.471;ORIGIN_MAPQ_MEAN=60;ORIGIN_EVENT_SIZE_STD_DEV=8.524;ORIGIN_EVENT_SIZE_MEDIAN=7232030;ORIGIN_EVENT_SIZE_MEAN=7232020;END_STARTS_STD_DEV=8.731;END_MAPQ_MEAN=60;END_EVENT_SIZE_STD_DEV=8.524;END_EVENT_SIZE_MEDIAN=7232030;END_EVENT_SIZE_MEAN=7232020;TUMOUR_DP_BEFORE=38,16;TUMOUR_DP_AT=33,20;TUMOUR_DP_AFTER=33,20;NORMAL_DP_BEFORE=22,26;NORMAL_DP_AT=22,26;NORMAL_DP_AFTER=22,26;TUMOUR_AF=0.091,0.15;NORMAL_AF=0,0;TUMOUR_TOTAL_HP_AT=15,18,0;NORMAL_TOTAL_HP_AT=11,10,1;TUMOUR_ALT_HP=1,0,2;TUMOUR_PS=98176815;NORMAL_ALT_HP=0,0,0;CLASS=PREDICTED_SOMATIC\tGT\t0/1";
        let variant: VcfVariant = row.parse()?;
        assert_eq!(AlterationCategory::DUP, variant.alteration_category());

        let row = "chr1\t9218455\tr_0\tC\t<DUP>\t.\tPASS\tSVTYPE=DUP;SVLEN=12572;END=9231027\tTR:VR\t37:9";

        let variant: VcfVariant = row.parse()?;
        println!("{:#?}", variant.alteration_category());
        println!("{:#?}", variant.bnd_desc());

        let row = "chr1\t47110295\tseverus_BND178_1\tN\t]chr1:47192169]N\t60.00\tPASS\tPRECISE;SVTYPE=BND;SVLEN=81874;MATE_ID=severus_BND178_2;STRANDS=+-;MAPQ=60\tGT:GQ:VAF:hVAF:DR:DV\t0/1:66:0.38000:0.38000,0.00000,0.00000:8:5";
        let variant: VcfVariant = row.parse()?;
        println!("{:#?}", variant.alteration_category());
        println!("{:#?}", variant.bnd_desc());

        Ok(())
    }

    #[test]
    fn variant_parse() -> anyhow::Result<()> {
        let row =
            "chr1\t1366\t.\tC\tCCCT\t8.2\tPASS\t.\tGT:GQ:DP:AD:VAF:PL\t1/1:4:6:1,4:0.66667:6,4,0";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);

        let row = "chr1\t1366\t.\tC\tCCCT\t8.2\tPASS\t.";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);

        let row = "chr1\t2628434\t.\tC\tT\t17.973\tPASS\tH;FAU=0;FCU=7;FGU=0;FTU=7;RAU=0;RCU=2;RGU=0;RTU=2\tGT:GQ:DP:AF:AD:NAF:NDP:NAD:AU:CU:GU:TU:NAU:NCU:NGU:NTU\t0/1:17:18:0.5:0,9:0:11:0,0:0:9:0:9:0:11:0:0";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);

        let row = "chr1\t52232\t.\tC\tCT\t18\t.\t.\tGT:GQ:DP:AD:AF\t1/.:1:24:3,5:0.208333";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);

        let row = "chr1\t52232\t.\tC\tCT\t18\t.\t.\tGT:GQ:DP:AD:AF\t1/1:1:24:3,5:0.208333";
        let variant_b: VcfVariant = row.parse()?;
        assert_eq!(variant, variant_b);
        let row = "chr1\t475157\t.\tA\tG\t12.301\tPASS\tH;FAU=2;FCU=0;FGU=2;FTU=0;RAU=3;RCU=0;RGU=3;RTU=0\tGT:GQ:DP:AF:AD:NAF:NDP:NAD:AU:CU:GU:TU:NAU:NCU:NGU:NTU\t0/1:12:10:0.5:0,5:0.0769:13:0,1:5:0:5:0:12:0:1:0";
        let variant: VcfVariant = row.parse()?;
        let var_string = variant.into_vcf_row();
        assert_eq!(row, &var_string);

        let row = "chr1\t161417408\tr_10_0\tT\t[chr1:161417447[TTGGCAGGTTCC\t.\tPASS\tSVTYPE=BND;MATEID=r_10_1;SVINSLEN=11;SVINSSEQ=TTGGCAGGTTC\tTR:VR\t22:3\t12:0";
        let variant: VcfVariant = row.parse()?;
        println!("{variant:#?}");
        let u = variant.bnd_desc();
        println!("{u:#?}");

        // Severus mates are not in RC
        let vcf = "chr7\t27304522\tseverus_BND6747_1\tN\t[chr6:32688062[N\t60\tPASS\tPRECISE;SVTYPE=BND;MATE_ID=severus_BND6747_2;STRANDS=--;MAPQ=60;CLUSTERID=severus_2\tGT:VAF:hVAF:DR:DV\t0/1:0.29:0.29,0,0:12:5";
        let variant: VcfVariant = vcf.parse()?;
        let bnd_a = variant.bnd_desc()?;

        let vcf = "chr6\t32688062\tseverus_BND6747_2\tN\t[chr7:27304522[N\t60\tPASS\tPRECISE;SVTYPE=BND;MATE_ID=severus_BND6747_1;STRANDS=--;MAPQ=60;CLUSTERID=severus_2 GT:VAF:hVAF:DR:DV\t0/1:0.29:0.29,0,0:12:5";
        let variant: VcfVariant = vcf.parse()?;
        let bnd_b = variant.bnd_desc()?;
        assert_eq!(bnd_a, bnd_b.rc());

        println!("{bnd_a}\n{bnd_b}");

        // Savana here each mate are in RC
        let vcf = "chr10\t102039096\tID_35957_2\tG\t]chr10:101973386]G\t.\tPASS\tSVTYPE=BND;MATEID=ID_35957_1;TUMOUR_READ_SUPPORT=7;TUMOUR_ALN_SUPPORT=7;NORMAL_READ_SUPPORT=0;NORMAL_ALN_SUPPORT=0;SVLEN=65710;BP_NOTATION=+-;SOURCE=SUPPLEMENTARY;CLUSTERED_READS_TUMOUR=7;CLUSTERED_READS_NORMAL=0;ORIGIN_STARTS_STD_DEV=0.35;ORIGIN_MAPQ_MEAN=60;ORIGIN_EVENT_SIZE_STD_DEV=7.248;ORIGIN_EVENT_SIZE_MEDIAN=65710;ORIGIN_EVENT_SIZE_MEAN=65705.4;END_STARTS_STD_DEV=7.007;END_MAPQ_MEAN=60;END_EVENT_SIZE_STD_DEV=7.248;END_EVENT_SIZE_MEDIAN=65710;END_EVENT_SIZE_MEAN=65705.4;TUMOUR_DP_BEFORE=38,29;TUMOUR_DP_AT=44,21;TUMOUR_DP_AFTER=44,21;NORMAL_DP_BEFORE=13,15;NORMAL_DP_AT=13,15;NORMAL_DP_AFTER=13,15;TUMOUR_AF=0.159,0.333;NORMAL_AF=0,0;TUMOUR_TOTAL_HP_AT=20,16,8;NORMAL_TOTAL_HP_AT=6,7,0;TUMOUR_ALT_HP=0,1,6;TUMOUR_PS=101917152;NORMAL_ALT_HP=0,0,0;CLASS=PREDICTED_SOMATIC\tGT\t0/1";
        let variant: VcfVariant = vcf.parse()?;
        let bnd_a = variant.bnd_desc()?;

        let vcf = "chr10\t101973386\tID_35957_1\tA\tA[chr10:102039096[\t.\tPASS\tSVTYPE=BND;MATEID=ID_35957_2;TUMOUR_READ_SUPPORT=7;TUMOUR_ALN_SUPPORT=7;NORMAL_READ_SUPPORT=0;NORMAL_ALN_SUPPORT=0;SVLEN=65710;BP_NOTATION=+-;SOURCE=SUPPLEMENTARY;CLUSTERED_READS_TUMOUR=7;CLUSTERED_READS_NORMAL=0;ORIGIN_STARTS_STD_DEV=0.35;ORIGIN_MAPQ_MEAN=60;ORIGIN_EVENT_SIZE_STD_DEV=7.248;ORIGIN_EVENT_SIZE_MEDIAN=65710;ORIGIN_EVENT_SIZE_MEAN=65705.4;END_STARTS_STD_DEV=7.007;END_MAPQ_MEAN=60;END_EVENT_SIZE_STD_DEV=7.248;END_EVENT_SIZE_MEDIAN=65710;END_EVENT_SIZE_MEAN=65705.4;TUMOUR_DP_BEFORE=29,38;TUMOUR_DP_AT=21,44;TUMOUR_DP_AFTER=21,44;NORMAL_DP_BEFORE=15,13;NORMAL_DP_AT=15,13;NORMAL_DP_AFTER=15,13;TUMOUR_AF=0.333,0.159;NORMAL_AF=0,0;TUMOUR_TOTAL_HP_AT=17,0,4;NORMAL_TOTAL_HP_AT=5,7,3;TUMOUR_ALT_HP=0,6,1;TUMOUR_PS=101917152;NORMAL_ALT_HP=0,0,0;CLASS=PREDICTED_SOMATIC\tGT\t0/1";
        let variant: VcfVariant = vcf.parse()?;
        let bnd_b = variant.bnd_desc()?;
        assert_eq!(bnd_a, bnd_b);

        println!("{bnd_a}\n{bnd_b}");

        // Deletions
        // Severus
        let vcf = "chr7\t143674704\tseverus_DEL7318\tN\t<DEL>\t60\tPASS≥\tPRECISE;SVTYPE=DEL;SVLEN=3642;END=143678346;STRANDS=+-;MAPQ=60\tGT:GQ:VAF:hVAF:DR:DV\t0/1:114:0.39:0.39,0,0:14:9";
        let variant: VcfVariant = vcf.parse()?;

        println!("{:?}", variant.infos);
        println!("{:?}", variant.formats);
        let del = variant.deletion_desc().unwrap();
        println!("{:?}", del);
        println!("{:?} {:?}", del.len(), variant.formats.n_alt_depth());

        assert_eq!(
            "chr7:143674704_143678346_del",
            variant.deletion_desc().unwrap().to_string()
        );

        println!("--\n");
        let vcf="chr7\t144003249\tr_106\tC\t<DEL>\t.\tPASS\tEND=144142733;SVTYPE=DEL;SVLEN=-139484;SVINSLEN=4;SVINSSEQ=GCCA\tTR:VR\t12:10\t51:0";
        let variant: VcfVariant = vcf.parse()?;

        println!("{:?}", variant.infos);
        println!("{:?}", variant.formats);
        let del = variant.deletion_desc().unwrap();
        println!("{:?}", del);
        println!("{:?} {:?}", del.len(), variant.formats.n_alt_depth());

        let path = "/data/ref/hs1/chm13v2.0_RefSeq_Liftoff_v5.1_Genes.bed";
        let r = read_bed(path)?;

        let deleted_genes = bedrow_overlaps_par(
            &r,
            &[&GenomeRange {
                contig: variant.position.contig,
                range: del.start..del.end,
            }],
        )
        .into_iter()
        .filter_map(|e| e.name)
        .collect::<Vec<String>>()
        .join(", ");
        println!("{deleted_genes}");

        Ok(())
    }

    #[test]
    fn parse_del_vep() -> anyhow::Result<()> {
        init();
        let config = Config::default();
        let vcf = "chr14\t48907506\t.\tCCAG\tC\t25.40\tPASS\t.\tGT:GQ:DP:AD:VAF:MID:PL\t0/1:25:58:34,24:0.41379:small_model:25,0,33";
        let variant: VcfVariant = vcf.parse()?;
        let annotations = Annotations::default();
        let ext_annot = ExternalAnnotation::init("test", &config)?;
        ext_annot.annotate_vep(std::slice::from_ref(&variant), &annotations)?;
        { 
            let annots = annotations.store
                .get(&variant.hash) // adapt if your API differs
                .expect("variant should have annotations");

            let vep = annots
                .iter()
                .find_map(|a| match a {
                    Annotation::VEP(v) => Some(v),
                    _ => None,
                })
                .expect("variant should have VEP annotation");

            let samd4a_nm = vep
                .iter()
                .find(|v| v.feature.as_deref() == Some("NM_001161576.2"))
                .expect("NM_001161576.2 transcript should be annotated");

            assert_eq!(samd4a_nm.gene.as_deref(), Some("SAMD4A"));
            assert_eq!(samd4a_nm.feature_type.as_deref(), Some("Transcript"));

            assert!(samd4a_nm
                .consequence
                .as_ref()
                .is_some_and(|c| c.contains(&annotation::vep::VepConsequence::InframeDeletion)));

            assert_eq!(samd4a_nm.cds_position.as_deref(), Some("468-470"));
            assert_eq!(samd4a_nm.protein_position.as_deref(), Some("156-157"));
            assert_eq!(samd4a_nm.amino_acids.as_deref(), Some("TS/T"));
            assert_eq!(samd4a_nm.codons.as_deref(), Some("acCAGc/acc"));

            assert_eq!(samd4a_nm.extra.impact, Some(annotation::vep::VepImpact::MODERATE));
            assert_eq!(samd4a_nm.extra.symbol.as_deref(), Some("SAMD4A"));
            assert_eq!(
                samd4a_nm.extra.hgvs_c.as_deref(),
                Some("NM_001161576.2:c.472_474del")
            );
            assert_eq!(
                samd4a_nm.extra.hgvs_p.as_deref(),
                Some("NP_001155048.2:p.Ser158del")
            );

            info!("VEP ok");
        }

        annotate_colorsdb(std::slice::from_ref(&variant), &annotations, &config.colors_db_smallvar)?;
        let annots = annotations.store
            .get(&variant.hash) // adapt if your API differs
            .expect("variant should have annotations");

        let colors = annots
            .iter()
            .find_map(|a| match a {
                Annotation::CoLoRSdb(v) => Some(v),
                _ => None,
            })
            .expect("variant should have CoLoRSdb annotation");

        assert_eq!(colors.ac, 2634);

        Ok(())
    }

    #[test]
    // fn variant_load_deepvariant() -> anyhow::Result<()> {
    //     init();
    //     let id = "ADJAGBA";
    //     let time = "diag";
    //     let dv = DeepVariant::initialize(id, time, Config::default())?;
    //     let annotations = Annotations::default();
    //     let variant_collection = dv.variants(&annotations)?;
    //     println!(
    //         "Deepvariant for {id} {time}: variants {} {}",
    //         variant_collection.variants.len(),
    //         variant_collection.vcf.n_variants
    //     );
    //     Ok(())
    // }

    // #[test]
    // fn variant_load_clairs() -> anyhow::Result<()> {
    //     init();
    //     let id = "ELKLIFI";
    //     let clairs = NanomonSV::initialize(id, Config::default())?;
    //     let u = clairs.should_run();
    //     info!("should_run {u:?}");
    //     // let annotations = Annotations::default();
    //     // let variant_collection = clairs.variants(&annotations)?;
    //     // println!("ClairS for {id}: variants {} {}", variant_collection.variants.len(), variant_collection.vcf.n_variants);
    //     Ok(())
    // }
    //
    // #[test]
    // fn variant_load_nanomonsv() -> anyhow::Result<()> {
    //     init();
    //     let id = "ADJAGBA";
    //     let nanomonsv = NanomonSV::initialize(id, Config::default())?;
    //     let annotations = Annotations::default();
    //     let variant_collection = nanomonsv.variants(&annotations)?;
    //     println!(
    //         "NanomonSV for {id}: variants {} {}",
    //         variant_collection.variants.len(),
    //         variant_collection.vcf.n_variants
    //     );
    //     println!("{:?}", variant_collection.variants.first());
    //     Ok(())
    // }

    // #[test]
    // fn variant_load_clairs_germline() -> anyhow::Result<()> {
    //     init();
    //     let id = "ADJAGBA";
    //     let clairs = ClairS::initialize(id, Config::default())?;
    //     let annotations = Annotations::default();
    //     let germline_variant_collection = clairs.germline(&annotations)?;
    //     println!("ClairS for {id}: variants {} {}", germline_variant_collection.variants.len(), germline_variant_collection.vcf.n_variants);
    //     Ok(())
    // }
    #[test]
    fn overlaps() {
        init();
        let positions = vec![
            &GenomePosition {
                contig: 1,
                position: 100,
            },
            &GenomePosition {
                contig: 1,
                position: 150,
            },
            &GenomePosition {
                contig: 1,
                position: 200,
            },
            &GenomePosition {
                contig: 2,
                position: 150,
            },
        ];

        let ranges = vec![
            &GenomeRange {
                contig: 1,
                range: 50..150,
            },
            &GenomeRange {
                contig: 2,
                range: 100..200,
            },
        ];

        let parallel_overlapping_indices = overlaps_par(&positions, &ranges);
        assert_eq!(parallel_overlapping_indices, vec![0, 3])
    }

    #[test]
    fn bed_read() -> anyhow::Result<()> {
        init();
        let path = &Config::default().mask_bed("ADJAGBA");
        let r = read_bed(path)?;
        println!("{}", r.len());
        Ok(())
    }

    #[test]
    fn test_read_dict() -> anyhow::Result<()> {
        init();
        let genome = read_dict(&Config::default().dict_file)?;
        let genome_length: usize = genome.into_iter().map(|(_, len)| len as usize).sum();
        println!("{genome_length}");
        Ok(())
    }

    #[test]
    fn bases_at() -> anyhow::Result<()> {
        init();
        let id = "ADJAGBA";
        let c = Config::default();
        let chr = "chr3";
        let position = 62416039; // 1-based
        let mut bam = rust_htslib::bam::IndexedReader::from_path(c.solo_bam(id, "diag"))?;
        let p = nt_pileup(&mut bam, chr, position - 1, false)?
            .iter()
            .map(|e| String::from_utf8(vec![*e]).unwrap())
            .collect::<Vec<_>>();
        let mut counts = HashMap::new();

        for item in p.iter() {
            *counts.entry(item.as_str()).or_insert(0) += 1;
        }

        for (key, value) in &counts {
            println!("{}: {}", key, value);
        }

        assert_eq!(8, *counts.get("C").unwrap());
        assert_eq!(13, *counts.get("G").unwrap());
        assert_eq!(6, *counts.get("D").unwrap());

        let chr = "chr1";
        let position = 3220; // 1-based

        let mut bam = rust_htslib::bam::IndexedReader::from_path(c.solo_bam(id, "mrd"))?;
        let p = counts_at(&mut bam, chr, position - 1)?;

        println!("{p:#?}");
        Ok(())
    }

    #[test]
    fn seq_at() -> anyhow::Result<()> {
        init();
        let c = Config::default();
        let chr = "chr1";
        let position = 16761;

        let mut fasta_reader =
            noodles_fasta::io::indexed_reader::Builder::default().build_from_path(c.reference)?;
        let r = io::fasta::sequence_at(&mut fasta_reader, chr, position, 3)?;
        println!(
            "{r} ({} {:.2})",
            r.len(),
            estimate_shannon_entropy(r.as_str())
        );
        Ok(())
    }

    #[test]
    fn ins_at() -> anyhow::Result<()> {
        init();
        let id = "PASSARD";
        let c = Config::default();
        let chr = "chr5";
        let position = 36122736; // 1-based like in vcf
        let mut bam = rust_htslib::bam::IndexedReader::from_path(c.solo_bam(id, "mrd"))?;
        // let p = ins_pileup(&mut bam, chr, position - 1, true)?.iter().map(|e| String::from_utf8(vec![*e]).unwrap()).collect::<Vec<_>>();
        let counts = counts_ins_at(&mut bam, chr, position - 1)?;

        println!("{counts:?}");
        for (key, value) in &counts {
            println!("{}: {}", key, value);
        }

        Ok(())
    }

    // #[test]
    // fn del_at() -> anyhow::Result<()> {
    //     let id = "PASSARD";
    //             let c = Config::default();
    //
    //     let mut bam = rust_htslib::bam::IndexedReader::from_path(c.solo_bam(id, "mrd"))?;
    //
    //     let pileup_start = crate::collection::bam::nt_pileup(
    //                                     &mut bam,
    //                                     "chr5",
    //                                     36122735,
    //                                     false,
    //                                 )?;
    //
    //     let uu = crate::collection::bam::nt_pileup_new(
    //                                     &mut bam,
    //                                     "chr5",
    //                                     36122735,
    //                                     false,
    //                                 )?;
    //
    //                                 let pileup_end = crate::collection::bam::nt_pileup_new(
    //                                     &mut bam,
    //                                     &var.position.contig(),
    //                                     del_repr.end.saturating_sub(1),
    //                                     false,
    //                                 )?;
    //
    //     println!("{pileup_start:?}");
    //     println!("{uu:?}");
    //     let depth = uu.len().max(pileup_end.len());
    //
    //     Ok(())
    //
    // }

    #[test]
    fn vep_line() -> anyhow::Result<()> {
        init();

        let line = "chr2_1922358_-/T\tchr2:1922357-1922358\tT\tMYT1L\tNM_001303052.2\tTranscript\tintron_variant\t-\t-\t-\t-\t-\t-\tIMPACT=MODIFIER;STRAND=-1;SYMBOL=MYT1L;SOURCE=chm13v2.0_RefSeq_Liftoff_v5.1_sorted.gff3.gz;HGVSc=NM_001303052.2:c.1619-613dup";
        let vep_line: VepLine = line.parse()?;
        println!("{vep_line:#?}");

        let vep: VEP = VEP::try_from(&vep_line)?;
        println!("{vep:#?}");

        Ok(())
    }

    // #[test]
    // fn savana_cn() -> anyhow::Result<()> {
    //     init();
    //     let id = "CAMARA";
    //     // let s = SavanaCopyNumber::load_id(id, Config::default())?;
    //     let s = SavanaReadCounts::load_id(id, Config::default())?;
    //     println!("tumoral reads: {}", s.n_tumoral_reads());
    //     println!("normal reads: {}", s.n_normal_reads());
    //     println!("tumoral:\n{:#?}", s.norm_chr_counts());
    //     Ok(())
    // }

    #[test]
    fn load_bam() -> anyhow::Result<()> {
        init();
        let id = "ADJAGBA";
        let time = "diag";
        let bam_path = Config::default().solo_bam(id, time);
        WGSBam::new(Path::new(&bam_path).to_path_buf())?;
        Ok(())
    }

    #[test]
    fn tar() -> anyhow::Result<()> {
        init();
        scan_archive("/data/lto/20241030_CUNVI-DG-N20_SEBZA-DG-N21.tar")?;
        Ok(())
    }

    // #[test]
    // fn somatic_cases() -> anyhow::Result<()> {
    //     init();
    //     let id = "PASSARD";
    //     let config = Config {
    //         somatic_pipe_force: true,
    //         ..Default::default()
    //     };
    //     match SomaticPipe::initialize(id, config)?.run() {
    //         Ok(_) => (),
    //         Err(e) => error!("{id} {e}"),
    //     };
    //     Ok(())
    // }

    #[test]
    fn load_variants() -> anyhow::Result<()> {
        init();
        let id = "ACHITE";
        let config = Config::default();
        let path = format!("{}/{id}/diag/{id}_somatic_variants.bit", config.result_dir);
        let variants = variant_collection::Variants::load_from_file(&path)?;
        println!("n variants {}", variants.data.len());

        let n_vep: usize = variants.data.iter().map(|v| v.vep().len()).sum();
        println!("VEP: {n_vep}");

        let translocations = variants.get_alteration_cat(AlterationCategory::TRL);
        println!("{} translocations", translocations.len());
        let threshold = 5;
        let res = group_variants_by_bnd_desc(&translocations, 5);

        let rres = group_variants_by_bnd_rc(&res, threshold);

        rres.iter().for_each(|group| {
            println!("{} {}", group.0.len(), group.1.len());
        });

        Ok(())
    }

    #[test]
    fn load_fc() -> anyhow::Result<()> {
        init();
        // FlowCells::load_archive_from_scan("/data/lto", "/data/archives.json")?;
        let r = FlowCells::load(
            "/home/prom/mnt/store",
            "/data/pandora_id_inputs.json",
            "/data/archives.json.gz",
        )?;
        println!("{r:#?}");
        Ok(())
    }

    #[test]
    fn alt_cat() -> anyhow::Result<()> {
        let id = "ADJAGBA";
        let config = Config::default();
        let path = format!("{}/{id}/diag/somatic_variants.json.gz", config.result_dir);
        let variants = variant_collection::Variants::load_from_json(&path)?;
        println!("n variants {}", variants.data.len());
        variants
            .data
            .iter()
            .filter(|v| v.alteration_category().contains(&AlterationCategory::TRL))
            .for_each(|v| {
                println!(
                    "{:?} {}",
                    v.vcf_variants
                        .iter()
                        .map(|v| v.bnd_desc())
                        .collect::<Vec<_>>(),
                    v.annotations
                        .iter()
                        .filter(|a| matches!(a, Annotation::Callers(..)))
                        .map(|a| a.to_string())
                        .collect::<Vec<String>>()
                        .join(";")
                )
            });
        Ok(())
    }

    #[test]
    fn variants_stats() -> anyhow::Result<()> {
        init();
        let config = Config::default();

        let all_variants_bit = find_files(&format!(
            "{}/*/diag/*_somatic_variants.bit",
            config.result_dir
        ))?;
        for v in all_variants_bit.into_iter() {
            let id = v
                .file_name()
                .unwrap()
                .to_str()
                .unwrap()
                .split("_somatic")
                .next()
                .unwrap();
            println!("{id}");

            let config = Config::default();
            let path = format!("{}/{id}/diag/{id}_somatic_variants.bit", config.result_dir);
            match variant_collection::Variants::load_from_file(&path) {
                Ok(mut variants) => {
                    let (mut high_depth_ranges, _) = somatic_depth_quality_ranges(id, &config)?;
                    high_depth_ranges.par_sort_by_key(|r| (r.contig, r.range.start));

                    let res = VariantsStats::new(&mut variants, id, &config, &high_depth_ranges)?
                        .save_to_json(&format!(
                            "{}/{id}/diag/{id}_somatic_variants_stats.json.gz",
                            config.result_dir
                        ));
                    if res.is_err() {
                        info!("{:#?}", res);
                    }
                }
                Err(e) => error!("{e}"),
            }
        }
        Ok(())
    }

    // #[test]
    // fn constit_stats() {
    //     init();
    //     let id = "ADJAGBA";
    //     let config = Config::default();
    //
    //     let _ = const_stats(id.to_string(), config);
    // }

    // #[test]
    // fn test_bnd() -> anyhow::Result<()> {
    //     init();
    //     let id = "COIFFET";
    //     let config = Config::default();
    //
    //     let annotations = Annotations::default();
    //     let s = Savana::initialize(id, config)?.variants(&annotations)?;
    //     s.variants.iter().for_each(|e| {
    //         if let Ok(bnd) = e.bnd_desc() {
    //             println!("{}\t{}\t{}", e.position, e.reference, e.alternative);
    //             println!("{:#?}", bnd);
    //         }
    //     });
    //     Ok(())
    // }

    // #[test]
    // fn parse_savana_seg() {
    //     init();
    //     let r = SavanaCN::parse_file("ADJAGBA", &config::Config::default())
    //         .unwrap()
    //         .segments;
    //     println!("{} lines", r.len());
    //     println!("{:#?}", r.first().unwrap());
    // }

    #[test]
    fn whole_scan() -> anyhow::Result<()> {
        init();
        let id = "CHENU";
        let mut config = Config::default();
        let u = config.solo_bam(id, "mrd");
        println!("{u}");

        config.somatic_scan_force = true;
        somatic_scan(id, &config)?;
        Ok(())
    }

    #[test]
    fn parse_gff() -> anyhow::Result<()> {
        init();
        let id = "ADJAGBA";
        let config = Config::default();
        let path = format!("{}/{id}/diag/somatic_variants.bit", config.result_dir);
        let exon_ranges = features_ranges("exon", &config)?;
        let exon_ranges = merge_overlapping_genome_ranges(&exon_ranges);

        let variants = variant_collection::Variants::load_from_file(&path)?;
        let full = variants_stats::somatic_rates(&variants.data, &exon_ranges, &config);
        info!("{full:#?}");

        // let restrained: Vec<variant_collection::Variant> = variants.data.iter().filter(|v| v.vcf_variants.len() >= 2)
        //     .cloned().collect();
        // let min_2 = variants_stats::somatic_rates(&restrained, &exon_ranges, &config);
        // info!("{min_2:#?}");
        //
        // let restrained: Vec<variant_collection::Variant> = restrained.iter().filter(|v| v.vcf_variants.len() >= 3)
        //     .cloned().collect();
        // let min_3 = variants_stats::somatic_rates(&restrained, &exon_ranges, &config);
        // info!("{min_3:#?}");

        // let mut high_depth_ranges = variants_stats::high_depth_somatic(id, &config)?;
        // high_depth_ranges.par_sort_by_key(|r| ( r.contig, r.range.start ));
        //
        // let exon_ranges_ref: Vec<&GenomeRange> = exon_ranges.iter().collect();
        // let exons_high_depth = range_intersection_par(&high_depth_ranges.iter().collect::<Vec<&GenomeRange>>(), &exon_ranges_ref);
        //
        // let full = variants_stats::somatic_rates(&variants.data, &exons_high_depth, &config);
        // info!("{full:#?}");
        //

        // info!("n variants loaded: {}", variants.data.len());
        //
        // let r = features_ranges("exon", &config::Config::default())?;
        // info!("n exon: {}", r.len());
        //
        // let merged = merge_overlapping_genome_ranges(&r);
        // info!("n merged exon: {}", merged.len());
        //
        // let ol = par_overlaps(&variants.data, &r);
        // info!("variants in exon {}", ol.len());
        //
        // let n_coding = ol.iter().filter_map(|i| variants.data[*i].best_vep().ok() ).filter_map(|bv| bv.impact()).filter(|impact| *impact <= VepImpact::MODERATE).count();
        // info!("coding variants {n_coding}");
        //
        // let n_bases_m = merged.par_iter().map(|gr| gr.length()).sum::<u32>();
        // info!("{n_bases_m}nt");
        //
        // let mega_base_m = n_bases_m as f64 / 10.0e6;
        //
        // let wgs_len = read_dict(&config.dict_file)?.iter().map(|(_, l)| *l).sum::<u32>();
        // info!("wgs len {wgs_len}");
        // let rate_wgs = variants.data.len() as f64 / (wgs_len as f64 / 10.0e6);
        // info!("somatic mutation rate {rate_wgs:.2}/mb");
        //
        // let n_exons_mb = ol.len() as f64 / mega_base_m;
        // info!("somatic mutation rate in the coding region {n_exons_mb:.2}/mb");
        //
        // let n_exons_mb = n_coding as f64 / mega_base_m;
        // info!("somatic non synonymous mutation rate in the coding region {n_exons_mb:.2}/mb");

        Ok(())
    }

    fn gr(contig: u8, start: u32, end: u32) -> GenomeRange {
        GenomeRange {
            contig,
            range: start..end,
        }
    }

    #[test]
    fn test_both_empty() {
        let a: Vec<&GenomeRange> = vec![];
        let b: Vec<&GenomeRange> = vec![];
        let result = range_intersection_par(&a, &b);
        assert!(result.is_empty());
    }

    #[test]
    fn test_one_empty() {
        let a = [gr(1, 0, 100)];
        let b: Vec<&GenomeRange> = vec![];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b);
        sort_ranges(&mut result);
        assert!(result.is_empty());
    }

    #[test]
    fn test_single_range_no_overlap() {
        let a = [gr(1, 0, 100)];
        let b = [gr(1, 100, 200)];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let b_refs: Vec<&GenomeRange> = b.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b_refs);
        sort_ranges(&mut result);
        assert!(result.is_empty());
    }

    #[test]
    fn test_single_range_full_overlap() {
        let a = [gr(1, 0, 100)];
        let b = [gr(1, 0, 100)];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let b_refs: Vec<&GenomeRange> = b.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b_refs);
        let mut expected = [gr(1, 0, 100)];
        sort_ranges(&mut result);
        sort_ranges(&mut expected);
        assert_eq!(result, expected);
    }

    #[test]
    fn test_different_contigs() {
        let a = [gr(1, 0, 100)];
        let b = [gr(2, 0, 100)];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let b_refs: Vec<&GenomeRange> = b.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b_refs);
        sort_ranges(&mut result);
        assert!(result.is_empty());
    }

    #[test]
    fn test_touching_ranges() {
        let a = [gr(1, 0, 100)];
        let b = [gr(1, 100, 200)];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let b_refs: Vec<&GenomeRange> = b.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b_refs);
        sort_ranges(&mut result);
        assert!(result.is_empty());
    }

    #[test]
    fn test_complete_subrange() {
        let a = [gr(1, 0, 200)];
        let b = [gr(1, 50, 150)];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let b_refs: Vec<&GenomeRange> = b.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b_refs);
        let mut expected = [gr(1, 50, 150)];
        sort_ranges(&mut result);
        sort_ranges(&mut expected);
        assert_eq!(result, expected);
    }

    #[test]
    fn test_multiple_overlaps_same_contig() {
        let a = [gr(1, 0, 50), gr(1, 75, 125)];
        let b = [gr(1, 25, 100), gr(1, 150, 200)];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let b_refs: Vec<&GenomeRange> = b.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b_refs);
        let mut expected = [gr(1, 25, 50), gr(1, 75, 100)];
        sort_ranges(&mut result);
        sort_ranges(&mut expected);
        assert_eq!(result, expected);
    }

    #[test]
    fn test_multiple_contigs() {
        let a = [gr(1, 0, 100), gr(2, 50, 150), gr(3, 200, 300)];
        let b = [gr(1, 50, 150), gr(2, 0, 100), gr(4, 0, 100)];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let b_refs: Vec<&GenomeRange> = b.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b_refs);
        let mut expected = [gr(1, 50, 100), gr(2, 50, 100)];
        sort_ranges(&mut result);
        sort_ranges(&mut expected);
        assert_eq!(result, expected);
    }

    #[test]
    fn test_adjacent_ranges() {
        let a = [gr(1, 0, 50), gr(1, 50, 100)];
        let b = [gr(1, 25, 75)];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let b_refs: Vec<&GenomeRange> = b.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b_refs);
        let mut expected = [gr(1, 25, 50), gr(1, 50, 75)];
        sort_ranges(&mut result);
        sort_ranges(&mut expected);
        assert_eq!(result, expected);
    }

    #[test]
    fn test_minimal_overlap() {
        let a = [gr(1, 0, 100)];
        let b = [gr(1, 99, 200)];
        let a_refs: Vec<&GenomeRange> = a.iter().collect();
        let b_refs: Vec<&GenomeRange> = b.iter().collect();
        let mut result = range_intersection_par(&a_refs, &b_refs);
        let mut expected = [gr(1, 99, 100)];
        sort_ranges(&mut result);
        sort_ranges(&mut expected);
        assert_eq!(result, expected);
    }

    /// helper to build a forward‑strand BND (same contig) where
    ///   A = pos,  B = pos + 5
    fn fwd(contig: &str, pos: u32) -> BNDDesc {
        BNDDesc {
            a_contig: contig.into(),
            a_position: pos,
            a_sens: true,
            b_contig: contig.into(),
            b_position: pos + 5,
            b_sens: true,
            added_nt: String::new(),
        }
    }

    /// Build a six‑node *forward* chain relying **only** on `auto_connect()`
    /// (no manual edges) and assert the Hamiltonian path spans all nodes.
    #[test]
    fn hamiltonian_chain_auto() {
        // positions 10,15   20,25   …  60,65  satisfy B(u) ≤ A(v)
        let bnds: Vec<BNDDesc> = (1..=6).map(|i| fwd("chr1", i * 10)).collect();
        let g: BNDGraph<()> = bnds.to_bnd_graph(); // trait uses auto_connect()
                                                   // ensure auto_connect produced 5 edges in a line
        assert_eq!(g.inner().edge_count(), 5);
        let path = g.hamiltonian_path().expect("chain should be Hamiltonian");
        assert_eq!(path.len(), 6);
    }

    /// Two disconnected "V" shapes -> auto_connect creates 2 edges in each
    /// component, no Hamiltonian path, components sorted by size.
    #[test]
    fn components_after_auto_connect() {
        // comp1: a->b<-c on reverse strand of chrX
        let a = BNDDesc {
            a_contig: "chrX".into(),
            a_position: 300,
            a_sens: false,
            b_contig: "chrX".into(),
            b_position: 250,
            b_sens: false,
            added_nt: String::new(),
        };
        let b = BNDDesc {
            a_contig: "chrX".into(),
            a_position: 200,
            a_sens: false,
            b_contig: "chrX".into(),
            b_position: 150,
            b_sens: false,
            added_nt: String::new(),
        };
        let c = BNDDesc {
            a_contig: "chrX".into(),
            a_position: 100,
            a_sens: false,
            b_contig: "chrX".into(),
            b_position: 50,
            b_sens: false,
            added_nt: String::new(),
        };
        // comp2: three‑node forward chain on chrY
        let d = fwd("chrY", 10);
        let e = fwd("chrY", 20);
        let f = fwd("chrY", 30);
        let g: BNDGraph<()> = vec![a, b, c, d, e, f].to_bnd_graph();
        assert!(g.hamiltonian_path().is_none());
        let comps = g.components_by_size();
        comps.iter().for_each(|a| println!("{}", g.fmt_path(a)));
        assert_eq!(comps.len(), 2);
        assert_eq!(comps[0].len(), 3);
        assert_eq!(comps[1].len(), 3);
    }

    #[test]
    fn bnd_connections() {
        // comp1: a->b<-c on reverse strand of chrX
        let a = BNDDesc {
            a_contig: "chrX".into(),
            a_position: 300,
            a_sens: true,
            b_contig: "chr14".into(),
            b_position: 250,
            b_sens: false,
            added_nt: String::new(),
        };
        let b = BNDDesc {
            a_contig: "chr14".into(),
            a_position: 200,
            a_sens: false,
            b_contig: "chrX".into(),
            b_position: 301,
            b_sens: true,
            added_nt: String::new(),
        };
        let c = BNDDesc {
            a_contig: "chrX".into(),
            a_position: 302,
            a_sens: true,
            b_contig: "chrZ".into(),
            b_position: 50,
            b_sens: false,
            added_nt: String::new(),
        };
        // comp2: three‑node forward chain on chrY
        let d = fwd("chrY", 10);
        let e = fwd("chrY", 20);
        let f = fwd("chrY", 30);
        let g: BNDGraph<()> = vec![a, b, c, d, e, f].to_bnd_graph();
        assert!(g.hamiltonian_path().is_none());
        let comps = g.components_by_size();
        comps.iter().for_each(|a| println!("{}", g.fmt_path(a)));
        assert_eq!(comps.len(), 2);
        assert_eq!(comps[0].len(), 3);
        assert_eq!(comps[1].len(), 3);
    }

    #[test]
    fn snv_hg38() -> anyhow::Result<()> {
        init();
        let id = "CHALO";
        let config = Config::default();
        let path = format!("{}/{id}/diag/{id}_somatic_variants.bit", config.result_dir);

        let mut variants = variant_collection::Variants::load_from_file(&path)?;
        info!("All: {}", variants.len());
        variants.retain(|v| {
            *v.alteration_category()
                .first()
                .unwrap_or(&AlterationCategory::Other)
                == AlterationCategory::SNV
        });
        info!("SNV: {}", variants.len());
        variants.in_place_merge();
        info!("SNV: {}", variants.len());

        // variants.write_vcf("/home/t_steimle/CHALO_hs1.vcf.gz", "/home/t_steimle/ref/hs1/chm13v2.0.dict", true)
        variants.save_into_lifted(
            "/home/t_steimle/CHALO_SNV_hg38.vcf",
            "./unmapped",
            "/home/t_steimle/ref/hs1/chm13v2-grch38.chain",
            "/home/t_steimle/ref/hg38/hg38.fa",
            100,
        )
    }
}