Comparison of variant allele frequency and number of mutant molecules as units of measurement for circulating tumor DNA

Here, we compare the number of mutant molecules and the variant allele frequency as units of measurement for circulating tumor DNA. We conclude that when using digital droplet PCR, both units of measurement yield greater agreement than when using next generation sequencing. Low frequent mutations and molecular coverage are important factors affecting agreement. These results will improve circulating tumor DNA analysis, contributing to better cancer patient management.

(Received 23 July 2020, revised 10 September 2020, accepted 15 October 2020, available online 31 October 2020) doi:10.1002/1878-0261.12827 [Correction added on 21 December 2020, after first online publication: Peer review history is not available for this article, so the peer review history statement has been removed.] Quantification of tumor-specific variants (TSVs) in cell-free DNA is rapidly evolving as a prognostic and predictive tool in patients with cancer. Currently, both variant allele frequency (VAF) and number of mutant molecules per mL plasma are used as units of measurement to report those TSVs. However, it is unknown to what extent both units of measurement agree and what are the factors underlying an existing disagreement. To study the agreement between VAF and mutant molecules in current clinical studies, we analyzed 1116 TSVs from 338 patients identified with next-generation sequencing (NGS) or digital droplet PCR (ddPCR). On different study cohorts, a Deming regression analysis was performed and its 95% prediction interval was used as surrogate for the limits of agreement between VAF and number of mutant molecules per mL and to identify outliers. VAF and number of mutant molecules per mL plasma yielded greater agreement when using ddPCR than NGS. In case of discordance between VAF and number of mutant molecules per mL, insufficient molecular coverage in NGS and high cell-free DNA concentration were the main responsible factors. We propose several optimization steps needed to bring monitoring of TSVs in cell-free DNA to its full potential.

Introduction
The genomic characteristics of solid tumors increasingly determine how patients are being treated. Although metastatic tissue can be obtained for this analysis, it is a cumbersome procedure and thereby limits repetitive sampling. Molecular profiling of cell-free DNA (cfDNA) in liquid biopsies from patients with cancer is evolving rapidly as a patient-friendly tool to measure tumor load as well as to gain insight into tumor characteristics [1].
Although isolation of cfDNA from plasma is an easy procedure, DNA fragments from nonmalignant cells (e.g., lymphoid and myeloid cells) hamper the subsequent detection of tumor-specific variants (TSVs) [2]. Generally, tumor-derived cfDNA fragments (circulating tumor DNA or ctDNA) represent a minority of all cfDNA fragments present in plasma [3]. The lower limit of detection (LOD) of TSVs has been improved by digital droplet PCR (ddPCR), enabling the detection of tumor-specific variants in a single DNA molecule [4]. In addition, the LOD of larger gene panels used for next-generation sequencing (NGS) has been optimized by the development of unique molecular identifiers (UMIs). UMIs are added to each molecule before amplification, allowing correction of sequencing errors and identification of individual mutated templates [5,6].
The quantity of a tumor-specific variant is typically reported as the ratio between the number of mutatedand wild-type DNA copies and is referred to as the variant allele frequency (VAF). During its clinical implementation, it has become clear that solely reporting the VAF does not suffice as it does not provide information on the concentration of a TSV. Use of a measurement that reports a concentration is common practice for other biomarkers, such as tumor antigens, as it is considered to reflect tumor load more adequately. Moreover, the cfDNA concentration is known to yield important prognostic value [7][8][9]. To this end, it might be preferable to use mutant molecules per mL plasma as a unit of measurement for monitoring of TSVs.
To what extent mutant copies per mL plasma relate to the VAF is currently unknown. To optimize characterization of cfDNA, including monitoring of TSVs in blood over time, information on agreement between VAF and mutant copies per mL plasma obtained by using real-life data is pivotal. Here, we report the agreement between VAF and mutant copies per mL plasma as units of measurement for 1116 TSVs quantified by NGS or ddPCR using current-day pipelines. Secondly, we identify pre-analytical and analytical factors that hamper agreement between both units of measurement.

