Identifying drug candidates for pancreatic ductal adenocarcinoma based on integrative multiomics analysis

Introduction

Pancreatic ductal adenocarcinoma (PDAC) is the most common exocrine pancreatic cancer, accounting for approximately 95% of all pancreatic cancers. PDAC, which has the poorest 5-year survival rate of approximately 10% in USA, is a highly invasive tumor formed by pancreatic tubule lining cells (1). The limited number of treatments and early detection methods are main reasons for the poor prognosis (2).

It is generally believed that PDAC is caused by the accumulation of mutations in genes, especially tumor suppressor genes KRAS, leading to the development of invasive cancer from pancreatic intraepithelial neoplasia (3). The genetic and phenotypic heterogeneity of PDAC makes it quite complicated to formulate a generally effective therapy. Through whole-genome expression profiling, PDAC can be divided into 3–4 molecular subtypes, which can be distinguished by unique molecular and phenotypic characteristics (4,5). Therefore, early detection of pancreatic cancer is urgently needed and can be accomplished by the identification of effective early diagnostic markers. Epigenetic research has confirmed that changes in DNA methylation have a critical relationship with the ontogenesis and progression of various diseases by regulating gene expression without altering the DNA sequence and early tumorigenesis, which provides us with an opportunity to screen reliable and effective markers for early diagnosis and prognosis evaluation (6). Hypermethylation in cytosine-phosphate-guanine (CpG) islands of tumor suppressor gene promoters may suppress gene transcription, which could cause dysfunction of the affected tumor suppressor genes and allow unrestricted cell growth and tumor occurrence (7). DNA methylation modification could arise at the early phase of tumorigenesis and last throughout tumor progression, which indicates that it could be a reliable biomarker for the early diagnosis of cancer and the prediction of patient prognosis (6). Copy number variations (CNVs) are also frequent genomic alteration events. A previous study has shown a close correlation between CNV and dysregulated gene expression in many cancer types (8). Integrative multiomics analysis makes it possible to explore diagnostic tumor markers.

The unsatisfactory therapeutic effect of drugs is one of the main reasons for the poor prognosis of PDAC. At present, gemcitabine plus nab-paclitaxel (GnP) and FOLFIRINOX are still the standard treatments for PDAC, but these regimens offer a limited advantage in overall survival due to drug resistance (9). The development of new anticancer drugs is time-consuming, expensive and low-yield, so screening and repurposing of existing drugs is an alternative and attractive approach (10). The Drug-Gene Interaction database (DGIdb, https://www.dgidb.org) can be used to identify plausible drug candidates through drug-gene interactions (11).

This study screened a group of differentially methylated regions (DMRs) using methylation expression data of PDAC from The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) datasets. Some significant DMRs were selected to build a diagnostic model for PDAC. Then, some hub genes, which might facilitate important molecular mechanisms in oncogenesis, were obtained by combining the transcript profile and CNV data. We further filtered out some potential therapeutic drugs for PDAC using interaction information related to existing drugs and genes. We present this article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-23-985/rc).

Methods

Data from the TCGA and ICGC cohorts

The RNA-sequencing (RNA-seq) profiles, methylation expression data, CNV data and corresponding clinical information of PDAC were downloaded from the Genomic Data Commons Application Programming Interface of TCGA (https://cancergenome.nih.gov/), containing 184 tumor and 10 normal cases as of March 2020. Methylation expression data of PDAC were also downloaded from ICGC (https://icgc.org/), which contained 261 tumor and 25 normal cases as of July 2013.

Data preprocessing

For methylation expression data from Human Methylation 450k, the missing value was replaced using the K-nearest neighbor (k-NN) algorithm. The probe was removed if the missing value of a probe exceeded 50% of the total samples; otherwise, it was filled with the mean value of ten adjacent probes. For RNA-seq data, genes with low expression levels, defined as an expression value of fragments per kilobase million (FPKM) <1 in more than 90% of samples, were removed. Samples with no survival time or missing event values were also removed.

Differential methylation analysis

Differential methylation analysis was performed using datasets from TCGA and ICGC, which had the same test platform. DMR was regarded as a false discovery rate (FDR) <0.05. To ensure the reliability of filtered methylation sites, only sites with fold change ≥2 for methylation difference level and average value of tumor sample ≥0.3 were filtered. Differences in methylation levels and similarity between normal and tumor tissues were both compared.

Correlation analysis of methylation and gene expression

Considering that gene expression and methylation and CNV had a causal relationship, genes associated with DMRs were filtered down using a regression model. As this method incorporates CNV and other potential factors that affect gene expression, the correlation analysis was more reasonable. Here, we established the following model: G_expr = α*G_methy + β*G_CNV + ε, in which G_expr represents the gene expression level, G_methy is the gene methylation level, G_CNV is the gene CNV status, and α and β are coefficients. ε was used as an error term because gene expression might be affected by the environment, quantitative trait loci, and upstream and downstream regulatory genes. According to the model, parameters were selected using least absolute shrinkage and selection operator (LASSO) regression, which used the L1-norm as the regular term of the loss function and could optimize the model by removing parameters with a coefficient of 0.

We filtered the results according to the following criteria: (I) if α≠0 and β=0, then preserved; (II) if α≠0 and β≠0, then the gene with coefficient α/β≥4 was preserved. Genes meeting both criteria were used for subsequent analysis. In addition, we also evaluated the correlation between gene methylation and its expression levels to further verify the results.

Construction of a diagnostic model with DMRs

To make the sample size larger, we merged the methylation data of TCGA and ICGC and then eliminated the batch effect. According to the categories of the samples, the random forest algorithm with the R package random forest and LASSO regression model were used to perform further dimensional reduction on the DMRs.

According to the random forest algorithm, the random variables of mtry and ntree parameters were set to 1–50 and 500, respectively, for each split. First, the mtry and ntree values with the lowest error rates were selected as the optimal values. Then, DMRs were sorted according to the importance of the results. For the LASSO regression model, 10-fold cross-validation was performed to determine the optimal lambda value, and DMRs with the lowest error rate were selected. DMRs shared by both algorithms were finally selected and then randomly divided into a training set and validation set in equal proportions. Then, a diagnostic model for PDAC was built with the logistic regression method. To avoid the bias effect of single sample grouping on the model, we conducted 100 random groupings and evaluated the area under the curve (AUC) value of the model under each grouping.

Marker gene characteristics test of the diagnostic model

The relationship of the methylation level of the selected DMRs and their corresponding gene expression levels was tested in tumorous and adjacent normal samples. The expression levels of PHF12, FXYD3, PRKCB and ZNF582, regulated by the four corresponding DMRs, were also analyzed in different tumor stages. Then, patients in TCGA were divided into high and low groups according to the median value of each gene expression, and prognostic survival analysis was performed to test their prognostic predictive values.

Coexpression network analysis

Pearson correlation coefficients of all selected genes and the four genes were calculated, and a correlation test was also conducted with preprocessed expression profile data in TCGA. Genes with FDR <0.001 were defined as coexpressed genes of the four genes. Weighted gene coexpression network analysis (WGCNA), which can be used to find highly relevant gene clusters and correlations among modules and samples, was used to further study the correlation patterns among the coexpressed genes. Before conducting WGCNA, abnormal samples were evaluated and removed in subsequent analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the genes was performed using the R package clusterProfiler.

Drug-gene interaction analysis

DGIdb, which provides information on the interactions between drugs and genes, was used to screen potential drugs for PDAC by searching with the four hub genes.

Immunohistochemical (IHC) validation of the protein expression of model-related genes

Thirty PDAC patients who underwent curative resection with no neoadjuvant therapy before surgery at The First Affiliated Hospital of Zhengzhou University (Zhengzhou, China) between February 2019 and November 2021 participated in the retrospective investigation. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013) and approved by ethics committee of The First Affiliated Hospital of Zhengzhou University (No. 2023-KY-0863-002). Informed consent was waived due to the retrospective nature of the study.

