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

Sample collection

All enrolled patients were consented to an institutional biospecimen banking protocol and MSK-IMPACT testing⁴², and all analyses were performed per a biospecimen research protocol. Eligible patients had International Federation of Gynecology and Obstetrics stage III to IV HGSOC with diagnosis verified through clinicopathological review. All protocols were approved by the Institutional Review Board of the Memorial Sloan Kettering (MSK) Cancer Center. Patients were consented following the Institutional Review Board-approved standard operating procedures for informed consent. Written informed consent was obtained from all patients before conducting any study-related procedures. The study was conducted in accordance with the Declaration of Helsinki and the Good Clinical Practice guidelines. We collected fresh tumour tissues at the time of upfront diagnostic laparoscopy or debulking surgery, as described in ref. ¹⁸. Blood collection was carried out longitudinally over a 5-year period (2019–2024). Two Streck tubes for cfDNA were collected in each visit. If possible, blood was collected in the Outpatient Clinic at the MSK Cancer Center. Alternatively, blood samples were collected in the operating room when patients were undergoing debulking surgery or laparoscopy.

Sample processing

Streck tubes were submitted to the MSK laboratory medicine facility after collection and processed for plasma and buffy coat separation, as well as DNA extraction¹³ (MagMAX cfDNA isolation kit).

Clinical data

In this cohort study, we extracted clinical annotations from electronic health records of 24 patients treated at the MSK Cancer Center for HGSOC. For these patients, we collected contemporaneous longitudinal data from their initial HGSOC diagnosis, as well as historical data, if available. Clinical data included laboratory measurements, surgical procedures and medications. CA-125 measurements were obtained as part of patients’ routine clinical care from blood samples collected at baseline, during therapy and subsequent follow-up visits. All dates are relative to the time of first surgery for each patient: that is, day 0 is the date of primary debulking or laparoscopic biopsy.

Recurrence data

Recurrence dates are defined by ‘progression of disease’, a patient without improvement after treatment or while on maintenance therapy based on CT scan. Improvement or lack thereof is determined on the basis of CT scan impressions (for example, an increase in a lymph node or unchanged tumour implants). We define patients as ‘alive with disease’ if they have not achieved remission but have also opted out of new treatment lines and/or are on observation.

DLP scWGS processing

The Mondrian scWGS suite of tools and pipeline was used for processing of the scWGS. This scWGS dataset is a subset of the dataset used in ref. ¹⁸, see this publication for full details of the data generation and processing. Here we describe the processing in brief. Sequencing reads were aligned to hg19 using BWA-MEM. Read counts were calculated in 500-kb bins across the genome and corrected for GC bias, then these values were input into HMMcopy to infer integer copy number (ranging from 0 to 11). We then applied the cell quality classifier described in ref. ² and removed any cells with quality less than 0.75. In addition, we removed replicating cells, multiplet cells and cells suspected to be the result of multipolar divisions, see ref. ¹⁸ for a detailed description of the filtering criteria. We then applied SIGNALS v.0.7.6 to infer haplotype specific copy number using default parameters.

SNV calling in DLP

To detect SNVs in each dataset, reads from all cells from a patient analysed using DLP+ were merged to form ‘pseudobulk’ bam files. SNV calling was performed on these libraries individually using Mutect2. A panel of normals was constructed by identifying normal cells from every patient, merging them and then running the mutect2 panel of normal option. Mutect2 filter was used to filter variants. We then ran Articull to remove artefacts that are specific to DLP+ due to the shorter than average insert size⁴³. This filtered set of variants were then genotyped in individual cells using cellSNP v.1.2.2 (ref. ⁴⁴).

SV calling in DLP

To detect SVs in each dataset, we also used the merged pseudobulk bam files. LUMPY⁴⁵ and deStruct⁴⁶ were run on these pseudobulk libraries. Events were retained if they were detected in deStruct and could be matched in the LUMPY calls. Breakpoint predictions were considered matched if the positions involved were each no more than 200 nucleotides apart on the genome and the orientation was consistent.