Sample collection
Samples were obtained from the following studies: IMPACT-CRC study on colon cancer (ClinicalTrials.gov, number NCT02117466) [10], REGORA study on colon cancer (ClinicalTrials.gov, number NCT02800330), START-TKI study on lung cancer (CCMO number, NL58664.078.16) [11], TAX-ESR1 study on breast cancer (trialregister.nl, number NL7280), and the CareMore-Trastuzumab study on breast cancer (trialregister.nl, number NL4977). The study was performed in accordance with the Declaration of Helsinki and approved by the medical ethics committee of the Erasmus MC. All patients gave written informed consent prior to study procedures.
Samples from patients with lung cancer were collected upon progression on current therapy for detection of primary activating and p.T790M EGFR mutations. Samples from patients with colon cancer were collected before start and during anti-EGFR therapy or before start and during regorafenib therapy. Samples from patients with breast cancer were collected before start of first-line taxane-based chemotherapy. For some patients with colon cancer, serially collected samples were available. Blood was collected in either EDTA or CellSave tubes. For blood collected in EDTA tubes, plasma was separated within 4 h as recommended by the latest guidelines [12,13].

Description of included cohorts: cfDNA isolation and quantification
The pre-analytical work-up of samples within a study cohort were consistent and are described below.

Cohort 1
Patients with metastatic colon cancer. CfDNA was isolated using the QiaSymphony Circulating DNA kit (Qiagen, Venlo, the Netherlands) from 950 to 4000 µL plasma and eluted in 70 µL elution buffer. Depending on the amount of cfDNA, 20 ng was used for NGS. Samples were concentrated using a speedvac concentrator if necessary. For a subset of patients, longitudinal ddPCR data were available. Generally, 7 µL eluate or less was used for ddPCR, depending on the eluate concentration. The maximum input was 50 ng, and the median input was 12 ng.

Cohort 2
Patients with lung cancer. CfDNA was isolated from 3 mL plasma with the QIAGEN QIAamp Circulating Nucleic Acid kit (Qiagen) and eluted in 50 µL elution buffer. For NGS, the samples were then concentrated to 25 µL using a Speedvac concentrator before cfDNA was quantified by Qubit fluorometer (Invitrogen, Carlsbad, CA, USA). Subsequently, 13 µL of concentrated sample with a maximum of 50 ng was used for NGS. For a subset of patients, ddPCR data were available. For ddPCR, 4 µL of the unconcentrated eluate was used regardless of cfDNA concentration.

