J-SPACE: a Julia package for the simulation of spatial models of cancer evolution and of sequencing experiments

Abstract

Background: The combined effects of biological variability and measurement-related errors on cancer sequencing data remain largely unexplored. However, the spatio-temporal simulation of multi-cellular systems provides a powerful instrument to address this issue. In particular, efficient algorithmic frameworks are needed to overcome the harsh trade-off between scalability and expressivity, so to allow one to simulate both realistic cancer evolution scenarios and the related sequencing experiments, which can then be used to benchmark downstream bioinformatics methods.

Result: We introduce a Julia package for SPAtial Cancer Evolution (J-SPACE), which allows one to model and simulate a broad set of experimental scenarios, phenomenological rules and sequencing settings.Specifically, J-SPACE simulates the spatial dynamics of cells as a continuous-time multi-type birth-death stochastic process on a arbitrary graph, employing different rules of interaction and an optimised Gillespie algorithm. The evolutionary dynamics of genomic alterations (single-nucleotide variants and indels) is simulated either under the Infinite Sites Assumption or several different substitution models, including one based on mutational signatures. After mimicking the spatial sampling of tumour cells, J-SPACE returns the related phylogenetic model, and allows one to generate synthetic reads from several Next-Generation Sequencing (NGS) platforms, via the ART read simulator. The results are finally returned in standard FASTA, FASTQ, SAM, ALN and Newick file formats.

Conclusion: J-SPACE is designed to efficiently simulate the heterogeneous behaviour of a large number of cancer cells and produces a rich set of outputs. Our framework is useful to investigate the emergent spatial dynamics of cancer subpopulations, as well as to assess the impact of incomplete sampling and of experiment-specific errors. Importantly, the output of J-SPACE is designed to allow the performance assessment of downstream bioinformatics pipelines processing NGS data. J-SPACE is freely available at: https://github.com/BIMIB-DISCo/J-Space.jl .

Keywords: Cancer Evolution; Next-generation sequencing; Spatial dynamics; Stochastic Simulation.

Fig. 1

The J-SPACE framework. Schematic representation of J-SPACE. A First, the algorithm simulates the spatial growth of the cells over an arbitrary graph. Then, J-SPACE simulates a spatial sampling (black circle) at a given time point. B J-SPACE reconstructs the phylogeny of the sampled cells (i.e., the leaves of the tree) and, given an ancestral genome, it generates the ground-truth sequence of the sampled cells using various substitution models. C A NGS experiment is simulated to return synthetic reads as outputs

Fig. 2

Phantom events and the reconstruction of phylogenetic trees. A Pictorial representation of the possible phantom events in a simulation with two different subpopulations. B Simplified scheme of the algorithm that generates the ground-truth phylogenetic tree from the list of birth events. First, the algorithm prunes the branches the leaves of which are not sampled (in red), then it removes the remaining edges that are not coalescent events

Fig. 3

Performances assessment. A The distribution of computational time in seconds to perform the simulation described in the text with respect to distinct sample size (over 50 simulation per configuration). B The distribution of computational time in seconds to generate the phylogenetic tree with respect to different sample size (over 50 simulation per configuration). C Distribution of computational time in seconds to generate the sequences for the phylogenetic trees above, with respect to distinct sample size (left) and genome length (right). In the top row, we present the results of the ISA-based model, in the bottom row we show the results of a finite-sites model (JC69) with indels (see the main text for further details)

Fig. 4

Analysis of cancer spatial dynamics and phylogenetic models. A The dynamics of the probability distribution of the number of cells is presented, divided by lattice dimensionality (2D or 3D). The dotted lines represent the expected values.B Box plots representing the distribution of the inferred steepness values of logistic growth are presented. C–D The distribution of the of the Sackin index and Beta-splitting statistic, evaluated on the trees divided by interaction rules and lattice dimensionality

Fig. 5

Variant calling with different mutational signatures. A An example dynamics of the number of cells for each subpopulation generated during the simulation. The bottom part of the panel presents the input driver mutational tree with the birth rate for each subpopulation. B At the top we present the phylogenetic tree generated by sampling 100 cells. We proceeded by simulating three different substitution models generated combinations of signatures SBS6 and SBS22 from the COSMIC database [71]. The difference between the three models consists in the time dynamics of the activation functions presented in this figure. C The count of the number of unique mutations simulated divided per class of substitution. The plot presents the result for the three different models. D The count of the number of unique mutations divided per class of substitution detected using the pipeline described in the main text