 Research article
 Open Access
 Published:
SeqClone: sequential Monte Carlo based inference of tumor subclones
BMC Bioinformatics volume 20, Article number: 6 (2019)
Abstract
Background
Tumor samples are heterogeneous. They consist of varying cell populations or subclones and each subclone is characterized with a distinct single nucleotide variant (SNV) profile. This explains the source of genetic heterogeneity observed in tumor sequencing data. To make precise prognosis and design effective therapy for cancer, ascertaining the subclonal composition of a tumor is of great importance.
Results
In this paper, we propose a statespace formulation of the feature allocation model. This model is interpreted as the blind deconvolution of the expected variant allele fractions (VAFs). VAFs are deconvolved into a binary matrix of genotypes and a matrix of genotype proportions in the samples. Specifically, we consider a sequential construction of the genotype matrix which we model by Indian buffet process (IBP). We describe an efficient sequential Monte Carlo (SMC) algorithm, SeqClone, that jointly estimates the genotypes of subclones and their proportions in the samples. When compared to other methods for resolving tumor heterogeneity, SeqClone provides comparable and sometimes, better estimates of model parameters. By design, SeqClone conveniently handles any number of probed SNVs in the samples. In particular, we can analyze VAFs from newly probed SNVs to improve existing estimates, an attribute not present in existing solutions.
Conclusions
We show that the SMC algorithm for deconvolving VAFs from tumor sequencing data is a robust and promising alternative for explaining the observed genetic heterogeneity in tumor samples.
Background
Tumor samples that are obtained temporally or spatially from a cancer patient are heterogeneous in nature [1, 2]. These samples contain genetically diverse subpopulation of cells often referred to as tumor subclones [1, 3, 4]. Each subclone harbors a distinct mutational profile that uniquely characterizes the genome of the cells in that particular subclone [5–7]. Mutational and evolutionary processes that drive tumor progression are partly responsible for the observed genetic differences that distinguish these subclones. For instance, somatic variations among the subclones are as a result of mutations that are acquired by chance in the cell during tumor progression [8, 9].
The advancements in highthroughput sequencing technologies over the last decade [1, 10] have put a searchlight on studies that are related to tumor heterogeneity. For instance, some methods concentrate on probing individual cell using fluorescent markers [11, 12] while others employ single cell sequencing [13–16]. However, these approaches have their downsides. As an example, the use of single cell sequencing to probe large number of cells remains too expensive. On the other hand, methods like whole genome sequencing (WGS) and whole exome sequencing (WES) of tumor samples allow for proper and adequate quantification of somatic mutations in the cells [17].
One way to resolve tumor heterogeneity is to computationally characterize and identify the tumor subclones in the samples, employing the datasets from WGS and WES. Generally, computional approach at resolving tumor heterogeneity is a very challenging task [18]. It involves an estimation of the distinct single nucleotide variant (SNV) profiles/genotypes and their respective proportions in the samples. The result from such task assists in the design of effective therapy in combating cancer, aids correct cancer prognosis [19] and minimizes chemotherapy resistance [20].
In the literature, various computational methods have been proposed to resolve tumor heterogeneity [18]. Most prominent among these methods model the SNV profiles/genotypes of subclones with a binary matrix. Each row of the genotype matrix corresponds to a locus/SNV and each column represents the SNV profile of a subclone. Further, computational approach can be viewed as either an indirect or a direct estimation problem, depending on how the genotype matrix is obtained. In the former, genotypes of subclones in the tumor samples are not directly inferred. Rather mutations with similar cellular prevalence are first grouped as mutation clusters. As a result, further analyses are often required in order to obtain the genotypes/SNV profiles of tumor subclones in the samples [21–25].
The direct approach employs the feature allocation model for the decomposition of the observed variant allele fractions (VAFs) into matrices of genotypes (Z) and proportions (W) [26–29]. In addition to the VAF dataset, some methods include copy number information in the analysis of tumor heterogeneity [24]. These methods simultaneously model the copy number variation and SNV datasets. A host of methods under the direct approach assume a fixed number of subclones, and model the genotypes of subclones with a binary matrix. Each column of the matrix corresponds to the SNV profile of a subclone: 0 and 1 denoting the absence and presence of a particular SNV in a subclone [26, 27]. However, in reality, the exact number of subclones is not known prior to the analysis of the samples. To estimate model parameters of the feature allocation model, [27] proposed an expectationmaximization (EM) algorithm [30] that returns point estimates of model parameters. Markov chain Monte Carlo (MCMC) [31, 32], which has been the gold standard algorithm in the literature [21, 24, 26, 28], returns point estimates and variabilities of model parameters. As noted in [26, 28], when the number of SNVs is large, MCMC algorithm is plagued with computational issues. With EM and MCMC algorithms, whenever more VAFs are available from newly called SNV(s), there is no provision for improvement of the existing parameter estimates with the new datasets.
In this paper, we propose a statespace formulation of the feature allocation modeling framework. Our work also describes a sequential Monte Carlo (SMC) algorithm [33, 34] for inferring all the unknown model parameters that explain tumor heterogeneity. These parameters include the binary matrix of genotypes and the proportions of tumor subclones in the samples. In particular, our statespace formulation considers the sequential construction of the binary genotype matrix by making use of Indian buffet process (IBP) [35–37]. IBP describes the prior distribution of a binary matrix with a fixed number of rows and an unknown number of columns. Other parameters of the feature allocation model, including the proportions of tumor subclones, are considered as the parameters of our statespace model. The observed VAF, which is the input data, is processed rowwise: this enables scalability to any number of rows. In the SMC framework, observed measurements are processed one at a time. At every instance of time, the posterior probability density function (PDF) of the state at that time is computed via approximation [38–41]. With extensive simulation, we compare SeqClone with other computational methods for resolving tumor heterogeneity. Overall, in terms of accuracy of the estimates of model parameters, SeqClone demonstrates comparable and sometimes superior performance to other methods.
The remainder of this paper is organized as follows. In the “Results” section, we investigate the performance of SeqClone, using simulated datasets and chronic lymphocytic leukemia (CLL) datasets, the real tumor samples obtained from three patients in [42]. In the “Discussion” section, we discuss the results obtained from the proposed algorithm. “Conclusions” section concludes the paper. Finally, the “Method” section details the description of system model and problem formulation.
Results
In this section, we report the performance of the proposed algorithm using simulated and real tumor datasets. We compared model estimates, matrices of genotypes and proportions, from the proposed algorithm to those obtained from other similar algorithms. In real tumor datasets, similar to the manual approach considered in [27], we hypothesized phylogenetic trees from the estimated matrix of genotypes. Particularly, we assumed that the set of mutations that are grouped together in a tumor subclone comprises of: all the mutations that belong to its ancestors on the tree and the mutations on the edge that connect the subclone to its parent subclone. With this simple rule, we were able to construct the possible phylogenetic trees that are consistent with the estimated matrix of genotypes. For the simulation experiments, we employed a reverse of the above rule to generate binary genotype matrices from phylogenetic trees. Finally, we compared the runtimes of the different algorithms for subclone inference.
Simulated datasets
We generated datasets for average sequencing depth r∈{50,200,1000} per locus, number of tumor subclones C∈{3,4,5}, number of tumor samples S∈{3,4,...,10} and number of genomic loci T∈{20,40,60,80,100,5000}. For a given number of tumor subclones C and number of genomic loci T, we simulated a phylogenetic tree from where the genotype matrix Z is obtained. For the phylogenetic tree simulation, we grouped the T mutations into C subclones uniformly at random. The mutations in each subclone are assumed to first appear in that particular subclone on the tree. One of the subclones is randomly selected as the root node and the rest C−1 subclones are iteratively connected to the tree. Specifically, an unattached subclone (child) and a parent subclone on the tree are randomly selected. The child subclone is attached to the parent subclone and the new set of mutations in the child subclone is a union of the mutations in the parent and the mutations in the child subclone. The mutational profiles of the subclones are the columns of the genotype matrix Z.
Given the genotype matrix, along with specific values of r and S, we generated the input data to the proposed algorithm, i.e., the matrices of variant count Y and total count V. We generated each entry of V, i.e., v_{ts} from Pois(r). We generated each entry of Y, i.e., y_{ts} as follows: sampled each column of the proportion matrix W independently from Dir([a_{0},a_{1},...,a_{4}]) (a_{0}=0.2 and a_{c}, c∈{1,...,4} randomly chosen from the set {2,4,5,6,7,8}), defined p=0.02, computed p_{ts} following (2) in the “Method” section, and sampled y_{ts} from binomial(v_{ts},p_{ts}).
The proposed algorithm, Clomial [27], BayClone [28], and Cloe [26] were run on the simulated datasets. We defined the following metrics to quantify the estimation strength of the algorithms: genotype error (e_{Z}), proportion error (e_{W}) and success probabilities error (\(e_{p_{ts}}\)). Mathematically, these errors are defined as
and
The problem of estimating genotype matrix and proportions matrix is a blind decomposition problem. This implies that after the analysis, we are unaware of the columns of the estimated genotype matrix that correspond to the columns of the true genotype matrix. We resolved this by computing the genotype error with every permutation of the columns of the estimated genotype matrix. We selected the permutation that resulted in the smallest error and we used the selected genotype in computing the other error values. All experiments were performed on Intel(R) Xeon(R) CPU @ 3.5GHz and a 24GB of RAM running a 64bit Windows 7.
In Tables 1, 2, 3 and 4 and Figs. 1, 2, 3, 4, 5, 6 and 7, we present the results obtained from analyses of simulated datasets. To compare the methods, we generated 20 datasets for every combination of number of genomic loci T, number of tumor subclones C, number of tumor samples S and average sequencing depth r. We computed the average and standard deviation of genotype error e_{Z} and proportion error e_{W} over all the 20 datasets. In Table 1, we present the average and standard deviation (in round parentheses) of the genotype and the proportion errors for all the methods when the number of loci T=20, number of subclones C=3, number of samples S=5 and average sequencing depth r∈{50,200,1000}. We excluded success probabilities error (\(e_{p_{ts}}\)) because not all the algorithms return an estimate of p in (2) (“Method” section). Similarly, in Table 2, we show, for all the methods, the average and the standard deviation of genotype and proportion errors when T=100, C=3, S=5 and r∈{50,200,1000}. The proposed algorithm demonstrates a comparable and sometimes, superior performance in terms of the accuracy of the estimated genotype and proportion matrices. It should be noted that, for BayClone, the ones in the true binary genotype matrices were changed to 0.5 before the simulation and the entries of the estimated genotype matrices greater than 0 were changed to ones before computing the errors.
In Figs. 1, 2, 3, 4, 5 and 6, for SeqClone, we present the errorbar plots for the average and standard deviation over 20 datasets for different combinations of the number of loci, sample size, number of subclones and average sequencing depth. The standard deviation is the vertical line above and below the average value in the errorbar plots. Figures 1, 2 and 3 show how the errors vary across different sample sizes for different subclones. There is an improvement, for all the subclones, in the estimates of all model parameters when the number of tumor samples increases. Similarly, in Figs. 4, 5 and 6, estimates of model parameters improves when the average sequencing depth increases. In the first row in Table 3, we present, for SeqClone, the result of the permutations of rows of the input data. For the dataset with T = 100,C = 3,r = 1000 and S=5, we ran SeqClone with randomly selected 100 permutations of the rows of the input data matrices and we computed the average and standard deviation of the errors (row one in Table 3). In row two in Table 3, we present results for higher number of genomic loci. In particular, we present the average and standard deviation of errors over 20 runs for the datasets with T=5000, C=3, r=1000 and S=5.
Lastly, we present the runtimes and memory consumption for all the methods when performing a section of the experiments in Table 1. For the proposed algorithm, we ran the algorithm 20 times with 1000 particles. For the MCMCbased algorithms (Cloe and BayClone), we ran 30,000 chains. For Clomial, we ran 2000 iterations. The runtimes for all the methods on a 3.5Ghz Intel 8 cores running MATLAB and the associated genotype and proportion errors for the dataset from T=20, C=3, r=1000, and S=5 are in Table 4. In addition, for this particular dataset, we report the estimated sample mean and sample standard deviation of the relative frequency of variant reads that are produced as error (parameter p in “Method” section) from SeqClone and BayClone. For SeqClone, the mean is 0.019 and the standard deviation is 0.0012. Likewise, for BayClone, the mean is 0.022 and the standard deviation is 0.0011. In Fig. 7, we present the memory consumption by all the algorithms for different genomic loci (T∈{20,40,60,80,100}). In general, Clomial is the most memory efficient of all the algorithms. However, SeqClone consumes lesser memory when compared to BayClone and Cloe.
Real biological tumor samples
Next, we present the results obtained from applying the proposed algorithm to real biological tumor datasets. Particularly, we analyzed the datasets of three patients with Bcell CLL namely: CLL077, CLL006, and CLL003 [42]. Complete datasets and the data preprocessing steps are in [42]. In Additional file 1, we include the analysis results with Clomial, BayClone and Cloe.
CLL077:
Here, we present the results obtained from analyzing the dataset from patient CLL077 with SeqClone. This dataset had 16 distinct loci probed for tumor heterogeneity. These are shown in the first row in Table 5. We present our analysis results in the main paper, and the estimates for other methods in Additional file 1. In concordance with other methods, SeqClone estimated 4 subclones as shown in Table 5, and also produced SNV profiles that are similar to those obtained from the three other methods. Also, the proportions of tumor subclones exhibit similar trend in all the 5 tumor samples across various methods. For instance, the abundance of sublone 1 in sample ‘a’ in Clomial, BayClone, Cloe and SeqClone are 0.27,0.21,0.16 and 0.27, respectively. This trend continues in all other samples except in sample ‘e’ where Clomial deviates from this normal trend, i.e., Clomial, BayClone, Cloe and SeqClone are 0.43,0.07,0.03, and 0.16, respectively. On this dataset, SeqClone produced a consistent result with other methods in estimating the SNV profiles of subclones and their proportions in all the samples (Tables 5, 6 and in Additional file 1). The constructed phylogenetic tree from the SNV profiles for CLL077 is shown Fig. 8.
CLL006:
This dataset comprises of 11 genomic loci. These are shown in the first row in Table 7. We analyzed the dataset with SeqClone, and the estimates of genotype and proportions matrices are in Tables 7 and 8. The constructed phylogenetic tree is shown Fig. 9. SeqClone and BayClone estimated 5 distinct subclones, Clomial had 4 subclones and Cloe recovered 6 subclones. Details of the estimates from Clomial, BayClone and Cloe are in Additional file 1.
CLL003:
The dataset from patient CLL003 has 20 distinct genomic loci. This is shown in the first row in Table 9. In this dataset, Clomial and Cloe produced 2 distinct subclones with considerably high proportions in the samples and 2 others with very small proportions across all samples. SeqClone and BayClone estimated the first 2 major subclones that dominate the 5 samples with proportions shown in Table 10 (and Additional file 1: Table S6). The constructed phylogenetic tree for CLL003 is shown in Fig. 10.
Finally, we investigated the behavior of the algorithms in terms of runtime and memory consumption, when applied to simulated and real datasets of similar size: T=20 and S=5. We present the results in Table 11. Runtimes (without parentheses) are in minutes and the consumed memory (in round parentheses) are in MB.
Discussion
Tumor heterogeneity describes a situation where bulk tumor samples have numerous subpopulations of cancer cells and each subpopulation has unique features that distinguish it from other subpopulations in the samples. It has been recognized as the major cause of relapse in cancer patients. One way to resolve tumor heterogeneity is by deconvolving the VAFs data from the nextgeneration sequencing to the genotypes and the proportions of subpopulations of cancer cells in the samples. In this paper, to resolve tumor heterogeneity, we interpreted the VAFs data using the feature allocation model [27, 28].
We developed the feature allocation model into a statespace framework so that VAFs with large number of genomic loci can be adequately modeled. We proposed a sequential algorithm, SeqClone, to infer all the parameters of our statestate model. The inferred parameters, which describe tumor heterogeneity, include: the genotypes of all the genomic loci in every subpopulation and their respective proportions in the tumor samples. With the statespace modeling framework and the sequential algorithm, computational problem that is often encountered by other methods for interpreting tumor heterogneity in the presence of large genomic loci is eliminated [26, 28]. It should be noted that, in this work, like some previous methods [27, 43], only somatic SNVs/mutations are modeled and we assume that these mutations are unaffected by copy number aberrations or rearrangements in the cancer genome. With this modeling assumption, extreme care must be taken when using SeqClone to interpret tumor heterogeneity.
In the “Results” section, we presented the results from running SeqClone and three other algorithms: Clomial, BayClone and Cloe, on simulated and real cancer datasets. For the simulation experiments, we generated several simulated datasets and compared the results from all the algorithms. SeqClone produced comparable, and sometimes better performance in the estimation of model parameters. On the real cancer datasets ([42]), SeqClone produced satisfying results that are comparable to other methods.
Also, because of the sequential nature by which the VAFs are processed by SeqClone, VAFs from previously unprobed genomic loci can be analyzed to improve the existing results, a feature that is absent in other algorithms.
Conclusions
Finally, we have demonstrated the efficacy of sequential Monte Carlo algorithm in the analysis of VAFs datasets that are obtained from heterogeneous tumor samples. The proposed method does not assume that the number of subclones is known/fixed prior to analysis and this allows the ‘correct’ number of subclones to be estimated from the tumor samples. Also, because of the sequential nature by which the proposed algorithm handles the VAFs datasets, the analysis can easily be scaled to a very large dataset. In addition, the current framework can be extended to a more general case that involves the estimation of mutation and the copy number profiles of the tumor subclones that are present in the tumor samples.
Method
System model and problem formulation
Before going to the details of our modeling approach, we define all the mathematical notations that are used in this paper. p(·) denotes a PDF, p(··) denotes a conditional PDF, P(·) denotes a probability mass function (PMF) and P(··) denotes a conditional PMF. Likewise, binomial(n,p) denotes a binomial distribution with n exact number of trials and p probability of success at each trial. Bern(p) denotes a Bernoulli distribution with success probability p and \(\mathcal {N}\left (\mu,\sigma ^{2}\right)\) denotes a univariate Gaussian distribution with mean μ and variance σ^{2}. Also, gamma(α_{0},β_{0}) denotes a gamma distribution (α_{0} is the shape parameter and β_{0} is the rate parameter) and beta(α_{1},β_{1}) denotes a beta distribution where α_{1} and β_{1} are the shape parameters. Pois(λ) denotes a Poisson distribution with mean parameter λ and Dir(α) denotes a Dirichlet distribution with a vector of concentration parameters α. Lastly, Γ(·) denotes the gamma function and \(\hat {x}\) denotes the estimate of variable x.
Two important quantities that are obtained from WGS and WES of tumor samples are the variant count and total count at each of the probed genomic locus. We denote the matrix of variant count by Y and the matrix of total count by V. Each of the matrices has a dimension T×S, where T is the number of genomic loci/SNVs and S is the total number of tumor samples. We denote the number of reads that bear the variant count at locus t in sample s as y_{ts}. Likewise, we denote the total number of reads at locus t in sample s as v_{ts}. In our formulation, we assume that the genomic loci are unaffected by copy number aberrations or rearrangement of the cancer genome [27, 43]. We employ the binomial sampling model [27, 28] in modeling the input data matrices, given as
p_{ts}, t=1,...,T, s=1,...,S are the success probabilities defined as [28]
where z_{tc}, a binary variable, represents the two possible states of an allelic genotype at locus t in subclone c and C represents the number of tumor subclones, an unknown variable. Under this framework, if z_{tc}=1, it implies that locus t in subclone c has reads that bear variant sequence. Likewise, if z_{tc}=0, there are no reads that bear variant sequence at that locus. We assume that if a mutation is present in a particular subclone, then at that genomic locus, the subclone is heterozygous with copy number equal to one.
The term \({\sum }_{c=1}^{C} z_{tc}w_{cs}\) in (2) defines p_{ts} as a weighted sum of effects of an unknown number of subclones in the tumor samples. Also, effects of experimental and data processing noises are captured by w_{0s}p in (2). In particular, p is the relative frequency of variant reads that are generated as a result of error during upstream data analysis [28]. For t=1,...,T, s=1,...,S, we can write (2) as
with \(\mathbf {Z}^{\prime } = \left [\mathbf {p} \hspace {1mm} \frac {1}{2}\mathbf {Z}\right ]\). P_{ts} is a T×S matrix of success probabilities, Z is a T×C binary matrix and p is a column vector with all its elements equal to p.
Each column of matrix Z represents the SNV profile of a tumor subclone and each column of matrix W represents the proportions of subclones in a sample. Thus, Z, W, C and p explain the inherent heterogeneity in the tumor samples. We perform a joint inference on all these variables by formulating the system model in a statespace framework and then derive an SMC algorithm to infer all the model parameters.
Statespace formulation
Here, we describe the state transition and the observation models of our statespace formulation of the feature allocation model (solution to (3)). Before going through the details of our formulation, we will briefly describe the prior distribution on a leftordered binary matrix Z that has a finite number of rows and an unknown number of columns [35, 36]. By leftordered, we mean that the columns of the binary matrix are ordered from left to right according to the magnitude of the binary in the columns and the first row is considered the most significant. Mathematically, the distribution is expressed as
where m_{c} represents the number of nonzero entries in the c^{th} column of matrix Z, T represents the finite number of rows in matrix Z, C_{+} represents the number of columns in matrix Z that do not sum to zero. \(H_{T} = {\sum }_{t = 1}^{T} 1/t\) represents the T^{th} harmonic number and C_{h} represents the number of columns in matrix Z that form a sequence of ones and zeros corresponding to the binary representation of the number h when read toptobottom.
Fortunately, the prior distribution described in (4) can be viewed as the outcome of IBP, a sequential generative process for the binary matrix. Given that in an Indian buffet restaurant, we have T customers who come into the restaurant one after the other. Assume that the first customer comes into the restaurant and fills her plate from the first c_{1}=Pois(α) distinct dishes. Then the t^{th} customer chooses a particular dish with probability m_{c}/t, m_{c} being the number of people that have chosen the c^{th} dish before her, and in addition, she adds Pois(α/t) new dishes. Following the dish serving rule, if we record the choices of the T customers on the different dishes as a binary matrix such that an entry is one if the customer chose the dish and zero otherwise, such a matrix is a draw from the distribution in (4) [36]. The IBP process is a sequential process in such a way that the choices of the t^{th} customer are only dependent on the customers that were in the restaurant before her.
In our statespace framework, we designate tumor subclones as the dishes, the genomic loci as the customers and the t^{th} customer as the observation at time t (t^{th} row of the input data). We write z_{t}=[z_{t1},z_{t2},...,z_{tC}], the t^{th} row of Z as the state at time t. Thus, according to the sequential process described by the IBP, our state transition model is written as
where Z_{t−1} represents a binary submatrix of the top t−1 rows in Z. We present the algorithm to draw a sample from (5) in Algorithm 1. In the algorithm, Z_{t} is obtained from Z_{t−1} and if in the process, additional nonzero column(s) is/are created in Z_{t}, i.e., Pois(α/t)>0, then additional row(s) will also be added to matrix W. We reparameterize matrix W to easily accommodate any possible change in its dimension by writing \(w_{cs} = \theta _{cs}/{\sum }_{c^{\prime }= 0}^{C} \theta _{c^{\prime }s}\). In other word, we estimate θ_{cs} instead of w_{cs} and we compute w_{cs} from the estimates of θ_{cs}.
For all the parameters of our statespace model, i.e., matrix W and the relative frequency of variant reads p, we employ random walk model to create artificial dynamics
Thus, (5)(6) describe the system state transition of our statespace framework.
The observation model that describes the measurement at time t in the statespace framework is given by
where y_{t} represents the measurement at time t, the t^{th} row of Y. Note that, given the state value at time t z_{t}), the measurement at this timestep is conditionally independent of all the past measurements Y_{t−1}. Thus, (7) details the observation model for the proposed statespace framework. Finally, (5)  (7) state the proposed statespace framework, comprising of the state transition and the observation models for resolving tumor heterogeneity. In summary, the framework described considers, at time t, the t^{th} row of the input data matrices (Y and V) as the observed measurement at time t. The t^{th} row of the binary genotype matrix Z is treated as the hidden state at time t. The proportions W and the relative frequency p are treated as the parameters of the model.
The SMC algorithm
Here, we present a brief description of the SMC filtering approach [33, 34] to make inference on the states (matrix Z) and the parameters (matrix W and p) of the proposed statespace framework. Assume that we have a dynamic system which has a hidden state variable x_{t}, a measurement variable y_{t}, an initial state model (state model when t=0) and a state transition model for other timesteps (∀t>0). In this paper, x_{t} comprises of two types of variables: continuous variables ϕ_{t}, \(\phi _{t} \in \left \{ p_{0}^{t}, \theta _{cs}^{t}, c = 0,1,...,C, s = 1,...,S\right \}\) and discrete variable z_{t}. Also, (5)  (6) describe the state transition model and (7) describes the observation model. At every timestep, given that we have the sequence of measurements up to the present timestep, i.e., Y_{t}={y_{1},y_{2},...,y_{t}}, we are interested in inferring the unobserved sequence X_{t}={x_{1},x_{2},...,x_{t}}.
If we can obtain samples (particles) from the posterior distribution p(X_{t}Y_{t}), then p(X_{t}Y_{t}) can be approximated by the drawn particles. But in most cases, obtaining these particles is not viable. One way to get an estimate is by obtaining weighted particles from a different distribution q(X_{t}Y_{t}) that has a support which incorporates the support of p(X_{t}Y_{t}). This distribution is known as importance distribution. Given that we sample N times from q(X_{t}Y_{t}), i.e., \( \left \{ \mathbf {X}_{i} \right \}_{i = 1}^{N}\), the associated weights are computed as
Thus, an approximation \(\hat {p}\left (\mathbf {X}_{t}\mathbf {Y}_{t}\right)\) of the original posterior distribution p(X_{t}Y_{t}) is by
This procedure is termed the importance sampling theory.
Next, we describe the sequential version of the importance sampling theory. The first step is to factorize the posterior distribution of state variables at time t, X_{t}, given all the measurements up to and including at time t Y_{t}, i.e.,
At time t, instead of sampling from the original distribution p(X_{t}Y_{t}) to approximate p(X_{t}Y_{t}), we obtain N weighted particles from the importance distribution q(X_{t}Y_{t}). We write the importance distribution as q(X_{t}Y_{t})=q(x_{t}X_{t−1},Y_{t})q(X_{t−1}Y_{t−1}), and we compute the associated unnormalized weights as
Imagine that at time t−1, we followed the description of the sequential version of importance sampling and we had N particles, \(\left \{ \mathbf {X}_{t1}^{i} \right \}_{i = 1}^{N}\), drawn from q(X_{t−1}Y_{t−1}), and the associated normalized weights given as
From the weighted particles at time t−1, we easily obtain weighted particles at time t, i.e., \(\left \{ \mathbf {X}_{t}^{i} \right \}_{i = 1}^{N} = \left \{ \mathbf {x}_{t}^{i},\mathbf {X}_{t1}^{i} \right \}_{i = 1}^{N}\), where \(\mathbf {x}_{t}^{i} \sim q\left (\mathbf {x}_{t}\mathbf {X}_{t1}^{i},\mathbf {Y}_{t}\right)\). By substituting (12) into (11), the associated unnormalized weights at time t satisfy the recursion
The weights are normalized to sum to one.
The optimal importance distribution that reduces variability due to one step reweighting is \( p\left (\mathbf {x}_{t}^{i}\mathbf {X}_{t1}^{i},\mathbf {Y}_{t}\right)\). This choice reduces the weights equation in (13) to \(\tilde {w}_{t}^{i} \propto w_{t1}^{i} p\left (\mathbf {y}_{t}\mathbf {X}_{t1}^{i},\mathbf {Y}_{t1}\right)\) [44, 45]. However, we only have closed form solutions for \(p\left (\mathbf {x}_{t}^{i}\mathbf {X}_{t1}^{i},\mathbf {Y}_{t}\right)\) and \(p\left (\mathbf {y}_{t}\mathbf {X}_{t1}^{i},\mathbf {Y}_{t1}\right)\) if and only if \(p\left (\mathbf {y}_{t}\mathbf {X}_{t}^{i},\mathbf {Y}_{t1}\right)\) and \(p\left (\mathbf {x}_{t}^{i}\mathbf {X}_{t1}^{i},\mathbf {Y}_{t1}\right)\) are conjugates. Such conjugacy does not exist in our statespace framework. An equally efficient solution is to choose \(p\left (\mathbf {x}_{t}^{i}\mathbf {X}_{t1}^{i}\right)\) in (5)(6) as the importance distribution [46–49]. Because of independence assumption in the model, i.e., \(p\left (\mathbf {x}_{t}^{i}\mathbf {X}_{t1}^{i},\mathbf {Y}_{t1}\right) = p\left (\mathbf {x}_{t}^{i}\mathbf {X}_{t1}^{i}\right)\) and \(p\left (\mathbf {y}_{t}\mathbf {X}_{t}^{i},\mathbf {Y}_{t1}\right) = p\left (\mathbf {y}_{t}\mathbf {x}_{t}^{i}\right)\), we rewrite (13) as
and then normalize the weights.
As time progresses, there is degeneracy, a condition where the variance of the weights increases [33]. To combat this, we perform resampling at every timestep [46–49]. The resampling procedure [38] is as follows : view each weight \(w_{t}^{i}\) as the probability of obtaining the particle index, draw N particles from the probability distribution \(\left \lbrace w_{t}^{i} \right \rbrace \), replace the old particles with the newly drawn particles and set the new weights to a constant value 1/N.
The proposed sequential algorithm, SeqClone, for estimating the states variables and the parameters of our statespace framework is highlighted in Algorithm 2. To initialize the algorithm, we assume the following prior distributions of the model parameters
In this way, we have \(w_{cs} = \theta _{cs}/{\sum }_{c^{\prime } = 0}^{C} \theta _{c^{\prime }s}\) and as a result, \({\sum }_{c^{\prime }= 0}^{C} w_{c^{\prime }s} = 1\). At every time step of the algorithm, we adaptively perturb the particles of the parameters in ϕ_{t} by choosing σ=2% of the value of the particle. We report the posterior estimates of all the state variables and model parameters using the method described in [50]. We detail this in Additional file 1.
Abbreviations
 CLL:

