Possible Molecular Markers for the Diagnosis of Pancreatic Ductal Adenocarcinoma.

BACKGROUND We aimed to identify pivotal genes and pathways involved in pancreatic ductal adenocarcinoma (PDAC), and explore possible molecular markers for the early diagnosis of the disease. MATERIAL AND METHODS The array data of GSE74629, including 34 PDAC samples and 16 healthy samples, was downloaded from GEO (Gene Expression Omnibus) database. Then, the DEGs (differentially expressed genes) in PDAC samples were compared with healthy samples using limma (linear models for microarray). Gene functional interaction networks were analyzed with Cytoscape and ReactomeFIViz. PPI networks were constructed with Cytoscape software. In addition, PPI (protein-protein interaction) network clustering modules were analyzed with ClusterONE, and the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses for modules were performed. RESULTS A total of 630 upregulated and 1,002 downregulated DEGs were identified in PDAC samples compared with healthy samples. Some ribosomal protein genes with higher average correlation in module 0 were enriched in the ribosome pathway. NUP107 (nucleoporin 107 kDa) and NUP160 (nucleoporin 160 kDa) were enriched in module 3. HNRNPU (heterogeneous nuclear ribonucleoprotein U) with higher average correlation in module 8 was enriched in the spliceosome pathway. The ribosome pathway and the spliceosome pathway were significantly enriched in cluster 1 and cluster 2, respectively. CONCLUSIONS Ribosomal protein genes Nup170, Nup160, and HNRNPU, and the ribosome pathway as well as the spliceosome pathway may play important roles in PDAC progression. In addition, ribosomal protein genes Nup170, Nup160, and HNRNPU may be used as possible molecular markers for the early diagnosis of the disease.

Background

Pancreatic ductal adenocarcinoma (PDAC), the most common pancreatic neoplasm, is expected to turn into the second most common cause of deaths associated with cancer by 2030 [1,2]. It is estimated that a total of 45,220 patients will be diagnosed with pancreatic cancer and 38,460 will die of this disease in the United States in 2013 [3]. Risk factors for PDAC include age, gender, ethnicity, cigarette smoking [4,5]. In addition, a very poor overall prognosis is one of the most important characteristic of PDAC, and the median survival time after diagnosis is only 3–5 months [6]. Thus, accurate early diagnosis for PDAC is needed and discovery of the molecular mechanisms may provide a new option for the early diagnosis of PDAC.

Recently, a wealth of previous studies has been applied to gain a better understanding of the molecular mechanisms of PDAC. One study showed that SEL1L (Sel-1-like) was downregulated by upregulated has-miR-155 in PDAC patients, which suggested that specific miRNAs played important roles in the pathogenesis of this disease [7]. Steele et al. suggested that microRNA (miRNA) such as miR-21 might serve as a screening tool for PDAC in the future [8]. The miR-10b is upregulated in PDAC and can be used as a diagnostic marker in suspicious pancreatic lesions [9]. Furthermore, the serum MMP7 (matrix metalloproteinase 7) level in PDAC patients correlates with metastatic disease and survival [10]. Oji et al. indicated that WT1 (Wilms Tumor 1) played a significant part in the tumorigenesis of PDAC and provided a new option for the treatment of this disease [11]. In addition, mutation of TP53 (tumor protein 53) contributes to rapid progression of premalignant cells to pancreatic tumor and FoxM1 (forkhead box M1) can be used as prognostic molecular marker for PDAC [12,13]. Besides, activation of the Jak/Stat3 (Janus kinase/signal transducer and activator of transcription 3) pathway is related with adverse outcome of PDAC after resection [14]. The Wnt pathway regulates the metastasis-promoting mucin 4 in PDAC, and then contributes to disease metastases and progression [15]. However, since these genes and pathways are not sufficient to clarify the molecular mechanisms of the PDAC, there is a need to identify more pivotal genes and pathways important in the progression of this disease.

Gene expression technologies are increasingly being used to identify DEGs (differentially expressed genes) in neoplastic tissues compared with their normal tissues [16]. In the present study, we downloaded the array data of GSE74629 and analyzed the DEGs associated with PDAC. Then, functional network analyses were performed for these DEGs. In addition, protein-protein interaction (PPI) network and important clustering modules were analyzed. We aimed to find pivotal genes and pathways involved in PDAC, and explore possible molecular markers for the early diagnosis of the disease.

Material and Methods

Microarray data

