Tumor Copy Number Deconvolution Integrating Bulk and Single-Cell Sequencing Data

2.1. NMF deconvolution model

As in previous work of Schwartz and Shackney (2010), we formalize the generic problem of genomic deconvolution in terms of a mixed membership model, but here relating bulk and SCS data. We focus here specifically on deconvolution of copy number data, as in Tolliver et al. (2010) and Zaccaria et al. (2018), which we assume is profiled on a set of m genomic regions. In the pure deconvolution problem, we assume a set of n bulk samples, which might correspond to measurements from distinct tumor sites or regions in one patient. These bulk samples are collectively encoded in an $m\times n$ matrix B, where element b_ij corresponds to the mean copy number of locus i in sample j. Our goal is to identify an $m\times k$ matrix of mixture components C, representing copy numbers of inferred common clones, and a $k\times n$ matrix of mixture fractions, F, describing the degree to which each column of C is represented in each column of B. B is presumed to be approximated by the product of C and F. We seek to minimize the deviation between B and C × F by some measure, such as the Frobenius norm. With the additional constraints that B, C, and F are non-negative, the problem is known as NMF (Wang and Zhang, 2013). More formally, we seek the following:

where $\left|\right|\cdot |{|}_{Fr}$ is the Frobenius norm of the matrix, $\left|\right|B-CF|{|}_{Fr}=\sqrt{{\sum }_{i=1}^m{\sum }_{j=1}^n{\left(b_{ij}-{\sum }_{\ell =1}^k{c}_{i\ell }\cdot f_{\ell j}\right)}^2}$, subject to the constraints $f_{\ell j}\ge 0,\forall \ell \in \left\{1,\dots ,k\right\},j\in \left\{1,\dots ,n\right\};{\sum }_{\ell =1}^k{f}_{\ell j}=1,\forall j\in \left\{1,\dots ,n\right\};c_{i\ell }\in {\mathbb{N}}_0,\forall i\in \left\{1,\dots ,m\right\},\ell \in \left\{1,\dots ,k\right\}.$

This optimization problem is nonconvex, but prior work showed that the Euclidean distance between B and CF is nonincreasing under the following multiplicative update rules in Lee and Seung (2001):

providing formulas for iterative local optimization by fixing C or F on alternate steps. In practice, we modify this process heuristically to renormalize columns of F after each iteration to ensure they add to 1. Since this heuristic might undermine the guarantee of monotonicity, we manually verify that $\left|\right|B-CF|{|}_{Fr}$ decreases on each iteration, terminating the optimization if it fails to yield continuing improvements. More details are provided in Section 5.

Figure 1 provides an illustrative example of the deconvolution model. Suppose we have a possible B, C, and F. The two data points B₁ and B₂ represent bulk tumor samples combining three mixture components C₁, C₂, and C₃. For ease of illustration, we assume data are assayed on the copy numbers of just two genomic loci, G₁ and G₂. The matrix B represents the average copy numbers of G₁ and G₂ in the bulk tumor samples B₁ and B₂. In each component of C, the copy numbers should be integers, but since the bulk tumors are weighted mixtures of components, the values in B need not be integers. The matrix F represents the fractional weights used to generate B₁ and B₂ from the pure components in C. For example, the first column of F indicates that B₁ is a mixture of equal parts of C₁ and C₂. This relationship can be expressed via matrix multiplication, $B=CF$, as shown in the right part of Figure 1.

2.2. Extending NMF with SCS data

The multiplicative update algorithm is a standard method for the pure NMF optimization problem, provided the number of samples n is large compared with the intrinsic dimension k of the mixture. We would, however, expect it to perform poorly for our problem, in part, because real tumor data generally include few samples per patient and, in part, because deconvolution of copy numbers is an underdetermined problem. We sought to improve the optimization by biasing the objective function to favor inferred clones similar to the observed SCS data via an auxiliary penalty in the objective function, similar to the approach of Malikic et al. (2017, 2018). Intuitively, we assume that the inferred clones (C) should be closely related to one or more of the observed single cells, which we call observed cell components ($C^{\left(observed\right)}$). While any given single cell may not exactly match a consensus clone, we propose that the method will be able to approximately infer mixture components reflecting dominant clones by balancing quality of deconvolution against similarity to observed single cells. We quantify this intuition using the Euclidean distance between the inferred clones and observed cells, introducing a regularization parameter $\alpha$ to balance the weight of this penalty relative to the prior cost based on deconvolution quality. The resulting combined objective appears as Equation (2):

which we optimize subject to the constraints $f_{\ell j}\ge 0,\forall \ell \in \left\{1,\dots ,k\right\},j\in \left\{1,\dots ,n\right\};$ ${\sum }_{\ell =1}^k{f}_{\ell j}=1,\forall j\in \left\{1,\dots ,n\right\};c_{i\ell }\in {\mathbb{N}}_0,\forall i\in \left\{1,\dots ,m\right\},\ell \in \left\{1,\dots ,k\right\}$.

We solve for the revised model through an extension of the iterative update algorithm in Lee and Seung (2001) and Berry et al. (2007):

adding the constraints on C and F to the update rules in Thurau et al. (2009). We further heuristically improve on the standard practice of random initialization by initializing the cell component matrix C with true SCS data. Pseudocode for the complete algorithm is provided in Section 5 as Algorithm 1. Collectively, these additions to the pure NMF iterative update algorithm constitute our first approach to integrating SCS data for improved deconvolution of CNAs from bulk DNA-seq, which we dub our phylogeny-free method.

2.3. Extending the NMF model with a single-cell phylogeny objective

We next developed an alternative phylogeny-based approach, seeking to deconvolve the bulk data into clonal subpopulations while simultaneously inferring a phylogeny on those deconvolved clones, similar to the SNV PHiSCS method of Malikic et al. (2018). Intuitively, evolutionary distance provides a more biologically motivated measure of what we mean in asserting that inferred single cells should be similar to observed single cells. As with the phylogeny-free method, we would expect that any small sample of single cells will not exactly reflect the spectrum of dominant clones, but that the method will be able to approximately infer dominant clones by balancing deconvolution quality against evolutionary distance of mixture components to observed single cells. This approach trades off a more principled measure of solution quality for a harder optimization problem.

