|
|
@@ -1,7 +1,5 @@
|
|
|
use crate::{
|
|
|
- create_should_run_normal_tumoral, init_solo_callers_normal_tumoral,
|
|
|
- scan::scan::SomaticScan,
|
|
|
- variant::variant::{run_if_required, ShouldRunBox},
|
|
|
+ annotation::is_gnomad_and_constit_alt, collection::ShouldRun, create_should_run_normal_tumoral, init_solo_callers_normal_tumoral, scan::scan::SomaticScan, variant::variant::{run_if_required, ShouldRunBox}
|
|
|
};
|
|
|
use anyhow::Context;
|
|
|
use itertools::Itertools;
|
|
|
@@ -11,7 +9,6 @@ use std::{
|
|
|
fs::{self, File},
|
|
|
io::Write,
|
|
|
path::Path,
|
|
|
- sync::Arc,
|
|
|
};
|
|
|
|
|
|
use crate::{
|
|
|
@@ -31,240 +28,172 @@ use crate::{
|
|
|
},
|
|
|
};
|
|
|
|
|
|
-pub struct Somatic {
|
|
|
+/// Runs the full somatic variant calling pipeline for a single sample (`id`).
|
|
|
+///
|
|
|
+/// This function orchestrates the entire somatic variant discovery process,
|
|
|
+/// starting from raw variant caller outputs (`PASSED VCF`) and applying multiple filtering
|
|
|
+/// and annotation steps to produce high-confidence somatic variants.
|
|
|
+///
|
|
|
+/// This function orchestrates the end-to-end somatic variant discovery process:
|
|
|
+/// - Executes and verifies upstream components if necessary
|
|
|
+/// - Loads variants from multiple callers (tumor and normal samples)
|
|
|
+/// - Applies several annotation and filtering steps (e.g. depth, population frequency, entropy)
|
|
|
+/// - Tracks filtering statistics at each step
|
|
|
+/// - Outputs high-confidence somatic variants in both `.json.gz` and `.bit` formats
|
|
|
+///
|
|
|
+/// ## Output Overview
|
|
|
+/// The final output includes:
|
|
|
+/// - `{tumoral_dir}/{id}_somatic_variants.json.gz`: annotated somatic variants
|
|
|
+/// - `{tumoral_dir}/{id}_somatic_variants.bit`: compact binary variant representation
|
|
|
+/// - `{stats_dir}/`: multiple intermediate JSON files with annotations and statistics
|
|
|
+///
|
|
|
+/// ## Steps
|
|
|
+///
|
|
|
+/// This pipeline performs the following high-level steps:
|
|
|
+///
|
|
|
+/// ### 1. Output Existence Check
|
|
|
+/// If the final JSON result already exists and [`Config::somatic_pipe_force`] is not set,
|
|
|
+/// the pipeline aborts early to avoid overwriting results.
|
|
|
+///
|
|
|
+/// ### 2. Initialization
|
|
|
+/// - Prepares statistics and output directories.
|
|
|
+/// - Initializes variant annotations.
|
|
|
+///
|
|
|
+/// ### 3. Pre-requisite Tool Execution
|
|
|
+/// Runs any required upstream components (e.g., alignment, basecalling, variant callers) if
|
|
|
+/// their outputs are missing, using the [`run_if_required`] logic.
|
|
|
+///
|
|
|
+/// ### 4. Load Variant Collections
|
|
|
+/// - Initializes the configured somatic variant callers and loads their output variants from PASSED VCF.
|
|
|
+/// - Also loads germline variants (from [`ClairS::germline`]) for comparative germline filtering.
|
|
|
+///
|
|
|
+/// ### 5. Statistics Initialization
|
|
|
+/// - Initializes a [`SomaticPipeStats`] object to track the number of variants at each step.
|
|
|
+/// - Captures initial variant counts before filtering.
|
|
|
+/// - Aggregates variant counts from all sources and stores initial annotations for quality control and
|
|
|
+/// comparison before filtering.
|
|
|
+///
|
|
|
+/// ### 6. Filter: Germline & Solo Constitutional Variants
|
|
|
+/// - Removes variants labeled as either Germline or Solo Constitutional, assuming they are unlikely to be somatic.
|
|
|
+/// - Records count of removed variants in [`SomaticPipeStats::n_constit_germline`].
|
|
|
+///
|
|
|
+/// ### 7. Annotation: BAM Read Depth and Alt Allele Counts
|
|
|
+/// - Uses the constitutional BAM file to annotate each variant with read depth and the number of alternate reads observed.
|
|
|
+/// - Flags variants with low depth or excessive alt reads for filtering as specified in [`Config`].
|
|
|
+///
|
|
|
+/// ### 8. Filtering: Low Depth / High Alt Alleles
|
|
|
+/// - Removes variants with low coverage in the constitutional sample, or excessive alt allele support (suggestive of germline origin).
|
|
|
+/// - Updates stats:
|
|
|
+/// - [`SomaticPipeStats::n_low_constit`]
|
|
|
+/// - [`SomaticPipeStats::n_high_alt_constit`]
|
|
|
+///
|
|
|
+/// ### 9. Annotation: Sequence Entropy
|
|
|
+/// Adds Shannon entropy annotation based on the reference sequence context
|
|
|
+/// around each variant (cf. [`Config::entropy_seq_len`]) to flag low-complexity regions (often repetitive).
|
|
|
+///
|
|
|
+/// ### 10. Annotation: External Databases (COSMIC, GnomAD):
|
|
|
+/// - Uses external resources to annotate variants with:
|
|
|
+/// - COSMIC hits (somatic mutation database)
|
|
|
+/// - GnomAD allele frequencies
|
|
|
+///
|
|
|
+/// ### 11. Filtering: GnomAD + Alt Support in Constitutional Sample
|
|
|
+/// - Removes variants that are both present in GnomAD **and** show
|
|
|
+/// alternate allele support in the constitutional BAM.
|
|
|
+/// These are highly likely to be non-somatic germline polymorphisms.
|
|
|
+/// - Updates [`SomaticPipeStats::n_high_alt_constit_gnomad`] stat.
|
|
|
+///
|
|
|
+/// ### 12. Filtering: Low Shannon Entropy:
|
|
|
+/// - Removes variants from low-complexity regions with entropy below the configured threshold
|
|
|
+/// (cf. [`Config::min_shannon_entropy`]).
|
|
|
+/// - Updates [`SomaticPipeStats::n_low_entropies`].
|
|
|
+///
|
|
|
+/// ### 13. Annotation: VEP (Variant Effect Predictor)
|
|
|
+/// Adds transcript-level annotations from Ensembl VEP, providing functional consequences,
|
|
|
+/// impacted genes, and regulatory features.
|
|
|
+///
|
|
|
+/// ### 14. Merging
|
|
|
+/// Merges variant collections into a unified [`Variants`] structure,
|
|
|
+/// preserving annotations and applying deduplication logic.
|
|
|
+///
|
|
|
+/// ### 15. Final Statistics and Saving
|
|
|
+/// - Saves final annotation stats, VEP summaries, and variant-level statistics.
|
|
|
+/// - Exports the final somatic variant set to both compressed JSON and `.bit` formats.
|
|
|
+///
|
|
|
+/// # Returns
|
|
|
+/// - `Ok(())` if all steps completed successfully.
|
|
|
+/// - `Err` if any tool fails, if file I/O fails, or if logical conditions are violated (e.g., pre-existing output).
|
|
|
+///
|
|
|
+/// # Errors
|
|
|
+/// - Returns early if the output file already exists and [`Config::somatic_pipe_force`] is `false`.
|
|
|
+/// - Wraps all component-specific errors using `anyhow::Result` with context.
|
|
|
+///
|
|
|
+/// # Side Effects
|
|
|
+/// - Runs external tools conditionally (e.g., [`ClairS`], [`DeepSomatic`]).
|
|
|
+/// - Creates intermediate directories and annotation JSONs for debugging/QC.
|
|
|
+/// - May consume significant compute time depending on the number of callers, annotations, and variants.
|
|
|
+///
|
|
|
+/// # TODOs
|
|
|
+/// - Support compressed intermediate files (`// TODO: GZ !!!`)
|
|
|
+/// - Improve filtering metrics reporting (currently not all filtered variants are tracked in final stats).
|
|
|
+///
|
|
|
+/// # Example Output Files
|
|
|
+/// - `tumoral_dir/sample123_somatic_variants.json.gz`
|
|
|
+/// - `tumoral_dir/sample123_somatic_variants.bit`
|
|
|
+/// - `stats_dir/sample123_annotations_*.json` (intermediate annotation snapshots)
|
|
|
+///
|
|
|
+/// # See Also
|
|
|
+/// - [`Annotations`] – core structure for managing per-variant metadata
|
|
|
+/// - [`Variants`] – the merged final variant structure
|
|
|
+/// - [`SomaticPipeStats`] – used for tracking variant counts throughout filtering
|
|
|
+///
|
|
|
+pub struct SomaticPipe {
|
|
|
+ /// Unique identifier for the sample.
|
|
|
pub id: String,
|
|
|
+ /// Configuration parameters for the pipeline.
|
|
|
pub config: Config,
|
|
|
- pub annotations: Annotations,
|
|
|
}
|
|
|
|
|
|
-impl Initialize for Somatic {
|
|
|
+impl Initialize for SomaticPipe {
|
|
|
+ /// Initializes a new `Somatic` instance with default annotations.
|
|
|
fn initialize(id: &str, config: crate::config::Config) -> anyhow::Result<Self> {
|
|
|
let id = id.to_string();
|
|
|
- Ok(Self {
|
|
|
- id,
|
|
|
- config,
|
|
|
- annotations: Annotations::default(),
|
|
|
- })
|
|
|
+ Ok(Self { id, config })
|
|
|
}
|
|
|
}
|
|
|
|
|
|
-#[derive(Debug, Default, Clone)]
|
|
|
-pub struct SomaticPipeStats {
|
|
|
- pub input: InputStats,
|
|
|
- pub n_constit_germline: usize,
|
|
|
- pub n_low_constit: usize,
|
|
|
- pub n_high_alt_constit: usize,
|
|
|
- pub n_high_alt_constit_gnomad: usize,
|
|
|
- pub n_low_entropies: usize,
|
|
|
-}
|
|
|
-
|
|
|
-#[derive(Debug, Default, Clone)]
|
|
|
-pub struct InputStats {
|
|
|
- pub solo_tumor: Vec<(Annotation, usize)>,
|
|
|
- pub solo_constit: Vec<(Annotation, usize)>,
|
|
|
- pub germline: Vec<(Annotation, usize)>,
|
|
|
- pub somatic: Vec<(Annotation, usize)>,
|
|
|
-}
|
|
|
-
|
|
|
-impl InputStats {
|
|
|
- pub fn from_collections(collections: &[VariantCollection]) -> Self {
|
|
|
- let mut stats = Self::default();
|
|
|
- for collection in collections.iter() {
|
|
|
- match collection.caller {
|
|
|
- Annotation::Callers(_, Sample::SoloTumor) => stats
|
|
|
- .solo_tumor
|
|
|
- .push((collection.caller.clone(), collection.variants.len())),
|
|
|
- Annotation::Callers(_, Sample::SoloConstit) => stats
|
|
|
- .solo_constit
|
|
|
- .push((collection.caller.clone(), collection.variants.len())),
|
|
|
- Annotation::Callers(_, Sample::Germline) => stats
|
|
|
- .germline
|
|
|
- .push((collection.caller.clone(), collection.variants.len())),
|
|
|
- Annotation::Callers(_, Sample::Somatic) => stats
|
|
|
- .somatic
|
|
|
- .push((collection.caller.clone(), collection.variants.len())),
|
|
|
- _ => (),
|
|
|
- };
|
|
|
- }
|
|
|
- stats
|
|
|
+impl ShouldRun for SomaticPipe {
|
|
|
+ fn should_run(&self) -> bool {
|
|
|
+ todo!()
|
|
|
}
|
|
|
}
|
|
|
|
|
|
-impl SomaticPipeStats {
|
|
|
- pub fn init(collections: &[VariantCollection]) -> Self {
|
|
|
- Self {
|
|
|
- input: InputStats::from_collections(collections),
|
|
|
- ..Default::default()
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- pub fn annot_init(&self, stats: &AnnotationsStats, json_path: &str) -> anyhow::Result<()> {
|
|
|
- let stats: Vec<(Vec<Annotation>, u64)> = stats
|
|
|
- .categorical
|
|
|
- .iter()
|
|
|
- .map(|e| {
|
|
|
- let anns = e
|
|
|
- .key()
|
|
|
- .split(" + ")
|
|
|
- .map(|k| k.parse())
|
|
|
- .collect::<anyhow::Result<Vec<Annotation>>>()
|
|
|
- .map_err(|err| {
|
|
|
- anyhow::anyhow!("Error while splitting key in AnnotationsStats.\n{err}")
|
|
|
- })?;
|
|
|
- Ok((anns, *e.value()))
|
|
|
- })
|
|
|
- .collect::<anyhow::Result<Vec<(Vec<Annotation>, u64)>>>()?;
|
|
|
-
|
|
|
- let callers_somatic_solo_tumor = [
|
|
|
- self.input
|
|
|
- .somatic
|
|
|
- .iter()
|
|
|
- .map(|(caller, _)| caller.clone())
|
|
|
- .collect::<Vec<Annotation>>(),
|
|
|
- self.input
|
|
|
- .solo_tumor
|
|
|
- .iter()
|
|
|
- .map(|(caller, _)| caller.clone())
|
|
|
- .collect(),
|
|
|
- ]
|
|
|
- .concat();
|
|
|
-
|
|
|
- let callers_germline_solo_constit = [
|
|
|
- self.input
|
|
|
- .germline
|
|
|
- .iter()
|
|
|
- .map(|(caller, _)| caller.clone())
|
|
|
- .collect::<Vec<Annotation>>(),
|
|
|
- self.input
|
|
|
- .solo_constit
|
|
|
- .iter()
|
|
|
- .map(|(caller, _)| caller.clone())
|
|
|
- .collect(),
|
|
|
- ]
|
|
|
- .concat();
|
|
|
-
|
|
|
- let mut with_germline: HashMap<String, HashMap<String, u64>> = HashMap::new();
|
|
|
- stats.iter().for_each(|(anns, v)| {
|
|
|
- if anns.iter().any(|a| {
|
|
|
- matches!(
|
|
|
- a,
|
|
|
- Annotation::Callers(_, Sample::SoloConstit)
|
|
|
- | Annotation::Callers(_, Sample::Germline)
|
|
|
- )
|
|
|
- }) {
|
|
|
- let n_by_tumor: Vec<(String, u64)> = callers_somatic_solo_tumor
|
|
|
- .iter()
|
|
|
- .flat_map(|tumor| {
|
|
|
- if anns.contains(tumor) {
|
|
|
- vec![(tumor.to_string(), *v)]
|
|
|
- } else {
|
|
|
- vec![]
|
|
|
- }
|
|
|
- })
|
|
|
- .collect();
|
|
|
-
|
|
|
- let mut germline_caller: Vec<String> = callers_germline_solo_constit
|
|
|
- .iter()
|
|
|
- .flat_map(|germ| {
|
|
|
- if anns.contains(germ) {
|
|
|
- vec![germ.to_string()]
|
|
|
- } else {
|
|
|
- vec![]
|
|
|
- }
|
|
|
- })
|
|
|
- .collect();
|
|
|
- germline_caller.sort();
|
|
|
- let germline_caller = germline_caller.join(" + ");
|
|
|
-
|
|
|
- n_by_tumor.iter().for_each(|(tumoral_caller, n)| {
|
|
|
- if let Some(row) = with_germline.get_mut(tumoral_caller) {
|
|
|
- if let Some(col) = row.get_mut(&germline_caller) {
|
|
|
- *col += *n;
|
|
|
- } else {
|
|
|
- row.insert(germline_caller.to_string(), *n);
|
|
|
- }
|
|
|
- } else {
|
|
|
- let mut row = HashMap::new();
|
|
|
- row.insert(germline_caller.to_string(), *n);
|
|
|
- with_germline.insert(tumoral_caller.to_string(), row);
|
|
|
- }
|
|
|
- });
|
|
|
- }
|
|
|
- });
|
|
|
-
|
|
|
- let mut germlines_callers: Vec<String> = with_germline
|
|
|
- .iter()
|
|
|
- .flat_map(|(_, r)| {
|
|
|
- r.iter()
|
|
|
- .map(|(k, _)| k.to_string())
|
|
|
- .collect::<Vec<String>>()
|
|
|
- })
|
|
|
- .collect();
|
|
|
- germlines_callers.sort();
|
|
|
- germlines_callers.dedup();
|
|
|
-
|
|
|
- let mut json = Vec::new();
|
|
|
- let mut lines: Vec<String> = with_germline
|
|
|
- .iter()
|
|
|
- .map(|(tumor, row)| {
|
|
|
- json.push(format!(
|
|
|
- "{{\"caller_name\": \"{tumor}\", \"germline\": [{}] }}",
|
|
|
- germlines_callers
|
|
|
- .iter()
|
|
|
- .map(|g| {
|
|
|
- let v = row.get(g).unwrap_or(&0);
|
|
|
- format!("{{\"{g}\": {v}}}")
|
|
|
- })
|
|
|
- .join(", ")
|
|
|
- ));
|
|
|
- format!(
|
|
|
- "{tumor}\t{}",
|
|
|
- germlines_callers
|
|
|
- .iter()
|
|
|
- .map(|g| {
|
|
|
- let v = row.get(g).unwrap_or(&0);
|
|
|
- format!("{g}: {v}")
|
|
|
- })
|
|
|
- .join("\t")
|
|
|
- )
|
|
|
- })
|
|
|
- .collect();
|
|
|
- lines.sort();
|
|
|
- println!("{}", lines.join("\n"));
|
|
|
-
|
|
|
- let json = format!("[{}]", json.join(", "));
|
|
|
- let mut file = File::create(json_path)?;
|
|
|
- file.write_all(json.as_bytes())?;
|
|
|
-
|
|
|
- Ok(())
|
|
|
- }
|
|
|
-}
|
|
|
-
|
|
|
-impl Run for Somatic {
|
|
|
+impl Run for SomaticPipe {
|
|
|
+ /// Executes the full somatic variant analysis pipeline.
|
|
|
fn run(&mut self) -> anyhow::Result<()> {
|
|
|
let config = self.config.clone();
|
|
|
let id = self.id.clone();
|
|
|
+ info!("Running somatic pipe for {id}.");
|
|
|
|
|
|
+ // Define output paths for the final somatic variants
|
|
|
let result_json = format!("{}/{id}_somatic_variants.json.gz", config.tumoral_dir(&id));
|
|
|
let result_bit = format!("{}/{id}_somatic_variants.bit", config.tumoral_dir(&id));
|
|
|
|
|
|
if Path::new(&result_json).exists() && !config.somatic_pipe_force {
|
|
|
- return Err(anyhow::anyhow!("already exists"));
|
|
|
+ return Err(anyhow::anyhow!(
|
|
|
+ "Somatic Pipe output already exists for {id}."
|
|
|
+ ));
|
|
|
}
|
|
|
|
|
|
- info!("Running somatic pipe for {id}.");
|
|
|
- let annotations = Arc::new(self.annotations.clone());
|
|
|
+ let mut annotations = Annotations::default();
|
|
|
|
|
|
- // Stats dir
|
|
|
+ // Create stats directory if it doesn't exist
|
|
|
let stats_dir = config.somatic_pipe_stats(&id);
|
|
|
if !Path::new(&stats_dir).exists() {
|
|
|
fs::create_dir(&stats_dir)?;
|
|
|
}
|
|
|
- // TODO: GZ !!!
|
|
|
- // LongphasePhase::initialize(&id, self.config.clone())?.run()?;
|
|
|
|
|
|
- // Initalize variants collections
|
|
|
- info!("Initialization prerequired pipe components...");
|
|
|
+ // Initialize and run any pre-required input if necessary
|
|
|
+ info!("Initialization prerequired pipe inputs...");
|
|
|
|
|
|
let mut to_run_if_req = create_should_run!(
|
|
|
&id,
|
|
|
@@ -276,13 +205,13 @@ impl Run for Somatic {
|
|
|
Savana,
|
|
|
DeepSomatic
|
|
|
);
|
|
|
- to_run_if_req.extend(create_should_run_normal_tumoral!(&id, &config, DeepVariant,));
|
|
|
+ to_run_if_req.extend(create_should_run_normal_tumoral!(&id, &config, DeepVariant));
|
|
|
|
|
|
info!("Running prerequired pipe components.");
|
|
|
-
|
|
|
run_if_required(&mut to_run_if_req)
|
|
|
.context("Failed to run a prerequired component of somatic pipe.")?;
|
|
|
|
|
|
+ // Initialize variant callers
|
|
|
let mut callers = init_somatic_callers!(
|
|
|
&id,
|
|
|
&config,
|
|
|
@@ -295,18 +224,21 @@ impl Run for Somatic {
|
|
|
|
|
|
callers.extend(init_solo_callers_normal_tumoral!(&id, &config, DeepVariant,));
|
|
|
|
|
|
+ // Load variants from each caller
|
|
|
info!("Loading variants.");
|
|
|
let mut variants_collections = load_variants(&mut callers, &annotations)
|
|
|
.map_err(|e| anyhow::anyhow!("Error while loading variants\n{e}"))?;
|
|
|
|
|
|
- info!("Loading Germline");
|
|
|
+ // Load germline variants using ClairS
|
|
|
+ info!("Loading germline variants.");
|
|
|
let clairs_germline =
|
|
|
ClairS::initialize(&id, self.config.clone())?.germline(&annotations)?;
|
|
|
variants_collections.push(clairs_germline);
|
|
|
|
|
|
+ // Initialize statistics
|
|
|
let mut somatic_stats = SomaticPipeStats::init(&variants_collections);
|
|
|
info!(
|
|
|
- "Variants collections from {} vcf ({} variants)",
|
|
|
+ "Variants collections loaded from {} vcf (total of {} variants loaded)",
|
|
|
variants_collections.len(),
|
|
|
variants_collections
|
|
|
.iter()
|
|
|
@@ -314,8 +246,7 @@ impl Run for Somatic {
|
|
|
.sum::<usize>()
|
|
|
);
|
|
|
|
|
|
- let mut annotations = Arc::try_unwrap(annotations)
|
|
|
- .map_err(|e| anyhow::anyhow!("Failed to unwrap Arc: {:?}", e))?;
|
|
|
+ // Initial annotation stats (caller annotations only)
|
|
|
let caller_cat_anns = |v: &Annotation| matches!(v, Annotation::Callers(_, _));
|
|
|
let annot_init = annotations.callers_stat(Some(Box::new(caller_cat_anns)));
|
|
|
somatic_stats.annot_init(
|
|
|
@@ -324,8 +255,11 @@ impl Run for Somatic {
|
|
|
)?;
|
|
|
annot_init.save_to_json(&format!("{stats_dir}/{id}_annotations_01.json"))?;
|
|
|
|
|
|
- // Filter: Variants neither Germline nor SoloConstit
|
|
|
- info!("Keeping somatic variants (variants neither in solo nor in germline).");
|
|
|
+ // Filter out germline and solo constitutional variants
|
|
|
+ info!(
|
|
|
+ "Keeping somatic variants (variants neither in solo {} nor in germline).",
|
|
|
+ config.normal_name
|
|
|
+ );
|
|
|
somatic_stats.n_constit_germline =
|
|
|
annotations.retain_variants(&mut variants_collections, |anns| {
|
|
|
!anns.iter().any(|ann| {
|
|
|
@@ -342,7 +276,7 @@ impl Run for Somatic {
|
|
|
"{stats_dir}/{id}_annotations_02_post_germline.json"
|
|
|
))?;
|
|
|
|
|
|
- // Annotation: BAM depth, n_alt
|
|
|
+ // Annotate with depth and number of alternate reads from constitutional BAM
|
|
|
info!("Reading Constit BAM file for depth and pileup annotation.");
|
|
|
variants_collections.iter().try_for_each(|c| {
|
|
|
c.annotate_with_constit_bam(&annotations, &self.config.normal_bam(&id), 150)
|
|
|
@@ -364,7 +298,7 @@ impl Run for Somatic {
|
|
|
})))
|
|
|
.save_to_json(&format!("{stats_dir}/{id}_annotations_03_bam.json"))?;
|
|
|
|
|
|
- // Filter: Remove LowConstitDepth from annotations and variants collections
|
|
|
+ // Filter based on low constitutional depth
|
|
|
info!(
|
|
|
"Removing variants when depth in constit bam < {}.",
|
|
|
self.config.somatic_min_constit_depth
|
|
|
@@ -403,17 +337,22 @@ impl Run for Somatic {
|
|
|
})))
|
|
|
.save_to_json(&format!("{stats_dir}/{id}_annotations_04_bam_filter.json"))?;
|
|
|
|
|
|
- // Annotation: Entropy
|
|
|
+ // Annotate variants with sequence entropy
|
|
|
info!(
|
|
|
"Entropy annotation from {} sequences.",
|
|
|
self.config.reference
|
|
|
);
|
|
|
variants_collections.iter().for_each(|c| {
|
|
|
- c.annotate_with_sequence_entropy(&annotations, &self.config.reference, 10, 150);
|
|
|
+ c.annotate_with_sequence_entropy(
|
|
|
+ &annotations,
|
|
|
+ &self.config.reference,
|
|
|
+ self.config.entropy_seq_len,
|
|
|
+ self.config.somatic_pipe_threads,
|
|
|
+ );
|
|
|
});
|
|
|
|
|
|
- // Annotation: Cosmic and GnomAD
|
|
|
- info!("Annotation with Cosmic and GnomAD.");
|
|
|
+ // Annotate with external databases like COSMIC and GnomAD
|
|
|
+ info!("Annotation with external databases like COSMIC and GnomAD.");
|
|
|
variants_collections
|
|
|
.iter()
|
|
|
.try_for_each(|c| -> anyhow::Result<()> {
|
|
|
@@ -430,31 +369,11 @@ impl Run for Somatic {
|
|
|
})))
|
|
|
.save_to_json(&format!("{stats_dir}/{id}_annotations_05_gnomad.json"))?;
|
|
|
|
|
|
- // Filter: Remove variants in Gnomad and in constit bam
|
|
|
+ // Filter: Remove variants present in GnomAD and have alt reads in constitutional sample
|
|
|
info!("Filtering out variants in GnomAD and in constit bam at low AF.");
|
|
|
- somatic_stats.n_high_alt_constit_gnomad =
|
|
|
- annotations.retain_variants(&mut variants_collections, |anns| {
|
|
|
- !anns
|
|
|
- .iter()
|
|
|
- .find_map(|a| {
|
|
|
- if let Annotation::GnomAD(gnomad) = a {
|
|
|
- Some(gnomad.gnomad_af > 0.0)
|
|
|
- } else {
|
|
|
- None
|
|
|
- }
|
|
|
- })
|
|
|
- .and_then(|gnomad_condition| {
|
|
|
- anns.iter()
|
|
|
- .find_map(|a| {
|
|
|
- if let Annotation::ConstitAlt(n_alt) = a {
|
|
|
- Some(*n_alt > 0)
|
|
|
- } else {
|
|
|
- None
|
|
|
- }
|
|
|
- })
|
|
|
- .map(|constit_alt_condition| gnomad_condition && constit_alt_condition)
|
|
|
- })
|
|
|
- .unwrap_or(false)
|
|
|
+ somatic_stats.n_high_alt_constit_gnomad = annotations
|
|
|
+ .retain_variants(&mut variants_collections, |anns| {
|
|
|
+ !is_gnomad_and_constit_alt(anns)
|
|
|
});
|
|
|
|
|
|
info!(
|
|
|
@@ -473,12 +392,13 @@ impl Run for Somatic {
|
|
|
"{stats_dir}/{id}_annotations_06_gnomad_filter.json"
|
|
|
))?;
|
|
|
|
|
|
- // Annotation low entropy
|
|
|
+ // Filter low entropy variants
|
|
|
annotations.low_shannon_entropy(self.config.min_shannon_entropy);
|
|
|
- // annotations.callers_stat();
|
|
|
|
|
|
- // Filtering low entropy for solo variants.
|
|
|
- info!("Filtering low entropies");
|
|
|
+ info!(
|
|
|
+ "Filtering out variants with low entropies ({})",
|
|
|
+ config.min_shannon_entropy
|
|
|
+ );
|
|
|
annotations
|
|
|
.callers_stat(Some(Box::new(|v| {
|
|
|
matches!(v, Annotation::Callers(_, _) | Annotation::LowEntropy)
|
|
|
@@ -489,14 +409,15 @@ impl Run for Somatic {
|
|
|
.retain_variants(&mut variants_collections, |anns| {
|
|
|
!anns.contains(&Annotation::LowEntropy)
|
|
|
});
|
|
|
+
|
|
|
annotations
|
|
|
.callers_stat(Some(Box::new(|v| matches!(v, Annotation::Callers(_, _)))))
|
|
|
.save_to_json(&format!(
|
|
|
"{stats_dir}/{id}_annotations_08_entropy_filter.json"
|
|
|
))?;
|
|
|
|
|
|
- // VEP
|
|
|
- info!("VEP annotation.");
|
|
|
+ // Final VEP annotation
|
|
|
+ info!("Annotation with VEP.");
|
|
|
variants_collections
|
|
|
.iter()
|
|
|
.try_for_each(|c| -> anyhow::Result<()> {
|
|
|
@@ -504,12 +425,14 @@ impl Run for Somatic {
|
|
|
ext_annot.annotate_vep(&c.variants, &annotations)?;
|
|
|
Ok(())
|
|
|
})?;
|
|
|
+
|
|
|
annotations
|
|
|
.callers_stat(Some(Box::new(caller_cat_anns)))
|
|
|
.save_to_json(&format!("{stats_dir}/{id}_annotations_09_vep.json"))?;
|
|
|
|
|
|
annotations.vep_stats()?;
|
|
|
|
|
|
+ // Merge all variants into a final collection
|
|
|
let variants = variants_collections.into_iter().fold(
|
|
|
Variants::default(),
|
|
|
|mut acc, variants_collection| {
|
|
|
@@ -550,3 +473,292 @@ pub fn const_stats(id: String, config: Config) -> anyhow::Result<()> {
|
|
|
|
|
|
Ok(())
|
|
|
}
|
|
|
+
|
|
|
+/// Holds statistical data for somatic variant pipeline processing,
|
|
|
+/// including summary counts and input categorization.
|
|
|
+#[derive(Debug, Default, Clone)]
|
|
|
+pub struct SomaticPipeStats {
|
|
|
+ /// Summary of input variant collections grouped by sample type.
|
|
|
+ pub input: InputStats,
|
|
|
+
|
|
|
+ /// Number of variants labeled as both constitutional and germline.
|
|
|
+ pub n_constit_germline: usize,
|
|
|
+
|
|
|
+ /// Number of variants in constitutional samples with low allele frequency.
|
|
|
+ pub n_low_constit: usize,
|
|
|
+
|
|
|
+ /// Number of variants in constitutional samples with high alternative allele count.
|
|
|
+ pub n_high_alt_constit: usize,
|
|
|
+
|
|
|
+ /// Number of high-alt constitutional variants that are also found in gnomAD.
|
|
|
+ pub n_high_alt_constit_gnomad: usize,
|
|
|
+
|
|
|
+ /// Number of variants filtered due to low entropy (indicative of low complexity regions).
|
|
|
+ pub n_low_entropies: usize,
|
|
|
+}
|
|
|
+
|
|
|
+impl SomaticPipeStats {
|
|
|
+ /// Initializes a `SomaticPipeStats` object with populated `InputStats` based on the
|
|
|
+ /// provided `VariantCollection`s.
|
|
|
+ pub fn init(collections: &[VariantCollection]) -> Self {
|
|
|
+ Self {
|
|
|
+ input: InputStats::from_collections(collections),
|
|
|
+ ..Default::default()
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ /// Generates a tumor-vs-germline annotation matrix and writes it to a JSON file.
|
|
|
+ ///
|
|
|
+ /// This method analyzes co-occurrence patterns between tumor and germline/constit
|
|
|
+ /// variant callers by iterating through annotation statistics from `AnnotationsStats`.
|
|
|
+ /// It builds a matrix where each row corresponds to a **tumor caller** and each column
|
|
|
+ /// corresponds to a **germline or constitutional caller**. Each cell contains the number
|
|
|
+ /// of variant calls that were annotated by both the tumor and the germline/constit caller.
|
|
|
+ ///
|
|
|
+ /// In addition to writing a structured JSON file to the specified path, the function also
|
|
|
+ /// prints a tab-separated human-readable matrix to stdout for convenience.
|
|
|
+ ///
|
|
|
+ /// # Parameters
|
|
|
+ ///
|
|
|
+ /// - `stats`: A reference to [`AnnotationsStats`] containing categorical annotations.
|
|
|
+ /// This object holds a frequency map where keys are combinations of `Annotation`
|
|
|
+ /// values (serialized as strings) and values are occurrence counts.
|
|
|
+ /// - `json_path`: The path where the resulting JSON output file should be written.
|
|
|
+ ///
|
|
|
+ /// # Output Formats
|
|
|
+ ///
|
|
|
+ /// ## JSON Output (`json_path`)
|
|
|
+ ///
|
|
|
+ /// The output JSON is an **array of tumor caller records**, each containing:
|
|
|
+ /// - `caller_name`: The name of the tumor caller (as a string).
|
|
|
+ /// - `germline`: An array of JSON objects, each with one key-value pair,
|
|
|
+ /// where the key is a germline/constit caller (or combination) and the value is the count.
|
|
|
+ ///
|
|
|
+ /// ### Example
|
|
|
+ /// ```json
|
|
|
+ /// [
|
|
|
+ /// {
|
|
|
+ /// "caller_name": "DeepVariant SoloTumor",
|
|
|
+ /// "germline": [
|
|
|
+ /// {"ClairS Germline": 99978},
|
|
|
+ /// {"ClairS Germline + DeepVariant SoloConstit": 4710570}
|
|
|
+ /// ]
|
|
|
+ /// },
|
|
|
+ /// {
|
|
|
+ /// "caller_name": "ClairS Somatic",
|
|
|
+ /// "germline": [
|
|
|
+ /// {"ClairS Germline": 944},
|
|
|
+ /// {"ClairS Germline + DeepVariant SoloConstit": 362}
|
|
|
+ /// ]
|
|
|
+ /// }
|
|
|
+ /// ]
|
|
|
+ /// ```
|
|
|
+ ///
|
|
|
+ /// ## Console Output (TSV)
|
|
|
+ ///
|
|
|
+ /// A tab-separated matrix is printed to stdout. Each row begins with a tumor caller,
|
|
|
+ /// followed by columns showing germline/constit caller combinations and their counts.
|
|
|
+ ///
|
|
|
+ /// # Notes
|
|
|
+ /// - Tumor callers are collected from `self.input.somatic` and `self.input.solo_tumor`.
|
|
|
+ /// - Germline/constit callers are collected from `self.input.germline` and `self.input.solo_constit`.
|
|
|
+ /// - Keys in the original `AnnotationsStats` map are split using `" + "` and parsed into `Annotation` enums.
|
|
|
+ /// - Germline keys are sorted and joined to form canonical column labels.
|
|
|
+ /// - Tumor and germline annotations are matched by exact `Annotation` equality.
|
|
|
+ ///
|
|
|
+ /// # Errors
|
|
|
+ /// Returns an error if:
|
|
|
+ /// - Any annotation key in `AnnotationsStats` fails to parse into a valid `Annotation`.
|
|
|
+ /// - The JSON output file cannot be created or written to.
|
|
|
+ ///
|
|
|
+ /// # Dependencies
|
|
|
+ /// Requires that `Annotation` implements `FromStr`, `ToString`, `PartialEq`, and `Clone`.
|
|
|
+ ///
|
|
|
+ /// # Example Usage
|
|
|
+ /// ```rust
|
|
|
+ /// let mut stats = SomaticPipeStats::init(&collections);
|
|
|
+ /// stats.annot_init(&annotation_stats, "output/matrix.json")?;
|
|
|
+ /// ```
|
|
|
+ pub fn annot_init(&self, stats: &AnnotationsStats, json_path: &str) -> anyhow::Result<()> {
|
|
|
+ // Parse annotations from stats
|
|
|
+ let stats: Vec<(Vec<Annotation>, u64)> = stats
|
|
|
+ .categorical
|
|
|
+ .iter()
|
|
|
+ .map(|e| {
|
|
|
+ let anns = e
|
|
|
+ .key()
|
|
|
+ .split(" + ")
|
|
|
+ .map(|k| k.parse())
|
|
|
+ .collect::<anyhow::Result<Vec<Annotation>>>()
|
|
|
+ .map_err(|err| {
|
|
|
+ anyhow::anyhow!("Error while splitting key in AnnotationsStats.\n{err}")
|
|
|
+ })?;
|
|
|
+ Ok((anns, *e.value()))
|
|
|
+ })
|
|
|
+ .collect::<anyhow::Result<Vec<(Vec<Annotation>, u64)>>>()?;
|
|
|
+
|
|
|
+ // Collect tumor and germline callers from input stats
|
|
|
+ let callers_somatic_solo_tumor = [
|
|
|
+ self.input
|
|
|
+ .somatic
|
|
|
+ .iter()
|
|
|
+ .map(|(caller, _)| caller.clone())
|
|
|
+ .collect::<Vec<Annotation>>(),
|
|
|
+ self.input
|
|
|
+ .solo_tumor
|
|
|
+ .iter()
|
|
|
+ .map(|(caller, _)| caller.clone())
|
|
|
+ .collect(),
|
|
|
+ ]
|
|
|
+ .concat();
|
|
|
+
|
|
|
+ let callers_germline_solo_constit = [
|
|
|
+ self.input
|
|
|
+ .germline
|
|
|
+ .iter()
|
|
|
+ .map(|(caller, _)| caller.clone())
|
|
|
+ .collect::<Vec<Annotation>>(),
|
|
|
+ self.input
|
|
|
+ .solo_constit
|
|
|
+ .iter()
|
|
|
+ .map(|(caller, _)| caller.clone())
|
|
|
+ .collect(),
|
|
|
+ ]
|
|
|
+ .concat();
|
|
|
+
|
|
|
+ // Build a matrix of tumor vs germline hits
|
|
|
+ let mut with_germline: HashMap<String, HashMap<String, u64>> = HashMap::new();
|
|
|
+ stats.iter().for_each(|(anns, v)| {
|
|
|
+ // Only proceed if this annotation includes a germline/constit sample
|
|
|
+ if anns.iter().any(|a| {
|
|
|
+ matches!(
|
|
|
+ a,
|
|
|
+ Annotation::Callers(_, Sample::SoloConstit)
|
|
|
+ | Annotation::Callers(_, Sample::Germline)
|
|
|
+ )
|
|
|
+ }) {
|
|
|
+ // Find all tumor callers present in this annotation set
|
|
|
+ let n_by_tumor: Vec<(String, u64)> = callers_somatic_solo_tumor
|
|
|
+ .iter()
|
|
|
+ .flat_map(|tumor| {
|
|
|
+ if anns.contains(tumor) {
|
|
|
+ vec![(tumor.to_string(), *v)]
|
|
|
+ } else {
|
|
|
+ vec![]
|
|
|
+ }
|
|
|
+ })
|
|
|
+ .collect();
|
|
|
+
|
|
|
+ // Build a normalized germline key
|
|
|
+ let mut germline_caller: Vec<String> = callers_germline_solo_constit
|
|
|
+ .iter()
|
|
|
+ .flat_map(|germ| {
|
|
|
+ if anns.contains(germ) {
|
|
|
+ vec![germ.to_string()]
|
|
|
+ } else {
|
|
|
+ vec![]
|
|
|
+ }
|
|
|
+ })
|
|
|
+ .collect();
|
|
|
+ germline_caller.sort();
|
|
|
+ let germline_caller = germline_caller.join(" + ");
|
|
|
+
|
|
|
+ // Update matrix: tumor -> germline -> count
|
|
|
+ n_by_tumor.iter().for_each(|(tumoral_caller, n)| {
|
|
|
+ if let Some(row) = with_germline.get_mut(tumoral_caller) {
|
|
|
+ if let Some(col) = row.get_mut(&germline_caller) {
|
|
|
+ *col += *n;
|
|
|
+ } else {
|
|
|
+ row.insert(germline_caller.to_string(), *n);
|
|
|
+ }
|
|
|
+ } else {
|
|
|
+ let mut row = HashMap::new();
|
|
|
+ row.insert(germline_caller.to_string(), *n);
|
|
|
+ with_germline.insert(tumoral_caller.to_string(), row);
|
|
|
+ }
|
|
|
+ });
|
|
|
+ }
|
|
|
+ });
|
|
|
+
|
|
|
+ // Extract all unique germline caller labels
|
|
|
+ let mut germlines_callers: Vec<String> = with_germline
|
|
|
+ .iter()
|
|
|
+ .flat_map(|(_, r)| {
|
|
|
+ r.iter()
|
|
|
+ .map(|(k, _)| k.to_string())
|
|
|
+ .collect::<Vec<String>>()
|
|
|
+ })
|
|
|
+ .collect();
|
|
|
+ germlines_callers.sort();
|
|
|
+ germlines_callers.dedup();
|
|
|
+
|
|
|
+ // Print a readable tab-separated matrix
|
|
|
+ let mut json = Vec::new();
|
|
|
+ let mut lines: Vec<String> = with_germline
|
|
|
+ .iter()
|
|
|
+ .map(|(tumor, row)| {
|
|
|
+ json.push(format!(
|
|
|
+ "{{\"caller_name\": \"{tumor}\", \"germline\": [{}] }}",
|
|
|
+ germlines_callers
|
|
|
+ .iter()
|
|
|
+ .map(|g| {
|
|
|
+ let v = row.get(g).unwrap_or(&0);
|
|
|
+ format!("{{\"{g}\": {v}}}")
|
|
|
+ })
|
|
|
+ .join(", ")
|
|
|
+ ));
|
|
|
+ format!(
|
|
|
+ "{tumor}\t{}",
|
|
|
+ germlines_callers
|
|
|
+ .iter()
|
|
|
+ .map(|g| {
|
|
|
+ let v = row.get(g).unwrap_or(&0);
|
|
|
+ format!("{g}: {v}")
|
|
|
+ })
|
|
|
+ .join("\t")
|
|
|
+ )
|
|
|
+ })
|
|
|
+ .collect();
|
|
|
+ lines.sort();
|
|
|
+ println!("{}", lines.join("\n"));
|
|
|
+
|
|
|
+ // Write JSON to file
|
|
|
+ let json = format!("[{}]", json.join(", "));
|
|
|
+ let mut file = File::create(json_path)?;
|
|
|
+ file.write_all(json.as_bytes())?;
|
|
|
+
|
|
|
+ Ok(())
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+#[derive(Debug, Default, Clone)]
|
|
|
+pub struct InputStats {
|
|
|
+ pub solo_tumor: Vec<(Annotation, usize)>,
|
|
|
+ pub solo_constit: Vec<(Annotation, usize)>,
|
|
|
+ pub germline: Vec<(Annotation, usize)>,
|
|
|
+ pub somatic: Vec<(Annotation, usize)>,
|
|
|
+}
|
|
|
+
|
|
|
+impl InputStats {
|
|
|
+ pub fn from_collections(collections: &[VariantCollection]) -> Self {
|
|
|
+ let mut stats = Self::default();
|
|
|
+ for collection in collections.iter() {
|
|
|
+ match collection.caller {
|
|
|
+ Annotation::Callers(_, Sample::SoloTumor) => stats
|
|
|
+ .solo_tumor
|
|
|
+ .push((collection.caller.clone(), collection.variants.len())),
|
|
|
+ Annotation::Callers(_, Sample::SoloConstit) => stats
|
|
|
+ .solo_constit
|
|
|
+ .push((collection.caller.clone(), collection.variants.len())),
|
|
|
+ Annotation::Callers(_, Sample::Germline) => stats
|
|
|
+ .germline
|
|
|
+ .push((collection.caller.clone(), collection.variants.len())),
|
|
|
+ Annotation::Callers(_, Sample::Somatic) => stats
|
|
|
+ .somatic
|
|
|
+ .push((collection.caller.clone(), collection.variants.len())),
|
|
|
+ _ => (),
|
|
|
+ };
|
|
|
+ }
|
|
|
+ stats
|
|
|
+ }
|
|
|
+}
|