The array data of GSE74629 was downloaded from GEO (Gene Expression Omnibus, ) database with the platform of GPL10558 (Illumina HumanHT-12 V4.0 expression beadchip, Jaen, Spain). A total of 34 PDAC samples and 16 healthy samples were included in this study.

Data preprocessing

The downloaded non-normalized gene expression data were preprocessed by limma (linear models for microarray) [17] in R package, including background correcting, normalizing and calculating expression. Then, the probe ID was transformed to gene symbol with illuminaHumanv4.db [18] in R annotation package and the probe that was not matched with gene symbol was eliminated. For one gene symbol mapped by several probes, the mean value of these probes was set as the final gene expression level. Finally, totally 21,318 gene expression values were obtained.

Screening of DEGs

The DEGs in PDAC samples compared with healthy samples were analyzed with limma. The p-values were calculated by t-test [19] in limma package. BH (Benjamini-Hochberg) [20] was used to adjust p-values into FDR (false discovery rate) values. FDR <0.05 and |log2FC (fold change)| ≥0.58 were set as cutoff criterion for DEGs.

Functional network analysis

The ReactomeFIViz app [21] can calculate correlations (as weights for edges in the whole functional interaction network) among genes involved in the same functional interactions, apply Monte Carlo localization graph clustering algorithm to the weighted functional interaction network, and form a sub-network for selected network modules on the basis of module size and average correlation.

In the present study, Cytoscape [22] and ReactomeFIViz app were used to analyze gene functional interaction network. Modules size ≥6 and average correlation ≥0.25 were set as criteria. Subsequently, pathway enrichment analyses were performed for each function module, and biological pathways involved by each module gene were observed with FDR <0.05.

PPI network analysis

STRING (Search Tool for the Retrieval of Interacting Genes) [23] is a database used to provide information of interaction of proteins, and neighborhood, gene fusion, co-occurrence, co-expression experiments, databases and text mining were prediction method of this database. In this study, the input gene sets were gene modules and species were Homo sapiens. PPI score >0.9 was set as the cutoff value and at least one protein interaction nodes was a module gene. PPI networks were constructed with Cytoscape software.

PPI network clustering modules analyses

KEGG (Kyoto Encyclopedia of Genes and Genomes) can be used to reveal functions of genes or other molecules [24]. The subnetwork modules in the PPI network tended to take part in common biological processes. ClusterONE [25] was used to analyze network modules, and the KEGG pathway enrichment analyses for modules were performed. Modules that p-value <3E-4 were set as significant modules.

Results

Data processing and DEGs analysis

As shown in Figure 1, a total of 630 upregulated and 1,002 downregulated genes were differently expressed in PDAC samples compared with healthy samples.

Figure 1

Heat map of differentially expressed genes (DEGs). The green represents higher expression levels; the red represents lower expression levels.

A total of 13 subnetwork modules including 146 nodes and 684 gene pairs were included in the network (Figure 2). Particularly, the average correlation of module 0, module 3, and module 8 was greater than 0.6, which indicated strong correlation (Table 1).

Figure 2

Functional network of differentially expressed genes (DEGs). Different colors represent different modules; diamond-shaped nodes represent downregulated genes; circular nodes represent upregulated genes.

Table 1

Information list of 13 sub network modules.

