Patient sets and study design

We used the same reference samples as in our previous publication, containing 205 primary PAAD, 144 primary iCCA, and 50 normal bile duct samples from 7 different studies [1].

For developing the anomaly detection layer, we put together the anomaly detection training samples, which included 20% of the samples from 10 different carcinomas from TCGA. These 10 carcinomas formed the negative class and were selected because of their high metastatic potential and their prevalence as possible differential diagnoses of PAAD in a metastatic setting. Briefly, we included: breast invasive carcinoma (BRCA, n = 159), esophageal carcinoma (ESCA, n = 35), lung adenocarcinoma (LUAD, n = 97), stomach adenocarcinoma (STAD, n = 79), liver hepatocellular carcinoma (LIHC, n = 78), colon adenocarcinoma (COAD, n = 63), rectum adenocarcinoma (READ, n = 21), uterine corpus endometrial carcinoma (UCEC, n = 93), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC, n = 62), and prostate adenocarcinoma (PRAD, n = 100). The reference samples were added to the anomaly detection layer to represent the positive class.

The biological validation cohort (n = 3579) consisted of 252 primary PAAD, 151 primary iCCA, and 20 normal bile duct samples partially representing the validation set for the initial classifier development [1], to which we added 16 primary PAAD samples, 20 PAAD liver metastases (PAAD met.Liv.) and 36 primary iCCAs (GSE217384) that represented our previous testing cohort, 26 previously published formalin fixed paraffin embedded (FFPE) primary PAAD from Benhamida et al. [13] and the additional 80% of the samples from the 10 mimicker carcinomas (n = 3120).

The technical validation cohort consisted of 6 not otherwise specified lung cancers and 5 BRCA samples that were analyzed using the Illumina Infinium MethylationEPIC v2.0 BeadChip (EPICv2) (Illumina, CA, USA) array (GEO GSE222919), as well as 4 in-house PAAD metastases: peritoneal carcinomatosis (PAAD met.PC), lymph node metastasis (PAAD met.LN), spleen metastasis (PAAD met.Spleen) and PAAD met.Liv..

The performance of the classifier was tested on a test set. As positive control we used 16 confirmed or clinically highly suspected PAAD metastases: 11 PAAD met.PC, 2 PAAD met.LN, 2 PAAD lung metastases (PAAD met.Lung), and 1 PAAD met.Liv..

The negative control, used to test the classifier’s performance in excluding potential mimickers, contained:

(i)

a unique set of 96 BRCA metastasis from 18 different anatomical sites (adrenal gland n = 3, bone n = 3, brain n = 12, chest n = 7, gastrointestinal tract (GI tract) n = 2, kidney n = 1, liver n = 24, lung n = 11, lymph node n = 12, ovary n = 3, pancreas n = 1, pericardium and pleura n = 3, peritoneum n = 2, skin n = 2, soft tissue n = 7, spleen n = 1, thyroid n = 1, uterus n = 1) from the AURORA US network resource. This set contained BRCA of all different molecular subtypes, samples were both FFPE and fresh frozen and were obtained both from autopsies and pathology specimens. Methylation was performed using the Illumina Infinium MethylationEPIC BeadChip (EPICv1) (Illumina, CA, USA) array [14].

(ii)

An external set of 13 carcinoma brain metastases (GSE249157): 1 BRCA, 7 LUAD, 4 COAD and 1 PRAD. All samples were fresh frozen and the methylation was performed using the EPICv1 array.

(iii)

An internal set of 15 carcinoma brain metastases: BRCA, LUAD, STAD, UCEC, PRAD, mucinous ovarian cancer (MOC) and CUP. A detailed presentation of these samples can be found in Additional file 1: Table S1.

The study was approved by the ethics commissions of Charité, Universitätsmedizin Berlin (EA1/079/22).

Tissue microarray construction and immunohistochemistry

