Simulation of tumor clonal data (advanced users)

Tumor model

The tumor model generates a clonal tree T and an associated matrix B, together with the clone proportions c and tumor blend at the moment of sampling. Briefly, the clonal tree T that represents the development of a tumor with a set of mutations M so that |M| = n, is a rooted tree on an n-sized vertex set V_n = {v₁, …, v_n} and an n − 1-sized edge set E_T = {e_ij}, where v_i corresponds to the vertex or clone i and an edge e_ij between two vertices v_i and v_j represents a direct ancestral relationship.

In the first place, an n-sized tree T with a random topology is created in the following manner. First, the root node of T, ℛ(T), is set and a random mutation M_i ∈ M is assigned to it. Then, for each of the remaining M_j ∈ M − {M_i} mutations, a new node v_j is created; then the mutation M_j is assigned to it and the node is attached as a child of an existing node in T. In order to meet the ISA model, each of the newly added nodes inherits all the mutations present in its parent node.

The attachment of nodes to the tree is not uniformly at random; instead, the nodes in the T that is being built have different probabilities of being chosen as parent nodes for the new nodes, depending on the number of ascendants, 𝒜(v_i), they have. Mathematically speaking, ∀M_j ≠ ℛ(T), its parent node 𝒫(M_j) is sampled from V_n following a multinomial distribution:

where

k ∈ (0, +∞) is the topology parameter that determines if the topology is more branched, with a decreasing probability for increasing numbers of ascendants (k < 1), or more linear, with an increasing probability for increasing numbers of ascendants (k > 1). At the same time, setting k to 1 assigns equal probabilities to all the nodes and generates a completely random topology. The relationship between k and the topology can graphically be seen in the following figure:

Function that determines the frequency of selection of an existing node in the tree T as the parent node of new children. The function is dependant on the number of ascendants and the topology parameter k.

This procedure to obtain a tree T and one of the possible B matrices associated to it is implemented in the function create_B. This function has no parameters other than n and k described in the section Basic instantiation so we do not recommend changing it; we suggest to just adjust the value of k depending on the desired topology.

Once B is obtained, the proportions of the clones in the tumor c are simulated. We stress that these proportions c_i refer to the moment of the sampling because, due to the dynamic nature of tumor growth and evolution, these proportions need not be the same at any two time instants. It is also worth mentioning here that these proportions are not the ones that show in the U matrix, which are the clone proportions, but the ones that do not depend on any sampling.

These clone proportions c are calculated by sampling a Dirichlet distribution in each multifurcation of T. For instance, for a node v_i with children 𝒦(v_i) = {v_j, v_k}, we draw one sample {x_i, x_j, x_k} that represents the proportions of the parent clone and its two children from Dir(α_i, α_j, α_k), respectively. When this sampling is performed in an internal node of the tree, i.e., not in the root node, these proportions are scaled to the original proportion of the parent clone so that once all the multifurcations have been visited, the proportions of all the clones in T sum up to one.

The parameters of the Dirichlet distribution depend on the evolution model of the tumor. In GeRnika we consider two basic cases: positive selection-driven evolution or neutral evolution. According to positive selection-driven evolution, some mutations provide growth advantage whereas most of them do not, so the former clones get over other existing clones and take up the biggest part of the tumor. This results in tumors dominated by few clones, while the rest of clones are observed in minor proportions. Under neutral evolution, instead, there is not a significant number of mutations that provide fitness advantage and clones accumulate just out of tumor progression, so all clones are present in similar proportions

Based on this, the α parameters for the Dirichlet distribution for positive selection-driven evolution have been set to α = 0.3, and for the neutral evolution to α_p = 5 and α_c = 10 . We use different alpha values for parent and children nodes in the neutral evolution so that clones that result from multiple multifurcations do not end up with smaller proportions just because of their position in the topology, ruining the expected clone proportion distribution for this type of evolution model. These values have been selected empirically and their effect is depicted in the following figure, which shows how 5,000 random samples from the mentioned Dirichlet distributions (for the particular case of 3 dimensions, i.e., one parent and two children nodes) would be distributed on a 2-simplex. As can be seen, for positive selection (α = 0.3), the {x_i, x_j, x_k} values are pushed towards the corners of the simplex (α_i is < 1) in an uniform manner (all α_i are equal); in other words, the samples tend to be sparse; usually one component has a large value and the rest have values close to 0. Instead, when the neutral selection is adopted (α = [5, 10, 10]), the values in {x_i, x_j, x_k} concentrate close to the middle of the simplex (all α_i are > 1) but tend to deviate to those components with larger α. This means that samples {x_i, x_j, x_k} are less sparse, with larger values for x₂ and x₃ in this case, which represent the children nodes.