Module	Nodes in module	Node percentage	Average correlation	Node list
0	25	0.1701	0.6255	CAND1, EEF1B2, EIF3E, RPL13, RPL17, RPL21, RPL22, RPL23, RPL26, RPL31, RPL35A, RPL36A, RPL37, RPL39, RPL7, RPS17, RPS23, RPS3A, RPS6, RPS7, SEC11C, SEC61G, SRP9, SRPRB, UPF3A
1	18	0.1224	0.4769	ZFP3, ZFP82, ZNF10, ZNF146, ZNF195, ZNF260, ZNF266, ZNF286B, ZNF30, ZNF302, ZNF329, ZNF483, ZNF514, ZNF529, ZNF544, ZNF696, ZNF789, ZNF91
2	15	0.102	0.4808	ICT1, MRPL1, MRPL12, MRPL24, MRPL3, MRPL38, MRPL39, MRPL48, MRPL50, MRPS17, MRPS23, MRPS26, MRPS30, MRPS33, MRPS35
3	11	0.0748	0.602	CENPT, GSR, KNTC1, KPNA3, MIS12, NUP107, NUP160, NUP205, NUP54, NUP88, PSIP1
4	11	0.0748	0.4621	AKR1C3, CA4, CARS, EIF2B3, FAIM3, MYC, PAICS, PRMT6, SUB1, TFB2M, TNS3
5	11	0.0748	0.4021	BAZ2A, DDAH2, DNMT1, EID2, FKBP3, HDAC2, PADI4, SAP30, TAF6, TBC1D8, WDTC1
6	9	0.0612	0.5257	CD40LG, DUSP12, DUSP14, DUSP3, EIF2AK3, KMT2C, MAP3K4, MAPK14, MAPKAPK5
7	9	0.0612	0.4102	ACTA2, ITGA2B, ITGB5, MYBPC3, MYL9, SORBS3, TLN1, TPM1, TPM2
8	8	0.0544	0.6486	CPSF2, HNRNPA0, HNRNPH1, HNRNPU, SNRPA1, SRSF2, SRSF7, ZRANB2
9	8	0.0544	0.5985	AKR1B1, CFLAR, CTDSPL, DROSHA, SMAD3, TFE3, TGIF2, ZMYND11
10	8	0.0544	0.4305	AK3, CDK5RAP2, CEP78, OFD1, TUBA4A, TUBB4B, TUBGCP5, VBP1
11	7	0.0476	0.4346	CARD16, DNAJC10, ERLEC1, HSPE1, SLX4, USP10, VCP
12	7	0.0476	0.3656	FES, FGF9, GSN, PELP1, PIK3R1, PIK3R5, PXK

Ribosomal protein genes, such as RPL13 (ribosomal protein L13), RPL17 (ribosomal protein L17), RPL21 (ribosomal protein L21), RPL22 (ribosomal protein L22), RPL23 (ribosomal protein L23), RPL26 (ribosomal protein L26), RPL31 (ribosomal protein L31), RPL35A (ribosomal protein L35A), RPL36A (ribosomal protein L36A), RPL37 (ribosomal protein L37), RPL39 (ribosomal protein L39), RPL7 (ribosomal protein L7), RPS17 (ribosomal protein S17), RPS23 (ribosomal protein S23), RPS3A (ribosomal protein S3A), RPS6 (ribosomal protein S6) and RPS7 (ribosomal protein S7) were enriched in module 0. NUP107 (nucleoporin 107 kDa) and NUP160 (nucleoporin 160 kDa) were enriched in module 3. HNRNPU (heterogeneous nuclear ribonucleoprotein U) was enriched in module 8.

Pathway analysis of functional network module and PPI network analysis

Pathways significantly enriched by each module are shown in Table 2. Ribosomal protein genes were mainly enriched in three pathways including SRP-dependent co-translational protein targeting to membrane, eukaryotic translation termination, and ribosome. NUP107 and NUP160 were mainly enriched in two pathways including hexose transport and metabolism of non-coding RNA. HNRNPU was mainly enriched in the spliceosome pathway.

Table 2

Pathways significantly enriched by each module.

