Characterization of tumor heterogeneity by latent haplotypes: a sequential Monte Carlo approach

State-space formulation

Our state-space formulation of the problem exploits the sequential construction of Z. Specifically, we consider the t^th row of the data matrix Y and V as the new observation at time t of our state-space model, treat the t^th row of the binary matrix Z as the hidden state at time t, and W and p as the model parameters. Before explicitly stating the state transition and the observation models, we succinctly describe the prior distribution on a “left-ordered” binary matrix (i.e., ordering the columns of the binary matrix from left to right by the magnitude of the binary expressed by that column, taking the first row as the most significant bit) with a finite number of rows and an unknown number of columns. The prior distribution as detailed in (Griffiths & Ghahramani, 2011; Ghahramani & Griffiths, 2006) is given by: (3)${\displaystyle \begin{array}{c}P\left(\mathbf{Z}\right)=\frac{{\alpha }^{C_{+}}}{\prod _{h=1}^{2^T-1}{C}_h!}exp\left\{-\alpha H_T\right\}\prod _{c=1}^{C^{+}}\frac{\left(T-m_c\right)!\left(m_c-1\right)!}{T!},\hfill \end{array}}$ where C₊ denotes the number of columns of Z with non-zero entries, m_c denotes the number of 1’s in column c, T denotes the number of rows in Z, $H_T={\sum }_{t=1}^{T}1∕t$ denotes the T^th harmonic number, and C_h denotes the number of columns in Z that when read top-to-bottom form a sequence of 1’s and 0’s corresponding to the binary representation of the number h.

The distribution in (3) can be derived as the outcome of a sequential generative process called the Indian buffet process (Griffiths & Ghahramani, 2011; Doshi-Velez, 2009). Imagine that in an Indian buffet restaurant, we have T customers who arrive at the restaurant sequentially, one after the other. The first customer walks into the restaurant and loads her plate from the first c₁ dishes, where c₁ = Pois(α) (α is similar to the dispersion parameter in the Chinese Restaurant Process (Zhang, 2008)). The t^th customer will choose a particular dish according to the popularity of the dish, i.e., choosing a dish with probability m_c∕t, where m_c denotes the number of people who have previously chosen the c^th dish, and in addition, chooses Pois(α∕t) new dishes as well. Now, if we record the choices of each customer on each row of a matrix, where each column corresponds to a dish on the buffet (1 if the dish is chosen, and 0 if not), then such a binary matrix is a draw from the distribution in (3) (Ghahramani & Griffiths, 2006). The entire process is sequential because the choices made by the t^th customer are dependent only on the choices made by the t − 1 preceeding customers and not on the remaining T − t customers.

In our case, the dishes in the IBP are the haplotypes in the tumor samples, the SNVs are the customers and more importantly, the t^th customer is the observation at timet in our state-space model. Moreover, if we consider z_t = [z_t1, z_t2, ..., z_tC] in (2), which is equivalently the t^th row of Z as the state at time t, then we can write our state transition model, following the sequential process described by the IBP as follows: (4)${\displaystyle \begin{array}{c}P\left({\mathbf{z}}_t|{\mathbf{Z}}_{t-1},\alpha \right),\hfill \end{array}}$ where Z_t−1 denotes the previous t − 1 rows in Z. The algorithm to sample from (4) is presented in Algorithm 1 in the Supplemental Information 1. Note that in the algorithm, Z_t is implicitly constructed from Z_t−1 and if in the process, new non-zero column(s) is/are introduced in Z_t (Pois(α∕t) > 0), then new row(s) will be added to W as well. On the other hand, if the numbers of non-zero columns in Z_t−1 and Z_t are the same, then the number of rows in W does not change between t − 1 and t. To account for any possible change of dimension in W, we re-parameterize matrix W. Specifically, we rewrite $w_{cs}={\theta }_{cs}∕{\sum }_{c^{\prime }=0}^C{\theta }_{c^{\prime }s}$, which implies that we estimate θ_cs and compute w_cs from the estimates of θ_cs. This procedure ensures that each column of W sums to unity at any point in time during the process.

