colorsdb update

pandora-config.example.toml

@@ -54,6 +54,7 @@ db_snp = "/home/t_steimle/ref/hs1/chm13v2.0_dbSNPv155.vcf.gz"
 
 # CoLoRSdb path (should be indexed, tbi)
 colors_db = "/home/t_steimle/ref/hs1/CoLoRSdb.CHM13.v1.2.0.pbsv.jasmine.vcf.gz"
+colors_db_smallvar = "/home/t_steimle/ref/hs1/CoLoRSdb.CHM13.v1.0.0.deepvariant.glnexus.vcf.gz"
 
 # BED with genes on the 4th column should be sorted
 genes_bed = "/home/t_steimle/ref/hs1/chm13v2.0_RefSeq_Liftoff_v5.1_Genes.bed"

src/annotation/colorscb.rs

@@ -41,7 +41,7 @@ impl TryFrom<&Infos> for ColorsDbEntry {
         let mut exchet = None;
         let mut nhomalt = None;
 
-        for info in &infos.0 {
+        for info in infos.0.iter() {
             match info {
                 Info::AC(v) => ac = Some(*v),
                 Info::AN(v) => an = Some(*v),

src/config.rs

@@ -72,6 +72,7 @@ pub struct Config {
 
     /// CoLoRSdb jasmine (should be indexed, tbi)
     pub colors_db: String,
+    pub colors_db_smallvar: String,
 
     /// BED with genes on the 4th column
     pub genes_bed: String,

src/helpers.rs

@@ -1074,8 +1074,13 @@ impl Drop for TempFileGuard {
     fn drop(&mut self) {
         if !self.disarmed.load(std::sync::atomic::Ordering::SeqCst) {
             warn!(
-                "Pipeline failed or panicked, cleaning up {} temporary files",
-                self.files.len()
+                "Cleaning up {} temporary files :{}",
+                self.files.len(),
+                self.files
+                    .iter()
+                    .map(|f| format!("\n\t- {}", f.display()))
+                    .collect::<Vec<String>>()
+                    .join("")
             );
             self.cleanup();
         }

src/lib.rs

@@ -182,10 +182,10 @@ mod tests {
     use self::{/* collection::pod5::{FlowCellCase, Pod5Collection}, */ config::Config};
     use super::*;
     use crate::{
-        annotation::Annotation,
+        annotation::{Annotation, colorscb::annotate_colorsdb},
         collection::{
             bam::{self},
-            flowcells::{scan_archive, FlowCells},
+            flowcells::{FlowCells, scan_archive},
             vcf::VcfCollection,
         },
         helpers::find_files,
@@ -194,9 +194,9 @@ mod tests {
         scan::scan::somatic_scan,
         variant::{
             variant_collection::{
-                group_variants_by_bnd_desc, group_variants_by_bnd_rc, VariantCollection,
+                ExternalAnnotation, VariantCollection, group_variants_by_bnd_desc, group_variants_by_bnd_rc
             },
-            variants_stats::{self, somatic_depth_quality_ranges, VariantsStats},
+            variants_stats::{self, VariantsStats, somatic_depth_quality_ranges},
             vcf_variant::{AlterationCategory, BNDDesc, BNDGraph, ToBNDGraph},
         },
     };
@@ -771,6 +771,78 @@ mod tests {
         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();

src/pipes/somatic.rs

@@ -60,7 +60,7 @@
 //! - Annotation matrix generation utilities
 use crate::{
     annotation::{
-        dbsnp::annotate_dbsnp, is_dbsnp_and_constit_alt, is_gnomad_and_constit_alt, Caller,
+        Caller, colorscb::annotate_colorsdb, dbsnp::annotate_dbsnp, is_dbsnp_and_constit_alt, is_gnomad_and_constit_alt
     },
     create_should_run_normal_tumoral, init_solo_callers_normal_tumoral,
     io::bed::read_bed,
@@ -69,7 +69,7 @@ use crate::{
     scan::scan::SomaticScan,
     variant::{
         variants_stats::somatic_depth_quality_ranges,
-        vcf_variant::{run_if_required, ShouldRunBox},
+        vcf_variant::{ShouldRunBox, run_if_required},
     },
 };
 use log::info;
@@ -633,7 +633,16 @@ impl Run for SomaticPipe {
                 )
             })
             .try_for_each(|c| -> anyhow::Result<()> {
-                annotate_dbsnp(&c.variants, &annotations, &config.db_snp)
+                if matches!(
+                    c.caller,
+                    Annotation::Callers(Caller::Savana, Sample::Somatic)
+                        | Annotation::Callers(Caller::Severus, Sample::Somatic)
+                        | Annotation::Callers(Caller::NanomonSV, Sample::Somatic)
+                ) {
+                    annotate_colorsdb(&c.variants, &annotations, &config.colors_db)
+                } else {
+                    annotate_colorsdb(&c.variants, &annotations, &config.colors_db_smallvar)
+                }
             })?;
 
         annotations

src/variant/vcf_variant.rs

@@ -974,16 +974,16 @@ fn recode_gt(gt: &str, allele_idx: usize) -> String {
         .join(&sep.to_string())
 }
 
-fn pick_nth(values: &str, idx: usize) -> &str {
-    values.split(',').nth(idx).unwrap_or(".")
-}
-
-fn pick_r(values: &str, allele_idx: usize) -> String {
-    let mut it = values.split(',');
-    let ref_val = it.next().unwrap_or(".");
-    let alt_val = it.nth(allele_idx).unwrap_or(".");
-    format!("{ref_val},{alt_val}")
-}
+// fn pick_nth(values: &str, idx: usize) -> &str {
+//     values.split(',').nth(idx).unwrap_or(".")
+// }
+//
+// fn pick_r(values: &str, allele_idx: usize) -> String {
+//     let mut it = values.split(',');
+//     let ref_val = it.next().unwrap_or(".");
+//     let alt_val = it.nth(allele_idx).unwrap_or(".");
+//     format!("{ref_val},{alt_val}")
+// }
 
 // ---------------------- VcfVariant SV graph
 
@@ -1574,15 +1574,15 @@ pub enum Info {
     RMSK_INFO(String),
     TSD(String),
     // CoLoRSdb
-    AC(u32), // Allele count in genotypes
-    AN(u32), // Total number of alleles in called genotypes
-    NS(u32), // Number of samples with data
-    AC_Hom(u32), // Allele counts in homozygous genotypes
-    AC_Het(u32), // Allele counts in heterozygous genotypes
+    AC(u32),      // Allele count in genotypes
+    AN(u32),      // Total number of alleles in called genotypes
+    NS(u32),      // Number of samples with data
+    AC_Hom(u32),  // Allele counts in homozygous genotypes
+    AC_Het(u32),  // Allele counts in heterozygous genotypes
     AC_Hemi(u32), // Allele counts in hemizygous genotypes
-    AF(f32), // Allele frequency
-    HWE(f32), // HWE test (PMID:15789306); 1=good, 0=bad
-    ExcHet(f32), // Test excess heterozygosity; 1=good, 0=bad
+    AF(f32),      // Allele frequency
+    HWE(f32),     // HWE test (PMID:15789306); 1=good, 0=bad
+    ExcHet(f32),  // Test excess heterozygosity; 1=good, 0=bad
     Nhomalt(u32), // The number of individuals that are called homozygous for the alternate allele.
 }
 
@@ -1674,6 +1674,17 @@ impl FromStr for Info {
                 "RMSK_INFO" => Info::RMSK_INFO(value.to_string()),
                 "TSD" => Info::TSD(value.to_string()),
 
+                "AC" => Info::AC(parse_value(value, key)?),
+                "AN" => Info::AN(parse_value(value, key)?),
+                "NS" => Info::NS(parse_value(value, key)?),
+                "AC_Hom" => Info::AC_Hom(parse_value(value, key)?),
+                "AC_Het" => Info::AC_Het(parse_value(value, key)?),
+                "AC_Hemi" => Info::AC_Hemi(parse_value(value, key)?),
+                "AF" => Info::AF(parse_value(value, key)?),
+                "HWE" => Info::HWE(parse_value(value, key)?),
+                "ExcHet" => Info::ExcHet(parse_value(value, key)?),
+                "nhomalt" => Info::Nhomalt(parse_value(value, key)?),
+
                 _ => Info::Empty,
             })
         } else {