Genome-wide cell-free DNA mutational integration enables ultra-sensitive cancer monitoring

Abstract

In many areas of oncology, we lack sensitive tools to track low-burden disease. Although cell-free DNA (cfDNA) shows promise in detecting cancer mutations, we found that the combination of low tumor fraction (TF) and limited number of DNA fragments restricts low-disease-burden monitoring through the prevailing deep targeted sequencing paradigm. We reasoned that breadth may supplant depth of sequencing to overcome the barrier of cfDNA abundance. Whole-genome sequencing (WGS) of cfDNA allowed ultra-sensitive detection, capitalizing on the cumulative signal of thousands of somatic mutations observed in solid malignancies, with TF detection sensitivity as low as 10-5. The WGS approach enabled dynamic tumor burden tracking and postoperative residual disease detection, associated with adverse outcome. Thus, we present an orthogonal framework for cfDNA cancer monitoring via genome-wide mutational integration, enabling ultra-sensitive detection, overcoming the limitation of cfDNA abundance and empowering treatment optimization in low-disease-burden oncology care.

Figures

Extended Data Fig. 1. Deep targeted sequencing analysis and DNA extraction optimization

(a-b) constitute a re-analysis of previously published data. (a) Histogram of samples detected as a function of the maximum variant allele fraction (VAF) in stage I-IV cancer. Maximum VAF serves as an estimation of tumor fraction (TF). Patient histogram (solid bars) and empirical cumulative distribution (solid line) are shown, displaying stage dependent detection probability from 50% sensitivity in stage I cancer to >90% sensitivity in stage IV. Color coding represents different cancer stages. Number of patients considered per cancer stage is indicated in the figure. (b) Histogram of variants (mutations) detected as a function of the variant allele fraction (VAF) for different cancer stages. Higher cancer stage exhibits an increase in the mutation VAF, associated with higher probability of detection. Number of mutations included in the analysis per cancer stage is indicated in the figure. (c) Comparison of cfDNA yields across various leading commercial extraction kits. Multiple extractions (n = 3) for each kit were performed over 1mL aliquots of the same input plasma sample (from plasmapheresis of a normal donor), \and total mass (per 1mL) was determined with Qubit (ThermoFisher, Waltham, MA). Mean value and confidence interval (standard deviation) is shown. (d) Extracted cfDNA was evaluated with Nanodrop (ThermoFisher, Waltham, MA) to detect impurities such as salt and genomic DNA. Omega Bio-Tek (Norcross, GA) had the lowest amount of contaminant carry over compared to other kits. Multiple tests (n = 3, solid lines) for each kit were performed over 1mL aliquots of the same input plasma sample

Extended Data Fig. 2. Patient-specific genome-wide SNV integration and error suppression