Moreover, since we are interested in the final estimates of the model parameters W and p, we create artificial dynamics for these parameters using the random walk model, i.e., (5)${\displaystyle \begin{array}{c}{\varphi }_t\sim p\left({\varphi }_t|{\varphi }_{t-1}\right)=\mathcal{N}\left({\varphi }_{t-1},{\sigma }^2\right),\hfill \\ {\varphi }_t\in \left\{p,{\theta }_{cs},c=0,1,...,C,s=1,...,S\right\},\hfill \end{array}}$ where σ denotes the standard deviation. Hence, (4)–(5) fully describe the system state transition.

Similarly, the observation at time t is given by: (6)${\displaystyle \begin{array}{cc}{\mathbf{y}}_t\sim P\left({\mathbf{y}}_t|{\mathbf{Z}}_{1:t},\mathbf{W},p\right)\hfill & =P\left({\mathbf{y}}_t|{\mathbf{z}}_t,\mathbf{W},p\right)\hfill \\ \hfill & =\prod _{s=1}^S\text{Binomial}\left(y_{ts}|v_{ts},p_{ts}\right),\hfill \end{array}}$ where y_t denotes the observation at time t (which is conditionally independent of the previous observations Y_t−1 given the state z_t), i.e., the t^th row of Y. (6) fully describes the measurement model for the system. Finally, (4)–(6) completely describe our proposed state-space model for estimating the mutational profile and the proportion of each haplotype, and the total number of haplotypes in the tumor samples.

The SMC algorithm

In this section, we briefly describe the SMC filtering framework that will be employed to estimate the states and the parameters of our state-space model (Doucet, Godsill & Andrieu, 2000; Doucet, De Freitas & Gordon, 2001). Consider the general dynamic system with hidden state variable x_t, in our case, consisting of discrete variables z_t and continuous variables ϕ_t, ${\varphi }_t\in \left\{p_0^t,{\theta }_{cs}^t,c=0,1,...,C,s=1,...,S\right\}$, and measurement variable y_t, where there is an initial state model p(x₀), and ∀t ≥ 1, a state transition model given in (4)–(5) and an observation model given in (6). The sequence X_t = {x₁, x₂, ..., x_t} is not observed and we want to estimate it for each time t, given that the we have the observations Y_t = {y₁, y₂, ..., y_t}.

Our goal is to approximate the posterior distribution of states p(X_t|Y_t) using particles drawn from it. However, getting such particles from p(X_t|Y_t) is usually not feasible. We can still implement an estimate using N particles, ${\left\{{\mathbf{X}}_t^i\right\}}_{i=1}^N$, taken from another distribution, q(X_t|Y_t), whose support includes the support of p(X_t|Y_t) (importance sampling theorem). For the approximation, the weights associated with the particles are calculated as follows: (7)${\displaystyle \begin{array}{c}{\tilde{w}}_t^i=\frac{p\left({\mathbf{X}}_t|{\mathbf{Y}}_t\right)}{q\left({\mathbf{X}}_t|{\mathbf{Y}}_t\right)}\text{and}{w}_t^i=\frac{{\tilde{w}}_t^i}{\sum _{m=1}^N{\tilde{w}}_t^m},i=1,...,N.\hfill \end{array}}$ Thus, the pair ${\left\{{\mathbf{X}}_t^i,w_{1:t}^i\right\}}_{i=1}^N$ is said to be properly weighted with respect to the distribution p(X_t|Y_t), and the approximation $\hat{p}\left({\mathbf{X}}_t|{\mathbf{Y}}_t\right)$ is then given by: (8)${\displaystyle \begin{array}{c}\hat{p}\left({\mathbf{X}}_t|{\mathbf{Y}}_t\right)=\sum _{i=1}^N{w}_t^i\delta \left({\mathbf{X}}_t-{\mathbf{X}}_t^i\right),\text{where}\delta \left(\mathbf{u}\right)=\left\{\begin{array}{cc}1,\hfill & \text{if}\mathbf{u}=\underset{¯}{\mathbf{0}}\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}\right.\hfill \end{array}}$

 
__________________________________________________________________________________ 
Algorithm 1 SMC Algorithm for Characterizing Tumor Heterogeneity 
__________________________________________________________________________________ 
Input: Y, V. 
  1:  Initialize N particles {zi0,pi0,Wi0}Ni=1 
  2:  for  t = 1,...,T do 
  3:     for i = 1,...,N do 
  4:         Sample zit from Zit−1 using Algorithm 1 in the Supplementary Mate- 
     rial. 
  5:         n1 ← number of columns in Zit−1 
  6:         n2 ← length of zit 
  7:         d ← (n2 − n1) 
  8:         if d = 0 then 
  9: 
      Zit ←[Zit−1 
              zit] 
 10:            Sample Wit using (5) 
 11:         else 
 12: 
    Zit ←[Zit−1 0 
             zit] 
 13:            Sample Wit using (5) 
 14:            Sample new rows of Wit from the priors in (14) 
 15:         end if 
 16:         Calculate  ˜ wit using (13) 
 17:     end for 
 18:     Normalize the weights 
 19:     Perform resampling 
 20:  end for 
 21:  Approximations of posterior estimates of all the unknown variables are ob- 
     tained from the final particles and weights, using the procedures highlighted 
     in (Lee et al., 2016) and discussed in the Supplementary Material. 
