Identification of DEGs between ICC tissues and adjacent normal tissues
To identify the presence of DEmRNAs, DEmiRs and DElncRNAs in ICC, the SRP126672 dataset of RNA-seq sequencing data, composed of 30 ICC tissues and 27 adjacent normal tissues, was downloaded from SRA database. FastQC and Trimmomatic applications were used to quality-control and filter the original sequencing data, and the DESeq2 package was used to screen DElncRNAs. A total of 912 DElncRNAs were obtained, specifically with 524 upregulated and 388 downregulated DElncRNAs (Fig. 1A, B). The different expression of DElncRNAs could help distinguish ICC tissues and adjacent normal tissues based on principal component analysis (Fig. 1C). DEmiR and DEmRNA sequencing data (counts) of ICC and adjacent normal tissues were obtained from the TCGA database, followed by differential analysis using the edgeR package. Finally, 66 DEmiRs (38 upregulated and 28 downregulated) (Fig. 1D, E) and 5,522 DEmRNAs (3158 upregulated and 2364 downregulated) (Fig. 1F, G) were acquired.
The lncRNA-miR-mRNA ceRNA network in ICC
To better understand the role of lncRNAs and miRs in the ceRNA network of ICC tissues, a lncRNA-miR-mRNA-ceRNA network was established. Initially, 66 DElncRNAs targeted by 66 DEmiRs was retrieved from miRcode database to obtain 31,452 lncRNA-miR pairs. We then screened the significantly downregulated lncRNAs targeted by upregulated miRs and obviously upregulated lncRNAs targeted by downregulated miRs in ICC. Then, only DEGs with the correct trends and targeting relationships served as candidate genes. The selected candidate genes contain 66 pairs of lncRNA-miR, including 18 markedly up-regulated DElncRNAs, 19 down-regulated DElncRNAs, 10 significantly up-regulated miRs, and 2 significantly down-regulated miRs.
In addition, 16457, 7052, and 2081 target mRNAs were predicted from the following databases: TargetScan, miRDB, and miRTarBase respectively. The 10 downregulated candidate DEmiRs were intersected with 2364 downregulated DEmRNAs, yielding 43 shared genes (Fig. 2A), and 2 downregulated candidate DEmiRs were intersected with 3158 upregulated DEmRNAs, yielding 17 share genes (Fig. 2B). The shared genes served as candidate genes and were applied to generate 136 miR-mRNA pairs including 43 significantly down-regulated and 17 significantly up-regulated DEmRNAs. A total of 37 DElncRNAs, 12 DEmiRs, and 60 DEmRNAs were used to construct the ceRNA network. Based on the aforementioned lncRNA-miR pairs and miR-mRNA pairs, a ceRNA network consisting of 37 lncRNA nodes, 12 miR nodes, and 60 mRNA nodes in ICC was constructed (Fig. 2C).
GO and KEGG pathway enrichment analysis on DEmRNAs
To identify the biological functions and pathways of the 60 DEmRNAs in the ceRNA network, we used the DAVID database to perform a GO function enrichment analysis and the KOBAS database for KEGG pathway enrichment analysis. GO enrichment analysis indicated that DEmRNAs related to biological processes were mainly enriched in GO terms including the following: insulin receptor signaling pathway, positive regulation of cell proliferation, cell respiration and other items (p < 0.05). DEmRNAs that were related to cellular components had a close resemblance with mitochondrial inner membrane (p < 0.05). DEmRNAs that were related to molecular functions were mainly enriched in GO terms of chemoattractant activity, growth factor activity, and transcriptional coactivator binding (p < 0.05) (Fig. 3A). The KEGG analysis revealed that DEmRNAs were significantly enriched in pathways such as cancer pathways, calcium signaling pathways, cytokin-cytokine receptor interactions, MAPK singling pathway, and proteoglycans in cancer (p < 0.001) (Fig. 3B). Coherently, 60 significantly different DEmRNAs exerted a pivotal role in the occurrence and development of ICC by the mediation of the above biological processes.
Identification of 6 hub genes in PPI network
The PPI network was constructed based on DEmRNAs, comprising of 60 nodes and 41 edges (Fig. 4A). To identify the top six hub genes in the PPI network, the MCC network topology in the cytoHubba plug-in was used to determine the relationship of DEmRNAs. Then, 6 hub genes were obtained, including Fos proto-oncogene, AP-1 transcription factor subunit (FOS), insulin like growth factor 1 receptor (IGF1R), hepatocyte growth factor (HGF), insulin like growth factor 2 (IGF2), forkhead box O1 (FOXO1), and neurotrophin 3 (NTF3) (Fig. 4B). Finally, the lncRNA-miR-hub gene subnetwork was constructed (Fig. 4C), including 86 ceRNA regulatory modules.
Association of hub gene-associated ceRNA network in ICC with OS of patients
To identify the relationships of 6 reliable hub genes (FOS, IGF1R, HGF, IGF2, FOXO1, and NTF3), and associated DElncRNAs and DEmiRNAs with OS in patients with ICC, a Kaplan-Meier curve analysis was performed. The results indicated that only DElncRNA MME-AS1 and DEmiRNA hsa-miR-182 were correlated with OS of patients with ICC (Fig. 5A, B). Moreover, the high expression of MME-AS1 was inversely correlated with the OS rate of ICC patients, whereas the high expression of hsa-miR-182 was positively correlated with the OS rate of ICC patients, and MME-AS1 could target and orchestrate hsa-miR-182 (Fig. 5C). MME-AS1 and miR-182 have not been reported in ICC according to literature, but miR-182 can inhibit cell proliferation of various cancers, such as ovarian cancer and colorectal cancer [13, 14], suggesting that MME-AS1 may affect the prognosis of patients with ICC by targeting hsa-miR-182.
Validation of FOS, HGF, IGF2, FOXO1, NTF3 and IGF1R expression in ICC
To better validate our analysis, the expression of 6 hub genes were evaluated in 33 ICC tissues and 8 adjacent normal tissues. The results (Fig. 6A) demonstrated that FOS, HGF, IGF2, FOXO1, and NTF3 were remarkably downregulated while IGF1R expression was notably elevated in ICC tissues (|logFC| > 2, FDR < 0.01) compared to that of adjacent normal tissues. To verify the accuracy of the analysis results, we also validated the expression of FOS, HGF, IGF2, FOXO1, NTF3, and IGF1R in the GSE45001 ICC dataset, and found that FOS, IGF2, FOXO1, and NTF3 were markedly poorly expressed, IGF1R was highly expressed, and HGF was slightly poorly expressed in ICC (Fig. 6B-F). Thus, HGF was excluded from the following analysis. The above findings showed that the results of our analysis have certain reliability.
GSEA prediction of biological pathways related to hub genes in ICC
To better understand the biological pathways related to FOS, IGF1R, IGF2, FOXO1, and NTF3, we classified 33 ICC samples from the TCGA database into the high expression group and the low expression group. A GSEA was conducted along with the annotated c2.cp.kegg.v7.0.symbols.gmt gene set in the MSigDB as the reference gene set. Based on the obtained results, the ICC tissues in the FOS, IGF1R, and IGF2 high expression groups were significantly enriched in “the TGF-β signaling pathway”, “the hedgehog signaling pathway”, and “the retinol metabolism”, respectively (Fig. 7A-C); the ICC tissues in the FOXO1 and NTF3 high expression groups were notably enriched on the “type 2 diabetes mellitus” (Fig. 7D, E). As previously reported, TGF-β is one of the main signaling pathways promoting the progression of cholangiocarcinoma. The inhibition of the TGF-β signaling pathway can induce the anti-proliferation properties of cholangiocarcinoma cells [15, 16]. Additionally, it was also documented that the inhibition of the hedgehog signaling pathway, a potential therapeutic target for human cholangiocarcinoma, attenuated the progression of carcinogenesis in vitro and subsequently increased the necrosis rate of cholangiocarcinoma [17, 18]. Retinol metabolism has been less studied in cholangiocarcinoma, but there is literature suggesting that retinol metabolism is a key pathway in the development of cholangiocarcinoma [19]. Moreover, diabetes is associated with an increased risk of ICC, showing similar etiologies [20-22]. However, other prospective studies are still required to validate the aforementioned associations. In summary, it was suggested that five key genes in ICC may exert effects through the TGF-β signaling pathway, the hedgehog signaling pathway, retinol metabolism, and type 2 diabetes mellitus.