This article is mentioned in:

Introduction

Cholangiocarcinoma (CCA) is the most common biliary malignancy and the second most common hepatic malignancy, following hepatocellular carcinoma (1). CCA accounts for 10–25% of primary liver carcinomas (2). CCA may be classified as intrahepatic (ICC), perihilar or distal CCA depending on the anatomical location (3). CCA is more prevalent in Asia compared with western countries. This is primarily attributed to the increased prevalence of established risk factors, including parasitic infections, bile duct cysts and hepatolithiasis (2). However, the incidence of ICC in the USA has increased from 0.44 to 1.18 cases/100,000 over the past three decades (4). The majority of CCA cases have the characteristics of insidious early and atypical clinical symptoms, rapid progression and poor prognosis. Surgery is the only curative treatment for patients with CCA; however, 50–95% of cases are not surgical candidates (5). The current 5-year survival rate for CCA following surgery and chemotherapy is <20% (6). In addition, CCA is difficult to diagnose, and existing CCA classification systems do not provide insights into the mechanisms of CCA tumorigenesis or potential targets for therapy (7). Therefore, a better understanding of the biology and molecular pathogenesis of CCA may provide the basis to target these markers for tumor prevention or therapy.

Genomic profiling studies have highlighted differing patterns of CCA, helping to stratify patients for targeted therapies (8). Previous studies have investigated the roles of genetic, epigenetic and transcriptomic alterations, in tumor suppressor genes and oncogenes, in the pathogenesis of CCA. Using integrative molecular analysis, Sia et al (9) described two distinct gene signature classes: A proliferation and an inflammatory class. The proliferation class has specific copy number alterations, activation of oncogenic pathways, and is associated with worse outcome. Based on microarray analysis, Jusakul et al (7) additionally revealed four distinct clusters characterized by different clinical features and genomic alterations. These previous results exemplify how genetics, epigenetics and environmental carcinogens may interplay across different geographies to generate distinct molecular subtypes of cancer. In addition, there are studies comparing gene expression profiles between ICC and hepatocellular carcinoma (HCC), in order to identify differences in their carcinogenic mechanisms (10,11). These studies identified genetic alterations in CCA that potentially render early diagnosis and precision treatment a possibility.

However, relatively small sample sizes, and differences in control design and platforms, has led to inconsistencies in terms of the identified genes. Additionally, certain studies have used HCC tissue as a control to identify differentially expressed genes (DEGs) in ICC (10,11). All these aspects increase the heterogeneity of the results. Therefore, an integrated analysis of multiple microarray studies may be helpful to define common DEGs and provide additional evidence for understanding the regulatory mechanism of CCA.

In the present study, an integrated analysis to identify DEGs between CCA and non-tumor tissues was performed by integrating gene expression files in the Gene Expression Omnibus (GEO) database using the web-based tool, NetworkAnalyst. The protein-protein interaction (PPI) network of these genes was subsequently constructed and visualized. In addition, significantly enriched functions of these DEGs were screened and analyzed to identify CCA-associated biological processes and pathways.

Materials and methods

Dataset collection and data processing

Gene expression profiles of CCAs were obtained from the GEO database (www.ncbi.nlm.nih.gov/geo). The following key words were used: ‘Homo sapiens’ and ‘cholangiocarcinoma’. Datasets containing gene expression profiles of CCA and non-tumor tissues or cultured cells were included in the present study. Studies with a sample number <10 were excluded. A total of seven datasets were included in this systematic review. The GEO IDs of these seven datasets were GSE26566, GSE32225, GSE32879, GSE89749, GSE22633, GSE45001 and GSE57555 (Table I) (7,9,11–15).

Table I. Characteristics of the individual studies for integrated analysis.

Table I.

Characteristics of the individual studies for integrated analysis.

			Sample type