Chronic lymphocytic leukemia
 EM:

Expectation maximization
 IBP:

Indian buffet process
 MCMC:

Markov chain Monte Carlo
 PDF:

Probability density function
 PMF:

Probability mass function
 SMC:

Sequential Monte Carlo
 SNV:

Single nucleotide variant
 VAF:

Variant allele fraction
 WES:

Whole exome sequencing
 WGS:

Whole genome sequencing
References
 1
Gerlinger M, Rowan AJ, Horswell S, Larkin J, Endesfelder D, Gronroos E, Martinez P, Matthews N, Stewart A, Tarpey P, et al. Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. N Engl J Med. 2012; 366(10):883–92.
 2
NikZainal S, Van Loo P, Wedge DC, Alexandrov LB, Greenman CD, Lau KW, Raine K, Jones D, Marshall J, Ramakrishna M, et al. The life history of 21 breast cancers. Cell. 2012; 149(5):994–1007.
 3
Hughes AE, Magrini V, Demeter R, Miller CA, Fulton R, Fulton LL, Eades WC, Elliott K, Heath S, Westervelt P, et al. Clonal architecture of secondary acute myeloid leukemia defined by singlecell sequencing. PLoS Genet. 2014; 10(7):1004462.
 4
Nowell PC. The clonal evolution of tumor cell populations. Science. 1976; 194(4260):23–8.
 5
