Whole-Transcriptome Sequencing Identifies Key Differentially Expressed mRNAs, miRNAs, lncRNAs, and circRNAs Associated with CHOL.

To systematically evaluate the whole-transcriptome sequencing data of cholangiocarcinoma (CHOL) to gain more insights into the transcriptomic landscape and molecular mechanism of this cancer, we performed whole-transcriptome sequencing based on the tumorous (C) and their corresponding non-tumorous adjacent to the tumors (CP) from eight CHOL patients. Subsequently, differential expression analysis was performed on the C and CP groups, followed by functional interaction prediction analysis to investigate gene-regulatory circuits in CHOL. In addition, The Cancer Genome Atlas (TCGA) for CHOL data was used to validate the results. In total, 2,895 differentially expressed messenger RNAs (dif-mRNAs), 56 differentially expressed microRNAs (dif-miRNAs), 151 differentially expressed long non-coding RNAs (dif-lncRNAs), and 110 differentially expressed circular RNAs (dif-circRNAs) were found in CHOL samples compared with controls. Enrichment analysis on those differentially expressed genes (DEGs) related to miRNA, lncRNA, and circRNA also identified the function of spliceosome. The downregulated hsa-miR-144-3p were significantly enriched in the competing endogenous RNA (ceRNA) complex network, which also included 7 upregulated and 13 downregulated circRNAs, 7 upregulated lncRNAs, and 90 upregulated and 40 downregulated mRNAs. Moreover, most of the DEGs and a few of the miRNAs (such as hsa-miR-144-3p) were successfully validated by TCGA data. The genes involved in RNA splicing and protein degradation processes and miR-144-3p may play fundamental roles in the pathogenesis of CHOL.

Introduction

Cholangiocarcinoma (CHOL) is a tumor arising from biliary tract epithelia with features of cholangiocyte differentiation.,,1, 2, 3 CHOL is commonly categorized by anatomic location: intrahepatic, perihilar, and distal CHOL.4, 5, 6 It is the subordinate most extensive hepatic malignancy after hepatocellular carcinoma, accounting for 15%–20% of primary hepatobiliary malignancies, and the overall incidence of CHOL seems to be increasing worldwide.7, 8, 9, 10 Most patients have advanced-stage disease at presentation, and the majority of patients with CHOL have a dismal prognosis; surgical resection is the preferred treatment option for patients who have early-stage disease, not for those at late or advanced stages because available systemic therapies are of limited effectiveness.11, 12, 13 Therefore, increased understanding of its molecular tumor biology is urgently required.

Recent advances in high-throughput next-generation sequencing (NGS) technologies enable the complete sequencing of entire genomes and allow us to analyze a number of cancer genome profiles.14, 15, 16, 17 Recently, fibroblast growth factor receptor 2 (FGFR2) kinase fusion was found to occur in 13.6% of intrahepatic CHOL by parallel whole-transcriptome sequencing in eight specimens. Competing endogenous RNAs (ceRNAs), natural decoys that compete for a common pool of microRNAs (miRNA, miR), represent a novel layer of gene regulation by systematically functionalizing miRNA response element (MRE)-harboring non-coding RNAs, such as long non-coding RNAs (lncRNAs), pseudogenes, and circular RNAs (circRNAs), and forming complex miRNA-mediated ceRNA networks.19, 20, 21, 22, 23 The perturbation of ceRNA crosstalk will disrupt the balance of cellular processes and functions, leading to development of diseases such as cancer.24, 25, 26 Recently, a study reported that lncRNAs H19 and HULC, stimulated by oxidative stress, regulate cell migration and invasion in CHOL by increasing the expression of IL-6 and CXCR4 in the inflammatory process using ceRNA methods of sponging let-7a/let-7b and miR-372/miR-373, respectively. However, few studies have utilized whole-transcriptome sequencing strategies, which allow accurate examination of global gene expression profile, in delineating the transcriptomic landscape of CHOL.

