Human Pancreatic Cancer Single-Cell Atlas Reveals Association of CXCL10⁺ Fibroblasts and Basal Subtype Tumor Cells

Pancreatic ductal adenocarcinoma (PDAC) is projected to become the second leading cause of cancer-related deaths in the United States by 2030 (1). More than 85% of patients are diagnosed at a metastatic disease stage, with a dismal 3% 5-year survival rate (2). The most effective treatment option for PDAC is surgical resection. However, within 2 years of surgery, 30% of patients will experience disease progression (2). PDAC lacks effective early biomarkers and is highly resistant to chemotherapy, radiation, and molecular/immune targeted therapies (3, 4), possibly because 80% of the tumor mass is made up of stromal components that protect the sparse population of tumor cells (4). Several groups have defined two major molecular subtypes of tumor cells based on their transcriptomic profiles, often coexisting in a single tumor, that differ in their response to therapies [for review, see (5)]. The classic subtype tends to have well-differentiated histology in contrast to the basal-like subtype, which is undifferentiated, more aggressive, more resistant to current therapies, and associated with worse clinical outcomes (6). This stratification of tumors based on transcriptional subtypes provides a molecular basis for treatment resistance for the basal-like tumors (7). To further confound this aspect of tumor biology in PDAC, single-cell RNA sequencing (scRNA-seq) studies have demonstrated that classic and basal gene expression programs coexist in individual tumors, sometimes in individual cells, suggesting that cellular plasticity is a key driver of treatment resistance (8).

Complex cellular cross-talk between tumor and stromal cells drives tumor growth and metastasis in PDAC. PDAC is characterized by an abundance of cancer-associated fibroblasts (CAF), extracellular matrix (ECM), and immune-suppressive cell types, including macrophages, regulatory T cells, and neutrophils (4). Recent work has revealed that within the CAF compartment, multiple populations of transcriptionally distinct CAFs exist, including types of inflammatory CAFs (iCAF) and myofibroblastic CAFs (myCAF), whose functional roles are still under active investigation (9–11). There remain inherent limitations of prior scRNA-seq studies, given the smaller sample sizes of those datasets. To address these limitations, we applied a “multitranscriptomics” approach, integrating data from an extensive collection of publicly available human pancreas scRNA-seq data (10–21), four independent PDAC bulk RNA-seq datasets (15, 16, 22, 23), and two spatial transcriptomics (ST; ref. 17) datasets to determine which compartment(s), if any, of the PDAC TME are enriched within the basal subtype tumors. As anticipated, we validated the existence of basal tumor cell populations in PDAC, the gene signatures of which correlate in bulk PDAC RNA-seq data with worse outcomes, compared with gene signatures derived from classic subtype tumor cells. In the current study, we reveal a novel enrichment and spatial association between basal tumor cells and CXCL10-enriched myofibroblasts in the PDAC TME. We determined that the CXCL10 receptor, CXCR3, was expressed by tumor cells in addition to immune cells. Our scRNA-seq atlas revealed that patients with basal-enriched tumors have fewer CD8 T cells, which may in part explain why they have a worse prognosis. Lastly, we show the utility of using scRNA-seq tumor cell gene signatures derived from a large number of primary tumor (PT) patient samples (n > 170) to map changes in epithelial cell plasticity in response to drug treatments from bulk RNA-seq data of in vitro models of PDAC.

The Institutional Review Board (IRB) at Henry Ford Health granted a non-human subjects determination (IRB 16150) for usage of deidentified samples and data collected through the Henry Ford Health Translational and Clinical Research Center Biorepository, IRB 13297. Studies were conducted in accordance with the U.S. Common Rule. Written informed consent was obtained from each patient prior to genomic profiling.

Raw scRNA-seq fastq files were downloaded from 12 studies who used 10X genomics sequencing. Using 10X Genomics CellRanger version 6.1.1, all samples were aligned to the human reference genome version GRCh38. Using Seurat version 3 (18), we preprocessed each dataset separately, removing cells expressing fewer than 200 genes, with mitochondrial percent greater than 25%, and total counts greater than the 95th percentile for the dataset. Doublet detection was then performed on a per sample basis using doubletFinder_V3 (19). All remaining cells were merged, normalized using Seurat’s NormalizeData function, and scaled using Seurat’s ScaleData, adjusting for number of features, total UMI counts, and mitochondrial percentage. Principal component analysis (PCA) was performed using RunPCA in Seurat for dimensionality reduction. Following established best practices (20, 21), batch correction was performed using harmony (24) to correct the PCA embedding space, treating each sample as a batch due to high variability in the data. Nearest neighbor analysis was performed with Seurat’s FindNeighbors. Clustering was performed in Seurat using FindClusters using the Leiden algorithm and projected into the Uniform Manifold Approximation and Projection (UMAP) space with RunUMAP with default parameters using corrected components that explained 90% of the variation of the data. Using published markers, we identified coarse cell types (e.g., myeloid and ductal). For each broad cell type, fine clustering was performed, in which we again normalized, scaled, harmonized, and clustered the data to identify fine labels. Differential gene expression analysis was performed, comparing each cluster with all other clusters, using a Wilcoxon rank-sum test. Genes that had an FDR-corrected P value less than 0.05 were considered differentially expressed. Violin plots were generated using Seurat’s VlnPlot function, UMAP expressions were generated using Seurat’s FeaturePlot function, and dot plots were generated using Seurat’s DotPlot function.

inferCNV (25) was used to identify cells with high proportions of somatic mutations to identify likely cancerous cells. For each sample, we extracted ductal cells as potentially cancerous and used myeloid, epithelial, and T and NK cells as healthy reference cells. Default parameters were used for inferCNV with random forest clustering. To measure copy-number scores (CNS), we calculated the quadratic mean of the copy-number alteration counts on a per cell basis. Cells above the 95th percentile of the reference healthy cells per sample were considered putative tumor.

ReactomePA (26) was used to determine functional properties of cell types. For each cell type, differential expression was performed using the Wilcoxon rank-sum test to identify differentially expressed genes. Prior to pathway analysis, mitochondrial genes and ribosomal proteins were removed. The remaining genes that met an FDR threshold of 0.05 were included in pathway analysis. ClusterProfiler (27) was used to visualize differences in enriched pathways between cell types.

Within each finely labeled cluster group, differential gene expression analysis was performed to identify marker genes that were cell type–specific. For each cell type, we selected the top 25 differential gene markers by highest average log₂ fold change and scored four bulk RNA-seq datasets [The Cancer Genome Atlas Pancreatic Adenocarcinoma] (TCGA-PAAD), PanCan, International Cancer Genome Consortium, and GSE71729] using single-sample gene set enrichment analysis (28). Resultant scores were converted to Z-scores to allow for easier comparison between cell types. Signature scores were used as a surrogate for cell type proportions in the bulk RNA data.