Ternary density plots of 5,000 samples drawn from two 3-dimensional Dirichlet distributions. The parameters of the Dirichlet distribution on the left are α = 0.3 and the distribution is used to represent positive-driven evolution. The distribution on the right has parameters α = [5, 10, 10] and it is used to represent neutral evolution.

In short, we have that the proportion of the clone v_i ∈ V_T|v_i ≠ ℛ(T) in the tumor at the moment of sampling, given by c_i, follows this distribution:

where

and

In turn, if v_i = ℛ(T):

The calculation of the clone proportions is implemented in the functions .distribute_freqs and calc_clone_proportions. The former function samples the Dirichlet distribution in the node corresponding to the clone designed by the clone_idx parameter and scales the resulting values to the original proportion of the parent clone, unless clone_idx is the root clone. The calc_clone_proportions function is responsible for calculating the clone proportions throughout the tree T, so it calls the .distribute_freqs function at each multifurcation.

Although we have seen that the parameter values used in the Dirichlet distributions are fairly good representatives of the underlying evolution models, it is true that, as we have said, the determination of these values has been done in an empirical way. If the user so wishes, they can try to vary these values by simply changing the values of the dirich_params variable in the .distribute_freqs function (the first row of the tibble in that variable corresponds to the clone clone_idx, whereas the second row corresponds to the children clones of clone_idx:

.distribute_freqs
#> function (B, clone_idx, clone_proportions, selection) 
#> {
#>     dirich_params <- dplyr::tibble(positive = rep(0.3, 2), neutral = c(5, 
#>         10))
#>     children_clone_idx <- get_children_idx(B = B, node_idx = clone_idx)
#>     if (!is.null(children_clone_idx)) {
#>         params <- dplyr::pull(dirich_params, eval(selection))
#>         dirich_values <- rdir(1, c(params[1], rep(params[2], 
#>             length(children_clone_idx))))
#>         norm_dirich_values <- dirich_values * clone_proportions[clone_idx]
#>         clone_proportions[c(clone_idx, children_clone_idx)] <- norm_dirich_values
#>     }
#>     return(clone_proportions)
#> }
#> <bytecode: 0x55af15a44e08>
#> <environment: namespace:GeRnika>

To finish with the tumor model, the tumor blend is simulated, i.e. how physically mixed are the clones between them. For doing so, we simplify the spatial distribution to one dimension and we model the tumor as a Gaussian mixture model of n components, where each component G_i represents a tumor clone and the mixture weights are given by c. The variance for all components is set to 1, while their mean value is a random variable. Specifically, a random clone is selected in the first place and the mean value of its component is set to 0. Then, the mean values of the n − 1 remaining components are calculated in a sequential manner by summing d units to the mean value of the immediately prior component. d ∈ {0, 0.1, …, 4} follows a Multinomial distribution so that:

For d = 0, the two clones are completely mixed; for d = 4, they are physically far from each other. This upper limit for d has been set empirically, after observing that with this value, the overlapping area between the two clones is already minimal, as shown in the following figure:

Overlapping area of pairs of Gaussian probability density functions with sd=1 and mean differences (d) ranging between 0 and 4. The functions represent 1-dimensional projections of tumor clone masses.

The tumor blend simulation is implemented in the place_clones_space function:

place_clones_space
#> function (B) 
#> {
#>     n <- nrow(B)
#>     . <- NULL
#>     max_sep <- 4
#>     mean_diffs <- seq(from = 0, to = max_sep, by = 0.1)
#>     clone_order <- sample(1:n, n)
#>     clone_mean_diff <- sample(x = mean_diffs, size = n - 1, prob = stats::dbeta(x = seq(0, 
#>         1, length.out = length(mean_diffs)), shape1 = 1, shape2 = 5), 
#>         replace = TRUE)
#>     clone_means <- c(0, cumsum(clone_mean_diff))
#>     clone_sd <- 1
#>     x <- seq(min(clone_means - 3 * clone_sd), max(clone_means + 
#>         3 * clone_sd), 0.01)
#>     y <- purrr::map(1:n, function(z) stats::dnorm(x, mean = clone_means[z], 
#>         sd = clone_sd)) %>% do.call(cbind, .) %>% dplyr::as_tibble(.name_repair = ~vctrs::vec_as_names(..., 
#>         repair = "unique", quiet = TRUE)) %>% magrittr::set_colnames(paste0("clon", 
#>         clone_order))
#>     spatial_coords <- y %>% dplyr::mutate(idx = x)
#>     return(list(spatial_coords = spatial_coords, x = x))
#> }
#> <bytecode: 0x55af14914138>
#> <environment: namespace:GeRnika>

The user can vary, e.g., increase or decrease, fix to a certain value, the range of values d can take to try other blending patterns. This can easily be done by changing the values of the vector mean_diffs, by modifying the minimum and maximum values of the sequence, or even the increment of the sequence.

Likewise, the user can try to vary the parameters in the Beta distribution with the same purpose. The current parameters α = 1 and β = 5 have been adjusted manually so that small d values are favoured and thus the samples obtained in the following steps do not become too sparse. If sparsity if desired, instead, they can, as said, increase the values in the mean_diffs variable or, equivalently, adjust the parameter values and give β a value larger than α in th clone_mean_diff variable. In any case, once the role of the Beta distribution regarding the tumor blend is clear, the user can make informed decisions on its parameter values.

Sampling model

So far we have explained how the clones of a tumor are modelled by the tumor model. However, in real practice there is no easy way of observing these tumor global properties; instead, we only have the information that samples or biopsies can provide. This means, among many other things, that the real clone proportions c can not be obtained. Instead, we can only get the clone proportions, which depend on certain sampling procedure; most of the time, the sampled clone proportions will not match the real ones. These sampled clone proportions are, in fact, the u_i. elements in the matrix U.

The sampling model in Gernika models the physical sampling of the tumor and, at the same time, allow us build the U matrix of the problem. This model actuates over the data simulated with the tumor model; specifically, it simulates some sampling performed in a grid-manner over the tumor Gaussian mixture model described in the previous subsection. Let G_i and G_n be the components with the lowest and largest mean values, respectively. Then, the domain is limited to [μ_Gi − 3 ⋅ σ_Gi, μ_Gn + 3 ⋅ σ_Gn] with σ_Gi = σ_Gn= 1, and it is divided into m+1 equal-sized bins by m cutpoints that represent the sampling sites. The 1^st and m^th cutpoints are always set to μ_Gi − 3 ⋅ σ_Gi + 20 and μ_Gn + 3 ⋅ σ_Gn − 20, and the remaining m-2 cutpoints are calculated accordingly.

The probability density of the components of the mixture model on those sites is normalized so that their sum sums up to 1 and the resulting values p_i are taken as the tumor clone composition values on the sampling spot i:

Finally, the effect of the cell count on the samples is considered so that the final tumor clone composition values in each sample, i.e., u_i., are modeled using a multinomial distribution:

The sampling model is implemented in the and create_U function:

create_U
#> function (B, clone_proportions, density_coords, m, x) 
#> {
#>     n_cells <- 100
#>     idx <- value <- proportion <- weights <- sampled_weight <- scaled_weight <- . <- NULL
#>     clone_order <- setdiff(colnames(density_coords), "idx")
#>     cutpoints_idx <- seq(20, length(x) - 20, length.out = m)
#>     cutpoints <- x[cutpoints_idx]
#>     U <- purrr::map(cutpoints, function(j) {
#>         reshape2::melt(density_coords, id = "idx") %>% dplyr::filter(idx == 
#>             j) %>% dplyr::mutate(weights = value/sum(value)) %>% 
#>             dplyr::left_join(clone_proportions, by = c(variable = "clone_idx")) %>% 
#>             dplyr::mutate(scaled_weight = weights * proportion/sum(weights * 
#>                 proportion)) %>% dplyr::mutate(sampled_weight = stats::rmultinom(n = 1, 
#>             size = n_cells, prob = scaled_weight)/n_cells) %>% 
#>             dplyr::select(sampled_weight)
#>     }) %>% do.call(cbind, .) %>% magrittr::set_rownames(clone_order) %>% 
#>         magrittr::set_colnames(as.character(cutpoints)) %>% t()
#>     U <- U[, order(as.numeric(gsub("clon", "", colnames(U))))]
#>     return(U)
#> }
#> <bytecode: 0x55af17563148>
#> <environment: namespace:GeRnika>

The user can again play with this function according to their needs. They could, for instance, vary the sampling sites; for doing so, they would just need to change the variables cutpoints_idx and/or cutpoints in the create_U function, as these variables contain the sampling sites in the Gaussian mixture model domain, which spans, it is worth recalling, from μ_Gi − 3 ⋅ σ_Gi to μ_Gn + 3 ⋅ σ_Gn.

It may also be interesting to adjust the n value in the multinomial distribution, which accounts for the number of cells sampled. Choosing a relatively low value for this parameter reduces the number of decimals in u_i., so that clones with frequencies several orders of magnitude lower than 1 can be modeled as absent in the sample, i.e., with composition values equal to 0. This is indeed way more realistic than truly observing them in so low frequencies. We set the value of n to 100 because it seemed a reasonable value and the resulting VAF values are in line with VAF values obtained in real samples; however should the user wish to try any other value, they could simply change the value of the variable n_cells in the first line of the create_U function.

Sequencing noise model

Up to this point, we have explained how the B and U matrices of an instance are simulated. In case we are simulating noise-free data, the instance is complete once Equation (1) has been applied on them in order to obtain both F_true and F matrices.

As a brief reminder, each element f_t − ij in F_true denotes the frequency or VAF of the mutation M_j on sample i or, in other words, the proportion of sequencing reads that harbour the mutation M_j on that sample. This also means that 1 - f_t − ij of the reads on that sample do not observe the mutation but the normal or reference nucleotide.

As the user may know, there exist numerous factors that alter the VAF value in an artificial manner so that this value deviates from the real ratio between the mutation and the reference nucleotides, including the noise that is introduced by the DNA sequencing process itself. The ways noise is introduced in this process are mainly two: through an incorrect nucleotide read of a fragment due to limitations of the sequencing instrument (e.g., a position that contains an A is read as a T) and through the production of a biased number of reads for a site that arises due to chemical reaction particularities or simply because not all the fragments get to be sequenced. This limitation can however be mitigated up to a certain level, as it has been shown that the higher the depth of coverage (μ_sd), i.e., the average number of reads that cover each position, the most accurate is the VAF value. In order to account for the effect of such noise in the data instances, we have designed a sequencing noise model. This model adds noise to the F_true matrix and builds a noisy F matrix, so that F ≠ U · B. Specifically, the model simulates the noise at the level of the sequencing reads and then recalculates the new f_ij values, as follows.

The sequencing depth r at the genomic position where M_j occurs in sample i is distributed as a negative binomial distribution so that:

The number of reads supporting the alternate allele r_ij^a is modeled by a binomial distribution:

The use of probability distributions for sampling the sequencing depth and the number of reads instead of using fixed values is the manner in which stochasticity is introduced in the read numbers and thus in the f_ij values.

Illumina mismatch errors are known to be around 0.1%. For simulating its effect on the VAF values, the number of reads r_ij^a′ that, despite originally supporting the alternate allele, contain a different allele as a result of a sequencing error are simulated as:

We also need to consider the situation where the reads contain the reference nucleotide but are read with the alternate allele as a result of this error. This is better understood with an example. Let us imagine that there is one genomic position where the normal cells have a T but in some cells there is a mutation and the T has changed to an A. For the normal cells, with probability 0.1% a sequencing error will occur and instead of a T, one of C, G or A will be read, all with equal chance. Thus, $\frac{0.001}{3}$% of the time, reads with the mutation of interest will arise from normal reads:

Thus, taking all these into consideration, the final noisy VAF values f_ijⁿ are simulated as:

In the following figure we have depicted for a collection of noisy instances, the density of the mean absolute error between the F and its corresponding noise-free F_true matrix. As it can be seen, the higher μ_sd, the lower the error we introduce to the F matrix. This is an expected behaviour since r_ij^a values follow a binomial distribution (Equation 14) where the number of trials is conditioned by μ_sd (Equation 13) and the event probability is the f_ij value.

Density of the mean absolute error in noisy F matrices for different μ_sd values that correspond to different noise levels.

The addition of sequencing noise to the data is controlled by the boolean noisy parameter in the create_instance function, and its quantity by the depth parameter in the same function. As explained in the previous paragraph, noise will increase with decreasing depth values. The default value for this parameter is 30 (30X); this value was chosen as most whole genome sequencing (WGS) data is sequenced at this depth. Again, the user can choose any other value that caters their needs.

The negative binomial distribution is a distribution commonly used to model sequencing depth. In such context, the mean of the distribution μ represents the expected number of reads for a position, whereas the overdispersion parameter α controls the degree of variability in the number of reads across different samples of the distribution. It is then easy to see that we can use this parameter to control for the noise introduced in the data. If we increase the value of α, the range of values of the samples drawn from the distribution will increase too, with some samples having much higher or lower depth values than the mean sequencing depth. In other words, increasing the value of α will increase the noise in the data. Conversely, if we decrease the value of α, there will be less variability across samples and the sequencing depth values will be more consistent across samples, i.e., data will be less noisy. The user can play with this value in the create_instance function, by changing the value of the overdispersion variable as we have just explained.

create_instance
#> function (n, m, k, selection, noisy = TRUE, depth = 30, seed = Sys.time()) 
#> {
#>     if (!is.character(selection)) {
#>         stop("\n selection must be a character")
#>     }
#>     if (selection != "neutral" & selection != "positive") {
#>         stop("\n selection must be neutral or positive")
#>     }
#>     if (depth < 1 || !identical(round(depth), depth)) {
#>         stop("\n depth must be a natural number greater than or equal to 1 ")
#>     }
#>     set.seed(seed)
#>     overdispersion <- 5
#>     B <- create_B(n, k)
#>     clone_proportions <- calc_clone_proportions(B, selection)
#>     clones_space <- place_clones_space(B)
#>     density_coords <- clones_space$spatial_coords
#>     domain <- clones_space$x
#>     U <- create_U(B = B, clone_proportions = clone_proportions, 
#>         density_coords = density_coords, m = m, x = domain)
#>     F_true <- calc_F(U = U, B = B, heterozygous = FALSE)
#>     if (noisy) {
#>         F_matrix <- add_noise(F_matrix = F_true, depth = depth, 
#>             overdispersion = overdispersion)
#>     }
#>     else {
#>         F_matrix <- F_true
#>     }
#>     clone_names <- purrr::map(1:n, function(x) paste0("clone", 
#>         x))
#>     mutation_names <- purrr::map(1:n, function(x) paste0("mut", 
#>         x))
#>     sample_names <- purrr::map(1:m, function(x) paste0("sample", 
#>         x))
#>     rownames(F_matrix) <- sample_names
#>     colnames(F_matrix) <- mutation_names
#>     rownames(B) <- clone_names
#>     colnames(B) <- mutation_names
#>     rownames(U) <- sample_names
#>     colnames(U) <- clone_names
#>     rownames(F_true) <- sample_names
#>     colnames(F_true) <- mutation_names
#>     return(list(F = F_matrix, B = B, U = U, F_true = F_true))
#> }
#> <bytecode: 0x55af161a0fb0>
#> <environment: namespace:GeRnika>

We would like to make two final notes regarding the sequencing noise model. The first one is that, even if our simulated data procedure follows the ISA, it is possible that due to the process of adding noise, resulting data breaks this assumption. The second one is a reminder that the sequencing noise model is an optional model only to be used if noisy data is to be generated.