For all positive control samples, for which material was available (n = 11), we constructed one tissue microarray (TMA) containing three 1.5 mm cores for each tumor. For the other samples (n = 5) that were too thin to be included into the TMA we performed whole slide staining. Next, the FFPE TMA and tumor blocks were cut into 2.5 μm sections. One section was used for hematoxylin and eosin (H&E) staining and ten others for p53, SMAD4, GATA6, Ki-67, CK7, CK20, Annexin 1 (ANXA1), Annexin 10 (ANXA10), CD3 and CD20 immunohistochemistry staining. Additionally, the PAAD met.Liv. samples from the biological validation set were stained for p53, SMAD4, GATA6, CK7, CK20, CD3 and CD20. Other immunohistochemistry (IHC) data for the PAAD met.Liv. from the biological validation samples were previously generated [1].

For the immunohistochemical staining, a BenchMark XT immunostainer (Ventana Medical Systems, Tucson, AZ) was used. For antigen retrieval, sections were incubated in CC1 mild buffer (Ventana Medical Systems, Tucson, AZ) for 30 min at 100 °C, or were incubated in protease 1 for 8 min. The sections were stained with anti-Ki-67 antibody (M7240, Dako, 1:50, CC1 mild buffer), anti-p53 (M7001, Dako, 1:50, CC1 buffer), anti-SMAD4 (Ab40759, Abcam, 1:200, CC1 buffer), anti-GATA6 (Q92908, R&D Systems, 1:100), anti-CK7 (M7018, Dako, 1:1000, protease 1), anti-CK20 (M7019, Dako, 1:1000, protease 1), anti-Annexin A10 (PA5-52151, Invitrogen, 1:2000), anti-Annexin I (610066, BD Biosciences, 1:5000), anti-CD3 (A045201-2, Dako, 1:100), and anti-CD8 (M7103, Dako, 1:100) for 60 min at room temperature, and visualized using the avidin–biotin complex method and DAB. We stained the cell nuclei by additionally incubating for 12 min with hematoxylin and bluing reagent (Ventana Medical Systems, Tucson, AZ). Histological images were acquired with the digital slide scanner PANNORAMIC 1000 (3DHISTECH).

Histological analysis and immunohistochemistry scoring

The ANXA1/10 immunohistochemistry score was proposed as a potential tool for detecting metastatic PAAD. For this purpose, we used the scoring and classification system proposed by Padden et al. [15]. The intensity [0 (none), 1 (weak), 2 (intermediate), or 3 (strong)] and percentage of positive tumor cells [0, 1 (≤ 5%), 2 (6–10%), 3 (11–50%), or 4 (> 50%)] for each tumor was scored separately and the two scores were multiplied, resulting the immunoreactive score (IRS). The IRS thresholds proposed by Padden et al. and validated by us in a previous paper [1] were used also in this study. Hence, an IRS of 5 or higher for Annexin 1 and an IRS of 0.5 or higher for Annexin 10 was suggestive for PAAD. According to the previous studies, only one of the two markers needed to be equal or higher to the IRS cut-off.

For Ki-67 the percentage of positive tumor cells was estimated in representative hot spots. For p53 complete loss or intense nuclear staining were considered to be specific for a mutated pattern. For SMAD4 complete loss was considered to be specific for a mutated pattern. For CK7 and CK20 any degree of cytoplasmic positivity was scored as positive. GATA6 was scored as previously described [16]. Briefly, semiquantitative scoring from 0 (negative) to 4 (intense nuclear) was performed. The samples with scores from 0 to 2 were considered GATA6 low, and the ones with score 3 and 4 were defined as GATA6 high.

We considered an IHC pattern to be specific for PAAD if a tumor showed an ANXA1/10 score that supported the diagnosis of PAAD, and additionally SMAD4 loss and/or CK7 positivity. We considered an IHC pattern to be inconclusive if a tumor showed an ANXA1/10 score that did not support the diagnosis of PAAD, and the tumor showed both SMAD4 loss and CK7 expression, or if the tumor showed an ANXA1/10 score that supported the diagnosis of PAAD, and SMAD4 was expressed and CK7 was negative. We considered an IHC pattern to be unspecific if a tumor showed an ANXA1/10 score that did not support the diagnosis of PAAD, and the tumor showed only SMAD4 loss or CK7 was expressed or none.

