From b9db19f30480bb66c6c6e80cfbf11c3d6761ef66 Mon Sep 17 00:00:00 2001 From: Francesco Terenzi Date: Tue, 5 Mar 2024 10:51:39 +0000 Subject: [PATCH] v4.1.0 --- CHANGELOG.md | 3 + Cargo.lock | 2 +- Cargo.toml | 2 +- src/clap_app.rs | 15 +++- src/genotype.rs | 2 +- src/main.rs | 35 +++++--- src/process.rs | 108 ++++++++++------------ src/proliferation.rs | 209 +++++++++++++++++-------------------------- src/stemcell.rs | 22 ++++- src/subclone.rs | 8 +- 10 files changed, 198 insertions(+), 208 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 10a41031..856d6a8f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,9 @@ # Changelog The semantic versioning is kind of random. +## 4.1.0 +- introduce asymmetric divisions and refactor the neutral mutations + ## 4.0.0 - Add new clap app cmds and organisation diff --git a/Cargo.lock b/Cargo.lock index d04f0a1e..1790369b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -323,7 +323,7 @@ checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" [[package]] name = "hsc" -version = "4.0.0" +version = "4.1.0" dependencies = [ "anyhow", "chrono", diff --git a/Cargo.toml b/Cargo.toml index 124d1e92..8755e364 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "hsc" -version = "4.0.0" +version = "4.1.0" edition = "2021" [dependencies] diff --git a/src/clap_app.rs b/src/clap_app.rs index e1478e7a..89253cf6 100644 --- a/src/clap_app.rs +++ b/src/clap_app.rs @@ -75,15 +75,18 @@ pub struct ExponentialThenMoran { #[derive(Args, Debug, Clone)] pub struct MoranPhase { - /// Rate of accumulation of neutral background mutations per year + /// Rate of accumulation of neutral background mutations per year in the fixed-size population phase #[arg(long, default_value_t = 14.)] pub mu_background: f32, - /// Rate of accumulation of neutral mutations upon division per year + /// Rate of accumulation of neutral mutations upon division per year in the fixed-size population phase #[arg(long, default_value_t = 1.14)] pub mu_division: f32, - /// Inter-division time for the wild-type cells + /// Inter-division time for the wild-type cells in the fixed-size population phase #[arg(long, default_value_t = 1.)] pub tau: f32, + /// Probability of asymmetric division in the fixed-size population phase + #[arg(long, default_value_t = 0.)] + pub asymmetric: f32, } #[derive(Args, Debug, Clone)] @@ -98,6 +101,9 @@ pub struct ExponentialPhase { /// Inter-division time for the wild-type cells #[arg(long, default_value_t = 0.065)] pub tau_exp: f32, + /// Probability of asymmetric division + #[arg(long, default_value_t = 0.)] + pub asymmetric_exp: f32, } #[derive(Debug, Parser)] // requires `derive` feature @@ -295,6 +301,7 @@ impl Cli { tau: moran.tau, mu_background: moran.mu_background, mu_division: moran.mu_division, + asymmetric_prob: moran.asymmetric, }; (options_moran, None) } @@ -315,6 +322,7 @@ impl Cli { tau: exp_moran.moran.tau, mu_background: exp_moran.moran.mu_background, mu_division: exp_moran.moran.mu_division, + asymmetric_prob: exp_moran.moran.asymmetric, }; let options_exponential = SimulationOptionsExp { gillespie_options: Options { @@ -329,6 +337,7 @@ impl Cli { tau: exp_moran.exponential.tau_exp, mu_background: exp_moran.exponential.mu_background_exp, mu_division: exp_moran.exponential.mu_division_exp, + asymmetric_prob: exp_moran.exponential.asymmetric_exp, }; (options_moran, Some(options_exponential)) } diff --git a/src/genotype.rs b/src/genotype.rs index 30aa26fa..dd6d2b87 100644 --- a/src/genotype.rs +++ b/src/genotype.rs @@ -125,7 +125,7 @@ impl Sfs { //! Compute the SFS from the stem cell population. if verbosity > 0 { println!("computing the SFS from {} cells", cells.len()); - if verbosity > 1 { + if verbosity > 2 { println!("computing the SFS from {:#?}", &cells); } } diff --git a/src/main.rs b/src/main.rs index 06c64eca..3acd0a87 100644 --- a/src/main.rs +++ b/src/main.rs @@ -4,10 +4,7 @@ use chrono::Utc; use clap_app::Parallel; use hsc::{ process::{Exponential, Moran, ProcessOptions, SavingCells, SavingOptions, Snapshot}, - proliferation::{ - DivisionAndBackgroundMutationsProliferation, DivisionMutationsProliferation, - MutateUponDivision, - }, + proliferation::{Division, NeutralMutations, Proliferation}, stemcell::StemCell, subclone::{Distributions, Fitness, SubClones, Variants}, write2file, @@ -15,6 +12,7 @@ use hsc::{ use indicatif::ParallelProgressIterator; use rand::SeedableRng; use rand_chacha::ChaCha8Rng; +use rand_distr::Bernoulli; use rayon::prelude::{IntoParallelIterator, ParallelIterator}; use sosa::{simulate, CurrentState, Options, StopReason}; use std::collections::VecDeque; @@ -27,6 +25,7 @@ pub struct SimulationOptionsMoran { tau: f32, mu_background: f32, mu_division: f32, + asymmetric_prob: f32, process_options: ProcessOptions, gillespie_options: Options, save_sfs_only: bool, @@ -39,6 +38,7 @@ pub struct SimulationOptionsExp { gillespie_options: Options, mu_background: f32, mu_division: f32, + asymmetric_prob: f32, } #[derive(Clone, Debug)] @@ -118,18 +118,21 @@ fn main() { ) .replace('.', "dot") .into(); - - let prolfieration = if app.background { - MutateUponDivision::DivisionAndBackgroundMutations( - DivisionAndBackgroundMutationsProliferation, - ) + let neutral_mutation = if app.background { + NeutralMutations::UponDivisionAndBackground } else { - MutateUponDivision::DivisionMutations(DivisionMutationsProliferation) + NeutralMutations::UponDivision }; let mut moran = if let Some(options) = app.options_exponential.as_ref() { // use 1 as we dont care about time during the exp growth phase let rates = subclones.gillespie_rates(&app.fitness, 1.0 / options.tau, rng); + let division = if (options.asymmetric_prob - 0f32).abs() > f32::EPSILON { + Division::Asymmetric(Bernoulli::new(options.asymmetric_prob as f64).unwrap()) + } else { + Division::Symmetric + }; + let proliferation = Proliferation::new(neutral_mutation, division); let mut exp = Exponential::new( subclones, Distributions::new( @@ -140,7 +143,7 @@ fn main() { options.mu_division, app.options_moran.gillespie_options.verbosity, ), - prolfieration, + proliferation, options.gillespie_options.verbosity, ); @@ -174,6 +177,14 @@ fn main() { rng, ) } else { + let division = if (app.options_moran.asymmetric_prob - 0f32).abs() > f32::EPSILON { + Division::Asymmetric( + Bernoulli::new(app.options_moran.asymmetric_prob as f64).unwrap(), + ) + } else { + Division::Symmetric + }; + let proliferation = Proliferation::new(neutral_mutation, division); Moran::new( app.options_moran.process_options.clone(), subclones, @@ -184,7 +195,7 @@ fn main() { save_population: app.options_moran.save_population, }, distributions, - prolfieration, + proliferation, app.options_moran.gillespie_options.verbosity, ) }; diff --git a/src/process.rs b/src/process.rs index 714fc723..8f778168 100644 --- a/src/process.rs +++ b/src/process.rs @@ -1,5 +1,5 @@ use crate::genotype::{MutationalBurden, Sfs}; -use crate::proliferation::{MutateUponDivision, Proliferation}; +use crate::proliferation::{NeutralMutations, Proliferation}; use crate::stemcell::{assign_background_mutations, StemCell}; use crate::subclone::{save_variant_fraction, CloneId, Distributions, SubClones, Variants}; use crate::MAX_SUBCLONES; @@ -65,7 +65,7 @@ pub struct Exponential { pub counter_divisions: usize, pub verbosity: u8, pub distributions: Distributions, - pub proliferation: MutateUponDivision, + pub proliferation: Proliferation, pub time: f32, } @@ -73,7 +73,7 @@ impl Exponential { pub fn new( initial_subclones: SubClones, distributions: Distributions, - proliferation: MutateUponDivision, + proliferation: Proliferation, verbosity: u8, ) -> Exponential { let hsc = Exponential { @@ -115,30 +115,27 @@ impl Exponential { save_population, proliferation: self.proliferation, }; - match moran.proliferation { - MutateUponDivision::DivisionAndBackgroundMutations(_) => { - // this is important: we update all background mutations at this - // time such that all cells are all on the same page. - // Since background mutations are implemented at each division, and - // cells do not proliferate at the same rate, we need to correct and - // update the background mutations at the timepoint corresponding - // to sampling step, i.e. before saving. - if self.verbosity > 0 { - println!("updating the neutral background mutations for all cells"); - } - for stem_cell in moran.subclones.get_mut_cells() { - assign_background_mutations( - stem_cell, - self.time, - &self.distributions.neutral_poisson, - rng, - self.verbosity, - ); - // this is required as we are restarting the time - stem_cell.last_division_t = 0f32; - } + if let NeutralMutations::UponDivisionAndBackground = moran.proliferation.neutral_mutation { + // this is important: we update all background mutations at this + // time such that all cells are all on the same page. + // Since background mutations are implemented at each division, and + // cells do not proliferate at the same rate, we need to correct and + // update the background mutations at the timepoint corresponding + // to sampling step, i.e. before saving. + if self.verbosity > 0 { + println!("updating the neutral background mutations for all cells"); + } + for stem_cell in moran.subclones.get_mut_cells() { + assign_background_mutations( + stem_cell, + self.time, + &self.distributions.neutral_poisson, + rng, + self.verbosity, + ); + // this is required as we are restarting the time + stem_cell.last_division_t = 0f32; } - MutateUponDivision::DivisionMutations(_) => {} } // restart the time moran.time = 0.; @@ -179,7 +176,7 @@ impl AdvanceStep for Exponential { // removes the cells from the clone with id `reaction.event`.** self.counter_divisions += 1; self.time += reaction.time; - self.proliferation.duplicate_cell_assign_mutations( + self.proliferation.proliferate( &mut self.subclones, self.time, reaction.event, @@ -211,7 +208,7 @@ pub struct Moran { pub filename: PathBuf, pub save_sfs_only: bool, pub save_population: bool, - pub proliferation: MutateUponDivision, + pub proliferation: Proliferation, } impl Default for Moran { @@ -231,7 +228,7 @@ impl Default for Moran { save_population: true, }, Distributions::default(), - MutateUponDivision::default(), + Proliferation::default(), 1, ) } @@ -246,7 +243,7 @@ impl Moran { time: f32, saving_options: SavingOptions, distributions: Distributions, - proliferation: MutateUponDivision, + proliferation: Proliferation, verbosity: u8, ) -> Moran { let hsc = Moran { @@ -388,28 +385,25 @@ impl Moran { if self.verbosity > 0 { println!("saving process at time {}", time); } - match self.proliferation { - MutateUponDivision::DivisionAndBackgroundMutations(_) => { - // this is important: we update all background mutations at this - // time such that all cells are all on the same page. - // Since background mutations are implemented at each division, and - // cells do not proliferate at the same rate, we need to correct and - // update the background mutations at the timepoint corresponding - // to sampling step, i.e. before saving. - if self.verbosity > 0 { - println!("updating the neutral background mutations for all cells"); - } - for stem_cell in self.subclones.get_mut_cells() { - assign_background_mutations( - stem_cell, - self.time, - &self.distributions.neutral_poisson, - rng, - self.verbosity, - ); - } + if let NeutralMutations::UponDivisionAndBackground = self.proliferation.neutral_mutation { + // this is important: we update all background mutations at this + // time such that all cells are all on the same page. + // Since background mutations are implemented at each division, and + // cells do not proliferate at the same rate, we need to correct and + // update the background mutations at the timepoint corresponding + // to sampling step, i.e. before saving. + if self.verbosity > 0 { + println!("updating the neutral background mutations for all cells"); + } + for stem_cell in self.subclones.get_mut_cells() { + assign_background_mutations( + stem_cell, + self.time, + &self.distributions.neutral_poisson, + rng, + self.verbosity, + ); } - MutateUponDivision::DivisionMutations(_) => {} } let population = self.subclones.get_cells_with_clones_idx(); if self.verbosity > 0 { @@ -486,7 +480,7 @@ impl AdvanceStep for Moran { // id `reaction.event`.** self.time += reaction.time; self.counter_divisions += 1; - self.proliferation.duplicate_cell_assign_mutations( + self.proliferation.proliferate( &mut self.subclones, self.time, reaction.event, @@ -516,7 +510,7 @@ mod tests { use rand::SeedableRng; use rand_chacha::ChaCha8Rng; - use crate::{genotype, proliferation::DivisionAndBackgroundMutationsProliferation}; + use crate::genotype; use super::*; @@ -562,9 +556,7 @@ mod tests { save_population: true, }, Distributions::new(u, 10f32, 1f32, 0), - MutateUponDivision::DivisionAndBackgroundMutations( - DivisionAndBackgroundMutationsProliferation, - ), + Proliferation::default(), 0, ) } @@ -650,9 +642,7 @@ mod tests { Exponential::new( SubClones::new(vec![StemCell::new(); cells.get() as usize], 3), Distributions::new(u, 10f32, 1f32, 0), - MutateUponDivision::DivisionAndBackgroundMutations( - DivisionAndBackgroundMutationsProliferation, - ), + Proliferation::default(), 0, ) } diff --git a/src/proliferation.rs b/src/proliferation.rs index 9e6b7c24..600fee78 100644 --- a/src/proliferation.rs +++ b/src/proliferation.rs @@ -1,68 +1,26 @@ -use enum_dispatch::enum_dispatch; use rand::Rng; +use rand_distr::{Bernoulli, Distribution}; use crate::{ - genotype::Variant, - stemcell::{assign_background_mutations, mutate, StemCell}, - subclone::{assign, proliferating_cell, CloneId, Distributions, SubClones}, + stemcell::{assign_background_mutations, assign_divisional_mutations}, + subclone::{assign_fit_mutations, proliferating_cell, CloneId, Distributions, SubClones}, }; -fn assign_mutations(stem_cell: &mut StemCell, mutations: Option>, verbosity: u8) { - if let Some(mutations) = mutations { - if verbosity > 2 { - println!("assigning {:#?} to cell {:#?}", mutations, stem_cell); - } - mutate(stem_cell, mutations); - } else if verbosity > 2 { - println!("no mutations to assign to cell {:#?}", stem_cell); - } +#[derive(Debug, Clone, Default)] +/// Stem cells can either self-renew by symmetric division or differentiate by +/// asymmetric division. +pub enum Division { + /// A stem cell proliferates and gives rise to a differentiate cell (which + /// we discard in our simulations) + Asymmetric(Bernoulli), + /// A stem cell proliferates and gives to a new stem cell with the same + /// genetic material as the mother cell plus some additional mutations + /// acquired upon division (fit divisions and neutral ones) + #[default] + Symmetric, } -fn proliferation( - stem_cell: StemCell, - subclones: &mut SubClones, - proliferating_subclone: CloneId, - distributions: &Distributions, - rng: &mut impl Rng, - verbosity: u8, -) { - if verbosity > 1 { - println!("cell {:#?} is dividing", stem_cell); - } - - let new_cell = stem_cell.clone(); - for mut cell in [new_cell, stem_cell] { - let division = distributions.neutral_poisson.new_muts_upon_division(rng); - if verbosity > 2 { - println!( - "assigning {} mutations upon cell division to cell {:#?}", - division.as_ref().unwrap_or(&vec![]).len(), - cell - ); - } - assign_mutations(&mut cell, division, verbosity); - assign( - subclones, - proliferating_subclone, - cell, - distributions, - rng, - verbosity, - ); - } -} - -#[derive(Clone, Debug)] -/// Simulate fit mutations, neutral background mutations and neutral mutations -/// upon division. -pub struct DivisionAndBackgroundMutationsProliferation; - -#[derive(Clone, Debug)] -/// Do not simulate neutral background mutations but only neutral mutations -/// and fit mutations upon division -pub struct DivisionMutationsProliferation; - -#[enum_dispatch] +#[derive(Clone, Debug, Default)] /// All the updates and changes in the system upon proliferation. /// /// The proliferation step is implemented as following: @@ -81,64 +39,21 @@ pub struct DivisionMutationsProliferation; /// old clone which `c` belonged to. /// /// 5. clone the proliferating cell `c` into `c1` and repeat step 3, 4 -/// with `c1`. -pub trait Proliferation { - fn duplicate_cell_assign_mutations( - &self, - subclones: &mut SubClones, - time: f32, - proliferating_subclone: CloneId, - distributions: &Distributions, - rng: &mut impl Rng, - verbosity: u8, - ); -} - -#[derive(Clone, Debug)] -#[enum_dispatch(Proliferation)] -/// Duplicate cells and deal with subclones. -/// Also specifies whether to simulate background neutral mutations or only -/// mutations upon divisions, see [`Proliferation`]. -pub enum MutateUponDivision { - DivisionAndBackgroundMutations(DivisionAndBackgroundMutationsProliferation), - DivisionMutations(DivisionMutationsProliferation), -} - -impl Default for MutateUponDivision { - fn default() -> Self { - MutateUponDivision::DivisionAndBackgroundMutations( - DivisionAndBackgroundMutationsProliferation, - ) - } +/// with `c1` only if we simulated a symmetric division, see [`Division`]. +pub struct Proliferation { + pub neutral_mutation: NeutralMutations, + pub division: Division, } -impl Proliferation for DivisionMutationsProliferation { - fn duplicate_cell_assign_mutations( - &self, - subclones: &mut SubClones, - time: f32, - proliferating_subclone: CloneId, - distributions: &Distributions, - rng: &mut impl Rng, - verbosity: u8, - ) { - let stem_cell = proliferating_cell(subclones, proliferating_subclone, verbosity, rng); - if verbosity > 1 { - println!("proliferation at time {}", time); +impl Proliferation { + pub fn new(neutral_mutation: NeutralMutations, division: Division) -> Self { + Proliferation { + neutral_mutation, + division, } - proliferation( - stem_cell, - subclones, - proliferating_subclone, - distributions, - rng, - verbosity, - ) } -} -impl Proliferation for DivisionAndBackgroundMutationsProliferation { - fn duplicate_cell_assign_mutations( + pub fn proliferate( &self, subclones: &mut SubClones, time: f32, @@ -150,21 +65,63 @@ impl Proliferation for DivisionAndBackgroundMutationsProliferation { let mut stem_cell = proliferating_cell(subclones, proliferating_subclone, verbosity, rng); if verbosity > 1 { println!("proliferation at time {}", time); + println!( + "cell with last division time {} is dividing", + stem_cell.last_division_t + ); + if verbosity > 2 { + println!("cell {:#?} is dividing", stem_cell); + } + } + if let NeutralMutations::UponDivisionAndBackground = self.neutral_mutation { + assign_background_mutations( + &mut stem_cell, + time, + &distributions.neutral_poisson, + rng, + verbosity, + ); + } + let cells = match self.division { + Division::Asymmetric(p) => { + // true means asymmetric + if p.sample(rng) { + if verbosity > 1 { + println!("asymmetric division"); + } + vec![stem_cell] + } else { + vec![stem_cell.clone(), stem_cell] + } + } + Division::Symmetric => vec![stem_cell.clone(), stem_cell], + }; + for mut stem_cell in cells { + assign_divisional_mutations( + &mut stem_cell, + &distributions.neutral_poisson, + rng, + verbosity, + ); + assign_fit_mutations( + subclones, + proliferating_subclone, + stem_cell, + distributions, + rng, + verbosity, + ); } - assign_background_mutations( - &mut stem_cell, - time, - &distributions.neutral_poisson, - rng, - verbosity, - ); - proliferation( - stem_cell, - subclones, - proliferating_subclone, - distributions, - rng, - verbosity, - ) } } + +#[derive(Clone, Debug, Default)] +/// Specifies the kind of neutral mutations to simulate. +pub enum NeutralMutations { + /// Neutral mutations upon division due to errors in the DNA duplication + UponDivision, + /// Neutral mutations upon division and background mutations (mutations + /// appearing during the lifetime of the cell) + #[default] + UponDivisionAndBackground, +} diff --git a/src/stemcell.rs b/src/stemcell.rs index 2dcf0159..a82b5b76 100644 --- a/src/stemcell.rs +++ b/src/stemcell.rs @@ -46,10 +46,27 @@ impl StemCell { } } -pub fn mutate(cell: &mut StemCell, mut mutations: Vec) { +fn mutate(cell: &mut StemCell, mut mutations: Vec) { cell.variants.append(&mut mutations); } +pub fn assign_divisional_mutations( + stem_cell: &mut StemCell, + neutral_poisson: &NeutralMutationPoisson, + rng: &mut impl Rng, + verbosity: u8, +) { + let mutations = neutral_poisson.new_muts_upon_division(rng); + if let Some(mutations) = mutations { + if verbosity > 2 { + println!("assigning {:#?} to cell {:#?}", mutations, stem_cell); + } + mutate(stem_cell, mutations); + } else if verbosity > 2 { + println!("no mutations to assign to cell {:#?}", stem_cell); + } +} + pub fn assign_background_mutations( stem_cell: &mut StemCell, time: f32, @@ -61,6 +78,9 @@ pub fn assign_background_mutations( //! current simulation time. //! //! This updates also the time of the last division for the cell. + if verbosity > 1 { + println!("assigning background mutations"); + } let interdivison_time = time - stem_cell.last_division_t; // 2. draw background mutations and assign them to `c` if let Some(background) = neutral_poisson.new_muts_background(interdivison_time, rng, verbosity) diff --git a/src/subclone.rs b/src/subclone.rs index 4d539035..e4227064 100644 --- a/src/subclone.rs +++ b/src/subclone.rs @@ -155,7 +155,7 @@ impl SubClone { } } -pub fn assign( +pub fn assign_fit_mutations( subclones: &mut SubClones, old_subclone_id: CloneId, cell: StemCell, @@ -501,7 +501,7 @@ mod tests { let cell2assign = StemCell::new(); let mut subclones = SubClones::new(cells, cells_present as usize + 1); - assign(&mut subclones, 0, cell2assign, &distr, &mut rng, 0); + assign_fit_mutations(&mut subclones, 0, cell2assign, &distr, &mut rng, 0); subclones.0[0].cell_count() as usize == before_assignment + 1 } @@ -515,7 +515,7 @@ mod tests { let cell2assign = StemCell::new(); let mut subclones = SubClones::new(cells, cells_present as usize + 1); - assign(&mut subclones, 0, cell2assign, &distr, &mut rng, 0); + assign_fit_mutations(&mut subclones, 0, cell2assign, &distr, &mut rng, 0); subclones.0[0].cell_count() as usize == before_assignment } @@ -528,7 +528,7 @@ mod tests { let mut subclones = SubClones::default(); let cell = StemCell::new(); - assign(&mut subclones, 0, cell, &distr, &mut rng, 0); + assign_fit_mutations(&mut subclones, 0, cell, &distr, &mut rng, 0); } #[should_panic]