a) Base quality (BQ) signal for high quality germline single nucleotide polymorphisms (SNPs; n = 9,142,326) vs. low VAF (single supporting read) artifactual variants (n = 407,061) from representative PBMC sample (Pat.01), showing distinct separation between SNP and the low VAF artifacts BQ distributions (P value < 10−100, two-sample t-test), and supporting effective filtration of sequencing artifacts by BQ filtration. (b) Receiver-operating-curve (ROC) analysis for BQ filtering including germline SNPs (true labels) and low VAF artifactual variants (false labels), using the same data set as in (a), and showing high filter performance for this simple quality metric (AUC = 0.9, n = 9,142,326). (c) Variant position-in-read (PIR) shows association between low VAF artifactual variants (n = 407,061) and position at the 3’ of the sequencing read, while germline SNPs (n = 9,142,326) show uniform spread across the sequencing read length. (d) Support vector machine (SVM) classification performance between germline SNPs (random subsampling, n = 100,000) vs. low VAF artifactual variants (random subsampling, n = 100,000) from all PBMC WGS data (n = 8). Performance of SVM and random forest classification was compared over the same sample set with 10-fold cross validation. (e) Box plot of error rate estimations before error suppression (blue) and after SVM-based error suppression (red) over four cancer types and 40 PBMC-derived replicates (2 patients per cancer type, 20 replicates per patient’s PBMC sample). Error rate was calculated as the number of mismatches detected over the number of bp checked. Showing a uniform error reduction (median 14 fold-change reduction, range 11–17). (f-h) Patient-specific SNV signal-to-noise quantification over a range of TFs (10−5-10−2) compared to basal noise signal detected in control (TF = 0, subsampled PBMC DNA fragments) samples (left column). Signal-to-noise was estimated by calculating the log difference between the number of detections in each plasma-like admixture (TF > 0) and the mean number of detections in the controls (TF = 0). Analysis was done separately using tumor and matched germline (PBMC) WGS from lung (f, Pat.04), breast (g, Pat.05) and osteosarcoma (h, Pat.08) patients. Inset panel shows discrimination of tumor and control samples down to tumor fraction 10−5 after utilizing machine-learning-based sequencing error suppression (red) vs. reduced sensitivity with the raw unfiltered data (blue). (i) Benchmarking of mutation detection performance for mutation centric method vs. read-centric method (MRDetect). Patient-specific SNV signal-to-noise quantification over a range of TFs (10−5-2X10−1) compared to basal noise signal detected in control (TF = 0, subsampled PBMC DNA fragments) samples (right column). Signal-to-noise was estimated by calculating the log difference between the number of detections in each plasma-like admixture (TF > 0) and the mean number of detections in the controls (TF = 0). (j) Single nucleotide variant (SNV) point mutation detection in plasma mixtures with different tumor fractions (TF > 0) and controls (TF = 0) is shown. Y-axis shows the number of detections (variants observed in tumor WGS and also detected in plasma synthetic admixture) as a function of TF (x-axis). Red line constitutes the number of detections predicted for each TF based on the mutation load, coverage, noise model. Gray area represents the area under the background noise model threshold (1.5std), showing robust discrimination from noise for TF > 10-3. Analysis was done on 35X coverage lung cancer (Pat.04) admixture cohort. Centre values represent mean and error bars represent standard deviation. In (f-j), n = 11 independent admixture samples for TFs > 0 and n = 20 independently down-sampled PBMC replicates for the control (TF=0) of each patient. Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR.

Extended Data Fig. 3. Patient-specific genome-wide SNV integration provides accurate read out of tumor fraction

(a-h) Tumor fraction (TF) inference plots using MRDetect genome-wide SNV integration for two melanomas (a: Pat.01, b: Pat.02), two lung cancers (c: Pat.03, d: Pat.04), two breast cancers (e: Pat.05, f: Pat.06), and two osteosarcomas (g: Pat.07, h: Pat.08). Each plot was generated using in silico admixtures of varying TF (range 10−5-0.2) with 35X depth of coverage by randomly downsampling and mixing tumor reads (mean coverage 97X, range 85X-110X) and germline reads (mean coverage 49X, range 43X-56X) from WGS data (see Methods; Supplementary Table 1). For TFs > 0, n = 11 independent admixture samples. For the control (TF = 0), n = 20 independently down-sampled PBMC replicates. We observed accurate TF estimation as low as 10−5, discriminated from control (TF = 0) samples (left box-plot), with high Pearson correlation (two-sided test) between the input TF mixture (x-axis) and the SNV-based estimated TF prediction, confirming accurate inference based on genome-wide mutational integration. Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR.

Extended Data Fig. 4. CNA load across tumor types

Histogram of CNA load in WGS samples across cancer types from the TCGA cohort. Measured as a function of the size of genome altered by CNA (in log10Mb). Dashed lines represent the percentage of samples that have CNA load of over 10 Mb and 1 Gb respectively, for each cancer type. Cancer types include LUSC: Lung squamous cell carcinoma (n=50), HNSC: Head and Neck squamous cell carcinoma (n=50), CESC: Cervical squamous cell carcinoma and endocervical adenocarcinoma (n=18), OV: Ovarian serous cystadenocarcinoma (n=50), KICH: Kidney Chromophobe (n=50), COAD: Colon adenocarcinoma (n=53), THCA: Thyroid carcinoma (n=50), LUAD: Lung adenocarcinoma (N=152), ESCA: Esophageal carcinoma (n=19).

Extended Data Fig. 5. CNA based detection of ctDNA tumor fraction

(a-c) Tumor fraction (TF) inference using MRDetect genome-wide CNA integration for representative patients, including melanoma (a: Pat.01) and lung cancer (b: Pat.03, c: Pat.04). Each plot was generated using in silico admixtures of varying TF (range 10−5 – 0.2) in 18X coverage, by randomly downsampling and mixing tumor reads (mean coverage 97X, range 85X-110X) and germline reads (mean coverage 49X, range 43X-56X) from WGS data (see Methods; Supplementary Table 1). Twenty replicates used for each TF > 0 sample, and for the control (TF = 0) samples, showing accurate TF estimation as low as 5*10−5, discriminated from control (TF = 0) samples (left box-plot), with high Pearson correlation between the input TF mixture (x-axis) and the CNA-based estimated TF prediction (two-sided test). (d-f) Tumor fraction (TF) inference in neutral regions (no copy number gain or loss in the tumor WGS data) for the same in silico admixtures (d: melanoma, Pat.01; e: lung, Pat.03; f: lung, Pat.04), shows the expected low Pearson correlation between input admixture TF and the signal (two-sided test), consistent with no expected coverage changes in the plasma admixtures in these regions. Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR.