Author, year	Datasets	Platforms	CCA, no.	Control, no.	(Refs.)
Jusakul et al, 2017	GSE89749	Illumina HumanHT-12 V4.0 expression beadchip	118	2	(7)
Sia et al, 2013	GSE32225	Illumina HumanRef-8 WG-DASL v3.0	149	6	(9)
Murakami et al, 2015	GSE57555	Agilent-039494 SurePrint G3 Human GE v2 8×60K Microarray	11	11	(11)
Andersen et al, 2012	GSE26566	Illumina humanRef-8 v2.0 expression beadchip	104	6	(12)
Oishi et al, 2012	GSE32879	Affymetrix Human Gene 1.0 ST Array	16	7	(13)
Seol et al, 2011	GSE22633	Illumina human-6 v2.0 expression beadchip	20	4	(14)
Sulpice et al, 2016	GSE45001	Agilent-028004 SurePrint G3 Human GE 8×60K Microarray	10	10	(15)

[i] CCA, cholangiocarcinoma.

Processed data in matrix form were collected from each qualifying microarray study. A global meta-analysis for identifying DEGs in CCA was conducted using the rank product algorithm (RankProd package in R statistical software; www.r-project.org) implemented in the web-based tool NetworkAnalyst (integrative meta-analysis of expression data; www.networkanalyst.ca) (16,17). Normalized gene expression datasets were uploaded into NetworkAnalyst. The datasets were subsequently processed and annotated to adjust the data format and class labels to a consistent style. Following an integrity check, the random effects model was used to calculate the pooled effect size according to the result of Cochran's Q tests (16).

Functional enrichment analysis of DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID; david.abcc.ncifcrf.gov/knowledgebase) (18,19) is a comprehensive set of functional annotation tools. In the present study, gene ontology (GO) enrichment analysis (including biological process, cellular component and molecular function categories) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for DEGs was performed using the DAVID tool. A P-value of 0.05 was selected as the cutoff criterion.

Construction of PPI network

Network analyst uses a comprehensive, high-quality PPI database downloaded from InnateDB (20), which is part of the International Molecular Exchange (IMEx) consortium (21). It additionally contains manually curated protein interaction data from published literature and experimental data from a number of PPI databases, including IntAct (22), MINT (23), DIP (24), BIND (25) and BioGRID (26). To construct the gene coexpression network of the DEGs, the DEGs were mapped on the protein interaction network tool in NetworkAnalyst. As the total nodes were >2,000, the network was switched to zero-order interactions, which composed only of the seed nodes and the edges that interconnect them (16). Visualization and functional analyses were performed. Furthermore, functional modules from the PPI network in NetworkAnalyst were identified and extracted, using enrichment analysis.

Results

Data processing and DEG identification

From the GEO database of the National Center for Biotechnology Information, seven GEO datasets associated with CCA that met our criteria for meta-analysis were extracted (Table I). Among these datasets, a total of 428 CCA cases and 46 controls were included in the integrated analysis.

An overview of the meta-analysis approach is outlined in Fig. 1. The meta-analysis page of the NetworkAnalyst website presents five common approaches for meta-analysis. The present study was performed based on combining effect sizes. According to the result of Cochran's Q test (data not shown), the random effects model was chosen for statistical meta-analysis. DEGs with P<0.05 were selected. A total of 12,081 genes were identified by integrated analysis, and 1,080 DEGs were identified from this meta-analysis, including 710 upregulated and 370 downregulated genes. The 10 most significantly upregulated genes (P<0.05) are secreted phosphoprotein 1 (SPP1), matrix metallopeptidase 11 (MMP11), collagen type I α1 chain (COL1A1), thymosin β10 (TMSB10), agrin (AGRN), collagen type IV α1 chain (COL4A1), Collagen type X α1 chain (COL10A1), minichromosome maintenance complex component 3 (MCM3), collagen type IV α2 chain (COL4A2) and solute carrier family 39 member 1 (SLC39A1; Table IIA). The 10 most significantly downregulated genes (P<0.05) are FXYD domain containing ion transport regulator 1 (FXYD1), cytochrome P450 family 2 subfamily A member 13 (CYP2A13), cystathionine gamma-lyase (CTH), apolipoprotein F (APOF), ornithine carbamoyltransferase (OTC), hydroxyacid oxidase 2 (HAO2), glycine-N-acyltransferase (GLYAT), phosphoenolpyruvate carboxykinase 2 (PCK2), microsomal triglyceride transfer protein (MTTP) and cytochrome P450 family 4 subfamily A member 22 (CYP4A22; Table IIB). The heat map of the top 100 upregulated and downregulated DEGs is presented in Fig. 2.

