Multidimensional analyses identify genes of high priority for pancreatic cancer research

Abstract

Pancreatic ductal adenocarcinoma (PDAC) is a drug-resistant and lethal cancer. Identification of the genes that consistently show altered expression across patient cohorts can expose effective therapeutic targets and strategies. To identify such genes, we separately analyzed 5 human PDAC microarray datasets. We defined genes as "consistent" if upregulated or downregulated in 4 or more datasets (adjusted P < 0.05). The genes were subsequently queried in additional datasets, including single-cell RNA-sequencing data, and we analyzed their pathway enrichment, tissue specificity, essentiality for cell viability, and association with cancer features, e.g., tumor subtype, proliferation, metastasis, and poor survival outcome. We identified 2,010 consistently upregulated and 1,928 downregulated genes, of which more than 50% to our knowledge were uncharacterized in PDAC. These genes spanned multiple processes, including cell cycle, immunity, transport, metabolism, signaling, and transcriptional/epigenetic regulation - cell cycle and glycolysis being the most altered. Several upregulated genes correlated with cancer features, and their suppression impaired PDAC cell viability in prior CRISPR/Cas9 and RNA interference screens. Furthermore, the upregulated genes predicted sensitivity to bromodomain and extraterminal (epigenetic) protein inhibition, which, in combination with gemcitabine, disrupted amino acid metabolism and in vivo tumor growth. Our results highlight genes for further studies in the quest for PDAC mechanisms, therapeutic targets, and biomarkers.

Keywords: Cancer; Gastroenterology; Glucose metabolism; Molecular genetics; Oncology.

Figure 1. Consistently upregulated or downregulated genes in human pancreatic tumors.

(A) Schematic overview illustrating the identification and potential utility of the consistent genes. Five PDAC microarray datasets were used for the identification of the genes and “consistent” was defined as genes upregulated or downregulated in at least 4 datasets (adjusted P < 0.05). CUGs, consistently upregulated genes; CDGs, consistently downregulated genes. See Methods (Identification of the consistent genes) or Supplemental Figure 1B for the sample size of each dataset. *Not the focus of this study. (B) Topmost 20 highly and lowly expressed genes in PDAC (i.e., CUGs and CDGs, respectively) based on the sum of expression rank in the 5 datasets. (C) Number of CUGs and CDGs that showed the same high or low expression pattern, respectively, in the independent microarray datasets GSE19279, GSE32676, GSE19650, and GSE62165. *Dataset of premalignant tumor stages. See Supplemental Methods for sample sizes/types of the independent datasets. (D) Pie chart indicating the distribution of consistent genes (i.e., CUGs and CDGs) in PDAC relative to the number (n) of prior publications on the genes (0, ≤5, and >5) as observed via PubMed search. Highlighted in red are topmost CUGs with no prior publication; in green are topmost CDGs with no prior publication. (E) Gene set enrichment analysis (GSEA) plots of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways or Hallmark associated with the consistent genes. The plots were generated using only the 3,938 genes that are consistent (adjusted P < 0.05 in tumor versus nontumor comparison in at least 4 of 5 datasets). The gene list used for GSEA was ranked by the sum of the expression score across the 5 datasets. NES, normalized enrichment score.

Figure 2. Consistent genes are components of multiple pathways.

(A) Heatmap showing cell cycle components with consistent expression pattern. (B and C) Gene set enrichment analysis (GSEA) plots indicating upregulated (B) IFN-α and (C) IFN-γ signatures. The listed genes are the most highly upregulated in the pathways shown. (D–F) Schematic depictions of consistent genes in the (D) glycolytic pathway and the tricarboxylic acid (TCA) cycle, (E) fatty acid (FA) biosynthesis and breakdown (β-oxidation), and (F) cholesterol metabolism. (G) Heatmap showing the most consistently expressed nutrient and ion transporters (top 20 up-/downregulated). (H) Heatmap showing the most consistently expressed transcription factors (top 20 up-/downregulated). (I and J) GSEA plots of transcriptional signature, indicating (I) upregulated BACH1_01 signature and (J) downregulated AR_01 and HNF_Q6 signatures. The genes in red are among the topmost of the enriched BACH1_01 signature, whereas genes in green are among the topmost downregulated in the AR_01 and HNF_Q6 plots, respectively. In D and E, *indicates genes that were consistent in 3 of 5 datasets; all other genes displayed are consistent in 4 or more datasets (in red – CUGs, in green – CDGs). Genes in gray were not captured as consistent and are shown based on their known position/role in the respective pathways. NES, normalized enrichment score.