__________________________________________________________________________________ 
    

Similar to the above importance sampling theory, a sequential algorithm can be obtained as follows. First, we express the full posterior distribution of states X_t given the observations Y_t as follows: (9)${\displaystyle \begin{array}{cc}p\left({\mathbf{X}}_t|{\mathbf{Y}}_t\right)\hfill & \propto p\left({\mathbf{y}}_t|{\mathbf{X}}_t,{\mathbf{Y}}_{t-1}\right)p\left({\mathbf{X}}_t|{\mathbf{Y}}_{t-1}\right)\hfill \\ \hfill & =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).\hfill \end{array}}$ At time t, we desire to obtain N weighted particles from p(X_t|Y_t), which is not feasible. Instead, we define an importance distribution q(X_t|Y_t) = q(x_t|X_t−1, Y_t)q(X_t−1|Y_t−1), where particles can be obtained from, and then calculate the associated unnormalized importance weights as follows: (10)${\displaystyle \begin{array}{c}{\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)}.\hfill \end{array}}$ Assuming that at time t − 1, we have already drawn the particles ${\left\{{\mathbf{X}}_{t-1}^i\right\}}_{i=1}^N$ from the importance distribution q(X_t−1|Y_t−1) and the corresponding normalized weights written as follows: (11)${\displaystyle \begin{array}{c}{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)},i=1,...,N,\hfill \end{array}}$ we can now draw particles ${\left\{{\mathbf{X}}_t^i\right\}}_{i=1}^N$ from the importance distribution q(X_t|Y_t) by drawing the new state particles for the time step t as ${\mathbf{x}}_t^i\sim q\left({\mathbf{x}}_t|{\mathbf{X}}_{t-1}^i,{\mathbf{Y}}_t\right)$, and write ${\left\{{\mathbf{X}}_t^i\right\}}_{i=1}^N={\left\{{\mathbf{x}}_t^i,{\mathbf{X}}_{t-1}^i\right\}}_{i=1}^N$. If we substitute (11) into (10), the weights at time t satisfy the recursion: (12)${\displaystyle \begin{array}{c}{\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)},i=1,...,N,\hfill \end{array}}$ and then the weights are normalized to sum to unity.

So far, we have presented a generic sequential sampling algorithm. We obtain the optimal importance distribution by setting $q\left({\mathbf{x}}_t^i|{\mathbf{X}}_{t-1}^i,{\mathbf{Y}}_t\right)=p\left({\mathbf{x}}_t^i|{\mathbf{X}}_{t-1}^i,{\mathbf{Y}}_t\right)$, and the weights in (12) become ${\tilde{w}}_t^i\propto w_{t-1}^{i}p\left({\mathbf{y}}_t|{\mathbf{X}}_{t-1}^i,{\mathbf{Y}}_{t-1}\right)$ (Ristic, Arulampalam & Gordon, 2004) i.e., if the distributions $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, then closed form solutions can be obtained for $p\left({\mathbf{x}}_t^i|{\mathbf{X}}_{t-1}^i,{\mathbf{Y}}_t\right)$, and hence, $p\left({\mathbf{y}}_t|{\mathbf{X}}_{t-1}^i,{\mathbf{Y}}_{t-1}\right)$. However, if no such conjugacy exists, which is the case for our state-space model, the most popular choice and equally efficient solution (Van Der Merwe, 2004) is to set $q\left({\mathbf{x}}_t^i|{\mathbf{X}}_{t-1}^i,{\mathbf{Y}}_t\right)=p\left({\mathbf{x}}_t^i|{\mathbf{X}}_{t-1}^i\right)$ (in (4)–(5)) (Wood & Griffiths, 2007; Särkkä, 2013). Considering the assumed independence in our 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)$, then (12) becomes: (13)${\displaystyle \begin{array}{c}{\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),\hfill \end{array}}$ and the weights are normalized. Such implementation is commonly referred to as a bootstrap filter in the literature (Särkkä, 2013).

