Bayesian Inference for Tumor Subclones Accounting for Sequencing and Structural Variants

Tumor samples are heterogeneous. They consist of different subclones that are characterized by differences in DNA nucleotide sequences and copy numbers on multiple loci. Heterogeneity can be measured through the identification of the subclonal copy number and sequence at a selected set of loci. Understanding that the accurate identification of variant allele fractions greatly depends on a precise determination of copy numbers, we develop a Bayesian feature allocation model for jointly calling subclonal copy numbers and the corresponding allele sequences for the same loci. The proposed method utilizes three random matrices, L, Z and w to represent subclonal copy numbers (L), numbers of subclonal variant alleles (Z) and cellular fractions of subclones in samples (w), respectively. The unknown number of subclones implies a random number of columns for these matrices. We use next-generation sequencing data to estimate the subclonal structures through inference on these three matrices. Using simulation studies and a real data analysis, we demonstrate how posterior inference on the subclonal structure is enhanced with the joint modeling of both structure and sequencing variants on subclonal genomes. Software is available at http://compgenome.org/BayClone2.


Biological background and motivation
Understanding tumor heterogeneity (TH) is critical for precise cancer prognosis. Not all tumor cells have the same genome and respond to the same treatment. TH arises when somatic mutations occur in only a fraction of tumor cells, and results in the observed spatial and temporal heterogeneity of tumor samples (Russnes et al., 2011;Greaves and Maley, 2012;Frank and Nowak, 2004;Biesecker and Spinner, 2013;Frank and Nowak, 2003;De, 2011;Bedard et al., 2013;Navin et al., 2011;Ding et al., 2012). In other words, a tumor sample is composed of different subclones of cells with each subclone being defined by a unique genome. Figure 1(a) illustrates this process with a hypothetical case in which accumulation of variants over the lifetime of a tumor gives rise to different subpopulations of tumor cells.
Researchers have recently started to recognize the importance of TH and realize the mistake of treating cancer using a one-size-fits-all approach. Instead, precision medicine now aims to focus on targeted treatment of individual tumors based on their molecular characteristics, including TH.
Rapid progress has been made in the development of computational tools for clonal inference in the past year (Oesper et al., 2013;Miller et al., 2014;Strino et al., 2013;Jiao et al., 2014;Zare et al., 2014). New methods continue to set new and higher standards in the statistical inference for TH that mimic the underlying biology ever more closely.
However, the current literature still lacks effective methods, computational or experimental, for assessing differences between subclonal genomes in terms of both structure variants, such as copy number variants (CNVs), and in terms of sequence variants, such as single nucleotide variants (SNVs). More importantly, current methods lack computational models that could jointly estimate copy numbers and variant allele counts within each subclone. Recent work by Li and Li (2014) adjusts the estimation of subclonal cellular fractions for both CNVs and SNVs, but still stops short of directly inferring subclonal copy numbers or variant allele counts. Figure 1(b) shows a stylized example of DNA-Seq data for a sample taken on day 360 of the process shown in Figure 1(a). The sample is a result of the underlying tumor evolution.
The sample has three tumor subclones. If the sample is sequenced and short reads are mapped, the total number of reads mapped to each locus will be affected by the copy numbers of all the subclones. In Figure 1(b), due to the copy number gains in subclones 2 and 3, we expect that there will be additional reads with sequence A at locus 1 and additional reads with sequence G at locus 3 (both marked by brown letters). In addition, On days 90, 180, and 360, four somatic mutations (represented by red letters) and three somatic copy number gains (represented by brown letters) result in three tumor subclones. (b) Observed short reads (some with variants) are results of heterogeneous subclonal genomes. In particular, the formula at the bottom shows that subclonal alleles are mixed in proportions to produce short reads, which are mapped to different loci.  Figure 2: Three matrices for inference to describe the subclonal structure in Figure 1. L describes the subclonal copy numbers, Z describes the numbers of subclonal variant alleles, and w describes the cellular fractions of subclones.

Life History of a Tumor
variant short reads will be generated due to the subclonal mutations in the sample, such as short reads with the red letters mapped to loci 1 and 2. Using NGS data we aim to recover the subclonal sequences at these loci and cellular fractions at the bottom of Figure   1(b) that explains the true biology in (a). In particular, we aim to provide three matrices as shown in Figure 2 to describe the subclonal genomes and sample heterogeneity. For illustration, Figure 2 fills in the (biological) truth for these three matrices corresponding to the hypothetical tumor heterogeneity described in Figure 1(a). In an actual data analysis, all three matrices are latent and must be estimated.

