Modeling clonal structure over narrow time frames via circulating tumor DNA in metastatic breast cancer

Abstract

Background: Circulating tumor DNA (ctDNA) offers minimally invasive means to repeatedly interrogate tumor genomes, providing opportunities to monitor clonal dynamics induced by metastasis and therapeutic selective pressures. In metastatic cancers, ctDNA profiling allows for simultaneous analysis of both local and distant sites of recurrence. Despite the promise of ctDNA sampling, its utility in real-time genetic monitoring remains largely unexplored.

Methods: In this exploratory analysis, we characterize high-frequency ctDNA sample series collected over narrow time frames from seven patients with metastatic triple-negative breast cancer, each undergoing treatment with Cabozantinib, a multi-tyrosine kinase inhibitor (NCT01738438, https://clinicaltrials.gov/ct2/show/NCT01738438 ). Applying orthogonal whole exome sequencing, ultra-low pass whole genome sequencing, and 396-gene targeted panel sequencing, we analyzed 42 plasma-derived ctDNA libraries, representing 4-8 samples per patient with 6-42 days between samples. Integrating tumor fraction, copy number, and somatic variant information, we model tumor clonal dynamics, predict neoantigens, and evaluate consistency of genomic information from orthogonal assays.

Results: We measured considerable variation in ctDNA tumor faction in each patient, often conflicting with RECIST imaging response metrics. In orthogonal sequencing, we found high concordance between targeted panel and whole exome sequencing in both variant detection and variant allele frequency estimation (specificity = 95.5%, VAF correlation, r = 0.949), Copy number remained generally stable, despite resolution limitations posed by low tumor fraction. Through modeling, we inferred and tracked distinct clonal populations specific to each patient and built phylogenetic trees revealing alterations in hallmark breast cancer drivers, including TP53, PIK3CA, CDK4, and PTEN. Our modeling revealed varied responses to therapy, with some individuals displaying stable clonal profiles, while others showed signs of substantial expansion or reduction in prevalence, with characteristic alterations of varied literature annotation in relation to the study drug. Finally, we predicted and tracked neoantigen-producing alterations across time, exposing translationally relevant detection patterns.

Conclusions: Despite technical challenges arising from low tumor content, metastatic ctDNA monitoring can aid our understanding of response and progression, while minimizing patient risk and discomfort. In this study, we demonstrate the potential for high-frequency monitoring of evolving genomic features, providing an important step toward scalable, translational genomics for clinical decision making.

Keywords: Circulating tumor DNA; Liquid biopsy; Neoantigens; Serial sequencing; Targeted panel sequencing; Tumor evolution; Ultra-low pass whole genome sequencing; ctDNA.

Fig. 1

Study design and sampling dynamics. a Schematic diagram of the analysis workflow from patient selection, sample capture, and sequencing to downstream analyses. We leveraged the Terra Genomics/FireCloud platform for data storage and high-performance computing tasks. b Schematic representation of sampling density for each of the seven cohort members on study, also specifying whether whole exome sequencing and/or targeted panel sequencing was performed on that sample. All samples received ultra-low-pass whole genome sequencing. c Tumor fraction dynamics colored by individual. Tumor fraction was measured on study using ultra-low-pass whole genome sequencing and the ichorCNA algorithm. d Tumor fraction dynamics recolored by RECIST v1.1 response by imaging categories. RECIST v1.1 bucket response type into several categories: complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD)

Fig. 2

Orthogonal ctDNA sequencing approaches are highly concordant. Somatic SNV and INDEL calling of whole exome sequencing (WES; average depth 150X) and targeted panel sequencing (TPS; nominal sequencing depth 10,000X) were completed on the Terra/Firecloud platform using gatk-Mutect2 pipelines (McKenna et al., 2010). a Variant recall assessment of TPS on somatic variants discovered in one or more WES assays. Only variants intersecting theoretical capture regions of TPS were considered. Variants used in assessment were those called in WES at any point, which also overlapped in genomic position with target or bait regions included in the TPS. X’s indicate a lack of adequate sequencing depth in the TPS. Center and right panels compare variant allele frequency (VAF) data from each assay. b Scatter plot comparing estimated VAF in TPS and WES sequencing across all individuals and time points. 1:1 line drawn for reference. c WES and ULP-WGS based algorithmic estimates of sample purity (a.k.a. tumor fraction) across samples and time points with high tumor fraction (TFx > 10%). d Algorithm estimation of ploidy (averaged copy number state across genome) across WES and ULP-WGS-based methods at time points with high tumor fraction. ABSOLUTE Soln.1 and Soln.2 represent the top two proposed solutions by model likelihood (Included here, as ABSOLUTE often suggests manual curation and/or override of the top solution)

Fig. 3.

Copy number profiles are stable. Ultra-low pass whole genome sequencing (ULP-WGS) was performed on all 42 ctDNA samples and tumor fraction and copy number data derived using ichorCNA. a Genome-wide copy profile of patient RP-466, derived from ULP-WGS on liquid biopsy ctDNA, showing changes in focal event resolution resulting from shifts in tumor fraction. Dark green segments represent a copy number of 1; blue represent neutral or 2 copies, brown and red represent 3 and 4+, respectively. b Scatter plot of computed log-ratios in ULP-WGS, compared to those derived from WES or TPS data using binned read-count of on and off target bins. c Discrete copy number confusion matrix for ULP-WGS based calls at first and last time points. All samples had tumor fraction ≥10%. Genomic positions assayed between first and last time points were uniformly and randomly sampled, and discrete copy number states were capped between one and seven during initial ichorCNA analyses

Fig. 4

Tumor subclonal dynamics vary across patients. Models of clonal and subclonal populations which make up the cancers of metastatic patients, derived using PyClone [34]. Variant inputs include union of filter-passing alterations from each sampled time point delivered by the commercially available liquid-biopsy targeted panel-sequencing pipeline at the Broad Institute. Copy number information and purity were derived from ichorCNA. a, b Clonal prevalence dynamics, clustering, and inferred phylogenetic tree structure for patient RP-466, revealing generally unchanging populations in the tumor, with important drivers occupying early positions in cell lineages. c, d RP-527 clonal dynamics profile and inferred tree structure showing statistically significant clonal expansion of cell lineage marked by non-synonymous DDR2 and RNF43 variants. e, f RP-557 profile and tree showing the opposite trend as RP-527, with a decreasing cell population marked by RB1 mutation

Fig. 5

Whole exome sequencing uncovers driver mutations and allows neoantigen prediction. Whole exome sequencing results from 31 total samples with tumor fraction ≥10% using short variant and INDEL calling tools from gatk-Mutect2 pipelines (McKenna et al., 2010), with subsequent neoantigen binding predictions for known MHC molecules from NetMHCpan 4.0 (Reynisson et al., 2020). a Driver mutations found via whole exome sequencing across time points. Variant data visualized are those whose genes have been previously annotated in literature as breast cancer drivers or pan cancer drivers. b Trends in predicted neoantigens among cohort members. Strong binders are denoted as those peptide sequences with NetMHCpan ranks <0.5%, and weak binders are those with ranks <2%. Neoantigen Generating sSNV are alterations whose changes to peptide structure are predicted to produce neoantigens capable of strong or weak binding to known MHC molecules. c, d Neoantigen dynamics from patient RP-527 and RP-535, showing proportions of detected neoantigens and dropout over time. Strong, weak, and ND labels correspond to binding affinity of predicted neoantigens, as well as a non-detected category to capture dropout. Threads are colored by their state at the final sequencing time point