Diagnostic model of combined ceRNA and DNA methylation related genes in esophageal carcinoma

Data source

The data analyzed in this study involved miRNA, mRNA, DNA methylation and clinical data, all of which were derived from the GEO (Barrett et al., 2007) database at NCBI (https://www.ncbi.nlm.nih.gov/geo/) and humans Cancer Genome Database TCGA (https://www.cancer.gov/) (Akbani et al., 2014). Table 1 listed the specific GEO data numbers used in this study (Li et al., 2014).

Screening of differentially expressed miRNAs

The R package (Affy, version 1.52.0) was used for background expression value correction and data normalization of the original data in each dataset (Yuen, Zhu & Leung, 2018). The probes in each file were then annotated according to the appropriate platform annotation file. Probes without matching genetic symbols were removed. When different probes map to the same gene, the average value of all the probes mapped to the gene was the final expression value of the gene. Through literature research and analysis, the miRNA datasets of ESCC in this study were GSE59973 and GSE114110 (Wen et al., 2018). Differential expression analysis was performed and the differential miRNA selection should meet p value < 0.05 and log2—FC—>1.

Differential methylation analysis

In order to confirm the results, we also downloaded the methylation (platform: Illumina HumanMethylation450 BeadChip) and expression (IlluminaHiSeq) microarray data from TCGA database for validation. Through DNA research and analysis, the DNA methylation dataset in this study was determined to be GSE52826 (Li et al., 2014), and the DNA methylation difference analysis was performed with GEO2R (Cao et al., 2019a). The adjusted P value less than 0.05, and delta expression value either greater than 1 (up-regulated gene) or less than −1 (down-regulated gene) were as the cut-off value of the expression chip data.

Target prediction

The 12 miRNAs selected in the above steps were used for target prediction, and the target prediction software RNA22 was used (Jabbar et al., 2019; Wang et al., 2019), which is target prediction software for predicting microRNA binding sites based on sequence characteristics and can be a good theory of ceRNA hypothesis mechanism.

Hypergeometric test

Hypergeometric testing is the most commonly used prediction method in the ceRNA mechanism (Tay, Rinn & Pandolfi, 2014). After the hypergeometric test calculation (Doi, Takahashi & Kawasaki, 2017), the result p value less than 0.05 is the potential ceRNA pair. The specific formula is as follows: ${\displaystyle P=1-F\left(x∕N,K,M\right)=1-\sum _{t=0}^{x-1}\frac{\left(\frac{K}{t}\right)\left(\frac{N-K}{M-t}\right)}{\left(\frac{N}{M}\right)}}$

Pearson correlation coefficient calculation

The gene expression matrix of GSE120356 was used to calculate the Pearson correlation coefficient (Ma et al., 2019c). Based on the ceRNA mechanism, the two RNAs of the ceRNA pair have a co-expression effect, which can be well supported for subsequent analysis. Therefore, Pearson correlation was taken here and 0.7 had a highly credible co-expression trend (Li et al., 2019a; Li et al., 2019b). Finally, the ceRNA that meets the hypergeometric test threshold was crossed with the ceRNA under the Pearson threshold to obtain the final ceRNA pairs.

Construction of ceRNA network

The selected differential mRNA, lncRNA, and differential miRNA were paired in the miRcode database. Target gene prediction was performed on the selected differential miRNAs using the starBase (Han et al., 2020) online software (http://starbase.sysu.edu.cn/), and the target genes were predicted by miRDB (Rigden & Fernández, 2020), miRTarBase (Liu, Huang & Luo, 2020), and TargetScan (Dong et al., 2020) databases, thereby obtaining ceRNA regulation of lncRNA-miRNA-mRNA network, using Cytoscape v3.5.1 software (Qi et al., 2019) for mapping.

Functional and pathway enrichment analyses

Gene Ontology (GO) is a comprehensive information tool providing gene function of individual genomic products, including three aspects: molecular function (MF), biological process (BP) and cellular component (CC). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for advanced biological functions and utilities. Both analyses and annotations are based on the DAVID database (https://DAVID.ncifcrf.gov/) (Kumar et al., 2019) and are used to explore and understand the biological significance of specific gene lists. In this study, GO and KEGG analysis of DEGs were based on standard error discovery rate (FDR) <0.05.

Identification of prognostic lncRNAs from ceRNA network

With the survival time of all esophageal cancers in the TCGA database, the Survival package in the R language was used to perform survival analysis on the differential lncRNA (Shi et al., 2020b). The difference was statistically significant with P ≤ 0.05.

Samples collection and ethics statement

In this study, we collected 10 pairs of ESCC specimens from patients in the First Affiliated Hospital of China Medical University from 2015 to 2018. All patients were confirmed histologically. All patients participating in the study signed the informed consent form and the study was approved by the Institutional Ethics Committee of China Medical University.

RNA extraction and reverse transcription- quantitative polymerase chain reaction analysis

Total RNA of tissues was extracted, and the sample was fully ground with liquid nitrogen, added 1ml Trizol solution, mixed, and placed at room temperature for 5 min to fully crack. Then, 14,000 g centrifuged at 4 °C for 15 min, and the RNA was divided into three layers, in the upper water phase, and transferred to another new RNase free EP tube. Later, isopropanol was added, gently and thoroughly mixed, and allowed to stand at room temperature for 10 min. At 4 °C, the samples were centrifuged at 14,000 g for 10 min, RNA precipitation was collected and washed twice with 75% ethanol. According to the amount of precipitation, an appropriate amount of DEPC water was added to dissolve the precipitation.

RNA samples were reverse transcribed into cDNA using TransScript One Step gDNA Removal and cDNA Synthesis SuperMix. RT qPCR was performed using TB Green Premix Ex Taq. The PCR conditions included an initial step at 95 °C for 10 min, followed by 40 cycles of amplification and quantification (95 °C for 15 s, 60 °C for 1 min, and 60 °C for 1 min). GAPDH was used as an endogenous control for normalization. The sequences of the primers used for RT qPCR were as follows:

HCG22 forward, 5′ CCTGGGGAGAGGTGTCATTT 3′  and reverse, 5′ TGGTCTCTGGGTGCTTAGTG 3′; RNF217-AS1 forward, 5′  TGGGAATGACAGCAGAAAG 3′  and reverse, 5′  TCCGCAGAGTGAACAAGAA 3′; HCP5 forward, 5′  GATGACTATGGGGTGAGGGG 3′  and reverse, 5′  TATGGAGATGAGGTGTGCCG 3′; BBOX1-AS1 forward, 5′  GGCACATTTGGAAGTT 3′  and reverse, 5′ TCAGGGTAACCGTAGC 3′.

Analysis of differentially expressed miRNAs

The miRNA datasets used in this study were GSE59973 and GSE114110, respectively. The threshold of selecting differential miRNAs was p.value < 0.05 and log2—FC—>1. In the volcano maps of differential miRNAs of GSE59973 and GSE114110, the black regions in the graph represented the distribution of non-differentiated miRNAs, the red fractions showed significantly differentially up-regulated miRNAs, and the green fractions represented significantly differentially down-regulated miRNAs (Figs. 1A, 1B). At the same time, the differentially up-regulated and down-regulated miRNAs obtained from the two datasets were separately taken for intersection. Finally, three differentially expressed down-regulated miRNAs were obtained, including hsa-mir- 139-5p, hsa-mir-1 and hsa-mir-133b. There were nine up-regulated miRNAs, including hsa-mir-1246, hsa-mir-34c-5p, hsa-mir-455-5p, hsa-mir-455-3p, hsa-mir-146b-5p, hsa-mir-3651 and hsa-mir-429, as shown in Figs. 1C and 1D. These 12 important miRNA differences will be used for subsequent target prediction and ceRNA correlation analysis.

CeRNA network construction

The gene expression matrix of GSE120356 was calculated by Pearson correlation coefficient. The ceRNA pair with r value more than 0.7 and the ceRNA pair with hypergeometric test value less than 0.05 were selected to determine the final ceRNA relationship for a total of 11,696 pairs. It can be seen that there were complex ceRNAs in ESCC and the network, as shown in Fig. 2, was the ceRNA network constructed in this study. The red dots represent lncRNA, the green dots represent mRNA, and the blue dots represent other types of RNA. The size of the dots indicates the degree of node degrees. It can be seen that the ceRNA network constructed in this study has a clear trend of ceRNA regulation between lncRNA.

Analysis of ceRNA combined with DNA methylation

The mRNA and lncRNA in the ceRNA network were combined with DNA methylation analysis to find ceRNAs with differential DNA methylation in ESCC. Figure 3A showed the intersection of ceRNA and differentially methylated genes and it could be seen that in the ceRNA network, there were 643 lncRNAs with methylation differences and 147 mRNAs with methylation differences. At the same time, we extracted 64 of the 147 mRNAs that form a ceRNA pair with lncRNA, such as KCNA3, USP44, OPLAH, SMTN, TTC6, COL27A1, SYNE2, LHX1, NRG1 and XKR4. These 64 mRNAs were subjected to GO-BP and KEGG pathway enrichment analysis to discover the biological regulation processes in which these lncRNAs were mainly involved in the ceRNA network. Figures 3B and 3C showed the KEGG pathway bubble diagram and the GO-BP pathway strip diagram of the 64 ceRNAs, respectively. The GO-BP pathway was only visualized by selecting the pathway with pvalue <0.01. It can be seen that for the KEGG pathway results, these significantly enriched pathways were mainly concentrated in the FoxO signaling pathway and related pathways such as glutamate synapses and EGFR tyrosine kinases, and these pathways have been reported to play an important role in ESCC (Zhang et al., 2017; Kashyap & Abdel-Rahman, 2018; Hong et al., 2015). For the bar graph of GO-BP, these mRNAs were mainly involved in the biological processes of the corresponding cascade of calcineurin. The above analysis indirectly indicates that lncRNA in the ESCC ceRNA network indirectly regulates the biological activities and life processes of ESCC.

Prognosis related lncRNAs

The intersection of 643 lncRNAs in the ESCC ceRNA network and the differentially expressed lncRNAs in GSE120356 screened 10 important lncRNAs, namely CASC9, SOX21-AS1, HCP5, HCG22, RNF217-AS1, CALML3-AS1, LINC00491, BBOX1-AS1, C5orf66, ZFPM2-AS1, shown in Fig. 4A. To explore the ceRNA regulatory network in which these 10 lncRNAs were involved in ESCC, a subset of the ceRNA networks from which 10 important lncRNAs were extracted, shown in Fig. 4B. The orange dots in the figure indicate these 10 important lncRNAs, the green dots represent mRNAs, and the blue dots represent important lncRNAs other than these 10 important lncRNAs. It can be seen that these 10 important lncRNAs are mainly related to other lncRNAs in this study. There is a clear trend of ceRNA regulation between ncRNA-interacting lncRNAs.

It was found by analysis that four lncRNAs in these lncRNAs have important effects on the survival prognosis of ESCC. The four lncRNAs were HCG22, RNF217-AS1, HCP5, BBOX1-AS1, and the survival curves were shown in shown in Figs. 5A–5D. At the same time, the survival model was analyzed by ROC curve, and the survival analysis of the next three years was predicted, in which the AUC value reached 0.908 (Figs. 5E, 5F).

Validation of differential lncRNAs by RT-qPCR

The present study aimed to determine whether the differentially expressed lncRNAs identified in the microarray analysis were up-regulated or down-regulated in clinical esophageal cancer patients. Esophageal cancer specimens and non tumor epithelial tissues were obtained, and the differential lncRNAs were validated with RT-qPCR (Fig. 6A). The experimental results showed that RNF217-AS1, HCP5 and ZFPM2-AS1 expression was significantly up-regulated in ESCC samples compared with healthy tissues, while HCG22 expression was lower than negative controls. These results were consistent with the result of TCGA analysis, as shown in Fig. 6B. Taking together, it was revealed that these differential lncRNAs may be involved in the mechanism of ESCC.