v4.1.0

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
 

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "hsc"
-version = "4.0.0"
+version = "4.1.0"
 edition = "2021"
 
 [dependencies]

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))
             }

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);
             }
         }

src/main.rs

@@ -4,17 +4,15 @@ 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,
 };
 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,
             )
         };

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,15 +65,15 @@ pub struct Exponential {
     pub counter_divisions: usize,
     pub verbosity: u8,
     pub distributions: Distributions,
-    pub proliferation: MutateUponDivision,
+    pub proliferation: Proliferation,
     pub time: f32,
 }
 
 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<MAX_SUBCLONES> 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<MAX_SUBCLONES> 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,
         )
     }