Skip to main content

SeqClone: sequential Monte Carlo based inference of tumor subclones



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.


In this paper, we propose a state-space 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.


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.


Tumor samples that are obtained temporally or spatially from a cancer patient are heterogeneous in nature [1, 2]. These samples contain genetically diverse sub-population 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 [57]. 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 high-throughput 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 [1316]. 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 [2125].

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) [2629]. 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 expectation-maximization (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 state-space 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 state-space formulation considers the sequential construction of the binary genotype matrix by making use of Indian buffet process (IBP) [3537]. 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 state-space 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 [3841]. 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.


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., vts from Pois(r). We generated each entry of Y, i.e., yts as follows: sampled each column of the proportion matrix W independently from Dir([a0,a1,...,a4]) (a0=0.2 and ac, c{1,...,4} randomly chosen from the set {2,4,5,6,7,8}), defined p=0.02, computed pts following (2) in the “Method” section, and sampled yts from binomial(vts,pts).

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 (eZ), proportion error (eW) and success probabilities error (\(e_{p_{ts}}\)). Mathematically, these errors are defined as

$$ \begin{aligned} e_{Z} \,=\, \frac{1}{TC} \sum\limits_{t = 1}^{T} \sum\limits_{c = 1}^{C} \left| \hat{z}_{tc} \,-\, z_{tc} \right|, \hspace{1mm} e_{W} = \frac{1}{CS} \sum\limits_{c = 0}^{C} \sum\limits_{s = 1}^{S} \left| \hat{w}_{cs} - w_{cs} \right|, \end{aligned} $$


$$ \begin{aligned} e_{p_{ts}} \,=\, \frac{1}{TS} \sum\limits_{t= 1}^{T} \sum\limits_{s = 1}^{S} \left| \hat{p}_{ts} \,-\, p_{ts} \right|, \hspace{1mm} \text{where} \hspace{1mm} \hat{p}_{ts} \,=\, \hat{p} \hat{w}_{0s} \,+\, \sum\limits_{c = 1}^{C} \hat{z}_{tc} \hat{w}_{cs}. \end{aligned} $$

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 64-bit 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 eZ and proportion error eW 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.

Fig. 1

Plot of genotype error eZ versus sample size S for T=20 loci, average sequencing depth r=1000 and C{3,4,5} subclones

Fig. 2

Plot of proportion error eW versus sample size S for T=20 loci, average sequencing depth r=1000 and C{3,4,5} subclones

Fig. 3

Plot of the error of success probability \(e_{p_{ts}}\) versus sample size S for T=20 loci, average sequencing depth r=1000 and C{3,4,5} subclones

Fig. 4

Plot of genotype error eZ versus sample size S for T=20 loci, C=3 subclones and average sequencing depth r{50,200,1000}

Fig. 5

Plot of proportion error eW versus sample size S for T=20 loci, C=3 subclones and average sequencing depth r{50,200,1000}

Fig. 6

Plot of the error of success probability \(e_{p_{ts}}\) versus sample size S for T=20 loci, C=3 subclones and average sequencing depth r{50,200,1000}

Fig. 7

Plot of consumed memory versus number of genomic loci T{20,40,60,80,100}

Table 1 Genotype error (eZ) and proportion error (eW) computed for SeqClone, Clomial, BayClone and Cloe for T=20, C=3, S=5 and r{50,200,1000}
Table 2 Genotype error (eZ) and proportion error (eW) computed for SeqClone, Clomial, BayClone and Cloe for T=100, C=3, S=5 and r{50,200,1000}
Table 3 Average and standard deviation of \(e_{p_{ts}}\), eZ and eW for T{100,5000}, S=5, C=3, and r=1000
Table 4 Runtimes, eZ and eW for T=20, S=5, C=3, and r=1000

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 MCMC-based 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 B-cell CLL namely: CLL077, CLL006, and CLL003 [42]. Complete datasets and the data pre-processing steps are in [42]. In Additional file 1, we include the analysis results with Clomial, BayClone and Cloe.


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.

Fig. 8

Phylogenetic tree from CLL077

Table 5 CLL077: estimate of genotype matrix/mutational profile
Table 6 CLL077: estimate of the proportions of subclones in the samples


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.

Fig. 9

Phylogenetic tree from CLL006

Table 7 CLL006: estimate of genotype matrix/mutational profile
Table 8 CLL006: estimate of the proportions of subclones in the samples


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.

Fig. 10

Phylogenetic tree from CLL003

Table 9 CLL003: estimate of genotype matrix/mutational profile
Table 10 CLL003: estimate of the proportions of subclones in the samples

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.

Table 11 Runtimes and memory consumption for simulated and real biological dataset


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 next-generation 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 state-space 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 state-state 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 state-space 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.


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.


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 yts. Likewise, we denote the total number of reads at locus t in sample s as vts. 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

$$ \begin{aligned} y_{ts} \stackrel{{ind.}}{\sim} \text{binomial}\left(v_{ts},p_{ts}\right), \hspace{2mm} t = 1,...,T, \hspace{1mm} s = 1,...,S, \end{aligned} $$

pts, t=1,...,T, s=1,...,S are the success probabilities defined as [28]

$$ \begin{aligned} p_{ts} = w_{0s}p + \frac{1}{2}\sum\limits_{c = 1}^{C} z_{tc}w_{cs}, \end{aligned} $$

where ztc, 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 ztc=1, it implies that locus t in subclone c has reads that bear variant sequence. Likewise, if ztc=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 pts 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 w0sp 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

$$ \begin{aligned} \mathbf{P}_{ts} = \mathbf{Z}^{\prime} \cdot \mathbf{W}, \end{aligned} $$

with \(\mathbf {Z}^{\prime } = \left [\mathbf {p} \hspace {1mm} \frac {1}{2}\mathbf {Z}\right ]\). Pts 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 state-space framework and then derive an SMC algorithm to infer all the model parameters.

State-space formulation

Here, we describe the state transition and the observation models of our state-space 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 left-ordered binary matrix Z that has a finite number of rows and an unknown number of columns [35, 36]. By left-ordered, 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

$$ {{}\begin{aligned} P(\mathbf{Z}) = \frac{\alpha^{C_{+}}}{\prod_{h=1}^{2^{T}-1}C_{h}!} \exp \left\{ -\alpha H_{T} \right\} \prod\limits_{c=1}^{C^{+}} \frac{ \left(T - m_{c}\right)! \left(m_{c} - 1\right)!}{T!}, \end{aligned}} $$

where mc represents the number of non-zero entries in the cth 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 Tth harmonic number and Ch 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 top-to-bottom.

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 c1=Pois(α) distinct dishes. Then the tth customer chooses a particular dish with probability mc/t, mc being the number of people that have chosen the cth 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 tth customer are only dependent on the customers that were in the restaurant before her.

In our state-space framework, we designate tumor subclones as the dishes, the genomic loci as the customers and the tth customer as the observation at time t (tth row of the input data). We write zt=[zt1,zt2,...,ztC], the tth 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

$$ \begin{aligned} P\left(\mathbf{z}_{t}|\mathbf{Z}_{t-1},\alpha\right), \end{aligned} $$

where Zt−1 represents a binary sub-matrix of the top t−1 rows in Z. We present the algorithm to draw a sample from (5) in Algorithm 1. In the algorithm, Zt is obtained from Zt−1 and if in the process, additional non-zero column(s) is/are created in Zt, i.e., Pois(α/t)>0, then additional row(s) will also be added to matrix W. We re-parameterize 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 wcs and we compute wcs from the estimates of θcs.

For all the parameters of our state-space model, i.e., matrix W and the relative frequency of variant reads p, we employ random walk model to create artificial dynamics

$$ \begin{aligned} \phi_{t} \sim p\left(\phi_{t}|\phi_{t-1}\right) = \mathcal{N}\left(\phi_{t-1},\sigma^{2}\right), \\ \phi_{t} \in \left\{ p, \theta_{cs}, c = 0,1,...,C, s = 1,...,S\right\}. \end{aligned} $$

Thus, (5)-(6) describe the system state transition of our state-space framework.

The observation model that describes the measurement at time t in the state-space framework is given by

$$ \begin{aligned} \mathbf{y}_{t} \sim P\left(\mathbf{y}_{t}|\mathbf{Z}_{1:t},\mathbf{W},p\right) &= P(\mathbf{y}_{t}|\mathbf{z}_{t},\mathbf{W},p) \\ & = \prod\limits_{s=1}^{S}\text{binomial}\left(y_{ts}|v_{ts},p_{ts}\right), \end{aligned} $$

where yt represents the measurement at time t, the tth row of Y. Note that, given the state value at time t zt), the measurement at this time-step is conditionally independent of all the past measurements Yt−1. Thus, (7) details the observation model for the proposed state-space framework. Finally, (5) - (7) state the proposed state-space framework, comprising of the state transition and the observation models for resolving tumor heterogeneity. In summary, the framework described considers, at time t, the tth row of the input data matrices (Y and V) as the observed measurement at time t. The tth 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 state-space framework. Assume that we have a dynamic system which has a hidden state variable xt, a measurement variable yt, an initial state model (state model when t=0) and a state transition model for other time-steps (t>0). In this paper, xt 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 zt. Also, (5) - (6) describe the state transition model and (7) describes the observation model. At every time-step, given that we have the sequence of measurements up to the present time-step, i.e., Yt={y1,y2,...,yt}, we are interested in inferring the unobserved sequence Xt={x1,x2,...,xt}.