Organoids

Organoids were established from a primary tumor and matched metastases (peritoneal carcinomatosis and liver metastases) from surgical specimens in accordance with the ethics approval EA1/157/21. The tissue was cut into small pieces using scalpels and digested with 100 µg/ml DNAse (STEMCELL Technologies, Vancouver, Canada), 125µg/ml Collagenase II (Sigma-Aldrich, Merck, Darmstadt, Germany), 1:2000 Rock-Inhibitor (Abmole Bioscience, Houston, TX, USA) and 1:200 Amphotericin B (Sigma-Aldrich, Merck, Darmstadt, Germany). The specimens were then incubated for 2 to 3 h. Cells were filtered through a sterile 100µm filter. Red blood cell lysis (Miltenyi Biotec, Bergisch-Gladbach, Germany) was performed if necessary and cells were plated in Cultrex (R&D Systems, Minneapolis, MN, USA). Culture medium as described by Broutier et al. [17] was added after solidification of domes. Amphotericin B was added for the first 7 days of culture to prevent fungal contamination. The culture medium was exchanged every 3 to 4 days and regularly checked for Mycoplasma contamination using the Mycoplasma detection kit (Applied biological Materials, Richmond, Canada). Organoids were split when they reached a size of 200 µm using TrypLE Express (Thermo Fisher Scientific, Waltham, MA, USA) and plated in a ratio of 1:2.

For histologic embedding, organoids were incubated PFA (4% in PBS) at 4 °C. After detachment they were transferred into prewarmed histogel (Thermo Fisher Scientific, Waltham, MA, USA). The organoids were then FFPE. The blocks were cut into 3 μm sections. The slides were stained with H&E using Tissue-Tek Prisma® Plus Automated Slide Stainer (SAKURA). For the immunohistochemical staining, BenchMark XT immunostainer (Ventana Medical Systems, Tucson, AZ) was used. The sections were stained with anti-Ki-67, anti-CA 19-9 (1116-NS-19-9, Dako, 1:500), anti-p53 antibody, and anti-GATA6 antibody. Finally, representative images were acquired with the digital slide scanner PANNORAMIC 1000 (3DHISTECH).

DNA extraction

For all samples, tumor areas were marked and the tumor cell content was determined using a light microscope (Olympus, BX46). Based on this information we determined the number of necessary slides for DNA extraction. Depending on the tumor purity and tumor surface we used between 7 and 20, 5 μm thick slides per sample from which the tumor contour was scratched for DNA extraction (Additional file 1: Table S2). Semi-automated DNA extraction was performed according to the manufacturer’s instructions (Maxwell RSC FFPE Plus DNA Purification Kit, Custom, Promega). DNA quantities were measured using Qubit HS DNA assay (Thermo Fisher Scientific).

DNA methylation

Whenever possible we used 500 ng of DNA for the DNA methylation analysis as input. For samples where there was not sufficient material available, we decreased the DNA input to as low as 182.8 ng of DNA. We used the Illumina Infinium HD FFPE DNA Restore Kit (Illumina, CA, USA) for DNA restoration from FFPE samples. Following this step, the EpiTect Bisulfite Kit (Qiagen) was used for bisulfite conversion. For the organoid models we used 100 ng of DNA extracted from fresh tissue. The EPICv1 array was used according to the manufacturer’s instructions for the DNA methylation analysis of the positive and in-house negative control samples. The EPICv2 array was used for the hybridization of the 4 PAAD metastases of the technical validation cohort.

Methylation array processing

Methylation data preprocessing was performed in R using various packages implemented in ChAMP [18]. Raw signals from all the IDAT files are loaded using the minfi package. In the training set, the EPICv1 and the Illumina Infinium HumanMethylation450 BeadChip (Illumina, CA, USA) samples were merged.