Module	GeneSet	FDR	Nodes
0	SRP-dependent cotranslational protein targeting to membrane(R)	<1.667e-04	RPL17, RPL36A, RPL13, RPL37, RPL39, RPL7, RPS3A, RPL31, RPS23, RPL35A, SEC11C, RPL26, SRPRB, RPS6, RPS7, RPL23, RPS17, RPL22, RPL21, SEC61G, SRP9
0	Eukaryotic Translation Termination(R)	<1.667e-04	RPL17, RPL36A, RPL13, RPL37, RPL39, RPL7, RPS3A, RPL31, RPS23, RPL35A, RPL26, RPS6, RPS7, RPL23, RPS17, RPL22, RPL21
0	Ribosome(K)	<1.667e-04	RPL17, RPL36A, RPL13, RPL37, RPL39, RPL7, RPS3A, RPL31, RPS23, RPL35A, RPL26, RPS6, RPS7, RPL23, RPS17, RPL22, RPL21
2	Mitochondrial translation(R)	<1.000e-03	MRPS35, MRPS17, MRPS26, MRPL1, MRPS33, MRPL50, MRPL3, MRPS23, MRPS30, MRPL24, MRPL12, ICT1, MRPL38, MRPL39, MRPL48
2	Ribosome(K)	<5.000e-04	MRPS17, MRPL1, MRPL3, MRPL24, MRPL12
3	ISG15 antiviral mechanism(R)	<1.000e-03	NUP160, NUP88, NUP205, NUP107, NUP54, KPNA3
3	Hexose transport(R)	<5.000e-04	NUP160, NUP88, NUP205, NUP107, NUP54
3	Metabolism of non-coding RNA(R)	<3.333e-04	NUP160, NUP88, NUP205, NUP107, NUP54
5	NoRC negatively regulates rRNA expression(R)	<1.000e-03	SAP30, HDAC2, DNMT1, BAZ2A
5	Signaling events mediated by HDAC Class I(N)	<5.000e-04	SAP30, HDAC2, FKBP3
5	Hedgehog signaling events mediated by Gli proteins(N)	3.83e-02	SAP30, HDAC2
6	Oxidative stress response(P)	<1.000e-03	DUSP3, MAP3K4, DUSP14, MAPK14, MAPKAPK5, DUSP12
6	MAPK signaling pathway(K)	8.00e-03	DUSP3, MAP3K4, MAPK14, MAPKAPK5
6	p38 MAPK signaling pathway(N)	2.17e-02	MAP3K4, MAPK14
7	Muscle contraction(R)	<1.000e-03	SORBS3, TLN1, ACTA2, MYBPC3, ITGB5, TPM2, TPM1, MYL9
7	Hypertrophic cardiomyopathy (HCM)(K)	<5.000e-04	MYBPC3, ITGB5, TPM2, TPM1, ITGA2B
7	Dilated cardiomyopathy(K)	<3.333e-04	MYBPC3, ITGB5, TPM2, TPM1, ITGA2B
8	Processing of Capped Intron-Containing Pre-mRNA(R)	<1.000e-03	SRSF2, SNRPA1, SRSF7, HNRNPH1, CPSF2, HNRNPA0, HNRNPU
8	Spliceosome(K)	<5.000e-04	SRSF2, SNRPA1, SRSF7, HNRNPU
8	spliceosomal assembly(B)	1.33e-03	SRSF2, SNRPA1
9	Regulation of nuclear SMAD2/3 signaling(N)	<1.000e-03	TFE3, SMAD3, TGIF2
9	Regulation of cytoplasmic and nuclear SMAD2/3 signaling(N)	7.00e-03	CTDSPL, SMAD3
10	Mitotic G2-G2/M phases(R)	<1.000e-03	OFD1, TUBGCP5, CEP78, TUBA4A, CDK5RAP2, TUBB4B
10	Assembly of the primary cilium(R)	<5.000e-04	OFD1, CEP78, TUBA4A, CDK5RAP2, TUBB4B
10	Protein folding(R)	<3.333e-04	VBP1, TUBA4A, TUBB4B
11	Protein processing in endoplasmic reticulum(K)	2.00e-03	VCP, DNAJC10, ERLEC1
11	Hedgehog ligand biogenesis(R)	2.15e-02	VCP, ERLEC1
12	Signaling by SCF-KIT(R)	7.00e-03	FGF9, FES, PIK3R1
12	Nongenotropic Androgen signaling(N)	9.00e-03	PELP1, PIK3R1
12	Osteopontin-mediated events(N)	7.67e-03	GSN, PIK3R1

FDR – false discovery rate.

The PPI network for other DEGs interacted with 13 module genes was shown in Figure 3. A total of 243 nodes and 676 interaction pairs were included in this network.

Figure 3

The PPI network of other differentially expressed genes (DEGs) interacted with 13 module genes. White represents other DEGs; the color of the modules was in line with Figure 2.

We obtained two subnetwork modules with ClusterOne (Figure 4), and the GO (Gene Ontology) and KEGG pathways significantly enriched in two subnetwork modules are shown in Figure 5.

Figure 4

Two significantly enriched subnetwork modules (cluster 1 and cluster 2).

Figure 5

GO and KEGG pathways significantly enriched by two subnetwork modules.

Cluster 1 was enriched by most of the ribosomal protein genes, and cluster 2 was enriched by NUP160 and HNRNPU. The ribosome pathway and the spliceosome pathway were significantly enriched in cluster 1 and cluster 2, respectively.

In this study, modules and clusters had similar meanings. The “modules” were obtained using ReactomeFI, and the “clusters” were obtained by ClusterOne. Sometimes, the results obtained by ClusterOne also could be called “modules”. The main difference between modules and clusters was that the former was function modules and the latter was cluster modules. In the present study, the results obtained by ClusterOne were presented as “clusters”.

Discussion

