This article is mentioned in:

Introduction

Intrahepatic cholangiocarcinoma (ICC) is a malignant tumor that originates from the epithelial cells of the secondary bile duct and higher order intrahepatic bile duct branches (1). Among primary malignant liver tumors, the incidence of ICC is second to only hepatocellular carcinoma, with increasing incidence yearly (2,3). Because the pathological mechanism and tumor characteristics of ICC differ from those of hepatocellular carcinoma and extrahepatic cholangiocarcinoma, the American Joint Committee on Cancer (AJCC) has formally listed ICC as an independent hepatobiliary system tumor for staging and analysis in the 7th edition of its TNM staging system published in 2010 (4). Although ICC is a rare malignant tumor (accounting for 3% of malignant digestive system tumors and 15% of malignant liver tumors), the global incidence of ICC has increased from 0.32/100,000 to 0.85/100,000 over the past three decades (an increase of 165%). At present, ICC risk factors include primary sclerosing cholangitis, congenital structural abnormalities of the bile duct, inflammatory bowel disease, intrahepatic bile duct stones and parasitic infection. In a recent study, 8% of primary sclerosing cholangitis patients developed cholangiocarcinoma during a 5-year follow-up period. For patients with choledochal cysts or Caroli's disease and other congenital anatomical abnormalities of the biliary system, the risk of malignant transformation in patients 20 years of age or older is 15% (5). In addition, chronic viral hepatitis infection is a high risk factor for ICC. A clinical report on 317 ICC patients by Zhou et al (6) showed that the HBV infection rate of ICC patients is significantly higher than that of the healthy population. However, although many ICC risk factors have been identified, most ICC patients are not exposed to these risk factors, which is an issue for early diagnosis and treatment of ICC.

Currently, the treatment for ICC involves surgical excision as the major measure. For ICC patients who have undergone surgical excision for early to middle stage tumors, the post-operation 3-year survival rate is only 40–50% (7). For patients with middle to late stage tumors who have lost the opportunity for surgery, the average life expectancy without treatment is only 5–8 months (8). Other treatment measures include local ablation, hepatic artery embolization, systemic chemotherapy and targeted therapy. Unfortunately, ICC has a low sensitivity to chemotherapeutic drugs; therefore, there is at present no generalized chemotherapy plan. However, there is preliminary evidence that cisplatin, used alone or in combination with gemcitabine and fluoropyrimidine (S-1), can extend the life expectancy of patients who cannot have surgery. There is currently no specific targeted therapeutic drug for ICC; however, VEGR, EGFR, RAF kinase and multiple receptors involved in the onset of ICC are emerging as targets.

In recent years, modular bioinformatics analysis has been frequently applied in disease and drug studies. This approach not only elucidates the characteristics of gene regulatory networks but also comprehensively and systematically analyzes signaling pathways. Here, we used a modular bioinformatics approach to analyze the expression of ICC genes and determine their roles in ICC progression. We then used a small molecule drug database to identify drugs involved in the regulation of ICC core functional groups. We identified many drugs that regulated the core functional groups that were not previously identified for the treatment of ICC or other cancers. This approach represents a new drug discovery concept for cancer research.

Materials and methods

Screening for differentially expressed mRNAs and miRNAs in cholangiocarcinoma

We analyzed two sets of expression data for cholangiocarcinoma. These data included disease and control samples, with paired mRNAs and miRNAs. Data were acquired from the Gene Expression Omnibus (GEO) database. The series accession numbers for the mRNA and miRNA expression profiles were GSE32879 and GSE32957, respectively (Table I). We organized the data in matrix format to maintain the cholangiocarcinoma samples and normal samples from the two datasets for subsequent analysis. Of these samples, 16 were cholangiocarcinomas and 5 were normal controls. To increase the reliability of target pair prediction and reduce the false-positive rate, we used miRNA-target pairs that were supported by clip-seq data from the Starbase database. To obtain additional potential targets, we selected target pairs predicted by at least one of the five prediction databases (a total of 423,976 target pairs). We pretreated the mRNA and miRNA expression profiles with the k-nearest neighbors algorithm to fill in missing values (define, k=5). We then statistically tested the fold-change values by t-test to screen differentially expressed mRNAs and miRNAs. Fold-change values >1.5 or <1/1.5 (FDR value <0.5) were used as thresholds (Fig. 1).