All analyses were performed in R 4.0.2. We tested for associations with survival using an unadjusted Cox proportional hazards model. Additionally, overall survival (OS) times were dichotomized at the 67th percentile of scores and visualized using Kaplan–Meier curves with a log-rank test for differences. To identify coexisting cell types, pairwise Pearson correlations were calculated between cell types of interest and clustered using hierarchical clustering. Differences in cell type proportions were tested using a hypergeometric test.

Formalin-fixed paraffin-embedded sections from pathologist-selected PDAC regions were processed onto barcoded slides or sent to the University of Michigan’s Advanced Genomics Core facility to be run with the CytAssist platform for 10X Genomics Visium ST. Data were aligned to human reference panel GRCh38 using SpaceRanger V2.1.0. In Seurat V3, we removed spots with low overall expression, few total counts, and high mitochondrial percent. Spots were deconvolved with scRNA-seq atlas coarse labels using Cottrazm (29). Spots at greater than 40% purity for a ductal cells, fibroblasts, and myeloid cells were further deconvolved for finely labeled cell types.

We investigated spatial interactions using a nearest neighbor approach previously developed for other spatial technologies (30). Due to the hexagonal layout of Visium slides and the possibility of multiple subtypes per spot, we calculated the 12 nearest neighbors for each spot, limiting the distance counted only to include spots within proximity to each other. The average number of cell type neighbors per spot was then calculated as the average number of cell type A adjacent to cell type B divided by the total number of cell type A that has a neighbor of cell type B.

