Our bootstrap approach samples sequencing reads with and without variants at genomic positions (Figures 2A, B). The total number of reads sampled at a given position remains the same as in the original dataset, but reads are sampled with replacement. Since the same read can be sampled multiple times, some reads and even variants will be missing at that position in the bootstrap replicate dataset. Reads are resampled for each position, and a bootstrap replicate data is generated for each tumor sample. A pseudo-multi-tumor dataset is then generated by combining bootstrap replicates of tumor samples, and this pseudo-multi-tumor dataset is analyzed using the desired clone prediction method (e.g., CloneFinder) to infer clones (Figure 2C). Similarly, more pseudo-multi-tumor datasets are generated, and clones are inferred in many bootstrap replicates.

The bootstrap approach can be used to build consensus clone sequences. Predicted clones from all replicates are pooled together and (nearly) identical clone genotypes (for a user-supplied SNV count cutoff) are grouped (Figure 2D). Each clone group is then represented by a consensus clone and the proportion of bootstrap replicates in which a clone appears is the bootstrap support for detecting that consensus clone. To construct a consensus clone sequence, a base reconstructed in a greater proportion of bootstrap clones than the desired threshold is selected for each variant position. The base assignment is marked ambiguous when none of the bases have received the minimum desired bootstrap support.

We lastly describe an application of the bootstrap approach to infer consensus tumor evolution estimates of metastatic migration history using predicted clones. As described above, a pseudo-multi-tumor dataset is first generated by sampling reads, and a set of bootstrap clones is inferred (Figures 2A–C). These bootstrap clones are subsequently analyzed (Figure 2E). For example, we infer bootstrap cell migration history using the bootstrap clones from every bootstrap replicate dataset. Then, each migration history is expressed as a collection of individual migration paths. All the migration paths are pooled to build a consensus migration history, where paths with the highest frequency are connected with each other first, followed by others with lower frequencies until the minimum desired bootstrap support requirement is met (Figure 2F).

We implemented and tested the bootstrap approach in our CloneFinder method which is known to perform well for inferring clones using bulk-sequencing datasets from multi-tumor samples (Miura et al., 2018; Miura et al., 2020). We also took this opportunity to advance CloneFinder by adding a step to preprocess the sequencing datasets by FastClone (Xiao et al., 2020) software that clusters SNVs based on VAFs. This allows CloneFinder to start with a larger collection of tumor genotypes than the original approach. Thus, the advanced CloneFinder, CloneFinder+, begins with the analysis of VAFs using FastClone (Figure 3A). For a tumor sample, FastClone clusters SNVs based on their VAF similarities and predicts relationships of SNV clusters, i.e., ancestor-descendant, sibling, or monoclonal (Figure 3B). CloneFinder+ constructs candidate clone sequences by accumulating all predicted SNVs from the root cluster to a target cluster. These clones are the candidate clones for a given tumor sample. This analysis is done for every tumor sample individually, and candidate clones are inferred for each tumor sample. All candidate clones from all the tumor samples are then pooled and duplicate clones are removed. Also, potentially spurious candidate clone sequences are filtered, e.g., those with many ambiguous base assignments (Supplementary Note for details). Lastly, a candidate clone sequence matrix, M, is constructed (Figure 3C). This is a binary matrix, where M _ij = 0/1 represents the absence/presence of a SNV at the jth variant in the ith candidate clone.

Next, we apply CloneFinder to infer clone genotypes (Miura et al., 2018) (Figures 3C–G). A phylogeny of candidate clones is first built using the matrix M. This phylogeny is used to identify missing ancestral clones that have persisted without being replaced by their descendant clones. Ancestral clone sequences are inferred at internal nodes of the phylogeny, and all unique ancestral clones are added to the collection of candidate clones, M (Figure 3E). Then, the presence of these candidate clones in the bulk-sequencing dataset is assessed using VAFs. A matrix of VAFs, V, is constructed, where the number of rows is equal to the number of tumor samples and the number of columns is equal to the number of SNVs (Figure 3D). Next, tumor frequencies of all candidate clones are estimated through regression analysis. Here, V, M, and f, a two-dimensional matrix of estimated clone frequencies of the tumor samples, have the following relationship, 1 / 2 f × M = V . (1)