Cohort 3
Patients with metastatic breast cancer. CfDNA was isolated from 600 to 4000 µL plasma using the QIA-GEN QIAamp Circulating Nucleic Acid kit (Qiagen) and eluted in 50 µL elution buffer. Depending on the amount of cfDNA 10 ng was used for NGS. Samples were concentrated using a speedvac concentrator if necessary. No ddPCR data were available for this cohort.
For all samples across cohorts, cfDNA concentrations were measured using the Quant-iT dsDNA highsensitivity assay (Invitrogen, Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions, and the Qubit fluorometer (Invitrogen) was used as read out.

Next-generation sequencing
All samples were sequenced with the Ion Torrent TM Oncomine TM cfDNA Assay for breast, colon or lung cancer for the respective cancer type, on the Ion Torrent S5XL Next Generation Sequencing (NGS) prime system, all according to protocols and consumables provided by the manufacturer [10,11] (Life Technologies, Thermo Fisher Scientific, Carlsbad, CA, USA). NGS panels were equipped with UMIs to enable detection of unique mutated copies. Molecular coverage was defined as the count of unique molecules and known hotspot variants were analyzed if they were detected in at least three independent molecules. Additionally, variants detected in cohort 1 and cohort 3 were called as true variants either if the molecular coverage was sufficient given the input of cfDNA (a minimum of 500 unique molecules for each 10 ng of DNA sequenced), or if the variant was detected with sufficient coverage in an earlier sample from the same patient.

ddPCR analysis
Analysis of mutations was performed using uniplex ddPCR mutation assays from Bio-Rad Laboratories or Thermo Fisher as previously described [11,14]. Variants were designated true variants if they were detected in at least three independent molecules. The number of droplets positive for mutant or wild-type molecules was fitted into a Poisson distribution to determine the absolute number of copies per µL eluate, thereby correcting for droplets containing more than one molecule. The number of mutant molecules was then derived from the concentration of mutant copies per µL eluate. The VAF was reported by the Bio-Rad software, after correction for the Poisson distribution. Molecular coverage was calculated as the sum of wildtype and mutant copies.

Data collection and definitions
Data on sequencing input, isolation protocol, Qubit measurements, amount of plasma used for isolation, molecular coverage, and number of mutant copies were collected. VAF was calculated as (number of mutant copies/(number of wild-type copies + number of mutant copies)) 9 100% and the number of mutant copies per mL plasma was calculated follows, for NGS: (number of mutant copies/DNA input for sequencing (ng)) 9 (cfDNA concentration (ngÁmL À1 plasma)), and for ddPCR: (number of mutant copies/ input for analysis (lL)) 9 (total eluate (lL)/amount of plasma used for isolation (mL)). As a result, the multiplication factor by which the number of mutant copies are multiplied to calculate the number of mutant molecules per mL plasma was defined as the ratio between the cfDNA concentration and the cfDNA assay input: cfDNA concentration (ngÁmL À1 plasma)/cfDNA sequencing input (ng) (NGS) or Total eluate (lL)/ (Assay input (lL) 9 Total amount of plasma used for isolation (mL)) (ddPCR).

Statistics
To assess the agreement between VAF and mutant molecules per mL plasma, a Deming regression analysis was performed and its 95% prediction interval (PI) was calculated. The average width of the 95% prediction interval was used as a surrogate for the 95% limits of agreement. This method is adopted from Bland and Altman [15], which is used to study agreement between measurements with different units. The average width of the 95% PI is an important measure that indicates whether it is clinically acceptable to replace one method by another. TSVs outside the 95% PI were considered outliers, in which TSVs that lie above the 95% PI upper limit were considered as upper limit outliers, whereas TSVs that lie below the 95% lower limit were considered lower limit outliers, for x = VAF and y = mutant copies per mL plasma.
To assess the association between cfDNA concentration and molecular coverage for both VAF and mutant copies per mL, we calculated Pearson correlation using linear regression.
For NGS, we analyzed the within-sample read coverage and molecular coverage of the different hotspots among the breast panel in cohort 3. Those samples were sequenced in different sequencing runs. The coefficient of variation for both the read coverage and the molecular coverage of different hotspots within a unique sample was reported. The coverage was normalized for cfDNA input using the total number of reads within a sample. Additionally, we calculated the positive and negative predictive value (PPV and NPV, respectively) of the factor molecular coverage < 500 molecules for lower limit outliers.
To analyze the correlation between pre-analytical and analytical factors, Pearson correlation was calculated. All variables were log transformed before analysis.
Descriptive statistics were performed using IBM SPSS STATISTICS 25 (IBM Corp, Armonk, NY, USA). All regression computations and graphics were performed in R program language [16]. Because the different cohorts included in this analyses were heterogeneous with regard to the pre-analytical work-up, we performed statistical analyses on each cohort separately. Thereby, our estimates of agreement approximated the real-life setting.
Outliers were identified in samples with both single and multiple TSVs. In some cases, all TSVs within a sample were outliers, whereas in other cases, samples contained TSVs that were either outlier or showed agreement (Table 2A,B, Table S2).

Distribution of pre-analytical and analytical variables among outliers
For NGS, the cfDNA concentration and molecular coverage were higher in upper limit outliers compared to TSVs showing agreement (cfDNA concentration median, cohort 1: 218 vs 54 ngÁmL À1 ; cohort 2: 192 vs 14 ngÁmL À1 ; cohort 3: 58 vs 20 ngÁmL À1 ) (molecular coverage median, cohort 1: 8512 vs 3145X; cohort 2: 9202 vs 2090X; cohort 3: 2842 vs 1338X) (Table 2A). Only in samples from cohort 2, the nanogram cfDNA input was higher among upper limit outliers (47 vs 20 ng). However, cfDNA input strongly correlated with cfDNA concentration and molecular coverage in this cohort since a fixed eluate volume was used for sequencing (Pearson's r: 0.600 and 0.606) (Table S1A). Additionally, both molecular coverage and cfDNA concentration were lower among lower limit outliers than in TSVs with agreement for all cohorts (molecular coverage median; cohort 1: 424 vs 3145X; cohort 2: 203 vs 2090X; cohort 3: 269-338 vs 1338X) (cfDNA concentration median; cohort 1: 19 vs 54 ngÁmL À1 ; cohort 2: 4 vs 14 ngÁmL À1 ; cohort 3: 8 vs 20 ngÁmL À1 ). Although the sequencing input from samples in cohort 1 did not substantially differ between upper or lower limit outliers and samples that showed agreement, both cfDNA concentrations and molecular coverage were lower in lower limit outliers and higher in upper limit outliers compared to TSVs with agreement. For ddPCR, no lower limit outliers were detected (Table 2B). Nanogram input, molecular coverage, and cfDNA concentration were higher in upper limit outliers in all cohorts. However, ddPCR all these variables were however strongly correlated with each other (Table S1B).

Molecular coverage in NGS and its impact on agreement
In lower limit outliers, the molecular coverage was substantially lower compared to samples showing agreement. A molecular coverage < 500 molecules, set arbitrary, resulted in a PPV of 0.55-1 and a NPV of 0.99-1 for lower limit outliers among all NGS detected TSVs. The median coefficient of variation (CV) of the read coverage and molecular coverage among different amplicons within a sample were 21% and 18% for the breast panel, respectively. The maximum observed CV within a sample was 36% for the read coverage and 44% for the molecular coverage. The coverage of gene positions located on the same amplicon did not vary. However, the coverage of amplicons within genes and across genes was highly divergent (Fig. 2).

Discussion
The use of ctDNA as a biomarker for real-time monitoring of disease burden in a minimally invasive way is a promising tool to evaluate TSVs in plasma of patients with solid tumors. Although the importance of ctDNA load in plasma has been recognized before, its quantification is still in its infancy as demonstrated by the different units of measurements that are currently used to report those TSVs. In our analyses, we assessed agreement between VAF and mutant copies per mL and determined factors that affect this agreement. Here, we propose several optimization steps that result from these analyses.
Firstly, our analyses demonstrate that a low molecular coverage resulted in a severe underestimation of the absolute number of mutant molecules or an overestimation of the VAF. Assuming the sequencing efficiency for both mutant and wild-type copies is equal, insufficient molecular coverage would affect the VAF to a lesser extent. A molecular coverage < 500 molecules in NGS had a high NPV for lower limit outliers indicating that a molecular coverage of > 500 In addition, we demonstrate that the molecular coverage in NGS is highly variable among amplicons present in a panel. For some individual samples, the within-run CV for different amplicons was as high as 44%. For different variants within a sample with similar VAFs this would result in nearly a doubling of the number of mutant molecules per mL, solely based on the molecular coverage of certain positions. As the read coverage was impacted by a very similar variation, this observation is mainly attributable to variability in sequencing efficiency and less likely from UMI adapter ligation and subsequent read loss and Table 2. Distributions of (pre-)analytical factors among TSVs with agreement and outliers identified by (A) NGS and (B) ddPCR. Results are presented as median (range) unless indicated otherwise. When ≤ 2 outliers are present, only a range is given. Multiplication factor: cfDNA concentration (ngÁmL À1 )/Sequencing input (ng), TSV, tumor-specific variant; NP, not present, no lower limit outliers were present in ddPCR.

A Variable
Agreement Lower limit outlier Upper limit outlier  [17,18]. Lower limit outliers resulting from insufficient molecular coverage were identified in both samples with single and multiple TSVs. In some samples, only one amplicon was affected by low sequencing efficiency, whereas in other samples the molecular coverage was below 500 molecules for all amplicons. Amplicons in the genes TP53 and EGFR were most frequently affected by sequencing efficiency (Table S2). In this study, we have used sequencing data generated by the use of Oncomine amplicon-based panels designed for the IonTorrent sequencer as these panels are commonly used in the diagnostic facility of our Pathology department. For both amplicon and capture-based panels, the molecular coverage is subjected to variation. Capture-based panels can isolate neighboring regions that are not of interest, which will result in lower overall coverage in the regions of interest. Contrarily, amplicon-based panels are subjected to issues related to primer design. Single nucleotide polymorphisms and short indels in the primer template region might cause allelic dropout resulting in decreased molecular coverage, whereas amplification of genes with high guanine-cytosine content will less likely be effective [19]. To overcome this problem, NGS panels should be equipped with adequate quality controls to enable correction for molecular coverage per amplicon after sequencing [20]. Another factor affecting agreement between VAF and mutant molecules per mL plasma was the cfDNA concentration. A high cfDNA concentration was associated with upper limit outliers. However, not all TSVs in samples with a high cfDNA concentration were upper limit outliers, they were mainly TSVs with a frequency of < 0.01% (Table S2). Low frequent TSVs (i.e. TSVs detected in a low absolute number of mutant copies) are more frequently prone to stochastic errors when calculating the number of mutant molecules that originate from 1 mL plasma. Alternatively, stochastic errors that occur during amplification favored by the abundance of wild-type molecules might result in an underestimation of the VAF.
In our analysis, quantification of TSVs by ddPCR yielded less outliers and smaller limits of agreement for mutant copies per mL. In addition, the molecular coverage was always sufficient to estimate mutant molecules per mL plasma resulting in the absence of lower limit outliers. Our analysis does not substantiate whether a fixed input based on eluate volume is superior to input based on cfDNA amount. Although not directly substantiated by our data, it is well known that the molecular coverage is an important determinant of the limit of detection of a TSV [19,20]. To avoid false-negative results, the molecular coverage should be sufficient to detect a TSV of a given frequency. To this end, the molecular coverage is an important parameter and should be reported for proper interpretation of sequencing results.
Although our analyses were limited by the inability to analyze the accuracy of all pre-analytical and analytical steps and their impact on agreement separately, it does provide insights into the agreement of both units of measurement in the current real-life setting ( Table 3). Although our primary aim was to investigate agreement between VAF and mutant molecules for different sequencing platforms, we lacked power to investigate differences between pre-analytical and analytical factors among outliers. We therefore only used descriptive statistics to describe differences among preanalytical and analytical variables between outliers and TSVs that showed agreement. For molecular coverage however, we were able to identify a threshold that showed a high NPV. Additionally, we did not include DNA from leukocytes to exclude germline variants or variants resulting from clonal hematopoiesis [2]. Although clonal hematopoiesis is known to occur in specific genes and we reported predominantly activating, cancer-specific hotspot mutations, a nonmalignant origin of reported TSVs cannot be excluded. All samples from cohorts that were included in this study were collected for research purposes and pre-analytical and analytical methods were homogeneous within study cohorts. The methods used in all cohorts are in agreement with current guidelines and recent literature. To this end, results from our analyses reflect the accuracy of the real-life analyses pipeline as a whole.

Conclusions
From our study, we conclude that VAF and number of mutant molecules per mL plasma yielded greater agreement when using ddPCR than when using NGS. Given the higher costs of NGS, ddPCR might be preferred in case a specific variant can be tracked. For clinical purposes, the unit of measurement should reflect the concentration of cfDNA, alike other biomarkers. However, we demonstrate that quantification of absolute numbers of molecules per mL plasma is currently more heavily affected by pre-analytical and analytical variables than the VAF. To this end, standardization of methods is necessary. Based on these analyses, we propose several optimization steps to be taken with respect to the isolation and quantification of cfDNA but also for NGS panels, to bring longitudinal monitoring of TSVs in cfDNA to its full potential.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1. (A) Correlation among (pre)analytical factors in NGS. (B) Correlation among (pre)analytical factors in ddPCR  Table S2. Pre-analytical variables in samples containing outliers