Figure 1 The fold-change values by t-test were statistically tested to screen differentially expressed mRNAs and miRNAs. (A) Differential expression of mRNAs, (B) differential expression of miRNAs, (C) the screened differential expression of mRNAs, (D) the screened differential expression of miRNAs. Fold-change values >1.5 or <1/1.5 (FDR value <0.5) were used as thresholds.

Table I Gene expression for intrahepatic cholangiocarcinoma from the GEO database.

Table I

Gene expression for intrahepatic cholangiocarcinoma from the GEO database.

Series	Array platform	Samples	Experimental design	Group
GSE32879	GPL6244 Affymetrix Human Gene 1.0 ST Array	GSM13712	Gene expression profiling of 16 ICC and 7 NT tissues were perfomed.	Tumor
		GSm13713		Tumor
		GSM13714		Tumor
		GSM13715		Tumor
		GSM13716		Tumor
		GSM13717		Tumor
		GSM13718		Tumor
		GSM13719		Tumor
		GSM13720		Tumor
		GSM13721		Tumor
		GSM13722		Tumor
		GSM13723		Tumor
		GSM13724		Tumor
		GSM13725		Tumor
		GSM13726		Tumor
		GSM13727		Tumor
		GSM13728		Tumor
		GSM13742		Control
		GSM13743		Control
		GSM13744		Control
		GSM13745		Control
		GSM13746		Control
		GSM13747		Control
		GSM13748		Control
GSE32957	GPL14732 Nanostring nCounter Human microRNA Expression Platform	GSM15966	Gene expression profiling of 16 ICC and 5 NT liver tissues	Tumor
		GSm15967		Tumor
		GSM15968		Tumor
		GSM15969		Tumor
		GSM15970		Tumor
		GSM15971		Tumor
		GSM15972		Tumor
		GSM15973		Tumor
		GSM15974		Tumor
		GSM15975		Tumor
		GSM15976		Tumor
		GSM15977		Tumor
		GSM15978		Tumor
		GSM15979		Tumor
		GSM15980		Tumor
		GSM15981		Tumor
		GSM15996		Control
		GSM15997		Control
		GSM15998		Control
		GSM15999		Control
		GSM16000		Control

[i] ICC, intrahepatic cholangiocarcinoma; NT, non-tumor.

Construction and analysis of the miRNA regulatory networks

The miRNA-mRNA target pairs were required to meet two conditions: opposite differential expression and a predicted target regulation relation. Differential expression was mapped to the network and interacting genes were labeled. Online target gene prediction software was used to identify miRNAs that were differentially expressed relative to the predicted target genes. The predicted target genes were compared with downregulated (or upregulated) mRNAs from the whole-genome microarray data (Fig. 2A). The network conforms to the in-degree distribution (Fig. 2B). Target genes without a signal were eliminated. SPSS 15.0 statistical analysis software was used to analyze the expression values of the miRNA and the corresponding target genes in the normal control group and the tumor group. Pearson's correlation analysis, with P<0.05, was used as a determination standard to screen for target genes that were negatively correlated with miRNA expression. The screening results from each group were integrated into a regulatory network of differentially expressed miRNAs and their ICC target genes. Cytoscape 2.8.2 was used for integrative network analysis of miRNAs and target genes.

Figure 2 Among the two groups of ICC expression profiles, a total of 1,473 pair regulated relations were found to include 488 mRNAs and 37 miRNAs. (A) This differential expression of mRNAs and miRNAs was mapped to the network. (B) We also analyzed the degree distribution of the network nodes; the nodes into degrees obey power-low distribution.

Generation of differentially expressed gene modules and analysis of module function