This relationship is valid only when variants are not affected by copy number alterations (CNAs). Thus, variants that are affected by CNAs need to be excluded from the dataset; or VAFs should be adjusted before the analysis using estimated cancer cell fraction (CCF), i.e., VAF is CCF divided by two. When an estimated clone frequency within a given tumor sample is lower than the desired threshold, that clone is assumed to be absent from the tumor sample. Candidate clones predicted to be absent from all tumor samples are removed.

Next, CloneFinder + decomposes potential hybrid sequences using the strategy in CloneFinder (Miura et al., 2018) (Figure 3F). Briefly, SNVs with the same expected VAF are subclassified into two groups based on their similarity to observed VAFs. Alternate candidate clone sequences are constructed by combining SNV groups. All possible combinations of SNV groups are added to the M and are evaluated using the regression analysis (Eq. 1). All combinations estimated to be present in appreciable frequency are retained. Accordingly, M is updated, and this process of searching ancestral and hybrid clones is repeated until no new clone sequences are added. Lastly, CloneFinder+ finalizes the predicted clone sequences, where variants not assigned to any clone sequences are added to a clone based on their VAF similarities (Figure 3G).

We obtained previously simulated datasets from https://github.com/raphael-group/machina. These datasets were generated by modeling the evolution of primary and metastatic tumors (El-Kebir et al., 2018). In this simulation, metastatic tumors were founded by one or more than one cancer cells that migrated from another tumor site (primary or another metastatic tumor). Clones were defined as a group of cells with the same sequences. Thus, a new clone could arise even from a single mutation. Each clone phylogeny was unique, and the number of clones was 6–26. Each dataset contained 9–99 SNVs. We excluded five datasets that contained tumor samples with only one variant because this is unrealistic in empirical data. In total, we analyzed 75 datasets.

We obtained metastatic cancer datasets from Zhao et al. (2016), who performed bulk sequencing (exome sequencing) and identified SNVs. For each SNV, these datasets contained total sequencing read counts and counts of reads with altered bases (SNVs). In total, we obtained 40 datasets with three to eight tumors.