In our study, whole-transcriptome sequencing was performed based on the tumorous (C) and their corresponding non-tumorous adjacent to the tumors (CP) from eight CHOL patients. Subsequently, differential messenger RNA (mRNA), miRNA, lncRNA, and circRNA expression analysis was performed of the C and CP groups, followed by functional interaction prediction analysis. In addition, The Cancer Genome Atlas (TCGA) CHOL data were used to validate the results. This study may provide novel insights into the molecular basis of initiation and progression of CHOL, and provide potential biomarkers for distinct prognostic and therapeutic implications.

Results

Differential Expression Analysis

According to the screening criteria, a total of 2,895 differentially expressed mRNAs (dif-mRNAs) were obtained, of which 2,290 were upregulated and 605 were downregulated; 56 differentially expressed miRNAs (dif-miRNAs) were identified, of which 30 were upregulated and 26 were downregulated; a total of 151 differentially expressed lncRNAs (dif-lncRNAs) were obtained, including 98 upregulated and 53 downregulated; a total of 110 differentially expressed circRNAs (dif-circRNAs) were found, of which 30 were upregulated and 80 were downregulated. The two-way clustering heatmaps of dif-mRNA, dif-miRNA, dif-lncRNA, and dif-circRNA were shown in Figure 1, from which we found that the cancer samples can be significantly separated from the control samples, indicating that the results of the differential expression analysis were reliable.

Figure 1

Heatmaps of Differentially Expressed Molecules

(A–D) Heatmaps of differentially expressed mRNAs (A), differentially expressed miRNAs (B), differentially expressed lncRNAs (C), and differentially expressed circRNAs (D). Red indicates upregulation, and green indicates downregulation.

Functional Enrichment Analysis of dif-mRNAs

When functional enrichment analysis of dif-mRNAs was performed, upregulated genes were enriched in 36 pathways, and downregulated genes were enriched in 46 pathways. In addition, the upregulated genes were enriched in 279 Gene Ontology Biological Process (GO-BP) terms, whereas the downregulated genes were enriched in 174 GO-BP terms. Figure 2 displayed only the top 20 pathways or GO-BP terms enriched by upregulated or downregulated dif-mRNAs. From the results, we found that the upregulated genes were mainly involved in GO:0000398∼mRNA splicing, via spliceosome, GO:0006364∼rRNA processing, hsa03040:Spliceosome, and hsa03013:RNA transport. Additionally, the downregulated genes were significantly associated with GO:0055114∼oxidation-reduction process, GO:0002576∼platelet degranulation, hsa01100:Metabolic pathways, and hsa04610:Complement and coagulation cascades.

Figure 2

Analysis of Gene Ontology Biological Process (GO-BP) and KEGG Pathway

Top 20 pathways or GO-BP terms enriched by upregulated or downregulated mRNAs. (A) Top 20 GO-BP terms enriched by upregulated genes. (B) Top 20 GO-BP terms enriched by downregulated genes. (C) Top 20 pathways enriched by upregulated genes. (D) Top 20 pathways enriched by downregulated genes.

Protein-Protein Interaction (PPI) Network and Module Extraction

The PPI network based on dif-mRNA consisted of 421 nodes and 934 interaction pairs. The nodes with high topological score can be regarded as a key node of the network. Using the Cytoscape plug-in MCODE (score ≥ 10), four sub-network modules were aggregated and extracted from the PPI network, of which the genes were all upregulated (Figure 3A). ModuleA (score = 20.5) contained 21 nodes and 205 interaction pairs, in which the ribosomal proteins were most included, such as ribosomal protein L23a (RPL23A, degree = 24), ribosomal protein S11 (RPS11, degree = 23), and ribosomal protein L8 (RPL8, degree = 23). ModuleB (score = 12.462) contained 14 nodes and 81 interaction pair, including small nuclear ribonucleoprotein D2 polypeptide (SNRPD2, degree = 20), small nuclear ribonucleoprotein polypeptide E (SNRPE, degree = 19), and small nuclear ribonucleoprotein D1 polypeptide (SNRPD1, degree = 18). ModuleC (score = 12) contained 12 nodes and 66 interaction pairs, in which the genes belonged to proteasome 26S subunit family, such as proteasome 26S subunit, non-ATPase 3 (PSMD3, degree = 13) and proteasome 26S subunit, non-ATPase 13 (PSMD13, degree = 12). ModuleD (score = 12) contained 13 nodes and 7 interaction pairs, including M-phase phosphoprotein 10 (MPHOSPH10, degree = 16) and UTP6, small subunit processome component (UTP6, degree = 14).