Model-based Inference for Tumor Heterogeneity
We propose a new class of Bayesian feature allocation models (Broderick et al., 2013) to implement inference on these three matrices. We first construct an integer-valued matrix L to characterize subclonal copy numbers. Each column corresponds to a subclone and rows correspond to loci. We use the column vector c = ( 1c , . . . , Sc ) of integers to represent copy numbers across S loci for subclone c. For example, in Figure 2, sc = 3 for s = 1 and c = 2 since subclone 2 has three alleles at locus 1. As a prior distribution for L, p(L), we will define a finite version of a categorical Indian buffet process (Sengupta, 2013;Sengupta et al., 2015), a new feature allocation model.
Next, we introduce a second integer-valued matrix Z with the same dimensions as L.
We use Z to record SNV's. Denote by z c the c-th column of Z. Conditional on c , z c = (z 1c , . . . , z Sc ), z sc ≤ sc , represents the number of alleles that bear a mutant sequence different from the reference sequence at locus s, s = 1, . . . , S in subclone c. For example, in Figure 2, z sc = 1 for s = 2 and c = 1, indicating that one allele bears a variant sequence. By definition, the number of variant alleles z sc in a subclone cannot be larger than the copy number sc of the subclone, i.e., z sc ≤ sc . Jointly, the two random integer vectors, c and z c describe a subclone and its genetic architecture at the corresponding loci. Lastly, we introduce the w matrix. Each row w t = (w t1 , . . . , w tC ) represents the cellular fractions of the C subclones in each sample (and we will still add an additional subclone c = 0).
The remainder of the paper is organized as follows: Section 2 describes the proposed Bayesian feature allocation model. Section 3 describes simulation studies. Section 4 reports a data analysis for an in-house data set to illustrate intra-tumor heterogeneity. The last section concludes with a final discussion.
Here, M st is the sample copy number that represents an average copy number across subclones. We will formally define and model M st using subclonal copy numbers (L) next. φ t is the expected number of reads in sample t if there were no CNV (the sample copy number equals 2). That is, when M st = 2, the Poisson mean becomes φ t .
Conditional on N st we assume a binomial sampling model for n st conditional on N st ; Here p st is the success probability of observing a read with a variant sequence. It is interpreted as the expected variant allele fractions (VAFs) in the sample. In the following discussion we will represent p st in terms of the underlying matrices L and Z.

