Identification of characteristic genes and construction of regulatory network in gallbladder carcinoma
BMC Medical Genomics volume 16, Article number: 240 (2023)
Gallbladder carcinoma (GBC) is a highly malignant tumor with a poor overall prognosis. This study aimed to identify the characteristic microRNAs (miRNAs) of GBC and the competing endogenous RNA (ceRNA) regulatory mechanisms.
The microarray data of GBC tissue samples and normal gallbladder (NGB) tissue samples from the Gene Expression Omnibus (GEO) database was downloaded. GBC-related differentially expressed miRNAs (DE-miRNAs) were identified by inter-group differential expression analysis and weighted gene co-expression network analysis (WGCNA). Machine learning algorithms were used to screen the characteristic miRNA based on the intersect between least absolute shrinkage and selection operator (LASSO) and Support vector machine-recursive feature elimination (SVM-RFE). Based on the differential expression analysis of GEO database, the ceRNA network of characteristic miRNA was predicted and constructed. The biological functions of the ceRNA network were revealed by carrying out the gene enrichment analysis was implemented. We further screened the key genes of ceRNA network and constructed a protein-protein interaction (PPI) network, and predicted and generated the transcription factors (TFs) network of signature miRNAs. The expression of characteristic miRNA in clinical samples was verified by quantitative real-time polymerase chain reaction (qRT-PCR).
A total of 131 GBC-related DE-miRNAs were obtained. The hsa-miR-4770 was defined as characteristic miRNA for GBC. The ceRNA network containing 211 mRNAs, one miRNA, two lncRNAs, and 48 circRNAs was created. Gene enrichment analysis suggested that the downstream genes were mainly involved in actin filament organization, cell-substrate adhesion, cell-matrix adhesion, reactive oxygen species metabolic process, glutamine metabolic process and extracellular matrix (ECM)-receptor interaction pathway. 10 key genes in the network were found to be most correlated with disease, and involved in cell cycle-related processes, p53, and extrinsic apoptotic signaling pathways. The qRT-PCR result demonstrated that hsa-miR-4770 is down-regulated in GBC, and the expression trend is consistent with the public database.
We identified hsa-miR-4770 as the characteristic miRNA for GBC. The ceRNA network of hsa-miR-4770 may play key roles in GBC. This study provided some basis for potential pathogenesis of GBC.
As one of the common digestive tract tumors, Gallbladder carcinoma (GBC) has a high degree of malignancy, unsatisfactory survival, and prognosis . According to GLOBOCAN 2020, GBC accounts for 0.9% of all cancers diagnosed but 1.7% of all cancer-related deaths . Because of the absence of specific symptoms in the early stages of GBC, the vast majority of patients are diagnosed at advanced stages, and the five-year survival rate is only 10% . More and more studies are coming to explore the treatment modalities for GBC. The surgical resection is only potentially curative therapy for GBC. Unfortunately, most patients with this type of cancer have unresectable disease, only 10–30% of patients can be considered for surgery on presentation . Therefore, further identifying diagnostic biomarkers and potential therapeutic targets for GBC, and achieving early diagnosis and treatment to promote survival rate are the urgent task and the goal we are eager to achieve currently.
Noncoding RNAs (ncRNAs), a class of transcripts without protein-coding potential, produce noncoding transcripts that regulate gene expression and protein function . The ncRNAs are typically classified into two groups according to their size: small ncRNAs of less than 200 nucleotides (nt), which include microRNAs (miRNAs), circular RNAs (circRNAs), Piwi-interacting RNAs and small nucleolar RNAs, and long noncoding RNAs (lncRNAs, >200 nt) . miRNAs are typically 22 nt long and evolutionarily conserved single-stranded RNA molecules, negatively regulating the expression of protein‐coding genes . There have been widely reported that miRNAs played crucial roles in development, metastasis, epithelial-mesenchymal transition, and prognosis of GBC . Yi-Jun Shu et al. reported that overexpression of miR-29c-5p effectively inhibited cancer cell growth, DNA replication and colony formation . Jia-Nan Chen et al. also showed that miR-139-5p significantly inhibited GBC proliferation, invasion, and glucose metabolism by targeting Pyruvate kinase isozyme type M2 (PKM2) . Competing endogenous RNAs (ceRNAs) are transcripts that cross-regulate each other by competing for shared miRNAs . In the ceRNA network, ncRNAs with miRNA response elements (MREs) competitively bind miRNAs to affects miRNA-mRNA binding, thereby regulating the expression of RNAs with the same MREs .
In this study, we employed three transcriptomic datasets associated with GBC as GSE104165, GSE76633 and GSE100363, respectively, and identified the characteristic miRNA of GBC by bioinformatic analysis. Next, based on the prediction of target molecules for key miRNAs, we constructed a potential ceRNA network for GBC. Further, we established a transcription factors (TFs) network of characterized miRNAs. Finally, quantitative real-time polymerase chain reaction (qRT-PCR) was used to detect the differences in the expression of the characteristic miRNA between the GBC samples and control samples. This study provided a basis for understanding the potential ceRNA-related pathogenesis of GBC and identify early diagnostic biomarkers and potential therapeutic targets for GBC.
The work flow of this study is shown in Fig. 1.
Three GBC microarray datasets (GSE104165, GSE76633 and GSE100363) were screened in the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). GSE104165 contains the miRNA expression profiles of eight normal gallbladder (NGB) tissue samples and 40 GBC tissue samples. GSE76633 was comprised of the mRNA and lncRNA expression profiles of nine pairs of GBC tissue samples and paired NGB tissue samples. The circRNA expression profiles of four GBC tissue samples and four matched NGB tissue samples were included in GSE100363. The summary of this dataset is shown in Table 1.
Identification of differentially expressed mRNA (DE-mRNAs), miRNAs (DE-miRNAs), lncRNAs (DE-lncRNAs) and circRNAs (DE-circRNAs)
The DE-mRNAs, DE-miRNAs, DE-lncRNAs, and DE-circRNAs between GBC and NGB samples were identified by “limma” R package . The cut-off criterion was |log2 Fold Change (FC)|> and P-value < 0.05. The corresponding volcano maps and heatmaps were created using the “ggplot2” and “pheatmap” R package respectively.
Weighted gene co-expression network analysis (WGCNA)
The expression data of GSE104165 dataset was analyzed using the “WGCNA” R package  to construct the co-expression network which utilized the tumor and survival as clinical traits. The “good Samples Genes” function was used to perform sample clustering to identify and remove outliers. For making the co-expression network contented the distribution of scale-free network, a soft-thresholding power was computed with the “pick Soft Threshold” function. The dynamic tree cutting method was used to identify different modules with the minimum number of genes in each module was 30. Subsequently, a merging threshold of 0.45 was set to merge similar modules. The correlation between these modules and clinical traits (tumor and survival) was further analyzed. Finally, the module with the highest Pearson correlation coefficient with clinical traits (tumor and survival) was chosen as key module, in which miRNAs with |MM|>0.8 and |GS|>0.2 were identified as key module genes, namely hub genes . Next, the GBC-related DE-miRNAs were identified by overlapping the hub genes and DE-miRNAs, the function annotation for these overlapped miRNAs was conducted using miEAA (https://ccb-compute2.cs.uni-saarland.de/mieaa2/) .
Screening for characteristic miRNAs by machine learning
The “glmnet” R package was utilized to screen for characteristic miRNAs using the least absolute shrinkage and selection operator (LASSO) regression analysis . The area under the receiver operating characteristic curve (ROC) was used to assess the diagnostic sensitivity and specificity of the LASSO model. Support vector machine-recursive feature elimination (SVM-RFE) is a machine learning method in terms of support vector machine algorithms, which can effectively derive a subset of informative genes and make the classification more reliable . The characteristic miRNAs were filtered out by the “e1071” R package with the ten-fold cross-validation method. Eventually, characteristic miRNAs, namely, potential biomarkers, were determined by overlapping the miRNAs identified by LASSO and SVM-RFE. The classification performance of characteristic miRNAs between GBC and NGB samples was evaluated using the area under the ROC curve, which was drawn by “pROC” R package .
Construction of ceRNA nework
We performed target genes (mRNAs), upstream circRNAs, and upstream lncRNAs prediction for the characteristic miRNAs using the Starbasedatabase (http://starbase.sysu.edu.cn/). The predicted mRNAs, circRNAs, lncRNAs, and the corresponding DE-mRNAs, DE-circRNAs, DE-lncRNAs were intersected according to the ceRNA regulatory mechanism. Finally, the lncRNA (circRNA)-miRNA-mRNA network was created by “Cytoscape” software .
Functional enrichment analysis
Gene Ontology (GO) and KEGG enrichment analysis was implemented by “clusterProfiler” R package . The threshold for significance was P-value<0.05. GO enrichment analysis mainly described the biological processes (BP), cellular components (CC), and molecular functions (MF) correlated with genes.
Screening for key genes by varElect
We applied varElect tothe genes in the ceRNA network for disease linkage using the keywords “gallbladder carcinoma” and “cancer”. The varElect collected information based on relevant literature reports and databases for disease linkage detection [21, 22]. The higher the score value, the greater the probability of disease linkage.
Interaction network analysis by GeneMANIA
GeneMANIA (http://www.genemania.org) is a database focused on gene function and related genes, including co-expression, physical interactions, genetic interactions, co-localization, and pathways. We accessed GeneMANIA to search for functions and interactions associated with key genes, and enriched the protein-protein interaction (PPI) network.
Construction of the miRNA-Transcription factors (TFs) network by TransmiR
In order to investigate the regulatory mechanism of characteristic miRNAs, we utilized TransmiR (http://www.cuilab.cn/transmir) to predict the characteristic miRNAs associated TFs and generate the miRNA-TFs network .
Patients and specimens
The study was confirmed by the China Ethics Committee of Registering Clinical Trials (Ethical review number: Chi ECRCT20190065). After obtaining informed consent from all participants, 20 tissue specimens were selected from 10 patients with GBC, including 10 tumor tissues and 10 paired adjacent normal gallbladder tissues. The selected GBC patients underwent surgery at the Second Affiliated Hospital of Kunming Medical University from June 2020 to June 2021, and the pathological diagnosis was finally diagnosed as GBC.
RNA extraction and quantitative real-time polymerase chain reaction (qRT-PCR)
Total RNA from fore-mentioned 10 NGB tissue samples and 10 GBC tissue samples was isolated using the TRIzol Reagent following the manufacturer’s instructions (Ambion, USA). Next, total RNA was reverse transcribed into cDNA utilizing the SweScript-First-strand-cDNA-synthesis-kit (Servicebio, China), according to the manufacturers’ protocol. qRT-PCR was subsequently performed using the 2×Universal Blue SYBR Green qPCR Master Mix (Servicebio, China). The following thermocycling conditions were used for qPCR: one cycle at 95 °C for 60 s (initial denaturation), followed by 40 cycles of 20 s at 95 °C (denaturation), 20 s at 55 °C (annealing), and 30 s at 72 °C (extension). The relative expression level was normalized to the endogenous control U6 and calculated using the 2 − ΔΔCq method . The primers of the target genes were: hsa-miR-4770-F: TGAGATGACACTGTAGCT, hsa-miR-4770-R: GTCGTATCCAGTGCAGGGT; U6-F: CTCGCTTCGGCAGCACA, U6-R: AACGCTTCACGAATTTGCGT.
All statistical analyses were performed using R software (version 4.3.1) and its appropriate packages One-way ANOVA analysis was used to compare data between GBC and NGB groups. ROC curve was generated to evaluate the sensitivity and specificity in distinguishing GBC and NGB samples. If not specified above, a P-value less than 0.05 was considered statistically significant.
Identification of GBC-related DE-miRNAs
To authenticated the GBC-related DE-miRNAs, we conducted differential expression analysis firstly. Based on the miRNA sequencing data in the GSE104165 dataset, we detected 360 DE-miRNAs in GBC samples compared to NGB samples using the “limma” R package, including 182 up-regulated miRNAs and 178 down-regulated miRNAs (Fig. 2, Table S1). To further mine the miRNAs associated with GBC, we performed WGCNA using the data in the GSE104165 dataset. No obvious outliers were excluded by cluster analysis (Fig. 3A). 7.0 was selected as the optimal soft threshold with an R2 = 0.85 (Fig. 3B). Based on the optimal soft threshold, we divided the miRNAs into different modules according to dynamic tree cutting algorithm. After merging, a total of six modules were generated (Fig. 3C and D). Correlations between modules and clinical traits (survival and tumor) were calculated (Fig. 3E). The turquoise module was considered as the key module as the |correlation coefficient| between module and clinical traits (survival and tumor) was the highest (Fig. 3E). According to the criterion of |MM|>0.8 and |GS|>0.2, 175 miRNAs related to tumor in turquoise module, and 169 miRNAs related to survival in turquoise module were authenticated as hub miRNAs, namely GBC-related miRNAs (Fig. 3F and G). Hence, 131 GBC-related DE-miRNAs were obtained by taking the intersection of hub miRNAs related to tumor, hub miRNAs related to survival, up-regulated miRNAs, and down-regulated miRNAs (Fig. 4A, Table S2) The results from the miRNA enrichment analysis showed that those dysregulated miRNAs were mainly associated with VEGF signaling pathway and positive regulation of endothelial cell proliferation (Fig. 4B, Fig. S1 and Table S3), which were considered as important for GBC progression.
Characteristic miRNAs for GBC
To determine the characteristic miRNAs from the 131 GBC-related DE-miRNAs, we adopted LASSO and SVM-RFE algorithms in GSE104165 dataset. Six characteristic miRNAs (hsa-miR-154-5p, hsa-miR-204-5p, hsa-miR-320b, hsa-miR-365a-3p, hsa-miR-379-5p, and hsa-miR-4770) were discerned using the LASSO algorithm (Fig. 5A and B). ROC curve with area under curve (AUC) value was 1.0 indicated the LASSO model performed well (Fig. 5C). Meanwhile, eight characteristic miRNAs (hsa-miR-5701, hsa-miR-125a-5p, hsa-miR-148b-3p, hsa-miR-128, hsa-miR-4770, hsa-miR-654-3p, hsa-miR-34a-5p, and hsa-miR-374a-5p) were selected with the SVM-RFE algorithm (Fig. 5D). Hence, one miRNA, namely hsa-miR-4770, was defined as characteristic miRNA for GBC by overlapping the miRNAs derived from these two algorithms (Fig. 5E). In order to investigate the diagnostic ability of the characteristic miRNA, we plotted the ROC curve. The AUC value of ROC curve was 1.0, indicating that the characteristic miRNA could distinguish GBC and NGB samples powerfully (Fig. 5F). The result indicated that hsa-miR-4770 had potential diagnostic value in clinical practice.
The ceRNA network of characteristic miRNA
To explore the ceRNA regulation mechanisms of characteristic miRNA in GBC, we first screened for DE-mRNAs, DE-lncRNAs, and DE-circRNAs between GBC and NGB samples in GSE76633 and GSE100363 datasets. 2758 up-regulated mRNAs and 1547 down-regulated mRNAs in GBC were identified in the GSE76633 dataset (Fig. 6A and B). Meanwhile, 13 up-regulated lncRNAs and three down-regulated lncRNAs in GBC samples were mined (Fig. 6C and D). In the GSE100363 dataset, 585 up-regulated circRNAs and 443 down-regulated circRNAs in the GBC samples were determined (Fig. 6E and F). Using the Starbase database, we predicted 1166 mRNAs targeted by hsa-miR-4770, while 119 lncRNAs and 1572 circRNAs were predicted to target hsa-miR-4770. Since hsa-miR-4770 was down-regulated in the GBC group, we intersected the predicted targeting mRNAs, upstream lncRNAs, upstream circRNAs, and up-regulated mRNAs, lncRNAs, circRNAs respectively, to obtain a total of 211 mRNAs, two lncRNAs, and 48 circRNAs (Fig. 6G and I). Thus, a final ceRNA network containing 211 mRNAs, one miRNA, two lncRNAs, and 48 circRNAs was created (Fig. 7, Table S4).
To further understand the function and the pathways involved in targeting mRNAs in ceRNA regulatory network of the characteristic miRNA, we performed the corresponding functional enrichment analysis. A total of 80 biological process (BP) entries, 21 cellular components (CC) entries and five molecular functions (MF) entries, and one KEGG pathway were enriched (Tables S5, S6, S7, S8). The enriched GO terms with statistical differences (P<0.05) were listed (Fig. 8). We noted that the downstream genes regulated by hsa-miR-4770 were mainly involved in actin filament organization, cell-substrate adhesion, cell-matrix adhesion, regulation of reactive oxygen species metabolic process, glutamine metabolic process, and extracellular matrix (ECM) receptor interactions pathway-receptor interaction pathway.
Key genes in the ceRNA network
To explore the correlation between genes in the ceRNA network of hsa-miR-4770 and disease, we used varElect. Using the key words “gallbladder carcinoma” and “cancer”, we identified the top 10 genes, namely BRCA1, CHEK2, RB1, CASP8, PTGS2, CD44, KRT19, CDK1, PVT1, and MXRA5 were most correlated with the disease, in which PVT1 encoded lncRNA (Table 2). Therefore, we presented the ceRNA subnetwork composed of these 10 key genes and hsa-miR-4770 (Fig. 9A). In the network, PVT1 regulated the transcription of BRCA1, CHEK2, RB1, CASP8, PTGS2, CD44, CDK1, PVT1 and MXRA5 by targeting hsa-miR-4770. Besides, hsa-miR-4770-related TFs were further predicted according to the TransmiR database (Fig. 9B), suggesting the potential effects of EOMES, FOXA1, FOXH1, OTX2, REST on GBC by regulating hsa-miR-4770 (Table S9). In GBC, hsa-miR-4770 was down-regulated and the corresponding CD44, KRT19, RB1, BRCA1, PTGS2, CASP8, CDK1, CHEK2, MXRA5, and PVT1 were up-regulated (Fig. 9C and E). Then, we collected 10 NGB tissue samples and 10 GBC tissue samples from clinic for further validation of the expression of characteristic miRNA at transcriptional level. Consistent with the result of GSE104165 dataset, hsa-miR-4770 was significantly reduced in GBC samples compared with NGB samples (P=0.0147) (Fig. 9F).
In order to further investigate the interactions and functions of key genes in the ceRNA network of hsa-miR-4770, we performed a relevant analysis using the GeneMANIA database. The interactions of key genes and associated genes were shown. Among these genes, multiple genes were involved in cell cycle-related processes, p53 related pathway, and extrinsic apoptotic signaling pathway, which played an important role in the development of GBC (Fig. 10).
GBC originates from the gallbladder epithelial cells and stimulated by various factors such as gallbladder stones and chronic inflammation. GBC is the most common malignant tumor of hepato-biliary tract cancer, and has a high degree of malignancy, unsatisfactory survival, and prognosis . Because of the absence of specific symptoms in the early stages of GBC, the vast majority of patients are diagnosed at advanced stages, and the five-year survival rate is only 10% . Research on the diagnosis and treatment of GBC has never stopped, and miRNAs are receiving increasing attention. The miRNAs have been widely reported to play crucial roles in development, metastasis, epithelial-mesenchymal transition, and prognosis of GBC [9, 25,26,27,28].
Integrated bioinformatics approaches have helped to unravel the mechanisms and progression of cancer and have also facilitated research advances in cancer diagnosis, treatment and prognosis [29, 30]. In this study, we took the intersection of hub miRNAs identified by WGCNA and DE-miRNAs and obtained a total of 131 GBC-related DE-miRNAs. Those dysregulated miRNAs were mainly associated with VEGF signaling pathway and positive regulation of endothelial cell proliferation, which were considered as important for GBC progression. It has been shown that tumors can generate lymphatic vessels through the VEGF signaling pathway, thus affecting the distant metastasis of primary tumor cells . Next, we overlapped the miRNAs derived from LASSO and SVM-RFE and defined hsa-miR-4770 as the characteristic miRNA for GBC. The ROC curve showed that hsa-miR-4770 has prominent diagnostic efficacy, which can clearly distinguish GBC samples from NGB samples, it also further illustrates that hsa-miR-4770 can be used as the diagnostic factor.
It has been reported that hsa-miR-4770 expression significantly changed in a wide variety of tumors. Shi-Han Xiao et al. found that hsa-miR-4770 was upregulated in tumor sediments of colorectal cancer. Based on six miRNAs including hsa-miR-4770, they constructed the signature nomogram for preoperative prediction of tumor deposits in colorectal cancer . However, Xiao-Bing Wu et al. conducted a comparison of miRNAs between colorectal cancer at different stages and paired adjacent normal mucosa, and found that hsa-miR-4770 was significantly downregulated in stages II, II-III, and III . Coincidentally, another study of colorectal adenocarcinoma also found significant downregulation of hsa-miR-4770 . Yan-Yan Liu et al. identified that hsa-miR-4770 as one of the specific miRNAs in Helicobacter pylori (+) gastric cancer, showed a large number of interactions in the ceRNA network of Helicobacter pylori (+) gastric cancer, and downregulated significantly . Fang-Yun Bai et al. showed that hsa-miR-4770 was downregulated in gastrointestinal stromal tumor cells by qRT-PCR . Elena Vila Navarro et al. sequenced and validated the miRNAs used to detect pancreatic tumors and found that hsa-miR-4770 was significantly upregulated in pancreatic ductal adenocarcinoma, and identified hsa-miR-4770 as one of the novel markers with potential diagnostic or therapeutic value . In addition, the downregulation of hsa-miR-4770 expression was also reported in esophageal squamous cell carcinoma, endometric cancer, and laryngeal squamous cell carcinoma [38,39,40]. We uncovered that hsa-miR-4770 was significantly downregulated in GBC compared with normal tissues in the GSE104165 dataset. We further performed qRT-PCR validation on hsa-miR-4770 expression, which yielded the same trend of variation. Hence, after a series of analyses and validations, we identified hsa-miR-4770 for the first time as a characteristic miRNA of GBC.
Through prediction and differential analysis, we constructed a ceRNA network of hsa-miR-4770 containing 211 mRNAs, one miRNA, two lncRNAs, and 48 circRNAs. GO and KEGG pathway analysis has been widely used to evaluate the enriched biological functions. GO enrichment of mRNAs regulated by hsa-miR-4770 showed its involvement in regulating actin filament organization, cell-substrate adhesion, cell-matrix adhesion, regulation of reactive oxygen species metabolic process, glutamine metabolic process. Furthermore, KEGG pathway enrichment analysis suggested that mRNAs regulated by hsa-miR-4770 was involved in the ECM receptor interactions pathway. ECM is known to play a crucial role in regulating cell survival, proliferation, migration, and differentiation . These results revealed clues for the role of hsa-miR-4770 in GBC, and provided us with directions and ideas in subsequent studies.
We used GeneMANIA database to analyze the interaction and function of the above ceRNA sub-network, and found that hsa-miR-4770 target gene was involved in the cell cycle, DNA damage repair and p53 process. It is known that p53, as an oncogene, plays a key role in controlling cancer development and progression by regulating cell cycle arrest, apoptosis, senescence and DNA repair . The p53 pathway maintains genome stability by regulating cell cycle delays and DNA damage repair . It is thus clear that hsa-miR-4770 may participate in these biological processes by regulating its target genes, and then affect the occurrence and development of GBC, including the survival and prognosis of GBC patients.
We further constructed hsa-miR-4770 miRNA TFs network to predict the regulatory mechanism, and also paid attention to the research progress of TFs obtained. EOMES is mainly expressed in the nucleus of activated CD8+ T cells , which can enhance the proliferation and survival ability of tumor antigen-specific CD8+ T cells in immunotherapy, effectively control tumor growth or completely eliminate tumors . FOXA1 and FOXH1 are members of FOX family of transcription factors, FOXA1 is decreased in hepatocellular carcinoma , FOXH1 promotes cell proliferation, invasion and tumorigenesis by enhancing the Wnt/β-catenin pathway in lung cancer cells and tissues . OTX2 controls the regulatory landscape medulloblastoma through cooperative activity at enhancer elements and contributes to the expression of critical target genes . REST has been found to be an important negative transcriptional regulator in various cells. In non-nerve cells and tumors (such as lung, breast and colon), REST plays an anti-cancer role , but its research in GBC has not been reported. The above TFs provides a little reference value for us to preliminarily reveal the regulation mechanism of hsa-miR-4770.
In the aforementioned results, we have only preliminarily explored the expression changes, ceRNA network of has-miR-4770 and interaction mechanism. Due to insufficient clinical samples and funding, we only included 10 pairs of samples for qRT-PCR detection and validation of the characterized miRNA. We will conduct more experimental validation of the mechanisms involved in the regulation of hub genes by miRNAs and TFs in the future to further explore and clarify their potential mechanisms.
Taken together, we firstly identified hsa-miR-4770 as the characteristic miRNA for GBC based on the GEO database, and further constructed a ceRNA network for hsa-miR-4770. The regulatory pathways include cell cycle-related processes, p53, and extrinsic apoptotic signaling pathways. We preliminarily verified that hsa-miR-4770 expression was downregulated in GBC. Five key TFs involved in the regulation of hsa-miR-4770 were observed in the ceRNA network. This study provided some basis for understanding the underlying pathogenesis of GBC and identifying diagnostic biomarkers and potential therapeutic targets for GBC.
Availability of data and materials
Data analyzed in this study is publicly available from GEO database (Accession No.: GSE104165, GSE76633 and GSE100363).
Competing endogenous RNA
Differentially expressed miRNAs
Weighted gene co-expression network analysis
Least absolute shrinkage and selection operator
Support vector machine-recursive feature elimination
Ng noncoding RNAs
Quantitative real-time polymerase chain reaction
Gene Expression Omnibus
miRNA response elements
Differentially expressed mRNAs
Differentially expressed lncRNAs
Differentially expressed circRNAs
Receiver operating characteristic
Jiao K, Jiang W, Zhao C, Su D, Zhang H. Bmi-1 in gallbladder carcinoma: clinicopathology and mechanism of regulation of human gallbladder carcinoma proliferation. Oncol Lett. 2019;18(2):1365–71.
Sung H, Ferlay J, Siegel R, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global Cancer Statistics 2020: GLOBOCAN estimates of incidence and Mortality Worldwide for 36 cancers in 185 countries. Cancer J Clin. 2021;71(3):209–49.
Misra S, Chaturvedi A, Misra NC, Sharma ID. Carcinoma of the gallbladder. Lancet Oncol. 2003;4(3):167–76.
Cai Q, Wang S, Jin L, Weng M, Zhou D, Wang J, Tang Z, Quan Z. Long non-coding RNA GBCDRlnc1 induces chemoresistance of gallbladder cancer cells by activating autophagy. Mol Cancer. 2019;18(1):82.
Huang J, Peng J, Guo L. Non-coding RNA: a new tool for the diagnosis, prognosis, and therapy of small cell lung cancer. J Thorac Oncol. 2015;10(1):28–37.
Guo H, Ingolia NT, Weissman JS, Bartel DP. Mammalian microRNAs predominantly act to decrease target mRNA levels. Nature. 2010;466(7308):835–40.
Shu Y, Bao R, Jiang L, Wang Z, Wang X, Zhang F, Liang H, Li H, Ye Y, Xiang S, et al. MicroRNA-29c-5p suppresses gallbladder carcinoma progression by directly targeting CPEB4 and inhibiting the MAPK pathway. Cell Death Differ. 2017;24(3):445–57.
Shu Y-J, Bao R-F, Jiang L, Wang Z, Wang X-A, Zhang F, Liang H-B, Li H-F, Ye Y-Y, Xiang S-S. MicroRNA-29c-5p suppresses gallbladder carcinoma progression by directly targeting CPEB4 and inhibiting the MAPK pathway. Cell Death & Differentiation. 2017;24(3):445–57.
Chen J, Yu Y, Chen X, He Y, Hu Q, Li H, Han Q, Ren F, Li J, Li C, et al. MiR-139-5p is associated with poor prognosis and regulates glycolysis by repressing PKM2 in gallbladder carcinoma. Cell Prolif. 2018;51(6):e12510.
Qi X, Zhang DH, Wu N, Xiao JH, Wang X, Ma W. ceRNA in cancer: possible functions and clinical implications. J Med Genet. 2015;52(10):710–8.
Jiang J, Bi Y, Liu XP, Yu D, Yan X, Yao J, Liu T, Li S. To construct a ceRNA regulatory network as prognostic biomarkers for bladder cancer. J Cell Mol Med. 2020;24(9):5375–86.
Smyth G. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Kern F, Fehlmann T, Solomon J, Schwed L, Grammes N, Backes C, Van Keuren-Jensen K, Craig DW, Meese E, Keller A. miEAA 2.0: integrating multi-species microRNA enrichment analysis and workflow management systems. Nucleic Acids Res. 2020;48(W1):W521–w528.
Kern F, Fehlmann T, Solomon J, Schwed L, Grammes N, Backes C, Van Keuren-Jensen K, Craig DW, Meese E, Keller A. miEAA 2.0: integrating multi-species microRNA enrichment analysis and workflow management systems. Nucleic Acids Res. 2020;48(W1):W521–8.
Friedman J, Hastie T, Tibshirani R. Regularization Paths for generalized Linear Models via Coordinate Descent. J Stat Softw. 2010;33(1):1–22.
Huang M, Hung Y, Lee W, Li R, Jiang B. SVM-RFE based feature selection and Taguchi parameters optimization for multiclass SVM classifier. TheScientificWorldJournal. 2014;2014:795624.
Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez J, Müller M. pROC: an open-source package for R and S + to analyze and compare ROC curves. BMC Bioinformatics. 2011;12:77.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Yu G, Wang L, Han Y, He Q. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.
Stelzer G, Plaschkes I, Oz-Levi D, Alkelai A, Olender T, Zimmerman S, Twik M, Belinky F, Fishilevich S, Nudel R, et al. VarElect: the phenotype-based variation prioritizer of the GeneCards suite. BMC Genomics. 2016;17(2):195-206.
Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, Stein T, Nudel R, Lieder I, Mazor Y, et al. The GeneCards suite: from Gene Data Mining to Disease Genome sequence analyses. Curr Protocols Bioinf. 2016;54:13031–313033.
Tong Z, Cui Q, Wang J, Zhou Y. TransmiR v2.0: an updated transcription factor-microRNA regulation database. Nucleic Acids Res. 2019;47(D1):D253–d258.
Livak K, Schmittgen T. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods (San Diego Calif). 2001;25(4):402–8.
Shu YJ, Bao RF, Jiang L, Wang Z, Wang XA, Zhang F, Liang HB, Li HF, Ye YY, Xiang SS, et al. MicroRNA-29c-5p suppresses gallbladder carcinoma progression by directly targeting CPEB4 and inhibiting the MAPK pathway. Cell Death Differ. 2017;24(3):445–57.
Chang Y, Liu C, Yang J, Liu G, Feng F, Tang J, Hu L, Li L, Jiang F, Chen C, et al. MiR-20a triggers metastasis of gallbladder carcinoma. J Hepatol. 2013;59(3):518–27.
Chandra V, Kim JJ, Mittal B, Rai R. MicroRNA aberrations: an emerging field for gallbladder cancer management. World J Gastroenterol. 2016;22(5):1787.
Li Z, Yu X, Shen J, Law PT, Chan MT, Wu WK. MicroRNA expression and its implications for diagnosis and therapy of gallbladder cancer. Oncotarget. 2015;6(16):13914.
Zhou L, Du Y, Kong L, Zhang X, Chen Q. Identification of molecular target genes and key pathways in hepatocellular carcinoma by bioinformatics analysis. Onco Targets Ther. 2018;11:1861-9.
Chinnappan J, Akilandeswari R, Vidhya RV. Integrative Bioinformatics approaches to therapeutic gene target selection in various cancers for nitroglycerin. Sci Rep. 2021;11(1):22036.
Yan J, Xiao G, Yang C, Liu Q, Lv C, Yu X, Zhou Z, Lin S, Bai Z, Lin H, et al. Cancer-Associated Fibroblasts Promote Lymphatic Metastasis in Cholangiocarcinoma via the PDGF-BB/PDGFR-β Mediated Paracrine Signaling Network. Aging Dis. 2023;14(6).
Xiao S, Guo J, Zhang W, Hu X, Wang R, Chen Z, Lai C. A Six-microRNA signature Nomogram for Preoperative Prediction of Tumor deposits in Colorectal Cancer. Int J Gen Med. 2022;15:675.
Wu X, Li S, Xu X, Wu S, Chen R, Jiang Q, Li Y, Xu Y. The potential value of miR-1 and miR-374b as biomarkers for colorectal cancer. Int J Clin Exp Pathol. 2015;8(3):2840.
Wu X, Xu X, Li S, Wu S, Chen R, Jiang Q, Liu H, Sun Y, Li Y, Xu Y. Identification and validation of potential biomarkers for the detection of dysregulated microRNA by qPCR in patients with colorectal adenocarcinoma. PLoS ONE. 2015;10(3):e0120024.
Liu Y, Zhu J, Ma X, Han S, Xiao D, Jia Y, Wang Y. ceRNA network construction and comparison of gastric cancer with or without Helicobacter pylori infection. J Cell Physiol. 2019;234(5):7128–40.
Bai F, Zhang N, Fang W, He X, Zheng Y, Gu D. PCAT6 mediates cellular biological functions in gastrointestinal stromal tumor via upregulation of PRDX5 and activation of wnt pathway. Mol Carcinog. 2020;59(6):661–9.
Vila-Navarro E, Vila-Casadesús M, Moreira L, Duran-Sanchon S, Sinha R, Ginés À, Fernández-Esparrach G, Miquel R, Cuatrecasas M, Castells A. MicroRNAs for detection of pancreatic neoplasia: biomarker discovery by next-generation sequencing and validation in 2 independent cohorts. Ann Surg. 2017;265(6):1226.
Sang C, Chao C, Wang M, Zhang Y, Luo G, Zhang X. Identification and validation of hub microRNAs dysregulated in esophageal squamous cell carcinoma. Aging. 2020;12(10):9807–24.
Gao L, Nie X, Zhang W, Gou R, Hu Y, Qi Y, Li X, Liu Q, Liu J, Lin B. Identification of long noncoding RNA RP11-89K21. 1 and RP11-357H14. 17 as prognostic signature of endometrial carcinoma via integrated bioinformatics analysis. Cancer Cell Int. 2020;20:1.
Hui L, Wang J, Zhang J, Long J. lncRNA TMEM51-AS1 and RUSC1-AS1 function as ceRNAs for induction of laryngeal squamous cell carcinoma and prediction of prognosis. PeerJ. 2019;7:e7456.
Järveläinen H, Sainio A, Koulu M, Wight TN, Penttinen R. Extracellular matrix molecules: potential targets in pharmacotherapy. Pharmacol Rev. 2009;61(2):198–223.
Qin JJ, Li X, Hunt C, Wang W, Wang H, Zhang R. Natural products targeting the p53-MDM2 pathway and mutant p53: recent advances and implications in cancer medicine. Genes Dis. 2018;5(3):204–19.
Liu K, Zang Y, Guo X, Wei F, Yin J, Pang L, Chen D. The ∆133p53 isoform reduces Wtp53-induced Stimulation of DNA Pol γ Activity in the Presence and absence of D4T. Aging Dis. 2017;8(2):228–39.
Intlekofer AM, Takemoto N, Wherry EJ, Longworth SA, Northrup JT, Palanivel VR, Mullen AC, Gasink CR, Kaech SM, Miller JD. Effector and memory CD8 + T cell fate coupled by T-bet and eomesodermin. Nat Immunol. 2005;6(12):1236–44.
Furusawa A, Reiser J, Sadashivaiah K, Simpson H, Banerjee A. Eomesodermin increases survival and IL-2 responsiveness of tumor-specific CD8 + T cells in an adoptive transfer model of cancer immunotherapy. J Immunotherapy (Hagerstown, Md: 1997). 2018;41(2):53.
Dou C, Wang Y, Li C, Liu Z, Jia Y, Li Q, Yang W, Yao Y, Liu Q, Tu K. MicroRNA-212 suppresses tumor growth of human hepatocellular carcinoma by targeting FOXA1. Oncotarget. 2015;6(15):13216.
Zhang J, Zhang X, Yang S, Bao Y, Xu D, Liu L. FOXH1 promotes lung cancer progression by activating the Wnt/β-catenin signaling pathway. Cancer Cell Int. 2021;21(1):1–13.
Boulay G, Awad ME, Riggi N, Archer TC, Iyer S, Boonseng WE, Rossetti NE, Naigles B, Rengarajan S, Volorio A. OTX2 activity at distal regulatory elements shapes the chromatin landscape of group 3 medulloblastoma. Cancer Discov. 2017;7(3):288–301.
Li C, Zou H, Wang Z, Tang X, Fan X, Zhang K, Liu J, Li Z. REST, not REST4, is a risk factor associated with radiotherapy plus chemotherapy efficacy in glioma. Drug Des Devel Ther. 2018;12:1363–71.
Thank the public dataset providers for the data contribution.
This research was funded by the Yunnan Medical Discipline Leader Training Program (D-2019012), Yunnan Provincial Department of Science and Technology-Kunming Medical University Joint Funding for Applied Basic Research (202201AY070001-102) and Innovation Fund for Master’s Degree Students of the Second Affiliated Hospital of Kunming Medical University (2023S086).
Ethics approval and consent to participate
The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by the China Ethics Committee of Registering Clinical Trials (No: Chi ECRCT20190065) and informed consent was taken from all individual participants.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Heatmap plot for the enriched pathways by miEAA.
Total differentially expressed miRNAs in GBC and normal tissues.
The 131 GBC-related DE-miRNAs.
The miRNA enrichment analysis results of dysregulated miRNAs.
The ceRNA network of GBC containing 211 mRNAs, 1 miRNA, 2 lncRNAs, and 48 circRNAs.
Eighty biological process entries analysis of ceRNA regulatory network.
Twenty-one cellular components analysis of ceRNA regulatory network.
Five molecular functions entries analysis of ceRNA regulatory network.
One KEGG pathway analysis of ceRNA regulatory network.
The potential effects of EOMES, FOXA1, FOXH1, OTX2, REST on GBC by regulating hsa-miR-4770.
About this article
Cite this article
Shao, H., Zhu, J., Zhu, Y. et al. Identification of characteristic genes and construction of regulatory network in gallbladder carcinoma. BMC Med Genomics 16, 240 (2023). https://doi.org/10.1186/s12920-023-01663-z