Figure 3

Four Modules Extracted from Protein-Protein Interaction (PPI) Network and GO-BP Enrichment Analysis

Four modules extracted from PPI network and GO-BP enrichment analysis of genes in these four modules. (A) Significant clustered modules from the PPI network. Orange prism indicates upregulated genes. Node size corresponded to degree value; larger node has higher degree value. (B) Top 10 GO-BP terms enriched by genes in those four modules.

Moreover, genes in the modules were performed for GO-BP enrichment analysis. According to the significance order, the top 10 terms for each module were selected for display (Figure 3B). The genes in moduleA were significantly involved in GO:0006614∼SRP-dependent cotranslational protein targeting to membrane, GO:0019083∼viral transcription, and GO:0000184∼nuclear-transcribed mRNA catabolic process, nonsense-mediated decay; genes in moduleB were significantly associated with GO:0000398∼mRNA splicing, via spliceosome and GO:0000387∼spliceosomal small nuclear ribonucleoprotein particle (snRNP) assembly. Genes in moduleC were concerned with GO:0006521∼regulation of the cellular amino acid metabolic process and GO:0038061∼NIK/nuclear factor-kappa B (NF-κB) signaling. Genes in moduleD were enriched in GO:0006364∼rRNA processing and GO:0030490∼maturation of small subunit rRNA (SSU-rRNA).

Enrichment Analysis of miRNA, lncRNA, and circRNA-Related Target Genes

Based on the mRNAs involved in the dif-miRNA-dif-mRNA regulatory relationship and dif-lncRNA-dif-mRNA and dif-circRNA-dif-mRNA co-expression relationship, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed, and the results were shown by bubble diagram (Figure 4). The results showed that those genes related to miRNA, lncRNA, and circRNA were significantly enriched in endocytosis and spliceosome.

Figure 4

Bubble Maps of KEGG Pathway Enrichment

(A) Pathways enriched by genes regulated by differentially expressed miRNAs. (B) Pathways enriched by genes regulated by differentially expressed lncRNAs. (C) Pathways enriched by genes regulated by differentially expressed circRNAs.

ceRNA Network Construction

Based on the lncRNA-mRNA co-expression relationship and the regulation relationship of dif-miRNA-dif-mRNA and dif-miRNA-dif-lncRNA, the lncRNA and mRNA that were significantly differentially expressed and regulated by the same miRNA were screened. In total, 91 lncRNA-miRNA-mRNA interactions were finally obtained (Figure 5), including 7 up-lncRNAs, 71 upregulated and 12 downregulated mRNAs, and 1 upregulated and 1 downregulated miRNA.

Figure 5

The lncRNA-miRNA-mRNA Network

Orange prism represents the upregulated genes, green circle indicates downregulated genes, red triangle indicates upregulated miRNAs, green arrow shows the downregulated miRNAs, and white quadrilateral indicates the upregulated lncRNAs.

Based on the circRNA-mRNA co-expression relationship and the regulation relationship of dif-miRNA-dif-mRNA and dif-miRNA-dif-circRNA, dif-circRNA and mRNA regulated by the same miRNAs were screened, resulting in finally 494 interaction relationships of circRNA-miRNA-mRNA. There are 22 up-circRNAs and 42 down-circRNAs, 204 upregulated mRNA and 95 downregulated mRNAs, and 6 upregulated miRNAs and 5 downregulated miRNAs. The circRNA-miRNA-mRNA network was shown in Figure 6.

