v0.26.0

CHANGELOG.md

@@ -1,5 +1,8 @@
 # Changelog
 
+## 0.26.0
+- Subsample the ecDNA distribution at the end of the simulation
+
 ## 0.25.0
 - Save based on the number of cells not time
 

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "ecdna-evo"
-version = "0.25.0"
+version = "0.26.0"
 edition = "2021"
 repository = "https://github.com/fraterenz/ecdna-evo"
 description = "Evolutionary models of extra-chromosomal DNA (ecDNA)"

src/clap_app.rs

@@ -90,21 +90,24 @@ pub struct Cli {
     /// Number of independent realisations to simulate of the same stochastic process
     runs: usize,
     #[arg(long, value_delimiter = ',', require_equals = true, num_args = 0..)]
+    /// Subsample the ecDNA distribution at the end of the simulation
+    subsamples: Option<Vec<u64>>,
+    #[arg(long, value_delimiter = ',', require_equals = true, num_args = 0..)]
     /// Number of cells that will trigger the saving of the ecDNA distribution.
-    snapshots: Option<Vec<usize>>,
+    snapshots: Option<Vec<u64>>,
     #[arg(short, long, action = clap::ArgAction::Count, conflicts_with = "debug", default_value_t = 0)]
     verbosity: u8,
 }
 
 fn build_snapshots(
     cells: NbIndividuals,
-    snapshots: Option<Vec<usize>>,
+    snapshots: Option<Vec<NbIndividuals>>,
 ) -> anyhow::Result<VecDeque<SnapshotCells>> {
     let mut snapshots = match snapshots {
         Some(s) => VecDeque::from_iter(
             s.into_iter().map(|cells| SnapshotCells { cells }),
         ),
-        None => VecDeque::from(build_snapshots_from_cells(11, cells as usize)),
+        None => VecDeque::from(build_snapshots_from_cells(11, cells)),
     };
 
     snapshots.make_contiguous();
@@ -117,9 +120,9 @@ fn build_snapshots(
 
 fn build_snapshots_from_cells(
     n_snapshots: usize,
-    cells: usize,
+    cells: NbIndividuals,
 ) -> Vec<SnapshotCells> {
-    let dx = cells / (n_snapshots - 1);
+    let dx = cells / (n_snapshots as NbIndividuals - 1);
     let mut x = vec![1; n_snapshots];
     for i in 1..n_snapshots - 1 {
         x[i] = x[i - 1] + dx;
@@ -155,6 +158,7 @@ impl Cli {
 
         let snapshots = build_snapshots(cells, cli.snapshots).unwrap();
         let segregation = cli.segregation.into();
+        let subsamples = cli.subsamples;
 
         // if both d0, d1 are either unset or equal to 0, pure birth,
         // else birthdeath
@@ -214,6 +218,7 @@ impl Cli {
             segregation,
             snapshots,
             growth: cli.growth,
+            subsamples,
             runs,
         };
         if verbosity > 0 {

src/lib.rs

@@ -9,12 +9,13 @@ pub mod segregation;
 use std::{collections::VecDeque, path::PathBuf};
 
 pub use ecdna_lib::{distribution, DNACopy};
+use sosa::NbIndividuals;
 
 #[derive(Debug, Clone)]
 /// Save the ecDNA distribution when the number of cells in the system is equal
 /// to `cells`
 pub struct SnapshotCells {
-    pub cells: usize,
+    pub cells: NbIndividuals,
 }
 
 pub struct SavingOptions {

src/main.rs

@@ -13,7 +13,7 @@ use indicatif::ParallelProgressIterator;
 use rand::SeedableRng;
 use rand_chacha::ChaCha8Rng;
 use rayon::prelude::{IntoParallelIterator, ParallelIterator};
-use sosa::{simulate, CurrentState, Options, ReactionRates};
+use sosa::{simulate, CurrentState, NbIndividuals, Options, ReactionRates};
 
 use crate::clap_app::{Cli, Parallel, ProcessType};
 
@@ -40,6 +40,7 @@ pub struct SimulationOptions {
     path2dir: PathBuf,
     options: Options,
     snapshots: VecDeque<SnapshotCells>,
+    subsamples: Option<Vec<NbIndividuals>>,
 }
 
 fn main() {
@@ -106,6 +107,20 @@ fn main() {
                 .expect(
                     "cannot save the ecDNA distribution at the end of the sim",
                 );
+                if let Some(samples) = app.subsamples.as_ref() {
+                    for nb_cells in samples.iter() {
+                        save(
+                    &process.path2dir,
+                    &process.filename,
+                    process.time,
+                    &process.get_ecdna_distribution().into_subsampled(*nb_cells, &mut rng),
+                    process.verbosity,
+                )
+                .expect(
+                    "cannot save the subsampled ecDNA distribution at the end of the sim",
+                );
+                    }
+                }
                 (
                     stop_reason,
                     [initial_state.population[0], initial_state.population[1]],
@@ -166,6 +181,20 @@ fn main() {
                 .expect(
                     "cannot save the ecDNA distribution at the end of the sim",
                 );
+                if let Some(samples) = app.subsamples.as_ref() {
+                    for nb_cells in samples.iter() {
+                        save(
+                    &process.path2dir,
+                    &process.filename,
+                    process.time,
+                    &process.get_ecdna_distribution().into_subsampled(*nb_cells, &mut rng),
+                    process.verbosity,
+                )
+                .expect(
+                    "cannot save the subsampled ecDNA distribution at the end of the sim",
+                );
+                    }
+                }
                 (
                     stop_reason,
                     [initial_state.population[0], initial_state.population[1]],

src/process.rs

@@ -105,6 +105,10 @@ where
     pub fn get_ecdna_distribution(&self) -> &EcDNADistribution {
         &self.distribution
     }
+
+    pub fn ecdna_distribution(self) -> EcDNADistribution {
+        self.distribution
+    }
 }
 
 impl<P: EcDNAProliferation, S: Segregate> AdvanceStep<2> for PureBirth<P, S> {
@@ -119,7 +123,7 @@ impl<P: EcDNAProliferation, S: Segregate> AdvanceStep<2> for PureBirth<P, S> {
             && self.snapshots.iter().any(|s| {
                 *self.distribution.get_nminus()
                     + self.distribution.compute_nplus()
-                    == s.cells as u64
+                    == s.cells
             })
         {
             let snapshot = self.snapshots.pop_front().unwrap();
@@ -246,6 +250,10 @@ where
     pub fn get_ecdna_distribution(&self) -> &EcDNADistribution {
         &self.distribution
     }
+
+    pub fn ecdna_distribution(self) -> EcDNADistribution {
+        self.distribution
+    }
 }
 
 impl<P: EcDNAProliferation, S: Segregate> AdvanceStep<4> for BirthDeath<P, S> {
@@ -260,7 +268,7 @@ impl<P: EcDNAProliferation, S: Segregate> AdvanceStep<4> for BirthDeath<P, S> {
             && self.snapshots.iter().any(|s| {
                 *self.distribution.get_nminus()
                     + self.distribution.compute_nplus()
-                    == s.cells as u64
+                    == s.cells
             })
         {
             let snapshot = self.snapshots.pop_front().unwrap();