However, the variance of the weights increases over time, a condition referred to as degeneracy in the literature (Doucet, De Freitas & Gordon, 2001). To avoid this, we perform resampling, at every time step, owing to the choice of the importance distribution (Wood & Griffiths, 2007; Särkkä, 2013), discarding the ineffective particles and multiplying the effective ones. The resampling procedure (Särkkä, 2013) is described in the Supplemental Information 1.

Finally, our proposed SMC algorithm for estimating the mutational profiles and the proportions of the haplotypes in the tumor samples i.e., the states and the parameters of our state-space model, is presented in Algorithm 1. The algorithm is initialized by taking particles from the prior distributions of the parameters. We assume the following: (14)${\displaystyle \begin{array}{cc}{\theta }_{0s}\hfill & \stackrel{i.i.d}{\sim }\text{Gamma}\left(a_0,1\right),s=1,...,S,p\sim \text{Beta}\left(a_{00},b_{00}\right)\hfill \\ {\theta }_{cs}\hfill & \stackrel{i.i.d}{\sim }\text{Gamma}\left(a_1,1\right),s=1,...,S,c=1,...,C,\hfill \end{array}}$ such that $w_{cs}={\theta }_{cs}∕{\sum }_{c^{{}^{\prime }}=0}^C{\theta }_{c^{{}^{\prime }}s}$ and consequently, ${\sum }_{c^{{}^{\prime }}=0}^C{w}_{c^{{}^{\prime }}s}=1$ and assume that a₀₀ <  < b₀₀ to impose a small p. We report the posterior estimates of all the unknown variables using the procedure highlighted in Lee et al. (2016), with the details discussed in the Supplemental Information 1.

Simulated data

We produced simulated datasets with average sequencing depth r ∈ {20, 40, 50, 200, 100, 10,000} per locus. We fixed the number of haplotypes C = 4 and number of samples S ∈ {1, 3, 5, 10, 20}. For all combinations of r and S, we generated the variants count matrix Y and the total count matrix V for different number of SNVs T ∈ {20, 60}. Specifically, we generated each entry of V, i.e., v_ts from Pois(r) and to generate each entry of Y, i.e., y_ts, we did the following: (i) generate each column of W from Dir([a₀, a₁, ..., a₄]), where a₀ = 0.2, and a_c, c ∈ {1, ..., 4} is randomly chosen from the set {2, 4, 5, 6, 7, 8}, (ii) generate entries of Z independently from Bern(0.6), (iii) set p = 0.02, (iv) compute p_ts using (2), and finally, (v) generate y_ts as an independent sample from Binomial(v_ts, p_ts).

Next, we run the proposed SMC-based, MCMC-based and the MAP-based algorithms on the simulated Y and V datasets. To quantify the performance of the algorithms, we define the following metrics: haplotype error (e_Z), proportion error (e_W) and the error of the success probabilities (e_pts) as follows: ${\displaystyle \begin{array}{c}{e}_Z=\frac{1}{TC}\sum _{t=1}^T\sum _{c=1}^C|{\hat{z}}_{tc}-z_{tc}|,e_W=\frac{1}{CS}\sum _{c=0}^C\sum _{s=1}^S|{\hat{w}}_{cs}-w_{cs}|,\hfill \end{array}}$ and ${\displaystyle \begin{array}{c}{e}_{p_{ts}}=\frac{1}{TS}\sum _{t=1}^T\sum _{s=1}^S|{\hat{p}}_{ts}-p_{ts}|,\text{where}{\hat{p}}_{ts}=\hat{p}{\hat{w}}_{0s}+\sum _{c=1}^C{\hat{z}}_{tc}{\hat{w}}_{cs}.\hfill \end{array}}$ However, since this is a blind decomposition, one does not know a priori which column of $\hat{\mathbf{Z}}$ corresponds to which column of Z. To resolve this, we calculate e_Z with every permutation of the columns of $\hat{\mathbf{Z}}$ and then select the permutation that results in the smallest e_Z. The selected permutation is then used in computing e_W and e_pts.