Prior
Construction of M st . Let C denote the unknown number of subclones in T samples. We first relate M st to CNV at locus s for sample t. We construct a prior model for M st in two steps, using the notion that each sample is composed of a mixture of C subclones. Let w tc denote the proportion of subclone c, c = 1, . . . , C, in sample t and let sc ∈ {0, 1, 2, . . . , Q} denote the number of copies at locus s in subclone c where Q is a pre-specified maximum number of copies. Here Q is an arbitrary upper bound that is used as a mathematical device rather than having any biological meaning. The event sc = 2 means no copy number variant at locus s in subclone c, sc = 1 indicates one copy loss and sc = 3 indicates one copy gain.
Then the mean number of copies for sample t can be expressed as the weighted sum of the number of copies over C latent subclones where the weight w tc denotes the cellular fractions of subclone c in sample t. We assume The second term C c=1 w tc sc reflects the key assumption of decomposing the sample copy number into a weighted average of subclonal copy numbers. The first term, s0 w t0 denotes the expected copy number from a background subclone to account for potential noise and artifacts in the data, labeled as subclone c = 0. We assume no CNVs at any the locations for the background subclone, that is, s0 = 2 for all s.
Construction of p st and prior on Z. To model the expected VAF of the sample, p st , we construct another feature allocation model linking p st with sc . We introduce an S × C matrix, Z whose entries, z sc ∈ {0, . . . , sc } denote the number z sc ≤ sc of alleles bearing a variant sequence among the total of sc copies at locus s in subclone c. Assume z sc = 0 if sc = 0, and given sc > 0, where DU(·) indicates a discrete uniform distribution.
Next, we write p st in (2) as a ratio between the expected number of variant alleles and the expected sample copy number. In particular, the expected number of variant alleles is a weighted sum of subclonal variant allele counts over (C + 1) latent subclones including the background subclone, and the expected VAF is Similar to the previous argument for (3), the term C c=1 w tc z sc in (5) reflects the assumption that the sample-level variant allele count is a weighted average of subclonal variant allele counts. The first term of the numerator, p 0 z s0 w t0 describes the background subclone and experimental noise. Specifically, we let z s0 = 2 for all s to denote the number of variant alleles in the background subclone. We add a global parameter p 0 to account for artifacts and experimental noise that would produce variant reads even if no subclones were to possess variant alleles. Since p 0 does not depend on s or t, it can be estimated by pooling data from all loci and samples, and does not affect the identifiability of the model. We consider p 0 ∼ Be(a 00 , b 00 ) with a 00 b 00 to inform a small p 0 value a priori. Equation (5)  Prior for w. Next, we introduce a prior distribution for the weights w tc in (3) and (5). The subclones are common for all tumor samples, but the relative weights w tc vary across tumor samples. We assume independent Dirichlet priors as follows. Let θ tc denote an (unscaled) abundance level of subclone c in tissue sample t. We assume θ tc | C iid ∼ Gamma(d, 1) for c = 1, . . . , C and θ t0 iid ∼ Gamma(d 0 , 1). We then define as the relative weight of subclone c in sample t. This is equivalent to for t = 1, . . . , T . Using d 0 < d implies that the background subclone takes a smaller proportion in a sample.
Finally, we complete the model construction with a prior on the unknown number of latent subclones C. We use a geometric distribution, C ∼ Geom(r) where E(C) = 1/r.
Conditional on C, the two latent matrices, L and Z describes C latent tumor subclones that are thought of composing the observed samples and w t provides the relative proportions over those C subclones in sample t. Joint inference on C, L, Z and w t explains tumor heterogeneity.
The construction of the subclones, including the number of subclones, C, the subclonal copy number, sc , and the number of copies having SNV, z sc are latent. The subclones are not directly observed. They are only defined as the components of the assumed mixture that gives rise to the observed CNV and VAFs. The key terms, C c=1 w tc sc in (3) and Lastly, we take account of different average read counts in T samples through φ t . φ t represents the expected read count with two copies in sample t and assume φ t