Figure 1. An overview of the meta-analysis approach. A total of seven published transcriptomic microarray datasets were obtained and uploaded to NetworkAnalyst. Following annotation and normalization, the summary-level data (effect size) were extracted and integrated to identify significant differentially expressed genes. Finally, the selected genes were presented within functional enrichment analysis and PPI network analysis. PPI network visualization, hub and module analysis were further performed using NetworkAnalyst. PPI, protein-protein interaction.

Figure 2. Heat map of differentially expressed genes. Heat maps representing expression patterns for the 100 most upregulated and downregulated genes across different datasets by row-wise comparison. Sample type (the first row): Red indicates the primary tumors of CCA (n=428); grey indicates non-tumor tissue (n=46). For each gene, red is upregulated and blue is downregulated in the corresponding sample. CCA, cholangiocarcinoma.

Table II. Top 10 most significantly up- or downregulated differentially expressed genes.

Table II.

Top 10 most significantly up- or downregulated differentially expressed genes.

A, Upregulated genes
Entrez ID	Name	Combined ES	P-value
6696	SPP1	−2.2097	0.014356
4320	MMP11	−2.0235	0.000592
1277	COL1A1	−1.9385	0.040954
9168	TMSB10	−1.8021	0.014024
375790	AGRN	−1.7029	0.008481
1282	COL4A1	−1.5809	0.004244
1300	COL10A1	−1.5545	0.000309
4172	MCM3	−1.5505	0.001245
1284	COL4A2	−1.5413	0.010102
27173	SLC39A1	−1.5064	0.015808
B, Downregulated genes
Entrez ID	Name	Combined ES	P-value
5348	FXYD1	2.508	0.009848
1553	CYP2A13	2.5026	0.029454
1491	CTH	2.4578	0.000691
319	APOF	2.4369	0.048764
5009	OTC	2.4193	0.039134
51179	HAO2	2.4088	0.040748
10249	GLYAT	2.3862	0.037987
5106	PCK2	2.3849	0.002871
4547	MTTP	2.3174	0.038962
284541	CYP4A22	2.2502	0.013350

[i] DEG, differentially expressed gene; ES, effect size.

Functional enrichment analysis

Functional enrichment analysis was performed to further study these DEGs. Following GO enrichment analysis, three categories (biological process, cellular component and molecular function) were detected using the DAVID database. The 10 most significantly enriched terms (P<0.05) in each category are presented in Fig. 3. KEGG pathway analysis revealed 23 significantly enriched pathways, with ‘DNA replication’ being the most significantly enriched pathway. In addition, ‘influenza A’ and ‘lysosome’ were also significantly enriched pathways (Fig. 4).

Figure 3. GO functional analysis of DEGs. GO enrichment analysis of DEGs was retrieved using DAVID. The 10 most significantly (P<0.05) enriched GO terms in biological process, molecular function and cellular component branches are presented. All the adjusted statistically significant values of the terms were negative 10-base log transformed. DEGs, differentially expressed genes; GO, gene ontology.

Figure 4. Significantly enriched pathways in KEGG pathway analysis of differentially expressed genes in cholangiocarcinoma. Enriched biological pathways were identified using the Database for Annotation, Visualization and Integrated Discovery. Statistically significant (P<0.05) categories among the pathways identified in KEGG. The adjusted statistically significant values were negative 10-base log transformed. KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network construction

