Screening for genes and subnetworks associated with pancreatic cancer based on the gene expression profile.

The present study aimed to screen for potential genes and subnetworks associated with pancreatic cancer (PC) using the gene expression profile. The expression profile GSE 16515 was downloaded from the Gene Expression Omnibus database, which included 36 PC tissue samples and 16 normal samples. Limma package in R language was used to screen differentially expressed genes (DEGs), which were grouped as up‑ and downregulated genes. Then, PFSNet was applied to perform subnetwork analysis for all the DEGs. Moreover, Gene Ontology (GO) and REACTOME pathway enrichment analysis of up‑ and downregulated genes was performed, followed by protein‑protein interaction (PPI) network construction using Search Tool for the Retrieval of Interacting Genes Search Tool for the Retrieval of Interacting Genes. In total, 1,989 DEGs including 1,461 up‑ and 528 downregulated genes were screened out. Subnetworks including pancreatic cancer in PC tissue samples and intercellular adhesion in normal samples were identified, respectively. A total of 8 significant REACTOME pathways for upregulated DEGs, such as hemostasis and cell cycle, mitotic were identified. Moreover, 4 significant REACTOME pathways for downregulated DEGs, including regulation of β‑cell development and transmembrane transport of small molecules were screened out. Additionally, DEGs with high connectivity degrees, such as CCNA2 (cyclin A2) and PBK (PDZ binding kinase), of the module in the protein‑protein interaction network were mainly enriched with cell‑division cycle. CCNA2 and PBK of the module and their relative pathway cell‑division cycle, and two subnetworks (pancreatic cancer and intercellular adhesion subnetworks) may be pivotal for further understanding of the molecular mechanism of PC.

Introduction

Pancreatic cancer (PC), a digestive system tumor, is one of the most aggressive types of cancer. With a high degree of malignancy, rapid progression and poor prognosis, the 1- and 5-year survival rates of patients with PC are only 8 and 3%, respectively (1). American Cancer Society statistics reported that there were an estimated 36,800 related fatalities and 43,140 new cases of PC in 2010. PC remains the fourth leading cause of cancer-related mortality in the United States, despite advances in detection, chemotherapy and surgery (2). In developing countries, for example, in China, the incidence of PC has also been markedly increasing during the past several decades, and PC has been ranked the sixth leading cause of death from malignant disease (3).

The inability to detect PC in its early treatable stage may be the critical factor contributing to high mortality. PC is characterized by the lack of notable clinical symptoms and patients often present with symptoms, such as back pain, weight loss, and digestive problems (4). As many as 80% of newly diagnosed patients with PC are already in the metastatic stage of the disease, which limits the potential for therapeutic intervention (5). At this stage, several epigenetic as well as genetic changes have taken place and result in the silencing of tumor suppressors and overexpression of oncogenes, ultimately leading to tumor progression (6). In recent years, important advances have been made to understand the molecular biology of PC and genetic analyses have verified that the basis of this malignant disease is heterogenous and complex (7). The occurrence and pathogenesis of PC, however, is not yet completely understood.

Similar to the majority of tumors, the development and growth of PC is a multistep process including initiation, progression, invasion and ultimately metastasis. Each step in this process is considered to be driven by the accumulation of genetic alterations (8). Numerous studies involving PC have been conducted in order to identify cancer-causing genes over the past decade, and as a result several cancer-related genes have been identified (9,10). For instance, DPC4, which encodes SMAD family member 4 (SMAD4), is found to be inactivated in ~50% of all PCs (11). KRAS, an oncogene which is associated with cell survival, proliferation and differentiation, has been identified in >90% of patients with PC, with the majority of these being point mutations at codon 12. In addition, it has been demonstrated that the detection of the KRAS mutation may be useful in identifying patients at high risk for developing PC (12). The identification and characterization of cancer-associated genes have increased the understanding of PC development. However, the survival rate has not improved as much in the past years due to the lack of early diagnosis and effective chemotherapeutic treatments. Therefore, identification of genes associated with the development of PC is required.