Marusyk A, Polyak K. Tumor heterogeneity: causes and consequences. Biochim Biophys Acta (BBA)  Rev Cancer. 2010; 1805(1):105–17.
 6
Meacham CE, Morrison SJ. Tumor heterogeneity and cancer cell plasticity. Nature. 2013; 501(7467):328.
 7
Heppner GH. Tumor heterogeneity. Cancer Res. 1984; 44(6):2259–65.
 8
Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. cell. 2011; 144(5):646–74.
 9
Hanahan D, Weinberg RA. The hallmarks of cancer. cell. 2000; 100(1):57–70.
 10
Reuter JA, Spacek DV, Snyder MP. Highthroughput sequencing technologies. Mol Cell. 2015; 58(4):586–97.
 11
Navin N, Krasnitz A, Rodgers L, Cook K, Meth J, Kendall J, Riggs M, Eberling Y, Troge J, Grubor V, et al. Inferring tumor progression from genomic heterogeneity. Genome Res. 2010; 20(1):68–80.
 12
Irish JM, Hovland R, Krutzik PO, Perez OD, Bruserud Ø, Gjertsen BT, Nolan GP. Single cell profiling of potentiated phosphoprotein networks in cancer cells. Cell. 2004; 118(2):217–228.
 13
Xu X, Hou Y, Yin X, Bao L, Tang A, Song L, Li F, Tsang S, Wu K, Wu H, et al. Singlecell exome sequencing reveals singlenucleotide mutation characteristics of a kidney tumor. Cell. 2012; 148(5):886–95.
 14
