Tracking clonal evolution during treatment in ovarian cancer using cell-free DNA

Abstract

Emergence of drug resistance is the main cause of therapeutic failure in patients with high-grade serous ovarian cancer (HGSOC)¹. To study drug resistance in patients, we developed CloneSeq-SV, which combines single-cell whole-genome sequencing² with targeted deep sequencing of clone-specific genomic structural variants in time-series cell-free DNA. CloneSeq-SV exploits tumour clone-specific structural variants as highly sensitive endogenous cell-free DNA markers, enabling the relative abundance measurements and evolutionary analysis of co-existing clonal populations over the therapeutic time course. Here, using this approach, we studied 18 patients with HGSOC over a multi-year period from diagnosis to recurrence and showed that drug resistance typically arose from selective expansion of a single or small subset of clones present at diagnosis. Drug-resistant clones frequently showed interpretable and distinctive genomic features, including chromothripsis, whole-genome doubling, and high-level amplifications of oncogenes such as CCNE1, RAB25, MYC and NOTCH3. Phenotypic analysis of matched single-cell RNA sequencing data³ indicated pre-existing and clone-specific transcriptional states such as upregulation of epithelial-to-mesenchymal transition and VEGF pathways, linked to drug resistance. In one notable case, clone-specific ERBB2 amplification affected the efficacy of a secondary targeted therapy with a positive patient outcome. Together, our findings indicate that drug-resistant states in HGSOC pre-exist at diagnosis, leading to positive selection and reduced clonal complexity at relapse. We suggest these findings motivate investigation of evolution-informed adaptive treatment regimens to ablate drug resistance in future HGSOC studies.

Fig. 1. Clone-specific mutations and structural variations in scWGS.

a, scWGS-based copy-number heatmap for patient OV-004. Each row is the copy number of a cell, cells are ordered according to a MEDICC2-computed single-cell phylogeny (shown on the left). b, Clone pseudobulk copy number at 10-kb resolution on chr. 17 for clones A and B. Truncal variants (TP53 missense and deletion) are annotated in purple, clone-specific duplications and SNVs are annotated in red and yellow, respectively. c, Phylogenetic trees annotated with cells that have support for variants shown in b. d, Clone pseudobulk copy number at 10-kb resolution on chr. 8 for clones A and B showing divergent chromothripsis on chr. 8. In b and d, notable regions that are different between clones A and B are highlighted in grey, and the absolute difference (Δ) between copy-number states is shown in the top track.

Fig. 2. SVs as highly specific markers of tumour DNA in cfDNA.

a, Schematic of CloneSeq-SV workflow illustrated with a translocation between chr. 8 and chr. 19 identified in OV-107. Read pileups were generated from BAM files using svviz. b, Distribution of VAFs for truncal SVs and SNVs in baseline samples (n = 21). c, Tumour fraction computed from truncal SVs versus tumour fraction computed from TP53 VAF. Plot annotated with the Pearson correlation coefficient and associated two-sided P value (P < 10⁻¹⁰). d, Schematic showing how patient-specific error rates are calculated by applying probe sets to off-target patients. e, Average background error rates in duplex, simplex and uncorrected sequences. Each violin or boxplot is a distribution over SVs and SNVs where each data point is the error rate for an individual patient (n = 21). Stars indicate P values from a two-sided t-test: ***P < 0.001. P values per group: duplex, 1.5 × 10⁻⁵; simplex, P < 10⁻¹⁰; uncorrected, P < 10⁻¹⁰. f, Fraction of SNV and SVs that have 0 background, that is, no read support in an incorrect patient. Box plots indicate the median, first and third quartiles (hinges) and the most extreme data points no farther than 1.5× interquartile range (IQR) from the hinge (whiskers).

Fig. 3. Detecting clone-specific SVs in cfDNA.

a–d, Single-cell phylogeny with tips coloured by clone membership (left), copy-number (CN) profiles of chromosomes of interest with clone-specific SVs driven by a mutational process: chromothripsis in OV-083, breakage-fusion bridge in OV-045, tandem duplication towers in OV-081 and a complex intra-chromosomal event in OV-022. The location of SVs are indicated above the copy-number profiles (middle), CCF of the two clones of interest in DLP, the number of clone-specific SVs and the VAF of those clone-specific SVs in cfDNA at baseline (right). Patient numbers: 083 (a), 045 (b), 081 (c) and 002 (d). e, VAF of all SVs at baseline in cfDNA stratified by clonality. Stars indicate P values from a two-sided t-test: ***P < 0.001, **P < 0.01, *P < 0.05, NS > 0.05. Box plots indicate the median, first and third quartiles (hinges) and the most extreme data points no farther than 1.5× IQR from the hinge (whiskers).

Fig. 4. Clonal evolution of drug resistance in patients.