The results obtained from the analyses of the simulated datasets are presented in Table 1, Figs. 1–3 and Table 1 in the Supplemental Information 1. Table 1 shows e_pts, e_Z and e_W obtained for the datasets from T = 60 SNVs, C = 4 haplotypes and S = 10 samples for all the average sequencing depth r ∈ { 20, 40, 50, 200, 1,000, 10,000}. Similar results are presented in the Supplemental Information 1 with T = 20 SNVs and S = 5 samples. From the results obtained for the different number of average sequencing depth r and number of SNVs, the proposed SMC algorithm yields more accurate estimates of the model parameters when compared with the other two algorithms. Specifically, the SMC-based algorithm produced lower error values e_pts, e_Z and e_W in all the datasets analyzed. Moreover, we investigated the effect of sample size on the proposed SMC algorithm. As shown in Figs. 1–3, the results show slight improvement as the the number of samples increased. It can be observed that the results are less sensitive to sample size when S > 5. Also noticed is a slight improvement in the results when the average sequencing depth r is increased.

By construction, the proposed SMC algorithm can handle datasets with any number of loci since the VAF of each loci is processed as an observation at every time step. Apart from the ability to process any number of loci, this property allows the proposed algorithm to process VAFs data from newly sequenced candidate SNVs to improve existing estimates without re-analyzing the previous datasets. To validate this, we generated datasets for 1,000 and 2,000 loci respectively and these datasets are analyzed with the proposed SMC algorithm with the results presented in Table 2. In fact, the proposed SMC algorithm benefits from larger number of loci because the larger the number of loci, the better the estimate of the proportions. This result is evident from the proportion errors in Table 2.

Lastly, we record the runtime (t_r) for the two algorithms on a 3.5 Ghz Intel 8 processors running MATLAB when analyzing some of the datasets described in Table 1 (i.e., T = 60, C = 4, S = 10 and r = 20). Observed t_r was 311 seconds and 585 seconds for the proposed SMC and the MCMC-based algorithms, respectively. For the MAP-based algorithm, a single run is 5 seconds but for a set of input data, the algorithm requires different random initializations.

Real tumor samples: CLL datasets

We evaluate the proposed SMC algorithm on the datasets for the B-cell chronic lymphocytic leukemia (CLL), obtained from three patients: CLL003, CLL006, and CLL077 (Schuh et al., 2012). These datasets represent the molecular changes in pre-treatment, post-treatment, and relapse samples in the three selected patients, i.e., the samples were taken temporally (see the Supplemental Information 1 for the summary of data pre-processing). The datasets are analyzed with the proposed SMC algorithm and the two other algorithms.

CLL003

The CLL dataset obtained from patient CLL003 has 20 distinct loci, shown in the first column in Table 3, and the dataset is analyzed with the proposed SMC algorithm. In Table 3, we present the posterior point estimate of the mutational profiles of the haplotypes in each of the 5 samples, where 1 and 0 denote the variant and the reference sequence, respectively. Moreover, in Fig. 4, we present a graphical representation of how the haplotypes are distributed across the samples. For instance, haplotype C2 with approximately 40 percent abundance in the first sample has reduced to approximately 3 percent after the last treatment. In the Supplemental Information 1, we present the table of proportions. The first row on the table and equivalently C0 in Fig. 4 comprises of the proportion of the background haplotype, which accounts for experimental noise in each sample. From Table 3, we found that each sample consists of at least 2 dominant haplotypes. For instance, tumor sample a is dominated by haplotypes C2 and C6, each with a proportion of approximately 0.4. Also, we analyzed the same dataset with the other two algorithms and the results are in the Supplemental Information 1.

Table 3:

CLL003: estimates of the mutational profiles of haplotypes, Z in the samples.

Gene	C1	C2	C3	C4	C5	C6
ADAD1	1	1	1	0	0	0
AMTN	0	1	0	0	0	0
APBB2	0	1	0	0	0	0
ASXL1	1	0	0	1	0	0
ATM	0	1	0	0	1	0
BPIL2	0	1	0	0	0	0
CHRNB2	1	0	0	1	0	0
CHTF8	1	1	1	0	0	0
FAT3	1	0	0	1	0	0
HERC2	1	1	1	0	0	0
IL11RA	1	1	1	0	0	0
MTUS1	0	1	0	0	0	0
MUSK	1	0	0	1	0	0
NPY	1	0	0	1	0	0
NRG3	1	0	0	1	0	0
PLEKHG5	0	1	0	0	0	0
SEMA3E	1	0	0	1	0	0
SF3B1	1	1	1	0	0	0
SHROOM1	1	1	1	0	0	0
SPTAN1	0	1	0	0	0	0