All hematoxylin and eosin-stained slides of tumor and adjacent nontumor samples were collected and confirmed by two experienced pathologists. Three genes, PHF12, FXYD3 and PRKCB, were selected for further validation by IHC method as previously described (12). All formalin-fixed paraffin-embedded PDAC samples were chosen to test the protein levels of the three genes.

Slides with 4 µm sections were prepared for deparaffinizing and blocking endogenous peroxidase activity. After washing the slides with phosphate buffered saline and incubating them with bovine serum, further incubation with primary and secondary antibodies were followed. All sections were counterstained with hematoxylin after incubating with 3,3'-diaminobenzidine solution. The IHC slides were analyzed by Image-Pro Plus 6.0 software (Media Cybernetics, Rockville, MD, USA). The results were recorded by two pathologists. For PHF12, FXYD3 and PRKCB expression analysis, relevant primary antibodies (PHF12, SC-514864, 200 µg/mL, Santa Cruz, CA, USA; FXYD3, AB205534, 1:1,000, Abcam, Cambridge, UK; PRKCB, BS-0267R, 1:1,000, Bioss, Beijing, China) were used. The results are expressed as the mean areal density (MAD).

Western blotting evaluation of the protein expression of model-related genes

To further evaluate the gene expression of PHF12, FXYD3 and PRKCB, another six pairs of fresh-frozen PDAC and adjacent normal samples were collected and tested using Western blotting.

Proteins were isolated from fresh-frozen tissues and processed as described earlier (12). The proteins were treated with primary rabbit polyclonal antibodies (PHF12, SC-514864, 200 µg/mL, Santa Cruz; FXYD3, AB205534, 1:1,000, Abcam; PRKCB, BS-0267R, 1:1,000, Bioss) and GAPDH (GB12002, 1:1,000, Servicebio, Wuhan, China), then with a horseradish peroxidase-conjugated goat anti-rabbit secondary antibody (GB23303, 1:3,000, Servicebio). Finally, the bound immunocomplexes were detected by enhanced chemiluminescence methods.

Statistical analysis