If we can obtain samples (particles) from the posterior distribution p(Xt|Yt), then p(Xt|Yt) 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(Xt|Yt) that has a support which incorporates the support of p(Xt|Yt). This distribution is known as importance distribution. Given that we sample N times from q(Xt|Yt), i.e., \( \left \{ \mathbf {X}_{i} \right \}_{i = 1}^{N}\), the associated weights are computed as

$$ {{}\begin{aligned} \tilde{w}_{t}^{i} = \frac{p\left(\mathbf{X}_{t}|\mathbf{Y}_{t}\right)}{q\left(\mathbf{X}_{t}|\mathbf{Y}_{t}\right)} \hspace{2mm} \text{and} \hspace{2mm} w_{t}^{i} = \frac{\tilde{w}_{t}^{i}}{{\sum}_{m = 1}^{N}\tilde{w}_{t}^{m}}, \hspace{1mm} i = 1,...,N. \end{aligned}} $$

Thus, an approximation \(\hat {p}\left (\mathbf {X}_{t}|\mathbf {Y}_{t}\right)\) of the original posterior distribution p(Xt|Yt) is by

$$ \hat{p}\left(\mathbf{X}_{t}|\mathbf{Y}_{t}\right) \,=\, \!\sum\limits_{i = 1}^{N} w_{t}^{i} \delta\left(\mathbf{X}_{t} \,-\, \mathbf{X}_{t}^{i}\right), \hspace{1mm} \!\text{where} \hspace{1mm} \delta(\mathbf{u}) \,=\, \left\{\begin{array}{ll} \!\!1, & \text{if}\ \mathbf{u} \,=\, \underline{\mathbf{0}} \\ \!\!0, & \text{otherwise}. \end{array}\right. $$

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, Xt, given all the measurements up to and including at time t Yt, i.e.,

$$ \begin{aligned} p\left(\mathbf{X}_{t}|\mathbf{Y}_{t}\right) &\propto p\left(\mathbf{y}_{t}|\mathbf{X}_{t},\mathbf{Y}_{t-1}\right) p\left(\mathbf{X}_{t}|\mathbf{Y}_{t-1}\right) \\ & \,=\, p\left(\mathbf{y}_{t}|\mathbf{X}_{t},\mathbf{Y}_{t-1}\right) p\left(\mathbf{x}_{t}|\mathbf{X}_{t\,-\,1},\mathbf{Y}_{t\,-\,1}\right) p\left(\mathbf{X}_{t-1}|\mathbf{Y}_{t-1}\right). \end{aligned} $$

At time t, instead of sampling from the original distribution p(Xt|Yt) to approximate p(Xt|Yt), we obtain N weighted particles from the importance distribution q(Xt|Yt). We write the importance distribution as q(Xt|Yt)=q(xt|Xt−1,Yt)q(Xt−1|Yt−1), and we compute the associated unnormalized weights as

$$ \begin{aligned} \tilde{w}_{t}^{i} = \frac{p\left(\mathbf{y}_{t}|\mathbf{X}_{t}^{i},\mathbf{Y}_{t-1}\right) p\left(\mathbf{x}_{t}^{i}|\mathbf{X}_{t-1}^{i},\mathbf{Y}_{t-1}\right)}{q\left(\mathbf{x}_{t}^{i}|\mathbf{X}_{t}^{i},\mathbf{Y}_{t}\right)} \frac{p\left(\mathbf{X}_{t-1}^{i}|\mathbf{Y}_{t-1}\right)}{q\left(\mathbf{X}_{t-1}^{i}|\mathbf{Y}_{t-1}\right)}. \end{aligned} $$

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}_{t-1}^{i} \right \}_{i = 1}^{N}\), drawn from q(Xt−1|Yt−1), and the associated normalized weights given as