Hou Y, Song L, Zhu P, Zhang B, Tao Y, Xu X, Li F, Wu K, Liang J, Shao D, et al. Singlecell exome sequencing and monoclonal evolution of a jak2negative myeloproliferative neoplasm. Cell. 2012; 148(5):873–85.
 15
Navin N, Kendall J, Troge J, Andrews P, Rodgers L, McIndoo J, Cook K, Stepansky A, Levy D, Esposito D, et al. Tumor evolution inferred by single cell sequencing. Nature. 2011; 472(7341):90.
 16
Potter NE, Ermini L, Papaemmanuil E, Cazzaniga G, Vijayaraghavan G, Titley I, Ford A, Campbell P, Kearney L, Greaves M. Singlecell mutational profiling and clonal phylogeny in cancer. Genome Res. 2013; 23(12):2115–25.
 17
Marusyk A, Almendro V, Polyak K. Intratumour heterogeneity: a looking glass for cancer?. Nat Rev Cancer. 2012; 12(5):323.
 18
Beerenwinkel N, Schwarz RF, Gerstung M, Markowetz F. Cancer evolution: mathematical models and computational inference. Syst Biol. 2014; 64(1):1–25.
 19
Schwarz RF, Ng CK, Cooke SL, Newman S, Temple J, Piskorz AM, Gale D, Sayal K, Murtaza M, Baldwin PJ, et al. Spatial and temporal heterogeneity in highgrade serous ovarian cancer: a phylogenetic analysis. PLoS Med. 2015; 12(2):1001789.
 20
