- Proceedings
- Open Access

# MixClone: a mixture model for inferring tumor subclonal populations

- Yi Li
^{1}and - Xiaohui Xie
^{1, 2, 3}Email author

**16 (Suppl 2)**:S1

https://doi.org/10.1186/1471-2164-16-S2-S1

© Li and Xie; licensee BioMed Central Ltd. 2015

**Published:**21 January 2015

## Abstract

### Background

Tumor genomes are often highly heterogeneous, consisting of genomes from multiple subclonal types. Complete characterization of all subclonal types is a fundamental need in tumor genome analysis. With the advancement of next-generation sequencing, computational methods have recently been developed to infer tumor subclonal populations directly from cancer genome sequencing data. Most of these methods are based on sequence information from somatic point mutations, However, the accuracy of these algorithms depends crucially on the quality of the somatic mutations returned by variant calling algorithms, and usually requires a deep coverage to achieve a reasonable level of accuracy.

### Results

We describe a novel probabilistic mixture model, MixClone, for inferring the cellular prevalences of subclonal populations directly from whole genome sequencing of paired normal-tumor samples. MixClone integrates sequence information of somatic copy number alterations and allele frequencies within a unified probabilistic framework. We demonstrate the utility of the method using both simulated and real cancer sequencing datasets, and show that it significantly outperforms existing methods for inferring tumor subclonal populations. The MixClone package is written in Python and is publicly available at https://github.com/uci-cbcl/MixClone.

### Conclusions

The probabilistic mixture model proposed here provides a new framework for subclonal analysis based on cancer genome sequencing data. By applying the method to both simulated and real cancer sequencing data, we show that integrating sequence information from both somatic copy number alterations and allele frequencies can significantly improve the accuracy of inferring tumor subclonal populations.

## Keywords

- Cancer genomics
- Subclonal inference
- Whole genome sequencing
- Somatic copy number alteration
- Allele frequency
- Mixture model

## Background

Tumor genomes have been shown to present extensive cellular heterogeneity for decades since Nowell's original clonal theory for tumor progression [1]. Identifying tumor subclonal populations is important for both understanding the evolution of tumor cells, and for designing more effective treatments as pre-existing mutations occurring in some subclones could lead to drug resistance [2]. For example, a research in lymphocytic leukemia has shown links between the presences of driver mutations within subclones and adverse clinical outcomes [3].

With the advancement of next-generation sequencing (NGS) and launch of large-scale cancer genome sequencing projects [4], computational methods have recently been developed to infer tumor subclonal populations based on cancer genome sequencing data [5–9].

Most of these methods rely on sequence information from somatic point mutations, such as PyClone [5], EXPANDS [6], PhyloSub [7] and rec-BTP [8]. Methods in this category leverage the cluster pattern of allele frequencies at somatic point mutations to detect distinct subclonal populations. However, as the determination of somatic point mutations is imperfect and the inclusion of false-positives is unavoidable [10], deep sequencing with more than 100X coverage is often required for subclonal inferences with high sensitivity and specificity [5, 7, 8].

Other approaches utilizing the read depth information from genomic segments with somatic copy number alterations (SCNAs) to infer the cellular prevalences of subclonal populations have also been developed, such as THetA [9]. THetA explores all combinations of copy number changes across all segments to infer the most likely collection of subclonal populations [9]. However, with the copy number information alone, THetA suffers from the "identifiability problem", where distinct combinations of tumor purity and ploidy are able to explain the read depth information from SCNAs equally well [9]. Additionally, the running time of THetA scales exponentially with the number of genomic segments [9], and often takes a prohibitively long time to run under certain parameter settings.

In this article, we present a novel probabilistic mixture model, MixClone, to infer the cellular prevalences of subclonal populations. MixClone integrates both read depth information from genomic segments with SCNAs and allele frequency information from heterozygous single-nucleotide polymorphism (SNP) sites within a unified probabilistic framework. Such integrative framework has been shown to significantly improve the accuracy of tumor purity estimation in our previous work [11]. Here, we present that MixClone achieves two major advantages compared to the existing methods that (i) it does not require deep sequencing data, (ii) it resolves the identifiability problem. To demonstrate MixClone's utility, we conducted simulation studies and showed that it outperforms existing methods. We also applied MixClone on a breast cancer sequencing dataset [12], and showed that it was able to discover subclonal events not reported before.

