v3.0.5

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

Cargo.toml

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

src/clap_app.rs

@@ -60,10 +60,15 @@ fn fitness_in_range(s: &str) -> Result<f32, String> {
 #[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<f32>,
     /// 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<f32>,
+    pub mu_division_exp: Option<f32>,
     /// 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(),

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<u16, u64>);
 
 impl MutationalBurden {
@@ -175,6 +179,30 @@ impl MutationalBurden {
         Ok(MutationalBurden(burden))
     }
 
+    pub fn from_moran(moran: &Moran, verbosity: u8) -> anyhow::Result<Self> {
+        //! 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<Self> {
+        //! 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 =

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

src/process.rs

@@ -495,3 +495,201 @@ impl AdvanceStep<MAX_SUBCLONES> 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::<u16>() < burden_after.0.keys().sum::<u16>());
+    }
+
+    #[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::<u16>() < burden_after.0.keys().sum::<u16>());
+    }
+    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::<u16>() < burden_after.0.keys().sum::<u16>());
+    }
+
+    #[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::<u16>() < burden_after.0.keys().sum::<u16>());
+    }
+}

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