Based on the IMEx database, the PPI network of the DEGs was constructed using NetworkAnalyst. The interaction network included 471 nodes and 896 edges (Fig. 5). In the PPI network, degrees were defined to determine the number of neighbors a node is directly connected to, and nodes with higher degrees were considered to be hub proteins. The five most significant hub proteins were histone deacetylase 1 (HDAC1; degree=50), cullin-associated NEDD8-dissociated protein 1 (CAND1; degree=47), ubiquitin D (UBD; degree=44), early growth response protein 1 (EGR1; degree=33) and glycogen synthase kinase 3β (GSK3B; degree=23). KEGG pathways were subsequently extracted from the PPI network. The five most significantly enriched KEGG pathways in the PPI network were ‘influenza A’, ‘cell cycle’, ‘pathways in cancer’, ‘colorectal cancer’ and ‘T cell receptor signaling pathway’ (Fig. 6). Numerous hub genes were associated with these pathways.

Figure 5. An overview of the PPI network. The PPI network was generated using NetworkAnalyst. Red and green color nodes represent upregulated and downregulated differentially expressed genes, respectively. PPI, protein-protein interaction.

Figure 6. The five most significantly enriched Kyoto Encyclopedia of Genes and Genomes pathways extracted from the protein-protein interaction network. (A) Influenza A; (B) cell cycle; (C) pathways in cancer; (D) colorectal cancer; (E) TCR signaling pathway. Red and green colored nodes represent upregulated and downregulated genes, respectively. TCR, T cell receptor.

Discussion

CCA is the most common primary malignancy of the biliary tract. The prognosis of this malignancy is dismal owing to its silent clinical character, difficulties in early diagnosis and limited therapeutic approaches; median survival is less than 24 months (27). There are a number of established risk factors for the development of CCA, although the majority of patients have no identifiable risks (2). Therefore, the identification of novel tumor biomarkers for the early diagnosis and effective treatment of patients with CCA is an important future direction. The transcriptional regulatory network, screened and analyzed with advanced technologies including transcriptomic and proteomic analysis, may be informative to understand the underlying regulatory mechanisms and provide additional evidence for therapeutic applications.

In the present study, a total of 1,080 DEGs were identified based on an integrated analysis. Through a PubMed literature search, it was identified that six out of the 10 most upregulated genes have been associated with CCA in biological or clinical experiments: SPP1 (28,29), MMP11 (30), COL1A1 (31), TMSB10 (32), AGRN (33), and COL4A1 (29). By microarray analysis and reverse transcription-quantitative polymerase chain reaction, Hass et al (28) determined that SPP1 is the most overexpressed gene in ICC. Another study demonstrated that SPP1 expression in the stroma of ICC is significantly associated with the overall patient survival (29). A previous study illustrated the role of MMP11 in cancer progression, with positive MMP11 expression in CCA indicating poor prognosis (30). COL1A1 is a component of type I collagen, which has been reported to be involved in tumor invasion and progression. COL1A1 is significantly upregulated in CCA compared with non-tumor tissues (31). Tissue microarray analysis by Sulpice et al (29) demonstrated increased expression of COL4A1 in the stroma of ICC. Abnormal expression of TMSB10 may contribute to the malignant progression of HCC, and high expression of TMSB10 predicts poor prognosis in patients with HCC following hepatectomy (34). High TMSB10 expression is significantly associated with clinicopathological features, poor prognosis, and distant metastases in patients with breast cancer (35). AGRN is a multidomain heparan sulfate proteoglycan, with different modules homologous to domains present in basement membrane proteins. AGRN expression and deposits are increased in CCA compared with HCC and nontumorous livers, which implies multiple roles in the pathogenesis and progression of CCA (33). However, few downregulated DEGs have been reported to be associated with CCA. Certain genes may be associated with CCA indirectly. For example, the metabolic gene HAO2 is downregulated in HCC, and HAO2 expression levels are inversely correlated with grading, overall survival and metastatic ability (36). Cystathionine-γ-lyase expression is regulated by the Wnt pathway at the transcriptional level and is involved in colon cancer (37); it additionally leads to the development of breast cancer in association with the STAT3 signaling pathway (38). The low expression level of these genes may be associated with the development of CCA, although there are no specific experiments investigating these genes in CCA.

