From 7b4d5ebbf56c55ea55d35290270e388814212dc6 Mon Sep 17 00:00:00 2001 From: Francesco Terenzi Date: Wed, 28 Feb 2024 14:47:44 +0000 Subject: [PATCH] v3.0.5 --- CHANGELOG.md | 3 + Cargo.lock | 2 +- Cargo.toml | 2 +- src/clap_app.rs | 11 ++- src/genotype.rs | 34 ++++++++- src/main.rs | 15 +++- src/process.rs | 198 ++++++++++++++++++++++++++++++++++++++++++++++++ src/subclone.rs | 4 + 8 files changed, 258 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0abb7bbe..f64d430e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,9 @@ # Changelog The semantic versioning is kind of random. +## 3.0.5 +- specfify the neutral background mutation rate for the exponential phase + ## 3.0.4 ### BugFix - When cloning a cell upon division, **do not** set `last_division_t` of the daughter cell to 0 diff --git a/Cargo.lock b/Cargo.lock index f323390e..a1a8232b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -323,7 +323,7 @@ checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" [[package]] name = "hsc" -version = "3.0.4" +version = "3.0.5" dependencies = [ "anyhow", "chrono", diff --git a/Cargo.toml b/Cargo.toml index 88e52e80..1624f7f7 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "hsc" -version = "3.0.4" +version = "3.0.5" edition = "2021" [dependencies] diff --git a/src/clap_app.rs b/src/clap_app.rs index 1e0067a5..a362563a 100644 --- a/src/clap_app.rs +++ b/src/clap_app.rs @@ -60,10 +60,15 @@ fn fitness_in_range(s: &str) -> Result { #[derive(Args, Debug, Clone)] #[group(required = false, multiple = true)] pub struct NeutralMutationRate { + /// Poisson rate for the exponential growing phase for the background + /// mutations, leave empty to simulate just a Moran process with fixed + /// population size + #[arg(long, requires = "mu_division_exp")] + pub mu_background_exp: Option, /// Poisson rate for the exponential growing phase, leave empty to simulate /// just a Moran process with fixed population size #[arg(long)] - pub mu_exp: Option, + pub mu_division_exp: Option, /// Background mutation rate (all neutral mutations **not** occuring in the /// mitotic phase) for for the constant population phase /// units: [mut/year] @@ -271,8 +276,8 @@ impl Cli { }; // Exp - let options_exponential = cli.neutral_rate.mu_exp.map(|_| { - assert!(cli.neutral_rate.mu_exp.is_some()); + let options_exponential = cli.neutral_rate.mu_division_exp.map(|_| { + assert!(cli.neutral_rate.mu_division_exp.is_some()); let process_options = ProcessOptions { path: cli.path, snapshots: snapshots.clone(), diff --git a/src/genotype.rs b/src/genotype.rs index 4bf05bf2..6a6e17cb 100644 --- a/src/genotype.rs +++ b/src/genotype.rs @@ -5,7 +5,10 @@ use rustc_hash::FxHashMap; use std::{fs, path::Path}; use uuid::Uuid; -use crate::stemcell::StemCell; +use crate::{ + process::{Exponential, Moran}, + stemcell::StemCell, +}; pub type Variant = Uuid; @@ -39,7 +42,8 @@ impl NeutralMutationPoisson { //! mutations acquired upon cell-division and the other modelling the //! acquisition of neutral background mutations, i.e. all mutations //! occuring not during cell-division. - ensure!(lambda_division > 0., "invalid value of lambda_background"); + ensure!(lambda_division > 0., "invalid value of lambda_division"); + ensure!(lambda_background > 0., "invalid value of lambda_background"); Ok(Self { background: lambda_background, division: Poisson::new(lambda_division) @@ -154,7 +158,7 @@ impl Sfs { /// mutations. /// /// To plot it, plot on the x-axis the keys (number of mutations) and on the -/// y-axis the keys (the number of cells with those mutations). +/// y-axis the values (the number of cells with those mutations). pub struct MutationalBurden(pub FxHashMap); impl MutationalBurden { @@ -175,6 +179,30 @@ impl MutationalBurden { Ok(MutationalBurden(burden)) } + pub fn from_moran(moran: &Moran, verbosity: u8) -> anyhow::Result { + //! Create the single-cell mutational burden from all cells in the Moran + //! process + let cells: Vec<&StemCell> = moran + .subclones + .get_cells_with_clones_idx() + .iter() + .map(|ele| ele.0) + .collect(); + MutationalBurden::from_cells(&cells, verbosity) + } + + pub fn from_exp(exp: &Exponential, verbosity: u8) -> anyhow::Result { + //! Create the single-cell mutational burden from all cells in the Moran + //! process + let cells: Vec<&StemCell> = exp + .subclones + .get_cells_with_clones_idx() + .iter() + .map(|ele| ele.0) + .collect(); + MutationalBurden::from_cells(&cells, verbosity) + } + pub fn save(&self, path2file: &Path, verbosity: u8) -> anyhow::Result<()> { let path2file = path2file.with_extension("json"); let burden = diff --git a/src/main.rs b/src/main.rs index 403e0412..0403ff69 100644 --- a/src/main.rs +++ b/src/main.rs @@ -124,16 +124,25 @@ fn main() { }; let mut moran = if let Some(options) = app.options_exponential.as_ref() { - let rate = app + let mu_div_exp = app .neutral_rate - .mu_exp + .mu_division_exp + .expect("find no exp neutral rate but options_exp"); + let mu_back_exp = app + .neutral_rate + .mu_division_exp .expect("find no exp neutral rate but options_exp"); // use 1 as we dont care about time during the exp growth phase let rates = subclones.gillespie_rates(&app.fitness, 1.0, rng); // we assume no background mutation for the exponential growing phase let mut exp = Exponential::new( subclones, - Distributions::new(u, rate, rate, app.options_moran.gillespie_options.verbosity), + Distributions::new( + u, + mu_back_exp, + mu_div_exp, + app.options_moran.gillespie_options.verbosity, + ), prolfieration, options.gillespie_options.verbosity, ); diff --git a/src/process.rs b/src/process.rs index 90c9fa15..d0685efd 100644 --- a/src/process.rs +++ b/src/process.rs @@ -495,3 +495,201 @@ impl AdvanceStep for Moran { state.population = Variants::variant_counts(&self.subclones); } } + +#[cfg(test)] +mod tests { + use std::num::{NonZeroU64, NonZeroU8}; + + use quickcheck_macros::quickcheck; + use rand::SeedableRng; + use rand_chacha::ChaCha8Rng; + + use crate::{genotype, proliferation::DivisionAndBackgroundMutationsProliferation}; + + use super::*; + + struct MoranFitVariant { + tot_cells: u64, + process: Moran, + } + + struct MoranNoFitVariant { + tot_cells: u64, + process: Moran, + } + + fn create_moran_fit_variants(cells: NonZeroU64) -> MoranFitVariant { + MoranFitVariant { + tot_cells: cells.get(), + process: create_moran(true, cells), + } + } + + fn create_moran_no_fit_variants(cells: NonZeroU64) -> MoranNoFitVariant { + MoranNoFitVariant { + tot_cells: cells.get(), + process: create_moran(false, cells), + } + } + + fn create_moran(fit_variants: bool, cells: NonZeroU64) -> Moran { + let u = if fit_variants { 1. } else { 0. }; + Moran::new( + ProcessOptions { + path: PathBuf::default(), + snapshots: VecDeque::from([Snapshot { + cells2sample: 2, + time: 10., + }]), + }, + SubClones::new(vec![StemCell::new(); cells.get() as usize], 3), + 0f32, + SavingOptions { + filename: PathBuf::default(), + save_sfs_only: true, + save_population: true, + }, + Distributions::new(u, 10f32, 1f32, 0), + MutateUponDivision::DivisionAndBackgroundMutations( + DivisionAndBackgroundMutationsProliferation, + ), + 0, + ) + } + + fn assert_pop_const(moran: &Moran, cells: u64) { + assert_eq!(moran.subclones.compute_tot_cells(), cells); + } + + fn assert_not_all_cells_in_wild_type(moran: &Moran, cells: u64) { + assert_ne!(moran.subclones.get_clone_unchecked(0).cell_count(), cells); + } + + fn assert_all_cells_in_wild_type(moran: &Moran, cells: u64) { + assert_eq!(moran.subclones.get_clone_unchecked(0).cell_count(), cells); + assert_eq!(moran.subclones.the_only_one_subclone_present().unwrap(), 0); + } + + #[quickcheck] + fn advance_moran_no_fit_variant_test(cells: NonZeroU8, seed: u64) { + let mut moran = create_moran_no_fit_variants(NonZeroU64::new(cells.get() as u64).unwrap()); + let reaction = NextReaction { + time: 12.2, + event: 0, + }; + let burden_before = genotype::MutationalBurden::from_moran(&moran.process, 0).unwrap(); + assert_pop_const(&moran.process, moran.tot_cells); + assert_all_cells_in_wild_type(&moran.process, moran.tot_cells); + let rng = &mut ChaCha8Rng::seed_from_u64(seed); + + moran.process.advance_step(reaction, rng); + + assert_pop_const(&moran.process, moran.tot_cells); + assert_all_cells_in_wild_type(&moran.process, moran.tot_cells); + let burden_after = genotype::MutationalBurden::from_moran(&moran.process, 0).unwrap(); + assert!(burden_before.0.keys().sum::() < burden_after.0.keys().sum::()); + } + + #[quickcheck] + fn advance_moran_test(cells: NonZeroU8, seed: u64) { + let mut moran = create_moran_fit_variants(NonZeroU64::new(cells.get() as u64).unwrap()); + let reaction = NextReaction { + time: 12.2, + event: 0, + }; + let burden_before = genotype::MutationalBurden::from_moran(&moran.process, 0).unwrap(); + assert_pop_const(&moran.process, moran.tot_cells); + assert_all_cells_in_wild_type(&moran.process, moran.tot_cells); + let rng = &mut ChaCha8Rng::seed_from_u64(seed); + + moran.process.advance_step(reaction, rng); + + assert_pop_const(&moran.process, moran.tot_cells); + assert_not_all_cells_in_wild_type(&moran.process, moran.tot_cells); + let burden_after = genotype::MutationalBurden::from_moran(&moran.process, 0).unwrap(); + assert!(burden_before.0.keys().sum::() < burden_after.0.keys().sum::()); + } + struct ExpFitVariant { + tot_cells: u64, + process: Exponential, + } + + struct ExpNoFitVariant { + tot_cells: u64, + process: Exponential, + } + + fn create_exp_fit_variants(cells: NonZeroU64) -> ExpFitVariant { + ExpFitVariant { + tot_cells: cells.get(), + process: create_exp(true, cells), + } + } + + fn create_exp_no_fit_variants(cells: NonZeroU64) -> ExpNoFitVariant { + ExpNoFitVariant { + tot_cells: cells.get(), + process: create_exp(false, cells), + } + } + + fn create_exp(fit_variants: bool, cells: NonZeroU64) -> Exponential { + let u = if fit_variants { 1. } else { 0. }; + Exponential::new( + SubClones::new(vec![StemCell::new(); cells.get() as usize], 3), + Distributions::new(u, 10f32, 1f32, 0), + MutateUponDivision::DivisionAndBackgroundMutations( + DivisionAndBackgroundMutationsProliferation, + ), + 0, + ) + } + + fn assert_not_all_cells_in_wild_type_exp(exp: &Exponential, cells: u64) { + assert_ne!(exp.subclones.get_clone_unchecked(0).cell_count(), cells); + } + + fn assert_all_cells_in_wild_type_exp(exp: &Exponential, cells: u64) { + assert_eq!(exp.subclones.get_clone_unchecked(0).cell_count(), cells); + assert_eq!(exp.subclones.the_only_one_subclone_present().unwrap(), 0); + } + + #[quickcheck] + fn advance_exp_no_fit_variant_test(cells: NonZeroU8, seed: u64) { + let mut exp = create_exp_no_fit_variants(NonZeroU64::new(cells.get() as u64).unwrap()); + let reaction = NextReaction { + time: 12.2, + event: 0, + }; + let burden_before = genotype::MutationalBurden::from_exp(&exp.process, 0).unwrap(); + let old_cells = exp.process.subclones.compute_tot_cells(); + assert_all_cells_in_wild_type_exp(&exp.process, exp.tot_cells); + let rng = &mut ChaCha8Rng::seed_from_u64(seed); + + exp.process.advance_step(reaction, rng); + + assert_eq!(exp.process.subclones.compute_tot_cells(), old_cells + 1); + let burden_after = genotype::MutationalBurden::from_exp(&exp.process, 0).unwrap(); + assert!(burden_before.0.keys().sum::() < burden_after.0.keys().sum::()); + } + + #[quickcheck] + fn advance_exp_test(cells: NonZeroU8, seed: u64) { + let mut exp = create_exp_fit_variants(NonZeroU64::new(cells.get() as u64).unwrap()); + let reaction = NextReaction { + time: 12.2, + event: 0, + }; + let burden_before = genotype::MutationalBurden::from_exp(&exp.process, 0).unwrap(); + let old_cells = exp.process.subclones.compute_tot_cells(); + assert_all_cells_in_wild_type_exp(&exp.process, exp.tot_cells); + let rng = &mut ChaCha8Rng::seed_from_u64(seed); + + exp.process.advance_step(reaction, rng); + + assert_eq!(exp.process.subclones.compute_tot_cells(), old_cells + 1); + assert_not_all_cells_in_wild_type_exp(&exp.process, exp.tot_cells); + let burden_after = genotype::MutationalBurden::from_exp(&exp.process, 0).unwrap(); + assert!(burden_before.0.keys().sum::() < burden_after.0.keys().sum::()); + } +} diff --git a/src/subclone.rs b/src/subclone.rs index 57c4969c..2af6944f 100644 --- a/src/subclone.rs +++ b/src/subclone.rs @@ -250,6 +250,10 @@ impl SubClones { self.0.get(id) } + pub fn get_clone_unchecked(&self, id: usize) -> &SubClone { + &self.0[id] + } + pub fn get_mut_clone_unchecked(&mut self, id: usize) -> &mut SubClone { &mut self.0[id] }