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. For instance, if we want to simulate a tumor with 10 clones with a rather linear topology, we could use the following code:

n <- 10 
k <- 8

B_mat <- create_B(n, k)
B_mat
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]    1    0    0    0    0    0    0    0    1     1
#>  [2,]    1    1    1    1    1    1    1    1    1     1
#>  [3,]    1    0    1    0    0    0    1    1    1     1
#>  [4,]    1    0    1    1    1    1    1    1    1     1
#>  [5,]    1    0    1    0    1    0    1    1    1     1
#>  [6,]    1    0    1    0    1    1    1    1    1     1
#>  [7,]    1    0    0    0    0    0    1    1    1     1
#>  [8,]    1    0    0    0    0    0    0    1    1     1
#>  [9,]    0    0    0    0    0    0    0    0    1     1
#> [10,]    0    0    0    0    0    0    0    0    0     1

Once B is obtained, 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 clone proportions is managed internally by the functions .distribute_freqs and calc_clone_proportions. The .distribute_freqs function samples the Dirichlet distribution for the node corresponding to the designated clone and scales the resulting values to the original proportion of the parent clone, unless the node represents the root clone. The calc_clone_proportions function oversees the calculation of clone proportions throughout the tree T, calling .distribute_freqs at each multifurcation to distribute proportions accordingly.

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 figure below:

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 internally in the place_clones_space function, which currently does not allow modification of the blending parameters. The values in the mean_diffs vector and the Beta distribution parameters (α = 1 and β = 5) have been set to achieve a specific blending pattern, favoring smaller d values to prevent excessive sparsity in the resulting samples.

Sampling model

So far we have explained how the clones of a tumor are modeled 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, along with the proportions of the clones in the tumor and the tumor blend described in the tumor model section, is implemented in the create_U function. Within this sampling model, the function sets parameters internally, including the sampling sites (cutpoints_idx and cutpoints) within the Gaussian mixture model domain [μ_Gi − 3 ⋅ σ_Gi, μ_Gn + 3 ⋅ σ_Gn], and uses a user-defined number of sampled cells, n_cells, as the parameter in the multinomial distribution.

To simulate and sample clone proportions from a B matrix and construct the U matrix, we can use the following code. In this example, we assume neutral evolution to simulate clone proportions and generate 6 samples, each with 100 cells:

m <- 6  
selection <- "neutral"

U_mat <- create_U(B = B_mat, m = m, selection = selection, n_cells = 100)
U_mat
#>       clon1 clon2 clon3 clon4 clon5 clon6 clon7 clon8 clon9 clon10
#> -2.81  0.00  0.00  0.00  0.00  1.00  0.00  0.00  0.00  0.00   0.00
#> 0.23   0.00  0.00  0.00  0.00  0.99  0.00  0.00  0.00  0.01   0.00
#> 3.27   0.00  0.11  0.00  0.11  0.00  0.00  0.01  0.13  0.64   0.00
#> 6.31   0.18  0.07  0.00  0.00  0.00  0.08  0.25  0.39  0.00   0.03
#> 9.35   0.00  0.00  0.09  0.00  0.00  0.02  0.00  0.00  0.00   0.89
#> 12.4   0.00  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.00   0.98

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_noisy matrices. In GeRnika, Equation (1) is encoded in the create_F function.

The create_F function, in addition to the two matrices to be multiplied, includes the boolean parameter heterozygous. This parameter determines whether to adjust the clone proportions for heterozygous states. If heterozygous is set to TRUE, the clone proportions are halved to reflect heterozygous states. Instead, if heterozygous is set to FALSE, the clone proportions remain unchanged. The function can be used as follows:

F_mat <- create_F(U = U_mat, B = B_mat, heterozygous = FALSE)

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_noisy matrix, so that F_noisy ≠ U · B. Specifically, the model simulates the noise at the level of the sequencing reads and then recalculates the new f_n − 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_n − 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_n − 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_noisy 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_noisy 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 general instantiation create_instance function. Essentially, when this parameter is set to TRUE, the add_noise function is called, which is who actually manages the addition of noise to the matrix F_true. The create_instance function also includes the depth parameter, which quantifies the amount of noise introduced. As explained in the previous paragraph, noise increases as depth values decrease. 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. However, users may choose any other value that meets 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 from the distribution. This parameter can therefore be used to control the level of noise in the data. Increasing α will widen the range of values sampled from the distribution, resulting in some samples with much higher or lower depth than the mean sequencing depth, thereby increasing noise in the data. Conversely, decreasing α will reduce variability across samples, leading to more consistent sequencing depth values and, thus, less noisy data.

The α parameter is encoded as overdispersion within the functions in GeRnika. The create_instance function does not include this parameter, but it is available in the add_noise function. Users who wish to modify it will need to first create a noise-free F matrix (either by using the sequence of create_B, create_U, and create_F, or by extracting F_true from the result of calling create_instance). They can then cal add_noise on that F matrix, where the value of the overdispersion parameter can be adjusted as described above.

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.