To test the hypothesis that the FK506-binding protein 51 (FKBP51) may function as a tumor suppressor, Pei et al (13) performed microarray analysis and submitted the expression profile, including 36 pancreatic cancer tissue samples and 16 normal samples, to the Gene Expression Omnibus database (GEO). This previous study was predominantly focused on the functional mechanism of the single gene FKBP51. Based on the gene expression profile submitted by Pei et al (13) and bioinformatics methods, differentially expressed genes (DEGs) between PC tissue samples and normal samples were determined in the present study. Furthermore, functional annotation of DEGs was conducted, followed by the construction of the protein-protein interaction (PPI) network. This study aimed to increase the understanding of the mechanism underlying PC development.

Materials and methods

Affymetrix microarray data

The expression profiles were accessible at the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo) using the series accession number GSE16515, which was deposited by Pei et al (13). This data set was based on the GPL570 platform of [Affymetrix Human Genome U133 Plus 2.0 Array (HG-U133_Plus_2); Affymetrix, Santa Clara, CA, USA] and updated on Aug 22, 2014. A total of 52 chips were divided into 2 groups: PC tissue samples (T-group, n=36) and normal samples (N-group, n=16).

Data processing and DEG identification

The probe-level data were firstly transformed into gene expression data. Then background corrections and quartile data normalization were conducted using the robust multiarray average (RMA) in the affy package (Fred Hutchinson Cancer Research Center, Seattle, WA, USA) with default parameters (14).

To screen DEGs between the T-group and N-group, the Limma package (Linear Models for Microarray Data) in R language was used (15). The raw P-value was adjusted to the false discovery rate based on the Benjamini-Hochberg approach (16,17) using the Limma package (version 3.22.1; Fred Hutchinson Cancer Research Center). DEGs were identified with the cutoff value of FDR<0.05 and |log (fold change)|>1 (18,19).

Subnetwork analysis

Paired Fuzzy SNet (PFSNet) (20) is a powerful method to identify smaller parts of pathways termed subnetworks. Comparison with previously published methods shows that significant subnetworks (and the genes therein) identified by PFSNet are up to 51% (64%) more consistent across independent datasets of the same disease phenotypes, even for datasets based on different platforms (20). In order to obtain the genes and subnetworks that may be associated with the biological characteristics of a sample, PFSNet was used in this study to analyze the subnetworks of the genes from the T-group as well as the N-group based on pathways from PathwayAPI (21), which integrated Wikipathways (22), Kyoto Encyclopedia of Genes and Genomes (KEGG) (23) and Ingenuity (24). Steps were conducted as follows: i) The pathways were divided into several subnetworks according to the genes with high expression level; ii) the subnetwork of each group was scored as 1 or 2, based on the equation as previously described (20); iii) the difference of scores in each group was evaluated by t-test and the subnetworks with significantly different scores were screened out. Parameters in the PFSNet were set as b= 0.5, t1=0.95, and t2= 0.8.

Gene ontology (GO) function annotation and pathway enrichment analysis

The Database for Annotation Visualization and Integrated Discovery (DAVID) provides a comprehensive set of functional annotation tools to determined the biological meaning of a large list of genes (25). DAVID was used for GO function annotation and REACTOME pathway enrichment analysis of the up- and downregulated genes, respectively. P<0.05 was selected as the cut-off criterion.

PPI network construction

Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org/) is an online database which includes experimental as well as predicted interaction information and comprises >1,100 completely sequenced organisms (26). The protein interactions in the STRING database were shown with a confidence score. To identify the interactive associations between the target genes and other genes, the up- and downregulated genes between the T-group and N-group were inputted into STRING and protein pairs with a confidence score ≥0.7 were considered to be significant. Cytoscape (National Institute of General Medical Sciences of the National Institutes of Health, Bethesda, MD, USA) was performed to visualize the PPI network.