We quantify phylogenetic distance as the minimum over evolutionary trees incorporating both observed single cells and inferred clones of the L₁ distance between copy number vectors describing each tree edge. Let $C^{\ast }=\left[C,C^{\left(observed\right)}\right]$ be an $m\times k^{\ast }$ matrix consisting of columns representing inferred clonal copy numbers followed by columns representing the copy numbers of the observed cells. Let $c_u^{\ast }$ denote column u of $C^{\ast }$. We introduce a $k^{\ast }\times k^{\ast }$ matrix of binary variables S. A value of $s_{uv}=1$ indicates the existence of a directed edge from node u to node v, and a value $s_{uv}=0$ indicates the absence of such a edge; we set $s_{uu}=0$ to avoid self-loops. In other words, S is an adjacency matrix for a directed graph; in the full formulation [Appendix 1 (Supplementary Methods)], we introduce constraints that ensure the graph is a tree. We define our measure of tree cost to be

Intuitively, $J\left(S,C,C^{\left(observed\right)}\right)$ is a form of minimum evolution model on a phylogeny defined by S. While there are more sophisticated and realistic models for CNA distance (Chowdhury et al., 2014, 2015; El-Kebir et al. 2017), we favored L₁ distance here as a tractable approximation easily incorporated into the overall integer linear programming (ILP) framework. Similarly, while there are now a number of sophisticated methods available specifically for phylogenetics of single-cell sequences (c.f., Kuipers et al., 2017), these are largely focused on SNV rather than CNA phylogenetics (Jahn et al., 2016; Ross and Markowetz, 2016; Zafar et al., 2017) with limited exceptions (Wang et al., 2014; Subramanian and Schwartz, 2015).

More specifically, we modify the NMF objective function as follows:

where $\left|\right|B-CF|{|}_1={\sum }_{i=1}^m{\sum }_{j=1}^n\left|b_{ij}-{\sum }_{\ell =1}^k{c}_{i\ell }\cdot f_{\ell j}\right|$, and $\beta$ is a regularization parameter to balance deconvolution quality against parsimony of the evolutionary model. The norm $\left|\right|\cdot |{|}_1$ is the element-wise L₁ matrix norm, that is, the sum of the absolute values of matrix elements, rather than the induced L₁ matrix norm for which the same notation is sometimes used. These are optimized subject to the same constraints as in the previous formulations: $f_{\ell j}\ge 0,\forall \ell \in \left\{1,\dots ,k\right\},j\in \left\{1,\dots ,n\right\};{\sum }_{\ell =1}^k{f}_{\ell j}=1,\forall j\in \left\{1,\dots ,n\right\};c_{i\ell }\in \mathbb{N},\forall i\in \left\{1,\dots ,m\right\},\ell \in \left\{1,\dots ,k\right\}.$

The discrete tree optimization term lacks an analytic expression and hence does not lend itself to the prior iterative update strategy. We therefore use a different computational strategy based on ILP to replace the linear algebra steps of the method in Lee and Seung (2001), similar to other recent works in joint deconvolution and phylogenetics in Zaccaria et al. (2018) and Eaton et al. (2018).

For this optimization problem, we use an iterative coordinate descent approach. There are three sets of variables over which to optimize: the weight matrix F, the tree structure S, and the inferred copy numbers C. We solve for variables F, S, and C alternately, in this order, while holding all other variables as constant. The iterative coordinate descent continues until the decrease between successive values of C falls below some threshold. To initialize C, we used observed single-cell data. Whenever two of the three sets of variables are held constant, the resulting optimization problems can each be expressed as either a linear program (LP) or an ILP.

When certain subsets of the variables are fixed, the resulting LP or ILP may be simplified. When solving for F with fixed values of S and C, the term $J\left(S,C,C^{\left(observed\right)}\right)$ is constant and the value of S is irrelevant. Similarly, when solving for S for fixed values of F and C, the term $\left|\right|B-CF|{|}_1$ is constant and therefore F is irrelevant. The optimal value of C for fixed values of F and S, however, depends both on F and S. In the limit of using no single-cell data, our problem statement and method are similar to that of Zaccaria et al. (2018) for incorporating tree mixtures into purely bulk CNA deconvolution.

We solve for S via an ILP that uses a flow model to constrain solutions to a minimum evolution tree, adapting a similar ILP method originally developed for finding maximum parsimony character-based phylogenies (Sridhar et al., 2007). Intuitively, the model forces a tree structure by setting up a flow from an arbitrary root to each other clone in the tree and minimizing the cost of edges needed to accommodate all such flows. The full ILP is described in Section 5.

2.4. Validation via observed single-cell data

To validate the method, we require bulk data for which clone copy number vectors and frequencies are known. As this is unavailable for any real data set, we use semisimulated data generated from copy number calls (Baslan et al., 2012) from real SCS data from two human glioblastoma cases of Wu et al. (2016). The full single-cell data set consists of low-depth SCS DNA-seq used to establish mean copy numbers at 9934 genomic positions throughout the genome, at intervals of ∼40 kbp. Each tumor was subdivided into three regions (i.e., samples), with each single cell labeled by its region (1, 2, or 3) of origin. We used these true SCS CNA data to generate a series of synthetic bulk data sets, simulating either one, two, or three bulk samples from each region for a total of three, six, or nine bulk samples per trial.

Each simulated sample is generated by sampling two dominant cells from a region to represent major clones, 23 other cells from the same region to represent minor clones, and 50 cells from the other regions to represent contamination, which are mixed with Dirichlet-sampled proportions with weight parameters for dominant, minor, and contaminant clones in the ratio 10 to 0.1 to 0.01. We then assessed our ability to deconvolve the bulk data across a range of regularization parameter values and random replicates of the chosen single cells. We assessed accuracy by the fraction of genomic positions assigned correct copy number and by the root mean square deviation (RMSD) between true and inferred clones and mixture fractions. Figure 2 summarizes the overall experimental design, which is described in more detail in Section 5.4. This design treats observed SCS as the ground truth, allowing us to ignore the problem of doublet cells that typically must be addressed with SCS data. We would normally require that likely doublets be removed from SCS data in preprocessing before applying our method. This design also does not explicitly include calling CNA markers on bulk data, itself a hard problem that would need to be performed in preprocessing before applying our method.

