Defining the Comprehensive Genomic Landscapes of Pancreatic Ductal Adenocarcinoma Using Real-World Endoscopic Aspiration Samples

Pancreatic ductal adenocarcinoma (PDAC) is a disease with a dire prognosis and one of the few cancers with a rising incidence, leading to estimates of it becoming the second leading cause of cancer-related death in the United States within the next decade (1). A series of studies on large-scale molecular characterization of PDAC have elucidated the comprehensive genomic landscape of this neoplasm, including the so-called “long tail” of potentially actionable mutations (many at as low prevalence rates as ∼1%), which may form the basis for precision oncology and clinical trial inclusion (2–6). Most of these prior “discovery” studies have used archival surgically resected specimens with optimal tumor cellularity (e.g., >60% for The Cancer Genome Atlas cohort), limiting much of our current knowledge to localized, and therefore resectable tumors (3, 4). Nonetheless, the majority of patients with PDAC (∼80%–85%) are diagnosed with locally advanced or metastatic disease, precluding surgical options, which, in turn, restricts many patients from the opportunity for in-depth genomic analyses from resected specimens (1, 7).

All patients with PDAC, irrespective of disease stage, are required to have histologic or cytological confirmation of their underlying diagnosis prior to onset of therapy (7). The two most common avenues for obtaining diagnostic biospecimens in advanced disease include either endoscopic ultrasound–guided fine-needle aspiration (EUS-FNA) of the primary tumor, performed by a gastroenterologist, or a percutaneous core biopsy of a biopsy-amenable metastatic lesion, when possible, the latter performed by an interventional radiologist. Overall, with the gastroenterology clinic being the typical “portal of entry” for symptomatic individuals, an EUS-FNA represents one of the most common, if not the most common, sample type obtained in patients with PDAC at the time of de novo presentation. Although both sources of biopsies are comprised of minimal tissue material, the challenges of using EUS-FNA from PDAC for next-generation sequencing (NGS) studies are unique. In contrast to the generally high cellularity of metastatic samples obtained by core biopsy, EUS-FNA samples of desmoplastic primary PDAC tend to be hypocellular, and potentially contaminated with nonneoplastic gastric and duodenal mucosa, replete with mucosal inflammatory cells, which can “drown” cancer-associated alterations (8). Further, the limited neoplastic DNA yield from EUS-FNA samples often precludes the use of even the currently available targeted NGS panels, let alone more comprehensive whole-exome sequencing (WES). Consequently, most published studies have been limited to assessing a few common “hotspot” mutations in PDAC, such as KRAS or BRAF (9). At the same time, it is worth noting that the currently available targeted NGS panels are geared toward common actionable genes, and not included are the aforementioned low-frequency alterations (so-called “long tail”). Third, the most commonly used first-line, and by some measures, most efficacious therapy in advanced PDAC is FOLFIRINOX, a platinum-containing multidrug regimen (10). Although preclinical data support that patients with PDAC with homologous recombination repair defects (HRD) are likely to bear greatest susceptibility to platinum agents like oxaliplatin (an integral component of FOLFIRINOX), the targeted NGS panels only interrogate a fraction of HRD, such as the approximately 5% of patients with deleterious BRCA1/2 mutations. Yet, prior studies have suggested that the full compendium of HRD in PDAC, as reflected in the “unstable genome” phenotype comprised of multiple structural abnormalities, might be several fold higher in prevalence (3). As a result, most first-line therapy decisions in the advanced disease setting are made empirically, warranting strategies that can provide improved treatment prediction information from the PDAC genome.

In this study, we use real-world limited biospecimens obtained via EUS-FNA, almost all in the diagnostic “first encounter” that a patient with PDAC has with their gastroenterologist, for in-depth genomic analysis of PDAC. Applying a relatively facile enrichment strategy, we were able to significantly increase neoplastic cell fraction, which allowed us to detect actionable mutations using molecularly barcoded WES for clinical decision-making. Further, we demonstrate the ability to use the WES information paired with bioinformatics pipelines to determine potential HRD status in the cases, which can help guide the choice of first-line regimen in a clinically feasible timeline. Our data establish the feasibility for WES in the most commonly obtained diagnostic biospecimen in PDAC, EUS-FNA samples, with consistent results compared with previous studies.

Investigators obtained written informed consent under MD Anderson protocol Lab00–396 from each patient prior to tissue and blood sampling. The study was performed in accordance with standard ethical guidelines approved by the Institutional Review Board and in accordance with the Declaration of Helsinki. All patients had clinically and histologically confirmed localized or metastatic PDAC. In the pilot phase, resected tissue specimens were only included if tumor cellularity was below 10% (Supplementary Fig. S1B). All biopsies (CT-guided core or EUS-FNA) were taken within routine clinical procedures. No EUS-FNB (fine-needle biopsy) samples have been included. All research passes harvested were sampled after routine diagnostic EUS-FNA passes and ROSE-rapid on-site cytology assessment which showed adequate cellularity in diagnostic specimens.

No complications have been reported in any tissue sampling.