## Methods

In this section, we introduce the generative mixture model of MixClone, which is an extension of our previous work on tumor purity estimation[11]. First, we introduce the notations for input data. Then, we describe the probabilistic models for sequence information of both SCNAs and allele frequencies. Finally, we combine these two types of data into a single likelihood model, and describe an algorithm to solve the model.

### Basic notations

The raw input data for MixClone are two aligned whole genome sequencing read sets of paired normal-tumor samples and a genome segmentation file based on the tumor sample. Following the notations from our previous work [11], we assume the tumor genome has been partitioned into *J* segments. We also assume there are *I*_{
j
} heterozygous SNP sites within segment *j* in the corresponding normal genome, and use (*i*, *j*) to index SNP site *i* within segment *j*. For each SNP site (*i*, *j*) we define the A allele to be the reference allele and the B to be the alternative allele, with respect to the reference genome. We also use a superscript N to denote data from normal samples and superscript T to denote data from tumor samples. Overall, the observed data are summarized in the following notations [13]:

${b}_{i\phantom{\rule{2.77695pt}{0ex}}j}^{N}$ = number of reads mapped to the B allele in the normal sample at site (*i*, *j*).

${d}_{i\phantom{\rule{2.77695pt}{0ex}}j}^{N}$ = reads depth of the normal sample at site (*i*, *j*).

${D}_{j}^{N}$ = total number of reads mapped to segment *j* of the normal sample.

The notations for the observed data from tumor samples are similarly defined, e.g. ${D}_{j}^{T}$ denotes total number of reads mapped to segment *j* of the tumor sample.

### Modeling SCNAs

*j*, we define an allelic configuration

*H*

_{ j }to represent its underlying allele-specific copy number status. For example, if the absolute copy number of segment

*j*is 2, then the compatible allelic configurations are PP, MM and PM, where P and M denotes the paternal and maternal allele of the tumor genome, respectively. Since PP and MM are not distinguishable based on sequence information alone as the reference human genome is not phased, we define the set of all possible allelic configuration as

*K*subclonal populations within the tumor sample, each of which has an associated cellular prevalence

*ϕ*

_{ k }∈ 0[1]. The subclonal type of each segment

*j*is denoted as

*K*possible subclonal populations. Given the allelic configuration

*H*

_{ j }=

*h*and the subclonal type

*Z*

_{ j }=

*k*, the average copy number of segment

*j*within the tumor sample, taking into account the subclonal cellular prevalence

*ϕ*

_{ k }, is

*j*depends on three main factors: 1) its copy number, 2) its total genomic length, and 3) its mappability, which depends on factors such as repetitive sequence and GC content [9]. For each segment

*j*, we associate a coefficient

*θ*

_{ j }to account for the effect of its mappability and genomic length. Thus the expected read counts mapped to segment

*j*, which is denoted as

*λ*

_{ j }, is proportional to ${\stackrel{\u0304}{C}}_{j}{\theta}_{j}$. For example, for segment

*x*and segment

*y*, we have