Garraway LA, Lander ES. Lessons from the cancer genome. Cell. 2013; 153(1):17–37.
 21
Roth A, Khattra J, Yap D, Wan A, Laks E, Biele J, Ha G, Aparicio S, BouchardCôté A, Shah SP. Pyclone: statistical inference of clonal population structure in cancer. Nat Methods. 2014; 11(4):396–8.
 22
Jiao W, Vembu S, Deshwar AG, Stein L, Morris Q. Inferring clonal evolution of tumors from single nucleotide somatic mutations. BMC Bioinformatics. 2014; 15(1):35.
 23
Miller CA, White BS, Dees ND, Griffith M, Welch JS, Griffith OL, Vij R, Tomasson MH, Graubert TA, Walter MJ, et al. Sciclone: inferring clonal architecture and tracking the spatial and temporal patterns of tumor evolution. PLoS Comput Biol. 2014; 10(8):1003665.
 24
Phylowgs: reconstructing subclonal composition and evolution from wholegenome sequencing of tumors. Genome Biol. 2015; 16(1):35.
 25
Yuan K, Sakoparnig T, Markowetz F, Beerenwinkel Nq. Bitphylogeny: a probabilistic framework for reconstructing intratumor phylogenies. Genome Biol. 2015; 16(1):36.
 26