We were unable to identify any competitive tool for bulk deconvolution of purely CNA data applicable to small numbers of bulk samples and for which software is publicly available. We therefore compare our methods to standard NMF, as implemented by our code with zero regularization parameters.

2.5. Real bulk sequencing data

In addition to the semisynthetic data, we applied our methods to true bulk tumor data from the same human glioblastoma cases described in Section 2.4. One bulk sample was available for each region for each case. Copy number estimation was performed by the Complete Genomics Standard pipeline, developed by the Beijing Genomics Institute (BGI), to establish mean copy numbers for 931 regions throughout the genome. Single-cell CNA estimation had been performed on the same tumor regions to obtain estimated copy numbers at 9934 genomic locations. We treated the most extreme locations on each of the 24 chromosomes (includes Y) as chromosome ends, and hence, we interpret the data as copy numbers at 9910 ( = 9934–24) intervals. These SCS copy number intervals were always smaller than the genomic intervals obtained from the bulk data.

Because the intervals in the SCS and bulk data did not have their endpoints in common, we subdivided the bulk data into 9910 intervals aligned to the same endpoints as in the SCS intervals by applying two rules: (1) for each SCS interval that was wholly contained in a bulk interval, we created a new interval that had the same coordinates as the SCS interval, and assigned this new interval the copy number of the bulk interval; (2) for each SCS interval that overlapped two bulk intervals, A and B, we created a new interval for the bulk data, again having the same endpoints as the SCS interval, but with a copy number that is a weighted sum of the copy number of interval A and the copy number of interval B. The weights were chosen to be proportional to the number of nucleotides in the overlap between the new interval and A and B. This process yielded SCS and bulk data with copy numbers assigned to common genomic intervals (Section 5.7). We applied our deconvolution and phylogeny methods as with the semisynthetic data to the collection of three bulk genomic samples corresponding to the three tumor regions.

2.6. Implementation

The methods described in Section 2 (Methods) and refined below were all implemented in Python3, using Gurobi. One practical change from the formulation above is that we replaced the theoretical $f_{\ell j}\ge 0$ with $f_{\ell j}>10^{-4}$ to avoid having the f values trapped at 0. The observed human subjects' data cannot be redistributed, but code for the methods is available along with artificial data on GitHub (https://github.com/CMUSchwartzLab/SCS_deconvolution).

3.1. Phylogeny-free method

We first assessed the accuracy of the phylogeny-free method relative to pure NMF and simple heuristic improvements. Figure 3 provides a summary of accuracy and RMSD for inference of true SCS components via the method of Section 2.2 for variations in the number of tumor samples (3, 6, 9) and regularization parameter $\alpha \left(0-1\right)$ over 40 replicates per condition. To provide a baseline for comparison, each plot provides equivalent accuracy measures for NMF (Lee and Seung, 2001, i.e., Algorithm 1 with $\alpha =0$) with random initial integer valued C (red dashed line in Fig. 3) and with the proposed solution that all copy numbers have the normal value of 2, which we call the “all-diploid baseline” (black dashed line in Figs. 3 and 12). In each case, the bulk data are simulated from $k\prime =6$ fundamental cell components (2 out of a random 25 cells selected in each region).

Pure NMF with random initialization performed poorly, which is unsurprising since NMF on CNA data is an underdetermined problem, although the simple heuristic of biasing the search toward biologically plausible solutions by initializing with real SCS data improves accuracy. Bringing true SCS data into the objective function yielded modest improvements in accuracy over using SCS data solely for initialization for at least some values of the regularization parameter. The phylogeny-free method with $\alpha =0$ corresponds to pure NMF initialized with true SCS data, and this performed slightly worse than the all-diploid baseline solution. Modestly increasing $\alpha$ led to some improvement in accuracy, but above some value, $\alpha$ put too much weight on similarity to the observed SCS data and too little weight on quality of the deconvolution, giving worse overall results. The best value of $\alpha$ depended on sample size, which we attribute again to NMF being underdetermined if the number of desired components is larger than the number of samples. The plots suggest that the method is fairly robust to $\alpha$ if the number of samples exceeds the intrinsic dimension of the data (six), but that SCS data can overcome that limit for small numbers of samples with a well-tuned regularization term. Additional results in Section 6 show minimal additional improvement even with unrealistically large sample sizes (Fig. 13), and also show that the performance is consistent across individual inferred clones (Fig. 14).

Figure 4 provides an illustrative example of performance for a single selected clone inferred from three, six, or nine samples, intended to demonstrate kinds of errors the method tends to produce. We chose the one inferred clone with smallest RMSD for each sample size to simplify visual inspection. We see that at least in these high-quality cases, the distributions of copy numbers are similar for the inferred and true clones. For loci at or just above diploid, the modified NMF can usually infer the exact copy number. Where errors occur, they tend to be in loci with large (5–10) or smaller copy numbers (0–1).

3.2. Phylogeny-based method

We next examined results of the phylogeny-based method of Section 2.3 under the same conditions used to assess the phylogeny-free method. Figure 5 summarizes average accuracy and RMSD as a function of regularization parameter $\beta$. The figure compares the results of pure NMF with the all-diploid baseline. Setting $\beta =0$ provides poor performance, substantially below the all-diploid baseline solution. Making $\beta =0$ for Figure 5 represents the same optimization problem as $\alpha =0$ for Figure 3, but solved by the coordinate descent method we developed to accommodate the ILP phylogeny objective rather than by the modified iterative update algorithm with the simpler L₂ objective. Figure 5 thus suggests that the new coordinate descent method is less effective at pure NMF than is the prior iterative update algorithm. Despite that observation, the results on $\beta \ge 0.2$ show substantially better accuracy than was achieved by pure NMF or the phylogeny-free algorithm. Furthermore, the results appear robust to variation in $\beta$ across the range examined. Results in Section 6 distinguishing accuracy across cells (Fig. 15) support the robustness of the phylogeny-based method to a range of $\beta$ values in cell-to-cell inferences.