Biopsy samples were mechanically and enzymatically digested into single-cell suspension as described before (11). The cell suspension was then processed using the EasySepTM Human EpCAM Positive Selection Kit (Stemcell, cat#18356) following the manufacturer's protocol to enrich for epithelial cell adhesion molecule (EpCAM)-positive epithelial cells. In brief, cell suspensions were incubated with EpCAM antibodies and magnetic beads, followed by subsequent clean-up steps using magnetic separation (Supplementary Fig. S1C).

Bulk genomic DNA was extracted from enriched and nonenriched cell suspensions using the QIamp DNA Micro Kit (Qiagen, cat#56304) and DNeasy Blood & Tissue Kit for peripheral blood mononuclear cells (PBMC; Qiagen, cat#69506) following the manufacturer's protocol. At least two different quantification methods for DNA quality and quantity were performed in parallel using a QubitTM dsDNA BR Assay Kit (Thermo Fisher, cat#Q32853), a NanoDropTM 2000/2000c spectrophotometers (Thermo Fisher, cat#ND2000), and/or high sensitivity D1000 screentape (HSD1000) with a Tapestation 2200 system (Agilent, cat#5067–5584). Digital droplet PCR (ddPCR) was performed for validation of KRAS mutations and MYC amplification as described in the Supplementary Materials.

A median of 160 ng of enriched tumor DNA or 200 ng of matched PBMC DNA was diluted in a total volume of 52 μL low TE buffer and fragmented to 150 to 200 bq using a Covaris LE 220 ultrasonicators system (Covaris). The following optimized settings were used: Duty Factor 30%, Peak Incident Power (PIP) W 450, cycles per burst 200, time 300 seconds, temperature 7°C, and water level 6. Adequate fragmentation was documented using HSD1000. Molecular-barcoded libraries were constructed following the SureSelect XT HT–targeted enrichment protocol (Version A1, July 2017). In brief, end repair and dA-tailing were followed by ligating individual molecular barcoding to each strand, PCR amplification, and bead-based cleaning. Libraries were then hybridized and incubated with a whole-exome capture library (SureSelect Clinical Research Exome V2; Agilent, cat#5190–9492), captured to streptavidin-coated beats, washed, and amplified. Libraries were multiplexed, denatured, and diluted to a final concentration of 1.8 pmol/L for sequencing and cluster generation as per the manufacturer's recommendation. Clustered flow cells were sequenced on the Illumina NextSeq 500 instrument (Illumina) using standard Illumina paired-end primers and chemistry. A median of 74,991,096 reads/PBMC (range, 25,067,034–251,278,724 reads) and 144,308,769 reads/enriched tumor (range, 25,451,842–236,573,016 reads) was detected. Median on-target coverage reached 56x (range, 23–157x) for matched PBMC and 108x (range, 14–177x) for enriched tumor samples for family size of 1 library.

Sequencing data were processed as detailed in the Supplementary Materials. Three algorithms were used for single-nucleotide variant (SNV) calling, and these were filtered as detailed in the Supplementary Materials. Of note, the median mutational burden of non–PBMC-paired samples (n = 5) was higher than in paired samples but did not reach statistical significance (4.89 mut/Mb; range, 0.72–23.91 vs. 1.63 mut/Mb; range, 0.29–21.46; P = 0.06) and seemed to had no impact on our ability to identify patients with increased tumor mutational burden (TMB; Supplementary Fig. S3C).

When available, Clinical Laboratory Improvement Amendments (CLIA) reports were used for SNV validation. These reports were generated from resected or core biopsied tissue samples during the course of treatment. Reports were available for 1 core biopsy patient (pilot phase) and for 4 FNA patients (second phase). Four of the CLIA reports were generated in-house at MD Anderson Cancer Center (MDACC) with a PCR-based sequencing platform and assembled using GRCh37/hg19 builds. MDACC CLIA reports included a panel of 134 individual genes (Solid Tumor Genomic Assay V1 report). The fifth patient report was performed by an outside CLIA-certified laboratory using a company-specific gene panel (Perthera).

Dana-Farber Cancer Institute (DFCI) cohort data were processed using the same bioinformatics pipeline used in our study with the exception of unique molecular identifiers (UMI)-specific processing.

We evaluated each sample for presence of 24 potentially targetable genes/alterations and SCNAs: BRCA1, BRCA2, PALP2, FGFR1, FGFR2, FGFR3, FGFR4, PDGFR, c-Kit/CD117, ROS1 (fusion), MET, NOTCH1, JAK1, JAK2, JAK3, mTOR, BRAF, RNF43, PI3K/PIK3CA, AKT, NTRK, MYC, KRAS wild-type, and HER2 amplification.

TMB was defined by the sum of all synonymous SNV, nonsynonymous SNV, stopgain, stoploss, frameshifts, or indels per sample that passed filter criteria as described above. For PBMC-matched samples, germline mutations called by HaplotypeCaller were excluded but for mutations in important PDAC genes like, e.g., BRCA1/2, PALB2, or ATM (12). Non–PBMC-matched samples are indicated in the analysis workflow and throughout results. TMB was then calculated by dividing the total number of mutations/patient by the coding region target of the SureSelect Clinical Research Exome V2 panel used (67.3 Mb).

Mutational signature analysis of FNAs was performed using MutationalPatterns following the steps outlined for cancer/Catalog of Somatic Mutations in Cancer (COSMIC) signature analysis (13). Only FNA samples with a matched germline control (PBMC sample) and those from treatment-naïve tumors were included in the mutational analysis.

Two algorithms were used for detection and classification of somatic copy-number alterations [SCNA; amplifications, deletions, and copy-neutral LOH (cnLOH)]. HapLOHseq was used for detection of genomic regions exhibiting allelic imbalance (AI). Results from this algorithm were combined with output from standard log₂ copy ratio segmentation data as described in the Supplementary Materials.

This approach detects B-allele frequency (BAF) shifts at germline heterozygous sites, indicative of AI, allowing for detection of chromosomal alterations in low mutant cell fraction settings (14, 15). Germline heterozygous sites with a depth greater than or equal to 10 were included as input for AI analysis. The hidden Markov model of hapLOHseq was used to compute the probability that a set of adjacent markers span a region of AI. An AI event was defined as a continuous set of markers with posterior probabilities exceeding the threshold of 0.85.

Genome Analysis Toolkit (GATK; ref. 16) was used for segmentation of log₂ copy ratio data. SCNAs were called by overlaying HapLOHseq AI and GATK segmentation calls (Supplementary Fig. S5A). For SCNA calls exclusive to GATK, only those with a log₂ copy ratio below −1 and above 0.58 for deletion/amplifications were included in the final call set. This approach allows for the inclusion of balanced amplifications/losses, which do not result in AI, as well as focal SCNAs spanning a modest number of germline heterozygous sites. Genomic regions with posterior probability of AI > 0.85 were included in the final SCNA call set and were classified using GATK segmentation values, and log₂ copy ratios below −0.41 and above 0.32 were called as deleted and amplified, respectively. cnLOH was defined as genomic regions with a log₂ copy ratio between ±0.15 and a BAF deviation >0.1. Events with a posterior probability of >0.85 and log₂ copy ratio between ±0.15 and −0.41 or +0.32 were defined as undetermined (Supplementary Fig. S4B). An aneuploidy score was calculated for each sample. It was defined as the number of chromosome arms (out of 39) with arm-level aneuploidy, following principles previously used in the analysis of genomic instability (17). A chromosome arm (except for the short arm of acrocentric chromosomes 13, 14, 15, 21, and 22) was categorized as aneuploid when an amplification, deletion, cnLOH, or undeterminable SCNA(s) spanned more than 75% of the chromosome arm. To accommodate for high-confidence AI events, with subtle log₂ copy ratio deviations, an arm-level SCNA was also called in cases where one type of SCNA (amplification, deletion, or cnLOH) plus undeterminable SCNA(s) spanned more than 75% of the arm, but a single type of SCNA did not reach the 75% threshold.

For identification of SCNA patterns suggestive of DNA damage repair deficiencies, we calculated a new and easy to compute score. This complexity score (CS) was calculated by summing the number of chromosome arms with two SCNAs with opposing classification (i.e., gain and loss, or gain cnLOH). Undetermined events and those smaller than 1 Mb were excluded. This chromosome arm–level signature is inferred to have arisen from double-strand breaks, thus higher scores are present in patients with deficiencies in DNA repair.

Progression-free survival (PFS) was defined as time period from the start of any treatment (chemotherapy, radiotherapy, or surgery) to date of disease progression (defined by RECIST 1.1 guidelines; ref. 18) or death. Overall survival (OS) was defined as time from tissue sampling (surgery, core biopsy, or EUS-FNA) to death for any reason. Survival curves were estimated using the Kaplan–Meier method. Statistical analyses were performed with SPSS statistical software, version 24 (IBM) or Prism 8 (GraphPad Software, Inc.). All tests were two-sided. Statistical significance was defined as a P value of <0.05. We used bootstrapping to estimate the sampling distribution of the CS score in the DFCI cohort.

To analyze the feasibility and impact of epithelial cell enrichment, we utilized samples from either surgically resected primary PDAC (N = 5) or CT-guided metastatic core biopsies (N = 4, 3 patients) in a pilot “first phase” of this study. The 9 pilot phase samples were obtained from 8 independent patients. Subsequently, for the main (“second phase”) study, we assessed 23 cytology samples from 23 patients with PDAC undergoing endoscopic FNAs between April 2017 and November 2018 at MDACC. Patient details and processing workflow for both phases are summarized in Supplementary Table S1 and Supplementary Fig. S1A. Matched blood samples were obtained in 24 of 31 patients (77.4%).

A major obstacle for genomic characterization of small biopsies, especially using WES, is the limited amount of starting material and the generally low tumor cellularity of pancreatic cancer. To overcome this known hurdle, we used a magnetic EpCAM-pulldown approach to enrich for epithelial tumor cells. In the pilot phase, we used five fresh resected tissues samples with a low cellularity (<10%) and four CT- or ultrasound-guided core biopsies (liver and lung metastases, N = 2 each; Supplementary Fig. S1A and S1B). All samples were mechanically and enzymatically dissociated into a single-cell suspension before further processing. EpCAM enrichment was performed using an EpCAM-based magnetic pulldown followed by DNA isolation (Supplementary Fig. S1C). ddPCR showed a significant increase in KRAS mutant allele frequency (MAF; median, 4.3% vs. 21.7%; P = 0.049; Fig. 1A and B). In 4 patients with a nonenriched KRAS MAF falling below the cutoff for inclusion in CLIA reports (<5%), enrichment increased the MAF above this threshold (Fig. 1A and B). In addition, we compared performance for WES SNV and SCNA calls between EpCAM-enriched samples from resected PDAC and matched archived formalin-fixed paraffin-embedded (FFPE) slides of the matched tumor at resection (n = 2). Although this comparison showed only a few alterations in FFPE tissues, a greater number of genomic alterations were observed in the enriched samples, underscoring the value of EpCAM enrichment (Supplementary Table S2; Supplementary Fig. S2A and S2B). For example, in patient 13, sequencing of the FFPE block did not identify any mutations in PDAC driver genes, whereas they were detected in the EpCAM-enriched sample (KRAS, SMAD4, and RNF43; Supplementary Table S2). In addition, no major AI events were found in the matched FFPE-derived samples compared with the enriched counterpart in both patients (Supplementary Fig. S2A and S2B).

All but one sample (1/9) of this pilot phase (patient 16, resected tumor) passed sequencing quality control (QC), which was therefore excluded from further analysis, with the exception of ddPCR profiling. The genomic landscape of mutations in known PDAC-associated genes in this pilot phase generally matches those previously reported (Fig. 2; refs. 3, 4, 19, 20).

In the second phase of this study, we applied our EpCAM-enrichment strategy to 23 independent real-world EUS-FNA samples collected during routine diagnostic procedures (18/23 with paired PBMC, median of one FNA pass, range 1–2). Comparable with the results of the pilot phase, EpCAM enrichment followed by ddPCR showed a significant increase in KRAS MAF (median, 13.9% vs. 28.0%; P = 0.03; Fig. 1C and D). These results were also seen by sequencing as we split up one pooled FNA with two research passes (patient 18) to conduct WES on matched enriched and nonenriched DNA. The median MAF of all overlapping mutations and SCNAs, in the enriched versus the nonenriched samples, increased significantly (P < 0.001, Supplementary Table S2; Supplementary Fig. S2C). Only three cases showed either no demonstrable KRAS or an MAF below 5% after EpCAM enrichment (range KRAS MAF, 0%–4.39%). Nondiagnostic EUS-FNA passes cannot be excluded in these cases as our single research passes were acquired after routine clinical passes and may therefore harbor lower levels of cancer cells (21).

All 23 enriched FNA samples had enough DNA to proceed with library preparation and passed sequencing QC. The genomic landscape of mutations in known PDAC-associated genes generally matched those from previous reports (3, 4, 10, 15). In our cohort, significantly recurrent mutations (>10%) were identified in KRAS, TP53, CDKN2A, SMAD4, GNAS, RNF43, and ARID1A (Fig. 2). Prevalence of KRAS mutations using sequencing alone was 82.6% (19/23 patients), whereas combination of both sequencing and ddPCR increased the KRAS mutation detection rate to 91.3% (21/23 patients). In order to confirm our methods for calling SNVs and CNVs in patients that underwent a subsequent resection or core biopsy, we used pregenerated targeted sequencing CLIA laboratory reports as the gold standard for validation. All but one SNV (patient 6: TP53) showed a concordance (92.3%) between the CLIA report and our sequencing results (Supplementary Table S3). In addition, comparison of sequencing with ddPCR results for KRAS MAF and GNAS p.R201C showed significant correlation (R² = 0.93; P < 0.0001; Supplementary Fig. S3A and S3B; Supplementary Table S4).

Although not present in a majority of patients with PDAC, low-frequency, potentially actionable alterations were found in individual cases. Based on knowledge of recently published actionable mutations (2) and by including only highly deleterious SNVs from the COSMIC database and amplifications with more than 3 copies, a total of 5 of 23 EUS-FNA patients (21.7%) showed potentially actionable alterations such as MYC amplifications = 1 (SCNA with a level of 5.5), MTOR = 1 (nonsense mutation, p.A835S), BRAF = 1 (nonsense mutation, p.R509L), MSH6 = 1 (nonsense mutation, p.P1087R), and POLE = 1 (frameshift insertion, p.F699Vfs*11; Fig. 3A). These results match prevalent actionable mutations previously reported in 17% to 48% of PDAC cases (19, 20, 22). Potential treatment options for these mutations include GSK525762 for MYC-dependent carcinomas; everolimus, VS5584, or LY3023414 for MTOR-mutant tumors; or vemurafenib for BRAF-mutant PDACs and agents targeting tumors with mismatch repair defects caused by MSH6 and POLE such as pembrolizumab and nivolumab which have shown response rates of over 30% in patients with noncolorectal cancer (2, 23). Nonetheless, it should be emphasized that estimates for the proportion of actionable mutations in PDAC vary widely depending on the inclusion criterion.

Taken together, we were able to show that EpCAM enrichment of epithelial tumor cells significantly increases MAF and facilitates genomic profiling in PDAC, especially in challenging sample types like EUS-FNA and samples with low overall tumor cellularity.

PDAC is known for a relatively modest mutational burden, compared with melanoma or lung cancer (1–4). Matching previous reports, our median mutational burden for (non-) synonymous SNVs (18/23 PBMC paired samples) was 1.63 mutations/Mb (range, 0.29 mut/Mb–21.46 mut/Mb; ref. 24). TMB was not affected by the presence of matching germline samples, stage, or tissue origin (Supplementary Fig. S3C–S3E).

Two of 23 individuals in our EUS-FNA cohort (8.7%) showed an increased mutational burden (TMB_high), defined as carrying more than 10 mutations/Mb by recent Pan-cancer analysis (Fig. 3B; ref. 25). Both TMB_high patients harbored SNVs in DNA-damage response (DDR) genes or mismatch repair defect genes, whereas the majority of TMB_low displayed no mutations in these genes (P < 0.001; Fig. 3C), which is in line with previously reported results (26). Patient 29 harbored a MSH6 germline mutation (Lynch syndrome) that has been confirmed by clinical CLIA germline testing as a variant of unknown significance and an additional nonsynonymous mutation (Supplementary Table S3). This patient had an extensive family history of multiple cancers including a mother with breast cancer, a sister with PDAC, a parental cousin with an aggressive melanoma, and two parental cousins with brain tumors. Patient 8 had a COSMIC-annotated frameshift insertion in POLE (COSM2001733, p.F699Vfs*11, no clinical CLIA testing performed) and also reported a family history for cancer (gastrointestinal cancer: mother and grandfather). In our limited EUS cohort, TMB_high patients demonstrated a trend toward improved PFS (median PFS, 674 days vs. 191 days; P = 0.09; Fig. 3D) but not OS (417 days vs. undefined; P = 0.19) compared with non-TMB_high (Supplementary Fig. S4A).

To ensure a rigorous standard for annotating SCNA, we used two independent algorithms: HapLOHseq (14), a powerful tool to analyze the presence of AI in low purity samples, and the GATK standard segmentation pipeline (Supplementary Fig. S4B). The results from these tools were merged to obtain a high-confidence SNCA call set, which was validated by a significant correlation of the MYC amplification estimates between sequencing and ddPCR results (n = 7, R² = 0.91; P < 0.0001; Supplementary Fig. S5A–S5C).

Consistent with previous studies, deleterious SCNAs were mainly detected in known tumor-suppressor genes such as CDKN2A (9/23, 39.1%), SMAD4 (7/23, 30.4%), TP53 (7/23, 21.7%), and ARID1A (6/23, 26.1%; all P < 0.05; refs. 3, 4, 19, 20). Interestingly, CDKN2A showed a large number of focal SCNA events (5/9, 55.5%), whereas most other loci showed a mixture of both focal and larger events. Not unexpectedly, amplifications were most commonly centered on established oncogenes such as GATA6 (5/23, 21.7%) or MYC (3/23, 13.0%; Figs. 2 and 4; Supplementary Fig. S6A). Mirroring results from previous PDAC publications, at least one arm-level SCNA occurred in 78.3% of patients. The majority of patients harbored multiple SCNAs, most frequently spanning chromosome arms 6q (9/23, 39.1%), 8p (8/23, 34.8%), 9p (7/23, 30.4%), 17p (11/23, 47.8%), 18p (7/23, 30.4%), 18q (11/23, 47.8%), and 19p (9/23, 39.1%; refs. 4, 17; Fig. 4; Supplementary Fig. S6).

Genomic instability from SCNA data across cancer types can be estimated using a number of different approaches. For PDAC, structural rearrangements contributing to an “unstable genome” (>200 rearrangements) have been used as a measure of genomic instability (3, 27). We were unable to derive this previously described “unstable genome” phenotype from the exome capture platform (SureSelect, Agilent), as it only covers 71 possible breakpoints. Therefore, other previously proposed quantitative metrics, such as HRD-LST (large-scale state transitions), HRD-LOH, the total SCNA burden, as well as aneuploidy were analyzed (17, 28, 29), and all of these quantitative assessments were readily feasible on EUS-FNA WES data. All scores were significantly correlated with each other, which is why we used aneuploidy in the subsequent analyses (all P < 0.015, Supplementary Table S5). In our cohort, aneuploidy score was not affected by patient's gender, primary tumor location, tissue source, or age (Supplementary Fig. S7A–S7D). Somatic TP53 mutation carriers showed a significantly higher aneuploidy score than patients without a TP53 mutation (median, 0; range, 0–8 vs. 11; range, 3–22; P = 0.0002), as has been previously reported (Supplementary Fig. S7E; ref. 17). The magnitude of aneuploidy measured by the number of aneuploid chromosome arms on WES has previously been used as an indirect measure of genomic instability, the prognostic value of which is an area of active research (17). Correspondingly, studies in many solid tumors (30) showed a correlation of high aneuploidy with later-stage disease and poorer prognosis [e.g., colorectal cancer (31) and esophageal cancer (32)]. We therefore evaluated the prognostic potential of the aneuploidy score. Classification of patients based on the median aneuploidy level into a low-level aneuploidy (aneuploidy_low: bottom half = 0–7 arm-level events) and high-level aneuploidy (aneuploidy_high: top half ≥ 8 arm-level events) revealed a significant and deleterious impact on prognosis (Fig. 5A–C). Aneuploidy_high tumors showed a median PFS of 104.5 days, whereas tumors with low-level aneuploidy experienced a significantly longer PFS of 365 days (P = 0.009; Fig. 5B). In addition, aneuploidy_high showed a trend for worse OS (225 days vs. undefined; P = 0.06; Fig. 5C). Importantly, the aneuploidy score maintained its prognostic trends also in localized patients (Fig. 5D and E), which suggested that aneuploidy's prognostic value is not based on advanced tumor stage. Correspondingly, after restricting the analysis to nonsurgically treated patients (n = 19) in this group, prognostic impact of aneuploidy remained significant (OS, 286.6 ± 58.8 days vs. 466.9 ± 69.5 days; P = 0.039 and PFS, 130.0 ± 51.4 days vs. 317.0 ± 24.4 days; P = 0.005). A potential confounder for our analysis may be the inclusion of TMB_high samples as these samples demonstrated longer PFS and therefore results have to be interpreted with caution. To further validate the prognostic impact of the aneuploidy score, we used an independent cohort from the DFCI in Boston, MA (20). There are several details that need to be highlighted when comparing the two datasets. In the DFCI cohort: (i) most of the specimens (96%) were taken using core needles or intraoperative biopsies, primarily from metastatic lesions, whereas only 4% of cases were sampled by EUS-FNA with 9% coming from primary pancreatic lesions; (ii) no enrichment of the biopsies was performed and the authors report a median cellularity of 40% (range, 5%–80%); (iii) most of the cohort (91%) and all of the patients with a cellularity >20% were diagnosed with metastatic disease; (iv) no UMI-enhanced sequencing was used for WES.

To minimize bias in the aneuploidy score related to low cellularity samples in the DFCI cohort, we restricted the analysis to specimens with a tumor cellularity of >20% (n = 52/73). The median aneuploidy score of the DFCI cohort was 15 (range, 2–33), and only 8 of 52 patients showed an aneuploidy score below 8. This may be due to the advanced stage in the DFCI cohort and mirrors our findings that patients with metastasis show a trend for higher aneuploidy score compared with localized patients [median of 5.5 (range, 0–18) vs. 11 (range, 0–22); P = 0.07]. PFS for most patients in the DFCI cohort was unavailable; therefore, we used OS data. Using our cutoffs for aneuploidy_high (≥8) and aneuploidy_low (≤7), there was no prognostic value for patients with aneuploidy_high (median OS in days: 625 ± 35.8 vs. 635 ± 74.8; P = 0.60). Because of the differences between the cohorts as detailed above, we then evaluated median aneuploidy score as a cohort-specific cutoff value and found a significantly worse OS survival for aneuploidy_high (≥15) patients (313.2 ± 38.4 vs. 792.7 ± 160.8; P = 0.0013; Supplementary Fig. S7F).

Collectively, we speculate that increased aneuploidy in localized patients may be associated with more aggressive disease and earlier metastatic spread.

Therapeutic decisions on first-line therapy in PDAC are currently mostly made on a clinical basis (e.g., performance status), without using biomarkers, unlike established markers for other cancer types, like colon cancer (33). There have been ongoing efforts using genomic or transcriptomic data to predict response to platinum-based therapy. In particular, cancers with BRCA1/2 mutations or HRD show superior response to platinum-based therapies in multiple cancers, including PDAC (22, 34–36). We therefore evaluated the potential of WES data within our unique EUS-FNA set to predict response to platinum regimens.

Eight patients (4 localized and 4 metastatic) whose EUS-FNA samples were sequenced were treated with platinum-based therapies. Three of these showed a response to FOLF(IRIN)OX—defined as stable disease or partial response based on RECIST 1.1 criteria—whereas 5 patients progressed. None of the patients harbored deleterious SNVs in the following genes: BRCA1/2, PALB2, FANCM, XRCC4/6, CHEK2, BRIP1, or BARD1, except for 2 patients with ATM variants (patients 26 and 29). There was significance between patients with higher TMB and platinum response (8.86 mut/Mb; range, 4.89–21.37 vs. 1.1 mut/Mb; range, 0.79–3.36; P = 0.04; Fig. 6A), which is in line with former reports (3). Previously proposed parameters including aneuploidy, HRD-LST, and HRD-LOH (17, 29) did not show a significant correlation (Fig. 6B; Supplementary Table S5) with platinum response, whereas signature 3 showed a marginal positive association with response (P = 0.06), a finding that is in line with previous results (37).

We therefore evaluated our SCNA data for new approaches to predict platinum response which might correlate with HRD or DDR deficiency. Based on the assumption that DDR deficiency and especially an impaired DNA double-strand break repair contribute to HRD, we analyzed the number of chromosomal arms with opposing SCNA segments (e.g., gain and deletion present in the same chromosome arm; Fig. 6C) and use the sum of arms exhibiting this pattern to calculate a complexity score (CS). Our assumption is that a higher CS is inversely correlated with DNA double-strand repair capabilities. The CS is novel with regards to the integration of copy-number (segmentation data) with AI calls derived using a haplotype-based approach. However, this approach is similar in concept to using large-scale state transitions as a signature for HRD, which identifies adjacent chromosome segments with different copy-number states (29). Patients with platinum therapy response showed a significantly higher CS than those with progressive disease (CS: 3, range 2–3 vs. 0, range 0–1; P = 0.04; Fig. 6D). In addition, a high CS showed a trend for better PFS (674 days vs. 175 days; P = 0.11; Fig. 6E), OS (417 days vs. undefined; P = 0.14; Supplementary Fig. S7G), and a positive TMB correlation with tumors harboring a mutation in DDR genes with highest CS (R² = 0.37; P = 0.002; Fig. 6F). In addition to the 8 patients treated with a platinum-based therapy, 8 patients of the FNA cohort received a gemcitabine-based regimen. Seven patients were excluded from this analysis due to trial treatment, upfront surgery, side effects, fragility, or patient's choice not to undergo any therapy. When examining CS for chemotherapy prediction, we observe a marginally significant association between the CS score and response: median 1.5, range 0–3 for responders and median 1, range 0–1 for nonresponders (P = 0.07; Supplementary Fig. S7H). We also attempted to test the predictive value in the DFCI cohort (see above). There were three responders in the dataset (with 1 having an estimated tumor purity of only 4%); this severely restricted power to test for an association between the CS score and response to platinum therapy. However, we do note that patients with BRCA1/2 mutations (4 germline and 1 somatic) in the DFCI cohort have significantly higher CS scores than patients without these mutations (P = 0.008, mean = 3.6 CS for those with a mutation and mean = 0.7 CS for those without). In addition, we performed mutational signature analysis as described previously (13) on FNA samples from treatment-naïve tumors (Supplementary Fig. S8). Herein, signature 3, which is associated with double-strand repair deficiency, showed a positive association with the CS (R² = 0.49; P = 0.07).

Taken together, CS shows potential as a predictive biomarker to estimate platinum responsiveness based on an indirect DDR evaluation even from limited FNA samples.

To demonstrate the clinical significance and feasibility of WES from limited clinical biopsies, we interrogated two consecutive endoscopic lung biopsies in a patient (patient 17) with recurrent metastatic PDAC to the lung (both biopsies were included in the pilot phase). The first biopsy [timepoint 1 (T1)] was taken prior to immunotherapy with a PD-L1 antibody (durvalumab) and a STAT3 inhibitor (AZD9150; ref. 38). The second endoscopic biopsy [timepoint 2 (T2)] was taken after progression on aforementioned therapy and prior to initiation of a hypoxia-activated prodrug regime with evofosfamide (39). Before immunotherapy initiation, this patient was diagnosed with a multifocal pulmonary recurrence 2 years after upfront distal pancreatectomy. The recurrent tumor was then treated with capecitabine and oxaliplatin but progressed within a year (Fig. 7A).

Comparative analysis of the biopsy taken before immunotherapy treatment (T1) and after progression (T2) revealed new-onset focal deletions of STK11, TSC2, and the CD274 locus (PD-L1; Fig. 7B). These loci have been previously reported to have a great impact on immunotherapy resistance. STK11 loss, for example, has been repeatedly linked to the development of immunotherapy resistance in KRAS-mutant lung cancer (40), and PD-L1 loss is associated with adaptive immune resistance in Hodgkin lymphoma and other neoplastic diseases (41). In addition to these focal deletions, chromosome 3q, 10q, 17q, 6p, and 6q also developed new arm-wide SCNAs, indicating LOH at all MHC genes, as evident by high confidence AI calls. Excluding these focal changes, T1 and T2 had overall congruous AI profiles. T1 and T2 biopsies shared focal deletions in major PDAC tumor suppressors such as SMAD4 and CDKN2A, as well as amplifications in GATA6 and MYC. T1 and T2 SNVs in established drivers, KRAS, SMAD4, and TP53, were also shared (Fig. 7C). Only two of the nonoverlapping deleterious SNVs had a COSMIC annotation (CCDC42: p.R104W, COSM4443319 and RCBTB2: p.V347M, COSM6445879); upon inspection and literature review, these do not have a known mechanistic link that would explain the development of immunotherapy resistance. Overall, mutational burden remained roughly the same between T1 and T2 (2.1 mut/Mb vs. 1.6 mut/Mb). As illustrated by this example, paired longitudinal endoscopic biopsies obtained over the course of therapy in patients with PDAC might elucidate genomic alterations that could explain the development of resistance to immunotherapy (or targeted therapies).

Over the past decade, precision oncology has improved clinical outcomes for many patients with solid cancers. Patients with melanoma and breast cancer are now routinely screened for their V600E proto-oncogene B-Raf (BRAF; ref. 42) and HER2 status, respectively (43). Targeting these aberrations has alleviated prognosis in both entities. Despite significant efforts, PDAC has proven less susceptible to targeted or immune therapies, mainly due its stroma-dense, immune “cold” tumors that additionally show a remarkable genomic heterogeneity (1–6).

Previous WES, whole-genome sequencing (WGS), and whole transcriptomic approaches in molecular profiling of PDAC have used tissue from multiple imaging-guided core biopsies for DNA and RNA profiling. The amount of tissue from these multiple cores is substantially greater than our “one-pass” EUS-FNA starting material (11, 13, 20, 22). It is also important to mention that multiple EUS-FNA passes are not necessarily associated with an increased amount of neoplastic cells harvest (44, 45) nor with an increased sensitivity of detection (21), and indeed, might increase the associated risk of complications (21, 46). Other sampling techniques such as EUS-FNB have shown promising results for in-depth genetic analysis, but their applicability and widespread use can be constrained by their significant expense (47, 48). In addition, a major advantage of EUS-FNA compared with EUS-FNB is that the former allows for ROSE-rapid on-site cytology assessment that ensures that there is adequate cellularity and a diagnostic specimen. Here, we report our efforts to characterize EpCAM-enriched real-world, limited EUS-FNA. Despite the fact that most patients with PDAC are diagnosed using EUS-FNA, these samples have so far been deemed not suitable for in-depth molecular analysis such as WES (49, 50). To increase the tumor to normal cell proportion, many studies have either enriched their specimens, e.g., laser microdissection (time- and labor-intensive), or excluded low tumor cellularity samples from analysis altogether (2–6, 20). In contrast, we used a simple and fast enrichment method to successfully isolate tumor cells resulting in an improved detection of deleterious mutations.

There is evidence that evaluation of aneuploidy as a marker for genomic instability may be a valid surrogate in metastatic PDAC (22, 51). Aneuploidy demonstrated prognostic significance independent of disease stage within the current study and an independent dataset of patients with pancreatic cancer. This finding is in line with results showing a worse prognosis and later disease stage correlated with increased aneuploidy (3, 24, 32). In this regard, tumors with high aneuploidy have shown increased tumor progression and chemotherapy resistance due to a higher intratumoral heterogeneity and ultimately an elevated capacity for adaptation (52). It is important to mention that aneuploidy seems to be context-dependent and therefore might differently affect distinct tumor types (32).

One of the important clinical questions in PDAC is the susceptibility to platinum-based therapies as one of two first-line regimens, which should ideally be discernible in WES data on limited samples material like EUS-FNA. In light of recent findings, up to 25% of patients with PDAC harbor a genomic footprint of HRD (3), and many of these potential responders might be missed using only “classic” DDR mutations for identification, like BRCA1/2 or PALB2 that are present in 5% to 7% of patients with PDAC (53). Other HRD-causing alterations like BRCA1 promoter hypermethylation or biallelic loss of PALB2 are even more challenging, if not impossible, to detect by most commonly used targeted sequencing panels (2). In a cohort of patients with breast cancer, over 30% of HRD-positive tumors had no underlying genetic alteration associated with HRD and yet demonstrated platinum susceptibility (35). Within this aforementioned study, HRD has been estimated by a WGS algorithm called “HRDetect,” but the algorithm performed suboptimally with WES data because essential predictor components such as rearrangement signatures and indels analysis are restricted (35, 36).

We identified a relatively facile WES-based scoring system of the “complexity” of SCNA events within a chromosomal arm that predicts platinum response in this cohort of patients with PDAC. This CS analyzes the presence of alternating SCNA events (gains, losses, and cnLOHs > 1M) within a single chromosomal arm boundary. This kind of event pattern can only occur upon a severe disruption of double-strand DNA repair and results in missegregation of “shattered” arm fragments in mitosis. It is relatively easy to calculate using WES data and standard GATK segmentation pipeline which makes it an attractive possible biomarker. We therefore hypothesize that CS may be a marker for increased genomic instability indicating susceptibility to platinum-based treatment. Interestingly, CS positively correlated with increased TMB that by itself also showed a predictive value but did not show a correlation with other SCNA markers, like aneuploidy.

We acknowledge important drawbacks in our study. One downside of our approach is that it requires a freshly acquired aspiration biopsy with dissociated cells as starting material, as we did not attempt WES on archival cytology samples. On the other hand, dissociation of cells from FNAs is a relatively simple procedure when compared with the dissociation of core biopsies which is time consuming. We were also not able to conduct a transcriptional analysis to validate loss or gain of function of genes because EUS-FNA provided sparse material to conduct additional assays. In addition, our WES library design refrained us from analyzing important RNA fusion patterns which are accessible primarily using WGS.

Although significant EpCAM expression is seen in epithelial-derived carcinoma cells when compared with normal epithelial, up to one third of pancreatic cancers may have low to no EpCAM expression (54, 55). This underlines a potential limitation of the EpCAM-enrichment approach of tumor cells and may warrant alternative strategies such as negative selection by targeting cells expressing CD45, CD31, and FAP to clean the sample of immune/blood cells, endothelial cells, and fibroblasts. Based on the minimal amount of cellular material acquired from a single EUS-FNA pass, we were not able to process each pass with multiple selection methods to compare intrapatient accuracy and performance. In addition, with advances and improvements in sampling strategies (like EUS-FNB), enrichment in the context of adequate tissue acquisition may become redundant. Finally, due to the small number of overall patients, the non-CLIA setting of our assay, and the inherent cohort heterogeneity, our conclusions need to be validated in a larger prospective EUS-FNA cohort.

In conclusion, in this proof-of-concept study, we show the importance and ease of EpCAM enrichment for genomic analysis of real-world limited aspiration (EUS-FNA) samples. Using enriched EUS-FNA samples, we were able to perform an in-depth whole high-quality molecular-barcoded WES analysis, reproducing genomic patterns consistent with previous results. We show the prognostic potential of SCNA (namely aneuploidy) for PDAC diagnosis and also identify a novel predictive biomarker (CS) for platinum responsiveness in our cohort. CS may be a relatively facile biomarker for predicting platinum responsiveness in newly diagnosed PDAC not harboring overt mutations in canonical DDR genes.

A. Semaan reports grants from DFG during the conduct of the study. A. Maitra reports other from Cosmos Wisdom Biotech, Thrive Earlier Detection, and JNJ outside the submitted work. Y.A. Jakubek reports other from Invitae outside the submitted work. No disclosures were reported by the other authors.

A. Semaan: Conceptualization, data curation, investigation, writing-original draft. V. Bernard: Conceptualization, methodology. J.J. Lee: Conceptualization. J.W. Wong: Data curation, formal analysis, methodology. J. Huang: Data curation, formal analysis. D.B. Swartzlander: Data curation, investigation. B.M. Stephens: Data curation. M.E. Monberg: Data curation. B.R. Weston: Resources. M.S. Bhutani: Resources. K. Chang: Data curation, formal analysis. P.A. Scheet: Supervision, funding acquisition. A. Maitra: Supervision, funding acquisition, project administration. Y.A. Jakubek: Data curation. P.A. Guerrero: Supervision.

We thank A.J. Aguirre, B.M. Wolpin, and the Hale Family Research Center at Dana-Farber Cancer Institute for providing an independent dataset used for replication of our results. Most of all we thank the patients and their families. A. Maitra is supported by the MD Anderson Pancreatic Cancer Moon Shot Program, the Khalifa Bin Zayed Al-Nahyan Foundation, and the NIH (U01CA196403, U01CA200468, and P50CA221707). A. Semaan is supported by the German Research Foundation (SE-2616/2–1). V. Bernard is supported by the NIH (U54CA096300, U54CA096297, and T32CA217789). J.J. Lee is supported by the NIH (T32CA009599).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.