m

src/io/fasta.rs

@@ -0,0 +1,26 @@
+use std::fs::File;
+
+// 0-based position
+pub fn sequence_at(
+    fasta_reader: &mut noodles_fasta::IndexedReader<noodles_fasta::io::BufReader<File>>,
+    contig: &str,
+    position: usize,
+    len: usize,
+) -> anyhow::Result<String> {
+    // convert to 1-based
+    let position = position + 1;
+
+    let start = position.saturating_sub(len / 2).max(1);
+    let end = start + len - 1;
+    // debug!("Region {contig}:{start}-{end} (1-based inclusive)");
+
+    let start = noodles_core::Position::try_from(start)?;
+    let end = noodles_core::Position::try_from(end)?;
+    let interval = noodles_core::region::interval::Interval::from(start..=end);
+
+    let r = noodles_core::Region::new(contig.to_string(), interval);
+    let record = fasta_reader.query(&r)?;
+    let s = String::from_utf8(record.sequence().as_ref().to_vec())?.to_uppercase();
+
+    Ok(s)
+}

src/io/mod.rs

@@ -3,3 +3,4 @@ pub mod readers;
 pub mod vcf;
 pub mod bed;
 pub mod dict;
+pub mod fasta;

src/pipes/somatic.rs

@@ -286,9 +286,9 @@ impl Run for Somatic {
             .map_err(|e| anyhow::anyhow!("Error while loading variants\n{e}"))?;
 
         info!("Loading Germline");
-        // let clairs_germline =
-        //     ClairS::initialize(&id, self.config.clone())?.germline(&annotations)?;
-        // variants_collections.push(clairs_germline);
+        let clairs_germline =
+            ClairS::initialize(&id, self.config.clone())?.germline(&annotations)?;
+        variants_collections.push(clairs_germline);
 
         let mut somatic_stats = SomaticStats::init(&variants_collections);
         info!(
@@ -496,16 +496,6 @@ impl Run for Somatic {
 
         annotations.vep_stats()?;
 
-        // Annotate echtvar
-        info!("GnomAD and Cosmic annotation.");
-        variants_collections
-            .iter()
-            .try_for_each(|c| -> anyhow::Result<()> {
-                let ext_annot = ExternalAnnotation::init()?;
-                ext_annot.annotate(&c.variants, &annotations)?;
-                Ok(())
-            })?;
-
         let variants = variants_collections.into_iter().fold(
             Variants::default(),
             |mut acc, variants_collection| {
@@ -525,7 +515,7 @@ impl Run for Somatic {
 }
 
 // 0-based position
-pub fn sequence_at(
+fn sequence_at(
     fasta_reader: &mut noodles_fasta::IndexedReader<noodles_fasta::io::BufReader<File>>,
     contig: &str,
     position: usize,

src/variant/variant.rs

@@ -259,6 +259,7 @@ impl VcfVariant {
         }
     }
 }
+
 #[derive(Debug, Clone, Serialize, Deserialize, Eq, PartialEq, Hash)]
 pub struct BNDDesc {
     pub a_contig: String,

src/variant/variant_collection.rs

@@ -27,8 +27,7 @@ use crate::{
         vcf::Vcf,
     },
     helpers::{app_storage_dir, estimate_shannon_entropy, temp_file_path, Hash128},
-    io::{readers::get_reader, vcf::vcf_header},
-    pipes::somatic::sequence_at,
+    io::{fasta::sequence_at, readers::get_reader, vcf::vcf_header},
     positions::GenomePosition,
 };