Figure 6 shows copy numbers for a single minimum-RMSD pair for inferred and true clones for each number of samples, again to visualize the nature of inference errors. The results again show exact fitting for most loci, as well as better fitting for both large (5–10) and small (0–1) copy numbers than the phylogeny-free method of Figure 4. There is no evident pattern to the smaller number of errors that do occur for the phylogeny-based versus phylogeny-free method, which are observed for a range of low and high copy number values.

Figure 7 compares the two methods at their optimal regularization parameters for three, six, and nine tumor samples. The phylogeny-based method outperforms the phylogeny-free method in accuracy and copy number RMSD in all cases. It is slightly better in mixture fraction RMSD for three samples, but worse for six and nine samples. Figure 16 in Section 6 shows comparative performance of the two methods in individual cell components. Given the poorer performance at pure NMF of the phylogeny-based method's algorithm versus the phylogeny-free method's algorithm, we tentatively attribute the phylogeny-based method's better overall performance to better evolutionary distance estimates and not to a better optimization algorithm.

The phylogeny-based method also provides as output the phylogeny. While we cannot exhaustively show trees across all replicates, we provide three representative examples in Figure 8. Because we use true SCS data to generate our synthetic mixtures, we do not know the full ground truth trees for the data and do not attach any biological meaning to the inferred trees. We can partially validate correctness of the trees using the fact that the cells were gathered from distinct tumor regions, and while we would not expect clonal ancestry to segregate perfectly by region, we should see a trend toward closer evolutionary relationship among cells in spatial proximity. We tested whether pairs of cells from distinct regions cluster together in disjoint subtrees (a kind of partial-information quartet distance); we found that a significant majority of pairs-of-pairs do (79% for three-sample data, 74% each for six- or nine-sample data) providing some support for the biological relevance of the trees.

3.3. Application to real data

We then applied our method to true bulk data from the same glioblastoma data set used in the semisimulated validation. For both patients (GBM07 and GBM33) and each region of the tumor, a single run of bulk sequencing and separate SCS were performed to estimate the copy number profile. We focus on the analysis on the patient GBM07 by using the phylogeny-based method since that method performed best with semisimulated data (Section 3.2).

For bulk data, we do not know the ground truth composition of clones, but may compare our deconvolution results to genomic alterations that have previously been reported in glioblastoma (Fig. 9, detailed analysis methods in Section 5.7). A total of 92 genes have been reported to be frequently altered in glioblastoma (Brennan et al., 2013; Körber et al., 2019), and 89 of these 92 genes were found to have altered copy numbers in the mixture components (clones) inferred by our methods. Among these known mutated genes, EGFR, CDK4, and PDGFRA have been reported to have significant recurrent amplifications, while CDKN2A/B and PTEN have been reported to have significant recurrent deletions according to the Cancer Genome Atlas Research Network et al. (2008).

Our analysis showed similar results. We observed copy number gains of PDGFRA in all but two components (Fig. 9C, blue text), and this gene generally exhibited substantial copy number gain in bulk sequencing (Section 6.6). EGFR was gained in all but one component. PTEN showed copy number loss in eight components and CDKNA2/B was lost in 10 components (Fig. 9C, orange text). We also found ERBB2 gained in three components and RB1 hemizygous or homozygous deletion in at least some glioblastoma samples (Cancer Genome Atlas Research Network et al., 2008; Abou-El-Ardat et al., 2017). These results are consistent with CNAs influencing the RTK/PI3K and RB signaling pathways frequently altered in glioblastoma (Fig. 9B). However, we could not identify CNAs affecting the p53 signaling pathway. We cannot tell from CNA inferences alone whether or not TP53 activity was normal in this particular patient; there are several common TP53-inactivating SNVs that would not show up in a CNA analysis (Magali et al., 2010). Figure 9C shows a more complete profile of inferred gains and losses in potential glioblastoma driver genes.

At the chromosome level, chromosomes 7, 9p, and 10 (10q) frequently display alterations in glioblastoma (Crespo et al., 2011; Davis et al., 2015; Abou-El-Ardat et al., 2017). We therefore looked at the copy numbers at a coarser scale for each inferred component, following the rules in Section 5.7 to decide whether or not the chromosome has a gain or loss. We found that in 11 components, almost all of the genomic loci showed gains in chromosome 7, although some specific loci did not display gain in some components, for example, component 10. All components show loss in chromosome 9p and eight components show loss in chromosome 10. Again, these observations are consistent with aneuploidy or large structural variants in these chromosomes as has been previously observed. Those alterations contain some key genes, for example, EFGR and MET are located on chromosome 7, CDKNA2/B on chromosome 9p, and PTEN on chromosome 10.

Abou-El-Ardat et al. (2017) and Körber et al. (2019) suggested that gain of chromosome 7, loss of chromosome 9p, and loss of chromosome 10 are likely to be early events in glioblastoma, followed by the focal loss of CDKN2A/B, and then later events such as EGFR and PDGFR amplification. Our inferred phylogeny (Fig. 9C) is largely consistent with, although more complicated than, this coarse model. The inferred phylogeny shows substantial branching rather than a simple linear structure, but the relative order of different components reflects early gain of chromosome 7 and loss of chromosome 9p. The loss of chromosome 10 is also inferred to be an early event, but this event could not be detected in all components. Loss of chromosome 10 is not apparent in all bulk samples, either, and the average copy number from the raw SCS data only shows slight loss at chromosome 10 (Section 6.6). The phylogeny suggests reversion of the chromosome 10 loss in some lineages, although that might be an error due to the difficulty of accurately calling CNAs that involve small, subclonal changes in copy number. Nevertheless, our inferred clonal mixture fractions (Fig. 9A) indicate that the first component (C0) has a large proportion in all three regions of the tumor sample, again consistent with the notion that gain of chromosome 7, and loss of chromosome 9p and chromosome 10 are early events in glioblastoma evolution.

We would expect that the model of early versus late events suggested by the prior work is a coarse consensus of what is a stochastic process for any individual tumor. It would require larger cohorts to test whether deviations reflect the variability of the process patient-to-patient or some systematic discrepancy in the overall model.