$$ \begin{aligned} w_{t-1}^{i} \propto \frac{p\left(\mathbf{X}_{t-1}^{i}|\mathbf{Y}_{t-1}\right)}{q\left(\mathbf{X}_{t-1}^{i}|\mathbf{Y}_{t-1}\right)}, \hspace{2mm} i = 1,...,N. \end{aligned} $$

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}_{t-1}^{i} \right \}_{i = 1}^{N}\), where \(\mathbf {x}_{t}^{i} \sim q\left (\mathbf {x}_{t}|\mathbf {X}_{t-1}^{i},\mathbf {Y}_{t}\right)\). By substituting (12) into (11), the associated unnormalized weights at time t satisfy the recursion

$$ \begin{aligned} \tilde{w}_{t}^{i} \propto w_{t-1}^{i} \frac{p\left(\mathbf{y}_{t}|\mathbf{X}_{t}^{i},\mathbf{Y}_{t-1}\right) p\left(\mathbf{x}_{t}^{i}|\mathbf{X}_{t-1}^{i},\mathbf{Y}_{t-1}\right)}{q\left(\mathbf{x}_{t}^{i}|\mathbf{X}_{t}^{i},\mathbf{Y}_{t}\right)}, \hspace{2mm} i = 1,...,N. \end{aligned} $$

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}_{t-1}^{i},\mathbf {Y}_{t}\right)\). This choice reduces the weights equation in (13) to \(\tilde {w}_{t}^{i} \propto w_{t-1}^{i} p\left (\mathbf {y}_{t}|\mathbf {X}_{t-1}^{i},\mathbf {Y}_{t-1}\right)\) [44, 45]. However, we only have closed form solutions for \(p\left (\mathbf {x}_{t}^{i}|\mathbf {X}_{t-1}^{i},\mathbf {Y}_{t}\right)\) and \(p\left (\mathbf {y}_{t}|\mathbf {X}_{t-1}^{i},\mathbf {Y}_{t-1}\right)\) if and only if \(p\left (\mathbf {y}_{t}|\mathbf {X}_{t}^{i},\mathbf {Y}_{t-1}\right)\) and \(p\left (\mathbf {x}_{t}^{i}|\mathbf {X}_{t-1}^{i},\mathbf {Y}_{t-1}\right)\) are conjugates. Such conjugacy does not exist in our state-space framework. An equally efficient solution is to choose \(p\left (\mathbf {x}_{t}^{i}|\mathbf {X}_{t-1}^{i}\right)\) in (5)-(6) as the importance distribution [4649]. Because of independence assumption in the model, i.e., \(p\left (\mathbf {x}_{t}^{i}|\mathbf {X}_{t-1}^{i},\mathbf {Y}_{t-1}\right) = p\left (\mathbf {x}_{t}^{i}|\mathbf {X}_{t-1}^{i}\right)\) and \(p\left (\mathbf {y}_{t}|\mathbf {X}_{t}^{i},\mathbf {Y}_{t-1}\right) = p\left (\mathbf {y}_{t}|\mathbf {x}_{t}^{i}\right)\), we rewrite (13) as

