Comprehensive analysis of the molecular mechanism for gastric cancer based on competitive endogenous RNA network
Hong-Jin Wu, Wei-Wei Dai, Li-Bo Wang, Jie Zhang, Cheng-Long Wang
Central Laboratory for Science and Technology, Longhua Hospital Shanghai University of Traditional Chinese Medicine, Shanghai, China
Correspondence Address:
Dr. Hong-Jin Wu
Central Laboratory of Science and Technology, Longhua Hospital, Shanghai University of Traditional Chinese Medicine, Shanghai 200032
China
Source of Support: None, Conflict of Interest: None
DOI: 10.4103/2311-8571.355010
Objective: To explore the regulatory mechanism of competitive endogenous RNAs (ceRNA) in gastric cancer (GC) and to predict the prognosis of GC. Materials and Methods: Expression profiles of long noncoding RNAs (lncRNAs), microRNAs (miRNAs), and mRNAs were obtained from The Cancer Genome Atlas platform. Differentially expressed RNAs (DERNAs) were screened to construct a lncRNA-miRNA-mRNA ceRNA network. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses were performed on the ceRNA network-related differentially expressed mRNAs (DEmRNAs). Next, the DERNAs were subjected to Cox regression and survival analyses to identify crucial prognostic factors for patients with GC. Results: We detected 1029 differentially expressed lncRNAs, 104 differentially expressed miRNAs, and 1659 DEmRNAs in patients with GC. Next, we performed bioinformatic analysis to construct the lncRNA-miRNA-mRNA ceRNA network, which included 10 miRNAs, 65 lncRNAs, and 10 mRNAs. Subsequently, KaplanMeier (K-M) analysis showed that the survival rate of the high-risk group was significantly lower than that of the low-risk group, and the area under the curve value of the receiver operating characteristic curve revealed that the polygenic model had good predictive ability. The results indicated that ADAMTS9-AS1, ATAD2, and CADM2 might be potential therapeutic targets and prognostic biomarkers for GC. Conclusions: Our study has implications for predicting prognosis and monitoring surveillance of GC and provides a new theoretical and experimental basis for the clinical prognosis of GC.
Keywords: Competitive endogenous RNAs network, gastric cancer, prognostic biomarkers, The Cancer Genome Atlas
Gastric cancer (GC) is a malignant tumor derived from the gastric mucosal epithelium, and its morbidity and mortality rank fifth and third, respectively, among lung, breast, colorectal, and prostate cancers.[1] GC has a higher incidence and mortality rate in East Asia.[2] Approximately 42% of all patients with GC are from China.[3] Due to atypical early symptoms, most patients with GC are already in an advanced stage of the disease when they are diagnosed.[4] Although some patients have received radical surgical treatment, they still have high recurrence and metastasis after surgery, with some of them developing advanced GC with poor prognosis.[5] Chemotherapy can prolong advanced GC patients' survival time and improve their quality of life, but the median survival time cannot exceed 1 year.[6] Therefore, it is still necessary to explore the mechanisms underlying the pathogenesis and prognosis of GC and to identify genes that affect GC development and have diagnostic and predictive value.
Bioinformatic analysis can be used to explore genes related to the occurrence, development, and prognosis of tumors.[7] Bioinformatics analysis can also resolve the inconsistent results caused by sample sizes or different microarray platforms.[8] Many studies have used bioinformatics to predict biomarkers for cancer treatment.[9],[10],[11],[12] Advances in postgenomic technology have allowed molecular biologists to study DNA (genome), mRNA (transcriptome), and protein sequences (proteome) in greater detail and to combine them for analysis in novel ways. Identification of molecular markers and expression profiles is used to predict tumor classification, diagnosis, and clinical outcomes. Although several cancer-related studies have been conducted, the biological functions of most mRNAs remain unclear. The genetic characteristics, individual differences, and complexity of the molecular mechanism of GC determine that its occurrence and development need to be summarized by gene groups or gene clusters, including the accurate assessment of prognosis of GC patients.[13]
To this end, we performed a bioinformatics analysis to furtherscreen some molecular indicators with higher specificity, to realize the early diagnosis of GC and individualized treatment based on more targets, and to improve the prognosis of GC. Long noncoding RNAs (lncRNAs) can regulate gene expression by indirectly interacting with microRNAs (miRNAs). However, the role of competitive endogenous RNAs (ceRNA) networks related to mRNAs, lncRNAs, and miRNAs in GC has not been fully elucidated.
Growing evidence shows that lncRNAs are involved in regulating the occurrence and development of GC through the mechanism of ceRNAs; however, there is less research on the function of lncRNAs as ceRNAs. Therefore, lncRNAs that function as ceRNAs require further exploration. In this study, we performed bioinformatics analysis to construct a lncRNA-miRNA-mRNA ceRNA network of GC. Then, through survival analysis based on the ceRNA network, we identified some independent prognostic indicators, including ADAMTS9-AS1, ATAD2, and CADM2, which are essential for survival prediction. The results of our study may provide more insights into the molecular mechanisms of GC.
Materials and MethodsPatients and The Cancer Genome Atlas More Details data retrieval
The RNA sequence data of 413 GC samples were retrieved from The Cancer Genome Atlas (TCGA) database (http://cancergenome.nih.gov) (up to September 20, 2020), which were divided into 381 tumor and 32 normal samples. The miRNA sequencing data of 497 GC cases were collected, including 452 GC tissue samples and 45 adjacent tissue samples. Using the ensemble database as a reference, 14166 lncRNAs and 19646 mRNAs were identified from the RNA SEQ expression profile data, and 1881 miRNAs were identified from the miRNA SEQ expression profile data for subsequent analysis. This study was conducted in accordance with the publication guidelines provided by TCGA (https://cancergenome.nih.gov/publication/publication guideline).
Identification of differentially expressed long noncoding RNAs, differentially expressed mRNAs, and differentially expressed microRNAs
We analyzed the differentially expressed lncRNAs (DElncRNAs), differentially expressed miRNAs (DEmiRNAs), and differentially expressed mRNAs (DEmRNAs) using the “edgeR” package of the R software (version 4.0.2). The false discovery rate (FDR) was used for all P values to correct the statistical significance of the multiple tests. |Fold change| ≥2 and FDR <0.05 were considered significant. The sampling hierarchical cluster analysis method was used to cluster the differentially expressed genes, and clustering results are shown by sampling heat map and volcano map.
Establishment of a competitive endogenous RNA regulatory network
We established the ceRNA control network for GC using the following steps: First, the interactions of DElncRNA-DEmiRNA were predicted using the miRcode database (http://www.mircode.org) and starBase v2.0 (http://starbase.sysu.edu.cn/). Next, miRTarBase (http://mirtarbase.mbc.nctu.edu.tw), miRDB (http://www.mirdb.org), and TargetScan (http://www.TargetScan.org) databases were used to retrieve DEmiRNA-targeted DEmRNAs. The target gene was intersected with the DEmRNAs. The regulatory relationships among DEmiRNAs, DElncRNAs, and DEmRNAs were obtained, and the interaction network diagram was generated using Cytoscape v3.6.0.
Differentially expressed mRNAs functional enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to identify LncRNAs with potential ceRNA function using the R software (version 4.0.2, The Bell Labs, New Jersey, USA). The critical value of GO analysis was P < 0.05. The results of the enrichment analysis were visualized using the GOplot and ggplot2 package in R software.
Construct the prognostic risk score
Univariate and multivariate Cox regression analyses were used to establish the prognostic risk assessment model using R software, and a predictive signature model of GC prognosis was established. We assessed the risk score of patients and classified them into high- and low-risk groups.
Survival analysis
Using the clinicopathological data of patients with GC in TCGA, the overall survival curve of differentially expressed RNAs (DERNAs) was generated by combining differential gene expression data. All data were processed by logarithmic transformation and median centralization and analyzed using the Kaplan–Meier survival method. P < 0.05 was considered to indicate that gene expression is related to prognosis.
Validation of prognostic long noncoding RNAs and differentially expressed mRNAs
The immunohistochemistry (IHC) method was used to validate the protein expression of prognostic genes based on the Human Protein Atlas database (HPA) (https://www.Proteinatlas.org) software.
Statistical processing
The threshold for differentially expressed genes was defined as | fold change| ≥2 and FDR < 0.05. The survival package of the R software was used for survival analysis, and differences were statistically significant at P < 0.05.
ResultsPatient characteristics
We obtained a complete dataset by deleting the samples with missing data. We then described the clinicopathological features of patients with GC [Table 1]. We analyzed the clinicopathological data of patients with GC, including sex, age, stage, cancer site, histological grade, lymph node stage, metastasis, and vital status.
Identification differentially expressed long noncoding RNAs, differentially expressed mRNAs and differentially expressed microRNAsys in gastric cancer
According to the threshold for differentially expressed genes (|Fold change| ≥2 and FDR <0.05), 1029 DElncRNAs of GC samples were identified compared to those in normal tissues, including 815 upregulated lncRNAs (79.2%) and 214 downregulated lncRNAs (20.8%). Furthermore, we identified 1659 DEmRNAs, including 916 (55.2%) upregulated mRNAs and 743 (44.8%) downregulated mRNAs using the edge R package. We also examined the differential expression profiles of miRNAs between the tumor and normal tissues. Finally, 104 DEmiRNAs were identified in patients with GC, including 87 upregulated (83.7%) and 17 downregulated miRNAs (16.3%). Using a volcano map and heatmap, the DElncRNAs, DEmRNAs, and DEmiRNAs are presented in [Figure 1]a, [Figure 1]b, [Figure 1]c. Hierarchical clustering analysis revealed differences in lncRNA, miRNA, and mRNA expression [Figure 1]d, [Figure 1]e, [Figure 1]f]. The top 10 upregulated and downregulated DElncRNAs, DEmiRNAs, and DEmRNAs between GC and normal tissues are listed in [Table 2].
Construction of a competitive endogenous RNA regulatory network in gastric cancer
To understand the role of DERNAs in GC, we established the DElncRNA-DEmiRNA-DEmRNA-related ceRNA regulatory network of GC. First, we obtained 1029 DElncRNAs from the miRcode database and then identified 168 pairs of interacting lncRNAs and miRNAs using the Perl program. We predicted that 10 of the retrieved DEmiRNAs could interact with the 65 DElncRNAs. Next, we identified miRNA-targeted mRNAs, including six miRNAs and 10 mRNAs, using the MiRTarBase, TargetScan, and miRDB databases. The ceRNA network, including 65 lncRNAs, 10 miRNAs, and 10 mRNAs, was constructed based on the interactions between RNAs [Figure 2]a and [Figure 2]b. Representative interactions among the genes are presented (lncRNAs, miRNAs, and mRNAs) [Table 3] and [Table 4].
Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis of Differentially expressed mRNAs
To further observe the functional characteristics of DEmRNAs involved in the ceRNA network, GO terms and KEGG pathway analyses were performed using R software [Table 5]. The results revealed that biological processes were primarily associated with cornification, keratinocyte differentiation, epidermal development, antimicrobial humoral response, epidermal cell differentiation, antimicrobial humoral immune response mediated by antimicrobials, whereas molecular functions were associated with receptor ligand activity, signaling receptor activator activity, hormone activity, and serine-type endopeptidase activity. KEGG pathway analysis showed that the genes were dominant in the neuroactive ligand-receptor interaction, protein digestion and absorption, and metabolism of xenobiotics by cytochrome P450 [Figure 3]a, [Figure 3]b, [Figure 3]c, [Figure 3]d.
Survival model construction and predictive accuracy evaluation
To identify the potential prognostic characteristics of DERNAs, the expression levels of 65 DElncRNAs, 10 DEmiRNAs, and 10 DEmRNAs in the ceRNA network were analyzed using univariate Cox proportional hazards regression. Nine DElncRNAs and three DEmRNAs were identified using univariate Cox analysis (P < 0.05) [Figure 4]a and [Figure 5]a and [Table 6], [Table 7]. The most significant survival-associated DElncRNAs and DEmRNAs were further analyzed by multivariate Cox regression, and the independent prognostic predictors for GC were identified [Figure 4]c and [Figure 5]c. The expression heat maps of five prognostic genes in the low-risk and high-risk groups are shown in [Figure 4]b and [Figure 5]b. Furthermore, the prognostic indices (PIs) for survival prediction were determined based on five lncRNAs and two mRNAs: ERVMER61-1 (P = 0.083), ADAMTS9-AS1 (P = 0.008), LINC00326 (P = 0.118), AL139002.1 (P = 0.043), AL391152.1 (P = 0.086), SERPINE1 (P < 0.001), and CADM2 (P = 0.02). The formulae for the two prognostic signatures were as follows: lncRNA-based prognostic index (PI) = (expression of ERVMER61-1 × 0.071908) + (expression of ADAMTS9-AS1 × 0.175152) + (expression of LINC00326 × 0.094621) + (expression of AL139002.1 × 0.103898) + (expression of AL391152.1 × 0.090519); As shown in [Figure 4]c, the concordance index of this prognostic model was 0.63; and the mRNA-based PI = (expression of SERPINE1 × 0.215545) + (expression of CADM2 × 0.08513); As shown in [Figure 5]c, the concordance index of this prognostic model was 0.61. The results of the K-M survival curve suggested significant differences in GC between the two groups [Figure 4]d and [Figure 5]d. Both projected PIs showed moderate prognostic assessment ability, with area under the curve values of 0.64 and 0.631 [Figure 4]e and [Figure 5]e. As shown in [Figure 4]f and [Figure 4]g and [Figure 5]f and [Figure 5]g, the scores assigned to each patient provided a good prognosis assessment. The expression heat map of the five prognostic genes in the low-and high-risk groups is shown in [Figure 4]b and [Figure 5]b.
The construction of prognostic predictors, including different types of RNA, could improve the accuracy of predicting the prognosis of GC. We performed K-M curve analysis on the above 9-lncRNA, 3-mRNA between the low-and high-expression groups [Figure 4]a and [Figure 5]a. Finally, we observed that four lncRNAs (ADAMTS9-AS1, LINC00326, AL139002.1, and AL391152.1) and two mRNAs (ATAD2 and SERPINE1) showed significant differences in survival prediction between the high and low expression groups (P < 0.05) [Figure 6]a, [Figure 6]b, [Figure 6]c, [Figure 6]d, [Figure 6]e, [Figure 6]f. Notably, low expression of ADAMTS9-AS1, LINC00326, AL139002.1, and SERPINE1 was associated with better prognosis of patients with GC. When the expression of the above genes was high, the survival time was relatively short in patients with GC, and the expression profiles of the above three lncRNAs (LINC00326, P = 0.125) and two mRNAs in tumor samples were significantly different from those in normal samples [Figure 7]a, [Figure 7]b, [Figure 7]c, [Figure 7]d, [Figure 7]e, [Figure 7]f.
Correlations between prognostic long noncoding RNAs s and differentially expressed mRNAs
The constructed ceRNA network showed potential interactions between DElncRNAs and mRNAs in GC. Based on multivariate Cox regression analysis, we found that ADAMTS9-AS1 was associated with the prognosis of patients with GC (hazard ratio = 1.9; 95% CI: 1.19–3.1; P = 0.008; [Figure 4]c). Theoretically, lncRNAs involved in lncRNA-miRNA-mRNA relationships positively regulate mRNA expression levels. To validate this mechanism in the context of GC, we performed a linear regression analysis on ATAMTS9-AS1 and 10 mRNAs that were related to the ceRNA subnetwork. We found three positive correlations between lncRNA-mRNA pairs and two negative correlations between lncRNA-mRNA pairs (r > 0.40, P < 0.05) [Figure 8]a, [Figure 8]b, [Figure 8]c, [Figure 8]d, [Figure 8]e.
Validation of prognostic RNA expression
To determine the clinical significance of GC-related prognostic genes in patients, ADAMTS9-AS1, ATAD2, and CADM2 were selected. The expression of these three genes was identified by HPA analysis. The results of IHC showed that ADAMTS9-AS1, ATAD2, and CADM2 were significantly upregulated in GC compared to normal tissues [Figure 9]a, [Figure 9]b, [Figure 9]c.
GC is a malignant tumor that poses a serious threat to global health. It ranks fifth in terms of incidence and third in terms of mortality worldwide.[14] Gastroscopy is the primary choice for examining and treating early stage GC, and the standard treatment for locally advanced GC is D2 gastrectomy combined with perioperative chemotherapy. However, most GC patients have a confirmed diagnosis at a relatively advanced stage due to lack of symptoms, which results in a less favorable prognosis.[15] Although cancer treatment has entered the era of precision medicine, the molecular mechanisms of GC pathogenesis and development have not been fully elucidated. Compared with other tumors, the research progress in accurate diagnosis and individualized treatment of GC is slow, and there is no ideal curative approach or prognosis prediction index.[16],[17] This may also be due to the lack of sensitive and accurate markers for early GC diagnosis. It is particularly urgent to identify highly sensitive and specific prognostic indicators for GC screening and early diagnosis. Bioinformatics and other technologies have been the most widely applicable and successful means of transforming molecular and genomic data into clinical practice. Recent studies have shown that lncRNAs play an important role in the development of stomach tumors.[18],[19] Increasing evidence shows that ceRNAs are related to the initiation, progression, invasion, and metastasis of many types of tumors.[20] However, the regulatory network of ceRNAs in GC remains unclear. It has been reported that ceRNAs include various types of RNA transcripts involved in the complex transcriptional regulation network in organisms, including mRNA, lncRNA, pseudogenes, and circRNAs. The region in which miRNA can bind is called miRNA response elements (MRE), and ceRNA can bind miRNAs through MREs to affect miRNA-induced gene silencing.[21] The construction of a GC-specific ceRNA regulatory network can provide ideas and directions for identifying potential targets for the diagnosis and treatment of GC. Furthermore, it may provide new biomarkers and potential effective therapeutic targets for GC. Recent studies have shown that lncRNAs are significantly associated with the occurrence and development of various malignant tumors, and may act as potential tumor suppressors, including GC[22] and pancreatic cancer.[23]
In the ceRNA network constructed in this study, we obtained 1029 lncRNAs, 104 miRNAs, and 1659 mRNAs and found that their expression was significantly different between gastric tumor and noncancerous tissues. Furthermore, we constructed a ceRNA network using bioinformatics tools combined with the results of DEmRNAs GO and KEGG functional analysis. According to the ceRNA mechanism, lncRNAs can regulate mRNA expression by directly interacting with miRNAs, which suggests that the biological function of lncRNAs may be similar to that of mRNAs. GO enrichment analysis showed that DEmRNA are involved in biological processes, including regulation of cornification, keratinocyte differentiation, epidermis development, digestion, and antimicrobial humoral immune response mediated by antimicrobials. The enriched KEGG pathways of mRNAs included neuroactive ligand-receptor interaction, bile secretion, protein digestion and absorption, and metabolism of xenobiotics by cytochrome P450, which were previously reported to be related to malignant tumors.[24],[25]
Recently, many scholars have proposed prognostic models based on lncRNAs,[26] miRNAs[27],[28] and mRNAs.[29] However, a prognostic model that includes different types of RNA for patients with GC is still missing. In this study, we demonstrated two different PIs, based on lncRNAs and mRNAs that could predict the overall survival and mortality of GC patients. It is worth noting that several genes were differentially expressed and related to the prognosis of GC. Particularly, four lncRNAs (ADAMTS9-AS1, LINC00326, AL139002.1, and AL391152.1) and two mRNAs (SERPINE1 and CADM2) were associated with the overall survival of GC patients. No miRNAs related to the prognosis of GC were found in the network, which may be due to the relatively small number of patients and insufficient database comparison.
Through univariate and multivariate Cox regression analyses, five lncRNAs (ERVMER61-1, ADAMTS9-AS1, AL139002.1, LINC00326, and AL391152.1) and two mRNAs (SERPINE1 and CADM2) were identified and used to establish a risk score model. The prediction efficiency of the model for the 5-year survival rate of patients with GC was 0.64 and 0.63, respectively, which has a medium degree of prediction efficiency according to the evaluation criteria. Moreover, screening hub genes in GC based on differentially expressed genes related to overall survival may be considered an independent prognostic index with high sensitivity and specificity, which has certain clinical significance. The results were also verified by analyzing the expression of these genes based on TCGA GC clinical database and by investigating the correlation of prognostic lncRNAs and DEmRNAs. The expression of ADAMTS9-AS1 and ATAD2 was significantly different between the normal and tumor tissues. Furthermore, we found that ADAMTS9-AS1 was most closely associated with ATAD2 and CADM2. The findings from this study indicate that ADAMTS9-AS1, ATAD2 and CADM2 could be considered as potential prognostic biomarkers in GC.
ADAMTS9-AS1 is an important anti-angiogenic factor that regulates tumor progression and metastasis by regulating fibroblast growth factor and vascular endothelial growth factor.[30] ADAMTS9 plays an important role in the inhibition of cancer genes in esophageal, gastric, breast, and other tumors via DNA methylation.[31] Recently, many studies have shown that the lncRNA ADAMTS9-AS1 plays a role in promoting cancer as a miR-206 sponge, and it could be used as a novel biomarker of GC.[32] Several studies have verified that promoter methylation of CADM2 is one of the reasons for the hypoexpression of CADM2 in tumor tissues.[33],[34] ATAD2 has AAA and bromine domains.[35] Studies have shown that ATAD2 and C-myc form a co-expression complex that is involved in the development of tumors and plays an important role in tumor proliferation.[36] HPA is a useful database for the analysis of protein expression between cancer and normal tissues based on IHC.[37] In the present study, HPA was used to determine the expression of prognostic genes. The results showed that the expression of ADAMTS9, CADM2, and ATAD2 was significantly upregulated in GC tissues compared to that in normal tissues.
Above all, these studies provided a powerful basis for the prediction model, and we hypothesized that the abnormal expression of these RNAs might be related to the development of gastric carcinomas and might impact the prognosis of patients with GC, which has potential clinical research value in the diagnosis and prognosis evaluation of GC. However, this study has several limitations. First, the potential mechanisms of these RNAs in GC should be further verified to determine their feasibility as diagnostic biomarkers. Second, despite the limited power of detecting individual events, the proposed PI has promising clinical implications. In summary, we performed this study to demonstrate that the ceRNA network and other important biomarkers may improve the prediction of prognosis of GC. This regulatory network can provide a reference for studying the mechanisms of GC.
ConclusionsWe established a ceRNA network by analyzing lncRNAs, miRNAs, and mRNAs in GC tissues from TCGA database. Importantly, the high expression of ADAMTS9-AS1, ATAD2, and CADM2 is closely related to the prognosis of GC. It can be used as a potential biomarker to assess the survival and prognosis of patients with GC.
Financial support and sponsorship
This study was supported by the Shanghai Natural Science Foundation of China (Grant No. 16ZR1447300).
Conflicts of interest
There are no conflicts of interest.
References
Comments (0)