Besides the qualitative concordance with previous work, we also find some potential driver genes are gained in subsets of components (e.g., KIT, KDR, and SIGLEC8). Such subclonal events may support tumor survival or progression and thus are possible therapeutic targets for inhibitors (Gomes et al., 2007; Nobusawa et al., 2011; Kiwamoto et al., 2012; Pearson and Regad, 2017). Gain of chromosome 19 was found in two components, an event, previously identified in brain tumor-initiating cells (Davis et al., 2015), that in combination with chromosome 7 gain may help a tumor survive radiation therapy (Huhn et al., 1999). We did not see evidence of loss of heterozygosity on chromosome 19, which has been frequently found in secondary glioblastoma in Nakamura et al. (2000) and Ohgaki and Kleihues (2007). The gain of chromosome 17 was observed in two components, potentially a mechanism for amplification of the ERBB2 oncogene. Another interesting gene is ATRX, which is located on chromosome X. This gene has a loss in some components. Mutations of ATRX are mostly found in younger patients diagnosed with glioblastoma but rare in adults, and loss of the orthologous Atrx accelerated glioblastoma growth in mouse (Koschmann et al., 2016; Nandakumar et al., 2017).

The mixture fraction reconstructions together with the phylogeny suggest an overall picture of early evolution of clone 0, which seeded all three tumor regions, followed by diversification most prominently into clones 1 and 2, which are dominant in two different tumor regions and represent two distinct early divergence events in the phylogeny. Most other clones are inferred to be comparatively rare. Given the uncertainty in mixture fraction analysis, we would not attach any significance to inferences of relatively rare clones in a given tumor region. We would thus suggest that the model is largely consistent with distinct tumor regions evolving independently without clear evidence for migration or reseeding between regions.

5.1. Phylogeny-free method for integrating single-cell sequencing data into non-negative matrix factorization

We solve for the phylogeny-free variant of single-cell sequencing (SCS)-assisted non-negative matrix factorization (NMF) [Eq. (2)] via an implementation of the iterative update algorithm, for which we present pseudocode as Algorithm 1. The core of the algorithm consists of the modified update rules from Section 2.2. Additional heuristic modifications are used to enforce non-negativity and integrality of solutions and make use of SCS data to bias initialization toward biologically plausible solutions, as well as to better handle possible numerical errors arising from finite machine precision.

Algorithm 1: Modified Multiplicative Update Algorithm for NMF

Open in a new tab

5.2. Integer linear programming for phylogeny inference in the phylogeny-based method

The major change in the phylogeny-based method of Section 2.3 is the introduction of a minimum-evolution phylogeny cost, $J\left(S,C,C^{\left(observed\right)}\right)$, to the objective function [Eq. (4)]. We solve for the phylogeny on each pass of the coordinate descent algorithm via a multicommodity flow integer linear programming (ILP) in Sridhar et al. (2007). Let the vertex set $T=\left\{1,\dots ,k^{\ast }\right\}$ represent the set of all cells in $C^{\ast }$. Let r be the unique predetermined root of T, which in the present practice is a purely diploid node. Furthermore, let w_u,v be the L₁ distance between the copy number vectors corresponding to nodes $u,v\in T$. For $t,u,v\in T$, introduce the binary variables $g_{v,u}^t$ representing the amount of flow along edge $\left(u,v\right)$ with destination $t\in T$. The full ILP is then as follows:

Intuitively, the method defines a flow from a single root to every other vertex and requires the graph T to contain edges $\left(u,v\right)$, indicated by $s_{uv}=1$, such that all such flows can be accommodated. This forces the graph to be connected. We can further establish that the resulting graph is acyclic. For purposes of contradiction, assume the optimal T contains a cycle. We can then remove any edge $\left(u,v\right)$ on the cycle and reroute any flow using that edge through the cycle in the other direction, reducing the cost of the tree by $w_{uv}$. Since $w_{uv}$ is a non-negative L₁ distance, then this must reduce the cost of the tree provided $u\ne v$, showing that T was nonoptimal. This establishes by contradiction that the optimal T is connected and acyclic, that is, a tree. It is specifically a tree of minimum cost over the complete graph of observed single cells and inferred clones implied by C and $C^{\left(observed\right)}$.

At any stage of the algorithm, if there are multiple minimum-cost solutions for F, S, or C, any solution might be chosen.

5.3. Generation of single-cell sequence glioblastoma data

Data for this analysis are provided from a study of single-cell genomics in two glioblastoma patients in Wu et al. (2016). Each patient's primary tumor was divided into three tumor regions, with 59–82 single cells extracted from each region for sequencing, for a total of 448 cells. Nuclei of 432 of these cells were amplified by multiple displacement amplification and sequenced to a coverage of 0.17 × . We used 393 among the 410 cells that passed quality control. These cells were called for CNVs by modified variable binning in Baslan et al. (2012). The resulting calls provided the input for single-cell analysis and for construction of semisynthetic bulk data in the present work.

5.4. Generating semisynthetic data from SCS samples

This section provides additional details on the generation of semisynthetic data from SCS samples for use in validating the methods. For each simulated tumor, we generate a set of experiments in which we simulated one, two, or three bulk samples from each of three tumor regions, for a total of three, six, or nine bulk samples, respectively. To generate the simulated bulk samples, we randomly chose 25 single cells from each region, for a total of 75 cells. These 75 cells define the copy number vectors that make a nonzero contribution to any of the simulated bulk samples. For each region, we chose two cells from among the 25 to represent codominant clones. We refer to the two chosen cells as dominant cells. To model the noisy nature of bulk tumor data, the two dominant cells in the region, the 23 remaining cells from the same region, and the 50 cells from the other two regions all contribute to the copy number of each bulk sample, but at different mixture fractions. The two dominant cells from each region make by far the greatest expected contribution to the simulated bulk samples for that region. This design is intended to approximate clonal structure of real tumor samples, where one might see a small number of dominant clones and a long tail of rarer cell populations (Heselmeyer-Haddad et al., 2012). We assess the methods on their ability to infer the dominant clones, with the rare clones effectively serving as a source of noise in the analysis.