Marass F, Mouliere F, Yuan K, Rosenfeld N, Markowetz F, et al. A phylogenetic latent feature model for clonal deconvolution. Ann Appl Stat. 2016; 10(4):2377–404.
 27
Zare H, Wang J, Hu A, Weber K, Smith J, Nickerson D, Song C, Witten D, Blau CA, Noble WS. Inferring clonal composition from multiple sections of a breast cancer. PLoS Comput Biol. 2014; 10(7):1003703.
 28
Sengupta S, Wang J, Lee J, Müller P, Gulukota K, Banerjee A, Ji Y. Bayclone: Bayesian nonparametric inference of tumor subclones using ngs data. Pac Symp Biocomput. 2015; 20:467–78.
 29
Fischer A, VázquezGarcía I, Illingworth CJ, Mustonen V. Highdefinition reconstruction of clonal composition in cancer. Cell Rep. 2014; 7(5):1740–52.
 30
Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the em algorithm. J R Stat Soc Ser B (Methodol). 1977; 39:1–38.
 31
Green PJ. Reversible jump markov chain monte carlo computation and bayesian model determination. Biometrika. 1995; 82(4):711–32.
 32
Hastings WK. Monte carlo sampling methods using markov chains and their applications. Biometrika. 1970; 57(1):97–109.
 33