Figure 6

circRNA-miRNA-mRNA Network

Orange prism represents the upregulated genes, green circle represents downregulated genes, red triangle indicates the upregulated miRNAs, green arrow indicates the downregulated miRNAs, white quadrilateral indicates the upregulated circRNA, and white hexagon indicates downregulated circRNA.

Further, according to the lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks, differentially expressed circRNAs, lncRNAs, and mRNAs that were regulated by the same miRNA were further screened. Finally, 158 interaction pairs were obtained (Figure 7), of which 7 upregulated and 13 downregulated circRNAs, 7 upregulated lncRNAs, 90 upregulated and 40 downregulated mRNAs, and 1 upregulated and 1 downregulated miRNA (hsa-miR-144-3p, downregulated; hsa-miR-135a-5p, upregulated) were included.

Figure 7

The Competing Endogenous RNA (ceRNA) Network

Orange prism indicates the upregulated gene, green circle indicates the downregulated gene, red triangle indicates the upregulated miRNA, green arrow indicates the downregulated miRNA, yellow square indicates the upregulated lncRNA, white square indicates upregulated circRNA, and white hexagon indicates downregulated circRNA.

TCGA CHOL Data Verification

Regarding the TCGA CHOL mRNA and miRNA data, the screening threshold was set as p < 0.01 and a total of 9,446 dif-mRNAs and 178 dif-miRNAs were obtained, which were compared with the dif-mRNA and dif-miRNA obtained by our analysis. The Venn plot was presented in Figure S1. As shown in the Venn map, 2,444/2,895 (84.42%) dif-mRNAs were found differentially expressed in the TCGA data and 18/56 (32.14%) dif-miRNAs were also found differentially expressed in the TCGA data. The difference of the other results may have a certain relationship with sample differences or different threshold selection. However, hsa-miR-144-3p was included in the 18 common dif-miRNAs. Genes with a higher degree in the enriched modules were also included in the 2,444 common dif-mRNAs, including RPS11, RPL8, SNRPD2, SNRPE, SNRPD1, PSMD3, PSMD13, MPHOSPH10, and UTP6.

Discussion

NGS has numerous advantages over traditional sequencing technology and has enabled many researchers to identify tumor driver genes in several types of cancer.28, 29, 30, 31 In this study, by applying whole-transcriptome sequencing, we found 2,895 dif-mRNAs, 56 dif-mRNAs, 151 dif-lncRNAs, and 110 dif-circRNAs in CHOL samples compared with controls. Direct function enrichment analysis on the dif-mRNAs showed that these genes were mainly involved in RNA processing and transport, spliceosome, oxidation-reduction process, and platelet degranulation. PPI network analysis identified several hub genes, including ribosomal proteins, ribonucleoprotein, and Proteasome 26S Subunit family. Enrichment analysis on those genes related to miRNA, lncRNA, and circRNA also identified the function of spliceosome. The downregulated hsa-miR-144-3p and upregulated hsa-miR-135a-5p were significantly enriched in the ceRNA complex network. Moreover, most of the differentially expressed genes (DEGs) and a few miRNAs were successfully validated by TCGA data; the difference between the sequencing data and TCGA data may have a certain relationship with sample differences or different threshold selection.