SVs called in the pseudobulk library were then genotyped in single cells. To do this we used a modified version of SVtyper⁴⁷. One key modification was rounding the read count up rather down, the read count computation internally in SVtyper are MAPQ scores rather than counts so are non-integers before rounding and outputting to a vcf file. This change is necessary in single cells as typically we observe only a single read supporting a SV. SVtyper computes the number of reads that support the reference (these are reads that directly span the genome reference at the breakpoint locations) and the number of reads that support the alternate allele. Alternate allele counts are either split reads that directly sequence the breakpoint or discordant reads that have larger than expected insert sizes or align to different chromosomes in the case of translocations. Clipped reads that support the breakpoint are also computed, to be more conservative we did not include these reads in the total of SV supporting reads. We made an extra modification requiring split reads to match both sides of the breakpoint to contribute to read counts, in the default version, a split read aligning to one side of the breakpoint would contribute 0.5 counts. This option is available through the –both-sides command line option.

Phylogenetic inference and clone assignments

MEDICC2 was used to infer phylogenetic trees using haplotype specific copy number as input, see ref. ¹⁸ for further details. We then manually identified clades in the tree that were the ancestor of clade specific genomic features of interest. These included whole-genome doubling, whole chromosome and chromosome arm gains or losses and focal amplifications. Clones were then defined as the set of cells that were descendants of each clade of interest. Clones with their genomic features of interest can be found in Supplementary Table 5.

Clone-level 10-kb resolution copy-number calling