Extended Data Fig. 6. CNA and SNV MRDetect correlation and further error suppression

(a) Spearman correlation between SNV and CNA TF estimation across TF admixtures for a lung tumor (Pat.03) shows high correlation (two-sided test) between the two orthogonal inference methods. Red dots correspond to cancer plasma (TF > 0) samples and blue dots correspond to control plasma (TF = 0) samples. Detection threshold (dashed lines) were set on TF < 5*10−5 for both methods. Eleven replicates used for each TF > 0 sample, and 20 replicates were used for the control (TF = 0) samples. (b) Comparison to an orthogonal CNA-based TF method- ichor-CNA. Analyzing the same cohort of breast cancer (Pat.05) in silico synthetic admixture samples shows concordance in TF estimation for high TF (TF > 5*10−3), with extension of detection for MRDetect to lower TFs. The same 20 replicates were used for both MRDetect and ichor-CNA for each of the TF > 0 and control (TF = 0) samples. (c) Proportion of variant concordant (brown) vs. discordant (gray) read pairs (R1 and R2) detected in germline SNPs and artifactual variants. Analysis was done across 10 control (benign lung lesions) plasma samples, comparing read-pairs associated with germline SNPs (right bar) vs. read-pairs associated with artifactual variants (left bar) per plasma sample. The artifactual variants were defined by read pairs with variants overlapping the union of all patient somatic SNVs compendia across all LUAD patients, that were observed with the same variant in the control plasma sample. The number of read-pairs used in the analysis is indicated above each bar. (d) Median genome-wide normalize (divided by mean coverage) coverage from matched germline PBMC WGS samples from patients with LUAD (cyan, n = 15) and control plasma WGS samples from patients with benign lung lesions (red, n = 11) before and after robust Z-score normalization. (e) Median absolute deviation (MAD) calculated over normalized-coverage (i.e., divided by mean coverage) from matched germline PBMC and control plasma WGS samples as in (d) before and after robust Z-score normalization. Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR.

Extended Data Fig. 7. Serial application of MRDetect to monitor melanoma response to immunotherapy

Melanoma treatment response (Patient MEL02) during immunotherapy (Nivolumab) is monitored by blood samples. Upper panels- Treatment monitoring by computed tomography (CT) shows response to therapy but residual disease after 3 months of therapy. Middle panel- MRDetect Z-scores effectively track tumor responses, matching radiographic changes, in higher temporal resolution than that feasible with imaging. Lower panel- ichor-CNA captures treatment response dynamics but showing lower signal to noise ratio compared to the MRDetect method. For both MRDetect and ichor-CNA methods, Z-score is calculated from a single plasma sample for each timepoint compared to a panel of control samples (n = 30). Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR.

Extended Data Fig. 8. MRDetect performance in colorectal monitoring