The functional mechanisms of these DEGs using GO and KEGG pathway analyses were further investigated. A total of 482 significantly enriched terms in the biological process, 42 in the cellular component and 73 in the molecular function category were identified. All of the 10 most significantly enriched terms in the biological process category are associated with DNA and chromosomes. Active DNA synthesis means proliferation signaling pathways may be activated in CCA, and an activated cell cycle process may be a sign of proliferation or cancer progression. There were 141 DEGs enriched in the cell cycle process, including enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2), cyclin C and cyclin E1. EZH2 expression promotes the progression of CCA cells by regulating the cell cycle and is associated with poor CCA prognosis (39). In the top 10 cell cycle process terms, apart from terms associated with DNA and chromosomes, there were three GO terms associated with extracellular components: ‘Extracellular exosome’, ‘extracellular vesicle’ and ‘extracellular organelle’. Among the 10 most significantly upregulated DEGs, SPP1, MMP11, COL1A1, COL4A1, COL4A2 and COL10A1 are able to be secreted to the extracellular space (40,41). These results indicated enhanced communication between CCA and the tumor microenvironment. Thus, the tumor microenvironment may be a crucial component governing tumor development and progression and, by targeting it, may be a new direction for CCA treatment (42).

Following KEGG pathway analysis, 22 pathways were identified according to the P-value cut-off. There were five terms, in accordance with the results published by Huang et al (31). Certain pathways, including ‘alcoholism’ (hsa05034), ‘viral carcinogenesis’ (hsa05203) and ‘colorectal cancer’ (hsa05210) have a close association with CCA. Additionally, heavy alcohol use, hepatitis B virus infection, hepatitis C virus infection and inflammatory bowel disease are possible risk factors for CCA (2). Furthermore, a PPI network was established using the NetworkAnalyst visualization tool. PPI network analysis revealed the significant hub proteins, HDAC1, CAND1, UBD, EGR1 and GSK3B. HDACs serve an oncogenic role in the occurrence and development of ICC. Abnormal expression of HDAC1 is significantly associated with lymph node metastasis, high stage carcinoma, vascular invasion and poor prognosis of ICC (43).

It was hypothesized that analyzing multiple datasets may increase the accuracy of the findings, compared with conclusions raised from analyzing a single datasets. In addition, identifying cell markers that are overexpressed in the majority of CCAs may lead to novel drug targets with high specificity. The present results may help to identify combinations of treatments to target various signaling pathways that are altered in CCA. However, certain limitations remain with the present study. Firstly, DEGs were identified via integrated analysis of microarray data. Although differential expression of several of these genes was confirmed by previous biological or clinical research, further in vitro and in vivo validation of these results is required. Secondly, the DEGs identified only represent relative expression levels compared with CCA and non-tumor tissues. Gene mutation and methylation, which are also important for CCA tumorigenesis and progression, were not included in the present analysis. Thirdly, the integrative analysis method used in the present study is based on combining the effect sizes, although there are a number of other algorithms used to identify DEGs, including P-values and gene ranks (15). Different methods may increase the heterogeneity between different published analyses of DEGs. Additionally, the sample cases are much more than the controls (428 vs. 46), and this case-control ratio may be another limitation of the present study.

In summary, the present study identified potential novel markers in the pathogenesis of CCA. In addition, genes consistently differing in expression in CCA were identified through NetworkAnalyst tools. The 10 most significantly upregulated and downregulated genes may serve as potential diagnostic biomarkers. GO annotation and KEGG pathway analysis demonstrated that the identified candidates have a strong association with CCA. Furthermore, a number of novel CCA-associated genes were identified. Further experimental validation is required to fully understand the mechanism of these novel genes in CCA, and to investigate any therapeutic and diagnostic potential these genes may hold.

Acknowledgements

Not applicable.

Funding

The present study was supported by nursery grants from the 175th Hospital of People's Liberation Army, Zhang Zhou, China (grant no. 16Y006).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

The authors' contributions were as follows: WZ, conception and design, data collection, data analysis and manuscript writing; LD, data collection and data analysis; JL, conception and design and data analysis; SZ, conception and design and final approval of manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References