The GeneCodis tool (http://genecodis.cnb.csic.es/) was used for gene annotation, including Gene Ontology (GO), biological process (BP), cellular component (CC), molecular function (MF) and KEGG pathway enrichment analyses. Hyp_c <0.05 was used as the threshold for the significant enrichment of GO terms and KEGG pathways. All analyses were performed using R (version 3.6.1) software, and chip data preprocessing and probe annotation were performed using Bioconductor related packages. A t-test was used for paired data, and a U test was used for unpaired data. Values with P<0.05 were regarded as statistically significant.

Results

Data preprocessing

The analysis procedure for this study is shown as a flowchart in Figure 1. After merging the TCGA and ICGC datasets and eliminating batch effects, differences between the two datasets of tumor samples and normal samples were eliminated, indicating that the merged dataset had low bias for subsequent analysis (Figure S1). A total of 46 DMRs with cumulative importance >95% were used for subsequent analysis.

Figure 1 Flowchart describing the procedure of building a diagnostic model and exploring potential therapeutic drugs by comprehensive multiomics analysis of PDAC. TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; DMRs, differentially methylated regions; LASSO, least absolute shrinkage and selection operator; RNA-seq, RNA-sequencing; CNV, copy number variation; WGCNA, weighted gene coexpression network analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; DGIdb, Drug-Gene Interaction database; PDAC, pancreatic ductal adenocarcinoma; IHC, immunohistochemical.

Identification of DMRs among TCGA and ICGC groups

The comparison of tumor and normal samples from the TCGA and ICGC yielded 7,628 and 9,463 DMRs, respectively (supplementary file 1 available at https://cdn.amegroups.cn/static/public/jgo-23-985-1.xls), with high internal consistency (Figure S2A,S2B). The methylation trend indicated that there were more hypermethylation sites than hypomethylation sites in the tumor group in both the TCGA and ICGC datasets (Figure 2A,2B). According to the distribution of these DMRs on the genome, they were mainly concentrated on the island, N shore, and S shore but were rarely distributed on the N shelf and S shelf (Figure S2C,S2D). There were 4,164 DMRs, including 3,624 hypermethylated regions and 540 hypomethylated regions, that were shared by the two datasets, and these were rarely distributed on the N shelf and S shelf (Figure 2C,2D).

Figure 2 Identification and distribution of DMRs in TCGA and ICGC cohorts, which had 7,628 and 9,463 DMRs, respectively. (A) Different methylation levels between tumor and adjacent normal tissues of the TCGA cohort. (B) Different methylation levels between tumor and adjacent normal tissues of the ICGC cohort. (C) 4,164 DMRs shared by the TCGA and ICGC cohorts, with 3,624 hypermethylated regions and 540 hypomethylated regions. (D) Distribution of DMRs in the genome, mainly concentrated on the island, N shore, and S shore. DMRs, differentially methylated regions; Hyper-M, hypermethylated; Hypo-M, hypomethylated; ICGC, International Cancer Genome Consortium; TCGA, The Cancer Genome Atlas; PDAC, pancreatic ductal adenocarcinoma.

Potential regulatory genes for DMRs

The TCGA dataset provided PDAC gene expression and CNV data. We used three omics datasets to analyze 4,164 genes possibly regulated by DMRs and finally obtained 285 genes involving 589 DMRs (supplementary file 2 available at https://cdn.amegroups.cn/static/public/jgo-23-985-2.xls). The expression levels of these genes were significantly negatively correlated with the methylation levels (Figure 3A,3B). The significance test of the correlation coefficient only had 34 genes with P values ≥0.05, and the other 555 genes had P values <0.05 (supplementary file 2 available at https://cdn.amegroups.cn/static/public/jgo-23-985-2.xls), which showed that the method had better performance. KEGG pathway enrichment analysis showed that the 285 potential regulatory genes were mainly involved in neuroactive ligand-receptor interactions, pathways in cancer, the PI3K-Akt signaling pathway, the MAPK signaling pathway, the calcium signaling pathway, and other components. These pathways were closely related to the cancer signaling process (Figure 3C).

Figure 3 The DMRs and potential regulatory genes in TCGA samples. (A) Expression levels of 589 DMRs in the TCGA samples. (B) Expression levels of 285 potential genes regulated by the DMRs. (C) KEGG pathway enrichment analysis of 285 potential regulated genes. *, P<0.05, **, P<0.01, ***, P<0.001. FPKM, fragments per kilobase million; KEGG, Kyoto Encyclopedia of Genes and Genomes; DMRs, differentially methylated regions; TCGA, The Cancer Genome Atlas.

Diagnostic model constructed with DMRs

After merging the TCGA and ICGC datasets, we randomly divided the tumor and normal samples into a training set (tumor/normal: 222/18) and a validation set (tumor/normal: 223/17). Random forest and LASSO methods were used to obtain 46 and 15 DMRs, respectively (Figure 4A,4B). cg10547050, cg08823209, cg03306374, and cg09568464, shared by both methods, were used to construct a PDAC diagnostic risk score model with the following formula: (−10.86 × expression level of cg10547050) + (−10.38 × expression level of cg08823209) + (7.18 × expression level of cg03306374) + (12.02 × expression level of cg09568464). The results showed that the risk scores of tumor samples were significantly higher than those of normal samples (Figure 4C). Furthermore, the ability of the model to classify tumor and normal samples was also evaluated. The results showed that the model could strongly classify the two groups of samples in both the training set and validation set, with AUC values reaching 0.984 and 0.988, respectively (Figure 4D).

Figure 4 Construction and validation of the pancreatic cancer diagnostic model. (A) Two-dimensional distribution diagram of adjacency matrix among samples by random forest analysis in the training set. (B) Relationship between the number of different variables and the classification error rate in the LASSO regression model. (C) Tumor samples (T) have a significantly higher risk score than normal samples (N) tested with the diagnostic model in both the TCGA dataset and the ICGC dataset. (D) ROC curve analysis of the diagnostic model for the training set and validation set. ***, P<0.001. Dim, dimension; ICGC, International Cancer Genome Consortium; TCGA, The Cancer Genome Atlas; AUC, area under the curve; CI, confidence interval; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

Analysis of marker genes regulated by DMRs in the diagnostic model

The expression of the four DMRs in tumor and normal tissues was analyzed. The results indicated that the levels of cg03306374 and cg09568464 were both significantly higher in tumorous tissue than in normal tissue, and the other two DMRs showed the opposite pattern (Figure 5A,5B). We further analyzed the methylation levels of all methylation sites within the upstream 2 kb transcription start site (TSS) of the four genes regulated by the DMRs. The results showed that these sites all had a high degree of methylation, which was consistent with the methylation trend of the four DMRs. Although only the expression of FXYD3 was significantly different between tumor and normal tissues (P<0.05), the expression patterns of the other three DMR-regulated genes, PHF12, PRKCB, and ZNF582, were also consistent with their methylation patterns (Figure 5C). Correlation tests for the levels of the four DMRs and the corresponding gene expression levels in tumorous tissue all revealed significant negative results (P<0.05) (Figure 5D). In the analysis of gene expression in different stages of PDAC, only the expression of ZNF582 was significantly lower in stage IV than in the early stages (P<0.05), and there was a significant correlation between ZNF582 expression and prognosis (Figure 5E). Survival analysis indicated that the group with low expression of ZNF582 had a worse prognosis (P<0.05), but the other three genes did not have significant prognostic values (P>0.05) (Figure 5F).

Figure 5 Expression of the 4 DMRs and corresponding regulatory genes in the diagnostic model. Expression levels of the 4 DMRs in tumor and normal samples from the TCGA dataset (A) and ICGC dataset (B). (C) Expression levels of 4 hub DMR-regulated genes in tumor samples and normal samples from GEPIA2 (http://gepia2.cancer-pku.cn/). (D) Correlation between the expression levels of 4 DMRs and the corresponding gene expression levels in tumorous tissue from the TCGA dataset. (E) Expression levels of 4 hub genes in different stages of PDAC. (F) Prognostic survival analysis of 4 hub genes by dividing samples into high and low groups according to the median value of each gene expression in TCGA. *, P<0.05; **, P<0.01; ***, P<0.001. TPM, transcripts per million; FPKM, fragments per kilobase million; DMRs, differentially methylated regions; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; PDAC, pancreatic ductal adenocarcinoma.

Coexpression network and functional analysis of the genes

As shown in the supplementary file 3 (available at https://cdn.amegroups.cn/static/public/jgo-23-985-3.xls), 5,696 genes obtained by coexpression analysis were significantly related to the four hub genes. All these genes were aggregated into nine modules (Figure 6A). PHF12, FXYD3, PRKCB and ZNF582 belonged to the green, yellow, brown and turquoise modules, respectively. The module-enriched KEGG pathway analysis indicated that green module genes were mainly related to herpes simplex virus 1 infection. The yellow module genes were mainly related to cytokine-cytokine receptor interactions and chemokine signaling pathways. The brown module and turquoise module genes were both mainly related to metabolic pathways (Figure 6B-6E).

Figure 6 Analysis of coexpression network and hub genes. (A) Cluster dendrogram of 5,696 coexpressed genes, which are significantly related to 4 hub genes. Nine modules were obtained. Enrichment analysis of the 4 modules, in which PHF12 belonged to the green module (B), FXYD3 to the yellow module (C), PRKCB to the brown module (D), and ZNF582 to the turquoise module (E). *, P<0.05; **, P<0.01; ***, P<0.001. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Potential drug analysis by PDAC-related signatures

To predict drugs with potential therapeutic effects for PDAC, the four hub genes were introduced into DGIdb. In the DGIdb, drugs that might interact with the four genes were analyzed. The results revealed that only 18 genes had potential interactions with PRKCB (supplementary file 4 available at https://cdn.amegroups.cn/static/public/jgo-23-985-4.xls). These drugs mainly acted as inhibitors of PRKCB, and only one drug acted as an activator (Figure 7A).

Figure 7 Potential drug analysis by PDAC-related signatures. (A) Potential drugs interacted with PRKCB, among which only ingenol mebutate acted as an activator. (B) Consistency cluster matrix of PDAC samples in TCGA. (C) Clustering analysis divided PDAC samples into C1 and C2 categories using the k-NN algorithm. (D) Overall survival rates of patients in the C1 and C2 groups. (E) Response of 27 PDAC cell lines to selected drugs. Activity area indicates sensitivity of the drugs to a certain cell line. A greater value indicates more sensitivity of the drug. Predicted prob is the probability that PDAC cell lines are predicted to be in the C2 category. PDAC, pancreatic ductal adenocarcinoma; TCGA, The Cancer Genome Atlas; t-SNE, t-distributed stochastic neighbor embedding; k-NN, k-nearest neighbor.

Based on the expression data of PHF12, FXYD3, PRKCB and ZNF582, PDAC samples in the TCGA could be further clustered into the optimal C1 and C2 groups (Figure 7B,7C). The two groups also had significantly different prognoses, in which patients in the C1 group had a better outcome (Figure 7D). Using the k-NN algorithm, 27 PDAC cell lines from the Cancer Cell Line Encyclopedia (CCLE, https://portals.broadinstitute.org/ccle) were all divided into the C2 group based on the four gene expression profiles (supplementary file 5 available at https://cdn.amegroups.cn/static/public/jgo-23-985-5.xls). According to the drug responses of the 27 cell lines, the drugs topotecan, PD-0325901, panobinostat, paclitaxel and tanespimycin (17-AAG) were predicted as the five most effective candidate drugs (Figure 7E).

Validation of protein expression of the model-related genes

To further validate the results, IHC analysis was conducted to determine PHF12, FXYD3 and PRKCB protein expression levels in PDAC. The results showed that PHF12, FXYD3 and PRKCB proteins were all upregulated in tumor tissues compared with adjacent normal pancreatic tissues (P<0.05) (Figure 8A). From the immunostaining, we observed that the positive immunoreactivity of PHF12 was primarily localized in the cell nucleus, and FXYD3 and PRKCB were both primarily localized in the cytoplasm and cytomembrane (Figure 8B). Western blotting was also applied to analyze PHF12, FXYD3 and PRKCB protein expression in fresh PDAC tissues and adjacent nontumor tissues. The same results were obtained with IHC analysis, which showed elevated protein expression of PHF12, FXYD3 and PRKCB in tumor tissues compared to adjacent normal tissues (P<0.05) (Figure 8C,8D; Figure S3).

Figure 8 Measurement of PHF12, FXYD3 and PRKCB protein levels in our cohort. (A) The mean protein expression levels of PHF12, FXYD3 and PRKCB in PDAC were all significantly higher than those in adjacent nontumor tissue by immunohistochemistry. (B) Representative immunohistochemical staining images of PHF12, FXYD3 and PRKCB, among which PHF12 is primarily localized in the cell nucleus, and FXYD3 and PRKCB are both primarily expressed in the cytoplasm and cytomembrane. Original magnification: ×200. (C) Representative Western blotting images show that the proteins PHF12, FXYD3 and PRKCB were all overexpressed in PDAC tissue compared with normal adjacent tissue. (D) Western blotting analysis demonstrates that the mean grayscale values of PHF12, FXYD3 and PRKCB are all higher in fresh-frozen PDAC tissues than in matched adjacent normal tissues. MAD, mean areal density; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; PDAC, pancreatic ductal adenocarcinoma.

Discussion

To provide a new strategy for the diagnosis and treatment of pancreatic cancer, it is necessary to perform in-depth exploration of its pathogenesis. In recent years, studies on the correlation between epigenetics and the occurrence and development of pancreatic cancer have become increasingly common (13,14). DNA methylation is a predominant epigenetic modification that affects gene transcription and expression and participates in many physiological activities of cells, such as the inactivation of the X chromosome, cell senescence, and tumorigenesis. Aberrant hypermethylation of CpG islands in the enhancer and promoter regions of tumor suppressor genes plays an important role in tumorigenesis, tumor development and invasion. Hypermethylation of tumor suppressor gene promoters would result in downregulation or silencing of gene expression, while hypomethylation in the promoter region of oncogenes would cause an increase in gene expression (15,16). DNA methylation often occurs at the early stage of tumorigenesis, which has been confirmed in colorectal (17), breast (18), pancreatic (5) and lung cancers (19).

In this study, we analyzed PDAC methylation data in the TCGA and ICGC databases. In the tumor group, the number of hypermethylated regions (Hyper-M) was greater than the number of hypomethylated (Hypo-M) regions. This showed that hypermethylation promotes the downregulation of genes, especially tumor suppressor genes, which might play an important role in the formation of PDAC. DNA methylation always precedes somatic cell mutation and occurs in the early stage of tumorigenesis. A study has shown that DNA methylation is virtually unregulated in every cancer, as cancer cells exhibit extensive differential methylation compared to normal cells (20). To construct a diagnostic predictive model for PDAC, we merged data from TCGA and ICGC and randomly assigned them to the training set and validation set. The four most important DMRs, cg10547050, cg08823209, cg03306374 and cg09568464, were selected, and a diagnostic model was constructed for PDAC, which could well classify the tumor and adjacent normal samples both in the training set and validation set, with AUC values reaching 0.984 and 0.988, respectively. This meant that the DNA methylation signatures could be used as potential early diagnostic biomarkers for pancreatic cancer.

Recently, a study has disclosed that genes with hypermethylation are always tumor suppressors and are usually silenced or expressed at low levels in pancreatic cancer (21). Furthermore, the most common genes with hypomethylation are oncogenes of which its activity or expression increases in the context of tumor progression (21). In this study, we also selectively analyzed the genes regulated by the model-related DMRs. The expression of FXYD3 in tumor tissues was significantly higher than that in normal tissues (P<0.05), while the expression of cg08823209, which regulates FXYD3, showed the opposite result. Another three genes and their corresponding DMRs all showed a significant negative correlation, which was consistent with the negative regulation mode of methylation on gene expression. In our external validation cohort, the protein expression levels of PHF12, FXYD3 and PRKCB were all significantly upregulated in tumor tissues. If we could screen abnormal genes or methylation sites at the early stage of intraepithelial neoplasia, then we could screen patients with early-stage tumors and give appropriate treatment. ZNF582 was more highly expressed in stage I PDAC than in other stages and had a significant relationship with prognosis, indicating that it might serve as a prognostic marker.

KEGG enrichment analysis of the coexpressed genes showed that they were mainly related to herpes simplex virus 1 infection, cytokine-cytokine receptor interactions, chemokine signaling pathways, and metabolic pathways, which indicated that all these selected genes might be involved in the formation and development of pancreatic cancer. Therefore, we analyzed the functions of the four genes to further comprehend their effects on tumorigenesis and tumor progression. PHF12, plant homeodomain (PHD) zinc finger protein 12, encodes a member of the PHD zinc finger family of proteins involved in the regulation of ribosomal biogenesis and senescence (22). FXYD, domain-containing ion transport regulator 3 (FXYD3), encodes a cell membrane protein that can regulate the function of ion pumps and ion channels and plays a role in tumor progression. FXYD3 has been reported to be highly expressed in several types of cancers, including bladder cancer and breast cancer, and is related to survival and metastasis (23,24). FXYD3 overexpression has also been reported in preclinical PDAC models (25). Protein kinase C beta (PRKCB) is a protein kinase C (PKC) gene family member coding PRKCB protein, and is involved in many different cellular functions, such as B-cell activation, apoptosis induction and endothelial cell proliferation (26). PRKCB1 was considered a suppressor of tumorigenic behavior both in vitro and in an in vivo pancreatic cancer model (27) and its overexpression could also promote gastric cancer cell proliferation and mobility (28). ZNF582 is widely regarded as a tumor suppressor gene. Hypermethylation of ZNF582 has been reported in esophageal squamous cell carcinoma and cervical cancer (29,30).

To search for potential drugs that might have therapeutic effects on pancreatic cancer, the DGIdb database was used to analyze drugs associated with the four genes. Only PRKCB had corresponding drug information in the database, and 28 drugs were found, among which 27 drugs all acted as inhibitors. We further verified whether the potential drugs had been investigated or tested in vitro or in vivo. First, the responses to these drugs were assessed in 27 pancreatic cancer cell lines, all clustered into the C2 category based on the expression of the four genes. The results indicated that the cell lines showed similar responses as patients in the C2 group in the TCGA, and this group of patients had shorter survival times than those in the C1 group.

According to the degree of response of the 27 pancreatic cancer cell lines to the drugs, topotecan, PD-0325901, panobinostat, paclitaxel and 17-AAG had the highest activity, indicating that they might have potential therapeutic roles for pancreatic cancer, especially for the C2 group. Topotecan is an inhibitor of topoisomerase I, which can form a complex with DNA and cause double-stranded DNA damage. Currently, topotecan is regarded as the first-line chemotherapy drug for refractory small cell lung cancer in many countries (31), and it is also used to treat recurrent ovarian cancer (32). However, a phase II trial of topotecan performed by Scher et al. demonstrated that it had limited activity in patients suffering from advanced or metastatic pancreatic cancer with a dosage of 1.5 mg/m²/d for 5 days intravenously and repeated every 21 days (33). Furthermore, although topotecan has demonstrated a wide range of antitumor activity in other preclinical and phase I studies, a phase II study was ineffective for patients with pancreatic carcinoma (34). PD-0325901 (mirdametinib) is a selective, non-ATP-competitive mitogen-activated extracellular signal-regulated kinase 1/2 (MEK1/2) inhibitor. MEK/ERK regulates cell proliferation, survival and differentiation stimulated by extracellular signals (35). PD-0325901 potently suppressed MEK/ERK pathways and displayed strong anti-proliferative, apoptotic, anti-angiogenesis and anti-tumor activity in head and neck squamous cell carcinoma cells and triple-negative breast cancer cells (35,36). Williams et al. explored the activation of PD-0325901 in response to radiation in multiple pancreatic tumor cell lines and xenografts and found that PD-0325901 could result in growth arrest, apoptosis, and radiosensitization (37). A phase 1 study performed by van Geel et al. revealed that PD-0325901 combined with the human epidermal growth factor receptor inhibitor dacomitinib did not show acceptable antitumor activity in patients with KRAS mutation-positive pancreatic cancer (38). These studies indicated that the effect of PD-0325901 in clinical models needs further verification, combining with other antitumor treatments might be the direction of future research (38).

Preclinical and clinical data suggested that panobinostat, a nonselective histone deacetylase inhibitor, had potential inhibitory activity in hepatocellular, pancreatic, colorectal, gastric, gastrointestinal stromal tumors, and myeloma (39,40). Panobinostat in combination with other antitumor drugs, such as gemcitabine, bortezomib, BEZ235, MK-1775, and IMC-RON8, could work more efficiently than monotherapy for pancreatic cancer (41). GnP and FOLFIRINOX, established as first-line treatment, have improved survival time for patients suffering from pancreatic cancer in the last few years. In addition, GnP after FOLFIRINOX are expected to be one of the second-line recommendations for patients with unresectable pancreatic cancer (42,43). A phase II study showed that 17-AAG, a heat shock protein 90 inhibitor, combined with gemcitabine appeared to decrease levels of checkpoint kinase 1 (Chk1), which prevented S-phase checkpoint inhibition, but clinical trials had not revealed similar results to support the clinical significance of such regimens (44,45). Another phase II multicenter study combining gemcitabine and 17-AAG also did not result in increased 6-month and overall survival times, as expected, in an interim analysis of 21 patients with stage IV pancreatic cancer (46).

A limitation of the study is that the results were mainly obtained by bioinformation analysis. The results of preclinical tumor cell experiments need to be applied quite cautiously to clinical treatment. More careful and rigorous clinical trials should be designed to verify the preclinical results.

Conclusions

In summary, we screened four significant DMRs using multiomics profiles of PDAC from TCGA and ICGC datasets. The diagnostic model constructed based on the DMRs showed excellent performance for distinguishing tumor and normal tissues. Furthermore, four genes, PHF12, FXYD3, PRKCB and ZNF582, regulated by the corresponding DMRs were identified. Moreover, drug target analysis revealed that PRKCB might be a potential gene target, as it had 28 potential interactions. Five potential drugs were screened out, and these might serve as promising therapeutics for particular PDAC patients.

Acknowledgments

We thank American Journal Experts (www.aje.com) for editing a draft of this manuscript.

Funding: The study was supported by grants from the Henan Province Medical Science and Technology Project of China (Nos. SBGJ202002109 and LHGJ20190136).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-23-985/rc

Data Sharing Statement: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-23-985/dss

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-23-985/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-23-985/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013) and was approved by ethics committee of The First Affiliated Hospital of Zhengzhou University (No. 2023-KY-0863-002). Informed consent was waived due to the retrospective nature of the study.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Mizrahi JD, Surana R, Valle JW, et al. Pancreatic cancer. Lancet 2020;395:2008-20. [Crossref] [PubMed]
2. Zaharia C, Roalsø M, Søreide K. The prognostic relevance of examined lymph nodes for accurate staging of resected pancreatic adenocarcinoma. Hepatobiliary Surg Nutr 2022;11:632-5. [Crossref] [PubMed]
3. Mahadevan KK, McAndrews KM, LeBleu VS, et al. KRASG12D inhibition reprograms the microenvironment of early and advanced pancreatic cancer to promote FAS-mediated killing by CD8+ T cells. Cancer Cell 2023;41:1606-1620.e8. [Crossref] [PubMed]
4. Collisson EA, Bailey P, Chang DK, et al. Molecular subtypes of pancreatic cancer. Nat Rev Gastroenterol Hepatol 2019;16:207-20. [Crossref] [PubMed]
5. Inoue Y, Ito H, Takahashi Y. Neoadjuvant therapy for resectable pancreatic cancers. Hepatobiliary Surg Nutr 2022;11:583-5. [Crossref] [PubMed]
6. Jones PA, Issa JP, Baylin S. Targeting the cancer epigenome for therapy. Nat Rev Genet 2016;17:630-41. [Crossref] [PubMed]
7. Liu B, Pilarsky C. Analysis of DNA Hypermethylation in Pancreatic Cancer Using Methylation-Specific PCR and Bisulfite Sequencing. Methods Mol Biol 2018;1856:269-82. [Crossref] [PubMed]
8. Scherer SW, Lee C, Birney E, et al. Challenges and standards in integrating surveys of structural variation. Nat Genet 2007;39:S7-15. [Crossref] [PubMed]
9. Franck C, Müller C, Rosania R, et al. Advanced Pancreatic Ductal Adenocarcinoma: Moving Forward. Cancers (Basel) 2020;12:1955. [Crossref] [PubMed]
10. Pushpakom S, Iorio F, Eyers PA, et al. Drug repurposing: progress, challenges and recommendations. Nat Rev Drug Discov 2019;18:41-58. [Crossref] [PubMed]
11. Wagner AH, Coffman AC, Ainscough BJ, et al. DGIdb 2.0: mining clinically relevant drug-gene interactions. Nucleic Acids Res 2016;44:D1036-44. [Crossref] [PubMed]
12. Ge PL, Li SF, Wang WW, et al. Prognostic values of immune scores and immune microenvironment-related genes for hepatocellular carcinoma. Aging (Albany NY) 2020;12:5479-99. [Crossref] [PubMed]
13. King G, Green S, Chiorean EG. Finding a role for cancer vaccines in pancreatic cancer: a model of resilience. Hepatobiliary Surg Nutr 2022;11:115-8. [Crossref] [PubMed]
14. Wang Y, Wang Y, Patel H, et al. Epigenetic modification of m(6)A regulator proteins in cancer. Mol Cancer 2023;22:102. [Crossref] [PubMed]
15. Vidal E, Sayols S, Moran S, et al. A DNA methylation map of human cancer at single base-pair resolution. Oncogene 2017;36:5648-57. [Crossref] [PubMed]
16. Joliat GR, Labgaa I, Sulzer J, et al. International assessment and validation of the prognostic role of lymph node ratio in patients with resected pancreatic head ductal adenocarcinoma. Hepatobiliary Surg Nutr 2022;11:822-33. [Crossref] [PubMed]
17. Schuebel KE, Chen W, Cope L, et al. Comparing the DNA hypermethylome with gene mutations in human colorectal cancer. PLoS Genet 2007;3:1709-23. [Crossref] [PubMed]
18. Jeschke J, Van Neste L, Glöckner SC, et al. Biomarkers for detection and prognosis of breast cancer identified by a functional hypermethylome screen. Epigenetics 2012;7:701-9. [Crossref] [PubMed]
19. Su Y, Fang HB, Jiang F. An epigenetic classifier for early stage lung cancer. Clin Epigenetics 2018;10:68. [Crossref] [PubMed]
20. Saghafinia S, Mina M, Riggi N, et al. Pan-Cancer Landscape of Aberrant DNA Methylation across Human Tumors. Cell Rep 2018;25:1066-1080.e8. [Crossref] [PubMed]
21. Gotoh M, Arai E, Wakai-Ushijima S, et al. Diagnosis and prognostication of ductal adenocarcinomas of the pancreas based on genome-wide DNA methylation profiling by bacterial artificial chromosome array-based methylated CpG island amplification. J Biomed Biotechnol 2011;2011:780836. [Crossref] [PubMed]
22. Graveline R, Marcinkiewicz K, Choi S, et al. The Chromatin-Associated Phf12 Protein Maintains Nucleolar Integrity and Prevents Premature Cellular Senescence. Mol Cell Biol. 2017;37:e00522-16. [Crossref] [PubMed]
23. Ramírez-Backhaus M, Fernández-Serra A, Rubio-Briones J, et al. External validation of FXYD3 and KRT20 as predictive biomarkers for the presence of micrometastasis in muscle invasive bladder cancer lymph nodes. Actas Urol Esp 2015;39:473-81. [PubMed]
24. Liu CC, Teh R, Mozar CA, et al. Silencing overexpression of FXYD3 protein in breast cancer cells amplifies effects of doxorubicin and γ-radiation on Na(+)/K(+)-ATPase and cell survival. Breast Cancer Res Treat 2016;155:203-13. [Crossref] [PubMed]
25. Kayed H, Kleeff J, Kolb A, et al. FXYD3 is overexpressed in pancreatic ductal adenocarcinoma and influences pancreatic cancer cell growth. Int J Cancer 2006;118:43-54. [Crossref] [PubMed]
26. Saba NS, Levy LS. Protein kinase C-beta inhibition induces apoptosis and inhibits cell cycle progression in acquired immunodeficiency syndrome-related non-hodgkin lymphoma cells. J Investig Med 2012;60:29-38. [Crossref] [PubMed]
27. Cirigliano SM, Mauro LV, Grossoni VC, et al. Modulation of pancreatic tumor potential by overexpression of protein kinase C β1. Pancreas 2013;42:1060-9. [Crossref] [PubMed]
28. Chen Z, Ju H, Zhao T, et al. hsa_circ_0092306 Targeting miR-197-3p Promotes Gastric Cancer Development by Regulating PRKCB in MKN-45 Cells. Mol Ther Nucleic Acids 2019;18:617-26. [Crossref] [PubMed]
29. Tang L, Liou YL, Wan ZR, et al. Aberrant DNA methylation of PAX1, SOX1 and ZNF582 genes as potential biomarkers for esophageal squamous cell carcinoma. Biomed Pharmacother 2019;120:109488. [Crossref] [PubMed]
30. Zhang C, Fu S, Wang L, et al. The Value and Clinical Significance of ZNF582 Gene Methylation in the Diagnosis of Cervical Cancer. Onco Targets Ther 2021;14:403-11. [Crossref] [PubMed]
31. DiBonaventura MD, Shah-Manek B, Higginbottom K, et al. Adherence to recommended clinical guidelines in extensive disease small-cell lung cancer across the US, Europe, and Japan. Ther Clin Risk Manag 2019;15:355-66. [Crossref] [PubMed]
32. Herzog TJ. Update on the role of topotecan in the treatment of recurrent ovarian cancer. Oncologist 2002;7:3-10. [Crossref] [PubMed]
33. Scher RM, Kosierowski R, Lusch C, et al. Phase II trial of topotecan in advanced or metastatic adenocarcinoma of the pancreas. Invest New Drugs 1996;13:347-54. [Crossref] [PubMed]
34. O'Reilly S, Donehower RC, Rowinsky EK, et al. A phase II trial of topotecan in patients with previously untreated pancreatic cancer. Anticancer Drugs 1996;7:410-4. [Crossref] [PubMed]
35. You KS, Yi YW, Cho J, et al. Dual Inhibition of AKT and MEK Pathways Potentiates the Anti-Cancer Effect of Gefitinib in Triple-Negative Breast Cancer Cells. Cancers (Basel) 2021;13:1205. [Crossref] [PubMed]
36. Mohan S, Vander Broek R, Shah S, et al. MEK Inhibitor PD-0325901 Overcomes Resistance to PI3K/mTOR Inhibitor PF-5212384 and Potentiates Antitumor Effects in Human Head and Neck Squamous Cell Carcinoma. Clin Cancer Res 2015;21:3946-56. [Crossref] [PubMed]
37. Williams TM, Flecha AR, Keller P, et al. Cotargeting MAPK and PI3K signaling with concurrent radiotherapy as a strategy for the treatment of pancreatic cancer. Mol Cancer Ther 2012;11:1193-202. [Crossref] [PubMed]
38. van Geel RMJM, van Brummelen EMJ, Eskens FALM, et al. Phase 1 study of the pan-HER inhibitor dacomitinib plus the MEK1/2 inhibitor PD-0325901 in patients with KRAS-mutation-positive colorectal, non-small-cell lung and pancreatic cancer. Br J Cancer 2020;122:1166-74. [Crossref] [PubMed]
39. Ma X, Ezzeldin HH, Diasio RB. Histone deacetylase inhibitors: current status and overview of recent clinical trials. Drugs 2009;69:1911-34. [Crossref] [PubMed]
40. Singh A, Patel P. The Safety, Efficacy and Therapeutic Potential of Histone Deacetylase Inhibitors with Special Reference to Panobinostat in Gastrointestinal Tumors: A Review of Preclinical and Clinical Studies. Curr Cancer Drug Targets 2018;18:720-36. [Crossref] [PubMed]
41. Singh A, Patel VK, Jain DK, et al. Panobinostat as Pan-deacetylase Inhibitor for the Treatment of Pancreatic Cancer: Recent Progress and Future Prospects. Oncol Ther 2016;4:73-89. [Crossref] [PubMed]
42. Von Hoff DD, Ervin T, Arena FP, et al. Increased survival in pancreatic cancer with nab-paclitaxel plus gemcitabine. N Engl J Med 2013;369:1691-703. [Crossref] [PubMed]
43. Hayuka K, Okuyama H, Murakami A, et al. Gemcitabine Plus Nab-Paclitaxel as Second-Line Chemotherapy following FOLFIRINOX in Patients with Unresectable Pancreatic Cancer: A Single-Institution, Retrospective Analysis. Chemotherapy 2021;66:58-64. [Crossref] [PubMed]
44. Hubbard J, Erlichman C, Toft DO, et al. Phase I study of 17-allylamino-17 demethoxygeldanamycin, gemcitabine and/or cisplatin in patients with refractory solid tumors. Invest New Drugs 2011;29:473-80. [Crossref] [PubMed]
45. Hendrickson AE, Oberg AL, Glaser G, et al. A phase II study of gemcitabine in combination with tanespimycin in advanced epithelial ovarian and primary peritoneal carcinoma. Gynecol Oncol 2012;124:210-5. [Crossref] [PubMed]
46. Pedersen KS, Kim GP, Foster NR, et al. Phase II trial of gemcitabine and tanespimycin (17AAG) in metastatic pancreatic cancer: a Mayo Clinic Phase II Consortium study. Invest New Drugs 2015;33:963-8. [Crossref] [PubMed]