Figure 3. Consistent genes show a tumor-specific or non–tumor-specific expression pattern.

(A) Venn diagram showing overlap between CUGs and genes with low expression in normal pancreas and/or highly expressed in epithelium compared with stroma from the laser microdissection dataset (GSE93326). Right: Venn diagram showing overlap between CDGs and genes highly expressed in normal pancreas and/or highly expressed in the stroma (low in epithelium). Gene expression in the normal pancreas was determined by comparing pancreatic tissue gene expression to that of other normal tissues in the Human Protein Atlas (HPA, n = 42 other tissue types), Genotype-Tissue Expression (GTEx, n = 33 other tissue types), and GSE71729 (Moffitt) datasets (normal tissues: pancreas, n = 46; liver, n = 27; lung, n = 19; lymph node, n = 10; spleen, n = 11). Genes that showed a similar expression pattern in at least 2 of HPA, GTEx, or GSE71729 were selected. (B and C) Heatmaps showing (B) topmost tumor-specific CUGs differentially expressed in the tumor epithelium samples, and (C) tumor-specific CDGs that show low expression in the tumor epithelium samples of the laser capture microdissection dataset GSE93326 (n = 65 pairs of tumor epithelium and stroma samples). (D) Pathway annotation of the 564 tumor-specific CUGs (high expression in the epithelium and low expression in normal pancreas). (E and F) Gene ontology for cellular component of the (E) 564 core CUGs and (F) the 250 core CDGs.

Figure 4. scRNA-seq data showing the expression of the consistent genes in tumor and surrounding cell population.

(A) Uniform manifold approximation and projections (UMAPs) of single-cell RNA-sequencing (scRNA-seq) data of human PDAC, depicting various cell populations and cell type marker plots. KRT19, keratin 19; KRT8, keratin 8 (pancreatic ductal epithelial markers); AMY2A, amylase α2A; CTRB2, chymotrypsinogen B2 (pancreatic acinar cell markers); FOXP3, forkhead box P3 (regulatory T cell [Treg] marker); GZMB, granzyme B (cytotoxic T/natural killer [NK] cell marker); ACTA2, actin α2 smooth muscle; COL1A1, collagen type I α1 chain (fibroblast markers); CD14, cluster of differentiation 14; APOE, apolipoprotein E; C1QA; complement C1q A chain (macrophage/myeloid cell markers); HBA2, hemoglobin subunit α2; HBB, hemoglobin subunit β (red blood cell markers). Sample size: n = 61 primary tumors. (B) UMAPs showing cell populations expressing CUGs that showed tumor-specific upregulation in microarrays and laser capture microdissection dataset (in Figure 3) and (C) CDGs that showed tumor-specific downregulation. (D) UMAPs showing consistently upregulated glycolysis genes that emerged as more tumor-specific (upper row, except PFKM) or ubiquitously expressed in most other cell types (lower row). On the right is a schematic summary showing the glycolysis steps associated with the displayed genes. *Consistent in 3 of 5 datasets shown in Figure 1.

Figure 5. Consistent genes correlate with poor cancer prognostic features.

