MEDICC2: whole-genome doubling aware copy-number phylogenies for cancer evolution

Abstract

Aneuploidy, chromosomal instability, somatic copy-number alterations, and whole-genome doubling (WGD) play key roles in cancer evolution and provide information for the complex task of phylogenetic inference. We present MEDICC2, a method for inferring evolutionary trees and WGD using haplotype-specific somatic copy-number alterations from single-cell or bulk data. MEDICC2 eschews simplifications such as the infinite sites assumption, allowing multiple mutations and parallel evolution, and does not treat adjacent loci as independent, allowing overlapping copy-number events. Using simulations and multiple data types from 2780 tumors, we use MEDICC2 to demonstrate accurate inference of phylogenies, clonal and subclonal WGD, and ancestral copy-number states.

Keywords: Aneuploidy; Cancer evolution; Chromosomal instability; Intratumor heterogeneity; Phylogenetic reconstruction; Single-cell sequencing; Somatic copy-number alterations; Whole-genome doubling.

Fig. 1

MEDICC2 algorithm. a MEDICC2 infers cancer phylogenies from SCNA data from single cells or bulk sequencing using a minimum-event distance (MED) and infers the ancestral genomes. It allows for backmutations, obeys biological constraints, and solves the phylogeny problem where ancestral genomes are not sampled. b Computing distances with WGD. Copy-number profiles are represented as vectors of positive integer copy numbers across chromosomes (here: two chromosomes with four segments each). To infer the correct MED, LOH events are considered first as lost segments cannot be re-gained by later events. WGD events span the full copy-number profile, whereas gain and loss events can affect an arbitrary number of segments within a chromosome. c Symmetric distance calculation. The MED from an ancestral state to a sample is asymmetric due to biological constraints. The final symmetric distance between two samples is computed as the sum of distances from an ancestral genome to both samples, while minimizing over all possible ancestors. d Schematic overview of the MEDICC2 workflow. Haplotype-specific copy-number profiles are either pre-phased or undergo evolutionary phasing (see e). Pairwise MEDs are computed between all genomes, followed by tree inference and ancestral reconstruction which determines the final branch lengths of the tree. Results are reported to the user as a patient summary and plot. e Evolutionary phasing. Copy-number profiles for both alleles are jointly encoded as an unweighted phasing FST P where both possible allele configurations are encoded at each position in the sequence. Evolutionary phasing then determines the optimal configuration (bold arrows) and extracts final haplotypes (orange and blue) by computing the MED between the phasing FST and two reference haplotypes. An example of major/minor copy number, phased copy number, and the MED from the diploid is shown at the bottom. Abbreviations: FST: Finite-state transducer, MED: Minimum-event distance, LOH: Loss-of-heterozygosity, WGD: Whole-genome doubling

Fig. 2

Algorithm performance and validation. a Runtime of different composition strategies. Copy-number profiles were simulated with increasing lengths from 20 to 200 segments. Computation time of the MEDs is linear with respect to the length of the input sequences. While MED-WGD took significantly longer to compute than the MED without WGD, the new lazy composition strategy reduced runtime by orders of magnitude. Shaded areas correspond to standard errors. b Performance on simulated data: Using an independent simulation routine we benchmarked MEDICC2 against a range of other methods. The reconstructed trees were compared to the simulated trees using the generalized Robinson-Foulds distance. As expected, the GRF distance rises with increasing tree size. MEDICC2 outperforms all other methods for all tree sizes. c Validation of MEDICC2 events with SVs. Pairs of MEDICC2 events and SVs were chosen based on an overlap of the starting segment. We assume MEDICC2 events to be supported by the SV if the ends also overlap. Shown here are the results using only duplications and deletions with size larger than 10Mbp. d The MEDICC2 WGD score for 2778 cancer genomes. Individual cancers are plotted based on their average ploidy and fraction of genome with LOH. The original separating line between WGD and non-WGD tumors was estimated by Dentro et al. as y = 2.9 − 2x. Correct “WGD” and “no WGD” predictions from MEDICC2 were marked in orange and blue while false predictions were marked in black and gray (latter if the PCAWG WGD status was “uncertain”). Abbreviations: NJ: Neighbor joining, Min. Ev.: Minimum Evolution

Fig. 3

Evolutionary history of tumor subclones from patient A31. a SNV-based phylogeny. Reproduction of the SNV-based phylogeny as described in Gundem et al. [40] for the multi-sample prostate cancer tumor case with one sample (C) from the primary tumor and four samples (A, D, E, and F) from distinct metastatic sites. Original reconstruction was performed using an n-dimensional Bayesian Dirichlet process to cluster estimated cancer cell fractions of the single-nucleotide variants (SNV) identified in the WGS across samples. Only the major subclone of each sample is shown (“Methods”). b MEDICC2 phylogeny. Using multi-sample phased copy-number profiles, MEDICC2 detected the presence of WGD in the metastatic samples and its absence in the primary sample from A31. The MEDICC2 analysis identifies multiple MSAI events as well as parallel LOH on 6 and 13 (purple arrows). Individual events are marked in the copy-number track where they occur: gains (red) and losses (blue). The gray number in each branch corresponds to its bootstrap-confidence score while the WGD events from the MEDICC2 event detection are marked in green

Fig. 4

Event detection for the Gundem et al. [40] cohort. a WGD detection. In the 10 patients, a total of 4 WGDs were detected, two of which were clonal, one subclonal and one in a terminal branch. b Distribution of arm-level events. Using the MEDICC2 event detection routine, we detected the number of times a whole chromosome arm was either gained or lost in a single branch. The gains and losses were aggregated over all patients and samples into a single score. This score was compared against the oncogene - tumor suppressor gene (OG-TSG) score derived by Davoli et al. [41]. A clear correlation between the gains/losses and the OG-TSG score (which is not based on copy numbers) is visible. c Distribution of gene-level events. The analysis was repeated on the basis of all 1729 individual genes present in the Davoli et al. dataset. On the x-axis, we plotted the base-10 logarithm of the genes’ p-values and flipped the sign for the oncogenes to create a single, continuous x-axis for both genesets. A small correlation is visible which becomes more pronounced when only considering the top 100 genes. Names are given for genes with p < 10⁻²⁰

Fig. 5

Inferred phylogeny for single-cell data with 1023 cells. Inferred phylogeny and allele-specific copy-number profiles for patient TN2 from Minussi et al. [46]. The diploid and most recent common ancestor to all cells are marked with green and blue circles, respectively. We manually selected clades from the phylogeny to match the superclones and subclones of the original publication. These are marked next to the tree in the colors of the original publication and with horizontal lines. The structure of the tree corresponds very clearly with distinct features of the copy-number profiles and matches the clonal structure derived in the original publication. Selected synapomorphies of the clone structure are highlighted with a yellow border and annotated on the figure