To generate random mixture fractions for the sampled cells, we sampled from Dirichlet distributions. Dirichlet distributions are the conjugate priors of multinomial distributions. Thus, a Dirichlet distribution is a distribution, with vector valued parameter γ, of probabilities for the multinomial distribution—in other words, of mixture fractions for copy numbers obtained by sampling. For each region, each cell i of the two dominant cells was assigned ${\gamma }_i=10$, each cell j of the 23 other selected cells from the same region was assigned ${\gamma }_j=0.1$, and each cell $\ell$ of the other 50 cells was assigned ${\gamma }_{\ell }=0.01$.

Figure 10 shows the overall experimental design, including different regions from which single cells are collected and DNA sequenced. We extracted the copy number for each genomic locus from the SCS results to compose the cell matrix, as shown in the heatmap in Figure 10.

More formally, we chose k total single-cell samples, $\widehat{k}$ dominant cell samples per region, and $\widehat{n}$ simulated bulk samples per region. Here, $k=75$, $\widehat{k}=2$, and $\widehat{n}$ is 1, 2, or 3. We drew $\widehat{n}$ column vectors of mixture fractions for each region from a Dirichlet distribution using the parameters $\gamma$, assigned as described above, for a total of $n=3\widehat{n}$ columns of mixture fractions. We use these n columns to form the $k\times n$ matrix $F^{\left(sel\right)}$. The $m=9934$ and $k=75$ columns of copy numbers from the single-cell data are used to form the $m\times k$ matrix $C^{\left(sel\right)},$ and the $m\times n$ simulated matrix B is the product $B=C^{\left(sel\right)}{F}^{\left(sel\right)}$. The entries of B are real values, not rounded to integers.

We also selected a number $k_{observed}=3\widehat{k}$ of single cells to form the matrix $C^{\left(observed\right)}$ for use in some algorithms. The number of single cells selected varied, but in each case, the same number of cells was selected from each region, and the cells were selected from those that were not chosen to be among the 75 used to simulate the bulk tumor.

Our goal is to recover the $k\prime =3\widehat{k}$ (here $k\prime =6$) dominant cells, representing the major clones from the simulated bulk samples. With this notation, the problem can be formally stated as follows:

Input: A matrix of bulk tumor $B_{m\times n}$ and a matrix of observed SCS samples $C_{m\times k_{observed}}^{\left(observed\right)}$.

Output: A set of inferred fundamental cell components $C_{m\times k\prime }^{\left(inferred\right)}$ and corresponding set of mixture fraction $F_{k\prime \times n}^{\left(inferred\right)}$.

5.5. Assessing solution quality

Given $B_{m\times n}$, $C_{m\times k_{observed}}^{\left(observed\right)}$, and $k\prime$, we would like to find $C_{m\times k\prime }^{\left(inferred\right)}$ and $F_{k\prime \times n}^{\left(inferred\right)}$ such that:

We repeated the above sampling, simulating, and deconvolving procedure for $N=40$ experiments to assess the performance. Figure 2 summarizes the design by which we simulated bulk tumor from the existing SCS data and how we integrated the SCS data into bulk tumor deconvolution.

We used two ways to estimate the performance. We calculated accuracy, measured as the fraction of genomic positions (among the 40 replicates), with a correctly inferred copy number relative to the true cell components. We also measured the accuracy of inferred clones and mixture fractions by the root mean square deviation (RMSD) between true and inferred data.

For each genomic position in each inferred clone, if the copy number is equal to the copy number in the true clone, we consider the value to be accurately inferred. Then, we calculate the copy number error for each genomic position as follows:

We define the accuracy to be the fraction of 0's in $C^{\left(acc\right)}$.

Following Schwartz and Shackney (2010), we also calculate the RMSD of $C_{m\times k\prime }^{\left(inferred\right)}$ and $C_{m\times k\prime }^{\left(true\right)}$, specified as the root mean square distance over all entries of all clones between the two matrices:

Similarly, we can measure the RMSD between $F_{k\prime \times n}^{\left(inferred\right)}$ and $F_{k\prime \times n}^{\left(true\right)}$ over all the loci and mixture fractions:

Following the same idea, when we assess the RMSD between inferred results and true data in each clone, we calculate the RMSD in the pairwise columns of the matrices:

In the inferred trees of 13 nodes each (six inferred, six observed, and a diploid root), we tested for clustering by tumor region as follows. Let the two-element sets $\left\{x_1,x_2\right\},\left\{y_1,y_2\right\},\left\{z_1,z_2\right\}$ be the nodes representing the two observed cells selected from the first region, second region, and third region, respectively. Considering the inferred tree as an undirected tree, we can define unique undirected paths between x₁ and x₂, between y₁ and y₂, and between z₁ and z₂. A pair of nodes/cells representing one region is considered to “cluster together” relative to another region if the path between the two nodes from the region of interest does not pass through either of the two nodes from the other region.

5.6. Fully simulated SCS data

We conducted additional tests on fully simulated data to provide some test case for which the ground truth is known and for which data could be distributed without restriction. We modeled these fully simulated data on the real data to match the true number of regions d ($d=3$ in our case), number of cells $c_i\left(i\in \left\{1,2,3\right\}\right)$ per region, estimated rate $r_{ai}$ of copy number variation a per region $\left(a\in \left\{0,1,2,\dots ,10\right\},i\in \left\{1,2,3\right\}\right)$, and probability $p_{mi}$ that each genomic position that has a nondiploid copy number $\left(m\in \left\{1,2,\dots ,9934\right\},i\in \left\{1,2,3\right\}\right)$.

For each region i, we created a root with diploid copy number in all genomic positions. We then created a binary tree where each node represents one clone. We created copy number vectors for the nodes of the tree by adding mutations according to a Poisson distribution with empirical rate $r_{ai}$ and $p_{mi}$ to mutate the copy number in different genomic positions so that the overall copy number distribution would be similar to that of the real SCS samples. We chose the depth D of the tree to be the minimum value that allows the number of nodes to exceed the number of cells sampled in the real SCS data in each region (in our case, $D=6$).

From the full tree, we constructed subtrees starting from the root with c_i nodes (excluding the root) by random walk with depth first order and with the equal possibility to go right or left (Fig. 11). We collected the chosen nodes to establish an artificial single-cell data set encoded as a matrix whose rows are the genomic positions, and columns are the simulated single-cell samples. We applied the methods described in Section 5.4 on these simulated SCS samples to simulate bulk tumor samples, which we call the fully simulated tumor data.