(a) Robust Z-score discrimination between signal detected across 20 random subsamplings (80% of reads per subsampling iteration) of LUAD patient pre-operative plasma (black, n = 36 patients) and the cohort of control plasma test set (gray, n = 30). The signal was measured on the subsampling set and control test set using the same patient-specific point mutation (SNV) compendium. Z-score was calculated using the noise parameters estimated in the control test cohort (see Methods). (b) Receiver-operating-curve (ROC) analysis was performed over all SNV-based Z-score values calculated on the patients’ pre-operative plasma and control plasma as in (a). (c) Cross patient noise evaluation. Robust Z-score discrimination between signal detected at 20 random subsamplings (80% of reads per subsampling iteration) of LUAD patient pre-operative plasma WGS (black, n = 36 patients), cross-patient noise estimation via application of the patient-specific compendium to all other patient pre-operative plasma (n = 35, gray). Z-score was calculated using the noise parameters estimated in the cross-patient cohort (see Methods). (d) Receiver-operating-curve (ROC) analysis was performed over all SNV-based Z-score values calculated on the matched patients and cross-patient plasma. (e) Z-score discrimination between MRDetect-CNA on LUAD patient pre-operative plasma (red, n = 36 patients) compared to signal detected in neutral regions (as a negative control, blue), control plasma test cohort (n = 30) and cross-patient cohort (n = 35). Cross-patient noise was estimated by applying the patient-specific CNA compendium to other patient plasma samples (n = 35, all other patients). Z-score was calculated using the noise parameters estimated by the control plasma cohort. (f) Receiver-operating-curve (ROC) analysis was performed over all CNA-based Z-score values calculated on the patients’ pre-operative plasma and control patients. (g) ROC analysis was performed over all ichor-CNA TF values calculated on the LUAD patients’ pre-operative plasma and control patients (n = 66). Interestingly the two patient plasma samples detected by ichor-CNA included events that do not appear in the tumor, one of them was found to be a PBMC specific somatic event (potentially from clonal hematopoiesis). Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR. (h) Z-score discrimination between signal detected in 20 random subsampling (80% of reads per subsampling iteration) of LUAD patient plasma WGS (n = 22 patients) collected at a median of 17 days after surgery and a cohort of control plasma test samples (gray, n = 30). The signal was measured on the matched plasma and control set using the same patient-specific point mutation (SNV) compendium. Z-score was calculated using the noise parameters estimated in the control cohort (see Methods). (i) Z-score discrimination between MRDetect-CNA on LUAD patient plasma (red, n = 22 patients) collected at a median of 17 days after surgery, compared to signal detected in neutral regions (as a negative control, blue) and control plasma cohort (n = 30). Z-score was calculated using the noise parameters estimated by the control plasma cohort (see Methods). Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR.

Extended Data Fig. 9. MRDetect performance in LUAD monitoring

(a) Robust Z-score discrimination between signal detected across 20 random subsamplings (80% of reads per subsampling iteration) of LUAD patient pre-operative plasma (black, n = 36 patients) and the cohort of control plasma test set (gray, n = 30). The signal was measured on the subsampling set and control test set using the same patient-specific point mutation (SNV) compendium. Z-score was calculated using the noise parameters estimated in the control test cohort (see Methods). (b) Receiver-operating-curve (ROC) analysis was performed over all SNV-based Z-score values calculated on the patients’ pre-operative plasma and control plasma as in (a). (c) Cross patient noise evaluation. Robust Z-score discrimination between signal detected at 20 random subsamplings (80% of reads per subsampling iteration) of LUAD patient pre-operative plasma WGS (black, n = 36 patients), cross-patient noise estimation via application of the patient-specific compendium to all other patient pre-operative plasma (n = 35, gray). Z-score was calculated using the noise parameters estimated in the cross-patient cohort (see Methods). (d) Receiver-operating-curve (ROC) analysis was performed over all SNV-based Z-score values calculated on the matched patients and cross-patient plasma. (e) Z-score discrimination between MRDetect-CNA on LUAD patient pre-operative plasma (red, n = 36 patients) compared to signal detected in neutral regions (as a negative control, blue), control plasma test cohort (n = 30) and cross-patient cohort (n = 35). Cross-patient noise was estimated by applying the patient-specific CNA compendium to other patient plasma samples (n = 35, all other patients). Z-score was calculated using the noise parameters estimated by the control plasma cohort. (f) Receiver-operating-curve (ROC) analysis was performed over all CNA-based Z-score values calculated on the patients’ pre-operative plasma and control patients. (g) ROC analysis was performed over all ichor-CNA TF values calculated on the LUAD patients’ pre-operative plasma and control patients (n = 66). Interestingly the two patient plasma samples detected by ichor-CNA included events that do not appear in the tumor, one of them was found to be a PBMC specific somatic event (potentially from clonal hematopoiesis). Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR. (h) Z-score discrimination between signal detected in 20 random subsampling (80% of reads per subsampling iteration) of LUAD patient plasma WGS (n = 22 patients) collected at a median of 17 days after surgery and a cohort of control plasma test samples (gray, n = 30). The signal was measured on the matched plasma and control set using the same patient-specific point mutation (SNV) compendium. Z-score was calculated using the noise parameters estimated in the control cohort (see Methods). (i) Z-score discrimination between MRDetect-CNA on LUAD patient plasma (red, n = 22 patients) collected at a median of 17 days after surgery, compared to signal detected in neutral regions (as a negative control, blue) and control plasma cohort (n = 30). Z-score was calculated using the noise parameters estimated by the control plasma cohort (see Methods). Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR.