Doucet A, De Freitas N, Gordon N. Sequential monte carlo methods in practice springer. New York:2001.
 34
Doucet A, Godsill S, Andrieu C. On sequential monte carlo sampling methods for bayesian filtering. Stat Comput. 2000; 10(3):197–208.
 35
Griffiths TL, Ghahramani Z. The indian buffet process: An introduction and review. J Mach Learn Res. 2011; 12(Apr):1185–224.
 36
Ghahramani Z, Griffiths TL. Infinite latent feature models and the indian buffet process. In: Advances in Neural Information Processing Systems. Cambridge: MIT Press: 2006. p. 475–82.
 37
Ayinde BO, Zurada JM. Deep learning of constrained autoencoders for enhanced understanding of data. arXiv preprint arXiv:1802.00003. 2018; 29:3969–79.
 38
Arulampalam MS, Maskell S, Gordon N, Clapp T. A tutorial on particle filters for online nonlinear/nongaussian bayesian tracking. IEEE Trans Sig Process. 2002; 50(2):174–88.
 39
Ogundijo OE, Elmas A, Wang X. Reverse engineering gene regulatory networks from measurement with missing values. EURASIP J Bioinforma Syst Biol. 2017; 2017(1):2.
 40
Ogundijo OE, Wang X. A sequential monte carlo approach to gene expression deconvolution. PloS ONE. 2017; 12(10):0186167.
 41