The PPI network was complicated; thus, further analysis was required to expose the enriched functional modules of the PPI network using ClusterONE (Clustering with overlapping neighborhood expansion) in Cytoscape (27). Then DAVID was used to annotate the function of genes in each module.

Results

DEG identification

After data preprocessing, the normalized expression profile data were analyzed using Limma package in R language. With FDR<0.05 and |log (fold change)|>1, 1,989 DEGs including 1,461 up- and 528 downregulated genes, were screened out in the T-group compared with the N-group.

In the N-group, 5 significant subnetworks were identified and were shown to be associated with glutathione metabolism (Fig. 1A), leucine and isoleucine metabolism (Fig. 1B), pancreatic cancer (Fig. 1C), calcium signaling pathway (Fig. 1D) and the mitogen-activated protein kinase pathway (Fig. 1E), respectively. Significant subnetworks of genes from the T-group were in association with galactose metabolism (Fig. 2A), alanine, aspartic acid and glutamic acid metabolism (Fig. 2B), intercellular cell adhesion (Fig. 2C) and contraction of vascular smooth muscle (Fig. 2D and E).

Figure 1

Five significant subnetworks of the N-group, associated with (A) glutathione metabolism (B) leucine and isoleucine metabolism (C) pancreatic cancer (D) calcium signaling pathway and (E) the mitogen-activated protein kinase pathway. Red nodes indicate upregulated genes and blue nodes downregulated genes.

Figure 2

Five significant subnetworks of genes from the T-group, associated with (A) galactose metabolism (B) alanine, (C) aspartic acid and glutamic acid metabolism, (D) intercellular junction and (E) contraction of vascular smooth muscle. Red nodes represent the upregulated genes and the blue nodes represent downregulated genes.

Function and pathway annotation

To gain further insight into the function of the identified DEGs, the online biological classification software DAVID was applied to annotate the DEGs. The upregulated genes were enriched in 14 GO subcategories with the most genes enriched in the cell adhesion pathway (Table I). The downregulated genes were enriched in 9 subcategories with the highest number of genes enriched in the proteolysis pathway (Table I). In addition, 8 significant REACTOME pathways for upregulated genes, such as metabolism of carbohydrates, hemostasis and cell cycle, mitotic were identified (Fig. 3A). Moreover, 4 significant REACTOME pathways for downregulated genes, regulation of β-cell development, transmembrane transport of small molecules, hemostasis and metabolism of amino acids were screened out (Fig. 3B).

Table I

Top gene ontology functional enrichment of up- and downregulated genes.

A, Upregulated
Term	Gene count	P-value
GO:0007155~cell adhesion	73	<0.001
GO:0022610~biological adhesion	73	<0.001
GO:0006955~immune response	66	<0.001
GO:0006952~defense response	59	0.000001
GO:0042127~regulation of cell proliferation	58	0.001423
GO:0007049~cell cycle	57	0.001696
GO:0009611~response to wounding	52	0.000002
GO:0010033~response to organic substance	51	0.006513
GO:0055114~oxidation reduction	47	0.004401
GO:0022402~cell cycle process	46	0.000678
GO:0008219~cell death	46	0.046477
GO:0016265~death	46	0.047681
GO:0008283~cell proliferation	43	0.000016
GO:0022403~cell cycle phase	41	0.000024

Term, Gene Ontology pathway name; gene counts, number of differentially expressed genes enriched in the Gene Ontology function.

Figure 3

Pathway enrichment for differentially expressed genes. (A) Upregulated genes and (B) downregulated genes.

PPI construction and module analysis