We then further perturbed the artificial single-cell data to evaluate sensitivity of the methods to noise in copy number alteration (CNA) calling. We randomly choose a desired fraction (0, 0.2, or 0.4) of total genomic positions in each cell sample of $C^{\left(observed\right)}$ to perturb, then sampled a set of uniformly random genomic loci without replacement, and increased or decreased (if the copy number is greater than 1) the copy number by 1. This modification resulted in noisy versions of $C^{\left(observed\right)}$. We then applied our methods and quality assessment as described in Sections 3.1, 3.2, and 5.5 to noiseless and noise-added variants of the fully simulated data.

5.7. Real data analysis

To analyze the effectiveness of the methods on bulk genomic data, we used the three regions of the real bulk data from the patient GBM07 as B in the objective, while we randomly chose 2 cells from single-cell data from each region, used as $C^{\left(observed\right)}$ $\left(k=6\right)$. We fixed the number of components to be six. Since we tried to take advantage of the resolution of single-cell data, we set the regularization parameter to be slightly larger ($\beta =1.0$). This setting would not be expected to hurt the inference since the accuracy of $\beta =1.0$ is only slightly worse than the optimized parameter (Fig. 5). We retrieved a list of known mutated genes that are frequently reported in previous work (Section 3.3) and matched their coordinates to the coordinates of intervals reported in the copy number data by following the matching rules mentioned in Section 2.5. We assert that a gene has a gain if the corresponding interval has copy number at least three, and has a loss if the corresponding interval has a copy number of at most one.

For purposes of visualizing mixture fractions and phylogenies on the fully real data, we randomly selected one result out of the 100 random trials. We present a single trial due to the difficulty of unambiguously matching mixture components and phylogenetic nodes across trials that involved distinct observed single cells and somewhat different inferred clones. We then plotted the fraction of each inferred clone in each tumor region for that one randomly selected trial to obtain Figure 9A and provided the inferred phylogenetic relationship among the inferred clones and observed single cells in that one trial as Figure 9C.

6.1. Inference quality on semisynthetic bulk data via pure NMF

In this subsection, we describe further analysis and additional experimental results on the use of pure NMF or simple heuristic extensions thereof. We assessed the methods with a series of experiments on semisynthetic data designed to assess the phylogeny-free and phylogeny-based methods in comparison with one another and with generic NMF as functions of various parameters of the simulated data and the algorithms. In general, NMF is used to decompose high-dimensional data $\left(m\times n\right)$ to a low rank basis $\left(m\times k\right)$ and fractions $\left(k\times n\right)$, where we usually have $k<<min\left(m,n\right)$. Therefore, we tested the pure NMF algorithm on different numbers of tumor samples to assess the data needs with respect to numbers of tumor samples. For these tests, we explored a wider range of sample numbers (3–99) than in the main article to get a sense of whether performance improvement saturates for unrealistically large numbers of samples. As described in the preceding section, we assume that we have $k\prime =6$ dominant clones (derived from 2 out of random 25 cells selected in each region) as well as 23 rare clones and assess our ability to infer the dominant clones, treating the remainder as noise.

The copy number is predominantly diploid (2) across genomic positions for these data. We therefore trivially achieve relatively high accuracy either in each clone or in all clones by guessing that all clones are purely diploid. As shown in the left plot in Figure 12, the accuracy in some clones can reach as high as 80%. The overall average accuracy resulting from setting all genetic positions to be diploid, indicated by the dashed line, means that over 70% of the entries are diploid. The right plot in Figure 12 shows that the overall RMSD in copy numbers is also small for an all-diploid baseline solution. The variance means that in each clone, there are some loci of which the copy numbers are far away from diploid, as shown in Figures 4 and 6. However, the all-diploid baseline is otherwise uninformative, since we are specifically interested in studying CNAs that change copy number for a subset of the genome. Nevertheless, we can use the accuracy of the all-diploid solution as a baseline against which to assess our algorithms.

We first applied NMF with random initialization to decompose the bulk tumor to infer the clonal cell components and the fraction matrix. We extended the test from small (3) to large (99) numbers of tumor samples with the regularization parameter $\alpha$ changing from 0 to 1 with the increment of 0.2. The $\alpha =0$ results in the top three plots in Figure 13 show that the pure NMF method performs poorly in resolving the problem. For some $\alpha$, we can achieve a good estimation of accuracy and RMSD in copy numbers in smaller numbers of tumor samples (3, 6, 9), but the improvement is generally minimal for larger numbers of tumor samples (33, 99), although some specific $\alpha$ values can return good estimates. These results indicate that adding the penalty to the objective has the potential to lead to a good local optimization. We also note that the performance of mixture fraction inference does not change much either across different numbers of tumor samples or among different $\alpha$.

The bottom three plots in Figure 13 show results of applying the full phylogeny-free method, with real SCS data for initialization and objective penalty, to a larger range of sample sizes. The results of smaller numbers of tumor samples (3, 6, 9) have been shown in the main article, but here we include the results from larger numbers of tumor samples (33, 99). As mentioned in the main article, using real SCS data in initialization and the objective function together can lead to good estimates of copy numbers and mixture fractions of clones. Increasing the number of tumor samples substantially above the intrinsic dimension of the mixture does not appreciably improve the peak accuracy, but does appear to make the method more robust to variation in $\alpha$.

6.2. Inference quality on semisynthetic bulk data via phylogeny-free augmentation with SCS data

Although we showed in the main article that the phylogeny-free method can yield reasonably good average results (Fig. 3), we were further interested in how performance of the method might vary across clones inferred. Thus, we calculated the accuracy and RMSD of inferred and true clones pairwise (Fig. 14) to see if the improvement is consistent in each clone. We also assessed this performance with the varying regularization parameter $\alpha$ from 0 to 1 in increments of 0.2. We only show the results for smaller tumor samples (3, 6, 9) since previous results have shown no significant difference for larger numbers of tumor samples. The top plot of Figure 14 shows that fine-tuned $\alpha$ can improve the accuracy of copy numbers in each clone, while continuously increasing $\alpha$ would result in worse performance in each clone. Similar results can also be observed in RMSDs of copy numbers (center, Fig. 14) and RMSDs of mixture fraction (bottom, Fig. 14). These results indicate that the modified NMF method affects the performance for each individual clone rather than only having the effect on a subset of the clones.