The differentially expressed genes were mapped to the network. Pairs in which both interacting genes were differentially expressed formed the sub-network. The MCL method was used to identify functional modules (Fig. 3). To functionally analyze the miRNAs and mRNAs that were differentially expressed in ICC, the upregulated and downregulated target genes regulated by differentially expressed miRNAs from each group were run through DAVID Bioinformatics Resources version 6.7 (http://david.Abcc.ncifcrf.gov/) for functional analysis of modules in the network.

Figure 3 As shown in Fig. 2, MCL was used to identify modules within the sub-networks. The modules were defined as modules 1–6.

Drug screening

The connectivity map (CMap) database contains the genome-wide transcriptional expression profile for human cells exposed to active small molecules. The CMap database includes 6,100 groups of small molecule intervention experiments (small molecule intervention group and normal control group), providing a total of 7,056 expression profiles. Of these profiles, 1,309 small molecules are represented. With these 1,309 small molecules from the CMap database, we performed differential gene analysis (drug administration vs. no drug) to obtain a gene set for each drug and its differential expression pattern (that is, how the genes were affected by each small molecule drug). We then compared the genes in each module with the differential genes corresponding to the abovementioned small molecules to determine the correspondence between the genes in each module and the small molecules. A Fisher's test was used to evaluate modular genes corresponding to the small molecules. Small molecules with p-values in the top 10 were considered the most valuable.

Literature review and drug components that connect regulatory network modules

To identify novel antitumor drugs, we analyzed 65 drug components related to the modules. Perl was used to write a program for literature mining. ActivePerl 5.16.2 was then used to mine literature information from Pubmed (NCBI). The mining scope was set to include titles and abstracts, along with drug names, modular gene names and 'cancer' as key words. This search provided articles related to each drug component and cancer. Integrating the drug component screening results, literature mining results and the modular genes, we identified the drugs with the greatest potential antitumor activity for the modular core function groups.

Results

Gene module screening

From comparison of the two sets of ICC and normal tissue data, the fold-change values were used to identify differentially expressed mRNA and miRNAs (detailed information in the Materials and methods). As shown in Table I, we identified 1,925 and 191 upregulated mRNAs and miRNAs, respectively. There were also 1,062 and 28 downregulated mRNAs and miRNAs, respectively. In addition, we performed cluster analysis on differentially expressed genes and miRNAs.

Screening of miRNA and mRNA regulatory networks and ICC modules

Although the mechanism of target gene mRNA degradation or translation inhibition caused by miRNAs is not completely clear, the complimentary pairing of a miRNA seed sequence with the 3′UTR sequence of the target gene is generally considered to be the main condition for a miRNA to act on a target gene. Accordingly, miRNA expression should be opposite to mRNA expression. We selected the upregulated (or downregulated) miRNAs from the miRNA microarray for each group. Targetscan 5.1 was used to predict the target genes from the differentially expressed miRNAs for each group. The target genes were then compared with the downregulated (or upregulated) mRNAs from the mRNA microarray to screen for target genes with differentially expressed miRNAs in each group. Subsequently, the differentially expressed miRNAs from each group and the screening results of their target genes were integrated to establish the regulatory network (Fig. 2A) of differentially expressed miRNAs for ICC and their target genes. The MCL method was used to explore the modules in the network. A total of six modules were identified (termed modules 1–6) (Fig. 3).

Functional analysis of gene modules

We analyzed the six modules in detail using Go_BP, with the arrangement based on the Go score [−log (t-test p-values)]. Each of the p-values of the 10 highest-scored Go analyses was <0.01. These modules were considered the most valuable and were used for subsequent analysis. Go analysis revealed the following: module 1 responded to hypoxia, module 2 responded to cell fraction, module 3 responded to organ development, module 4 responded to glycerophospholipid metabolism, module 5 responded to cytoplasm and module 6 responded to cytoplasm and cellular ketone metabolic process (Fig. 4).

Figure 4 Each module was run through DAVID Bioinformatics Resources version 6.7 for Go_BP analysis of modules in the network. Go term enrichment scores were calculated as the −log(p-value from t-test) and were arranged in descending order.

Drug component regulation of gene modules

Based on the network data, we calculated that there were 1,309 small molecules from the different gene expression experiments; however, the final results showed that a total of 1,108 small molecules induced the expression of the different genes. For the six modules identified by bioinformatic analysis, we compared the genes in each module with the differentially expressed genes that corresponded with the small molecules obtained from the CMap database. We determined the correspondence between the genes in each module and the small molecules. Using a Fisher's exact test on the raw data, we calculated the p-values of each drug component for each module. The p-values were arranged in descending order, with the first 10 drug components with the smallest p-values considered to be the most valuable. These drug components were selected for subsequent analysis. A total of 65 drug components regulated the 6 modules (Fig. 5).

Figure 5 Drugs that regulate the genes in the 6 core modules were identified through bioinformatics techniques. Using the CMap database, we calculated that genes were differentially expressed in 1,309 of the small molecule interference experiments examined. The final results showed that 1,221 small molecules caused these differential expression. Subsequently, the genes from each module were mapped with the differentially expressed genes found from the CMap database, allowing us to determine the relationships between the genes in each module and the tested small molecules. Fisher's exact test was used to calculate the p-value of each drug's regulation of each module, and the p-values were converted into statistical scores [statistical score = −log(Fisher's p-value)].

Discovery of novel potential anticancer drugs

To identify new anticancer drugs, we conducted a literature review and a retrospective study. The drug data were all from recent anticancer therapy studies. Interestingly, anticancer drug components that had not been well-studied showed potential as potent anticancer drugs (Fig. 6). Our review also identified many drug components that have been widely used in anticancer treatment, such as daunorubicin, irinotecan and mercaptopurine. However, many of the drug components we identified have not yet been well-recognized for anticancer therapy, for example, benfotiamine, telenzepine and chlorcyclizine. In addition, pimethixene, etamivan and mafenide have never been studied as anticancer drugs.

Figure 6 Potential new anticancer drugs that regulate the genes in the 6 core modules were identified through literature review. To obtain novel anticancer drugs, we conducted a literature review of the cancer-related research on the 65 drugs. The number of publications was an indicator of the extent of research on the drug's use in cancer treatment. In particular, fewer publications were associated with a higher likelihood that the drug could become a new anticancer drug. Fifteen potential new anticancer drugs were identified (the number of relevant publications was ≤10).

Discussion

The pathogenesis and progression of ICC is a complex, multifaceted process. During this process, various factors interact with each other. While ICC pathogenesis is an important area of study, the discovery of specific therapeutic drug targets for ICC is a more urgent issue that needs to be addressed. An in-depth understanding of drugs that can target and regulate ICC, specifically molecular-targeted therapeutic drugs, is very important for antitumor treatment.

Recently, many ICC tissues and cells from various stages have been collected as samples for whole genome microarray testing. In the present study, we integrated whole genome microarray information from ICC tissues of different origins using a modular bioinformatics approach. Modular analysis is an invaluable bioinformatics technology for functional analysis of single proteins, large-scale protein interactions and gene interactions. It has also been used to study the pathogenesis of many diseases and assess their treatments.

Among the six functional modules identified by modular analysis, module 1 contained the greatest number of genes, with the most complex interaction between the genes. Module 1 was involved in the response to hypoxia, response to drugs, protein tetramerization, response to oxidative stress and glucose metabolism. PTM (post-translational modification) is a key step in protein synthesis (9), involving the phosphorylation, glycosylation and methylation of proteins. One protein may have multiple different post-translational modifications (10). Although PTM has recently gained more appreciation, at present, PTM is rarely investigated in the context of antitumor drugs. The role of PTM has been appreciated in multiple myeloma, breast cancer and lung adenocarcinoma (11–13); however, important PTMs have not been reported for ICC or even CCA. The physiology of tumor cells is very different compared to normal cells, as highlighted by differences in cell growth and differentiation, cell metabolism, biochemical composition and gene regulation. Inhibition of tumor growth is an important facet of anticancer treatment. The change in SSTR1 (somatostatin receptor type 1) in module 2 further confirmed that the modules had the core functions of ICC. Interestingly, we identified significant differences in the ADI1 gene. Because ADI1 has only been previously implicated in prostate cancer, this provides a new direction for future studies (14). Due to its important function in many biological activities, the insulin-like growth factor (IGF) signaling system has attracted wide attention. IGF signaling is associated with the onset of many diseases, such as cretinism, diabetes and cirrhosis of the liver. In module 3, IGFBP5, which is closely associated with the pathogenesis of many cancers, was highly expressed. IGFBP5 suppresses cancer in human melanoma cells (15) and plays important roles in human uterine leiomyoma cells, breast cancer, urothelial carcinoma and papillary thyroid carcinoma (16–18). The role of tumor metabolism has also been appreciated in recent years. The canonical PI3K/AKT/mTOR pathway has been widely studied, as has BNIP3 (19,20). We hypothesize that ICC may inhibit the PI3K/AKT/mTOR pathway by upregulating BnIP3. However, this assertion must be verified by future experiments. The functions of modules 5 and 6 were related to cytoplasm regulation. The survival and prognosis of tumor patients are associated not only with tumor cells but also with the microenvironment of tumor cells. Thus, the pathogenesis of a tumor may be the result of the abnormal microenvironment of tumor cells. In these modules, PTMs, SSTR1, IGFBP5 and BNIP3 were most closely associated with the pathogenesis and progression of the cancer, indicating that these modules are tumor-related gene clusters, with these key genes at their cores. Moreover, in each gene cluster, multiple genes work synergistically to play an even more powerful role in tumor progression. Therefore, in antitumor therapy, the effect of a targeted drug treatment on a gene cluster is very likely better than the effect of targeting individual genes.

There has been a long history of traditional chemotherapy for ICC, with more and more drugs emerging with anticancer effects. However, because ICC itself is not sensitive to chemotherapy, there is no single effective drug, such as sorafenib for liver cancer. To screen for drug components with potential anticancer effects, we used the CMap small molecule database to identify drug components that regulated the core modules of ICC. Among the 65 drug components we identified, some drugs were found to have anticancer effects, as reported in extensive studies. For example, daunorubicin is effective for treating breast cancer (21), while irinotecan is a new drug for treating colon cancer that can strongly inhibit tumor growth and metastasis. Several studies have shown that irinotecan also has significant effects on colorectal liver metastases and metastatic colorectal cancers (22,23). Thus, many of the drug components identified through our screening have known anticancer effects, confirming the accuracy of the screening approach. Apart from the drugs that have been confirmed to have anticancer effects, literature mining identified 12 drug components with limited cancer-related data. For example, pimethixene, etamivan and mafenide have not been reported in any cancer-related studies. In our study, at least one of the drugs that regulated each module was from one of these 12 drug components. Through the Fisher's test, we determined the value of these drug components. For example, mercaptopurine, which regulated module 3, was scored first with daunorubicin, which regulated module 6, listed second. The Fisher's test indicated that the 12 identified drug components not previously reported in the anticancer literature had a strong ICC-regulating effect. This indicates that these drugs may have potent anticancer effects.

In summary, we analyzed two sets of mRNA and miRNA data for ICC. Using modular analysis, we identified six modules that were associated with ICC. Furthermore, based on the CMap database, we found 65 effective small molecule drug components that regulated the six modules. We accurately identified the most meaningful components using the Fisher's test. We then performed a detailed retrospective analysis to select 12 components from the small molecule drug components as new anticancer drug options. Our detailed target gene prediction via differential miRNA analysis and miRNA-mRNA correlation analysis revealed the regulatory roles of miRNAs in the progression of ICC. Thus, we provide a new approach for the selection of anticancer drug. Future experiments will verify the drugs we identified.

Acknowledgments

This study was supported by the Graduate Student Research Innovation fund in Heilongjiang Province of China (grant no. YJSCX2012-201HLJ), Changjiang Scholars and Innovative Research Team in University (grant no. IRT1122), the national natural Scientific Foundation of China (grant nos. 81272705, 81472322 and 81201878), China Postdoctoral Science Foundation Grant (2014M560271) and the natural Science Foundation of Heilongjiang Province (LC201437/H1617). The funders had no role in the study design, data collection and analysis, decision to publish, and preparation of the manuscript.

References