A Bayesian method to trace cell lineages from CRISPR barcode data with linked target sites

We designed GABI to consider GESTALT barcode alignments as data for Bayesian phylogenetic inference. We enable this by re-parameterizing and implementing a phylogenetic likelihood (i.e. the probability of the sequence alignment given the tree) under the model of GESTALT barcode evolution described in [37] (called GAPML model), in BEAST 2 [20]. In the following, we briefly introduce the notations, definitions and assumptions that were used to formulate this model of GESTALT sequence evolution and underline key assumptions and approximations used for the likelihood calculation under that model.

GESTALT relies on a Cas9-based technology to generate stochastic mutations introduced at specific sites within a transgene (referred to as barcode) built into the founding cell of the population whose lineage is traced. Lineage data consist of an alignment of such modified transgenes, generated through sequencing of GESTALT barcodes after cell populations of interest undergo proliferation. The following are general notations used to describe any of the several GESTALT barcode designs introduced in [16] but also hold for all subsequent GESTALT applications and related technologies [39–44] and can be directly modelled by only adjusting the available parameters.

The GESTALT barcode is a DNA sequence that consists of an array of M (typically, up to 10) subsequences (typically each in the range of 20 bp) each containing a CRISPR/Cas9 target site. Each of these subsequences contains a protospacer adjacent motif (PAM) sequence necessary for Cas9 nuclease binding and a spacer sequence complementary to single-guide RNA (sgRNA) for specific recognition. In the presence of the Cas9 enzyme and a sgRNA complementary to the spacer sequence at a target site, mutations of GESTALT barcode sequences occur as the result of a two-step process: (i) upon sgRNA binding at a target site k, Cas9 first generates a double-strand cut at a position 3 bp upstream of PAM sequences, denoted as c(k), and (ii) this then triggers an endogenous repair mechanism that results in the introduction of insertions and deletions into the barcode at the targeted sites (see electronic supplementary material, appendix figure S1, upper half). Introducing an insertion or deletion at each of the M target sites is irreversible as it breaks complementarity with its corresponding sgRNA. Therefore, such sets of insertions and/or deletions (together referred to as 'indels') are the building blocks of characters in GESTALT sequence alignments. For such GESTALT sequencing data, [37] suggested a model of this evolution, implemented in GAPML. It contains parameters that specify target-specific cut rates, weights for the relative frequency at which double-cut events occur and factors controling the frequency of long and short deletions left and right of these cuts (see table 1 for notation). In electronic supplementary material, appendix B, we describe the precise model assumptions and how the likelihood is calculated under that model. In our implementation, we additionally introduce a clock rate r, quantifying the rate per time unit of mutations per target site (also see table 1). This distinguishes GABI from GAPML, where no such rate is assumed, which results in the length of any branch i, τi being inferred relative to total tree height. Consequently, under a transition rate matrix Q built under the GAPML model, transition probabilities are expressed as elements of the transition probability matrix P=eQτi. Instead, in GABI, we additionally provide scaling in absolute time with the clock rate, where a branch i accumulates indels at a rate of ri indels per target site per unit of time. In BEAST 2, these branch rates can then either be specified as being constant for all branches, corresponding to a strict molecular clock [45] assumption or allowed to vary across branches or clades, corresponding to a relaxed clock [46,47]. Under this implementation, transition probabilities are obtained as described in electronic supplementary material, appendix B, but substituting for the absolute branch length ti instead P=eQriti, which is the core of our novel parameterization.

We re-parameterized and implemented the phylogenetic likelihood and the sequence simulation engine for the GESTALT model from [37] in Java for the Bayesian phylogenetic inference platform BEAST 2, as a package called GABI. This implementation enables its interaction with dozens of tree priors, molecular clock models and a growing body of tutorials to facilitate the specification of complex analyses. Its integration within BEAST 2 also ensures interoperability with software that provides visualization tools to summarize Markov chain Monte Carlo (MCMC) posterior samples and assess convergence.