Because the mappability coefficients (*θ*_{
j
}'s) matter only in a relative sense, we take ${\theta}_{x}\text{/}{\theta}_{y}={D}_{x}^{N}\text{/}{D}_{y}^{N}$, as these segments should have the same sequence properties between the normal and tumor samples.

*λ*

_{ j }, we curate a list of segments which contain no loss of heterozygosity according to their allele frequencies information. Based on the observed number of reads mapped to each segment, we further remove "outlier" segments from the list if their copy numbers are different from the bulk of the segments' copy numbers in the list. Finally, we call the remaining segments in the list as "baseline segments" and denote the set of these segments as

*S*. We assume the allelic configurations of all the baseline segments are PM with copy number

*n*

_{ s }= 2. Other possible allelic configurations for baseline segments, which have equal copy numbers for each allele (e.g.

*φ*, PPMM), are likely to be rare, and currently we do not model them. Then based on

*n*

_{ s }, we specify

*λ*

_{ j }as follows

where ${D}_{s}^{T}$ denotes the number of reads mapped to segment *s* of the tumor sample.

*j*in the tumor sample as a Poisson distribution, given

*H*

_{ j }and

*Z*

_{ j }

Details on curating the baseline segments are given in Supplementary, Additional file 1.

### Modeling allele frequencies

*i*within segment

*j*, we denote its tumor genotype as

*G*

_{ ij }, which is selected from the set of all possible tumor genotypes up to a maximum copy number alteration, e.g.

in which, *ε* ≪ 1 is a small random deviation accounting for general sequencing errors. We choose *E* = 0.01, which is equivalent to a Phred quality of 20 [15].

*G*

_{ ij }=

*g*, the allelic configuration

*H*

_{ j }=

*h*, and the subclonal type

*Z*

_{ j }=

*k*, the average BAF of site (

*i*,

*j*) within the tumor sample, taking into account the subclonal cellular prevalence

*φ*

_{ k }, is

*µ*

_{0}= 0.5 is the BAF of heterozygous SNP sites in the normal sample. Finally, we model the distribution of the B allele count ${b}_{ij}^{T}$ at site (

*i*,

*j*) as a binomial distribution, given

*G*

_{ ij },

*H*

_{ j }and

*Z*

_{ j }

### Combining SCNAs and allele frequencies

Now, we combine sequence information from both SCNAs and heterozygous SNP sites. For all the heterozygous SNP sites within the same segment, their genotypes should be consistent with the underlying allelic configuration of the segment. We model this consistency through a predefined conditional probability ${Q}_{gh}=\mathbb{P}\left({G}_{ij}=g|{H}_{j}=h\right)$. If the genotype *g* is inconsistent with the allelic configuration *h*, e.g. AA is inconsistent with PM, we assign a small probability *σ* as *Q*_{
gh
}, otherwise we assign equal probabilities to genotypes that are consistent with the allelic configuration.

*H*

_{ j }and subclonal type

*Z*

_{ j }, the probability of observing B allele read count ${b}_{ij}^{T}$ at site (

*i*,

*j*) is given as

*H*

_{ j }, the B allele read counts ${\left\{{b}_{ij}^{T}\right\}}_{i=1}^{{I}_{j}}$ at different sites within the same segment

*j*are independent of each other, and are also independent of the total read count ${D}_{j}^{T}$ of the segment. Then, the joint probability of observing the two types of read counts information of segment

*j*is

### Likelihood model

*j*. We then further model the allelic configuration

*H*

_{ j }and the subclonal type

*Z*

_{ j }of segment

*j*as random variables that follow categorical distributions

*ρ*

_{ j }= (

*ρ*

_{j∅}, ⋯,

*ρ*

_{jPPM/PMM}), where ${\rho}_{jh}=\mathbb{P}\left({H}_{j}=h\right)$ is the probability of observing

*h*as the allelic configuration of segment

*j*.

*π*= (

*π*

_{1}, ⋯,

*π*

_{ K }), where ${\pi}_{k}=\mathbb{P}\left({Z}_{j}=k\right)$ is the probability of observing subclonal type

*k*for all the segments. The model parameters Θ is defined as

We use Expectation-Maximization (EM) algorithm [16] to find the maximum likelihood estimation of Θ. The complete details of the EM updates are given in Supplementary, Additional file 1.

### Model selection

*K*. PyClone and PhyloSub use posterior sampling methods to estimate

*K*[5, 7], while THetA requires users to specify

*K*as an input [9]. Since the probabilistic model of MixClone is a generative mixture model, the model complexity and the corresponding log-likelihood increases as

*K*increases. Therefore, we use a criterion based on the increase of the log-likelihood to select

*K*. Practically, Mix-Clone allows the user to specify

*K*. If

*K*is not specified, MixClone runs the mixture model five times with different

*K*in range of 1 to 5. We denote the log-likelihoods under the five different settings as ${\left\{{L}_{K}\right\}}_{K=1}^{5}$, and the total log-likelihood increase as

*/L*

_{1}| < 0.01, which means the ratio of total log-likelihood increase is less than 0.01, MixClone predicts there is no subclonal event in the tumor sample and selects

*K*= 1 as the number of subclonal populations. If |Δ

*/L*

_{1}| ≥ 0.01, MixClone further calculates another quantity

which is the cumulative log-likelihood increase from *K* = 1 to *K* = *i* as a percentage regarding to the total increase Δ. If *δ*_{
i
} ≥ 0.9 and *δ*_{i−1 }< 0.9, MixClone selects *K* = *i* as the number of subclonal populations.

In practice, we suggest users use this criterion as a heuristic guide when analyzing real data, and determine the number of subclonal populations in conjunction with regard to other external information.

### MixClone software package

## Results

In this section, we evaluate the performance of MixClone on both simulated and real datasets and compare its performance with two published algorithms: (i) PyClone, a method based on somatic point mutations, and (ii) THetA, a method based on somatic copy number alterations.

### Results from simulated data

To generate simulation data, we simulated ten sets of NGS reads from chromosome 1 of artificial paired normal-tumor samples, each with 60X coverage. Heterozygous SNP sites from dbSNP [17] were inserted to the reference human genome to create the artificial normal genome. Both heterozygous SNP sites and somatic point mutations from [18] were inserted to the reference human genome to create artificial tumor genomes. Five of the artificial tumor genomes contain two subclonal populations and the other five contain three subclonal populations. Each artificial tumor genome was randomly assigned with segmentations, allelic configurations and subclonal cellular prevalences. We used segmentations based on both ground truth and BIC-seq [19] as the input for MixClone. We used ground truth somatic point mutation sites and copy numbers as the input for PyClone and THetA. Details on how reads were simulated and preprocessed are given in Supplementary, Additional file 1.

### Results from breast cancer sequencing data

We also applied MixClone on a whole-genome breast cancer sequencing dataset [12]. The details on data preprocessing are described in Supplementary, Additional file 1.

Figure 3b shows the five log-likelihoods of MB-116 under different numbers of sub-clonal populations. The magenta, red and yellow curves represent the log-likelihoods corresponding to number 1, 3, and 5, respectively. Because the distance between the magenta and red curves (the cumulative log-likelihood increase from 1 to 3) is greater than 0.9 of the distance between the magenta and yellow curves (the total log-likelihood increase from 1 to 5), MixClone selected *K* = 3 as the number of subclonal populations for MB-116.

^{−4}, which is less than the threshold of 0.01. Therefore, MixClone selected

*K*= 1 as the number of subclonal populations for MB-106. The estimated cellular prevalence of this single population is 83%, which is also consistent with the tumor purities estimated by PyLOH, ABSOLUTE and one result of THetA [11] (Figure 4a).

Besides MB-116, MixClone also detected significant subclonal events in MB-45 and MB-123. Results of MB-45 and MB-123 are given in Supplementary, Additional file 1.

## Discussion

In this article, we demonstrated MixClone's utility using whole genome sequencing data. However, most of the existing cancer genome sequencing data are from exome sequencing. An important future direction is to extend the current methodology to handle the exome sequencing data. Yet, extending MixClone to whole exome sequencing data is not trivial, as reads coverage on targeted exonic regions are no longer randomly distributed due to probe's variable efficiency [21]. Instead of Poisson distribution, using Gaussian distribution to model reads depth ratios between tumor and normal samples might be more appropriate to account for such additional variances, which has been demonstrated in whole exome sequencing based copy number analysis [21].

Another important future direction to extend MixClone is to implement joint analysis based on multiple samples, which is supported by PyClone and PhyloSub [5, 7]. Multiple samples have been obtained for a single heterogeneous tumor tissue both temporally and spatially, and joint analysis based on these samples may reveal additional patterns of the history of tumor progression [5].

Currently, MixClone runs the subclonal analysis five times with different number of subclonal populations in range of 1 to 5 by default. In reality, larger numbers of subclonal populations may coexist within one tumor sample, but in this case some of the populations are very likely to share similar cellular prevalences. Since Mix-Clone defines different subclonal populations based on distinct cellular prevalences, those populations with similar cellular prevalences may not be differentiated by MixClone. To achieve finer resolution of subclonal populations, subclonal lineages information would be necessary to further differentiate each population in addition to cellular prevalences. And phylogenetic methods may be possible solutions to explicitly incorporate subclonal lineages information [7].

## Conclusions

In summary, we have developed a new method for inferring tumor subclonal populations by integrating sequence information gathered from SCNAs and heterozygous SNP sites. We showed that our method outperforms existing ones on simulation data, and applying it to a real breast cancer dataset is able to reveal new subclonal events not discovered before. Compared with existing methods, our method requires no additional deep sequencing of somatic point mutation sites.

## Declarations

### Acknowledgements

The work was partly supported by National Institute of Health grant R01HG006870. The authors would also like to acknowledge dbGaP repository for providing the cancer sequencing datasets. The accession numbers for the breast cancer and prostate cancer datasets are phs000369.v1.p1 and phs000447.v1.p1, respectively.

This article has been published as part of *BMC Genomics* Volume 16 Supplement 2, 2015: Selected articles from the Thirteenth Asia Pacific Bioinformatics Conference (APBC 2015): Genomics. The full contents of the supplement are available online at http://0-www.biomedcentral.com.brum.beds.ac.uk/bmcgenomics/supplements/16/S2

## Authors’ Affiliations

## References

- Nowell PC: The clonal evolution of tumor cell populations. Science. 1976, 194 (4260): 23-28. 10.1126/science.959840.View ArticlePubMedGoogle Scholar
- Garraway LA, Lander ES: Lessons from the cancer genome. Cell. 2013, 153 (1): 17-37. 10.1016/j.cell.2013.03.002.View ArticlePubMedGoogle Scholar
- Landau DA, Carter SL, Stojanov P, McKenna A, Stevenson K, Lawrence MS, Sougnez C, Stewart C, Sivachenko A, Wang L, et al: Evolution and impact of subclonal mutations in chronic lymphocytic leukemia. Cell. 2013, 152 (4): 714-726. 10.1016/j.cell.2013.01.019.PubMed CentralView ArticlePubMedGoogle Scholar
- Hudson TJ, Anderson W, Aretz A, Barker AD, Bell C, Bernabé RR, Bhan M, Calvo F, Eerola I, Gerhard DS, et al: International network of cancer genome projects. Nature. 2010, 464 (7291): 993-998. 10.1038/nature08987.View ArticlePubMedGoogle Scholar
- 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. Nature methods. 2014, 11 (4): 396-398. 10.1038/nmeth.2883.View ArticlePubMedGoogle Scholar
- Andor N, Harness JV, M¨uller S, Mewes HW, Petritsch C: Expands: expanding ploidy and allele frequency on nested subpopulations. Bioinformatics. 2014, 30 (1): 50-60. 10.1093/bioinformatics/btt622.PubMed CentralView ArticlePubMedGoogle Scholar
- 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-10.1186/1471-2105-15-35.PubMed CentralView ArticlePubMedGoogle Scholar
- Hajirasouliha I, Mahmoody A, Raphael BJ: A combinatorial approach for analyzing intra-tumor heterogeneity from high-throughput sequencing data. Bioinformatics. 2014, 30 (12): 78-86. 10.1093/bioinformatics/btu284.View ArticleGoogle Scholar
- Oesper L, Mahmoody A, Raphael BJ: Theta: inferring intra-tumor heterogeneity from high-throughput dna sequencing data. Genome biology. 2013, 14 (7): 80-80. 10.1186/gb-2013-14-7-r80.View ArticleGoogle Scholar
- Roberts ND, Kortschak RD, Parker WT, Schreiber AW, Branford S, Scott HS, Glonek G, Adelson DL: A comparative analysis of algorithms for somatic snv detection in cancer. Bioinformatics. 2013, 29 (18): 2223-2230. 10.1093/bioinformatics/btt375.PubMed CentralView ArticlePubMedGoogle Scholar
- Li Y, Xie X: Deconvolving tumor purity and ploidy by integrating copy number alterations and loss of heterozygosity. Bioinformatics. 2014, 174:View ArticleGoogle Scholar
- Banerji S, Cibulskis K, Rangel-Escareno C, Brown KK, Carter SL, Frederick AM, Lawrence MS, Sivachenko AY, Sougnez C, Zou L, et al: Sequence analysis of mutations and translocations across breast cancer subtypes. Nature. 2012, 486 (7403): 405-409. 10.1038/nature11154.PubMed CentralView ArticlePubMedGoogle Scholar
- Roth A, Ding J, Morin R, Crisan A, Ha G, Giuliany R, Bashashati A, Hirst M, Turashvili G, Oloumi A, et al: Jointsnvmix: a probabilistic model for accurate detection of somatic mutations in normal/tumour paired next-generation sequencing data. Bioinformatics. 2012, 28 (7): 907-913. 10.1093/bioinformatics/bts053.PubMed CentralView ArticlePubMedGoogle Scholar
- Lander ES, Waterman MS: Genomic mapping by fingerprinting random clones: a mathematical analysis. Genomics. 1988, 2 (3): 231-239. 10.1016/0888-7543(88)90007-9.View ArticlePubMedGoogle Scholar
- Ewing B, Green P: Base-calling of automated sequencer traces using phred. ii. error probabilities. Genome research. 1998, 8 (3): 186-194.View ArticlePubMedGoogle Scholar
- Dempster AP, Laird NM, Rubin DB, et al: Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal statistical Society. 1977, 39 (1): 1-38.Google Scholar
- Sherry ST, Ward M-H, Kholodov M, Baker J, Phan L, Smigielski EM, Sirotkin K: dbsnp: the ncbi database of genetic variation. Nucleic acids research. 2001, 29 (1): 308-311. 10.1093/nar/29.1.308.PubMed CentralView ArticlePubMedGoogle Scholar
- Berger MF, Lawrence MS, Demichelis F, Drier Y, Cibulskis K, Sivachenko AY, Sboner A, Esgueva R, Pflueger D, Sougnez C, et al: The genomic complexity of primary human prostate cancer. Nature. 2011, 470 (7333): 214-220. 10.1038/nature09744.PubMed CentralView ArticlePubMedGoogle Scholar
- Xi R, Hadjipanayis AG, Luquette LJ, Lee E, Zhang J, Johnson MD, Muzny DM, Wheeler DA, Gibbs RA, et al: Copy number variation detection in whole-genome sequencing data using the bayesian information criterion. Proceedings of the National Academy of Sciences. 2011, 108 (46): 1128-1136. 10.1073/pnas.1110574108.View ArticleGoogle Scholar
- Carter SL, Cibulskis K, Helman E, McKenna A, Shen H, Zack T, Laird PW, Onofrio RC, Winckler W, Weir BA, et al: Absolute quantification of somatic dna alterations in human cancer. Nature biotechnology. 2012, 30 (5): 413-421. 10.1038/nbt.2203.PubMed CentralView ArticlePubMedGoogle Scholar
- Sathirapongsasuti JF, Lee H, Horst BA, Brunner G, Cochran AJ, Binder S, Quackenbush J, Nelson SF: Exome sequencing-based copy-number variation and loss of heterozygosity detection: Exomecnv. Bioinformatics. 2011, 27 (19): 2648-2654. 10.1093/bioinformatics/btr462.PubMed CentralView ArticlePubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.