Clonal evolution tracking in four patients. a, Anatomical sites sequenced with DLP, a phylogenetic tree of the clones, then clonal fractions, mean truncal SV VAF and TP53 VAF in cfDNA, CA-125 and treatment history over time for patient 044. Disease recurrences are annotated by grey boxes on the CA-125 track. b, Distribution of ERBB2 copy number across cells in clone B (n = 205) and E (n = 146) showing a two-sided t-test, P < 10⁻¹⁰. c, Pseudobulk copy number of clones B and E at 10 kb resolution in chr. 2 and chr. 17. A translocation specific to clone E and implicated in the ERBB2 amplification is highlighted. Bottom panel shows the read counts of this translocation across time points in cfDNA. d, CT scan images from day 0 and day 84 at indicated sites. e,g,i, Orange and white arrows indicate site of disease clonal tracking in patients 009 (e), 107 (g) and 045 (i). f, Diagram of mutations affecting the BRCA1 gene: location of frameshift deletion (FS del.) shown with a red dashed line; large 1.37 kb deletion shown in grey. h, NOTCH3 and CCNE1 single-cell copy-number distribution across clones (P < 10⁻¹⁰ for both comparisons, numbers of cells A = 294, B = 61, C = 24, D = 39) and CCNE1 DNA FISH image from patient 107. j, RAB25 and CCNE1 single-cell copy-number distribution across clones (P < 10⁻¹⁰ for both comparisons, numbers of cells A = 931, B = 308, C = 483, D = 59) and CCNE1 DNA FISH image from patient 045. Stars indicate P values from a two-sided t-test: ***P < 0.001. Box plots indicate the median, first and third quartiles (hinges) and the most extreme data points no farther than 1.5× IQR from the hinge (whiskers). In the FISH image keys (h,j), Cent. 19 refers to the centromere of chromosome 19. rev., reversion; del., deletion; tra., translocation; CN, copy number; R1, recurrence 1; R2, recurrence 2. Illustrations in a,e,g,i were created using BioRender (https://biorender.com).

Fig. 5. Clone-specific transcriptional programmes.

a, Hallmark pathway variability across genomically defined clones in scRNA-seq data from 20 patients. Each data point represents the maximal pathway score difference between clones in each patient. b, Clone frequencies inferred from cfDNA at baseline (B) and recurrence (R) for OV-107. c, UMAPs labelled by clone mapping (inferred using TreeAlign) and sites. d, Distribution of NOTCH3 expression, VEGF pathway, hypoxia and HIF1A across clones. Plots are annotated with the result of a two-sided t-test comparing distributions between the dominant clone at baseline (n = 2,735 cells) versus the dominant clone at recurrence (n = 132 cells), ***P < 0.001, P < 10⁻¹⁰ for all comparisons except for hypoxia (P = 3.2 × 10⁻⁶). e, Clone frequencies inferred from cfDNA at baseline (B) and recurrence (R) for OV-009. f, UMAPs labelled by clone mapping (inferred using TreeAlign) and sites. g, Distribution of EMT pathway, VIM expression, JAK-STAT pathway and fraction of cells in each cell cycle phase. Plots are annotated with the result of a two-sided t-test comparing distributions between dominant clone at baseline (n = 823 cells) versus dominant clone at recurrence (n = 3,738 cells), ***P < 0.001, P < 10⁻¹⁰ for all comparisons. Box plots indicate the median, first and third quartiles (hinges) and the most extreme data points no farther than 1.5× IQR from the hinge (whiskers). UPR, unfolded protein response; DN, down; V1, version 1; V2, version 2; UV, ultraviolet; RUQ, right upper quadrant; LUQ, left upper quadrant.

Fig. 6. Evolutionary modelling.

a, Summary of the approach used to accept or reject neutrality. Frequency of clones at baseline and changes in cancer cell population informed by CA-125 concentrations are used as input to a neutral Wright–Fisher (WF) model with varying population sizes. For each sample, 1,000 simulations temporally matched to clinical time course and cfDNA measurements are generated and then the distribution of frequencies at the final time point are compared with observed values. b, Example simulated trajectories and observed frequencies for three patients: 006, 009 and 045. Patients 009 and 045 have clones that deviate from the expectations in a neutral model, whereas clones in patient 006 are consistent with a neutral model. c, Summary of the results of the Wright–Fisher simulation-based test in 16 patients. From bottom to top, a change in clone frequencies between baseline and the final time point that had evidence of ctDNA (in most cases the final time-point samples), whether clones were observed in any time point postsurgery, adjusted P values (Holm–Bonferroni correction) from Wright–Fisher neutrality test (Methods) per clone, neutral or non-neutral classification based on a cut-off of P(adjusted) < 0.05.

Extended Data Fig. 1. Study cohort.

Swimmer plot showing clinical history of all 24 patients included in the study. Shown are survival status, therapies, surgeries time of first clinical recurrence and data generation timepoints. Days are relative to day of first surgery, ie Day 0 is the date of primary debulking or laparoscopic biopsy.

Extended Data Fig. 2. Data metrics.

a) Number of clonal and subclonal SVs per patient b) Total number of SVs called per patient by SV type c) Distribution of coverage per cell per patient d) Pseudobulk coverage per cell (summed coverage across all cells) e) Number of high quality cells per patient.