In the sample preprocessing for differential methylation analysis, several CpG sites were excluded: those on EPICv1 array not present in 450k arrays; any CpG sites with a detection p value greater than 0.01; low quality sites, defined as having fewer than 3 beads in at least 5% of the samples; all SNP-associated sites; multi-hit sites; and CpGs found on chromosomes X and Y.

While preprocessing samples for the classifier, no filtering was necessary, as only the 2048 CpGs which serve as features are selected at the end of the preprocessing pipeline.

Finally, the beta values were extracted and normalized using FunNorm and BMIQ, which together enhanced the process. Each cohort was pre-processed independently.

t-distributed stochastic neighbor embedding (t-SNE)

To generate the t-SNE plots, beta values of CpG sites were broken down into eigenvectors, and then handled using the R package Rtsne [19] using 5000 iterations. The count of eigenvectors (k) and the perplexity (p) were chosen individually for each plot to accommodate the varying number of samples.

The t-SNE in Fig. 2A and B were created using 30 eigenvectors (k = 30) and a perplexity of 15 (p = 15), while for t-SNE in Additional File 1: Fig. S4, we used k = 50 and p = 20, due to the larger number of samples plotted. The 2048 classifier features served as input for the eigenvector decomposition.

The t-SNE in Fig. 3A consists of 16 primary PAAD, 20 liver and 12 peritoneum metastases unmatched samples. The top 2000 CpGs with the highest standard deviation among these samples were selected and decomposed in 15 eigenvectors (k = 15) and then further reduced to two dimensions using t-SNE with a perplexity of 5 (p = 5).

Tumor purity estimation

We estimated the tumor purity using the InfiniumPurify R package [20]. For estimating the purity of iCCA and liver PAAD metastases samples, we selected “CHOL” as tumor type, and “PAAD” as tumor type for the primary PAAD samples and non-liver PAAD metastases.

Updated classification pipeline

The anomaly detection layer developed in the previous study could only differentiate between STAD, COAD tumors from PAAD, iCCA and normal bile tumors [1]. This layer has been replaced with an updated version that can separate PAAD, iCCA and normal bile tissue (positive class) from 10 different mimicker carcinomas (negative class).

The new model was built by training a neural network ensemble on the same 2048 CpGs identified in the previous study [1] and used by the classification layer. These are the top 2048 CpGs with the highest standard deviation among the reference samples in the training set. The models were trained using the anomaly detection training dataset consisting of all the reference samples (n = 399), which formed the positive class, and 20% of the samples of each of the 10 mimicker carcinomas, selected at random, forming the negative class (n = 787).

The model ensemble was built using the python library keras [21] using fourfold cross-validation, while python library optuna [22] was used to conduct hyperparameter optimization. The ensemble consists of 4 neural networks, each network having 5 layers (search space between 1 and 9 layers) with 2048 neurons per layer. A learning rate of 0.0088, a dropout rate of 0.2 (search space: 0, 0.1, 0.15, 0.2, 0.25, 0.3), L1 regularization of 0.00067 (search interval between 0 and 0.1), and 228 epochs (search space between 20 and 300) were found to be the optimal hyperparameters.

The other segments of the pipelines remained as they were in the previous study [1]. The python library reComBat [23] is used to fit a regularized empirical Bayes model to reduce sample storage material induced batch effects. The neural network model in the classification layer was developed in the previous study using python keras and optuna libraries [21, 22]. The optimal network was found to consist of 8 layers (search space 1 and 10 layers) with a starting width of 256 neurons (search space 64, 128, 256, 512, 1024, 2048, 4096) for the first hidden layer, incrementally decreasing to 16 neurons in the last hidden layer. It was trained with a learning rate of 0.00895 (search space between 0.0001 and 0.01). Finally, a dropout rate of 0 (search space 0, 0.1, 0.15, 0.2, 0.25, 0.3), L1 regularization of 0.00441 (search space between 0 and 0.1), and 191 epochs (search space between 20 and 300) were found to perform optimally. To further increase the accuracy and confidence in the model’s output, a threshold of 0.8 for the PAAD and iCCA classes and 0.5 for the normal bile were selected. Predictions which did not reach the threshold were put into the “No Match” class, together with the samples rejected by the anomaly detection layer described above.