Several key changes were made from the original Python framework with regard to the tree inference: GABI considers exclusively bifurcating trees and thus does not implement the scheme aimed at resolving multifurcations by creating spines and does not restrict the tree search to maximum parsimony topologies. Additionally, an advancement made in GABI is the increased speed of the likelihood calculation by caching partial likelihoods and ancestral states for unchanged parts of the lineage tree, leading to up to a 2.5-fold increase over 100 iterations of the MCMC chain with respect to a naive implementation. Overall runtime can also be further optimized by the use of threading natively supported by BEAST 2.

As outlined above, the core novel component of GABI is its ability to specify and infer a molecular clock rate, r, which denotes the rate that converts branch lengths into calendar time. This can further be extended by specifying branch or clade-specific rates when deemed appropriate. Building upon this component, our implementation of a sequence simulation engine, based on the BEAST 2 package feast [48] combined with the GAPML simulation algorithm, is the first fully integrated simulator allowing to generate GESTALT alignments along time-scaled lineage trees, while specifying branch specific mutation rates. Beyond the assessment of proper implementation as shown here, it also provides a framework for more investigations of the properties of GESTALT data with respect to its ability to encode signal for both lineage and cellular growth patterns.

While models of indel evolution in the context of classic sequence or language evolution have been implemented for BEAST 2 [49], it does not natively support indel sequences characteristic of GESTALT barcode alignments and their metadata. To enable their use, we defined a new BEAST 2 datatype format, which encodes sequences as lists of indels each characterized by their start position in the barcode, the deletion length, the target sites affected by the deletion and the insert sequence. This implementation also enables the specification of any GESTALT barcode design of interest as metadata associated with input alignments for GABI, containing the unedited barcode sequence and positions of the cut site relative to the 3′ end of each target site.

Input data compatible with GABI is generated as follows: GESTALT alleles (in the format made available in [16]) are created from consensus GESTALT barcode sequences built using unique molecular identifiers. After being aligned to the reference, unedited barcode, nucleotide sequences are converted to alleles that consist of lists of observed indels. From this starting format, additional processing of ambiguous indels is performed using the pipeline provided in [37], merging indels affecting the same site and ensuring that indels span target sites correctly. Finally, these alleles are built into an alignment compatible with GABI in BEAST 2 (see electronic supplementary material, appendix figure S2 for an example). Additional data needed for computing the likelihood, such as the unedited barcode sequence and site annotations, are provided as metadata for BEAST 2 along with the alignment. All main functions and data structures that enable the calculation of the tree likelihood in GABI interact as shown in electronic supplementary material, appendix figure S3.

We first validate the correct implementation of the phylogenetic likelihood in GABI by comparing its output values with the previous Python implementation of GAPML. These tests are all integrated as unit tests for GABI in BEAST 2 and consist of calculations of the likelihood output for 13 datasets that cover all mutation types (unedited barcode, focal deletion or long deletions) and a range of phylogenetic structures (single branch, 3-leaf tree and 4-leaf tree). For each of them, we additionally record the runtime needed for a single, initial evaluation of the likelihood in both implementations, recorded on a computer running OSX on Intel Core i5 dual-core at 2.30 GHz.

We also validate the correct statistical properties of our implementation of the likelihood using a well-calibrated simulation study [50]. We first draw 100 sets of mutation parameters from distributions described in electronic supplementary material, appendix table S1. These distributions were centred around previously inferred maximum likelihood values on V7 GESTALT barcode data [37]. For simplified interpretation (but without loss of generality), we simulate mutations of an unedited barcode identical to V7 of the GESTALT design [16], shortened to 4 positions instead of 10 (described in electronic supplementary material, appendix table S2). For every set of mutation parameters, we simulate a dataset along one of 100 simulated trees chosen to have between 50 and 100 tips (see electronic supplementary material, appendix table S3 for parameters used for tree simulation). We chose to simulate alignments under a molecular clock rate centred around approximately 0.02 indels per site per hour, and arbitrarily assumed an experimental duration of 25 days. Together, this is consistent with a relatively low mutation rate to reflect a lower bound in recording capacity, corresponding to an average of about 0.5 indels per barcode site in each simulated GESTALT alignment. Using these simulated alignments as data, we infer the mutation parameters with the priors being equal to the distributions that were used for simulation, fixing the tree to its true value. Each dataset was used for inference using MCMC for 106 iterations and a burn-in of 30% of the iterations was discarded. Convergence was defined by reaching an effective sample size (ESS) >200 for all parameters. Coverage at a given credibility level α is calculated for each of the parameters as the proportion of times where the true parameter is contained within the inferred α % highest posterior density (HPD) interval over all simulated datasets.