Once cells were assigned to clones, we also called integer copy number at 10-kb resolution at the clone level. Read counts were computed in 10-kb bins across the genome in every cell and then summed across cells assigned to each clone. Aggregated read counts were then normalized against the read counts from any normal diploid cells sequenced in the same library and then GC corrected using the same modal GC correction described in ref. ². These normalized GC corrected read counts were then adjusted for ploidy of the clone (inferred from the 500-kb data) and then we applied a HMM to compute integer read counts. The HMM model (code available at https://github.com/shahcompbio/HMMclone) uses a state space of 0–25 with each state assumed to be a normal distribution with standard deviation 0.2 and mean equal to the integer copy number. The standard deviation was determined empirically from the data. The viterbi algorithm was used to compute the most likely copy-number profile.

Identifying complex rearrangement processes

To identify complex structural events, we manually reviewed 10-kb clone-specific copy-number profiles overlayed with SVs for all patients and clones. We then assigned events as follows:

1. (1)

   chromothripsis: high density of rearrangement breakpoints and oscillating copy-number profiles (for example, 004 and 083)

2. (2)

   breakage-fusion bridge cycles: foldback inversion breakpoints coincident with copy-number changepoints and staircase such as copy-number profiles (for example, 045 and 051)

3. (3)

   Pyrgo: stacked tandem duplications resulting in tower-like copy-number profiles (for example, 081)

4. (4)

   complex inter-chromosomal rearrangements: high diversity of copy-number states with translocations between two or more chromosomes (for example, 002)

5. (5)

   ecDNA: focal high-level amplification with large copy-number diversity between cells (none detected).

In addition, oncogenic amplifications (whether they were part of an identifiable complex rearrangement process or not) were defined as high-level amplification (more than ten copies), amplifications (more than two times Ploidy). Homozygous deletions were identified by searching for bins with copy-number state = 0 in the 10-kb clone-level copy-number profiles. We therefore may miss homozygous deletion smaller than 10 kb, unless they were identified in MSK-IMPACT targeted sequencing, in which case they were probably truncal alterations.

Assigning SVs to clones

Assigning SVs to clones was done using the matrix of read counts per SV per cell, the cell to clone label mapping and the clone-level 10-kb copy-number profiles. First, we summed the SV supporting reads across clones giving an SV by clone matrix. Any SVs with non-zero read counts were assumed to be present in the clone. In addition, when SVs could be mapped to copy-number changepoints identified at 10-kb resolution, we checked whether there existed other clones that had the same copy-number changepoint but lacked read level support for the SV. In these cases, the SV was also assumed to present in that clone. This was to circumvent cases in which the total number of cells was too low to confidently assume the absence of a particular SV. In addition, we manually reviewed all SVs and copy-number profiles to produce a high confidence set of clone-specific SVs. Manual review involved inspecting the location of all SVs relative to copy-number profiles and ensuring absence of SV read support did not coincide with a proximal copy-number change. In some cases, no SVs could be found that were specific to a clone; this was largely due to clones being too small and consequently lacking the cumulative sequencing coverage to detect SVs in pseudobulks. In such cases, we used coarser clone definitions by merging sister clones resulting in clones with enough cells to call clone-specific alterations, typically at least 50 cells was needed.

Bulk WGS and MSK-IMPACT

The bulk WGS and MSK-IMPACT targeted sequencing was originally published in ref. ³. See this publication for data generation, processing and data access.

Probe design and synthesis

For most patients we aimed to include 1,000 genomic features encompassing SVs, SNVs and germline single-nucleotide polymorphisms (SNPs). For samples that constituted our pilot (patients 068, 065, 044 and 003) the number of features was lower, between 250 and 400 and included only a limited number of SNVs. Within the SV and SNV groups, these could be classed into truncal (present in every tumour cell) or subclonal (present in a fraction of tumour cells). The number of probes from each class was variable between patients due to differences in the number of SVs and SNVs called in each patient, as well as the clonal structure in each patient. We first required 200 truncal SVs and 200 truncal SNVs. The remaining 600 probes were split between subclonal SVs and SNVs. We ensured we had 200 subclonal SNVs and then the remaining slots were given to subclonal SVs; if there were still slots remaining then we included extra subclonal SNVs. Within the SNV class we included any SNV annotated as ‘High Impact’ in the MSK-IMPACT targeted sequencing. Probes were synthesized by Integrated DNA Technologies using the xGen MRD hybrid probes, from 120-bp sequences provided as FASTA files. We also included a small panel of germline SNPs to verify sample identity throughout the preparation process. No discrepancies were found.

cfDNA duplex sequencing analysis

We used the MSK-ACCESS protocol to generate the cfDNA sequencing data; this protocol is described in detail in ref. ¹³. The gene panel used in ref. ³ was replaced with the patient-specific probe sets. Patient probes from at least two patients were pooled together so that for each patient probe set we could estimate background error rates from the counts supporting SVs and SNVs in off-target patients.

To process the cfDNA sequencing we used a suite of tools developed by the Centre for Molecular Oncology informatics team at MSK for use with the MSK-ACCESS assay (https://github.com/msk-access). The nucleo pipeline was used to generate bam files from fastq files. The output of this pipeline is double strand error-corrected bam files (duplex), single strand (simplex) and uncorrected bam files that can then be used for downstream applications. Read counts of supporting and reference reads for SNVs and Indels were extracted using https://github.com/msk-access/GetBaseCountsMultiSample. This takes a MAF file as input and outputs a MAF file with extra columns for the read counts in duplex, simplex or uncorrected bam files. To extract read counts for SVs we used the same version of SVtyper modified for use with DLP+ described above. Discordant reads were only included in read count calculations if the SV size (end − start) was greater than 10⁴ or if it was a translocation. For a split read to be designated as supporting an SV, we required that alignments had evidence of both sides of the breakpoint to be included (option –both_sides). SVs in cfDNA sequencing were also visually inspected using IGV and svviz v.1.5.2 (ref. ⁴⁸).

Computing error rates in cfDNA

To compute error rates across sequencing types (duplex, simplex, raw uncorrected) and mutation types (SVs and SNVs) we applied the patient-specific probe set to many samples from off-target patients. We then summed the counts of reference supporting reads and variant-supporting reads for off-target variants across all patients, and defined the error rate as variants supporting reads divided by total number of reads.

We could then compare estimated error rates with the theoretical LOD on the basis of the number of mutations and coverage of the assay. We will define the LOD as the smallest tumour fraction that can be reliably detected; this can be defined as when the theoretical LOD (TLOD) is greater than the error rate of the assay. The TLOD is the LOD in the absence of any sequence error and depends on the number of genomic targets and the sequencing coverage. A framework for estimating this value is outlined in ref. ¹⁴ and given by:

where M is the number of detected mutations, N is the total number of mutations in the panel, TF is the tumour fraction and cov is the depth of coverage. As in ref. ¹⁴, when M < 1, we interpret it as a detection probability. We show this detection probability as a function of tumour fraction, coverage and the number of mutations in Extended Data Fig. 5a (this reproduces Fig. 1e from ref. ¹⁴).

We can then define the TLOD as the lowest TF with 99% probability of detection. So TLOD is given by (setting M = 0.99 and rearranging for TF):

Estimating clone frequencies

To estimate clone frequencies, we calculated variant allele frequencies (VAFs) for each clone by summing the total number of variant-supporting reads and dividing them by the total number of reads for all variants assigned to a clone. To mitigate against biases from variant copy number we developed an approach that leverages our ability to calculate accurate cancer cell fractions (CCFs) from our scWGS data. Using our single-cell copy-number phylogenies, we first identify the most-recent common ancestor (MRCA) of cells containing a SV and then calculate the tree-derived CCF, CCF_tree, by dividing the number of descendants of the MRCA (n_MRCA) by the total number of cells (n_total).

We then calculate the VAF of each SV from the pseudobulk BAM file. These two values give us a way to infer a correction factor, C, that we use to convert a given VAF measurement to a CCF value.

This is then applied to VAFs measured from cfDNA. To plot the changes in frequency over time we normalized corrected VAFs so that they summed to 1 at each time point, then applied a spline function to smooth values between time points. When no tumour DNA was detected, we allowed all clones to have VAF = 0. Smoothing was performed using the splinefun function in R with method = ‘monoH.FC’. This resulted in values that were greater than 1 or less than 0 in some cases, we therefore renormalized the data so that frequencies were positive and summed to 1 at each time point. In addition, when there were large periods of time preclinical recurrence without cfDNA samples we assumed tumour DNA was 0 (for example, in patient 107). We did not include clone frequency estimates when plasma tumour fractions were less than 10⁻⁴ (estimates based on truncal SVs), reasoning that clone frequency estimates at such low tumour fractions would be unreliable and suffer from dropout issues. Given the low error rates, we used the uncorrected raw sequencing for estimating clone frequencies using SVs. For estimating clone frequencies using SNVs we used the same approach but used duplex consensus sequences for read counting due to the higher error rates for SNVs. In patient 139, no SVs were identified in clone B, therefore we used the SNV based clone frequencies for this clone.

Identifying BRCA reversion mutations

For BRCA1/2 mutant cases, we also included probes that captured exonic regions within 200 bp of the mutation, enabling detection of proximal BRCA1/2 reversion mutations^49,50. We used revmut (https://github.com/inodb/revmut) to identify putative BRCA reversion mutations in the first instance. In addition, we inspected alignments in IGV around the BRCA mutations to look for any extra putative reversion mutations not identified by revmut. This is how we found the reversion mutation present in patient 009. This mutation was a large 1.37 kb deletion that excised the germline mutation, alignments with the same breakpoint sequence, aligning to the same locations were found in 3 postrecurrence samples. This mutation was probably not identified using revmut due to it being unusually large compared with previously reported BRCA reversion mutations.

Wright–Fisher modelling and hypothesis testing

To test for non-neutrality in clone frequencies over time we implemented a modelling and hypothesis testing framework based on a multi-species Wright–Fisher model with varying population size. Population size was assumed to be 10⁹ at the time of surgery (t = 0) and then varied according to CA-125 concentrations. We set the population size at the time point with the lowest CA-125 level N_low = 10⁴, assuming this was the period with the smallest tumour cell population. We then set the population size (N) to vary exponentially according to the following equation:

where \(A={N(0)\times {\rm{e}}}^{-b\times {\rm{CA}}125(0)}\) and \(b=\frac{\log ({N}_{{\rm{low}}})-N(0)}{{{\rm{CA}}125}_{{\rm{low}}}-{\rm{CA}}125(0)}\). We then used the multinomial distribution to simulate clone frequencies over time:

where X_1…k(i + 1) is the population size of each clone k in generation i + 1, N(i) is the total population size in generation i and p_1…k(i) are the clone frequencies in generation i for the k clones. For generation i = 1, N(i = 1) = 10⁹ and p_1…k(i = 1) are given by the clone frequencies estimated from cfDNA at t = 0. We then forward simulate this process for the clinical time course of each individual patient 1,000 times giving a distribution of clone frequencies at t_end. We assumed a generation time of 4 days⁵¹, and t_end was set to be the final cfDNA time point in each patient. We then calculated a z-score, comparing the observed clone frequency from data to the mean and standard deviation of the simulated frequencies to calculate a P value for each clone under the hypothesis of neutral evolution.

cfDNA WGS

Whole-genome libraries constructed during the duplex sequencing assay library prep were sequenced to 70× on an illumina NovaSeq using 100 bp reads. Reads were mapped to hg19 using BWA-MEM⁵². Read counts in 0.5-Mb bins across the genome were calculated and GC corrected using QDNAseq⁵³. To compare this to data from the duplex sequencing targeted assay, we used information from the DLP copy-number profiles and the clone fractions inferred from CloneSeq-SV to predict copy-number profiles measured with cfDNA WGS. Copy-number ratios (R) in bin i are given by:

where n is the normal fraction, c_i is the copy number in bin i and p is the ploidy of the tumour. We know n from the TP53 VAF in cfDNA, for c_i and p we took the weighted average across clones, with weights given by the estimated clone fractions at each time point.

Heterozygous SNPs identified in the scWGS pipeline were genotyped in the cfDNA WGS. Phasing information from SIGNALS⁹ was used to calculate phased B-allele frequencies (BAFs) in 0.5-Mb bins, resulting in BAF measurements with identical phasing between scWGS and cfDNA WGS. To assess similarity between cfDNA WGS and clone copy-number profiles, we calculated Pearson correlation coefficients between cfDNA WGS and all clones for both read-depth ratios and BAF values, the summed values of these two values were used to identify the most similar clone for each time point.

scRNAseq data generation and processing

The scRNAseq data were originally published in ref. ³; full details of the processing can be found in that study. Pathway scoring was performed with PROGENY⁵⁴ or the Seurat module scoring function using hallmark pathways.

TreeAlign

To match scRNAseq cells to clones identified in DLP we used TreeAlign³³. To do this, we genotyped the same set of heterozygous SNPs used to call allele-specific copy number in DLP+ in scRNAseq using cellSNP⁴⁴. The per-cell SNP count matrix was then input into TreeAlign along with clone assignments and 10-kb clone copy-number profiles derived from DLP. We used the CloneAlignClone method and used default parameter values apart from min_clone_assign_prob = 0.5. To compare transcriptional heterogeneity for each clone, we took the mean value of the per-cell Seurat-derived module score for all hallmark pathways, then for each patient calculated the maximum value minus the minimum value. These per-patient difference values were then plotted as violin plots ordered by the average difference across the cohort of patients. For patients with longitudinal cfDNA measurements, scRNAseq data were available for patients 107, 014, 045, 009, 002, 026, 083 and 037, however, following application of TreeAlign, in some patients, clones were represented minimally due to differences in data collection from different sites. In patients 107 and 009, all clones were represented with at least 100 cells present from each clone, we therefore focused on these cases when comparing drug-resistant to drug sensitive clones.

DNA FISH

DNA FISH was performed using locus-specific probes for CCNE1 (green) and a centromeric probe for chr. 19 (Cen19; orange) as an internal control (Empire Genomics). Representative regions for FISH imaging were selected on the basis of pan-cytokeratin immunofluorescence staining and on haematoxylin and eosin staining performed on adjacent serial sections, to enrich for epithelial tumour areas. Tumour histology annotations of immunofluorescence and haematoxylin and eosin images are described in refs. ³ and ¹⁸. For quality control, only cells with at least one signal from each probe were retained, ensuring the centromeric signal served as an internal control for ploidy and hybridization efficiency.

Statistical analysis

For between-group comparisons, we used two-sided t-tests. To assess correlations, we used Pearson correlation. All statistical analysis was done using R.

Data organization

To facilitate integration of data across several modalities we used the isabl platform⁵⁵. Isabl is a databasing and data access platform that allows users to straightforwardly link many datasets from the same patient and chain together pipelines across modalities.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.