Virtual microdissection identifies distinct tumor- and stroma-specific subtypes of pancreatic ductal adenocarcinoma

Virtual microdissection of PDAC

We used NMF to analyze gene expression in a cohort of microarray data from 145 primary and 61 metastatic PDAC tumors, 17 cell lines, 46 pancreas and 88 distant site adjacent normal samples using Agilent (Agilent Technologies) human whole genome 4x44K DNA microarrays (106 primary tumors were previously used in a separate bulk analysis of gene expression (GSE21501¹⁶). To validate our findings, RNA sequencing was performed on 15 primary tumors, 37 pancreatic cancer patient-derived xenografts (PDX), 3 cell lines, and 6 cancer associated fibroblast (CAF) lines derived from deidentified patients with pancreatic cancer. Histology of all available samples was reviewed by a single blinded pathologist (KEV). Table 1 summarizes the demographic and clinical characteristics of patients in our cohorts.

NMF distinguishes normal and tumor compartments

A key obstacle in the analysis of gene expression data, particularly in PDAC, is the removal of confounding normal or stroma gene expression from local and distant organ sites. Supplementary Figure 1 shows example histology of samples with both tumor, normal, and stromal tissue. We used NMF to identify gene expression which we attribute to normal pancreas, liver, lung, muscle, and immune tissues. Expression of exemplar genes from these factors, i.e. genes with distinctly large weights in a single column of G, as well as factor weights for the samples, i.e. rows of S, showed excellent agreement with known tissue labels (Fig. 1b, c, Supplementary Fig. 2). Investigation of the exemplar genes from these factors further confirmed their role as confounding normal tissue. For example, using the Kolmogorov–Smirnov test, the top-weighted genes from the liver factor show significant (p<10^-10) enrichment in the MSigDB term SU_LIVER (Supplementary table 1), and the highest weighted gene, fibrinogen beta (FGB) (Supplementary table 2), is specifically expressed in normal human liver tissue.

In addition to normal tissue from distant organs, we found two factors which were exclusive to pancreas tissue, but were differentiated from each other by their respective gene lists. One factor described endocrine function including expression of glucagon and insulin (GCG and INS) (Supplementary table 2), while the other factor described exocrine function including expression of digestive enzyme genes such as pancreatic lipase, PNLIP (Supplementary table 2). This unsupervised discovery of two molecularly distinct yet highly co-localized factors related to normal pancreatic function represents an important proof of concept in the use of NMF to identify novel features without pre-defined expression knowledge.

To validate our normal expression signatures, all available samples were reviewed by a single pathologist (KEV) to independently assess the amount of tumor, normal, and stroma cellularity. We found that many factor weights were correlated or anti-correlated to tumor cellularity (Supplementary Fig. 3). Among normal and metastatic liver samples, for example, tumor-specific factor weights were correlated with cellularity whereas the normal-specific liver factor weight was inversely related to the tumor content of a sample (Fig. 1d). These findings support our hypothesis that factor weights obtained from NMF are quantitatively indicative of underlying sample composition.

Identification of stroma-specific subtypes

Stroma is particularly important in PDAC. According to pathology assessments, stroma varies (Supplementary Fig. 1 c-e), and comprises on average 48% of our primary tumor samples with a standard deviation of 30%. Our analysis identified two factors, which describe gene expression from the stroma, which were distinctly different from the normal factors shown in Figure 1. Consensus clustering on exemplar genes from these two stroma factors divided tumor samples into two stromal subtypes, which we classified as “normal”, and “activated” (Fig. 2a). Patients with samples with an activated stroma subtype had worse median survival of 15 months and 60% 1-year survival when compared to patients with a normal stroma subtype (median 24 months, 1-year 82%, Fig. 2b). Both were notably absent in PDAC cell lines (Fig. 2c), which exhibited a distinct mitotic expression signature associated with mitotic checkpoints and DNA replication (Supplementary table 1)¹⁷. The fact that cell lines do not express these stromal factors, and that many metastatic samples express them at low levels (Supplementary Fig. 4), suggest that these genes are not expressed by the tumor epithelium. To further validate the stromal origin of these gene expression signatures, we isolated 6 CAF lines from primary tumors (Supplementary Fig. 5), and found that they robustly overexpressed our hypothesized stromal signatures compared to PDAC tumor cell lines which had no expression of the stromal signatures (Fig. 2c).

The vast majority of collagen gene expression was attributable to stromal compartments, with the lone exception being COL17A1, which was high in tumors. “Normal” stroma was characterized by relatively high expression of known markers for pancreatic stellate cells, smooth muscle actin, vimentin, and desmin, (ACTA2, VIM, and DES). Stellate cells have been shown to promote cancer cell survival in vitro¹⁸, but at the same time may restrain PDAC in mouse models^19,20. Targeting desmoplastic stroma by hedgehog pathway inhibition has been shown to both accelerate the development of disease²¹ and enhance delivery of chemotherapy²². In patients, the ratio of smooth muscle actin stained area to the collagen-stained area has been shown to be predictive of poor outcomes²³. We found that “activated” stroma was characterized by a more diverse set of genes associated with macrophages, such as the integrin ITGAM, and the chemokine ligands CCL13 and CCL18. “Activated” stroma also expressed other genes which point to its role in tumor promotion, including the secreted protein SPARC, WNT family members WNT2, and WNT5A, gelatinase B (MMP9), and stromelysin 3 (MMP11). The presence of fibroblast activation protein (FAP), which has previously been related to worse prognosis, in the activated stroma suggests that an activated fibroblast state may be partially responsible for the poor outcomes for these patients²⁴. This lead us to hypothesize that the “normal” stroma factor may describe a “good” version of stroma and that “activated” stroma factor may describe the activated inflammatory stromal response that has been seen in previous studies to be responsible for disease progression^25-27. Our factor analysis supports a complex, multi-gene model of stroma in PDAC, which may explain why single gene analysis has yielded mixed results.

Identification of tumor-specific subtypes

Independent of normal and stromal factors, we found that two tumor-specific factors define “classical” and “basal-like” subtypes of PDAC. When our samples were split into the two tumor subtypes (Fig. 3a), patients with basal-like subtype tumors had an overall worse median survival of 11 months and 44% 1-year survival compared to 19 months and 70% 1-year survival for those with classical subtype tumors (p=0.007, Fig. 3b). All cell lines assayed in this study (p < 0.001), as well as a majority of metastatic samples (p = 0.002), were classified as “basal-like”, suggesting that cell line models represent only one subset of PDAC. These subtypes as well as their prognostic value were independently validated within the recently published International Cancer Genome Consortium (ICGC) PDAC microarray data set (Fig. 3c, d)²⁸. Genes from the “basal-like” factor, including laminins and keratins, were also consistent with basal subtypes previously defined in bladder^14,28,29 and breast³⁰ cancers (Fig. 3e-h). Interestingly, genes from our “basal-like” subtype reproduced subtype calls (p<0.001) in breast cancer, had prognostic value in breast cancer samples (p<0.001) and reproduced previous subtype calls in bladder cancer (p<0.001). Given these promising results, we developed a single-sample cross-platform classifier of basal-like subtype which was trained on our microarray, TCGA bladder, and Perou breast cancer data, with a 93% cross validation accuracy, and was able to classify TCGA breast cancer data with 92% accuracy during external validation (Supplementary Fig. 6)

Potential subtypes of PDAC have previously been described by Collisson et al.⁴. We used the published exemplar genes for “exocrine-like”, “classical”, and “quasimesenchymal” subtypes to cluster normal pancreas, cell lines, and primary PDAC tumors from our cohort (Supplementary Fig. 7a). The three previous classifications were also observed in our data, but did not hold prognostic power either by cluster label or by supervised classification with PAM³¹(Supplementary Fig. 7b). Furthermore, inclusion of the Collisson et al. subtypes into a multivariate Cox regression with our proposed tumor subtypes did not remove the predictive power of our subtyping (p = 0.014). By cross-referencing the genes from Collisson et al.'s model with our NMF model, we observed three key findings. First, “exocrine-like” genes overlapped with genes from our exocrine pancreas factor (17/17). Tumors in this cluster had expression indistinguishable from adjacent normal samples from our data set. Second, Collisson et al.'s “classical” genes overlapped with our “classical” subtype genes (20/22), for which we retain the naming convention “classical”. Third, the gene set associated with “quasimesenchymal” subtype appeared to be a mixed collection of genes from our “basal-like” tumor (6/20) and stromal subtypes (6/20). Thus, the appearance of stromal factors in the Collisson et al. list of “quasimesenchymal” class genes may explain the apparent mesenchymal-like gene expression that was observed.

“Basal-like” and “classical” tumors were found within both “normal” and “activated” stroma subtypes (Fig. 4a). Differential prognosis among tumor and stroma subtypes was cumulative, as “classical” subtype tumors with “normal” stroma subtypes (n = 24) had the lowest hazard ratio of 0.39, with a 95% CI of [0.21, 0.73], while the “basal-like” subtype tumors with “activated” stroma subtypes (n = 26) had the highest hazard ratio of 2.28 with a 95% CI of [1.34, 3.87] (Fig. 4b). In a multivariate Cox regression model which included tumor subtypes, stromal subtypes, and clinical variables (gender, race, T stage, N stage, margin status, adjuvant therapy, histological grade, and age), both classifications were independently associated with survival (stroma subtypes: p = 0.037, tumor subtypes: p = 0.003).

Although basal-like subtype tumors have a worse prognosis, patients with basal-like subtype tumors showed a strong trend towards better response to adjuvant therapy (p=0.072, Fig. 4c,d). Among basal-like subtype patients, adjuvant therapy provided a hazard ratio of 0.38, (95% CI of [0.14, 1.09]), while in patients with classical subtype tumors, adjuvant therapy is associated with a hazard ratio of only 0.76 (95% CI [0.40, 1.43]). In our cohort, there was no association of most clinical variables (race, gender, T stage, N stage, differentiation, tumor cellularity) with survival, although positive nodal status trended towards significance, and positive margin status was significantly associated with worse survival (Table 1). Table 2 shows two-way associations of all subtype calls with clinical and pathological information from our cohort of PDAC patients. We found no association of tumor or stroma subtype with standard clinical or pathological variables, with the notable exception of mucinous features.

Tumor-specific subtypes in patient-derived xenografts

To assess the tumor or stromal specificity of our signatures, we performed RNAseq on a group of 37 PDX tumors. PDX tumors are composed of human tumor cells surrounded by mouse stroma (Supplementary Fig. 8)²⁹. Genes from both of our tumor signatures were expressed as human transcripts, whereas genes from both of our stromal signatures were expressed as mouse transcripts (Fig. 2d, Supplementary Fig. 9a). Furthermore, we found that PDX RNAseq expression divided PDX into both classical and basal-like groupings (Supplementary Fig. 9b) while predominantly expressing an activated stromal signature (Fig. 2d). We found that while tumor-specific subtype was not predictive of graft success (Fig. 5a), patient tumors with an activated stroma subtype had significantly higher graft success rates than those with normal stroma subtype or low amounts of stroma (Fig. 5b) (p=0.019). Basal-like subtype tumors also exhibited faster growth rates than classical tumors (p=0.032) as measured by the length of time that tumors took to grow to 200 mm³ (TT200, Fig. 5c-d), a previously used metric for PDX growth³⁰. Retrospective analysis of patients who had matched PDX tumors found that a shorter TT200 was associated with an unfavorable recurrence-free survival (p=0.035, Fig. 5e), suggesting that PDX tumor growth rate may reflect patient biology.

We also measured both mouse and human-specific expression of the Collisson et al. genes in our PDX models. We found that while genes from the “classical” subtype were expressed by human cells in PDX, “quasimesenchymal” transcripts were expressed by a mixture of human and mouse cells, and “exocrine-like” transcripts were infrequently expressed (Supplementary Fig. 7c). This supports our hypothesis that while the “classical” subtype is a bona fide group, the “quasimesenchymal” subtype is partially driven by non-tumor contributions of stroma and the “exocrine-like” subtype by normal pancreas.

KRAS codon mutations, tumor-specific subtypes, and race

Studies of KRAS codon mutations have demonstrated that different codon mutations may have differential functions^31,32 and in some clinical studies, have been shown to be associated with differential outcome. Because PDX tumors are enriched for human-specific tumor cells, we evaluated KRAS codon mutations in our PDX cohort using manually curated RNAseq data. While the overall frequency of KRAS codon mutations was similar to a recent study of PDAC⁷, we noted that the KRAS G12D mutation was significantly overrepresented in our basal-like subtype while G12V was isolated to the classical subtype (Figure 5f, p=0.030). Furthermore, we found an overrepresentation of KRAS G12V mutations in African Americans (Figure 5g, p<0.001). In contrast to basal-like breast cancers which occur most frequently in African American women and have a worse prognosis³³, African American patients in our cohort tended to have mainly classical subtype tumors (13 vs 2). We found that similar to other cancers, African Americans had a worse prognosis after adjusting for tumor subtype (Figure 4e, p=0.017). African American patients with classical subtype tumors had a median survival of 13 months compared to Caucasian patients with classical subtype tumors who had a median survival of 19 months.

Other commonly mutated genes and altered pathways in PDAC

Previously, loss of SMAD4 has been shown to promote tumor growth^34,35. We found that, similar to previous PDX studies of PDAC, loss of SMAD4 was associated with graft success in PDX models³⁶ (Fig. 5h, Supplementary Fig. 10, p=0.044). Furthermore, in our PDX cohort, we found that SMAD4 expression was significantly higher in classical compared to basal-like subtype PDX tumors (Fig. 5i, p=0.015), consistent with the observation that SMAD4 loss confers a more aggressive phenotype.

Using mutation, genomic subtype³, and gene expression²⁸ data from publically available ICGC data in which we show recapitulation of our subtypes and prognosis, we evaluated significantly mutated genes and pathways in PDAC, including ones recently identified through whole-exome sequencing of microdissected primary PDAC tumors^1-3,7. We found no significant associations between our expression subtypes and these mutationally altered pathways, i.e. TGFβ, RB, NOTCH, CTNNB1, SWI/SNF, and DNA repair (Supplementary Fig. 11). Furthermore, we found no overlap between our subtypes and recently identified genomic subtypes, or response to platinum therapy³. Consistent with this, a recent comprehensive study of somatic mutations in PDAC long-term survivors suggested that somatic mutations alone will not be sufficient to explain clinical outcome³⁷.

Given the overlap of our classical subtype with that of Collisson et al.⁵, it was not surprising to find that our classical subtype was also enriched for genes associated with GATA6 overexpression³⁸ (Supplementary Fig. 12a, Fig 4a). GATA6 has been found to promote epithelial cell differentiation^38,39. This prompted us to perform a more detailed histological markers of differentiation in our samples and found that samples with greater than 10% extracellular mucin, a marker of differentiation, comprised mostly of classical subtype tumors (88.5%, n=23) compared to only 11.5% (n=3) of basal-like subtype tumors (Supplementary Fig. 13, p=0.042, Table 2). Consistent with the increased presence of extracellular mucin, our classical subtype was enriched for genes upregulated in mucinous ovarian cancer (Supplementary Fig. 12b, WAMAUNYOKOLI_OVARIAN CANCER_GRADES_1_2_UP⁴⁰). Interestingly, our basal-like subtype was enriched for genes related to KRAS activation and STK11 loss in a lung cancer mouse model where STK11-deficient tumors demonstrated shorter latency and more frequent metastasis⁴¹ (Supplementary Fig. 12c). We found one sample with STK11 inactivation in the ICGC data; this sample was a basal-like subtype (Supplementary Fig. 11). Notably, our subtypes were not associated with other known signaling pathways in PDAC, including Fanconi anemia, DNA repair, chromatin remodeling, beta-catenin, RB, ARF, G1 (Fig. 4a). However, all of these pathways except for beta-catenin were considerably differentially expressed in cell lines compared to patient tumors, suggesting that gene expression in cell lines may be a deceptive representation of most tumors.

Low intra-patient heterogeneity in tumor-specific genes

We expect that only a subset of genes are relevant to the question of intra- and inter-patient heterogeneity in PDAC. Many methods exist to pre-select genes for supervised analysis³², but selection of the most differentially expressed genes is a common preprocessing step during unsupervised analysis³³. When clustering matched samples of metastatic and primary lesions using the 50 most differentially expressed genes among all matched samples, samples separated primarily by organ site instead of by patient (Fig. 6a, c). In contrast, when considering 25 top ranked exemplar genes each from the “basal-like” and “classical” factors, samples from the same patient clustered closer together, and were less dependent of organ site (Fig. 6b, d). This is further illustrated in a focused analysis of two patients (Fig. 6), whose tumor samples appear patient-specific when considering our tumor subtype gene list, but cluster by site when considering differentially expressed genes. Overall, we found that our tumor subtype gene list showed higher similarity (mean Pearson's ρ=0.53) between all other samples from the same patient than did the differentially expressed gene list (ρ=0.32, t-test p<0.001). Furthermore, our tumor subtype gene list produced much lower similarity among all other samples from the same organ site across different patients (ρ = 0.04) than the differentially expressed gene list (ρ=0.34, p<0.001). This observed similarity of tumor gene expression among tumors within the same patient suggests overall high inter-patient tumor heterogeneity and low heterogeneity between primary and metastatic sites. However, we did observe examples of intra-patient heterogeneity between metastatic sites. For example, lung metastases, even those from patients with “basal-like” tumors in other locations, clustered exclusively with the “classical” tumors, suggesting that some intra-patient heterogeneity may exist among metastatic sites, and supporting the previously reported divergent patterns of failure in PDAC³⁴.