We performed multiplex immunofluorescence (IF) and RNAScope in situ hybridizations on formalin-fixed paraffin-embedded sections from patients with PDAC under IRB 1714-01. Briefly, sections were deparaffinized and subjected to target retrieval and hydrogen peroxide blocking as per the manufacturer’s instructions (RNAScope Multiplex Fluorescent Detection Reagents v2, ACD Cell Diagnostics). Codetection antibody diluent was used to block (ACD Cell Diagnostics, #323160), and then primary antibodies [pan cytokeratin mAb (AE1/AE3), eBioscience, #53-9003-82, 1:100 and α smooth muscle actin (αSMA), Sigma, #A2547, 1:400, S100α2, AbCam #ab247877) were incubated]overnight at 4°C. Samples were subjected to postfixation with 10% neutral buffered formalin. Then the in situ protocol was followed with probes Hs-LAMC2-C2 #501371-C2 and Hs-LY6D #484681. RNA probes were detected with Akoya Biosciences kits (Cy 3.5 #NEL76001KT and Cy 5.5 #NEL766001KT). We selected blocks with tumor content in consultation with a board-certified gastrointestinal pathologist. For quantification of ductal subtypes from multiplex staining, custom segmentation steps were applied. In brief, the 4',6-diamidino-2-phenylindole (DAPI) expression images were first thresholded to identify the highest intensity pixels using Ostu algorithm, and then individual cells were identified using the watershed algorithm (Python v3.8, skimage). The other expression channels required additional presmoothing with Gaussian and post-threshold smoothing with median filtration. The mean intensity per cell was quantified for each expression channel. For quantification of proximity of CXCL10⁺ CAFs with ductal populations, we stained for LAMC2 (ductal 4), S100A2 (ductal 7), and CEACAM5 (ductal 1) in coordination with detection of CXCL10⁺ CAFs (αSMA). We took on average 15 randomly selected regions per patient tissue section to image. We quantified what percentage of the time (of 100%) we detected CXCL10⁺ CAFs within the high-power field adjacent to ductal 1– (n = 6), ductal 4– (n = 3), or ductal 7–expressing (n = 6) regions.

scRNA-seq files are from publicly available sources (EGAS00001002543, GSE154778, GSE155698, GSE156405, GSE158356, GSE194247, GSE202051, GSE205013, GSE211644, GSE229413, phs001840.v1.p1, and PRJCA001063). Processed data used to generate figures are available at https://zenodo.org/records/14199536 or through contacting the corresponding author. Bulk RNA-seq data are from publicly available sources (TCGA-PAAD, E-MTAB-6134, GSE71729, and PACA-CA). Raw ST files are available through dbGaP with accession number PRJNA1124001. Imaging files for ST are available at https://zenodo.org/records/13379726. Human Tumor Atlas Network ST data are available through Synapse (17). PDAC cell line data were obtained from European genomic archive and Gene Expression Omnibus (PRJEB25797, PRJEB25806, and GSE235027). Organoid bulk RNA-seq data were obtained from dbGaP (phs001611). Code used to perform the analysis and generate figures is available at https://github.com/PDAC-MULTIOMIC/PDAC_Atlas.

A limitation of previous scRNA-seq studies is the number of patients included in these datasets. It is well established that generally, increasing the number of patients analyzed will decrease false positive and false negative results. To test the hypothesis that aggregating a larger cohort of patients would reveal novel biology regarding the PDAC subtypes and associated TMEs, we combined 229 pancreas scRNA-seq samples from 12 different studies containing 172 PT, 25 metastatic lesions (Met), 26 adjacent normal (AdjN), and 6 healthy donor pancreata. The metastatic biopsy samples originated from lung (n = 1), liver (n = 21), vaginal apex (n = 1), peritoneum (n = 1), and omentum (n = 1), although it should be noted that most comparisons drawn between PT and Met samples will be predominantly liver metastases. After preprocessing and strict quality control parameters (detailed in “Materials and Methods”), 726,107 high-quality cells remained, with 562,321 from PT, 56,442 from Met, 74,035 from AdjN, and 33,309 from healthy donor samples. We included a large number of samples across different published studies (Supplementary Table S1). Our harmonization approach and stringent quality control (detailed in “Materials and Methods”) ameliorated any significant batch effects between datasets (Supplementary Fig. S1A–D). Using the Louvain algorithm, we clustered the cells in an unbiased manner. Major cell types were annotated using canonical markers (Supplementary Fig. S1E), and all cell types were identified in the four disease states (Fig. 1A and B). Of note, the metastatic samples had cells that grouped within the acinar cluster, but as acinar cells are absent in metastatic sites, these cells are negative for acinar lineage markers (Supplementary Fig. S1F and S1G), suggesting a similar transcriptional phenotype but distinct cell types (such as an epithelial/secretory population).

Acinar to ductal metaplasia (ADM) is a cellular process most evident in mouse models, resulting from injury or inflammation, in which differentiated pancreatic acinar cells transdifferentiate into duct-like cells that express common duct cell markers. These duct-like cells can then acquire mutations, transitioning first to pancreatic intraepithelial neoplasia (PanIN), which, with tumor-suppressor mutations, leads to ductal adenocarcinoma (31). PDAC has been divided into two molecular subtypes: the well-differentiated classic subtype and relatively undifferentiated basal-like subtype. These subtypes coexist in single tumors, emphasizing the profound heterogeneity of primary PDAC. Some PDAC tumor cells express markers of both subtypes, suggesting that the cells exist on a spectrum of subtype states, with inherent plasticity allowing them to interconvert (17, 32–34).

To investigate putative lineages and tumor subtypes in our large-scale atlas, we separated the KRT18⁺ /KRT19⁺ (encoding keratin 18 and 19) ductal population and clustered them independently. We found several populations of duct cells enriched in the donor and AdjN samples (Fig. 1C and D), including one transcriptionally similar to the PanIN cluster identified in a previous study (ductal E on Fig. 1C and F; expressing COL1A1 and COL1A2; ref. 17) and one producing acinar genes (AMY2A; Fig. 1C–F). We also identified several populations of duct cells highly enriched in the tumor samples (ductal A–C; Fig. 1C–F). CNSs were used to determine putative tumor status in the ductal clusters using inferCNV (25) with a numerical threshold of 0.15. As expected, cells enriched in the PT and Met disease states generally had higher inferred CNSs than the ductal clusters enriched in donor and AdjN, as well as cells displaying a PanIN signature (Fig. 1G; Supplementary Fig S1H). We noticed a cluster of cells expressing acinar cell markers that grouped with the duct cells (Fig. 1C–E), suggestive of ADM. To examine this in more detail, we combined the two distinct acinar clusters together and reclustered the data (Fig. 1H). We looked at the expression of acinar, ductal, ADM, and PanIN markers (Fig. 1I). Notably, the acinar cell markers, CTRB2, AMY2A, and PRSS1, are expressed in the majority of the cells in the UMAP, whereas duct cell markers (CA2, SLC4A4, and CFTR) are expressed primarily in cluster 1. Consistent with the population showing both acinar and ductal characters, ADM markers (AQP1, MUC6, and ONECUT2) also highlight the right side of the UMAP. We conclude that cluster 3 is primarily differentiated acinar cells, with cluster 1 being likely ADM. Bridging these populations in cluster 2 are cells that express PanIN markers, along with ADM markers, with these cells possibly representing low-grade PanIN cells.

We next agnostically identified distinct transcriptional subtypes by subclustering the tumor cell enriched ductal populations (ductal clusters A, B, and C Fig. 1C). The final number of clusters was determined using unbiased clustering analysis, resulting in nine transcriptionally distinct cell populations (Fig. 2A and B; Supplementary Fig. S2A). To further support the notion that the tumor enriched clusters indeed contained malignant cells, we confirmed that the inferred CNS values were higher for the nine clusters compared with cells isolated from AdjN and donor tissues (Fig. 2C). Figure 2D shows the distribution of ductal cell subtype proportions on a per sample basis. Clusters 2 and 8 were significantly enriched in donor and AdjN samples, whereas clusters 1, 4, and 7 were enriched in PT and Met samples (hypergeometric test, P < 0.001; Supplementary Fig. S2B and S2C). Whereas some tumors were predominantly composed of a single ductal subtype, typically population 1 or 2, most tumors were heterogeneous, containing multiple ductal populations, consistent with recent findings (Fig. 2D; ref. 8). These findings highlight the intratumoral heterogeneity identified in prior studies (35–37).

Using previously published scRNA-seq–derived classic and basal gene signatures (7) we scored each ductal population from our atlas and visualized the per cell scores (Fig. 2E). Populations 9, 3, and 1 scored the highest for the classic gene signature with mean cluster scores of 1.12, 0.71, and 0.51, respectively. Populations 7 and 4 had the highest scores for the basal signature (1.52 and 0.85, respectively; Fig. 2E). Clusters 5 and 6 scored positively for both signatures, suggesting they were intermediate or “hybrid” cell types. Ductal populations 2 and 8 scored relatively low or negatively for both gene signatures. We then evaluated the gene expression signatures that have been previously associated with classic and basal molecular subtypes (Fig. 2F; refs. 7, 15). We observed enrichment for classic genes (namely, TFF1, LYZ, and CTSE; left side of UMAPs, with the exception of ductal 8, which is characterized by a normal gene signature) and enrichment for basal genes (namely, S100A2, KRT17, and KRT6A; right side of UMAPs), consistent with the signature scoring (Fig. 2E). Using ReactomePA (26), we investigated pathways enriched in each of our clusters to determine potential functional distinctions (Fig. 2G). Ductal 1 and 3 were both enriched for pathways relating to metal transporters, and ductal 3 exhibited an association with matrix metalloprotease activity (Fig. 2G). Ductal 9 was enriched for pathways involving cell metabolism (e.g., ATP synthesis, citric acid cycle, and electron transport). The gene signature for ductal 6 was consistent with the endoplasmic reticulum stress response, and ductal 5 exhibited an IFN signaling signature (Fig. 2G). Ductal 7 was enriched for pathways relating to regulation of the cell cycle (Fig. 2G). In ductal 4, we observed enrichment in pathways relating to ECM organization and c-Met signaling (Fig. 2G). We noted that vimentin expression, a marker of an epithelial to mesenchymal transition, was highest in ductal 7 compared with all other ductal populations (Supplementary Fig. S2D). TP63 is known to be a driver of basal cell state (bioRxiv 2023.10.24.563848), so we evaluated its expression in the nine ductal populations and found it to be enriched in ductal 4 and 7 (Fig. 2H).

Using the Wilcoxon rank-sum test, we identified genes that were enriched in one cluster relative to the others by highest average log fold change. We determined whether these agnostic scRNA-seq–derived ductal gene signatures (Supplementary Table S2) correlated with survival of patients with PDAC from four independent bulk RNA-seq datasets containing PT samples (15, 16, 22, 23). As expected, both ductal 4 and 7 basal subpopulation gene signatures consistently correlated with worse OS using an unadjusted Cox proportional hazards model (Fig. 2I). Kaplan–Meier curves using the 67th percentile to dichotomize high versus low scores were used to further investigate this relationship (Supplementary Fig. S2E). Similarly, ductal 5 consistently correlated with worse OS, although none of the studies reached statistical significance. Additionally, ductal 1, 8, and 9 consistently correlated with better OS in the Puleo dataset, whereas ductal 1 correlated in TCGA-PAAD. The remaining populations showed mixed impact on the survival analysis, with only ductal 2 correlating with better survival in TCGA-PAAD (Fig. 2I). Overall, our findings are consistent with previously published studies showing that basal–enriched samples have worse OS, whereas classic enriched samples have better OS.

In the PDAC TME, CAFs and their associated ECM make up a large portion of the tumor mass (38). Interactions between tumor cells and CAFs can either restrict or promote tumor growth depending on the timing, context, and signaling pathways active in both cell types (9–11, 39–41). Three transcriptionally distinct fibroblastic populations have been described in PDAC, namely, iCAFs, myCAFs, and antigen-presenting CAFs that seem to have distinct functions and spatial distributions (42). To explore the relationships between subtypes of tumor cells and CAFs using our atlas, we separated fibroblasts (DCN⁺COL1A1⁺LUM⁺PDPN⁺) and clustered them alone, identifying eight transcriptionally distinct clusters (Fig. 3A and B; Supplementary Fig. S3A). Then, we evaluated whether those clusters were associated with an iCAF (defined by expression of C7 and DPT) or myCAF identity (defined by expression of ACTA2 and COL11A1; Fig. 3C; refs. 9, 11). We did not observe a transcriptionally unique antigen-presenting CAF population in our dataset. We identified five populations of myofibroblasts and three populations of inflammatory fibroblasts that were transcriptionally distinct (Fig. 3D). We labeled these populations as myFb and iFb, instead of myCAF and iCAF, respectively, as our atlas included nonmalignant pancreas tissues.

To gain insight into the potential functional relevance of these Fibroblast (Fb) clusters, we performed pathway analysis on the uniquely enriched genes for each population (Fig. 3E). As expected, inflammatory fibroblasts (iFb) 1 and iFb 2 were enriched for complement signaling. The iFb 3 population showed an association with both complement and ECM pathways and furthermore uniquely expressed high levels of the potent regulator of angiogenesis, VEGFA (Fig. 3D and E). Broadly, myofibroblast (myFb) 1 to 5 populations were highly enriched for pathways related to ECM signaling and collagen assembly (Fig. 3E). MMP11 and COL11A1 were highest expressed in myFb 1 compared with the other Fb populations (Fig. 3D). The myFb 2 population was enriched for pathways relating to ATP synthesis, citric acid cycle, and hypoxia signaling, suggesting these cells could be the previously described metabolic cancer associated fibroblasts (“meCAF”) population which was shown to correlate with increased risk of metastasis in PDAC but better responses to immunotherapy (Fig. 3E; ref. 43). CXCL14, a chemokine previously implicated in promoting tumor cell–invasive properties (44), was enriched in myFb 3 compared with other myFbs and iFbs (Fig. 3D). We identified a unique population of CXCL10 expressing myFbs (myFb 5) highly enriched for IFNγ genes, consistent with the notion that expression of this chemokine is downstream of IFNγ signaling (45). We noted expression of important chemokines/cytokines previously implicated in promoting tumor growth and immune-suppressive signaling in PDAC (46–50), including CXCL1, CXCL2 CXCL10, and CXCL14, throughout both fibroblast populations at varying levels, suggesting the majority of fibroblasts can express cytokines independent of classification (Fig. 3F).

We investigated the association of the fibroblast scores (Supplementary Table S3) with OS in the bulk RNA-seq datasets, as we did with the tumor clusters. We found that in general, the myFb clusters were associated with worse OS (Supplementary Fig. S3B). Overall, our scRNA-seq atlas reveals an enrichment in what would have been typically termed “inflammatory” fibroblasts in normal/adjacent pancreas while myofibroblasts expand during tumorigenesis (Fig. 3G; Supplementary Fig. S3A). However, data mined from our atlas suggest significant overlap in immune-modulatory functions (based on expression of cytokines and chemokines) shared between iFb and myFb populations. Furthermore, we find evidence that higher expression of a myofibroblast gene signature may correlate with worse overall outcomes in PDAC, consistent with previous reports (43, 51).

Myeloid cells are major cellular components of the PDAC TME and regulate tumor progression by communicating with tumor cells, CAFs, and antitumor T and NK cells (52). We identified several major myeloid lineages: neutrophils, dendritic cells (DC), macrophages, and monocytes (Fig. 3H), as classified according to known canonical markers (Fig. 3I). We compared the proportion of myeloid cells in each patient grouped by disease state (Fig. 3J; Supplementary Fig. S3C). Macrophages were the most commonly observed myeloid population, in line with previous studies (53, 54). Myeloid cell population scores (Supplementary Table S4) were evaluated for associations with OS in four bulk RNA-seq datasets, and we found that DCs consistently associated with better OS and neutrophils in general associated with worse OS (Supplementary Fig. S3D). Conversely, monocyte and macrophage scores did not have a clear correlation with survival. We concluded this may be due to antigen-presenting DCs priming more effective T-cell responses, whereas neutrophils are highly immune-suppressive in PDAC (immature neutrophils or myeloid-derived suppressor cells; refs. 50, 52, 55–57).

Given the worse OS observed in patients with PDAC with basal-enriched tumors, we investigated whether there were consistent components of the TME that correlated with basal tumor cells, and thus, with OS. For the putative tumor cells, fibroblasts, and myeloid cells, we generated gene signatures using the top 25 differentially expressed genes per cell type (mean log₂ fold change). Using single-sample gene set enrichment analysis (ref. 58), we scored the four PDAC bulk RNA-seq datasets to estimate cellular composition. To evaluate whether any cell types consistently coexist within patient tumors, we performed a pairwise Pearson correlation for all cell types followed by hierarchical clustering (Fig. 4A). Strikingly, three distinct clusters of commonly associated cell types emerged using the silhouette method for cluster determination. Of note, ductal 2 consistently correlated with ductal 1, 3, and 9, suggesting that it is likely a classic subtype as well. Ductal 4 and 7 were always detected in the same cluster, in line with our earlier scRNA-seq unbiased clustering analysis (Fig. 2E). Intriguingly, ductal 5, one of the two hybrid clusters, was also associated with ductal 4/7, reinforcing the possibility that this population is in transition to or from a basal state (Fig. 4A). Ductal 4 and 7 consistently clustered with myofibroblasts, specifically myFb 5 in all four bulk datasets. myFb 1, myFb 3, and myFb 2 clustered with both basal ductal populations in three of the four bulk datasets analyzed. myFb 5s were distinguished from the other Fb populations in part because of its expression of the chemokine CXCL10. We verified the correlations between our basal clusters and the CXCL10 expressing myofibroblast clusters in the scRNA-seq data (Fig. 4B). We found that the proportion of myFb 5s were positively correlated with both ductal 4 and 7 proportions (r = 0.247 and 0.517, respectively) and negatively correlated or lacking correlation with classic populations (ductal 2 r = −0.211 and ductal 9 r = −0.07). CXCL10 is a chemokine that has been linked to the regulation of immune cells (59), so we investigated the expression of its receptor CXCR3 in our TME populations. We found that in addition to expression in myeloid and T and NK cells, CXCR3 was expressed by tumor cells, with expression in all tumor subclusters apart from ductal 8 (Fig. 4C). Additionally, IFNγ, which regulates CXCL10 expression (60), was highly expressed in T and NK cells, whereas its receptors IFNγR1 and IFNγR2 were both highly expressed in fibroblasts (Fig. 4C). To determine whether there were compositional differences across studies and samples, gene signature scores were then minimum–maximum normalized and scaled so each sample’s scores sum to one (Fig. 4D). Importantly, the four studies show remarkably similar distributions of the different cell types, with some differences by study, (e.g., International Cancer Genome Consortium and TCGA have higher iFb 3, whereas Puleo has higher myFb 1). Overall, these data demonstrate the power of this multitranscriptomics approach to identify niches of putative cellular cross-talk that correlate to survival.

To determine which hypothetical cellular interactions were supported by spatial analysis in the highly complex and heterogenous TME of PDAC, we deconvolved 22 ST PDAC patient tumor sections [Human Tumor Atlas Network (HTAN) (17), and our own Henry Ford Pancreatic Cancer Center (HFPCC) cohorts; Supplementary Table S5]. Pathologist-annotated PDAC regions were selected for ST, and samples included both chemotherapy/radiation-treated and treatment-naïve PDAC. Using the Cottrazm pipeline, which maps putative malignant regions (using inferCNV) and overlays this with cell type–specific scRNA-seq–derived gene signature deconvolution (29), we mapped ductal, myeloid, and fibroblasts (the three most abundant components in the PDAC TME) subtypes to the Visium ST platform (spots are at 55-micron resolution, containing 1–30 cells; ref. 61). We identified regions with each cell type identified in the scRNA-seq data in all the spatial PDAC tissues (Fig. 5A), limiting concerns of the influence of dissociation-induced stress on our scRNA-seq gene signatures. Interestingly, the ductal 8 gene signature was found in the low- to high-grade PanIN region (Fig. 5B), coexisting with the ductal 1 gene program, suggesting retention of a more normal ductal phenotype along with upregulation of classic signature. Within the same patient tumor, we noted a well-differentiated tumor region (R1), which contains classic enriched ductal types 1 and 3 near two poorly differentiated tumor regions (R2 and R3) containing ductal types 4 and 7 (Fig. 5C). These data highlight both intratumoral heterogeneity as well as validation that our scRNA-seq–identified ductal populations are found in intact tissue. Contrary to our scRNA-seq analysis, in which ductal 1 was predominant, we found that, generally, ductal 4 was the most predominant in situ subtype whereas ductal 7 was rare (Fig. 5D). We quantified CAF and myeloid proportions separately per sample (Fig. 5E and F), with both showing high levels of intrapatient heterogeneity. T and NK cell content was evaluated in the spatial samples as well, but using current deconvolution methods, very few cells were detected, which may be a technical limitation of ST but corroborated the limited T and NK cell infiltration noted in PDAC tumors (Supplementary Fig. S4A; ref. 62).

Using our deconvolved ST analysis, we next interrogated whether the coenrichment of ductal 4/7 with myFBs observed in the bulk RNA-seq data was reflected by spatial proximity. We used a nearest neighbor approach previously developed for other spatial technologies (30, 63) to identify potentially interacting cells. Briefly, for each sample, neighbors of each ductal cell type were counted. We then randomly permuted cell labels for all cells 100 times and counted the number of neighbors for each ductal cell type, providing a distribution of counts under spatial randomness. We compared the true observed number of neighbors with this null distribution to derive empirical P values to determine whether cell types were consistently grouping together. Using a binomial test, we evaluated whether there were a consistent enrichment of spatial associations between myFB 5 and the ductal subtypes across samples. We identified basal-like ductal 4 as significantly strongly associated with myFbs, in particular those highly expressing the metalloproteinase MMP11 and the cytokine CXCL10, averaging more than two myFbs per ductal 4 cell (n = 13/15; P = 0.007; Fig. 5G), supporting our bulk RNA-seq correlation results (Fig. 4A). Additionally, we performed nearest neighbor testing between classic ductal cells and fibroblast subtypes (Fig. 5G) and saw that there was less proximity to fibroblast cells in general, with myFb 5 averaging one neighbor per ductal cell (n = 6/15; P = 0.607).

We used CIBERSORTx (64) to deconvolve the four bulk RNA-seq datasets to determine cellular content of the major labeled cell types (Supplementary Fig. S4B). As with the ST samples, very few T and NK cells were identified. Given the lack of reliability of interrogating T and NK cells from bulk RNA-seq and ST approaches, we turned to our scRNA-seq atlas. Previous work has demonstrated a positive association between infiltration of CD8⁺ T cells and improved OS in PDAC (60). We utilized our scRNA-seq atlas to determine the association between ductal 4 and 7 contents with abundance of CD8⁺ and CD4⁺ T cells (Supplementary Fig. S4C). In samples with higher ductal 4 content, we saw no significant association with CD4⁺ T cells; however, there was a significant correlation between ductal 7–high samples and CD4⁺ T cells. Interestingly, we found a negative correlation between high ductal 4/7 content and the presence of CD8⁺ T cells, suggesting reduced infiltration of potentially cytotoxic T cells in basal-enriched tumors, possibly partly explaining the poorer outcome in these patients. Thus, using our newly built scRNA-seq atlas and integrating with ST analysis, we uncovered a novel interaction between CXCL10-expressing myFbs and basal-like ductal cells (ductal 4) in PDAC tumors, highlighting the potential of using a large atlas and integrating with ST.

We wanted to test how representative 2D and 3D cell cultures are compared with in situ tumors, so we utilized bulk RNA-seq datasets from published in vitro patient-derived PDAC cell and organoid lines and applied CIBERSORTx (64) deconvolution to estimate the proportion of ductal clusters. We found that cluster 4 was highly represented in all 2D cell lines (Fig. 5H). Matsumoto and colleagues (65) found that treatment with an ERK inhibitor improved drug responses against gemcitabine in basal subtype patient-derived organoids. Our deconvolution showed that ERK inhibition or KRAS silencing in KRAS-mutant PDAC cell lines (66) decreased ductal 4 content and increased ductal 3 and 5 content, a transition from a basal-like to a classic phenotype (Fig. 5H) in all lines except PANC-1 and SW1990. Notably, both PANC-1 and SW1990 cells are ERK- and KRAS-independent (67, 68). We further investigated the maintenance of these ductal cell populations in two additional cell line datasets. As before, ductal 4 was the most well represented (Supplementary Fig. S4D and S4E). We noticed that in both datasets from Bryant and colleagues (66), the deconvolution results matched well, whereas in Zhou and colleagues (69), the ductal subtype composition differed in the nonbasal components, suggesting differences in culture, rather than an artifact of the computational method, consistent with findings that show that media formulation affects the transcriptional states of cell lines and organoids (7). We further evaluated a patient-derived organoid bulk RNA-seq dataset and identified one sample with ductal 9 and one with ductal 7, neither of which were present in the 2D cell line data (Fig. 5I), possibly highlighting differences between 2D and 3D cell cultures. Additionally, we identified ductal 8 consistently in the normal samples, reflecting the enrichment in the donor and AdjN samples. Similar to the cell line data, we observed that most organoid lines (70) were enriched for ductal 4, including the normal samples (names highlighted in blue; Fig. 5I). Overall, we noted that these data were consistent with the spatial deconvolution, which demonstrated that ductal 4 was the predominant subtype (Fig. 5D). Additionally, we used Visium ST data from Bell and colleagues (71) from AdjN samples, having computationally annotated cell types, to help delineate the differences in our basal and classic subtypes. Using Seurat’s AddModuleScore, six samples were scored for each of our nine ductal signatures. We saw enrichment for ductal 1, 2, 3, and 9 in PanIN regions, which they suggest are more similar to classic cells [in agreement with Bell and colleagues (71)], confirming that ductal 2 is likely classic (Supplementary Fig. S5). Furthermore, we saw an enrichment of ductal 8 in what they deemed normal ductal cells and low-grade PanIN (consistent with the enrichment we observed in the tumor ST data) with a degradation in high-grade PanIN, suggesting that ductal 8 could be a cellular component of both reactive duct cells and early PanIN (Supplementary Fig. S6).

Finally, to examine ductal subtypes in PDAC patient samples, we used IF to stain tissues for PANCK to locate ductal/epithelial cells and costained for markers associated with either ductal 4 (LAMC2) or ductal 7 (S100A2; Fig. 6A–E; Supplementary Fig. S7). Consistent with the scRNA-seq, ST, and bulk RNA-seq data, we often observed a high prevalence of the ductal 4 subtype (Fig. 6C–E) and less commonly ductal 7 subtype cells, although there were patients who were S100A2-enriched overall (Fig. 6C and D). We observed that within the same region, there was a stark difference in expressions of LAMC2 and S100A2 between normal ductal epithelium compared with tumor cells, as well as a shift in cell morphology to a more mesenchymal spindle shape suggestive of an epithelial to mesenchymal transition/basal state (representative image shown in Fig. 6C). Similar to what we see in the UMAP visualizations from the scRNA-seq data (Fig. 2A), we observed transitionary states, with varying degrees of single-marker expression of LAMC2 or S100A2 in many cells and coexpression of both markers in some cells (Fig. 6B; Supplementary Fig. S7). These data emphasize the high level of cellular plasticity between classic and basal populations but also distinct transcriptional cell states within the basal subtype.

To further investigate whether CXCL10-expressing myFbs exist in proximity to basal cells and less frequently in proximity to classic cells, as would be predicted from our transcriptomic data, we stained PDAC tissues for CEACAM5 to identify ductal 1 (classic), LAMC2 mRNA to identify ductal 4 (basal), S100A2 to identify ductal 7 (basal), and αSMA protein to identify myFbs in combination with an in situ hybridization probe against human CXCL10 (Supplementary Fig. S8A–S8C). We noticed a lack of proximity between ductal 1 and CXCL10⁺ fibroblasts, confirming what we saw in the ST, scRNA-seq, and bulk (Supplementary Fig. S8D). We observed ductal 4 and ductal 7, but not ductal 1, regions nearby CXCL10-expressing αSMA⁺ CAFs, now at individual cell resolution (Supplementary Fig. S8B). Interestingly, we found that regions with high myFb and ductal 4 contents also showed expression of CXCL10 within the tumor cells, reinforcing a potential for local TME responses to IFNγ signaling (Supplementary Fig. S8B).

Existing PDAC subtyping studies have been limited to analysis of bulk RNA-seq, which lacks granularity, or to scRNA-seq data that have been relatively underpowered due to limited sample sizes. By combining available scRNA-seq studies, we have substantially increased the statistical power needed to more accurately identify tumor subtypes and applied our findings to bulk RNA-seq data to determine what compartments of the PDAC TME are associated with these subtypes. With this approach, we have assigned duct cells as malignant by inferred copy-number variation. When these cells were analyzed in isolation, four classic (ductal 1, 2, 3, and 9), two intermediate/hybrid (ductal 4 and 5), one secretory (ductal 8, highly enriched in donor and AdjN pancreata), and two basal (ductal 4 and 7) populations were identified, of which only 4 and 7 correlate with poor survival. Although we examine the ductal clusters in finer granularity, given the larger sample size, these ductal gene signatures are in general alignment with previous efforts (7, 8, 33, 72, 73). We also leveraged scRNA-seq TME gene signatures combined with the ductal markers to deconvolve bulk RNA-seq datasets and found consistently enriched cell compartments in 744 patients with PDAC. Using this approach, we identified that the two basal and an intermediate population were consistently coenriched with myofibroblasts expressing the chemokine CXCL10. We validated this observation with ST and immunostaining to confirm this association in intact tissue. This finding demonstrates the powerful resource our atlas represents; when integrated with ST, it revealed a new putative signaling niche within human tumors between CXCL10⁺ fibroblasts and basal ductal cells.

Bulk RNA-seq studies are limited in that they struggle to accurately capture the high level of cellular heterogeneity that is a prime characteristic of PDAC. Moffit and colleagues (15) identified general classic and basal subtypes in PDAC using bulk RNA-seq data, associating with better and worse OS, respectively. In concordance with our results, they observed that cell lines are predominantly of the basal subtype; however, with our improved granularity, we identified a mixture of classic and intermediate subtypes in commonly used PDAC cell lines grown in 2D. Bailey and colleagues (32) identified four distinct subtypes: squamous, ADEX, pancreatic progenitor, and immunogenic, identifying additional subtypes, although this study was still limited by using bulk RNA-seq, and some of these subtypes may be influenced by the “contamination” of adjacent nontumor tissue. With the advent of scRNA-seq technology, Chan-Seng-Yue and colleagues (33) were able to push beyond the previously defined classifications, separating classic and basal into classic 1, classic 2, basal 1, and basal 2, respectively. In agreement with this prior study, our study also detected two unique basal populations, one of which was highly represented in nondissociated tissues. Raghavan and colleagues (7) further sought to understand the role of the TME in the differentiation of tumor cell states. They identified an intermediate tumor subtype and showed how the three tumor subtypes (classic, intermediate, and basal) vary in their associated TME. Whereas the general lack of cellular differentiation in basal-like tumors is presumed to be, in large part, the reason for increased aggressiveness of the tumors, in our study, the highest basal-scoring cells were negatively correlated with the presence of CD8⁺ T cells. Furthermore, our atlas reinforces the existence of two distinct intermediate/hybrid populations, one of which (ductal 5) is highly enriched with the two basal populations. Hwang and colleagues utilized snRNA-seq, to try and decouple the relationship between treated and untreated PDAC samples and their transcriptional subtype differences, identifying seven subtypes (the most distinct subtypes yet), some of which overlapped with previous subtypes; however, the sample size was still relatively underpowered. Williams and colleagues (74) further showed that the intratumoral heterogeneity identified in many of these previous studies and our own have significant association with patient outcomes.

Multitranscriptomic approaches have great power in defining the biology of numerous cancers. However, we found it important to validate our results with multiple technical methods in order to draw robust and reproducible conclusions. For instance, the consistent association between CXCL10-expressing CAFs and ductal 4 tumor cells suggested a possible physical association between these cell types, which we confirmed using ST and IF approaches. The CXCL10 receptor, CXCR3, is expressed on subsets of DCs, T and NK cells, myeloid cells, and tumor cells themselves. Future work is required to determine the functional relevance of this interaction. Interestingly, basal-enriched tumors are associated with fewer CD8⁺ T cells, suggesting a higher degree of immune suppression in the TME, an observation that may have implications with immune therapy approaches.

On a study of this scale, combining data from 12 published studies, the potential for unmanageable batch effect is high; thus, the choice of batch correction algorithm is critical to ensure the findings are meaningful. However, batch correction has the potential to introduce artifacts and remove biological signal. Many groups have compared existing batch correction algorithms for speed and removal of technical artifacts while maintaining biological relevance (20, 75). With these concerns in mind, we thoroughly considered the implications our chosen batch correction method would have on the interpretation of our results. To minimize the effects of batch correction, we chose to include only studies that were performed on the 10X Genomics scRNA-seq platform for which raw sequencing files were made available, allowing us to properly align the data with the same software and reference genome, removing two potential sources of batch effects. We compared batch correction methods, including rPCA, Seurat3 (18), and scVI (76). We elected for Harmony (24) in large part because other approaches continued to treat each epithelial cell from each sample as a separate cluster, preventing us from drawing the meaningful inference required to inform population-based interpretations. It is because of this that we also chose to treat each sample as a batch, rather than each study or disease state, although we compared the interpretation of our findings using study-based batch correction which yielded similar results. This may have limited our ability to discriminate rare PDAC subtypes. We were able to identify meaningful results that we could validate using a variety of independent methods, including mining publicly available bulk RNA-seq data, in situ hybridization, and ST.

We envision using the atlas signatures applied to cell lines to assess how various treatments cause changes in the ductal composition of tumors, allowing us in future studies to construct a predictive algorithm for patient treatment. With the identification of novel duct cell populations, deeper examination of these in primary versus metastatic stages or after therapy is needed. Regarding ductal 7, we noted that this population is extremely rare in in vitro models despite it being detected more commonly in scRNA-seq, bulk RNA-seq, and ST patient data. We do not know the lineage fate of these cells, but investigation of the relationships between these transcriptional states, and how this information can be leveraged to predict drug responses and resistance to targeted therapies, is warranted.

Overall, we show that our newly built scRNA-seq atlas containing data of more than 190 patients with pancreatic cancer and 700,000 cells, integrated with ST, is a powerful resource for the pancreatic cancer field. Although it corroborates several prior observations, it has also uncovered a novel relationship between CXCL10⁺ myofibroblasts and basal tumor cells, as one specific example of how it can be deployed. That it predicts validated spatial relationships between cell types and potential signaling niches that have not been identified in other studies reinforces the power of multiomic approaches in understanding this deadly disease.

B.Z. Stanger reports grants from Revolution Medicines and Boehringer Ingelheim and other support from iTeos Therapeutics outside the submitted work. A.M. Waters reports grants from the NCI during the conduct of the study and personal fees and nonfinancial support from Revolution Medicines outside the submitted work. E.J. Fertig reports grants from the NIH/NCI, the Lustgarten Foundation, and Break Through Cancer during the conduct of the study and grants from the NIH/NIA and Roche/Genetech and personal fees from Viosera/Resistance Bio, Mestag Therapeutics, and Merck outside the submitted work. No disclosures were reported by the other authors.

I.M. Loveless: Data curation, formal analysis, visualization, writing–original draft, writing–review and editing. S.B. Kemp: Formal analysis, writing–review and editing. K.M. Hartway: Formal analysis, writing–review and editing. J.T. Mitchell: Validation, writing–review and editing. Y. Wu: Formal analysis, writing–review and editing. S.D. Zwernik: Data curation, writing–review and editing. D.J. Salas-Escabillas: Writing–review and editing. S. Brender: Validation, visualization. M. George: Writing–review and editing. Y. Makinwa: Data curation, writing–review and editing. T. Stockdale: Data curation, writing–review and editing. K. Gartrelle: Writing–review and editing. R.G. Reddy: Writing–review and editing. D.W. Long: Writing–review and editing. A. Wombwell: Writing–review and editing. J.M. Clark: Writing–review and editing. A.M. Levin: Methodology, writing–review and editing. D. Kwon: Data curation, writing–review and editing. L. Huang: Writing–review and editing. R. Francescone: Writing–review and editing. D.B. Vendramini-Costa: Writing–review and editing. B.Z. Stanger: Writing–review and editing. A. Alessio: Methodology, writing–review and editing. A.M. Waters: Writing–review and editing. Y. Cui: Writing–review and editing. E.J. Fertig: Supervision, validation. L.T. Kagohara: Supervision, validation, writing–review and editing. B. Theisen: Data curation, writing–review and editing. H.C. Crawford: Conceptualization, supervision, writing–original draft, writing–review and editing. N.G. Steele: Conceptualization, resources, data curation, supervision, funding acquisition, writing–original draft, writing–review and editing.

We authentically and emphatically thank all patients who chose to consent to our study and donate their resected tissue for our research. We also thank our clinical research coordinator team in help with collection and processing of tissues.

Work in the N.S. laboratory was supported by NIH/NCI grant R00CA263154 and funding from the Sky Foundation. A.M. Waters was supported by NIH/NCI grant K22CA276632 and by the Steven Goldman Memorial Pancreatic Cancer Research Fund at the University of Cincinnati Cancer Center. S.B. Kemp was supported by a Cancer Research Institute Irvington Fellowship. Y. Cui was supported by American Heart Association Award 24TPA1288424. During the time of this study, D.J. Salas-Escabillas was supported by NIH/NCI grant F31CA271576 and the University of Michigan’s Rackham Merit Fellowship. D.B. Vendramini-Costa was supported by DOD PA220131P1. R. Francescone was supported by Pancreatic Cancer Action Network Career Development Award in memory of Skip Viragh 21-20-FRAN. H.C. Crawford is supported by NIH/NCI grants R01CA247516, 1R01CA281198, R01CA244931, R01CA248160, R37 CA23742101, R01 DK128102, U01 CA200468, and U01 CA274154 and funding from the Sky Foundation, the Lustgarten Foundation, and Henry Ford Health and Michigan State University Health Sciences Center.

This study was conducted using data provided by the Ontario Institute for Cancer Research which is funded by the Government of Ontario. The views expressed in the publication are those of the author and do not necessarily reflect those of the Ontario Institute for Cancer Research or the Government of Ontario.

This study utilizes data deposited in the European Nucleotide Archive at EMBL-EBI under accession number PRJEB25797 (https://www.ebi.ac.uk/ena/browser/view/PRJEB25797) and accession number PRJEB25806 (https://www.ebi.ac.uk/ena/browser/view/ PRJEB25806).

The data from the following studies were downloaded from dbGaP:

•  

  phs002071: Multimodal Mapping of the Immune Landscape in Human Pancreatic Cancer. This project was supported by NIH/NCI grants R01CA151588, R01CA198074, and U01CA224145 and the American Cancer Society. This work was also supported by the University of Michigan Cancer Center Support Grant (P30CA046592), including an Administrative Supplement.

•  

  phs001840: Single-Cell Analysis of Human Pancreatic Ductal Adenocarcinoma. For additional information, see Elyada and colleagues (9). This work was supported by the Lustgarten Foundation, in which D.A. Tuveson is a distinguished scholar and Director of the Lustgarten Foundation Dedicated Laboratory of Pancreatic Cancer Research. D.A. Tuveson is also supported by the Cold Spring Harbor Laboratory Association, the David Rubinstein Center for Pancreatic Cancer Research at Memorial Sloan Kettering Cancer Center, the V Foundation, the Thompson Foundation, and the Simons Foundation (552716). In addition, this work was supported by NIH grants P30CA045508, P50CA101955, P20CA192996, U10CA180944, U01CA168409, U01CA210240, R33CA206949, R01CA188134, and R01CA190092 to D.A. Tuveson and U01CA224013 to P. Robson and D.A. Tuveson. P. Robson is also supported by JAX Laboratory startup funds and the JAX Cancer Center Support Grant P30CA034196-30. A. Califano is supported by the NCI Outstanding Investigator Award (R35 CA197745), the NCI Cancer Systems Biology Consortium (U54 CA209997), and two NIH Shared Instrumentation Grants (S10 OD012351 and S10 OD0217640). A. Califano is also supported in part through the NCI Cancer Center Support Grant (P30 CA013696). E.M. Jaffee is supported by the James W. and Frances Gibson McGlothlin Foundation, the Skip Viragh Center for Pancreatic Cancer at Johns Hopkins, the Bloomberg∼Kimmel Institute for Cancer Immunotherapy at Johns Hopkins, the NCI (R01 CA18492-04 and R01 CA19729603), and a Stand Up To Cancer-Lustgarten Foundation Pancreatic Cancer Convergence Translational Research Grant (SU2C-AACR-DT14-14). C.L. Wolfgang is supported by the NIH (R01CA182076). R.A. Burkhart is supported by a Stand Up To Cancer-Lustgarten Foundation Pancreatic Cancer Interception Translational Cancer Research Grant (SU2C-AACR-DT26-17). Stand Up to Cancer is a division of the Entertainment Industry Foundation administered by the American Association for Cancer Research, the Scientific Partner of SU2C. We also thank the Cold Spring Harbor Laboratory shared resources: P. Moody and C. Kanzler at the flow cytometry facility, the animal facility, the histology core, the single-cell and bioinformatics cores, and the microscopy core. These shared resources are funded by the NIH Cancer Center Support Grant P30CA045508. In addition, we are grateful for support from the following: The Cold Spring Harbor Laboratory and Northwell Health Affiliation (for J. Preall and D.A. Tuveson), the Human Frontiers Science Program (LT000403/2014-L for E. Elyada and LT000195/2015-L for G. Biffi), European Molecular Biology Organization (EMBO) (ALTF 1203-2014 for G. Biffi), the NCI (R01 CA202762 for K.H. Yu and 5T32CA148056 and 5K99CA204725 for D.D. Engle), and the NIH (R50 CA211506-02 for Y. Park).

•  

  phs001611: Pancreas Cancer Organoid Profiling. We would like to thank the patients, families, and clinicians who contributed to this work. This work was supported by the Lustgarten Foundation. D.A. T. is a distinguished scholar of the Lustgarten Foundation and Director of the Lustgarten Foundation–designated Laboratory of Pancreatic Cancer Research. This project was also supported in part by a V Foundation Translational Grant to K.H.Y and D.A.T. We also acknowledge the Cold Spring Harbor Laboratory Tissue Culture and Next-Generation Sequencing Shared Resources, which are funded by the NIH Cancer Center Support Grant 5P30CA045508. D.A.T. is also supported by the Cold Spring Harbor Laboratory Association and the NIH (NIH 5P30CA45508-29, 5P50CA101955-07, P20CA192996-03, 1U10CA180944-04, 5U01CA168409-5, 1R01CA188134-01, and 1R01CA190092-04); 5T32CA148056 and 5K99CA204725 for D.D.E. In addition, we are grateful for the following support: Stand Up to Cancer/AACRPS09 (D.A.T.), Precision Medicine Research Associates (D.A.T.), SWOG ITSC 5U10CA180944-04 (D.A.T., H.T., and E.J.K.), NCI 5P01CA013106 (C.R.V), Pancreatic Cancer Action Network-AACR 16-20-25-VAKO (C.R.V.), State of New York C150158 (T.D.D.S.; The Opinions, results, findings and/or interpretation of data contained herein are the responsibility of the authors and do not necessarily represent the opinions, interpretations or policy of New York State), Concetta Greenberg in memory of Marvin S. Greenberg, M.D (J.R.B. and J.M.W.), 1R01CA212600-01 (J.R.B. and J.M.W.), a RAN grant from the American Association for Cancer Research Pancreatic Cancer Action Network (J.R.B. and J.M.W.), ASGE Endoscopic Research Award 71040 (J.M.B), Simons Foundation Award 415604 (E.L.), NCI P20 CA192994 (E.L.), R01CA202762 (K.H.Y.), P50CA127297 (UNMC RAP Program), U01CA210240 (UNMC RAP Program, D.A.T.), and 5R50CA211462 (University of Nebraska Medical Center's (UNMC) Research Apprenticeship Program (RAP) Program). We acknowledge the contributions of Production Sequencing, Genome Sequencing Informatics, and Tissue Portal (Diagnostic Development) at the Ontario Institute for Cancer Research, as well as the University Health Network BioBank. The pancreatic cancer patient study (International Cancer Genome Consortium-CA, COMPASS) was conducted with the support of the Ontario Institute for Cancer Research (PanCuRx Translational Research Initiative) through funding provided by the Government of Ontario. S.G. is supported by the Ontario Institute for Cancer Research, a Canadian Cancer Society Research Institute grant (702316), charitable donations from the Canadian Friends of the Hebrew University (Alex U. Soyka) and the Pancreatic Cancer Canada Foundation, and the Lebovich Chair in Hepatobiliary/Pancreatic Surgical Oncology.

We acknowledge and thank the Henry Ford Translational and Clinical Research Center Biorepository for providing resources and services in support of the research reported in this publication.