We showcase GABI’s implementation with a phylodynamic analysis of the growth of dome-stage transgenic D. rerio embryos, from [16]. In this experiment, zebrafish embryos are harvested for sequencing of GESTALT barcodes 4.33 hours post-fertilization (hpf) (also see figure 1A). This is a well-characterized developmental stage, where embryos are known to have undergone between 12 and 14 rounds of cell division and no cell death [51]. This makes the dataset ideal as a first, proof-of-concept application of our framework for the study of cell growth dynamics using GESTALT barcode data. Based on the sequencing data, we aim to infer the time-scaled phylogenetic trees together with the division rates (birth rates) and sampling intensity, employing GABI together with the birth–death and coalescent framework in BEAST 2 [7,20,38].

The GESTALT barcode design used for this experiment is V7, consisting of 10 target sites of 23 bp each (see electronic supplementary material, appendix table S6 for the full, unedited barcode sequence). Editing of these barcodes is achieved by injection of Cas9 editing reagents into transgenic fertilized D. rerio eggs and is thought to taper starting after dome stage [16]. To model these editing dynamics, we assume a constant, strict molecular clock for the duration of the experiment. Prior distributions on the GESTALT editing parameters around mean values estimated previously on the same barcode design (see electronic supplementary material, appendix table S7 and §2a). We fix the mean indel length parameters in this example, which were previously inferred most consistently over different experiments using barcode V7 in [37]. This dataset consisted of 10 replicates with a median of 2000 sequences available per embryo. Thus, it constitutes an ideal case study for the first application of phylodynamics to the quantification of zebrafish embryonic development. To model embryonic growth and the harvesting process, we apply the two foundational models of population growth available for phylodynamic inference: the birth–death-sampling model [52], and the coalescent with exponential growth [38]. In the birth–death-sampling model, birth rate β corresponds to the rate of cell divisions, while death rate δ (in our applications, fixed to 0) corresponds to the rate at which cells undergo apoptosis, and ρ corresponds to the fraction of cells used in the analysis over all cells at the end of the experiment. Under this model, the expected population size after time t is N(t)=e(β−δ)t.

The coalescent with exponential growth is parameterized with a population growth ratealso corresponding here to the rate of cell divisions (assuming no apoptosis happens). Here, the model assumes a population sizeat time. Under both models, we report estimates of the cell division rate and the sampling proportion. For the coalescent analyses on an alignment ofcells, the sampling proportion is reported as the ratio(see electronic supplementary material, appendix tables S8 and S9 for detailed priors). g N t ( ) = e g t t s s N t ( )

To present different use cases of GABI, we designed three analysis pipelines for this dataset, described in figure 1B–D. In the first one, we infer the lineage tree topology, its branch lengths and growth parameters for three of the dome alignments individually (figure 1B). For this analysis, we limit the number of sequences per dataset to 20 sequences each to remain within a short runtime (in the range of hours per dataset). As these datasets are replicates of the same experiment, we additionally perform a second analysis where parameters of the GESTALT barcode, and cell population growth parameters are inferred over the same three sets of sequences (overall, 60 sequences) pooled together (figure 1C). The third analysis pipeline uses a larger subsample of 100 sequences for the first replicate, where the tree topology is fixed to a maximum-likelihood estimate obtained with GAPML (figure 1D). For this analysis, only branch lengths in the lineage tree are estimated jointly with model parameters. Each of these analyses was run for approximately 10⁷ iterations, and approximately 20% of steps were discarded as burn-in. Convergence is confirmed for parameter estimates that reach an ESS larger than 200 (and assessed using software Tracer [53]). All trees shown in this study are maximum clade credibility (MCC) trees constructed to summarize the posterior tree sample.

All datasets analysed were obtained from the Gene Expression Omnibus repository, with query GSE81713. Specifically, the three datasets analysed are the following: GSM2171788 (dome 1), GSM2171791 (dome 2) and GSM2171791 (dome 3), where the provided GESTALT input sequences were aligned and processed with the pipeline described in [16] (see section above for processing specific to GABI).