For each SNV, we predicted if it was a driver mutation using the Cancer Genome Interpreter (CGI) web tool (http://www.cancergenomeinterpreter.org), which uses the OncodriveMUT method (Tamborero et al., 2018). We also used the CRAVAT web tool (http://www.cravat.us), which performs the CHASM prediction (Carter et al., 2009; Douville et al., 2013). Driver mutations were predicted without specifying a cancer type.

To map mutations at branches of the phylogeny, we analyzed predicted clone sequences and reconstructed ancestral clone sequences using MEGA (Kumar et al., 2012; Tamura et al., 2021).

Each simulated and empirical dataset was analyzed using CloneFinder+, and the reliability of the inferences was assessed using the bootstrap approach. In the bootstrap analysis, we generated 30 replicates because using more (100) replicates essentially produced the same result (Supplementary Figure S1). In the CloneFinder + analysis, we clustered variants without giving the tumor purity (the step of FastClone analysis) and used variants with at least 50 reference read counts and two mutant read counts to assess the quality of candidate clones. Note that CloneFinder+ does not require the value of tumor purity and the maximum number of clones to be inferred. During the analysis, candidate clone genotypes with <1% clone frequencies for all tumor samples were discarded. We grouped identical bootstrap clones to derive consensus clone sequences while allowing at most one base assignment difference. We selected a base in >90% of bootstrap clones for each variant position.

We inferred cell migration history using the PathFinder method (Kumar et al., 2020), and the reliability was assessed using the bootstrap approach. PathFinder was performed by providing the correct primary tumor sites, sequences of CloneFinder+ clones that were predicted with >5% clone frequencies, and tumor sites that contained each clone. Note that PathFinder does not require the value of tumor purity.

To evaluate the accuracy in inferring correct clones, we paired each simulated clone sequence with its most similar inferred clone sequence. We allowed an inferred clone to be paired with more than one simulated clone. We counted the number of sequence differences between inferred and simulated clones paired. We calculated the average when more than one inferred clone was paired with a given simulated clone, which was divided by the sequence length to estimate genotype error (GE) for a given simulated clone.

To evaluate the accuracy of inferred migration history, we counted the number of migration paths that were correctly inferred, those not identified, and incorrect paths following a previous study (Kumar et al., 2020).

The bootstrap approach is widely employed in molecular evolutionary and phylogenetic analyses to estimate variances and confidence limits (Efron and Tibshirani, 1994; Nei and Kumar, 2000). Our bootstrap approach samples sequencing reads. So, VAFs will be perturbed, with greater perturbations experienced by VAFs computed from a small number of reads (Figure 1). Since the number of reads is often highly variable among SNVs in empirical data [e.g., 8–1,436 for ATP401 from Zhao et al. (2016)], the variances of their VAFs are also variable. We first show the GEs of inferred clones and their bootstrap supports.

We found that the median GE of the clones with higher bootstrap supports (>50%) was zero (Figure 4A). Therefore, clones with >50% bootstrap support can be considered to be reliably inferred. We also found that many bootstrap consensus clones (77%) received low bootstrap supports (≤10%), and these clones contained many genotype errors (GE >5%; Figure 4A). Interestingly, these clones had a wide range of GEs, even overlapping with those from higher bootstrap supports. In an extreme case, a few consensus clones with <10% bootstrap support had correct clone sequences, i.e., GE = 0. Therefore, some simulated (correct) clones were not repeatedly detected among bootstrap replicates, reducing bootstrap support.

We next tested if different scenarios of the cell migration history affected the performance of the bootstrap approach. These simulated datasets were generated by modeling the evolution of primary and metastatic tumors, where metastatic tumors were seeded by clones that migrated from either a primary or a metastatic tumor. We classified the datasets into (1) those with metastatic tumors that received a single seeding clone from another tumor site, i.e., monoclonal seeding, and (2) those with more than one seeding event, i.e., polyclonal seeding. Thus, intra-tumor heterogeneity of tumors from polyclonal seeding is higher, while inter-tumor heterogeneity among tumor sites is smaller than those with monoclonal seeding.

We found that the bootstrap approach performed well on both types of datasets, as inferred clones with greater bootstrap confidence tended to have more accurately inferred clone sequences (lower GE) (Figures 4B, C). We also found that GE tended to be slightly better for datasets with monoclonal seeding, indicating that inferred clone sequences were slightly more accurate for data with monoclonal seeding, i.e., lower intra-tumor heterogeneity with higher-inter tumor heterogeneity. This pattern was consistent with previous studies (Miura et al., 2020).

We also tested the impact of the number of tumors on the performance of the bootstrap approach. These simulated datasets were composed of five or eight tumors, i.e., four or seven metastatic tumors with a primary tumor per dataset. For both datasets with five and eight tumors, we similarly found that clones with greater bootstrap confidence tended to have more accurate clone sequences, i.e., a lower GE (Figures 4D, E). Therefore, the number of tumors in a dataset did not affect the performance of the bootstrap approach to place the reliability, while inferred clone sequences were more accurate for datasets with a larger number of tumors (lower GE). This was consistent with previous studies (Miura et al., 2020). Overall, these results suggested that the bootstrap approach is useful for assessing the reliability of inferred clone sequences.

Next, we tested the performance of the bootstrap approach using 40 empirical datasets from various cancer types. As observed for the simulation study, most of the bootstrap consensus clones had low bootstrap supports (<10%), and few clones were identified with good (>50%) bootstrap supports (Figure 5A). We did not find a clear association between cancer types and bootstrap support values, as the distribution of bootstrap values was similar among different cancer types (Figure 5B). We also did not observe an association between tumor mutation burden (number of mutations) and bootstrap values for any cancer types (Figures 5C–F). Therefore, the performance of the bootstrap approach was not affected by either cancer type or tumor mutation burden.

Since the bootstrap approach performed well to place a confidence limit on inferred clones, we next tested if the bootstrap approach is also useful to assess the reliability of a downstream inference of predicted clones. As an example of a downstream analysis of predicted clones, we inferred metastatic cell migration histories using the same simulated datasets.

We found that correct paths often had a high bootstrap support (a median bootstrap support = 86.5%), while incorrect paths tended to have low bootstrap supports (a median of 10%) (Figure 6). However, as observed in the analysis of inferred clone sequences, a few correct paths were not well supported. Actually, bootstrap support for correct paths varied considerably, indicating that these paths were not repeatedly found in many bootstrap replicate datasets. Thus, some migration paths are difficult to reconstruct, which is consistent with previous findings (Kumar et al., 2020). Overall, these results indicated that the reliability of inferred migration path from predicted clones could be assessed using the bootstrap approach.

Since point estimates of clones and cell migration paths have limited accuracy, placing bootstrap confidence on inferences is essential in real empirical data analysis. As an example, we illustrate empirical data analysis for revealing the patterns of metastatic cell migrations and driver mutation occurrences. In this analysis, we excluded eight datasets without primary tumors because the information on the primary tumor site was necessary for the inference of migration history.

First, we show the result of a lung cancer patient with three metastatic tumors in the heart, liver, and GEJ (ATP401 patient). The inferred clone phylogeny using CloneFinder+ (without bootstrap reliability assessment) indicated that clone C1 originated from the root clone (the most recent common ancestor of all clones), and a metastatic tumor in the liver contained this clone (Figure 7A). The PathFinder analysis (without bootstrap reliability assessment) using these predicted clones and phylogeny produced by CloneFinder+ predicted that clone C1 migrated from the primary tumor site (lung) to the liver. While 128 mutations were mapped on the branch leading to the C1 (B1 branch in Figure 7A), the phylogeny does not tell us whether C1 acquired all, any, or a subset of new mutations in the lung. Nevertheless, seven B1 mutations were predicted to be drivers according to CGI (Tamborero et al., 2018), which may be important for metastasis. On the other hand, we also found that some migration events were not associated with any detected mutations as no mutations were mapped at corresponding branches (Figures 7A, B). For example, clone C2.1 from metastasis GEJ seeded a heart metastasis without associated mutations, as the same clone is found in both locations. Overall, we found a greater average number of driver mutations in clones that moved from the primary tumor than those from metastatic tumors for this patient (Figure 7B).

We next demonstrate how this observed pattern could be validated using the bootstrap approach. We found that all of the predicted cell migration paths in the single point PathFinder inference (Figure 7B) were supported with high bootstrap confidences (>60%), validating the PathFinder inference (Figure 7C). In the bootstrap analysis, we mapped mutations on each bootstrap migration path and calculated the average driver mutation count per path for each replicate. We found that the number of driver mutations per path was significantly greater for those from the primary tumor than those from a metastatic tumor (p < 0.01; t-test) (Figure 7D). Since migration paths with a larger number of associated mutations may simply result in a larger number of driver mutations, we normalized the count of driver mutation, i.e., we computed the driver mutation rate by dividing it by the total number of mutations for a path. We excluded paths without any associated mutations. Similarly, the driver mutation rate was significantly greater for those from the primary tumor than those from a metastatic tumor (p < 0.01 by t-test) (Figure 7E). Therefore, driver mutations occurred more frequently at migrations from the primary tumor than those from metastatic tumors for this patient.

To test if most of the patients similarly had higher driver mutation rates for migration paths from the primary tumor than metastatic tumors, we analyzed 32 datasets of metastatic cancer patients. Similar to the ATP401 patient, many paths were from metastatic tumors (Figure 8A), indicating that migration events from metastatic tumors were not rare, consistent with previous studies (Kumar et al., 2020; Chroni et al., 2022). We found that only six patients showed significant differences in driver mutation rates between those from primary tumors and from metastatic tumors (p < 0.01 by t-test for both driver prediction methods), and all of them had higher driver mutation rates for paths from the primary tumors, which was similar to the ATP401 patient (Figure 8B). However, most of the patients (26) did not show a significant difference in driver mutation rates between migration paths from primary and metastatic tumors (p > 0.01 by t-test for at least one driver prediction method). Therefore, the numbers of driver mutations were often not significantly different between metastatic cell migration events sourced from primary tumors and those from metastatic tumors.