Extended Data Fig. 10. Imaging of sample monitored LUAD cases and fragment length analysis

Positron emission tomography–computed tomography (PET-CT) of two patients (LUAD#3 and LUAD#6) confirms no radiographically observable metastatic spread at the time of surgery (a and c, respectively), while metastatic recurrence has been identified approximately six months post-operative by PET-CT (b and d, respectively). (e-h) Representative fragment size histograms showing the distribution of DNA fragments as a function of the fragment length. DNA fragments that are associated with tumor mutations (gray) are showing significantly shorter size compared to DNA fragments that are associated with artifactual non patient-specific detections (red, derived from applying cross-patient mutational compendia; median P value < 10−3, two-sample t-test). Plots include two pre-operative plasma samples (e: LUAD#18, tumor mutation detection n = 35, artefactual detection n = 18802 ; f: LUAD#31, tumor mutation detection n=82, artefactual detection n = 32530) and matching plasma samples after surgery (g: LUAD#18, tumor mutation detection n = 55, artefactual detection n = 22737; h: LUAD#31, tumor mutation detection n = 27, artefactual detection n = 14610). (i) Kernel-density-estimator (KDE) trained to discriminate between tumor-derived (human aligned reads from a patient derived xenograft model, see Methods, blue) and normal-derived (from control plasma samples, orange) cfDNA based on the fragment size signature. The log difference between the tumor and normal density functions (black solid line) was used as a score function that integrates the fragment size shift signal across the entire fragment sizes distribution (80bp-600bp). (j) Pre-operative cfDNA showed significant shift (two-sample t-test) in their tumor-specific mutation detections in the patient plasma (red, n=563), compared to non-tumor (cross-patient) detected mutations in the same samples (blue, n=4184). Violin plots depict kernel density estimates of the density distribution. Center dashed lines represent the median and dashed lines represent the interquartile range.

Figure 1:. Low cfDNA input material limits sensitive ctDNA mutation detection with deep targeted sequencing, but can be overcome by genome-wide mutational integration.

(a) Illustrative schematic showing that mutations captured by targeted panel deep sequencing are only a small fraction of genome-wide somatic tumor mutations (upper panel). The low input and random sampling nature of cfDNA collection by blood samples results in higher probability for mutation capture by WGS than by targeted panels (lower panel). (b-d) constitute a re-analysis of previously published data. (b) Histogram of samples detected as a function of the maximum variant allele fraction (VAF) in stage IV (left) and stage I (right) cancer. Maximum VAF serves as an estimation of tumor fraction (TF). Patient histogram (solid bars), empirical cumulative distribution (solid line) and predicted “optimal-sensitivity” detection (dashed line, see Methods) are shown. Lower limit of detection (red solid line) shows a residual of 44% of stage I patients that are not detected by the deep targeted sequencing panel (10% for stage IV). Most of these cancers are predicted to have been detected if the VAF lower-limit-of-detection was 10-5. Re-analysis across 64 stage I and 21 stage IV patient samples across tumor types. (c) Plasma input material (depth of coverage) correlates to the number of detected tumor mutations (correlation and statistical significance calculated with a two-sided Pearson correlation.) Upper panel: distribution (histogram) of unique genomic equivalents (GEs, unique coverage, x-axis) across all patient samples (n = 194 patients, across different tumor types and disease stages) sequenced with exhaustive deep sequencing (median 40,729X). A median of 6,000 (range 400–14,000) unique reads (GEs) was found in a median 4.2mL of plasma (range 1.6–9.8). Lower panel: the average number of mutations detected per patient shows a positive correlation with the number of unique GEs in the patient plasma, indicating that input material is a potential limiting factor for mutation detection with deep targeted sequencing (two-sided Pearson correlation). Dashed line shows regression line, shaded area represents 95% confidence interval. (d) Patient samples with higher number of unique GEs (coverage) have higher probability for cancer detection by deep targeted sequencing (P value = 7*10−4, two-sided Wilcoxon rank-sum test). Analysis includes all detected (n = 85) and non-detected (n = 53) stage I/II patients across tumor types. Boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR. (e) Predicted probability (analytic derivation from binomial distribution of Bernoulli trials, see Methods) of detecting at least one mutated (SNV supporting) read as a function of the number of unique GEs (coverage) and the fraction of ctDNA in cfDNA (tumor fraction, TF). Left: targeting only a single SNV (e.g., with digital droplet PCR), the probability for detection decreases rapidly with tumor fraction (e.g., probability to detect at 10−5 with 5000 GEs is <0.05, red line). Middle: targeting 10 patient-specific mutations (e.g., 1Mbp panel for patient with 10/Mbp mutation load or a typical exome-based patient specific panel) shows low probability of detection in low TF (e.g., probability to detect at 10−5 with 5000 GE is <0.4 , red line). Right: targeting 10,000 patient-specific mutations (e.g., whole genome sequencing for patient with >4/Mbp mutation load) shows significant probability of detection in low tumor fraction even with modest coverage (e.g., probability to detect at 10−5 is > 0.99 with 50X coverage, red line).

Figure 2:. Patient-specific genome-wide SNV integration provides ultra-sensitive ctDNA detection and precision tumor fraction (TF) estimation.

(a) Illustration of genome wide SNV integration for inference of plasma TF, and analytical validation scheme. Patient-specific SNV mutations are detected from matched tumor and peripheral blood mononuclear cells (PBMC) germline whole genome sequencing (WGS). Tumor somatic SNVs are then used to calculate the genome-wide tumor signal in the patient’s plasma sample. To test this approach, for each cancer, tumor and PBMC WGS reads were mixed to generate a range of TFs (10−5-0.2) in 35X coverage with multiple replicates (n = 11, see Methods). Across eight patient samples with four different tumor types (lung, melanoma, breast and osteosarcoma), we generated over 700 in silico admixture samples. Detection and filtering were applied on each patient sample to benchmark the tumor fraction detection sensitivity and accuracy. (b) Patient-specific SNV signal-to-noise quantification over a range of TFs (10−5-0.2) compared to basal noise signal detected in control (TF = 0) samples (left column), estimated using melanoma sample (Pat.01). Signal-to-noise was estimated by calculating the log difference between the number of detections in each plasma sample (TF > 0) and the mean number of detections in the controls (TF = 0). Inset panel shows discrimination of tumor and control samples down to tumor fraction 10−5 after utilizing machine-learning-based sequencing error suppression (red) vs. reduced sensitivity with the raw unfiltered data (blue). For TFs > 0, n = 11 independent admixture samples. For the control (TF = 0), n = 20 independently down-sampled PBMC replicates. (c) Correlation between the number of base-pairs evaluated for SNVs and the number of mismatches detected (artefactual detections), measured over all synthetic control plasma (no tumor DNA, TF = 0, n = 341) from eight patients and mutational compendia from four tumor types (lung, breast, melanoma and osteosarcoma). Results show a constant error rate, independent of tumor type, corresponding to previously published, Illumina sequencing error-rate estimates (~1/1000 bps) (correlation and statistical significance calculated with a two-sided Pearson correlation). (d-e) Tumor fraction (TF) inference using genome-wide SNV integration for lung cancer (Pat.03) (d) and melanoma (Pat.01) (e) samples, shows accurate TF estimation as low as 5*10−5 and 10−5, respectively, discriminated from control (TF = 0) samples (left box-plot). High Pearson correlation (two-sided test) between the input TF mixture (x-axis) and the SNV-based estimated TF prediction, confirms accurate inference based on genome-wide mutational integration. For TFs > 0, n = 11 independent admixture samples. For the control (TF = 0), n = 20 independently down-sampled PBMC replicates. (f) The lower-limit-of-detection (LLOD) was empirically measured as a function of the input mutational load and the WGS coverage showing that for high mutational loads (~60,000 mutations) achieving sensitivity nearing 10−6 is feasible. Analysis was done over 21,420 in silico admixtures, varying TF (10−3-10−6), coverage (10–120X), mutation load (2,000–63,000), and across 20 replicates (random read downsampling) for each admixture. Lower-limit-of-detection (LLOD) was defined by the lowest tumor fraction that show significant Z-score separation from the control (TF = 0) cohort, for the same context (mutation load and coverage depth, see Methods). Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR.

Figure 3:. Patient-specific genome-wide CNA integration provides ultra-sensitive ctDNA detection and precision tumor fraction (TF) estimation.

(a) Illustration of the genome-wide copy number alteration (CNA) integration for the inference of plasma tumor fraction (TF). Patient-specific CNA segments are detected from matched tumor and germline WGS. Tumor CNA segments are then used to calculate the genome-wide accumulated signal in the patient’s plasma sample. Amplification segments (right panels, pink) in synthetic plasma with TF of 1% and 0.1% show sparse TF correlated signal that can be integrated due to the positive read coverage bias. Deletion segments (left panel, blue) show similar sparse signal, which is negatively biased. Upper panel: normalized coverage difference between plasma and control (TF = 0) across a 10Kbp region. Lower panel: sum of coverage differences between plasma and control. (b) Utilizing patient-specific CNA mixture model allows for the integrated CNA signal to be directly converted to TF estimates in a breast cancer (Pat.05) in silico synthetic plasma sample, with accurate TF estimation as low as 5*10−5, discriminated from the control samples (TF = 0, left box-plot). High Pearson correlation (two-sided test) between the input TF mixture (x-axis) and the CNA-based estimated TF (y-axis) shows accurate inference based on CNA patterns. N = 20 for each TF > 0 sample, as well as the control (TF = 0) samples. (c) Utilizing the CNA mixture model on copy-neutral regions of the same tumor genome (i.e., regions that do not show differential read-coverage between tumor and germline WGS), do not show correlation with TF (two-sided Pearson correlation). N = 20 for each TF > 0 sample, as well as the control (TF = 0) samples. (d) CNA method lower-limit-of-detection (LLOD) was empirically measured as a function of the CNA load showing high sensitivity (5X10−5) for tumors with copy number variation footprint of 1Gb or more.Upper panel: ctDNA detection at various TF (10−5 - 5*10−1 and control (TF=0)) and CNA load (5Mb - 1050Mb) estimated in a breast cancer (Pat.05) in silico synthetic plasma samples. Detection is observed at TF=5*10−3 at 5Mb and as low as TF=5*10−5 at 1Gb (Z score >= 1). Lower panel: Smoothed histogram (kernel density estimate) of a representative cohort of WGS breast cancer (BRCA) samples from TCGA (n = 95). Majority (97.9%) of the samples have CNA load above 10Mb and 33.7% have CNA load above 1Gb, corresponding with ultra-high ctDNA sensitivity. CNA method lower-limit-of-detection (LLOD) was empirically measured as a function of the CNA load showing high sensitivity (5X10−5) for tumors with copy number variation footprint of 1Gb or more. Upper panel: ctDNA detection at various TF (10−5 - 5*10−1 and control (TF=0); 20 random admixtures each) and CNA load (5Mb - 1050Mb; 3 random subsampling each) estimated in a breast cancer (Pat.05) in silico synthetic plasma samples. Throughout the figure, boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR.

Figure 4:. Detection of ctDNA using MRDetect in Melanoma during immunotherapy and colon cancer post-operatively.

(a) Error rate estimation in a cohort the test control plasma samples (n = 30) with and without error suppression. Applying support vector machine (SVM) error suppression and paired-end read concordance allow sequencing error reduction by a median 21 fold and increase the uniformity between samples (2-fold decrease in coefficient of variation). Boxplots represent median, bottom and upper quartile; whiskers correspond to 1.5 x IQR. (b) Melanoma treatment response during immunotherapy (Nivolumab) is monitored by blood samples. Treatment monitoring by computed tomography (CT) shows response to therapy but residual disease after 3 months of therapy (upper panel). MRDetect Z-scores effectively track tumor responses, matching radiographic changes, in higher temporal resolution than that feasible with imaging (middle panel). Ichor-CNA sensitivity captures initial treatment response dynamics but does not detect residual disease after 3 months of treatment. Z-score is calculated from a single plasma sample for each timepoint compared to a panel of control samples (n = 30). (c) Illustration of the colon cancer clinical cohort and potential clinical utility for liquid biopsy residual disease detection post-operatively. Upon surgical resection of tumor, the primary tumor is sequenced to define the tumor mutational compendium. Pre-operative plasma is used to verify tumor mutational integration detection in cfDNA and establish detection performance metrics. Post-operative plasma (median 6 weeks) is used to establish the presence of residual disease. (d) Colorectal cancer (CRC) detection on pre-operative plasma samples (n = 19 patients), using SNV (red) and CNA (gray) MRDetect methods. Z-score estimation (bar plots) was used to calculate mutation signature detection in the patient plasma in comparison to detection in the control plasma test cohort (n = 30). Z-score (> 4 and >1.3) was used for the SNV or CNA methods, respectively, to define a positive tumor detection based on ROC optimization on each method (e) Receiver operating curve (ROC) analysis on a combined detection model of SNV and CNA mutational compendia. Pre-operative plasma samples (n = 19) were used as the positive label, and the panel of control plasma samples against all patient mutational compendia (n = 570, nineteen mutational compendia assessed across thirty control samples) was used as the negative label. Two operational points (Z-score thresholds) are shown that correspond with two different specificity regimes. (f) Detection of ctDNA post-operative (6 weeks after surgery) in plasma samples (n = 19) was performed using Z-score estimation in comparison to the cohort of control plasma samples (n = 30). Patients with clinical recurrence, are indicated by red dots. Post-operative ctDNA detection (lower panel) is associated with early disease recurrence as shown in the pie chart. None of patients with the non-detected post-operative ctDNA (upper panel) showed clinical recurrence after median clinical follow-up of 15 (range 4–32) months. (g) Kaplan Meier disease free survival analysis was done over all patients with detected (n = 7) and non-detected (n = 12) post-operative ctDNA. Post-operative ctDNA detection shows association with shorter recurrence-free survival (RFS) (two-sided log-rank test).

Figure 5:. Detection of ctDNA using MRDetect in lung adenocarcinoma (LUAD) pre- and post-operative.

(a) Illustration of the LUAD clinical cohort and potential clinical utility for liquid biopsy residual disease detection post-operatively. Upon surgical resection of tumor, the primary tumor is sequenced to define the tumor mutational compendium. Pre-operative plasma is used to verify tumor mutational integration detection in cfDNA and establish detection performance metrics. Post-operative plasma (median 2.5 weeks) is used to establish the presence of residual disease. Patients with detectable ctDNA in their post-operative plasma samples show higher probability for disease recurrence even when there is no other clinical or imaging confirmation at that time. (b) Early stage lung adenocarcinoma (LUAD) detection on pre-operative plasma samples (n = 36 patients), using SNV (red) and CNA (gray) MRDetect methods. Z-score estimation (bar plots) was used to calculate the mutation signature detection in the patient plasma in comparison to detection in the control plasma test cohort (n = 30). Z-score was used for the SNV or CNA methods to define a positive tumor detection, indicated by black dots. (c) Early stage lung adenocarcinoma (LUAD) tumor fraction estimation per stage. All pre-operative plasma samples were considered for tumor fraction estimation by MRDetect, non-detected samples were assigned a 0 value. Bar plot represent the median TF value per stage, black line represents the confidence intervals (iqr) and the specific patient TF values are overlaid for the detected (blue dots) and non-detected (red dots) patients. Estimated tumor fraction shows association with cancer stage, with stage III samples showing the highest TF values. (d) Receiver operating curve (ROC) analysis on a combined detection model of SNV and CNA mutational compendia. Pre-operative plasma samples (n = 36) were used as the positive label, and the panel of control plasma samples against all patient mutational compendia (n = 1,080, thirty-six mutational compendia assessed across thirty control samples) was used as the negative label. Two operational points (Z-score thresholds) are shown that correspond with two different specificity regimes. (e) Detection of ctDNA post-operative (median 2.5 weeks after surgery) in plasma samples (n = 22) was performed using Z-score estimation in comparison to the cohort of control plasma test cohort (n = 30). Z-score was used for the SNV or CNA methods to define a positive tumor detection, indicated by black dots. Patients with clinical recurrence are indicated by red dots. Post-operative ctDNA detection (lower panel) is associated with early disease recurrence, as shown in the pie chart. None of patients with the non-detected ctDNA (upper panel) showed clinical recurrence after median clinical follow-up of 18 (range 6–36) months. (f) Kaplan Meier disease free survival analysis was done over all patients with detected (n = 10) and non-detected (n = 12) post-operative ctDNA. Post-operative ctDNA detection shows association with shorter recurrence-free survival (RFS) (two-sided log-rank test). (g) Fragments containing tumor mutations show fragment size shift in comparison to non-tumor (cross-patient) mutations in the same patient plasma. The fragment size shift was estimated on detected post-operative patients and non-detected post-operative patients using a kernel-density-estimation method (see Methods). Patients that showed ctDNA detection post-operative also exhibit significant shift (two-sample t-test) in their tumor-specific mutation detections in the patient plasma (light red, n = 226), compared to non-tumor (cross-patient) detected mutations on the same samples (light blue, n = 1603). Patients that do not show ctDNA detection post-operative also do not exhibit significant fragment size shift in tumor-specific mutation detections (dark red, n = 160) compared to non-tumor (cross-patient) detected mutations on the same samples (dark blue, n = 2121). Violin plots depict kernel density estimates of the density distribution. Center dashed lines represent the median and dashed lines represent the interquartile range.