The up- and downregulated genes between the T-group and N-group were input into the STRING database to identify the significant interactions with a confidence score of ≥0.7. The PPI network reveals the molecular mechanisms of pancreatic cancer, but it contains too many nodes and interactions to select the useful information. Therefore, the functional modules in the PPI network were mined by ClusterONE. In the current study, the significant module with the lowest P-value for the upregulated DEGs was displayed in Fig. 4A, and 5 DEGs with higher connectivity degrees including cyclin-dependent kinase 1 (CDK1), maternal embryonic leucine zipper kinase (MELK), PDZ-binding kinase (PBK), Cyclin A2 (CCNA2) and nucleolar and spindle associated protein 1 (NUSAP1) were included in this module. GO analysis showed that the DEGs in this module were predominantly associated with cell-division cycle (Table II). While, 5 hub genes with the higher degrees including albumin (ALB), carboxypeptidase A1 (pancreatic) (CPA1), colipase, pancreatic (CLPS), epidermal growth factor (EGF) and complement component 5 (C5) were identified in the significant module with the lowest P-value for the downregulated DEGs with the lowest P-values, which is shown in Fig. 4B. Moreover, the DEGs in this module predominantly participated in biological processes, such as response to wounding, endogenous stimulus and regulation of cell proliferation (Table II).

Figure 4

Significant module of (A) upregulated genes and (B) downregulated genes. Red nodes represent the upregulated DEGs. While blue nodes represent downregulated genes.

Table II

Top Gene Ontology annotation of up- and downregulated genes in the significant module with the lowest P-value of the protein-protein interaction network.

A, Upregulated
Term	Gene count	P-value
GO:0007049~cell cycle	44	<0.001
GO:0022402~cell cycle process	38	<0.001
GO:0022403~cell cycle phase	37	<0.001
GO:0000279~M phase	36	<0.001
GO:0000278~mitotic cell cycle	32	<0.001
GO:0051301~cell division	29	<0.001

Discussion

PC is one of the leading causes of cancer-related mortality worldwide; however, the molecular mechanisms of PC progression remain unclear. With the rapid expansion of knowledge on genomics, emerging evidence suggests that the initiation, progression, invasion and metastasis of PC are generally caused by the differential expression of genes. In the present study, a total of 1,989 DEGs including 1,461 up- and 528 downregulated genes were screened out. In line with the results of the study by Pei et al (13), FKBP5 was identified as one of the downregulated genes in the PC samples. To understand the interaction of these DEGs, a PPI network was constructed and the significant module with the lowest P-value for upregulated genes with the top 5 nodes of CDK1, MELK, PBK, CCNA2 and NUSAP1 and the module with the lowest P-value for downregulated genes with the top 5 nodes of ALB, CPA1, CLPS, EGF and C5 were identified. Among all these proteins, CDK1, ALB, CPA1, CLPS and EGF were verified to be associated with PC (28–31). Moreover, the association of MELK and C5 with PC have been demonstrated in certain studies (32–34). However, according to the present results, CCNA2 and PBK, which have not previously been directly associated with PC, may be pivotal for the initiation and progression of PC. In addition, certain subnetworks may be important in PC via the differential expression of genes involved, such as the subnetwork directly associated with PC and the subnetwork associated with intercellular cell adhesion.

Cyclins are a family of proteins that control the progression of cells through the cell cycle by activating cyclin-dependent kinases (CDK) (35). As a member of the cyclin family, CCNA2 is produced at the onset of DNA synthesis in proliferating somatic cells and is critical in cell cycle progression by regulation of transition from G1 to S phase (36). Genetic variants of CCNA2, which may affect the function of the encoded protein by changing gene expression or by altering the protein structure, are found to significantly increase the risk of cancer development in a tissue-specific manner, such as colon, liver and lung cancer (37). In addition, Gao et al (38) reported that CCNA2 was a biomarker for the prognosis of breast cancer and a promising target for developing novel strategies to prevent or even reverse tamoxifen resistance. In addition, the expression of CCNA2 may aid in monitoring tamoxifen efficacy and directing personalized therapies in patients with breast cancer (38). Few previous studies have focused on the association of CCNA2 and PC, while high throughput bioinformatics analysis in the present study indicates that CCNA2 may be important for the initiation and development of PC. In the present study, CCNA2 was found to be upregulated in PC tissue samples, and functional analysis demonstrates that CCNA2 was predominantly enriched in the cell cycle pathway and participates in biological processes, such as regulation of cell proliferation, regulation of cell cycle, cell cycle checkpoint and mitosis. These findings were concordant with those of previous studies (39,40). Therefore, it was hypothesized that CCNA2 may be important in the pathogenesis of PC via regulation of the cell cycle and mitosis, which may further influence tumor occurrence.