Extended Data Fig. 3. Study summary.

a) Study summary, showing clinical history of a typical HGSOC patient and specimen sample collection protocol. b) Workflow showing clonal evolution tracking using structural variants identified in single-cell whole genome sequencing and assigned to clones using single-cell phylogenetics. These clone specific SVs are then followed in cfDNA using deep duplex error corrected sequencing. Illustrations were created using BioRender (https://biorender.com).

Extended Data Fig. 4. High resolution copy number profiles.

Copy number plots of chromosome 8 and 19 from OV-004 using 500 kb bins a) and 10 kb bins b). c) proportion of SVs that could be matched to copy number transitions using 10 kb and 500 kb bins.

Extended Data Fig. 5. Theoretical limit of detection in cfDNA.

a) Detection probability as a function of cfDNA tumour fraction, number of mutations assayed and coverage using formula provided by Zviran et al. b) Theoretical limit of detection (probability of detection > 0.99) as a function of the number of mutations (x-axis) and Coverage (colors). Error rates of SNVs and SVs estimated in duplex and uncorrected sequencing from our data are shown on the left and right respectively. Star indicates relevant parameters for this study (100 mutations and 1000X coverage).

Extended Data Fig. 6. Single cell copy number heatmaps and phylogenetics.

scWGS copy number heatmaps and phylogenetic trees for the 18 patients with longitudinal tracking data. Patient ID and the total number of cells are indicated at the top of each plot. Each row shows the copy number profile of a cells, rows are ordered by the MEDICC2 derived phylogenetic tree shown on the left of each plot. Trees are coloured by clone assignments. Only cells assigned to clones are shown.

Extended Data Fig. 7. Clonal evolution in patients.

a–k) Clonal evolution tracking in 11 patients. For each patient we show the anatomical sites sequenced with DLP, a phylogenetic tree of the clones, clonal fractions, mean truncal SV VAF and TP53 VAF from cfDNA, CA-125 and treatment history over time. l) Summary of the clonal composition at baseline and recurrence (final time point if more than one post-recurrence time point) for 18 patients. m) Distribution of shannon entropy at baseline and recurrence n) Number of clones detected at baseline and recurrence. Plots are annotated with p-values from a paired two-sided t-test. Box plots indicate the median, first and third quartiles (hinges) and the most extreme data points no farther than 1.5× IQR from the hinge (whiskers). IQR, interquartile range. Illustrations in a–k were created using BioRender (https://biorender.com).

Extended Data Fig. 8. Tumour suppressor mutations.

a) Oncoplot showing presence/absence of tumour suppressor mutations across patients and clones. Filled squares indicate instances where at least one sequencing read in a clone shows evidence for the mutation. Absence of reads supporting a mutation is inconclusive in some clones when the number of cells is low, these are annotated with black squares. Patients are ordered by mutational signature (HRD-Dup, HRD-Del, FBI, TD). b) VAF of somatic tumour suppressor mutations across time, shown are non-TP53 mutations from patients with time series.

Extended Data Fig. 9. Whole genome sequencing of cfDNA.

a) Top: clone frequencies at baseline and recurrence. Bottom: summed pearson correlation coefficient of normalised read counts and B-allele frequency values between cfDNA WGS and clone level consensus copy number profiles per clone. b) Normalized read counts at baseline and recurrence from whole-genome sequencing of cfDNA from 6 patients. Black dots are the data, red dots are predictions based on copy number profiles from DLP and inferred tumour and clone fractions from CloneSeq-SV. Patient ID and the tumour fraction (TF) based on TP53 mutation are indicated above each plot. cfDNA WGS data from baseline (left) and recurrence (right) are shown. c,d) Clone frequencies over time calculated from SVs (c) and SNVs (d) for patient OV-045. e) Scatter plot of all clone frequencies calculated using SNVs and SVs, dashed line indicates y = x line. Included in this plot are clone frequency estimates from samples with purity > 0.5% and clones with at least 5 SVs and SNVs. R and p show the Pearson correlation coefficient and associated p-value.

Extended Data Fig. 10. Single cell whole genome sequencing of second time points.

a) Clonal tracking in patient 046 b) patient 026 and c) patient 139. For each patient we show the anatomical sites sequenced with DLP, a phylogenetic tree of the clones, clonal fractions, mean truncal SV VAF and TP53 VAF from cfDNA, CA-125 and treatment history over time. For these patients, scWGS was generated from a second timepoint, anatomical site and time of second time point are annotated on the anotamical diagram and at the top of the clone fraction plot (S1 vs S2). Comparison of pre-and post treatment scWGS data for these 3 patients: d) 026 f) 139 and g) 046. For each panel we show the phylogenetic tree of all cells annotated by clone, timepoint (pre or post treatment) and whole genome doubled (WGD) state. Then the consensus total copy number heatmap for each clone based on pre-treatment cells and the total copy number heatmap for all post-treatment cells. e) Shows a pre-treatment WGD cell from clone C in 026 and the consensus copy number profile of all post-treatment cells. Illustrations in a-c were created using BioRender (https://biorender.com).