To show the correct implementation of the GESTALT phylogenetic likelihood computed in GABI, we first ensured that its output coincides with its previous implementation in Python. GABI implements a special case of the Python implementation of GAPML, where only bifurcating trees are considered. Correct calculation of transition probabilities alone is first assessed by calculating likelihoods along single branches, for all possible indel types. The full likelihood calculation, which recursively calculates and combines subtree likelihoods using Felsenstein’s pruning algorithm, is then tested on three small bifurcating trees with example sequences. Figure 2A shows the equality of the log-likelihood values calculated under both implementations. Additionally, we recorded the runtime needed to calculate these likelihood values under each implementation (see figure 2B) and found that GABI outperforms the original Python code, by up to an order of approximately 100 fold (which corresponds to up to seconds).

Next, we validated the correctness of the likelihood implementation within its Bayesian inference framework with a simulation-based calibration procedure. For all parameters, after ensuring convergence, consistency between the credibility level and the percentage of HPD intervals containing the true value is achieved, showing that the targeted posterior distribution is sampled by GABI (figure 2C). Additionally, we further investigated the inferential properties of the GABI framework at a credibility level of 95%. All parameters show the expected coverage at this level (within the 95% confidence interval (CI) of the binomial distribution with n = 100 trials and p=0.95, see electronic supplementary material, appendix table S4). Further quantification of the performance (root mean squared error, bias and relative HPD interval width) of the inference under these validation settings are also reported in electronic supplementary material, appendix table S5. Notably, these results are achieved with data generated under the assumption of a relatively low molecular clock rate. Therefore, we underline that the performance shown here represents a lower bound on the inferential performance of GABI as we studied its properties under a regime with low signal/little data. Higher mutation rates, leading to up to seven indels per barcode as observed in our example datasets, may provide more signal and improve GABI’s performance.

Estimates of the cell division rate and sampling proportion from dome-stage zebrafish GESTALT data under the three chosen analysis pipelines are shown in figure 3. Each compares estimates obtained under a constant birth–death model and the coalescent with exponential growth. While both model the growth of the studied population of cells as being exponential, their key difference is that the birth–death model assumes stochastic population dynamics, while the growth in the coalescent model is assumed to be deterministic. We consider both models to be potentially relevant in the context of this data. Overarching features of early embryonic development are thought to be largely deterministic, such as synchronicity until the tenth division [54,55]. However, many genetic and environmental factors, here including experimental conditions at play in generating replicates (e.g. growth conditions, handling and harvesting), are also likely to introduce some level of stochasticity in the growth of the cell populations under study [56]. As the comparative adequacy of these two models for the analysis of growing cell populations remains to be investigated (in the context of epidemiological data, see [57]), we consider estimates under each model here as equally capable of capturing dome-stage growth.

Individual division rate estimates under the birth–death-sampling model are consistent across all three replicates of the experiment, with median estimates around 2.2 h⁻¹, corresponding to cells until the dome stage dividing on average every 27 min (figure 3B). While overlapping, respective estimates of the division rate under the coalescent model suggest a faster division rate, with all three median division rates above 3.1 h⁻¹, corresponding to divisions occurring every 20 min or less on average. Our analysis infers a consistent sampling proportion for these three datasets and across both growth models used, where the alignments are estimated to contain between 0.1 and 0.2% of total cells in the embryos (figure 3F). Comparatively, median estimates of parameters of the GESTALT mutation process are estimated with slightly more variation across replicates than population parameters but do not suggest significant differences. For example, median molecular clock rate estimates range between 0.06 and 0.08 indels per site per hour across dome embryos (figure 4C). Likewise, while site-specific cut rates vary across replicates, sites 3, 8 and 10 are consistently inferred to be the least active sites. Together, these estimates show that the V7 barcode design used in this experiment carries reproducible properties across the three replicates (see electronic supplementary material, appendix figures S4 and S5, for all pipelines).