(A) Venn diagram showing overlap between CUGs or CDGs and genes expressed by PDAC basal-like or classical subtypes. The basal-like versus classical subtype genes included are statistically upregulated or downregulated (P < 0.05) in 2 or more datasets, namely, TCGA (n = 31 basal-like vs. 31 classical), GSE71729 (Moffitt dataset, n = 27 basal-like vs. 27 classical), and Puleo (n = 64 basal-like vs. 64 classical tumor samples). Genes in red/green are selected examples. (B) Heatmap depicting topmost upregulated or downregulated genes in proliferation-high PDAC. *Indicates proliferation markers used for the tumor stratification and included as positive controls. Genes included are statistically up-/downregulated in proliferation-high relative to proliferation-low tumors in 2 or more of TCGA (n = 64 proliferation-high vs. 86 proliferation-low), GSE71729 (n = 77 proliferation-high vs. 68 proliferation-low), and Puleo (n = 99 proliferation-high vs. 210 proliferation-low) tumors (P < 0.05). Pro, proliferation. (C) Pathway annotation of 783 CUGs and 588 CDGs that overlapped with genes differentially expressed in proliferation-high vs. -low tumors in at least 2 datasets. cAMP, cyclic adenosine monophosphate; cGMP, cyclic guanosine monophosphate; PKG, protein kinase G. (D) Venn diagrams showing CUGs and, below, CDGs overlapping with genes in liver metastasis compared with normal liver tissues from GSE71729 and GSE19279 datasets (P < 0.05). See Supplemental Methods for sample sizes/types. (E) Pathway annotation of the CUGs and CDGs that overlapped between liver metastasis compared with normal liver tissues. PI3K, phosphoinositide 3-kinase; ECM, extracellular matrix. (F) Kaplan-Meier (KM) overall survival (OS) plots (log-rank test, P < 0.05) of genes that predicted survival in the clinical cohorts analyzed. Tumor sample size: TCGA (n = 146), GSE71729 (n = 125), Puleo (n = 288); and ICGC (n = 267). (G) Topmost genes that predicted OS in at least 3 of the 4 pancreatic adenocarcinoma datasets, i.e., TCGA, Puleo et al., GSE71729/Moffitt, and International Cancer Genome Consortium (IGCG) datasets. High expression of the genes in green predicts “better” outcome, whereas those in red predict “worse” outcome.

Figure 6. High-priority consistent genes.

(A) Workflow for the analysis of the CRISPR/Cas9 and shRNA (collectively gene interference) screen data to identify CUGs that are potentially essential for PDAC cell survival or growth. (B) Venn diagram showing the overlapping essential genes in PDAC cell lines as derived from the Project Drive, Project Achilles, Behan et al., and GECKO screen data. In bold are the number of CUGs, in total 185, that overlapped as essential for survival in at least 3 of the 4 gene interference screen data. Right: Selected genes among the 185 high-priority targets. (C–F) Pathway annotation of genes that emerged as essential for PDAC survival/growth in (C) Project Achilles (n = 394 CUGs), (D) GECKO (n = 231 CUGs), (E) Behan (n = 211 CUGs), and (F) Project Drive (n = 200 CUGs). OM, oocyte maturation; HTLV, human T lymphotropic virus.

Figure 7. Expression of high-priority genes predicts therapy response.

(A) Heatmap showing the sensitivity of PDAC cell lines expressing high-priority genes to compounds tested in the Genomics of Drug Sensitivity in Cancer (GDSC2) data. In blue font are epigenetics drugs. EXP, expressing. (B) Viability assay of PDAC cell lines treated with BETi AZD5153 alone or in combination with gemcitabine (gem) for 72 hours, in quadruplicate. BETi, bromodomain and extraterminal motif (BET) protein inhibitor(s); DMSO, dimethyl sulfoxide. (C) Heatmaps showing the metabolomics profile of lysate (intracellular) and culture media (extracellular) extracts from human PDAC cell TU8902 treated with gemcitabine, AZD5153, or the combination, for 24 hours, in triplicate. AMP, adenosine monophosphate; NAD, nicotinamide adenine dinucleotide. (D) Viability assay of TU8902 treated with gemcitabine, BD-9136 (a bromodomain-containing protein 4 degrader), or both, for 48 hours. Represents more than 2 independent experiments. Below: Treatment with BD-9136 or AZD5153 alone or in combination with gemcitabine for 48 hours, in quadruplicate. (E) Workflow of the mouse tumor experiments. Mice were treated with gemcitabine at 100 mg/kg body weight 2 times per week and with BD-9136 at 20 mg/kg body weight 5 times per week. (F) Tumor weight of subcutaneously implanted TU8902 (n = 4 mice per group, injected on both flanks; n = 8 tumors per group except combination arm, n = 6). Right: Image depicting the size of tumors harvested at the end of the experiment. (G) Image depicting the size of KPC 7940b pancreatic orthotopic tumors following the indicated treatments (n = 6 mice per group). Below: Bar graph showing tumor weights. Statistical comparison of the indicated groups (F and G) was by 2-tailed t test with Welch’s correction. Comparison by 1-way ANOVA with Tukey’s post hoc correction was not statistically significant. Comparisons for B and D (below) were by unpaired t test and D was by 1-way ANOVA with Tukey’s post hoc test. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001 (B, D, F, and G).