This article has Open Peer Review reports available.
Prognostic relevance of molecular subtypes and master regulators in pancreatic ductal adenocarcinoma
- Rekin’s Janky†1Email author,
- Maria Mercedes Binda†2,
- Joke Allemeersch3,
- Anke Van den broeck2,
- Olivier Govaere4,
- Johannes V. Swinnen5,
- Tania Roskams4,
- Stein Aerts†1Email author and
- Baki Topal†2Email author
© The Author(s). 2016
Received: 21 February 2016
Accepted: 8 July 2016
Published: 12 August 2016
Pancreatic cancer is poorly characterized at genetic and non-genetic levels. The current study evaluates in a large cohort of patients the prognostic relevance of molecular subtypes and key transcription factors in pancreatic ductal adenocarcinoma (PDAC).
We performed gene expression analysis of whole-tumor tissue obtained from 118 surgically resected PDAC and 13 histologically normal pancreatic tissue samples. Cox regression models were used to study the effect on survival of molecular subtypes and 16 clinicopathological prognostic factors. In order to better understand the biology of PDAC we used iRegulon to identify transcription factors (TFs) as master regulators of PDAC and its subtypes.
We confirmed the PDAssign gene signature as classifier of PDAC in molecular subtypes with prognostic relevance. We found molecular subtypes, but not clinicopathological factors, as independent predictors of survival. Regulatory network analysis predicted that HNF1A/B are among thousand TFs the top enriched master regulators of the genes expressed in the normal pancreatic tissue compared to the PDAC regulatory network. On immunohistochemistry staining of PDAC samples, we observed low expression of HNF1B in well differentiated towards no expression in poorly differentiated PDAC samples. We predicted IRF/STAT, AP-1, and ETS-family members as key transcription factors in gene signatures downstream of mutated KRAS.
PDAC can be classified in molecular subtypes that independently predict survival. HNF1A/B seem to be good candidates as master regulators of pancreatic differentiation, which at the protein level loses its expression in malignant ductal cells of the pancreas, suggesting its putative role as tumor suppressor in pancreatic cancer.
The study was registered at ClinicalTrials.gov under the number NCT01116791 (May 3, 2010).
Pancreatic ductal adenocarcinoma (PDAC; also called pancreatic cancer) is one of the most aggressive cancers, associated with a poor prognosis . The lack of early diagnostic markers and efficient therapeutic modalities for PDAC results in extremely poor prognosis. For several decades, many efforts have been undertaken to better understand the pathogenesis and biology of PDAC, and to improve patient survival through early diagnosis and various therapeutic strategies. However, no substantial advances have been made to overcome its lethal destiny. Today, adequate surgical resection is the only chance for patients to be cured from PDAC, often in combination with peri- or post-operative chemo(radio)therapy [2, 3]. Unfortunately, only selected patients with localized disease are potential candidates for surgical management with curative intent. Even in the group of surgically treated curable patients, the majority will develop cancer recurrence and die within two years. Most patients with pancreatic cancer are not eligible for surgery as they present in advanced stages with distant organ metastases and/or locoregional extension. Systemic chemotherapy is the standard of care for patients with advanced inoperable PDAC, resulting in a median survival of about 8 months .
As currently available clinicopathological classification systems and treatment modalities fail to tailor patient management or improve survival substantially, molecular subtyping of PDAC may help unravel its mechanisms of carcinogenesis and progression, and help discover efficient therapeutic molecules. The quest to identify clinically relevant gene signatures of PDAC has been a rough journey resulting in a wide range of often non-reproducible or conflicting data. Recently, based on 27 microdissected surgical samples, three subtypes of PDAC (classical, quasimesenchymal, and exocrine-like) were identified and their gene signatures defined as PDAssign. Despite its small sample size the study presented a prognostic relevance for these subtypes . The aim of our study was to evaluate the prognostic relevance of molecular subtypes and identify key transcription factors as master regulators in a large cohort of PDAC patients. Hereto, in contrast to other studies, we analyzed also several relevant clinicopathological variables that have proven to influence survival significantly.
Between 1998 and 2010, tissue samples were collected, after written informed consent, from patients who underwent pancreatic resection for PDAC. Snap-frozen tissue samples were stored in liquid nitrogen and/or at −80 °C in RNALater (Qiagen) until further use. From the primary tumor of 171 patients and from surrounding non-tumoral pancreatic (control) tissue of 14 patients, total RNA was extracted using the RNeasy Mini kit (Qiagen) according the manufacturer’s instructions. Only samples with an RNA integrity number (RIN) of >7.0 were used for further analysis, i.e. 118 PDAC samples (male/female ratio: 65/53; age: 32–87 years with median of 64 years) and 13 control tissues (male/female ratio: 8/5; age: 51–78 years with median of 67 years). Two pathologists confirmed PDAC samples to contain at least 30 % cancer cells. Patients with pre-operative radio- or chemotherapy were excluded from the study.
RNA concentration and purity were determined spectrophotometrically using the Nanodrop ND-1000 (Nanodrop Technologies) and RNA integrity was assessed using a Bioanalyser 2100 (Agilent). Per sample, an amount of 100 ng of total RNA spiked with bacterial RNA transcript positive controls (Affymetrix) was amplified and labeled using the GeneChip 3′ IVT express kit (Affymetrix). All steps were carried out according to the manufacturers protocol (Affymetrix). A mixture of purified and fragmented biotinylated amplified RNA (aRNA) and hybridisation controls (Affymetrix) was hybridized on Affymetrix Human Genome U219 Array Plate followed by staining and washing in the GeneTitan® Instrument (Affymetrix) according to the manufacturer’s procedures. To assess the raw probe signal intensities, chips were scanned using the GeneTitan® HT Array Plate Scanner (Affymetrix).
Microarray data analysis
Analysis of the microarray data was performed with the Bioconductor/R packages  (http://www.bioconductor.org). The analysis was based on the Robust Multi-array Average (RMA) expression levels of the probe sets, computed with the package xps. Differential expression was assessed via the moderated t-statistic implemented in the limma package, described in . To control the false discovery rate, multiple testing correction was performed  and probe sets with a corrected p-value below 0.05 and an absolute fold change larger than two were selected.
Molecular subtype discovery
Intrinsically variable genes were first selected based on their expression variation over the 118 PDAC samples (2374 genes with s.d. > 0.8). The “PDAssign” genes were selected as the variable genes matching the published signature , i.e. 62 genes excluding 3 genes without probes in our microarray platform (CELA3B, PRSS2, SLC2A3) and 3 genes that are not variable (SLC16A1, GPM6B, SLC5A3).
Identification of subclasses using non-negative matrix factorization clustering
Subclasses of a data set consisting of unified expression data of 118 samples and variable genes were computed by reducing the dimensionality of the expression data from thousands of genes to a few metagenes by applying a consensus non-negative matrix factorization (NMF) clustering method (v5) [9, 10]. This method computes multiple k-factor factorization decompositions of the expression matrix and evaluates the stability of the solutions using a cophenetic coefficient. Consensus matrices and sample correlation matrices were calculated for 2 to 5 potential subtypes (k) using default parameters and Euclidian distance. The final subclasses were defined based on the most stable k-factor decomposition and visual inspection of sample-by-sample correlation matrices. For this we used the NMF clustering implemented from Gene Pattern software package .
Merging microarray data using DWD
Distance Weighted Discrimination (DWD) method  was applied for batch correction to the data of Collisson et al. and our expression data on variable genes after row median centering and column normalization according to the authors’ protocol . The Java version of DWD was used with default parameters (Standardized DWD, centered at zero).
Gene Set Enrichment Analysis (GSEA) was used to score how enriched the modules and regulons (identified above in the first section) were in the top differentially expressed genes for a given contrast . We performed the GSEA Preranked analysis using the list of the genes ranked by the signed p-value from each of the supervised and unsupervised biological contrasts (e.g. PDAC vs Control, k2.cl1 vs k2.cl2). This algorithm scores the positive or negative enrichment for all modules/regulons at the top or the bottom of the ranking. We also used WebGestalt , in which the hyper-geometric test was used for enrichment analysis and the Benjamini-Hochberg procedure was used to control the False Discovery Rate.
Top 250 KRAS dependency signature probes were extracted from Singh et al.  and provided a list of 187 genes, of which 165 genes were in our microarray data and 77 genes showed variable expression (sd > 0.8). The list of 77 genes was ranked according to their KRAS dependency and was used to make an expression heatmap of the 118 PDAC samples. Expression heatmaps are generated using R package heatmap. Hierarchical clustering based on a Spearman rank correlation as distance metric and an average linkage method (R function hclust) was used predicting 112 samples (95 %) as KRAS dependent samples (high level KRas activity) and 6 samples as KRAS independent (low level KRAS activity). The R function cutree automatically cut each dendrogram (from the top down) to form two groups of samples. KRAS expression levels are also significantly higher in KRAS dependent samples compared to other samples (p = 0.002).
Kaplan-Meier estimates were used for survival analysis. Overall survival (OS) was defined as time from surgery to death, irrespective of cause. Disease-free survival (DFS) was defined as time to tumor recurrence or death, irrespective of cause. Patients were followed up until death or until the date of study closure on November 2014. Together with the molecular subclasses the effect on survival of a set of 16 clinico-pathological prognostic factors was evaluated: patient age (years), gender (male/female), PDAC location (head/body or tail), tumor diameter (mm), differentiation grade (pG), depth of tumor invasion (pT), locoregional lymph node metastasis (pN), distant organ metastasis (pM), completeness of tumor resection (pR), magnitude of the surgical resection margin (pRM), perineural invasion (PNI), vascular invasion (VI), lymph vessel invasion (LVI), extra-capsular lymph node invasion (ECLNI), AJCC TNM Classification 7th Edition, adjuvant systemic chemotherapy (Yes/No). Log-rank tests and Cox regression models were used to verify the relation between a set of predictors and survival. A multivariable model was constructed combining the predictors with p < 0.10 in the univariable models, and p values less than 0.05 were considered significant.
Master regulator analysis
In order to characterize regulatory networks underlying the subtypes, we used iRegulon  to identify master regulators, i.e. transcription factors whose regulons (transcriptional target sets) are highly overlapping with the observed gene signatures. The master regulators are expected to be directly activated by signal transduction. In this approach, we use a large collection of transcription factor (TF) motifs (9713 motifs for 1191 TFs) and a large collection of ChIP-seq tracks (1120 tracks for 246 TFs).
Briefly, this method relies on a ranking-and-recovery strategy where the offline ranking aims at ranking 22284 genes of the human genome (hg19) scored by a motif discovery step integrating multiple cues, including the clustering of binding sites within cis-regulatory modules (CRMs), the potential conservation of CRMs across 10 vertebrate genomes, and the potential distal location of CRMs upstream or downstream of the transcription start site (TSS+/−10 kb). The recovery step calculates the TF enrichment for each set of genes, i.e. genes from co-expression modules, leading to the prediction of the TFs and their putative direct target genes in the module. An important advance of this method is that it can optimize the association of TFs to motifs using not only direct annotations, but also predictions of TF orthologs and motif similarity, allowing the discovery of more than 1191 TFs in human.
Samples (n = 6) showing top differential expression for HNF1B were selected for HNF1B immunohistochemistry staining (IHC). Five-micrometer-thick sections were prepared from formalin-fixed paraffin-embedded PDAC specimens. Stainings were made using the Benchmark Ultra (Ventana). Briefly, samples were deparaffinized at 72 °C and endogenous peroxidase activity was blocked using 0.3 % H2O2. Antigens were retrieved by heating the sections for 68 min at 91 °C in citrate buffer, pH6. Sections were incubated with the primary antibody against human HNF1B (Sigma, catalogue number HPA002083) dissolved 1:200 in Dako REAL antibody diluent at 37 °C for 32 min. The reaction product was developed using ultraView Universal DAB Detection Kit and sections were counterstained with hematoxylin. Sections were washed, dehydrated in progressively increasing concentration of ethanol and xylene, and mounted with xylene-based mounting medium. Normal human pancreas was used as a positive control. In order to check unspecific antibody binding, negative controls, in which the primary antibody was omitted, were also done. Samples were carefully analyzed by a pathologist. Slides were visualized using Leica DMR microscope (Leica Microsystems Ltd, Germany) and photographs were taken using Leica Application Suite v3.5,0 software (Leica Microsystems, Switzerland). HNF1B staining was scored based on intensity (on a scale from 0–3; 0, negative; 1, weak; 2, positive; 3, strong) and the proportion of reactive cells (0–100 %); histoscore was determined by multiplying both parameters (range 0–300) as published in Hoskins et al. . When more than one magnification area was available from a given tumor, the mean score was used.
Gene expression profiling
We applied gene expression profiling using microarrays on 118 tumor and 13 histologically normal pancreatic tissue samples (control) to investigate the molecular mechanisms driving PDAC and its different subtypes. Gene expression analysis of PDAC samples was performed on whole tumor tissue, i.e. cancer cells (at least 30 % of sample) and tumor stroma. Differential gene expression analysis using the contrast of all PDAC samples versus all control samples provided a large number (n = 6873) of genes that were differentially expressed (corrected p-value < 0.05; Additional file 1: Table S1). Our findings are in agreement with previously published pancreatic cancer gene expression data . When we compared the gene expression profile of each tumor sample against a published KRas dependent gene signature , we found 94 % of our samples (112/118) to be KRas-dependent, which is in agreement with the fact that more than 90 % of PDAC have a KRAS driver mutation (Additional file 2: Figure S1) [19, 20].
Molecular subtypes linked to survival
Results of univariable and multivariable Cox regression models for disease-free survival (DFS)
Number of Patients
Disease-free Survival Time (DFS; median (CI): months)
Hazard ratio (HR) (95% CI)
Hazard ratio (HR) (95% CI)
< 64 y.
10.2 (6.4 - 13.3)
0.756 (0.508 - 1.122)
> 64 y.
9.0 (7.3 - 10.9)
10.9 (8.1 - 13.5)
0.785 (0.528 - 1.161)
7.8 (6.0 - 11.1)
9.8 (7.4 - 12.4)
0.564 (0.357 - 0.921)
0.581 (0.340 - 1.015)
Body or Tail
8.3 (4.4 - 11.1)
< 2 cm
13.1 (6.7 - 17.0)
0.676 (0.411 - 1.064)
0.739 (0.399 - 1.306)
> 2 cm
7.8 (3.3 - 18.9)
13.1 (4.0 - 24.3)
10.2 (7.1 - 16.2)
7.7 (6.0 - 11.0)
7.7 (7.3 - 8.1)
11.6 (7.1 - 17.1)
9.0 (6.8 - 11.4)
8.4 (7.4 - 11.2)
10.3 (6.1 - 13.1)
0.839 (0.557 - 1.248)
8.8 (7.3 - 11.4)
10.0 (7.7 - 12.3)
0.411 (0.230 - 0.802)
0.441 (0.097 - 1.458)
4.7 (3.2 - 10.9)
10.0 (7.8 - 12.4)
0.712 (0.458 - 1.149)
7.7 (7.4 - 11.2)
< 1 mm
10.0 (6.3 - 13.1)
0.925 (0.621 - 1.382)
> 1 mm
9.8 (7.0 - 12.3)
17.1 (3.3 - 56.1)
0.528 (0.264 - 0.955)
0.516 (0.233 - 1.053)
9.8 (7.4 - 11.4)
11.1 (7.0 - 18.9)
0.629 (0.399 - 0.969)
0.805 (0.485 - 1.303)
8.0 (6.4 - 11.0)
10.2 (6.1 - 14.3)
0.833 (0.529 - 1.276)
8.8 (7.3 - 11.4)
10.4 (7.3 - 13.3)
0.821 (0.544 - 1.254)
8.4 (6.3 - 11.1)
AJCC TNM Stage 7th Ed.
10.9 (6.8 - 14.3)
0.730 (0.474 - 1.100)
8.5 (7.0 - 11.1)
N0 <T3 M0 (Early)
10.9 (6.8 - 14.3)
Early vs Adv 0.498 (0.273 - 0.950)
N1 <T3 M0 (LNM)
9.8 (7.3 - 12.6)
Early vs Adv 0.035
0.915 (0.295 - 4.022)
T4 or M1 (Advanced)
5 (3.3 - 10.4)
7.4 (4.6 - 11.0)
1.139 (0.733 - 1.728)
10.0 (7.7 - 12.6)
11.6 (7.4 - 16.2)
0.655 (0.440 - 0.976)
0.252 (0.092 - 0.888)
7.8 (6.7 - 10.0)
13.5 (10.9 - 17.1)
Cl1 vs Cl2 0.602 (0.382 - 0.940)
8.0 (7.0 - 10.0)
Cl1 vs Cl3 0.615 (0.363 - 1.066)
Cl1 vs Cl2 0.026
4.7 (3.3 - 11.2)
Cl1 vs Cl3 0.082
13.3 (9.8 - 16.7)
Cl1 vs Cl2 0.670 (0.418 - 1.065)
8.4 (6.7 - 10.0)
Cl1 vs Cl2 0.090
11.0 (4.8 - 17.1)
4.7 (3.3 - 11.2)
13.5 (10.9 - 17.0)
Cl1 vs Cl5 0.488 (0.251 - 1.021)
8.4 (7.0 - 10.2)
Cl1 vs Cl5 0.057
9.3 (4.8 - NA)
4.7 (3.3 - 11.2)
7.5 (4.3 - 11.0)
Results of univariable and multivariable Cox regression models for overall survival (OS)
Number of patients
Overall survival time (OS; median (CI): months)
Hazard ratio (HR) (95% CI)
Hazard ratio (HR) (95% CI)
< 64 y.
23.5 (12.6 - 33.0)
0.626 (0.422 - 0.924)
0.551 (0.350 - 0.859)
> 64 y.
13.7 (11.4 - 16.8)
20.5 (12.9 - 29.3)
0.865 (0.587 - 1.269)
12.6 (11.2 - 20.1)
19.5 (12.6 - 23.5)
0.600 (0.384 - 0.967)
0.714 (0.435 - 1.209)
Body or Tail
12.6 (10.0 - 26.4)
< 2 cm
20.5 (11.7 - 36.5)
0.746 (0.459 - 1.165)
> 2 cm
14.8 (11.9 - 20.8)
22.5 (1.5 - 33.2)
14.8 (11.2 - 29.3)
15.9 (11.5 - 23.5)
13.4 (12.3 - 14.6)
26.9 (9.4 - 38.8)
15.9 (11.8 - 21.0)
12.4 (1.3 - 33.0)
21.0 (14.6 - 29.3)
0.750 (0.499 - 1.111)
12.8 (11.7 - 17.8)
17.8 (12.9 - 23.5)
0.569 (0.323 - 1.097)
0.624 (0.166 - 1.928)
11.4 (5.8 - 12.4)
16.8 (12.9 - 25.6)
0.733 (0.474 - 1.174)
12.4 (7.0 - 23.4)
< 1 mm
15.4 (11.8 - 25.6)
1.096 (0.737 - 1.622)
> 1 mm
16.7 (12.3 - 23.5)
37.8 (10.6 - NA)
0.468 (0.227 - 0.860)
0.561 (0.252 - 1.115)
15.9 (12.4 - 20.8)
19.7 (11.9 - 33.2)
0.730 (0.466 - 1.115)
12.8 (11.5 - 23.4)
19.7 (10.6 - 33.0)
0.877 (0.561 - 1.334)
15 (12.3 - 21.7)
20.1 (12.9 - 30.5)
0.654 (0.437 - 0.987)
0.660 (0.398 - 1.089)
12.4 (10.2 - 20.8)
AJCC TNM Stage 7th Ed.
23.5 (16.8 - 31.7)
0.681 (0.443 - 1.023)
0.672 (0.166 - 2.254)
12.6 (11.4 - 16.7)
Early (pN=0,pT≤ 3,pM=0)
23.5 (16.8 - 31.7)
Early vs LNM 0.722 (0.463 - 1.106)
14.8 (11.2 - 21.7)
LNM vs Adv 0.736 (0.425 - 1.328)
Early vs LNM 0.136
Advanced (pT=4 or pM=1)
11.7 (6.6 - 12.4)
Early vs Adv 0.532 (0.295 - 0.997)
Early vs Adv 0.049
12.1 (7.0 - 16.8)
1.337 (0.874 - 2.002)
19.8 (13.9 - 25.6)
20.9 (12.9 - 29.3)
0.710 (0.482 - 1.048)
0.247 (0.092 - 0.860)
12.7 (11.2 - 16.7)
24.6 (16.8 - 33.2)
Cl1 vs Cl2 0.680 (0.437 - 1.050)
12.7 (11.2 - 16.7)
Cl2 vs Cl3 1.055 (0.641 - 1.788)
Cl1 vs Cl2 0.082
11.8 (6.6 - 20.5)
Cl1 vs Cl3 0.717 (0.426 - 1.236)
Cl1 vs Cl3 0.226
Cl1 vs Cl3 0.209 (0.057 - 0.809)
23.5 (14.8 - 33.0)
12.6 (10.6 - 16.7)
26.2 (9.3 - 38.8)
11.9 (6.6 - 21.0)
23.5 (14.8 - 33.2)
Cl1 vs Cl5 0.398 (0.210 - 0.808)
13.9 (11.4 - 21.7)
Cl2 vs Cl5 0.483 (0.251 - 0.988)
Cl1 vs Cl5 0.012
28.9 (9.3 - NA)
Cl3 vs Cl5 0.298 (0.067 - 0.949)
Cl2 vs Cl5 0.046
11.8 (6.6 - 20.5)
Cl3 vs Cl5 0.040
10 (6.9 - 16.6)
Functional analysis of molecular subtypes
PDAC subtypes are poorly characterized at the molecular level and little is known about the regulatory networks underlying the expression of the genes driving better or worse survival. As we could reproduce the three subtypes (NMF with k = 3, or briefly “k3”) and confirmed their prognostic relevance, we aimed to further characterize their gene expression profiles, functions, and pathways. Compared to normal tissue samples, all subtypes are enriched for “Neoplasms”, “invasiveness”, and “integrin family cell surface interactions”, and all subtypes are comparably enriched for typical pancreatic cancer gene signatures (FDR = 0.000, NES > =2.41).
When the k3 subtypes are compared directly against each other (Additional file 1: Table S1), we could define cluster-specific gene signatures as the genes that are specifically over- or under-expressed for a given subtype and missing PDAssign genes were added to these signatures to perform functional enrichment analysis (Additional file 4: Figure S3). For example, we found a specific gene signature with 148 genes over-expressed and 3 under-expressed in the predicted exocrine-like subtype that is enriched for processes related to the exocrine pancreas, such as pancreatic secretion and protease activity. For the QM subtype we identified 50 up-regulated genes specific for this subtype with 132 down-regulated genes, and this set of genes shows typical properties of epithelial and mesenchymal cancers. Focusing further on Epithelial-to-Mesenchymal Transition (EMT) properties, we found an enrichment of an EMT signature (NES = 2.38). Some EMT TFs, such as TWIST1 and SNAI2, show QM subtype specific expression. However, although this signature resembles some aspects of EMT, it does not capture the entire EMT signature, since there is limited gene overlap with a core mesenchymal transition signature derived by meta-analysis across cancer types . Notice that samples clustered by low and high expression of mesenchymal cancer attractors do not show a significant link with survival. Finally, the predicted classical subtype has very few specific genes compared to the other subtypes (only 14 genes), and lacks any specific biological pathway enrichment. Overall, despite a partial gene overlap with the published PDAssign genes (36.4 %, 20/55) (Additional file 4: Figure S3e), our larger cluster-specific gene signatures agree with the known description of the PDAC subtypes.
Master regulators of PDAC
We also found a ZEB1 motif (NES = 3.91) in the regulatory analysis of 1325 down-regulated genes (Fig. 3c; clustered with “LMO2” motifs) while ZEB1 is up-regulated in PDAC samples (log ratio = 2.17, p-value = 2.32 × 10−21). This finding is consistent with ZEB1 being a repressor . Expression of ZEB1 has been shown recently to be a strong predictor of survival in PDAC  and is a known TF inducing epithelial-mesenchymal transition (EMT) in cancer cells. Finally, we identified enriched GATA3 ENCODE tracks in the “classical” and “QM-PDA” specific gene signatures (NES ~ 4), but not in the contrasts of PDAC vs control (data not shown). Thus, besides the role of GATA6 in QM-PDA, as proposed by Collison et al., our data also suggests that GATA3 may be functional in the two other subtypes.
Within the set of 1325 down-regulated genes in PDAC versus control, the most strongly enriched TF motifs were those for HNF1A/B (NES = 5.036, Fig. 3c). The HNF1A/B regulon, defined by 320 predicted target genes, is furthermore differentially expressed between classical and QM subtypes (Additional file 5: Figure S4). HNF1A/B is also found as top enriched regulator (NES = 8.156) when using a gene signature specific for the classical subtype compared to the exocrine subtype (data not shown). Compared to HNF1A, HNF1B is the best candidate to bind to this motif because the HNF1B gene itself is also down-regulated in the tumor samples (log ratio = −1.34, p-value = 6.24 × 10−5) and its expression profile is strongly correlated with the predicted targets (Pearson correlation = 0.71, p-value < 2.2 × 10−16), although HFN1A is also strongly correlated with these genes (Pearson correlation = 0.52, p-value = 1.67 × 10−10).
HNF1B protein expression in PDAC
In a recent attempt to unravel the tumor biology of pancreatic ductal adenocarcinoma (PDAC), Collisson et al. reported the PDAssign gene signature to classify this lethal cancer into three molecular subtypes with prognostic relevance . The association of PDAssign with survival was based on gene expression data for 27 patients. In the current study, we evaluated the validity of PDAssign in a large cohort of 118 pancreatic cancer patients treated with surgery with or without adjuvant systemic chemotherapy. Apart from the sample sizes, another major difference between these two studies is the fact that we used whole-tumor samples including the micro-environment, whereas the former study used microdissection to enrich their samples for cancer cells. While microdissection of cells in fixed tissue could possibly be associated with higher levels of RNA degradation , we used high-quality samples with a pathologically proven minimum of 30 % cancer cells. By doing so, we kept the molecular information of the microenvironment, we have reduced RNA contamination and the large number of samples improves the signal-to-noise ratio. A future perspective may be to decipher the tumour specific response using single cell technology. Interestingly, a recent study shows that we can defined stroma and tumour specific subtypes by applying a similar NMF approach on a compendium of microarray expression data including 145 primary and 61 metastatic PDAC tumour samples . They also identify two stromal subtypes, normal and activated, with the latter showing the worse prognosis. Nonetheless, we could confirm PDAssign to be a reliable classifier of PDAC into three distinct molecular subtypes with prognostic relevance. With their approach, Moffit et al.  virtually dissected the samples to identify a ‘classical’ and a ‘basal-like’ tumor-specific subtypes showing similarities to our predicted clusters with the ‘basal-like’ subtype showing genes of the same family of the ‘quasi-mesenchymal’ subtype with the worse survival. In our survival analysis, the clinicopathological factors are not independent predictors of survival, including stage and grade features. This is consistent with the recent studies [5, 31] but not with the previous literature , which may be due to the larger size of the recent studies. Moreover, we found these molecular subtypes as independent predictors of both disease-free and overall survival. We confirmed the classical PDAC subtype to be associated with the best survival, though in contrast to Collisson et al., we showed the exocrine-like subtype to be associated with poor prognosis and comparable survival to that of the quasi-mesenchymal subtype. We therefore envision that future evaluation of these molecular subtypes in larger studies may provide new insights in novel treatment strategies, opening new perspectives in personalized targeted therapy for PDAC.
In order to better understand the biology of PDAC we used iRegulon to identify transcription factors as master regulators of PDAC and its subtypes. We found that HNF1A/B are among thousand TFs the top enriched master regulators of the genes expressed in the normal pancreatic tissue compared to the PDAC regulatory network. On immunohistochemistry staining of PDAC samples we confirmed low expression of HNF1B in well differentiated tumors and no expression in six poorly differentiated PDAC samples. Our IHC results are also confirmed in an independent study from Jiang X et al. , i.e. positive staining for HNF1β in acinar parenchyma and ducts from normal pancreas and negative or moderate staining for PDAC samples. HNF1β also plays an important role in human normal pancreas morphogenesis and terminal differentiation of pancreatic β-cells . Moreover, HNF1β is involved in regulating the β-cell transcription factor network and is necessary for glucose sensing or glycolytic signalling in the pancreatic β-cells . HNF1B was also found recently to be down-regulated in vitro in PDAC cells by a microRNA mechanism involving hsa-miR-24 and/or hsa-miR-23a . In the suggested mechanism, HNF1B deregulation in PDAC results in loss of the expression of the adhesion molecule E-cadherin, which induces epithelial-mesenchymal transition (EMT) and allows cells to detach from cell agglomerations and to migrate. In another study , HNF1B was also found deregulated in a mouse model of intraductal papillary mucinous neoplasm (IPMN) to PDAC progression with another duct-specific factor, SOX9, while the latter was not found deregulated in our data and was not predicted in our regulatory analysis. The authors highlighted the importance of these factors for the loss of mature ductal identity in tumor initiation. However, HNF1A has also been revealed as a specific key regulator of the transcriptome in pancreatic tumor tissues and was suggested as an important tumor suppressor in the pancreas . Hoskins et al. observed that inducible over-expression of HNF1A in pancreatic tumor-derived cells could generate growth inhibition, a G0/G1 cell cycle arrest and apoptosis. Taken together, these observations suggest that HNF1A and HNF1B can be co-expressed in normal pancreatic tissues and may act as tumor suppressors through their regulatory activity. These factors can dimerize as homo- or hetero-dimers and can present several tissue-specific and species-specific isoforms , which can explain why we can find their activity independently. Low expression of HNF1B in some PDAC samples could reflect the tumor origination from acinar cells with incomplete ductal reprogramming phenotype (as suggested by one of our peer reviewer). Additionally, we identified IRF/STAT, AP-1, and ETS-family members as key transcription factors in gene signatures downstream of mutated KRAS. However, this approach only captures a part of the regulatory network while the post-transcriptional regulation and the microRNA regulatory network were not taken into account. We believe these key TFs or master regulators represent a valuable set of molecules for further study in functional assays and in vivo experiments to assess their role in PDAC carcinogenesis, progression, and novel therapeutic strategies.
This is the first study describing in a large cohort of pancreatic cancer patients the prognostic relevance of molecular subtypes, which are driven by the PDAssign gene signature. Our results show molecular subtypes, but not clinicopathological factors, as independent predictors of survival. We have identified enriched transcription factors (TFs) as putative master regulators of PDAC, and their downstream networks, using iRegulon. Among them, the hepatocyte nuclear factor 1 homeobox A/B (HNF1A or HNF1B) and its predicted targets are globally down-regulated in PDAC. Immunohistochemistry for HNF1B shows a strong nuclear staining of normal pancreatic ductal cells, whereas its expression is low in malignant ductal cells of well differentiated and absent in poorly differentiated PDAC samples. As these TFs play a key role in PDAC, they may involve novel therapeutic targets to improve the survival of patients with PDAC.
CRM, cis-regulatory module; DFS, disease-free survival; DWD, distance weighted discrimination; EMT, epithelial-to-mesenchymal transition; GSEA, gene set enrichment analysis; HNF1A/B, hepatocyte nuclear factor 1 homeobox A/B; IHC, immunohistochemistry; NES, normalized enrichment score; NMF, non-negative matrix factorization; OS, overall survival; PDAC, pancreatic ductal adenocarcinoma; QM, quasi-mesenchymal; TF, transcription factor
We thank Frank Vanderhoydonc (Laboratory of Lipid Metabolism and Cancer, KU Leuven) and Kathleen Van den Eynde (Translational Cell & Tissue Research, KU Leuven) for their technical assistance.
RJ is supported by postdoc fellowships from Belspo, KU Leuven Research Fund (F+) and FWO Belgium.
BT is supported by a basic-clinical research mandate (Fundamenteel Klinisch Mandaat) from the FWO Belgium.
This study is supported by i) an unrestricted grant from Johnson & Johnson Medical Devices, Belgium; ii) the Special Research Fund (BOF) KU Leuven (http://www.kuleuven.be/research/funding/bof/) (grant PF/10/016 to SA), and iii) the Foundation Against Cancer (http://www.cancer.be) (grant 2012-F2 to SA).
Availability of data and materials
The dataset supporting the conclusions of this article is available in the NCBI Gene Expression Omnibus (GEO) repository, [GSE62165, http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62165].
RJ, MB, JS, SA and BT participated in the study design. RJ, MB, SA and BT drafted the manuscript. RJ carried out the bioinformatics, gathered the raw data and performed the statistical analysis. MB, AVdB and OG carried out the laboratory experiments. JA carried out the gene expression analyses. TR carried out the pathological examination. SA conceived of the study and supervised the bioinformatics. BT conceived of the study and its coordination, and supervised the clinical aspects. All authors reviewed and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
The study was approved by the UZ/KU Leuven Ethical Committee prior to sample analysis and was given study number ML6615. The study was registered at clinicaltrials.gov under the number NCT01116791. Tissue samples were collected, after written informed consent, from patients who underwent pancreatic resection for PDAC.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Stewart BW, Wild CP, editors. World Cancer Report 2014. Lyon: IARC Nonserial Publication; 2014.Google Scholar
- Mukherjee S, Hurt CN, Bridgewater J, Falk S, Cummins S, Wasan H, Crosby T, Jephcott C, Roy R, Radhakrishna G, McDonald A, Ray R, Joseph G, Staffurth J, Abrams RA, Griffiths G, Maughan T. Gemcitabine-based or capecitabine-based chemoradiotherapy for locally advanced pancreatic cancer (SCALOP): a multicentre, randomised, phase 2 trial. Lancet Oncol. 2013;14:317–26.View ArticlePubMedPubMed CentralGoogle Scholar
- Oettle H, Neuhaus P, Hochhaus A, Hartmann JT, Gellert K, Ridwelski K, Niedergethmann M, Zülke C, Fahlke J, Arning MB, Sinn M, Hinke A, Riess H. Adjuvant chemotherapy with gemcitabine and long-term outcomes among patients with resected pancreatic cancer: the CONKO-001 randomized trial. JAMA. 2013;310:1473–81.View ArticlePubMedGoogle Scholar
- Von Hoff DD, Ervin T, Arena FP, Chiorean EG, Infante J, Moore M, Seay T, Tjulandin S a, Ma WW, Saleh MN, Harris M, Reni M, Dowden S, Laheru D, Bahary N, Ramanathan RK, Tabernero J, Hidalgo M, Goldstein D, Van Cutsem E, Wei X, Iglesias J, Renschler MF. Increased survival in pancreatic cancer with nab-paclitaxel plus gemcitabine. N Engl J Med. 2013;369:1691–703.View ArticleGoogle Scholar
- Collisson EA, Sadanandam A, Olson P, Gibb WJ, Truitt M, Gu S, Cooc J, Weinkle J, Kim GE, Jakkula L, Feiler HS, Ko AH, Olshen AB, Danenberg KL, Tempero M a, Spellman PT, Hanahan D, Gray JW. Subtypes of pancreatic ductal adenocarcinoma and their differing responses to therapy. Nat Med. 2011;17:500–3.View ArticlePubMedPubMed CentralGoogle Scholar
- Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, Hornik K, Hothorn T, Huber W, Iacus S, Irizarry R, Leisch F, Li C, Maechler M, Rossini AJ, Sawitzki G, Smith C, Smyth G, Tierney L, Yang JYH, Zhang J. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5:R80.View ArticlePubMedPubMed CentralGoogle Scholar
- Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3.PubMedGoogle Scholar
- Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995;57:289–300.Google Scholar
- Brunet J-P, Tamayo P, Golub TR, Mesirov JP. Metagenes and molecular pattern discovery using matrix factorization. Proc Natl Acad Sci U S A. 2004;101:4164–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Gaujoux R, Seoighe C. A flexible R package for nonnegative matrix factorization. BMC Bioinformatics. 2010;11:367.View ArticlePubMedPubMed CentralGoogle Scholar
- Reich M, Liefeld T, Gould J, Lerner J, Tamayo P, Mesirov JP. GenePattern 2.0. Nat Genet. 2006;38:500–1.View ArticlePubMedGoogle Scholar
- Benito M, Parker J, Du Q, Wu J, Xiang D, Perou CM, Marron JS. Adjustment of systematic microarray data biases. Bioinformatics. 2003;20:105–14.View ArticleGoogle Scholar
- Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–50.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang J, Duncan D, Shi Z, Zhang B. WEB-based GEne SeT AnaLysis Toolkit (WebGestalt): update 2013. Nucleic Acids Res. 2013;41(Web Server issue):W77–83.View ArticlePubMedPubMed CentralGoogle Scholar
- Singh A, Greninger P, Rhodes D, Koopman L, Violette S, Bardeesy N, Settleman J. A gene expression signature associated with “K-Ras addiction” reveals regulators of EMT and tumor cell survival. Cancer Cell. 2009;15:489–500.View ArticlePubMedPubMed CentralGoogle Scholar
- Janky R, Verfaillie A, Imrichová H, Van de Sande B, Standaert L, Christiaens V, Hulselmans G, Herten K, Naval Sanchez M, Potier D, Svetlichnyy D, Kalender Atak Z, Fiers M, Marine J-C, Aerts S. iRegulon: from a gene list to a gene regulatory network using large motif and track collections. PLoS Comput Biol. 2014;10:e1003731.View ArticlePubMedPubMed CentralGoogle Scholar
- Hoskins JW, Jia J, Flandez M, Parikh H, Xiao W, Collins I, Emmanuel M a, Ibrahim A, Powell J, Zhang L, Malats N, Bamlet WR, Petersen GM, Real FX, Amundadottir LT. Transcriptome analysis of pancreatic cancer reveals a tumor suppressor function for HNF1A. Carcinogenesis. 2014;35:2670–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Nakamura T, Furukawa Y, Nakagawa H, Tsunoda T, Ohigashi H, Murata K, Ishikawa O, Ohgaki K, Kashimura N, Miyamoto M, Hirano S, Kondo S, Katoh H, Nakamura Y, Katagiri T. Genome-wide cDNA microarray analysis of gene expression profiles in pancreatic cancers using populations of tumor cells and normal ductal epithelial cells selected for purity by laser microdissection. Oncogene. 2004;23:2385–400.View ArticlePubMedGoogle Scholar
- Stathis A, Moore MJ. Advanced pancreatic carcinoma: current treatment and future challenges. Nat Rev Clin Oncol. 2010;7:163–72.View ArticlePubMedGoogle Scholar
- Eser S, Schnieke A, Schneider G, Saur D. Oncogenic KRAS signalling in pancreatic cancer. Br J Cancer. 2014;111:817–22.View ArticlePubMedPubMed CentralGoogle Scholar
- Cheng W-Y, Ou Yang T-H, Anastassiou D. Biomolecular events in cancer revealed by attractor metagenes. PLoS Comput Biol. 2013;9:e1002920.View ArticlePubMedPubMed CentralGoogle Scholar
- Fink K, Grandvaux N. STAT2 and IRF9: beyond ISGF3. Jak-Stat. 2013;2:e27521.View ArticlePubMedPubMed CentralGoogle Scholar
- Postigo AA, Dean DC. ZEB represses transcription through interaction with the corepressor CtBP. Proc Natl Acad Sci U S A. 1999;96:6683–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Bronsert P, Kohler I, Timme S, Kiefer S, Werner M, Schilling O, Vashist Y, Makowiec F, Brabletz T, Hopt UT, Bausch D, Kulemann B, Keck T, Wellner UF. Prognostic significance of Zinc finger E-box binding homeobox 1 (ZEB1) expression in cancer cells and cancer-associated fibroblasts in pancreatic head cancer. Surgery. 2014;156:97–108.View ArticlePubMedGoogle Scholar
- Grisanzio C, Werner L, Takeda D, Awoyemi BC, Pomerantz MM, Yamada H, Sooriakumaran P, Robinson BD, Leung R, Schinzel AC, Mills I, Ross-Adams H, Neal DE, Kido M, Yamamoto T, Petrozziello G, Stack EC, Lis R, Kantoff PW, Loda M, Sartor O, Egawa S, Tewari AK, Hahn WC, Freedman ML. Genetic and functional analyses implicate the NUDT11, HNF1B, and SLC22A3 genes in prostate cancer pathogenesis. Proc Natl Acad Sci U S A. 2012;109:11252–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Berndt SI, Sampson J, Yeager M, Jacobs KB, Wang Z, Hutchinson A, Chung C, Orr N, Wacholder S, Chatterjee N, Yu K, Kraft P, Feigelson HS, Thun MJ, Diver WR, Albanes D, Virtamo J, Weinstein S, Schumacher FR, Cancel-Tassin G, Cussenot O, Valeri A, Andriole GL, Crawford ED, Haiman C, Henderson B, Kolonel L, Le Marchand L, Siddiq A, Riboli E, et al. Large-scale fine mapping of the HNF1B locus and prostate cancer risk. Hum Mol Genet. 2011;20:3322–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Cuff J, Salari K, Clarke N, Esheba GE, Forster AD, Huang S, West RB, Higgins JP, Longacre TA, Pollack JR. Integrative bioinformatics links HNF1B with clear cell carcinoma and tumor-associated thrombosis. PLoS One. 2013;8:e74562.View ArticlePubMedPubMed CentralGoogle Scholar
- Worley MJ, Welch WR, Berkowitz RS, Ng S-W. Endometriosis-associated ovarian cancer: a review of pathogenesis. Int J Mol Sci. 2013;14:5367–79.View ArticlePubMedPubMed CentralGoogle Scholar
- Haumaitre C, Fabre M, Cormier S, Baumann C, Delezoide A-L, Cereghini S. Severe pancreas hypoplasia and multicystic renal dysplasia in two human fetuses carrying novel HNF1beta/MODY5 mutations. Hum Mol Genet. 2006;15:2363–75.View ArticlePubMedGoogle Scholar
- Harrell JC, Dye WW, Harvell DME, Sartorius CA, Horwitz KB. Contaminating cells alter gene signatures in whole organ versus laser capture microdissected tumors: a comparison of experimental breast cancers and their lymph node metastases. Clin Exp Metastasis. 2008;25:81–8.View ArticlePubMedGoogle Scholar
- Moffitt R a, Marayati R, Flate EL, Volmar KE, Loeza SGH, Hoadley K a, Rashid NU, Williams L a, Eaton SC, Chung AH, Smyla JK, Anderson JM, Kim HJ, Bentrem DJ, Talamonti MS, Iacobuzio-Donahue C a, Hollingsworth M a, Yeh JJ. Virtual microdissection identifies distinct tumor- and stroma-specific subtypes of pancreatic ductal adenocarcinoma. Nat Genet. 2015;47:1168–78.View ArticlePubMedPubMed CentralGoogle Scholar
- Donahue TR, Tran LM, Hill R, Li Y, Kovochich A, Hargan Calvopina J, Patel SG, Wu N, Hindoyan A, Farrell JJ, Li X, Dawson DW, Wu H. Integrative Survival-Based Molecular Profiling of Human Pancreatic Cancer. Clin Cancer Res. 2012;18:1352–1363.View ArticlePubMedGoogle Scholar
- Jiang X, Zhang W, Kayed H, Zheng P, Giese NA, Friess H, Kleeff J. Loss of ONECUT1 expression in human pancreatic cancer cells. Oncol Rep. 2008;19:157–63.PubMedGoogle Scholar
- Wang L, Coffinier C, Thomas MK, Gresh L, Eddu G, Manor T, Levitsky LL, Yaniv M, Rhoads DB. Selective deletion of the Hnf1beta (MODY5) gene in beta-cells leads to altered gene expression and defective insulin release. Endocrinology. 2004;145:3941–9.View ArticlePubMedGoogle Scholar
- Listing H, Mardin W a, Wohlfromm S, Mees ST, Haier J. MiR-23a/-24-induced gene silencing results in mesothelial cell integration of pancreatic cancer. Br J Cancer. 2014;112:131–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Roy N, Malik S, Villanueva KE, Urano A, Lu X, Figura V, Seeley ES, Dawson DW, Collisson E a, Hebrok M, Von Figura G, Seeley ES, Dawson DW, Collisson E a, Hebrok M. Brg1 promotes both tumor-suppressive and oncogenic activities at distinct stages of pancreatic cancer formation. Genes Dev. 2015;29:658–71.View ArticlePubMedPubMed CentralGoogle Scholar
- Harries LW, Brown JE, Gloyn AL. Species-specific differences in the expression of the HNF1A, HNF1B and HNF4A genes. PLoS One. 2009;4:e7855.View ArticlePubMedPubMed CentralGoogle Scholar