$$ \begin{aligned} \tilde{w}_{t}^{i} & \propto w_{t-1}^{i} p\left(\mathbf{y}_{t}|\mathbf{x}_{t}^{i}\right) \\ & = w_{t-1}^{i} p\left(\mathbf{y}_{t}|\mathbf{z}_{t}^{i}, \mathbf{W}_{t}^{i}\right), \end{aligned} $$

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 time-step [4649]. 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 state-space framework is highlighted in Algorithm 2. To initialize the algorithm, we assume the following prior distributions of the model parameters

$$ \begin{aligned} \theta_{0s} & \stackrel{{i.i.d}}{\sim} \text{gamma}\left(a_{0},1\right), \hspace{1mm} s = 1,...,S,\\ \theta_{cs} & \stackrel{i.i.d}{\sim} \text{gamma}\left(a_{1},1\right), \hspace{1mm} s = 1,...,S,c = 1,...,C, \hspace{1mm} \text{and}\\ p & \sim \text{beta}\left(a_{00},b_{00}\right). \end{aligned} $$

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.



Chronic lymphocytic leukemia


Expectation maximization


Indian buffet process


Markov chain Monte Carlo


Probability density function


Probability mass function


Sequential Monte Carlo


Single nucleotide variant


Variant allele fraction


Whole exome sequencing