Posterior Simulation
Let x = (L, Z, θ, φ, π, p 0 ) denote all unknown parameters, where θ = {θ tc } and π = {π cq }. For fixed C such MCMC simulation is straightforward. Gibbs sampling transition probabilities are used to update sc , z sc , π cq and φ t and Metropolis-Hastings transition probabilities are used to update θ and p 0 . It is possible to improve the mixing of the Markov chain by updating all columns in row s of the matrices L and Z jointly by means of a Metropolis-Hastings transition probability that proposes changes in the entire row vector z s and s .
The construction of transition probabilities that involves a change of C is more difficult, since the dimension of L, Z, π and θ changes as C varies. We use the approach proposed in Lee et al. (2014) for posterior simulation in a similar model. We split the data into a small training set (n , N ) with In the implementation we use b st generated from Be (25, 975) for the simulation studies and Be(30, 970) for the lung cancer data. We found that using a random b st worked better than a fixed fraction b across all samples and loci. Let p 1 (x | C) = p(x | N , n , C) denote the posterior distribution under C using the training sample.
We use p 1 in two instances. First, we replace the original prior p( and, second, we use p 1 (·) as proposal distribution q( x | C) = p 1 ( x | C) in a reversible jump (RJ) style transition probability where C is a proposed value of C. The test data is then used to evaluate the acceptance probability. The critical advantage of using the same p 1 (·) as prior and proposal distribution is that the normalization constant cancels out in the Metropolis-Hastings acceptance probability.

Simulation 1
We assess the proposed model via simulation. We generate hypothetical read counts for a set of S = 100 loci in T = 4 hypothetical samples. In the simulation truth, we assume two latent subclones (C TRUE = 2) as well as a background subclone (c = 0) with all SNVs bearing variant sequences with two copies. We use Q = 3. The simulation truth L TRUE is shown in To fit the proposed model, we fix the hyperparameters as r = 0.2, α = 2, γ q = 0.5 for   to the misspecification of 2 and z 2 . This simulation study illustrates that the proposed model reasonably recovers the simulation truth even with a small number of samples when the underlying structure is not complex.
For comparison, we implemented PyClone (Roth et al., 2014) with the same simulated data. We let the normal copy number, the minor parental copy number and the major parental copy number be 2, 0 and 3, respectively, at each locus. PyClone considers copy number changes and estimates the variant allelic prevalence (fraction of clonal population having a mutation) at a locus in a sample. The interpretation of variant allelic prevalences, referred to as "cellular prevalences" in PyClone, is similar to that of p st in the proposed The loci (rows) of the two heatmaps are re-arranged in the same order for easy comparison.
By comparing the two heatmaps, the cellular prevalence estimates under PyClone are close to p TRUE st and lead to a reasonable estimate of a clustering of the loci. However, PyClone does not attempt to construct a description of subclones with genomic variants.

Simulation 2
We carried out a second simulation study with a more complicated subclonal structure. We simulate read counts for a set of S = 100 loci in T = 25 hypothetical samples. In the simulation truth, we assume four latent subclones (C TRUE = 4) as well as a background subclone (c = 0) with all SNVs bearing variant sequences with two copies. We use Q = 3.
The simulation truths, L TRUE and Z TRUE are shown in Figure 6(a) and (b), respectively.
We generated φ TRUE t from Gamma(600, 3) for t = 1, . . . , 25 and then generated w TRUE t as follows. We let a TRUE = (13, 4, 2, 1) and for each t randomly permuted a TRUE . Let a TRUE π denote a random permutation of a TRUE . We generate w TRUE ∼ Dir(0.3, a TRUE π ). That is, the first parameter of the Dirichlet prior for the (C TRUE + 1)-dimensional weight vector was 0.3, and the remaining parameters were a permutation of a TRUE . The weights w TRUE are shown in Figure 6(c). The samples in the rows are rearranged for better display. From Figure 6(c), each sample has all the four subclones with its own cellular fractions, resulting in large heterogeneity within a sample. In addition, the random permutation of a TRUE induces heterogeneity among the samples. We observe that when the underlying subclonal structure is complicate and samples are heterogeneous, larger sample size is needed. In particular, T = 25 which is a large number compared to the typical sample size in real datasets is assumed for this simulation study. Using the assumed L TRUE , Z TRUE and w TRUE and letting  For comparison, we again applied PyClone (Roth et al., 2014) to the same simulated data. We used a similar setting for PyClone as in the previous simulation. Figure 8

Lung Cancer Data
We record whole-exome sequencing for four surgically dissected tumor samples taken from the same patient with lung cancer. We extracted genomic DNA from each tissue and constructed an exome library from these DNA using Agilent SureSelect capture probes. The exome  A total of nearly 115,000 SNVs and small indels were called within the exome coordinates.
The posterior estimates of prevalence and the estimated clustering of the loci are shown in Figure 11(a). The clustering identified five clusters of the loci. The mean prevalences within a locus cluster are similar across samples, which is similar to w in Figure 10(f). Panel (b) of Figure 11 is a heatmap of fractions of reads bearing mutation for each locus and sample.
Again, PyClone provides a reasonable estimate of a loci clustering based on the empirical fractions, but does not provide an inference on subclonal populations.

Conclusions
The proposed approach infers subclonal DNA copy numbers, subclonal variant allele counts and cellular fractions in a biological sample. By jointly modeling CNV and SNV, we provide the desired description of TH based on DNA variations in both, sequence and structure. Such inference will significantly impact downstream treatment of individual tumors, ultimately allowing personalized prognosis. For example, tumor with large proportions of cells bearing somatic mutations on tumor suppressor genes should be treated differently from those that have no or a small proportion of such cells. In addition, metastatic or recurrent tumors may possess very different compositions of cellular genomes and should be treated differently.
Inference on TH can be exploited for improved treatment strategies for relapsed cancer patients, and can spark significant improvement in cancer treatment in practice.
A number of extensions are possible for the present model. For example, sometimes additional sources of information on CNVs such as a SNP array may be available. We then extend the proposed model to incorporate this information into the modeling of L. Another meaningful extension is to cluster patients on the basis of the imputed TH, that is, we link a random partition and a feature allocation model. This extension may help clinicians assign different treatment strategies, and be the basis of adaptive clinical trial designs.
Inference for TH is a critical gap in the current literature. The ability to precisely break down a tumor into a set of subclones with distinct genetics would provide the opportunity for breakthroughs in cancer treatment by facilitating individualized treatment of the tumor that exploits TH. It would open the door for cocktail type of combinational treatments, with each treatment targeting a specific tumor subclone based on its genetic characteristics. We believe that the proposed model may provide a integrated view on subclones to explain TH that remains a mystery to scientists so far.