Spliceosome is a multimegadalton ribonucleoprotein machine composed of five snRNPs and numerous other proteins, and can accurately identify precursor mRNAs (pre-mRNAs) splice sites and then catalyze intron excision.32, 33, 34, 35 It should be noted that core spliceosome machinery has been demonstrated to be overexpressed in cancer and affect autophagy and cell proliferation, presenting an attractive therapeutic target for modulating RNA splicing in malignant solid tumors.36, 37, 38, 39 Recently, a study revealed ribosomal protein coding genes and proteasome genes involved in RNA splicing and the protein degradation processes across multiple cancer samples by RNA sequencing (RNA-seq). Moreover, evidence indicates that proteosome inhibition and accumulation of misfolded or unfolded proteins induces dysfunction of endoplasmic reticulum (ER) and caspase-independent cell death of human CHOL cells. In this study, a few observations were particularly noteworthy; for instance, the functional enrichment analysis of differential mRNA, miRNA, lncRNA, and circRNA demonstrated the involvement of spliceosome and RNA processing and transport. Moreover, the genes, including ribosomal proteins, ribonucleoprotein, and Proteasome 26S Subunit family, such as RPS11, RPL8, SNRPD2, SNRPE, SNRPD1, PSMD3, and PSMD13, were identified to be hub in the PPI network, which were also validated by TCGA data. These results strongly argue that the genes, such as RPS11, RPL8, SNRPD2, SNRPE, SNRPD1, PSMD3, and PSMD13, involved in RNA splicing and protein degradation processes may play fundamental roles in the pathogenesis of CHOL. Future studies are needed to consider them in the future targeted therapy development for this cancer.

Additionally, we found the significant role of downregulated hsa-miR-144-3p in the ceRNA network. Evidence has shown that miR-144-3p serves as tumor-suppressive miRNA in many types of cancer, such as glioblastoma, renal cell carcinoma, and esophageal squamous cell carcinoma.42, 43, 44, 45, 46 Wang et al. reported that lncRNA MALAT1 acted as a ceRNA of miR-144-3p for promoting cell proliferation in osteosarcoma cells. In this study, miR-144-3p was found to be in the center of the ceRNA network. Plausibly, miR-144-3p may play an important role in the pathogenesis of CHOL. Future mechanistic studies are warranted to determine the role of miR-144-3p in CHOL.

In conclusion, this study explored the molecular mechanism of CHOL using whole-transcriptome sequencing and through validation by TCGA data. The genes, such as RPS11, RPL8, SNRPD2, SNRPE, SNRPD1, PSMD3, and PSMD13, involved in RNA splicing and protein degradation processes may play fundamental roles in the pathogenesis of CHOL. miR-144-3p may play an important role in the pathogenesis of CHOL. These findings shed light on the genomic complexity of CHOL and also suggest potential new targets in CHOL characterization. However, this study is limited by small sample size in providing a comprehensive overview of CHOL transcriptome. Further in vivo validation using an animal model and functional characterization are needed to delineate the exact mechanistic details.

Materials and Methods

Patients and Tissue Collection

The tissue samples were from eight CHOL patients (four males and four females; age range: 47–75 years, mean age: 60 years) who underwent surgical resection at the Nantong Tumor Hospital, Jiangsu, between 2011 and 2016 and have no preoperative radiation or chemotherapy. All experimental procedures were approved by the Institutional Animal Care and Use Committee (IACUC) guidelines at Nantong Tumor Hospital (NT-19-215). The patients included in this study signed informed consents. Normal bile duct tissues adjacent to the tumors were used as controls. All of the specimens were examined by pathologists before preservation in liquid nitrogen.

RNA Extraction and Quality Control

Total RNA was extracted from frozen tumor tissues (named as C group) and matched adjacent non-cancerous tissues (named as CP group) using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s protocol. The RNA concentration and quality were evaluated with the NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies). Purified RNA was stored at −80°C.

Library Construction and Sequencing

Sequencing libraries were prepared using Illumina TruSeq stranded total RNA with Ribo-Zero gold kit according to the manufacturer’s protocol (Illumina, CA, USA).48, 49, 50, 51 First, rRNA-depleted RNA was fragmented. The cDNA library was constructed using the TruSeq RNA sample Prep Kit (Illumina, San Diego, CA, USA). The libraries were sequenced on an Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions, and 125-bp paired-end reads were generated. Similarly, Illumina’s TruSeq small RNA library preparation kit was used to prepare the miRNA library from samples. Libraries were pooled and sequenced in a single lane of Illumina HiSeq 2000 sequencing platform.