Whole genome sequencing


  1. 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.

    CAS  Article  Google Scholar 

  2. 2

    Nik-Zainal 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.

    CAS  Article  Google Scholar 

  3. 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 single-cell sequencing. PLoS Genet. 2014; 10(7):1004462.

    Article  Google Scholar 

  4. 4

    Nowell PC. The clonal evolution of tumor cell populations. Science. 1976; 194(4260):23–8.

    CAS  Article  Google Scholar 

  5. 5

    Marusyk A, Polyak K. Tumor heterogeneity: causes and consequences. Biochim Biophys Acta (BBA) - Rev Cancer. 2010; 1805(1):105–17.

    CAS  Article  Google Scholar 

  6. 6

    Meacham CE, Morrison SJ. Tumor heterogeneity and cancer cell plasticity. Nature. 2013; 501(7467):328.

    CAS  Article  Google Scholar 

  7. 7

    Heppner GH. Tumor heterogeneity. Cancer Res. 1984; 44(6):2259–65.

    CAS  PubMed  Google Scholar 

  8. 8

    Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. cell. 2011; 144(5):646–74.

    CAS  Article  Google Scholar 

  9. 9

    Hanahan D, Weinberg RA. The hallmarks of cancer. cell. 2000; 100(1):57–70.

    CAS  Article  Google Scholar 

  10. 10

    Reuter JA, Spacek DV, Snyder MP. High-throughput sequencing technologies. Mol Cell. 2015; 58(4):586–97.

    CAS  Article  Google Scholar 

  11. 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.

    CAS  Article  Google Scholar 

  12. 12

    Irish JM, Hovland R, Krutzik PO, Perez OD, Bruserud Ø, Gjertsen BT, Nolan GP. Single cell profiling of potentiated phospho-protein networks in cancer cells. Cell. 2004; 118(2):217–228.

    CAS  Article  Google Scholar 

  13. 13

    Xu X, Hou Y, Yin X, Bao L, Tang A, Song L, Li F, Tsang S, Wu K, Wu H, et al. Single-cell exome sequencing reveals single-nucleotide mutation characteristics of a kidney tumor. Cell. 2012; 148(5):886–95.

    CAS  Article  Google Scholar 

  14. 14

    Hou Y, Song L, Zhu P, Zhang B, Tao Y, Xu X, Li F, Wu K, Liang J, Shao D, et al. Single-cell exome sequencing and monoclonal evolution of a jak2-negative myeloproliferative neoplasm. Cell. 2012; 148(5):873–85.

    CAS  Article  Google Scholar 

  15. 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.

    CAS  Article  Google Scholar 

  16. 16

    Potter NE, Ermini L, Papaemmanuil E, Cazzaniga G, Vijayaraghavan G, Titley I, Ford A, Campbell P, Kearney L, Greaves M. Single-cell mutational profiling and clonal phylogeny in cancer. Genome Res. 2013; 23(12):2115–25.

    CAS  Article  Google Scholar 

  17. 17

    Marusyk A, Almendro V, Polyak K. Intra-tumour heterogeneity: a looking glass for cancer?. Nat Rev Cancer. 2012; 12(5):323.

    CAS  Article  Google Scholar 

  18. 18

    Beerenwinkel N, Schwarz RF, Gerstung M, Markowetz F. Cancer evolution: mathematical models and computational inference. Syst Biol. 2014; 64(1):1–25.

    Article  Google Scholar 

  19. 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 high-grade serous ovarian cancer: a phylogenetic analysis. PLoS Med. 2015; 12(2):1001789.

    Article  Google Scholar 

  20. 20

    Garraway LA, Lander ES. Lessons from the cancer genome. Cell. 2013; 153(1):17–37.

    CAS  Article  Google Scholar 

  21. 21

    Roth A, Khattra J, Yap D, Wan A, Laks E, Biele J, Ha G, Aparicio S, Bouchard-Côté A, Shah SP. Pyclone: statistical inference of clonal population structure in cancer. Nat Methods. 2014; 11(4):396–8.

    CAS  Article  Google Scholar 

  22. 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.

    Article  Google Scholar 

  23. 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.

    Article  Google Scholar 

  24. 24

    Phylowgs: reconstructing subclonal composition and evolution from whole-genome sequencing of tumors. Genome Biol. 2015; 16(1):35.

  25. 25

    Yuan K, Sakoparnig T, Markowetz F, Beerenwinkel Nq. Bitphylogeny: a probabilistic framework for reconstructing intra-tumor phylogenies. Genome Biol. 2015; 16(1):36.

    Article  Google Scholar 

  26. 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.

    Article  Google Scholar 

  27. 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.

    Article  Google Scholar 

  28. 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.

    Google Scholar 

  29. 29

    Fischer A, Vázquez-García I, Illingworth CJ, Mustonen V. High-definition reconstruction of clonal composition in cancer. Cell Rep. 2014; 7(5):1740–52.

    CAS  Article  Google Scholar 

  30. 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.

    Google Scholar 

  31. 31

    Green PJ. Reversible jump markov chain monte carlo computation and bayesian model determination. Biometrika. 1995; 82(4):711–32.

    Article  Google Scholar 

  32. 32

    Hastings WK. Monte carlo sampling methods using markov chains and their applications. Biometrika. 1970; 57(1):97–109.

    Article  Google Scholar 

  33. 33

    Doucet A, De Freitas N, Gordon N. Sequential monte carlo methods in practice springer. New York:2001.

  34. 34

    Doucet A, Godsill S, Andrieu C. On sequential monte carlo sampling methods for bayesian filtering. Stat Comput. 2000; 10(3):197–208.

    Article  Google Scholar 

  35. 35

    Griffiths TL, Ghahramani Z. The indian buffet process: An introduction and review. J Mach Learn Res. 2011; 12(Apr):1185–224.

    Google Scholar 

  36. 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.

    Google Scholar 

  37. 37

    Ayinde BO, Zurada JM. Deep learning of constrained autoencoders for enhanced understanding of data. arXiv preprint arXiv:1802.00003. 2018; 29:3969–79.

    Google Scholar 

  38. 38

    Arulampalam MS, Maskell S, Gordon N, Clapp T. A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking. IEEE Trans Sig Process. 2002; 50(2):174–88.

    Article  Google Scholar 

  39. 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.

    Article  Google Scholar 

  40. 40

    Ogundijo OE, Wang X. A sequential monte carlo approach to gene expression deconvolution. PloS ONE. 2017; 12(10):0186167.

    Article  Google Scholar 

  41. 41

    Ogundijo OE, Wang X. Characterization of tumor heterogeneity by latent haplotypes: a sequential monte carlo approach. PeerJ. 2018; 6:4838.

    Article  Google Scholar 

  42. 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.

    CAS  Article  Google Scholar 

  43. 43

    El-Kebir M, Oesper L, Acheson-Field H, Raphael BJ. Reconstruction of clonal trees and tumor composition from multi-sample sequencing data. Bioinformatics. 2015; 31(12):62–70.

    Article  Google Scholar 

  44. 44

    Jajamovich GH, Wang X, Arkin AP, Samoilov MS. Bayesian multiple-instance motif discovery with bambi: inference of recombinase and transcription factor binding sites. Nucleic Acids Res. 2011; 39(21):e146.

    CAS  Article  Google Scholar 

  45. 45

    Ristic B, Arulampalam S, Gordon N. Beyond the kalman filter. IEEE Aerosp Electron Syst Mag. 2004; 19(7):37–8.

    Article  Google Scholar 

  46. 46

    Wood F, Griffiths TL. Particle filtering for nonparametric bayesian matrix factorization. In: Advances in Neural Information Processing Systems.2007. p. 1513–20.

  47. 47

    Särkkä S, Vol. 3. Bayesian Filtering and Smoothing. Cambridge: Cambridge University Press; 2013.

    Google Scholar 

  48. 48

    Li P, Goodall R, Kadirkamanathan V. Estimation of parameters in a linear state space model using a rao-blackwellised particle filter. IEE Proc Control Theory Appl. 2004; 151(6):727–38.

    Article  Google Scholar 

  49. 49

    Li P, Goodall R, Kadirkamanathan V. Parameter estimation of railway vehicle dynamic model using rao-blackwellised particle filter. In: European Control Conference (ECC), 2003. IEEE: 2003. p. 2384–9.

  50. 50

    Lee J, Müller P, Sengupta S, Gulukota K, Ji Y. Bayesian feature allocation models for tumor heterogeneity. In: Statistical Analysis for High-Dimensional Data.Cham: 2016. p. 211–32.

Download references


We thank the Petroleum Technology Development Fund, the body responsible for the doctoral sponsorship of Oyetunji Ogundijo, one of the authors.


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: Similarly, MATLAB code is also available to download at:

Author information




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

Correspondence to Xiaodong Wang.

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(, 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( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Ogundijo, O.E., Wang, X. SeqClone: sequential Monte Carlo based inference of tumor subclones. BMC Bioinformatics 20, 6 (2019).

Download citation


  • Tumor heterogeneity
  • Bayesian model
  • Sequential Monte Carlo
  • Indian buffet process