DOI: 10.7717/peerj.4838/table-3

Notes:

The genes where the mutations are found are shown in the first column.

CLL077

The CLL dataset obtained from patient CLL077 has 16 distinct loci, shown in the first column in Table 4, and the dataset is analyzed with the proposed SMC algorithm. In Table 4, we present the posterior point estimate of mutational profiles of the haplotypes in each of the 5 samples. Moreover, in Fig. 5, we present a graphical representation of how the haplotypes are distributed across the samples, with a table of proportions presented in the Supplemental Information 1. From our analysis results, we found that each sample consists of at least two dominant haplotypes. Also, we analyzed the same dataset with the other two algorithms and the results are in the Supplemental Information 1.

Table 4:

CLL077: estimates of the mutational profiles of haplotypes, Z in the samples.

Gene	C1	C2	C3	C4	C5	C6	C7	C8	C9
BCL2L13	1	0	1	1	1	0	1	0	0
COL24A1	0	0	1	0	0	0	0	0	0
DAZAP1	0	0	0	1	1	0	0	1	0
EXOC6B	0	0	0	1	1	0	0	0	1
GHDC	0	0	0	1	1	0	0	0	1
GPR158	1	0	1	1	1	0	0	0	1
HMCN1	0	0	1	0	0	0	0	0	0
KLHDC2	0	0	1	0	0	0	0	0	0
LRRC16A	0	0	0	0	1	0	0	0	0
MAP2K1	0	0	1	0	0	0	0	0	0
NAMPT	1	0	1	1	1	0	1	0	0
NOD1	0	0	1	0	0	0	0	0	0
OCA2	0	0	0	1	1	0	0	0	1
PLA2G16	0	0	0	1	1	0	1	0	0
SAMHD1	0	1	1	1	1	0	1	0	0
SLC12A1	0	1	1	1	1	0	0	0	0

DOI: 10.7717/peerj.4838/table-4

CLL006

Here, we analyze the CLL dataset obtained from patient CLL006. The dataset comprises of 11 loci as shown in Table 5 in the first column, and is analyzed with the proposed SMC algorithm. Table 5 and Fig. 6 show the estimates of mutational profiles and proportions of the haplotypes, respectively. Also, in the Supplemental Information 1, we present the results obtained from analyzing the dataset with the two other algorithms.

Table 5:

CLL006: estimates of the mutational profiles of haplotypes, Z in the samples.

Gene	C1	C2	C3	C4	C5	C6	C7	C8	C9	C10
ARHGAP29	1	1	1	1	0	1	0	0	0	0
EGFR	1	1	1	1	1	0	0	0	0	0
IRF4	1	0	1	0	0	0	0	0	0	0
KIAA0182	1	1	1	0	1	0	1	0	0	0
KIAA0319L	1	0	1	1	0	0	0	1	0	0
KLHL4	1	1	1	0	1	1	1	1	0	1
MED12	1	1	1	1	1	1	1	0	1	0
PILRB	1	1	1	0	1	0	1	0	0	0
RBPJ	1	0	0	0	0	0	0	0	0	0
SIK1	1	1	1	1	1	0	0	0	0	0
U2AF1	1	1	1	0	1	0	0	0	0	0

DOI: 10.7717/peerj.4838/table-5

Notes:

The genes where the mutations are found are shown in the first column.

As presented in the Supplemental Information 1, the results obtained from the other two algorithms for all the patients are similar. When the estimated haplotypes are compared with some methods that estimate the mutational profiles of tumor subclones, Phylosub proposed in (Jiao et al., 2014) and the manual method proposed in (Schuh et al., 2012), we found that some of the haplotypes, specifically, C1, C2, C3 in CLL003; C3, C4, C5 in CLL077 and C1, C2, C3 in CLL006, carry the same set of unique mutations that are present in distinct genomes of subclones. The complete results of the clonal analyses of these methods are presented in the Supplemental Information 1.

Finally, the results presented so far indicate that each heterogeneous tumor sample is made up of more than two haplotypes: usually a few dominant haplotypes and other minor types. The multiple number of haplotypes in a tumor is an indication of the presence of heterogeneity in the sample.