With the gene expression patterns obtained from the GEO database, a total of 630 upregulated and 1,002 downregulated genes were differently expressed in PDAC samples compared with healthy samples in our present study. Our results showed some ribosomal protein genes with higher average correlation in module 0 were enriched in the ribosome pathway. NUP107 and NUP160 were enriched in module 3. HNRNPU, with higher average correlation in module 8, was enriched in the spliceosome pathway. The ribosome pathway and the spliceosome pathway were significantly enriched in cluster 1 and cluster 2, respectively.

RPS6 is an mTOR effector and the mTOR pathway is activated in PDAC [26,27]. Phosphorylation of RPS6 attenuates DNA damage and tumor suppression for pancreatic cancer progression [26]. RPS26 regulates p53 activity and p53 plays an important role in the development of PDAC [28,29]. Furthermore, the abundance of RPS8 can determine the susceptibility of pancreatic cancer cells to gemcitabine treatment [30]. In addition, a higher abundance of RPS8 is related to worse survival of PDAC, and RPS8 may be an important prognostic factor [31]. In addition, downregulation of RPL15 is related to tumor development in PDAC [32]. Thus, some ribosomal protein genes, such as RPS6, RPS26 and RPS8, have been associated with the progression of PDAC. In the present study, some ribosomal protein genes with higher average correlation in module 0 were also enriched in cluster 1. Although the significant roles of some other ribosomal protein genes have not been fully discussed, we speculated that ribosomal protein genes may be involved in the development of PDAC.

Nup170 and Nup160 are significant components of the Nup107–160 complex, and the Nup107–160 complex plays an important role in the regulation of NPC (nuclear pore complex) assembly [33]. On the one hand, the Nup107–160 complex regulates microtubule polymerization at kinetochores [34]. Nup107–160 complex is an essential component for the recruitment of RanGAP1 (Ran GTPase-activating protein 1), RanBP2 (Ran GTPase-activating protein 2), and Crm1 (chromosomal maintenance 1) to these structures at kinetochores [35]. Furthermore, Nup170 is involved in kinetochore function and the overexpression of kinetochore may contribute to the development of tumor via driving chromosome instability [36,37]. On the other hand, Ran is an essential protein in the formation of nuclear envelope, and this mechanism is dependent on importin-β associated with Nup107–Nup160 [38]. Furthermore, Ran plays a significant role in the metastatic progression of cancer [38]. One study indicated that Nup170 was upregulated in pancreatic cancer [33]. In our present study, NUP107 and NUP160 were enriched in module 3, and NUP160 also was enriched in cluster 2. Thus, we speculated that Nup170 and Nup160 may play important roles in the development of PDAC.

HNRNPU can modulate WT1 transcriptional activation by directly interacting with WT1 [39], and WT1 plays a significant part in the tumorigenesis of PDAC [11]. HNRNPU interacts with MDM2 (E3 ubiquitin-protein ligase Mdm2) involving the degradation of p53 [40]. Furthermore, p53 as tumor suppressor genes plays a part in cell cycle regulation, and mutation of this gene may promote the growth of tumor cell [41]. Thus, HNRNPU may play a role in PDAC progression indirectly. In this study, HNRNPU with higher average correlation in module 8 was enriched in cluster 2. Therefore, our results were similar to other research and suggest that HNRNPU may be an important gene involved in the development of PDAC.

One study indicated that Myc might be involved in the early neoplastic progression of PDAC [42]. Furthermore, the spliceosome is a treatment entry point for Myc-driven cancer [43]. Thus, spliceosome may play important roles in the PDAC progression. In this study, the spliceosome pathway was significantly enriched in cluster 2, and HNRNPU was enriched in the spliceosome pathway. Therefore, HNRNPU may play important roles in the development of PDAC via the spliceosome pathway. In addition, the ribosome pathway was significantly enriched in cluster 1, and the ribosomal protein genes mainly enriched in the ribosome pathway. In addition, as aforementioned, the ribosomal protein genes may be involved in the development of PDAC. Therefore, ribosomal protein genes may play significant parts in the development of PDAC via ribosome pathway.

Conclusions

Ribosomal protein genes, Nup170, Nup160, HNRNPU, and the ribosome pathway, as well as the spliceosome pathway, may play important roles in PDAC progression. In addition, ribosomal protein genes, Nup170, Nup160, and HNRNPU may be used as possible molecular markers for the early diagnosis of the disease. A limitation of this study is a lack of experimental verification, thus, further research is needed about genes that are used for early diagnosis of the disease and other possible molecular markers of PDAC.