Screening of dif-mRNAs, dif-miRNAs, dif-lncRNAs, and dif-circRNAs

Differential expression analysis was performed on C versus CP using paired t test and Benjamini-Hochberg method. A p value <0.01 and |log2 fold change (FC)| >2 were considered as the thresholds for screening significantly dif-mRNA, dif-miRNA, dif-lncRNA, and dif-circRNA. A heatmap was drawn on the dif-mRNA, dif-miRNA, dif-lncRNA, and dif-circRNA using the pheatmap R package (https://cran.r-project.org/web/packages/pheatmap/index.html).

Analysis of the dif-mRNA Enrichment Pathway

Analysis of the GO-BP functional terms and KEGG pathways in which dif-mRNA was involved was performed using the commonly used enrichment analysis tool DAVID (version 6.8; https://david-d.ncifcrf.gov/). The GO-BP terms and pathways with the enriched gene count ≥2 and the significance threshold p <0.05 were considered significant.

PPI Network and Module Analysis of dif-mRNA

The interaction between DEG-encoded proteins was analyzed by STRING (version 10.0; https://string-db.org/cgi/input.pl) database. The input gene set was set as all DEGs, and the species was set as human. To get the most closely related interaction pairs, we set the parameter PPI score to 0.99 (high confidence). The PPI network was built by Cytoscape software (version 3.2.0; https://cytoscape.org/). The CytoNCA plug-in (version 2.1.6; https://apps.cytoscape.org/apps/cytonca) was used to analyze the network topology properties of the nodes, and the parameter was set to “without weight.” Through the score ranking of each node, the important nodes of the PPI network are obtained. The Cytoscape’s plug-in MCODE (version 1.4.2; https://apps.cytoscape.org/apps/MCODE) was used to analyze the most significant clustering modules in the PPI network. The threshold was selected as score ≥10. The significant clustering module gene was also subjected to GO-BP enrichment analysis. The GO-BP terms with enriched gene count ≥2 and the significance threshold p <0.05 were considered significant.

Prediction of miRNA Regulation Relationship

The miRNA-gene predictive analysis of the dif-miRNA obtained was performed using the miRWalk 2.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) tool. The miRNA-gene regulatory relationship that existed in miRWalk, miRanda, miRDB, miRMap, RNA22, and TargetScan databases was predicted to screen for possible dif-miRNA-dif-mRNA regulatory relationships. The predicted miRNA-gene relationship was integrated with dif-mRNA to obtain dif-miRNA-dif-mRNA regulatory relationships.

The miRNA-lncRNA and miRNA-circRNA regulatory relationships of dif-miRNA were predicted using the StarBase (http://starbase.sysu.edu.cn/) database (v.2.0). The screening criteria were number of supporting experiments low stringency ≥1. The predicted miRNA-lncRNA and miRNA-circRNA regulatory relationships were integrated with dif-lncRNA and dif-circRNA to obtain the dif-miRNA-dif-lncRNA regulation relationship and dif-miRNA-dif-circRNA regulation relationship.

Co-expression Analysis of dif-lncRNA-dif-mRNA and dif-circRNA-dif-mRNA

Using the matrix data of dif-lncRNA, dif-circRNA, and dif-mRNA, the Pearson correlation coefficients of each dif-mRNA and dif-lncRNA and Pearson correlation coefficients of each dif-mRNA and dif-circRNA were calculated; a correlated test was performed to screen significant mRNA-lncRNA pairs and mRNA-circRNA pairs with the thresholds of |r| >0.9 and p <0.05. A list of mRNA-lncRNA and mRNA-circRNA co-expression relationships was further obtained.

KEGG Enrichment Analysis of dif-lncRNA, dif-circRNA, and dif-miRNA

According to the obtained dif-miRNA-dif-mRNA regulatory relationships and the obtained co-expression relationships of dif-mRNA-dif-lncRNA and dif-mRNA-dif-circRNA, the R package clusterprofiler was used for KEGG functional enrichment analysis on the related genes of miRNA, lncRNA, and circRNA in three pairs. The pathway results of count ≥2 and p <0.05 of lncRNA and circRNA were most significantly enriched; the miRNA results were relatively less, and the threshold only limited with p <0.05. Bubble maps are used to show pathways in which miRNA, lncRNA, and circRNA may be involved or affected.

miRNAs are known to cause gene silencing by binding to mRNA largely at the post-transcriptional level, while ceRNAs can act as decoys for miRNA binding and counteract the repressive activity of miRNA, forming a ceRNA network. We integrated the lncRNA-mRNA co-expression relationship and the dif-miRNA-dif-mRNA and dif-miRNA-dif-lncRNA regulatory relationships. According to the mechanism of ceRNA, we were more concerned with the positive correlation expression of lncRNA-mRNA, so the miRNAs capable of simultaneously regulating lncRNA and mRNA were focused on, and the positive correlation coexpression relationship between these miRNA-regulated mRNAs and lncRNA was further obtained. The lncRNA-miRNA-mRNA complex network was constructed by Cytoscape.

Similarly, we integrated the circRNA-mRNA co-expression relationship and the dif-miRNA-dif-mRNA and dif-miRNA-dif-circRNA regulatory relationships. We paid more attention to the positive correlation expression of lncRNA-mRNA, so we focused on miRNAs capable of simultaneously regulating circRNA and mRNA, and further obtained the positive correlation co-expression relationship between mRNA and circRNA regulated by these miRNAs. Cytoscape was used to construct a circRNA-miRNA-mRNA network.

These two complex networks of circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA were integrated, and we focused on screening miRNAs that can regulate circRNA, lncRNA, and mRNA. Positive correlation co-expression relationships of mRNA, lncRNA, and circRNA regulated by these miRNAs were further obtained. Cytoscape was applied to build a ceRNA network.

TCGA CHOL Data Validation

Genomic Data Commons (GDC) TCGA CHOL data provided by the University of California Santa Cruz (UCSC) Xena database (https://xenabrowser.net/datapages/) were used in this study for data validation. The data contained 36 tumor samples and 9 adjacent cancer normal samples. The mRNA and miRNA expression matrices and annotation information were separately downloaded, and the samples were subjected to differential expression analysis in tumor versus normal using unpaired t test in the limma package. A p value <0.01 was considered as a cutoff of significant differential expression in mRNA and miRNA. The differential expression results from TCGA data were compared with the dif-mRNA and dif-miRNA obtained in this analysis, and the Venn map was plotted.

Author Contributions

K.-J.C., Y.-S.M., J.-B.L., D.F., and X.-Q.J. designed the study. K.-J.C., Y.-S.M., X.-H.J., Z.-J.W., T.-M.W., Z.-Z.L., J.-H.W., Q.-X.G., B.Y. J.-B.L., D.F., and X.-Q.J. conducted the study. Y.S., H.-M.W., G.-R.W., J.-B.L., D.F., and X.-Q.J. collected data. K.-J.C., L.-P.G., S.-Q.Z., J.-B.L., D.F., and X.-Q.J. performed the statistical analyses and interpreted the data. K.-J.C., Y.-S.M., J.-B.L., D.F., and X.-Q.J. contributed to study materials and consumables. K.-J.C., Y.-S.M., J.-B.L., D.F., and X.-Q.J. wrote the manuscript. K.-J.C., Y.-S.M., X.-H.J., Z.-J.W., and T.-M.W. contributed equally to this work. All authors contributed to the final version of the manuscript and approved the final manuscript.

Conflicts of Interest

The authors declare no competing interests.