PBK, also known as PDZ-binding kinase, is a mitotic protein kinase and its encoding gene, PBK was found to be upregulated in PC tissue samples. Studying upregulated kinases in cancer may provide important clues as to the mechanism of malignant conversion (41,42). PBK is phosphorylated in vitro by Cdc2-cyclin B at a site in the amino terminus (Thr 9) which is implicated in the binding of α-tubulin, and then localizes to mitotic spindles and spindle poles during metaphase (43). Studies regarding PBK have demonstrated that the expression of PBK is regulated by cell cycle-specific transcription factors, such as E2F and CREB/ATF, and knockdown expression of PBK may lead to cytokinetic dysfunction in breast cancer (44,45). Ayllón et al (46) suggested that PBK is involved in DNA damage sensing and repair via phosphorylating c-H2AX. Nandi et al (47) confirmed that with western immunoblotting and immunoprecipation, and yeast two-hybrid analysis, PBK can directly interact with p53, downregulate its expression and attenuate G2/M checkpoint in fibrosarcoma cells, which was hypothesized to be a plausible explanation for the role of PBK in augmenting tumor cell growth. Similarly, the GO-biological process enrichment in the present study predicted that PBK was predominantly associated with nuclear division, cell division, M phase of mitotic cell cycle, and PBK with higher connectivity degree in the module with the lowest P-value of upregulated genes was enriched in cell division cycle. Based on these results, it was inferred that PBK may influence the occurrence of PC by regulating the mitotic cell cycle and other biological processes.

Intercellular cell adhesion determines the polarity of cells and participates in the maintenance of tissues (48). Several studies have shown that cell-cell adhesiveness is generally reduced in human cancer, which may result in influences as follows: Reduced intercellular adhesiveness allows loss and disruption of cell-cell adhesion, resulting in destruction of histological structure, which is the morphological hallmark of malignant tumors (48). Conversely, reduced intercellular adhesiveness is also indispensable for cancer invasion and metastasis (49). In line with the previous studies, subnetworks associated with intercellular cell adhesion were found to be significant in the T-group, and the majority of the genes in this subnetwork were identified to be upregulated in PC. Accordingly, intercellular cell adhesion may be important in the progression of PC.

The significant subnetwork directly associated with PC consisted of six genes, VEGFA, NFKB1, STAT3, PGF, RAC1 and ARHGEF6. Expression levels of the majority of these genes were identified to be significantly higher in PC samples and this subnetwork is directly involved in PC via the differential expression of genes involved. For instance, STAT3 is confirmed to be vital in anti-pancreatic cancer effects through its contributions to the positive feedback loop between reactive oxygen species and autophagy (50). The concentration of PGF is found to be significantly increased in pancreatic carcinoma compared with tumor-free tissue (51). Moreover, activation of RAC1-dependent superoxide generation leads to PC cell proliferation and inhibition of RAC1 may be a potential therapeutic strategy (52). Hence, as demonstrated, subnetworks directly associated with pancreatic cancer may be crucial in the pathogenesis of PC.

In conclusion, the results of this study may increase the understanding of the mechanism of the occurrence and development of PC. CCNA2 and PBK of the module and their relative pathway cell-division cycle may be pivotal for understanding the molecular mechanism of PC. In addition, two subnetworks (pancreatic cancer subnetwork and intercellular adhesion subnetwork) may be highly associated with PC. However, the whole study was conducted based on bioinformatics methods, and the conclusions have not been verified by corresponding experiments yet. Thus, further experiments are urgently required to confirm the results of this study.