The updated classification pipeline is therefore composed of three parts: (i) the anomaly detection layer that singles out PAAD, iCCA, and normal bile tissue from other carcinomas; (ii) the classification layer capable of differentiating between PAAD, iCCA, and normal bile samples; and (iii) a threshold-based filtering layer that weeds out samples with low confidence predictions. The result can therefore belong to one of four classes: PAAD, iCCA, normal bile tissue, or “No Match”. The “No Match” class contains all the samples rejected by the anomaly detection layer and the samples that passed the anomaly detection layer but did not reach the level required by the threshold-based filtering layer.

Copy number analysis

We calculated the copy number profiles from DNA methylation array data using the conumee package, version: 1.3.0. [24]. A set of 63 control samples derived from histologically confirmed normal pancreas tissue were used as baseline reference. The evaluation of copy number alterations was carried out manually with consideration of the tumor cell content for the evaluation of chromosomal gains or losses. In general, changes were considered relevant if the intensity ratio of a segment deviated from the baseline by at least more than 0.15 [25]. In addition, we created summary copy number profiles for three different groups: primary PAAD (n = 16), PAAD met.PC (n = 11), and PAAD met.Liv. (n = 21). This analysis was done using an adaptation of the conumee script (provided by Dr. Damian Stichel, Neuropathology Heidelberg). For the comparison of specific gene deletions and amplifications between the three groups we performed Fisher's exact test with Bonferroni correction for multiple testing.

Differentially methylated CpG probes and pathway analysis

The differentially methylated analysis (DMA) followed by pathway analysis was conducted on 48 samples from 3 groups, 16 primary PAAD tumors, 21 PAAD met.Liv. and 11 PAAD met.PC. The R package limma as implemented in ChAMP was used. Limma deploys a linear model alongside an empirical Bayes approach to gauge the mean methylation disparity between groups, following which it computes adjusted p values to accommodate multiple testing. An adjusted p value below 0.01 and an absolute logFC value exceeding 0.2 were chosen as thresholds to select the differentially methylated CpGs between groups.

We used the Illumina Infinium HumanMethylationEPIC manifest to annotate promoter and enhancer CpGs. Genes associated with differentially methylated promoter- and enhancer-associated CpGs between groups (primary PAAD, PAAD met.Liv., and PAAD met.PC) were analyzed for enriched pathways using Enrichr [26]. The Reactome (2022) Pathway Database was selected for the enrichment analysis. Volcano plots were created using VolcaNoseR [27]. From the top 10 enriched pathways we labeled the ones linked to epithelial and mesenchymal phenotypes. In addition, we used a second method for pathway analysis, methylGSA [28]. We used the same genes associated with differentially methylated promoter- and enhancer-associated CpGs between the three groups for this analysis. The KEGG pathway database was selected for the tested gene sets.

Classifier website

The website is a Vue.js application built with Nuxt.js running on a Node.js platform. It utilizes client-side rendering and leverages Google Firebase for secure user authentication. The application backend responsible for processing the raw data and running the prediction was developed in python using the FastAPI framework. The application is hosted on Google Cloud.

Statistical analysis

We performed statistical analyses and graphics using the GraphPad Prism 9 software. First, we determined whether the data followed a normal distribution, using the Shapiro–Wilk normality test. For the comparison between two groups, p‐values were determined with an unpaired t test if the data were normally distributed, while the nonparametric Mann–Whitney–Wilcoxon test was applied on data with a non‐normal distribution. For the comparison between multiple groups p‐values were determined using ordinary ANOVA test for normally distributed data, and the Kruskal–Wallis test for data with a non‐normal distribution. Correlations were performed using the Pearson correlation test. All tests were two‐sided, and a p value < 0.05 was considered statistically significant.