colorsdb annotation

pandora-config.example.toml

@@ -52,6 +52,9 @@ refseq_gtf = "/home/t_steimle/ref/hs1/chm13v2.0_RefSeq_Liftoff_v5.1_sorted.gtf"
 # dbSNP vcf.gz file (should be indexed)
 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"
+
 # 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

@@ -0,0 +1,157 @@
+use std::sync::atomic::{AtomicUsize, Ordering};
+
+use bitcode::{Decode, Encode};
+use log::{error, info, warn};
+use rayon::prelude::*;
+use serde::{Deserialize, Serialize};
+
+use crate::{
+    annotation::{Annotation, Annotations},
+    io::{readers::TabixReader, vcf::fetch_vcf_with_reader},
+    positions::GenomeRange,
+    variant::vcf_variant::{Info, Infos, VcfVariant},
+};
+
+#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode)]
+pub struct ColorsDbEntry {
+    pub ac: u32,
+    pub an: u32,
+    pub ns: u32,
+    pub ac_hom: u32,
+    pub ac_het: u32,
+    pub ac_hemi: u32,
+    pub af: f32,
+    pub hwe: f32,
+    pub exchet: f32,
+    pub nhomalt: u32,
+}
+
+impl TryFrom<&Infos> for ColorsDbEntry {
+    type Error = String;
+
+    fn try_from(infos: &Infos) -> Result<Self, Self::Error> {
+        let mut ac = None;
+        let mut an = None;
+        let mut ns = None;
+        let mut ac_hom = None;
+        let mut ac_het = None;
+        let mut ac_hemi = None;
+        let mut af = None;
+        let mut hwe = None;
+        let mut exchet = None;
+        let mut nhomalt = None;
+
+        for info in &infos.0 {
+            match info {
+                Info::AC(v) => ac = Some(*v),
+                Info::AN(v) => an = Some(*v),
+                Info::NS(v) => ns = Some(*v),
+                Info::AC_Hom(v) => ac_hom = Some(*v),
+                Info::AC_Het(v) => ac_het = Some(*v),
+                Info::AC_Hemi(v) => ac_hemi = Some(*v),
+                Info::AF(v) => af = Some(*v),
+                Info::HWE(v) => hwe = Some(*v),
+                Info::ExcHet(v) => exchet = Some(*v),
+                Info::Nhomalt(v) => nhomalt = Some(*v),
+                _ => {}
+            }
+        }
+
+        Ok(ColorsDbEntry {
+            ac: ac.ok_or("missing INFO/AC")?,
+            an: an.ok_or("missing INFO/AN")?,
+            ns: ns.ok_or("missing INFO/NS")?,
+            ac_hom: ac_hom.ok_or("missing INFO/AC_Hom")?,
+            ac_het: ac_het.ok_or("missing INFO/AC_Het")?,
+            ac_hemi: ac_hemi.ok_or("missing INFO/AC_Hemi")?,
+            af: af.ok_or("missing INFO/AF")?,
+            hwe: hwe.ok_or("missing INFO/HWE")?,
+            exchet: exchet.ok_or("missing INFO/ExcHet")?,
+            nhomalt: nhomalt.ok_or("missing INFO/Nhomalt")?,
+        })
+    }
+}
+
+pub fn annotate_colorsdb(
+    variants: &[VcfVariant],
+    annotations: &Annotations,
+    colorsdb_path: &str,
+) -> anyhow::Result<()> {
+    if variants.is_empty() {
+        return Ok(());
+    }
+
+    let n_threads = rayon::current_num_threads();
+    let chunk_size = variants.len().div_ceil(n_threads).max(500);
+
+    let annotated_count = AtomicUsize::new(0);
+
+    variants.par_chunks(chunk_size).for_each_init(
+        || TabixReader::open(colorsdb_path),
+        |reader_res, chunk| {
+            let Ok(ref mut reader) = reader_res else {
+                error!("Failed to open CoLoRSdb tabix reader for chunk");
+                return;
+            };
+
+            for variant in chunk {
+                if let Some(anns) = annotations.store.get(&variant.hash) {
+                    if anns
+                        .iter()
+                        .any(|a| matches!(a, Annotation::GnomAD(_) | Annotation::DbSnp(_, _)))
+                    {
+                        continue;
+                    }
+                }
+
+                let region = GenomeRange::new(
+                    &variant.position.contig(),
+                    variant.position.position,
+                    variant.position.position + 1,
+                );
+
+                let colorsdb_records = match fetch_vcf_with_reader(reader, &region) {
+                    Ok(v) => v,
+                    Err(e) => {
+                        warn!(
+                            "CoLoRSdb fetch failed at {}:{}: {e}",
+                            variant.position.contig(),
+                            variant.position.position
+                        );
+                        continue;
+                    }
+                };
+
+                for dv in &colorsdb_records {
+                    if dv != variant {
+                        continue;
+                    }
+
+                    match (&dv.infos).try_into() {
+                        Ok(colors) => {
+                            annotations
+                                .insert_update(variant.hash, &[Annotation::CoLoRSdb(colors)]);
+                            annotated_count.fetch_add(1, Ordering::Relaxed);
+                        }
+                        Err(e) => {
+                            warn!(
+                                "CoLoRSdb parse failed at {}:{}: {e}",
+                                variant.position.contig(),
+                                variant.position.position
+                            )
+                        }
+                    }
+
+                    break;
+                }
+            }
+        },
+    );
+
+    info!(
+        "Annotated {} variants with CoLoRSdb",
+        annotated_count.load(Ordering::Relaxed)
+    );
+
+    Ok(())
+}