Ogundijo OE, Wang X. Characterization of tumor heterogeneity by latent haplotypes: a sequential monte carlo approach. PeerJ. 2018; 6:4838.
 42
Schuh A, Becq J, Humphray S, Alexa A, Burns A, Clifford R, Feller SM, Grocock R, Henderson S, Khrebtukova I, et al. Monitoring chronic lymphocytic leukemia progression by whole genome sequencing reveals heterogeneous clonal evolution patterns. Blood. 2012; 120(20):4191–6.
 43
ElKebir M, Oesper L, AchesonField H, Raphael BJ. Reconstruction of clonal trees and tumor composition from multisample sequencing data. Bioinformatics. 2015; 31(12):62–70.
 44
Jajamovich GH, Wang X, Arkin AP, Samoilov MS. Bayesian multipleinstance motif discovery with bambi: inference of recombinase and transcription factor binding sites. Nucleic Acids Res. 2011; 39(21):e146.
 45
Ristic B, Arulampalam S, Gordon N. Beyond the kalman filter. IEEE Aerosp Electron Syst Mag. 2004; 19(7):37–8.
 46
Wood F, Griffiths TL. Particle filtering for nonparametric bayesian matrix factorization. In: Advances in Neural Information Processing Systems.2007. p. 1513–20.
 47
Särkkä S, Vol. 3. Bayesian Filtering and Smoothing. Cambridge: Cambridge University Press; 2013.
 48
Li P, Goodall R, Kadirkamanathan V. Estimation of parameters in a linear state space model using a raoblackwellised particle filter. IEE Proc Control Theory Appl. 2004; 151(6):727–38.
 49
Li P, Goodall R, Kadirkamanathan V. Parameter estimation of railway vehicle dynamic model using raoblackwellised particle filter. In: European Control Conference (ECC), 2003. IEEE: 2003. p. 2384–9.
 50
Lee J, Müller P, Sengupta S, Gulukota K, Ji Y. Bayesian feature allocation models for tumor heterogeneity. In: Statistical Analysis for HighDimensional Data.Cham: 2016. p. 211–32.
Acknowledgements
We thank the Petroleum Technology Development Fund, the body responsible for the doctoral sponsorship of Oyetunji Ogundijo, one of the authors.
Funding
No specific funding was received for this study.
Availability of data and materials
The datasets analyzed during the current study are available to download at: https://github.com/moyanre/seqclone. Similarly, MATLAB code is also available to download at: https://github.com/moyanre/seqclone.
Author information
Affiliations
Contributions
OE and XW conceived the project idea and the design of the methods. OE performed the computer experiments and contributed in the writing of the draft. XW reviewed the draft for submission. Both authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1
Supplementary Material for “SeqClone: Sequential Monte Carlo Based Inference of Tumor Subclones”. (PDF 203 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Ogundijo, O.E., Wang, X. SeqClone: sequential Monte Carlo based inference of tumor subclones. BMC Bioinformatics 20, 6 (2019). https://doi.org/10.1186/s128590182562y
Received:
Accepted:
Published:
Keywords
 Tumor heterogeneity
 Bayesian model
 Sequential Monte Carlo
 Indian buffet process