We also examine the shape of the tree topologies inferred for the three datasets (figure 4A), looking at the tree imbalance, which refers to the degree to which clade sizes, or groups of related cells, differ in the lineage tree. For a tree of size 20 under a birth–death model, the expected imbalance index [58], representing the sum of the differences in clade sizes across all internal nodes, is 38.5 (95% CI (36.4, 41.1)). Interestingly, while the MCC tree for replicate 1 is within this expected range (40), for replicates 2 and 3, the inferred MCC lineage trees have lower imbalance indices (respectively, 13 and 20). The MCC trees inferred under the coalescent model (see electronic supplementary material, appendix figure S6) show a similar pattern, where trees for the second and third replicate are more balanced than the first replicate, with indices of 38 and 23 (51 for the first replicate). Under the assumption of uniform sampling of the cell population, more balanced trees correspond to an underlying process where sub-populations of cells are relatively even. This may agree with some of the assumptions of the population process under study: zebrafish embryos do not experience cell death or apoptosis until later embryonic stages and early cell divisions characteristic of zebrafish development are known to be largely synchronous until at least 3 hpf [55]. Based on this knowledge, highly balanced trees are expected to be reconstructed from such an experiment. While neither the birth–death sampling nor the coalescent models specifically allow for synchronous divisions, the lower imbalance observed in the tree estimates for the second and third replicate may hint towards these alignments containing signal for such division events.

Our second analysis pipeline showcases how experimental replicates as analysed in our first example can be leveraged within a joint inference framework instead where cell division rates and mutation parameters are assumed to be the same across replicates (described in figure 1C). In BEAST 2, this is easily achieved by multiplying the individual likelihoods for each replicate in each step of the MCMC. In figure 3C,G, we show the division rate and sampling proportion estimated jointly over all three datasets, under the birth–death and coalescent models (see figure 4D for the clock rate). For the birth–death and coalescent models alike, the cell division rate inferred over this pooled data overlaps with individual estimates but is inferred with more confidence. For example, the estimated cell division rate from pooled data under the birth–death model shows a median 95%-level HPD interval width 26% narrower than the individual analyses. Pooled estimates under the exponential growth coalescent showcase a similar behaviour, where the division rate is estimated with lower uncertainty (10% narrower 95%-level HPD interval width). Interestingly, both analyses suggest slightly elevated division rates when compared to individual replicates, with median cell division rates elevated on average by 0.18 h⁻¹ (coalescent) and 0.19 h⁻¹ (birth–death) and suggest lower sampling proportions (with sampling proportions of 0.06 and 0.07% for the birth–death and the coalescent, respectively).

As a third example of an analysis pipeline using GABI for phylodynamic inference, we show parameter estimates obtained using as input a fixed GAPML tree topology for the first replicate of the experiment (dome 1) (figures 3D,H and 4E). In this scenario, only branch lengths are allowed to vary, where GABI provides estimates of the timing of individual cell divisions (see electronic supplementary material, appendix figure S7 for the inferred scaled lineage trees). This is consistent with many pipelines for phylodynamic inference in both fields of macroevolution and epidemiology (e.g. [59,60]), where pre-existing, or separately inferred, fixed tree topologies are used as input to infer population dynamic parameters. This is particularly relevant in the context of Bayesian inference, where the exploration of the exponentially growing lineage tree space quickly becomes intractable. Thus, using trees reconstructed with other, heuristic-based methods is often considered to enable the analysis of larger datasets. Here, we use a maximum-likelihood tree estimated under GAPML, for an alignment of 100 sequences as input, which is faster than GABI using their parsimony-guided tree search. The estimated division rate and sampling proportion under this pipeline are shown in figure 3D,H. Under the birth–death model, the division rate inferred remains mostly consistent with both previous analyses, but akin to the pooled analysis may suggest a slightly faster rate, 2.6 h⁻¹, with divisions occurring on average every 23 min. The division rate estimated for this dataset under the coalescent model is most shifted towards faster divisions, suggesting a median time until division of 16 min.

All three pipelines recover division rates that broadly agree with observed cell division lengths in zebrafish embryos: leading up to the mid-blastula transition (10th cycle post-fertilization), embryos undergo synchronous cell divisions every 18 min a most (from 6 min in early cell divisions to between 12.9 and 18 min starting from the sixth division [55]). After that stage, synchronicity is lost as individual cell cycles lengthen and are found to last for up to an hour. While we do not model a time-varying rate for the purpose of this analysis, our estimates, considering both posterior medians and their uncertainty, successfully capture these known average dynamics.