src/annotation/mod.rs

@@ -5,6 +5,7 @@ pub mod echtvar;
 pub mod gnomad;
 pub mod ncbi;
 pub mod vep;
+pub mod colorscb;
 
 use std::{
     collections::{HashMap, HashSet},
@@ -20,8 +21,7 @@ use std::{
 };
 
 use crate::{
-    helpers::{format_count, mean, Blake3BuildHasher, Hash128},
-    variant::{variant_collection::VariantCollection, vcf_variant::AlterationCategory},
+    annotation::colorscb::ColorsDbEntry, helpers::{Blake3BuildHasher, Hash128, format_count, mean}, variant::{variant_collection::VariantCollection, vcf_variant::AlterationCategory}
 };
 use bitcode::{Decode, Encode};
 use cosmic::Cosmic;
@@ -107,6 +107,9 @@ pub enum Annotation {
 
     /// query_start,query_end,class,class...
     RepeatMasker(String),
+
+    /// CoLoRSdb
+    CoLoRSdb(ColorsDbEntry),
 }
 
 /// Denotes the biological sample type associated with a variant call.
@@ -249,6 +252,7 @@ impl fmt::Display for Annotation {
             TSD(_) => "TSD".to_string(),
             InsertionType(v) => format!("InsertionType_{v}"),
             RepeatMasker(_) => "RepeatMasker".to_string(),
+            CoLoRSdb(_) => "CoLoRSdb".to_string(),
         };
         write!(f, "{}", s)
     }

src/config.rs

@@ -70,6 +70,9 @@ pub struct Config {
     ///  dbSNP vcf.gz file (should be indexed)
     pub db_snp: String,
 
+    /// CoLoRSdb jasmine (should be indexed, tbi)
+    pub colors_db: String,
+
     /// BED with genes on the 4th column
     pub genes_bed: String,
 

src/pipes/somatic.rs

@@ -621,6 +621,21 @@ impl Run for SomaticPipe {
                 Ok(())
             })?;
 
+        info!("Annotation of SV with CoLoRSdb.");
+        variants_collections
+            .iter()
+            .filter(|c| {
+                matches!(
+                    c.caller,
+                    Annotation::Callers(Caller::Savana, Sample::Somatic)
+                        | Annotation::Callers(Caller::Severus, Sample::Somatic)
+                        | Annotation::Callers(Caller::NanomonSV, Sample::Somatic)
+                )
+            })
+            .try_for_each(|c| -> anyhow::Result<()> {
+                annotate_dbsnp(&c.variants, &annotations, &config.db_snp)
+            })?;
+
         annotations
             .callers_stat(Some(Box::new(caller_cat_anns)))
             .save_to_json(&format!("{stats_dir}/{id}_annotations_09_vep.json"))?;

src/variant/vcf_variant.rs

@@ -1573,6 +1573,17 @@ pub enum Info {
     INSERT_TYPE(String),
     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_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
+    Nhomalt(u32), // The number of individuals that are called homozygous for the alternate allele.
 }
 
 impl FromStr for Info {
@@ -1683,6 +1694,7 @@ impl FromStr for Info {
 impl fmt::Display for Info {
     /// Converts the `Info` enum into a VCF-compliant string (key=value or flag).
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        use Info::*;
         match self {
             Info::Empty => write!(f, "."),
             Info::H => write!(f, "H"),
@@ -1769,6 +1781,18 @@ impl fmt::Display for Info {
             Info::RMSK_INFO(v) => write!(f, "RMSK_INFO={v}"),
             // Target site duplication sequence
             Info::TSD(v) => write!(f, "TSD={v}"),
+
+            // CoLoRSdb
+            AC(v) => write!(f, "AC={v}"),
+            AN(v) => write!(f, "AN={v}"),
+            NS(v) => write!(f, "NS={v}"),
+            AC_Hom(v) => write!(f, "AC_Hom={v}"),
+            AC_Het(v) => write!(f, "AC_Het={v}"),
+            AC_Hemi(v) => write!(f, "AC_Hemi={v}"),
+            AF(v) => write!(f, "AF={v}"),
+            HWE(v) => write!(f, "HWE={v}"),
+            ExcHet(v) => write!(f, "ExcHet={v}"),
+            Nhomalt(v) => write!(f, "nhomalt={v}"),
         }
     }
 }