6.3. Inference quality on semisynthetic bulk data via phylogeny-based augmentation with SCS data

This subsection elaborates on the results for the phylogeny-based method of Section 2.3, also specifically in evaluating accuracy and RMSD of clones and mixture fractions in each inferred clone as functions of the regularization parameter $\beta$. Figure 15 shows the assessment of performance in each inferred clone with the varying regularization parameter $\beta$ from 0 to 1 in increments of 0.2. The top and center rows in Figure 15 show that the accuracy and RMSDs in copy numbers were poor in each cell component when $\beta =0$, but were substantially improved by adding the phylogeny into the objective. Each cell component is fairly robust to $\beta$ variation and shows substantially better copy number accuracy and RMSD than observed with pure NMF ($\beta =0$). However, the RMSDs of mixture fractions (bottom row, Fig. 15) do not show appreciable improvements relative to pure NMF. Examination of the variances of the mixture fractions showed that adding the phylogeny increased the variances, which suggests that combining other objective function terms with the phylogeny might lead to better estimates of the fractions for each cell component.

6.4. Comparison of phylogeny-free and phylogeny-based methods

We next compared the performance in each inferred clone of the two methods for three, six, and nine tumor samples. We chose the best-performing regularization parameter, as indicated in Figure 7, for each method in each tumor sample, respectively. When comparing the results in each inferred clone, we found that the phylogeny-based method consistently gave better average accuracy and copy number RMSD for each inferred clone for three tumor samples (top and center row, Fig. 16). Also, the number of outliers from the phylogeny-based method is generally smaller than that from the phylogeny-free method in most of the inferred clones. In contrast, for mixture fraction comparison, the phylogeny-based method yields wider variances. This observation indicates that the performance is less consistent in mixture fraction inference, but, given the higher mean accuracy, can be interpreted as an ability to infer much better mixture fractions in a subset of cases (bottom row, Fig. 16).

6.5. Inference quality on fully synthetic data

In this section, we present the results of both phylogeny-free and phylogeny-based methods on fully synthetic tumor samples (Section 5.6) following procedures similar to those of Section 3.1 and 3.2. We varied numbers of tumor samples (3, 6, 9) for these tests and regularization parameters ($\alpha =0.0,0.2,0.4,0.8$ for the phylogeny-free method and $\beta =0.0,0.2,0.4,0.8$ for the phylogeny-based method). For the phylogeny-free method (top row, Fig. 17), the results are qualitatively similar to those on semisimulated data (Fig. 3) across the range of chosen regularization parameters. The phylogeny-free method shows notable improvement over pure NMF for at least some regularization parameters, with the range of effective regularization parameters expanding with increasing numbers of tumor samples. As with the semisimulated data, the results of the phylogeny-based method (bottom row, Fig. 17) for $\beta >0$ are substantially better than those for the phylogeny-free method. Furthermore, they show much less sensitivity to sample size and to regularization parameter across the range of nonzero values considered. The phylogeny-based results on fully simulated data are somewhat better than those seen on semisimulated data, suggesting that the more complicated error profiles of real SCS data do present some challenge to the methods.

Figure 18 shows sensitivity of the methods to noise in $C^{\left(reference\right)}$ at noise levels 0, 0.2, and 0.4. We conducted these tests for a single regularization value of 0.2 for $\alpha$ or $\beta$, chosen because it yielded relatively insensitive performance to sample size in the noise-free tests. For the phylogeny-free method (top row, Fig. 18), the plot shows that accuracy by each measure generally degrades with increasing noise, with copy number inference generally more sensitive to noise than mixture fraction inference. Robustness to noise by each measure increases noticeably with larger sample sizes. For the phylogeny-based method (bottom row, Fig. 18), the plot shows that accuracy is stable with increasing noise for all cases. RMSD results show some variability parameter-to-parameter, but no evident trend with increasing noise. These results show that the phylogeny-free method is robust to noise for a sufficiently large sample size, while the phylogeny-based method is robust to noise even for small sample sizes.

6.6. Average copy number in real bulk tumor and single-cell data

In this section, we examine a visualization of average copy numbers along the genomic loci in real bulk tumor samples and corresponding single-cell data. We matched the coordinates of our data to the hg19 human genomic reference and split the loci into different chromosomes (Fig. 19). We observe some variance in the copy number profile of bulk data (top, Fig. 19), in part, because the copy number changes are not identical across the three bulk samples. For example, two samples show the loss of chromosome 10, while the other does not. Also, not every bulk sample shows the gain of the full chromosome 7, with some only showing a gain of subset of the loci in chromosome 7. We further note that the bulk sequencing reflects only the average genomic profiles of the clones. In the single-cell data (center, Fig. 19), by contrast, we can clearly see the gain of chromosome 7 and the loss of chromosome 10. The loss of chromosome 10 is less pronounced, which may due to the coarse resolution of genomic breakpoints and uncertainty of measurement of calling copy number in individual cells. Some single cells exhibit the loss of full chromosome 10, while others do not. Similar reasons may account for why both bulk and single-cell data show the gain and loss of partial chromosome 11.

The copy number profile of inferred cell components (bottom, Fig. 19) not only reflects the common events in both bulk and single-cell data (e.g., the large copy number of PDGFR) but also shows the similarity with the available single-cell data. The figure shows the gain of chromosome 7 and the loss of chromosome arm 9p, while the loss of chromosome 10 is less pronounced. Comparing these three average copy number profiles, with the assistance of SCS, the variation of copy number changes in some genomic loci [chromosome 7, 9, 17, 19, X(23), Y(24) in this case] has been minimized, providing more confidence about the copy number aberrations in the genomic loci with which we are interested. We can also derive more confidence about the aberrations in chromosome regions that are less commonly altered (e.g., chromosome 19 and chromosome X, which contains ATRX). Most importantly, we show it is possible to accomplish this using only a very small subset of the single-cell data (six cells in this case) to take the advantage of the resolution of SCS technology.