Integrative analysis of the gastric cancer long non-coding RNA-associated competing endogenous RNA network.

Gastric cancer (GC) is a common type of cancer, and identification of novel diagnostic biomarkers associated with this disease is important. The present study aimed to identify novel diagnostic biomarkers associated with the prognosis of GC, using an integrated bioinformatics approach. Differentially expressed long non-coding RNAs (lncRNAs) associated with GC were identified using Gene Expression Omnibus datasets (GSE58828, GSE72305 and GSE99416) and The Cancer Genome Atlas database. A competing endogenous RNA network that incorporated five lncRNAs [long intergenic non-protein coding RNA 501 (LINC00501), LINC00365, SOX21 antisense divergent transcript 1 (SOX21-AS1), GK intronic transcript 1 (GK-IT1) and DLEU7 antisense RNA 1 (DLEU7-AS1)], 29 microRNAs and 114 mRNAs was constructed. Gene Ontology and protein-protein interaction network analyses revealed that these lncRNAs may be involved in 'biological regulation', 'metabolic process', 'cell communication', 'developmental process', 'cell proliferation', 'reproduction' and the 'cell cycle'. The results of receiver operating characteristic curve analysis demonstrated that LINC00501 (AUC=0.819), LINC00365 (AUC=0.580), SOX21-AS1 (AUC=0.736), GK-IT1 (AUC=0.823) and DLEU7-AS1 (AUC=0.932) had the potential to become valuable diagnostic biomarkers for GC. Associations with clinicopathological characteristics demonstrated that LINC00501 expression was significantly associated with sex (P=0.015) and tumor grade (P=0.022). Furthermore, LINC00365 expression was significantly associated with lymph node metastasis (P=0.025). Gene set enrichment analysis revealed that LINC00501, LINC00365 and SOX21-AS1 were enriched in signaling pathways associated with GC. Reverse transcription-quantitative PCR analysis demonstrated that LINC00501 expression (P=0.043) was significantly upregulated in GC tissues, whereas the expression levels of LINC00365 (P=0.033) and SOX21-AS1 (P=0.037) were significantly downregulated in GC tissues. Taken together, the results of the present study suggest that LINC00501, LINC00365, SOX21-AS1, GK-IT1 and DLEU7-AS1 may be used as novel diagnostic biomarkers for GC, and may be functionally associated with GC development and progression. Copyright: © Jiang et al.

Introduction

Among the most common types of cancer, gastric cancer (GC) was the second leading cause of cancer-associated mortality worldwide until 2018, with ~1 million new cases of GC and 782,685 GC-associated mortalities occurring in 2018 (1). Radical surgery remains the most reliable treatment for GC (2). Due to advancements in surgical techniques, chemotherapy and radiotherapy, the 5-year survival rate of patients with early stage GC is >95% (2). However, early stage GC is often asymptomatic and rarely detected, and thus, >70% of patients with GC eventually develop advanced GC and may be ineligible for surgical resection (2). Thus, it is important to identify novel biomarkers associated with GC to improve the diagnosis of this deadly disease.

Long non-coding RNAs (lncRNAs) are RNA molecules that are >200 nucleotides in length, lack protein coding potential, and can regulate the migration, survival and proliferation of cancer cells (3,4). Several studies have reported that specific lncRNAs are dysregulated in GC, including urothelial cancer-associated 1 (5), antisense RNA in the INK4 locus (6) and AC093818.1 (7). At the functional level, lncRNAs serve as competing endogenous RNAs (ceRNAs) that bind to specific microRNAs (miRNAs/miRs) in a competitive manner, thereby upregulating the expression levels of miRNA target genes (8). Several studies have demonstrated the interactions between lncRNAs, miRNAs and mRNAs, which influence the development and progression of different types of cancer. For example, lncRNA HNF1A antisense RNA 1 binds to miR-30b-3p in GC, which promotes activation of the PI3K/AKT signaling pathway (9). In addition, lncRNA DCST1 antisense RNA 1 regulates the survival, proliferation and invasion of GC cells by competitively binding and sequestering miR-605-3p (10). It has been reported that lncRNA metallothionein 1J pseudogene can control F-box and WD repeat domain containing 7 expression by competitively targeting miR-92a-3p in GC cells (11). Thus, broader analyses of lncRNA-miRNA-mRNA ceRNA networks may enable researchers to further understand complex GC-associated gene interactions. In addition, a ceRNA model can offer unique insights to understand the putative roles of several uncharacterized lncRNAs in GC development and progression.

In the present study, an integrated systematic analysis of lncRNAs and mRNA expression patterns in datasets of patients with GC was performed. Multiple publicly available Gene Expression Omnibus (GEO) datasets (GSE58828, GSE72305 and GSE99416) were normalized and integrated to identify genes that were differentially expressed between tumor tissues and adjacent normal tissues from patients with GC. In addition, The Cancer Genome Atlas (TCGA) database was used to confirm the differential expression profiles of lncRNAs in GC. A ceRNA network was constructed using these data. Gene Ontology (GO) and protein-protein interaction (PPI) analyses of mRNAs in the GC-related ceRNA network were performed to further investigate the association between the lncRNAs and GC in the ceRNA network. Receiver operating characteristic (ROC) curve analysis was performed to determine the diagnostic value of the central lncRNAs. Subsequently, associations with clinicopathological characteristics were assessed, and gene set enrichment analysis (GSEA) of the central lncRNAs was performed.

The present study aimed to identify novel biomarkers for the diagnosis of GC and provide insight into the mechanistic basis of GC development and progression.

Materials and methods

Microarray datasets

Studies from the GEO database (https://www.ncbi.nlm.nih.gov/geo) were considered eligible according to the following criteria: i) Studies with GC tissue samples; ii) studies with information on technology and platform used for studies. Based on these criteria, the GSE58828, GSE72305 (12) and GSE99416 (13) GC datasets were downloaded from the GEO database (14). Details of each microarray study are presented in Table I. GC RNA expression profile data were obtained from TCGA database (http://tcga-data.nci.nih.gov) (15). An overview of the bioinformatics analysis performed to assess the GEO and TCGA datasets is presented in Fig. S1.

Table I.

Details of gastric cancer studies and associated microarray datasets from the Gene Expression Omnibus database.

		Sample size, n						Sex, n
Series accession no.	Platform	Total	Tumor	Adjacent	DEGs, n	Upregulated, n	Downregulated, n	Female	Male
GSE58828	GPL15314	6	3	3	9,121	4,869	4,252	2	4
GSE72305	GPL15314	12	10	2	324	108	216	NA	NA
GSE99416	GPL16956	12	6	6	2,651	1,261	1,390	4	8

DEGs, differentially expressed genes.

Differential expression analysis

All three datasets were initially integrated to increase the overall sample numbers (19 tumor samples and 11 normal samples), and introduction of unreliable results was avoided by batch normalizing the three datasets, using the sva package (16) (version 3.38.0; http://www.bioconductor.org/packages/release/bioc/html/sva.html) and limma package (17) (version 3.46.0; http://bioconductor.org/packages/release/bioc/html/limma.html) in R software (version 3.5.1; www.R-project.org). Subsequently, tumor and normal tissue samples were compared using the limma package in R software (17), and differential RNA expression was detected using the following criteria: |log fold-change (FC)|>1 and adjusted P<0.05. A volcano map was constructed using the pheatmap package (version 1.0.12; http://cran.r-project.org/web/packages/pheatmap/index.html) in R software.

ceRNA network construction

A ceRNA network was constructed to visualize lncRNA-miRNA-mRNA interactions, based on the ceRNA model wherein lncRNAs can bind and sequester miRNAs, thereby altering their ability to influence mRNA translation (8). The miRcode database (version 11; http://www.mircode.org) (18) was used to identify intersecting lncRNA target miRNAs, while miRDB (version 5.0; http://mirdb.org), TargetScan 7.2 (version 7.2; http://www.targetscan.org) (19) and miRTarBase (version 7.0; http://mirtarbase.cuhk.edu.cn/php/download.php) (20) were used to predict interactions between miRNAs and mRNAs.

Genes identified in at least two databases were considered miRNA targets. The target mRNAs were further intersected with dysregulated mRNAs in GEO GC samples. Subsequently, the ceRNA network was constructed by combining the lncRNA-miRNA interactomes and miRNA-mRNA interactomes, which was visualized using Cytoscape software (version 3.6.1; http://cytoscape.org) (21).

GO and PPI analysis

To assess the functional relevance of differentially expressed lncRNAs in the present study, GO analysis of the molecular functions, cellular components and biological processes, for which these genes were enriched, was performed using Webgestalt 2019 (http://www.webgestalt.org) (22). In addition, PPI network data of the mRNAs were collected from the Search Tool for the Retrieval of Interacting Genes (STRING) database (version 10.0; http://string-db.org) (23), and the PPI network was constructed using Cytoscape software.

ROC curve analysis

ROC curve analysis was performed to determine the diagnostic values of the differentially expressed lncRNAs in the ceRNA network, using the survival ROC package (version 1.66.0; http://bioconductor.org/packages/ROC) in R software. An area under the curve (AUC) value of >0.50 indicated evidence for diagnosis.

Correlations between lncRNA expression and clinicopathological characteristics

To further evaluate the clinical value of the central lncRNAs in GC, the associations between lncRNA expression and clinicopathological characteristics, including age, sex, grade, stage, T classification, N classification, M classification of TCGA data were analyzed using the χ2 test. Unpaired Student's t-test was used to compare differences between two groups. The tumor grading system was based on the Goseki histological grading of GC (24).

GSEA

For GSEA, cut-off values for each lncRNA were determined based on the median expression levels in the TCGA-STAD database (https://portal.gdc.cancer.gov/projects/TCGA-STAD). The ‘c2.cp.kegg.v7.0.symbols.gmt’ dataset was used for 1,000 gene set permutations per analysis for each lncRNA, using GSEA software (version 4.0.1; http://software.broadinstitute.org/gsea/downloads.jsp). The present study focused on pathways with nominal P-value <0.05 and selected the most significantly enriched signaling pathways based on their size >50. The results were generated using the ggplot2 package (version 3.3.3; http://CRAN.R-project.org/package=ggplot2) in R software.

Tissue sample collection

A total of 15 GC tumor tissues and paired adjacent normal tissues (5 cm away from the tumor margin) were collected from patients with GC who underwent radical primary tumor excision at The First Affiliated Hospital of Chongqing Medical University (Chongqing, China) between December 2016 and January 2019. The inclusion criteria were as follows: i) Patients had received an initial diagnosis of GC according to the results of gastroscopy and pathological biopsy; ii) all patients with GC were evaluated according to TNM stage (25) and iii) there were no significant differences in blood routine test indexes, liver function indexes and renal function indexes of all patients with GC. The exclusion criteria were as follows: i) patients with poor general condition and unable to tolerate related examinations; ii) concurrently diagnosed with other tumors and patients with secondary GC and iii) patients with GC who had received radiotherapy, chemotherapy or other treatments prior to surgery. The patients included 9 men and 6 women, with a mean age of 64.6±4.9 years (age range, 50–75 years). Tissue samples were stored at −80°C until subsequent experimentation. The present study was approved by the Research Ethics Committee of the First Affiliated Hospital of Chongqing Medical University (approval no. 2016140) and written informed consent was provided by all patients prior to the study start.

RT-qPCR

A total of 15 pairs of tumor tissues and adjacent normal tissues from patients with GC were used for RT-qPCR validation. Total RNA was extracted from tissues using TRIzol® reagent (Tiangen Biotech Co., Ltd.). The concentration of RNA was measured using a Nanodrop 2000 instrument (Bio-Rad Laboratories, Inc.). Total RNA (500 ng RNA) was reverse transcribed into cDNA using the PrimeScript RT Reagent kit with gDNA Eraser (Takara Biotechnology Co., Ltd.). qPCR was subsequently performed using specific primers (TsingKe Biological Technology) and TB Green® Premix ExTaq™ II (Takara Biotechnology Co., Ltd.). The following primer sequences were used for qPCR: LINC00501 forward, 5′-GAACAATGACCGGGGAACAG-3′ and reverse, 5′-TTCTTCCTTTGTGCTTCCGC-3′; LINC00365 forward, 5′-AGCTGCTCATCCTTCCTCAG-3′ and reverse, 5′-ACACAGGTGCCAAAATCCAC-3′; SOX21-AS1 forward, 5′-GAGGTGCTGCAGGAGAGTTA-3′ and reverse, 5′-ACTCTCCACTCGCCTAAACC-3′; DLEU7-AS1 forward, 5′-AACAAATTTGGGGCACTGCT-3′ and reverse, 5′-CACCAAAGCACGGAAGGTAG-3′; and GK-IT1 forward, 5′-CTGAGGTTGGGAGTTCGAGAC-3′ and reverse, 5′-GGATTACAGGCATGAGCCAC-3′; and GAPDH forward, 5′-GGTCTCCTCTGACTTCAACA-3′ and reverse, 5′-GTGAGGGTCTCTCTCTTCCT-3′. The temperature protocol for RT was as follows: 37°C for 15 min, 85°C for 5 sec and 4°C storage. The resulting cDNA product was stored at −20°C. The following thermocycling conditions were used for qPCR: Initial denaturation at 95°C for 3 min, followed by 40 cycles of 95°C for 1 sec and 58°C for 30 sec. lncRNA expression levels were calculated using the 2−∆∆Cq method (26) and normalized to the internal reference gene GAPDH.

Statistical analysis

Statistical analysis was performed using R software (version 3.5.1; www.R-project.org). All experiments were performed in triplicate and data are presented as the mean ± standard deviation. Paired Student's t-test was used to compare differences between lncRNA expression levels in tumor tissues and adjacent normal tissues. P<0.05 was considered to indicate a statistically significant difference.

Results

Pooled analysis of GC-related gene expression profiles in the GEO datasets

Following reannotation, 25,476, 22,053 and 32,386 genes were obtained in the GSE58828 (GPL15314), GSE72305 (GPL15314) and GSE99416 (GPL16956) datasets, respectively. Genes that were differentially expressed between GC tissues and normal tissues (|log FC|>1 and false discovery rate <0.05) were identified using the limma package. In total, 9,121 genes were differentially expressed in the GSE58828 dataset (4,869 upregulated genes and 4,252 downregulated genes), while 324 genes were differentially expressed in the GSE72305 dataset (108 upregulated genes and 216 downregulated genes) and 2,651 genes were differentially expressed in the GSE99416 dataset (1,261 upregulated genes and 1,390 downregulated genes). These differentially expressed genes are presented in volcano plots in Fig. 1. Following batch normalization, a total of 465 genes were differentially expressed between GC tissues and normal tissues. Overall, the present study identified 48 lncRNAs and 175 mRNAs that were upregulated, and 45 lncRNAs and 197 mRNAs that were downregulated in tumor tissues. The top 20 upregulated and downregulated lncRNAs (Fig. 2A) and mRNAs (Fig. 2B) in these datasets were identified, and further information regarding all the lncRNAs and mRNAs are presented in Tables SI and SII.

Figure 1.

Differentially expressed gastric cancer-related genes in the Gene Expression Omnibus database. Red dots represent upregulated genes (adjusted P<0.05 and log FC>1), green dots represent downregulated genes (adjusted P<0.05 and log FC<-1); and black dots represent genes that were not identified to be significantly differentially expressed. FC, fold-change; adj.P.Val, adjusted P-value.

Figure 2.

Top 20 gastric cancer-related differentially expressed genes in the Gene Expression Omnibus database. (A) Differentially expressed long non-coding RNAs. (B) Differentially expressed mRNAs. Red represents upregulated genes, while green represents downregulated genes. The number indicates the fold-change in each dataset.

TCGA-based validation of differential lncRNA expression

GC sample RNA expression profiles were downloaded from TCGA database to validate the differential lncRNA expression profiles detected in the GEO dataset analysis. Comparisons of the GEO and TCGA datasets identified 25 differentially expressed lncRNAs (Table II).

Table II.

Differentially expressed long non-coding RNAs in the Gene Expression Omnibus and The Cancer Genome Atlas databases.

lncRNA	Regulation	logFC	Adjusted P-value
PGM5-AS1	Down	−3.601	<0.001
LHX5-AS1	Down	−4.367	<0.001
LINC02163	Up	5.976	<0.001
DUXAP8	Up	3.582	<0.001
LINC00982	Down	−2.251	<0.001
PCAT18	Down	−3.098	<0.001
DLEU7-AS1	Up	2.599	<0.001
LINC00582	Down	−2.196	<0.001
FEZF1-AS1	Up	4.654	<0.001
DLGAP1-AS2	Up	2.174	<0.001
LINC00365	Down	−2.210	<0.001
NCOA7-AS1	Down	−2.018	<0.001
LINC02404	Down	−2.802	<0.001
LINC02447	Down	−1.205	<0.001
SOX21-AS1	Down	−2.296	<0.001
CASC19	Up	2.354	<0.001
UNC5B-AS1	Down	−1.660	<0.001
LINC00501	Up	2.184	<0.001
LINC01133	Down	−1.296	<0.001
LINC00853	Up	1.141	<0.001
LINC01985	Down	−1.304	<0.001
GK-IT1	Up	1.337	<0.001
HRAT92	Up	1.366	<0.001
LINC01589	Down	−1.061	<0.001
PICSAR	Up	1.947	<0.001

lncRNA, long non-coding RNA; FC, fold-change.

ceRNA network

Based on the lncRNA-miRNA, circRNA-miRNA and miRNA-mRNA predicted interactions, a ceRNA network was constructed. The network incorporated five lncRNAs (LINC00501, LINC00365, SOX21-AS1, GK-IT1 and DLEU7-AS1), 29 miRNAs and 114 mRNAs (Fig. 3). Overall, 29 specific miRNAs were predicted to target five specific lncRNAs (Table III), and 26 specific miRNAs (hsa-miR-125b-5p, hsa-miR-129-5p and hsa-miR-455-5p did not target mRNAs from the intersecting specific mRNAs) were identified to interact with the 114 intersecting mRNAs (Table IV).

Figure 3.

Competing endogenous RNA network. Red arrows represent long non-coding RNAs, green diamonds represent microRNAs and blue circles represent mRNAs. miR, microRNAs.

Table III.

Specific lncRNAs that target specific miRNAs.

lncRNAs	miRNAs
LINC00501	miR-301b-3p, miR-4295, miR-3666, miR-206, miR-613, miR-429, miR-23b-3p, miR-24-3p, miR-363-3p, miR-338-3p, miR-455-5p, miR-129-5p, miR-490-3p
LINC00365	miR-17-5p, miR-20b-5p, miR-429, miR-590-5p, miR-761, miR-3619-5p, miR-216b-5p, miR-363-3p, miR-338-3p, miR-449c-5p, miR-125a-5p, miR-125b-5p, miR-129-5p
SOX21-AS1	miR-301b-3p, miR-4295, miR-3666, miR-212-3p, miR-761, miR-3619-5p, miR-107, miR-338-3p, miR-125a-5p, miR-125b-5p, miR-455-5p, miR-129-5p
GK-IT1	miR-4770, miR-24-3p, miR-129-5p
DLEU7-AS1	miR-507, miR-142-3p, miR-761, miR-3619-5p, miR-23b-3p, miR-27a-3p, miR-107, miR-338-3p, miR-375

lncRNA, long non-coding RNA; miRNA/miR, microRNA; LINC00, long intergenic non-protein coding RNA; DLEU7-AS1, DLEU7 antisense RNA 1; GK-IT1, GK intronic transcript 1; SOX21-AS1, SOX21 antisense divergent transcript 1.

Table IV.

Specific miRNAs that target specific mRNAs.

miRNAs	mRNAs
miR-107	CYP2C8, DEPDC1B, FOXC1, CREBRF, GPRC5A, SYNDIG1, SERPINB5, KIF23, RPS6KB1,
miR-125a-5p	SWSAP1, UBE2J1, SWSAP1, UBE2J1, AGTR1, DMRTA2, GK, ADHFE1, ITIH5, RPS6KB1, CAMK2N1, REP15, NALCN, PDE3A, CKAP2L, GIF, NEUROD2, FANCI, PRKCB, ZIC2
miR-142-3p	ADAMTS3, CDC6, FLVCR1
miR-17-5p	PRKCB, KIF23, MSTO1, TMEM100, NAA50, LPGAT1, PBK, SYT16, AUNIP, ATP1B3, RHOV, BAMBI, MUC17, MMP3, UBE2J1, AADAC, RND3, STIL, NETO2, KIF26B, FOXC1, IGFBP5, PON2, TGM2, FAM83D, NMUR2
miR-206	RSPO3, SLC2A12, NETO2, PABPC1L2B, NALCN, DGKH
miR-20b-5p	NETO2, KIF26B, KIF23, TMEM100, LPGAT1, BAMBI, PON2, UBE2J1
miR-212-3p	SOWAHA, TLN2, OCLN, MUC13, CCNB1, DUSP9
miR-216b-5p	ANKRD22, ZIC5, FADS1, SOX11, RND3, ANXA10, TRAM1L1
miR-23b-3p	VSNL1, NAA50, LPGAT1, ALDH1A2, ZIC5, CA2
miR-24-3p	CCNB1, RSPO3, FSCN1, SULT2A1, CDX2, CD3D, STIL, NETO2, CDKN2A, DSC2, TLN2, OSMR
miR-27a-3p	GRIN2D, PRKCB, HRK, RPS6KB1, CREBRF, WISP1, RND3, HOXA10, PAPSS2, NAA50, SOWAHA, DEPDC1B, ONECUT2
miR-301b-3p	GDA, LGALSL, ZIC5, IGFBP5, NAA50
miR-338-3p	ONECUT2, GPRC5A, LGALSL, SERPINA5, CEACAM6
miR-3619-5p	NAA50, LPGAT1, PLCB2, P2RX6, SCAMP4, LTF, NRG4, FSCN1
miR-363-3p	FZD10, OSMR, DSC2, KLHL29, SYNDIG1, PAPSS2, SOX11, TRIP13, ZIC5, SOSTDC1, SYNPO2
miR-3666	NAA50
miR-375	NETO2, CLDN1, TGM2, MOSPD3
miR-429	ADAMTS3, NOG, RND3, OCLN, TRAM1L1, NALCN, FLVCR1, FSCN1, KLHL29, PRKCB, ONECUT2, KCTD8, TLN2, RPS6KB1, PMAIP1, KYNU
miR-4295	NAA50
miR-449c-5p	PAPSS2, ACADL, NMUR2, TRIM73, ATP1B3, CARNS1, MLANA, ACAD8, SERPINB5, KLHL29, NETO2, RPS6KB1, DEPDC1B, SOX11, AGR3
miR-4770	IGFBP5
miR-490-3p	SOX11, GCNT3
miR-507	DGKH, DSC2, RPS6KB1, LPGAT1, WISP1, RABGEF1
miR-590-5p	FBXL13, VSNL1, IL1B, ALDH1A1
miR-613	DGKH, RSPO3, NALCN, NETO2, PABPC1L2B, SLC2A12
miR-761	SCAMP4, NAA50, FSCN1, P2RX6, NRG4, LTF, LPGAT1

miRNA/miR, microRNA.

Functional enrichment and interaction analyses

To further investigate the association between the five lncRNAs and GC, GO and PPI network analyses of the 114 mRNAs within the constructed ceRNA network were performed. The results demonstrated that the mRNAs were involved in biological processes, including ‘biological regulation’, ‘metabolic process’, ‘response to stimulus’, ‘multicellular organismal process’, ‘cell communication’, ‘developmental process’, ‘localization’, ‘cellular component organization’, ‘cell proliferation’, ‘reproduction’ and ‘multi-organism process growth’ (Fig. 4A). The mRNAs were primarily enriched for cellular component terms, including ‘membrane’, ‘cytosol’, ‘nucleus’, ‘endomembrane system’, ‘vesicle’, ‘extracellular space’, ‘membrane-enclosed lumen’ and ‘protein-containing complex’ (Fig. 4B). For molecular functions terms, the mRNAs were concentrated in ‘protein binding’, ‘ion binding’, ‘nucleotide binding’, ‘hydrolase activity’, ‘nucleic acid binding’, ‘transferase activity’, ‘enzyme regulator activity’, ‘molecular transducer activity’, ‘transporter activity’ and ‘lipid binding’ (Fig. 4C).

Figure 4.

Functional enrichment analysis of differentially expressed genes. Gene Ontology analysis of significantly enriched (A) biological processes, (B) cellular components and (C) molecular functions. (D) Protein-protein interaction network.

A PPI network incorporating the mRNAs in the ceRNA network was constructed using the STRING database. The results demonstrated that the lncRNAs in the ceRNA network may be involved in regulating the cell cycle via cyclin B1 (CCNB1), family with sequence similarity 83 member D (FAM83D) and cell division cycle 6 (CDC6) in GC (Fig. 4D).

Diagnostic values of the lncRNAs

ROC curve analysis was performed to determine the diagnostic values of the central lncRNAs in the ceRNA network. The results generated AUC values of 0.819, 0.932, 0.823, 0.580 and 0.736 for LINC00501, DLEU7-AS1, GK-IT1, LINC00365 and SOX21-AS1, respectively (Fig. 5), which suggests that the five lncRNAs have diagnostic value for patients with GC.

Figure 5.

ROC curve analysis of the five long non-coding RNAs. The red line represents the sensitivity curve, while the black line represents the identifying line. ROC, receiver operating characteristic; AUC, area under the curve; LINC00, long intergenic non-protein coding RNA; SOX21-AS1, SOX21 antisense divergent transcript 1; GK-IT1, GK intronic transcript 1; DLEU7-AS1, DLEU7 antisense RNA 1.

Associations between lncRNAs and clinicopathological characteristics

The present study further assessed the clinical value of the central lncRNAs in GC, and associations between lncRNA expression and clinicopathological characteristics were analyzed using TCGA database. The results demonstrated that LINC00501 expression was significantly associated with sex (P=0.015; Fig. 6A) and tumor grade (P=0.022; Fig. 6B). Furthermore, LINC00365 expression was significantly associated with lymph node metastasis (P=0.025; Fig. 6C). However, due to the limited sample size, data of associations between other lncRNAs expression and clinicopathological characteristics were non-significant (Table SIII).

Figure 6.

Associations between long non-coding RNAs and clinicopathological characteristics. Association between LINC00501 expression and (A) sex and (B) tumor grade. (C) Association between LINC00365 expression and lymph node metastasis. LINC00, long intergenic non-protein coding RNA.

Identification of lncRNA-associated signaling pathways

GSEA was performed to identify signaling pathways that are differentially activated in GC. The present study focused on signaling pathways that had nominal P<0.05 and selected the most significantly enriched signaling pathways based on their SIZE >50. The results demonstrated that LINC00501 was enriched in the TGF-beta, mTOR, MAPK and WNT signaling pathways (Fig. 7A), while LINC00365 was enriched in the MAPK signaling pathway (Fig. 7B). In addition, SOX21-AS1 was enriched in the cell cycle (Fig. 7C). Taken together, these results suggest that these lncRNAs may modulate these signaling pathways, thereby regulate the development and progression of GC.

Figure 7.

GSEA pathways of the long non-coding RNAs. GSEA pathways of (A) LINC00501, (B) LINC00365 and (C) SOX21-AS1. GSEA, gene set enrichment analysis; LINC00, long intergenic non-protein coding RNA; SOX21-AS1, SOX21 antisense divergent transcript 1.

RT-qPCR verification of differential lncRNA expression in GC tissues

To confirm that the central lncRNAs were differentially expressed in GC tissues, RT-qPCR analysis was performed to detect their expression levels in 15 pairs of GC tissues and adjacent normal tissues. Consistent with the results of the GEO and TCGA dataset analyses, the results demonstrated that the expression levels of LINC00501, DLEU7-AS1 and GK-IT1 were upregulated in GC tumor tissues, whereas the expression levels of LINC00365 and SOX21-AS1 were downregulated in GC tumor tissues (Fig. 8). Notably, LINC00501 expression was increased 5.47-fold (P=0.043), LINC00365 expression was decreased 3.67-fold (P=0.033) and SOX21-AS1 expression was decreased 3.24-fold (P=0.037). Although the differences in DLEU7-AS1 and GK-IT1 expression were not significant, they exhibited an upward trend consistent with the bioinformatics analyses.

Figure 8.

Reverse transcription-quantitative PCR-based validation of the expression levels of five long non-coding RNAs in gastric cancer tissues. LINC00, long intergenic non-protein coding RNA; DLEU7-AS1, DLEU7 antisense RNA 1; GK-IT1, GK intronic transcript 1; SOX21-AS1, SOX21 antisense divergent transcript 1.

Discussion

Increasing evidence suggest that lncRNAs are involved in the progression of different types of cancer (27–30). However, the etiology of GC remains partly unknown. In addition, the roles of lncRNAs in GC remain unclear. The present study analyzed multiple GEO and TCGA datasets to comprehensively identify lncRNAs associated with GC development or progression. By constructing a ceRNA network, performing functional enrichment analyses and ROC analysis, assessing associations with clinicopathological characteristics, GSEA and RT-qPCR verification, five lncRNAs (LINC00501, LINC00365, SOX21-AS1, GK-IT1 and DLEU7-AS1) were identified that may be functionally associated with GC development, and thus may be used as novel diagnostic biomarkers for patients with GC.

In the present study, five lncRNAs (LINC00501, LINC00365, SOX21-AS1, GK-IT1 and DLEU7-AS1), 29 miRNAs and 114 mRNAs were incorporated into a ceRNA network. Among the included lncRNAs, LINC00501 and LINC00365 were associated with the most miRNAs, suggesting that they may be key regulators of GC due to their ceRNA functionality. Previous studies have demonstrated that a number of the selected mRNAs and predicted miRNAs, including miR-301b-3p (31), miR-4295 (32), miR-206 (33), miR-613 (34), miR-429 (35), forkhead box C1 (36), G protein-coupled receptor class C group 5 member A (37) and kinesin family member 23 (38), are associated with GC. GO analysis revealed that the 114 mRNAs within this network may be involved in biological processes, including ‘biological regulation’, ‘metabolic process’, ‘response to stimulus’, ‘multicellular organismal process’, ‘cell communication’, ‘developmental process’, ‘localization’, ‘cellular component organization’, ‘cell proliferation’, ‘reproduction’ and ‘multi-organism process growth’. The mRNAs were primarily enriched for cellular component terms, including ‘membrane’, ‘cytosol’, ‘nucleus’, ‘endomembrane system’, ‘vesicle’, ‘extracellular space’, ‘membrane-enclosed lumen’ and ‘protein-containing complex’. For molecular functions terms, the mRNAs were concentrated in ‘protein binding’, ‘ion binding’, ‘nucleotide binding’, ‘hydrolase activity’, ‘nucleic acid binding’, ‘transferase activity’, ‘enzyme regulator activity’, ‘molecular transducer activity’, ‘ transporter activity’ and ‘lipid binding’. PPI analyses revealed that these lncRNAs may promote or inhibit the development of GC by targeting a range of cell cycle-related mRNAs, such as CCNB1 (39,40), FAM83D (41) and CDC6 (42). The ceRNA hypothesis posits that lncRNAs can competitively bind miRNAs, thereby altering their ability to influence mRNA translation (8). Recent studies have demonstrated that ceRNAs play important functional roles in GC (43,44). In the present study, a ceRNA network was constructed, GO and PPI analyses were performed to elucidate the mechanisms regarding how these lncRNAs exert their functions. The results suggest LINC00501 may be involved in cell proliferation by targeting mRNA fascin actin-bundling protein 1 (FSCN1) via miR-429. Notably, a previous study demonstrated that miR-429 acts as a tumor suppressor by targeting FSCN1 in GC (45), which was consistent with the results of the present study. The results of the present study also suggest that LINC00365 may be involved in the cell cycle by targeting mRNA FAM83D via miR-17-5p. In addition, SOX21-AS1 and GK-IT1 may be involved in the cell cycle by targeting mRNA CCNB1 via miR-212-3p and miR-24-3p, respectively. Furthermore, DLEU7-AS1 may play a role in carcinogenesis by regulating the cell cycle through miR-142-3p/CDC6. In the present study, ROC analysis revealed that AUC values for LINC00501, LINC00365, SOX21-AS1, GK-IT1 and DLEU7-AS1 were all >0.50, which suggests that these five lncRNAs have the potential to serve as valuable diagnostic biomarkers for GC. Notably, analysis of the associations with clinicopathological characteristics revealed that LINC00501 expression was associated with sex and tumor grade. Furthermore, LINC00365 expression was upregulated in N1 stage samples compared with that in N0 stage samples, suggesting that LINC00365 may be an indicator for early stage GC diagnosis without lymph node metastasis. GSEA of the three lncRNAs indicated that LINC00501 may modulate GC progression via the TGF-beta, MAPK, WNT and mTOR signaling pathways, LINC00365 via MAPK signaling and SOX21-AS1 via the cell cycle. Several studies have highlighted the important roles of TGF-beta signaling pathway (46), MAPK signaling pathway (47), WNT signaling pathway (48) and the cell cycle (49) in GC.

LINC00501 has been identified as a prognostic factor associated with the overall survival of patients with hepatocellular carcinoma (HCC) (50–52). LINC00365 influences the Wnt/β-catenin signaling pathway, which modulates colorectal cancer (CRC) progression (53), and breast cancer cell viability may be regulated by the LINC00365-secretoglobin family 2A member 1 (SCGB2A1) axis, which targets NF-κB signaling (54). In GC, the expression levels of LINC00365 and SCGB2A1 are downregulated in tumor tissues, which is associated with a shorter survival time (55). Notably, this LINC00365/SCGB2A1 axis and associated NF-κB suppression are associated with GC progression (55). In the present study, SOX21-AS1 was demonstrated to be associated with GC in the ceRNA network. Previous studies have reported that SOX21-AS1 is associated with different types of cancer, including cervical cancer (56,57), HCC (58), lung adenocarcinoma (59), CRC (60) and oral cancer (61). Previous studies have demonstrated that SOX21-AS1 can function as a ceRNA to sequester miRNAs, which in turn modulates gene expression (57,60). For example, SOX21-AS1 participates in cervical cancer progression by competitively binding miR-7/voltage dependent anion channel 1 (57), while it sequesters miR-145 in CRC, which promotes tumor progression via the enhanced expression of myosin VI (60). SOX21-AS1 knockdown in nephroblastoma cells in vitro disrupts the proliferation and colony formation of these cells through a mechanism associated with p57 upregulation, and results in cell cycle arrest (62). In addition, SOX21-AS1 is associated with HCC progression and patient prognosis via a mechanism associated with p21 epigenetic silencing (58).

Although the differences in DLEU7-AS1 and GK-IT1 expression were not significant, they exhibited an upward trend consistent with the present bioinformatics analyses. This lack of significance may be attributed to the limited sample size in the present study. Thus, future studies with larger sample sizes are required to confirm the results presented here. A previous study demonstrated that DLEU7-AS1 expression is upregulated in CRC, and is associated with CRC stage, metastasis and poor patient prognosis (63). From a mechanistic perspective, DLEU7-AS1 is considered to modulate the Wnt/β-catenin signaling pathway in CRC cells, which influences the ability of these cells to proliferate and invade proximal or distal tissues (63). GK-IT1 has been reported to be positively associated with overall survival in patients with esophageal adenocarcinoma (64).

Several studies have investigated ceRNA networks for GC, which have aimed to screen lncRNAs involved in GC via analysis of gene expression data obtained from TCGA database (65–67). Different original non-coding RNA microarray data were used for the analysis, thus different lncRNAs were identified as biomarkers for GC. The present study aimed to screen novel lncRNAs involved in GC via analysis of gene expression data obtained from both the GEO and TCGA databases, which makes the results more reliable and repeatable. However, the present study is not without limitations. Although the expression levels of the five lncRNAs were verified in GC tissues, the lack of validation of these biomarkers in serum is a major limitation of the present study. Furthermore, future studies are required to confirm the functions of the identified lncRNAs in GC.

In conclusion, integrated analysis of TCGA and GEO datasets and a series of analyses identified five differentially expressed lncRNAs (LINC00501, LINC00365, SOX21-AS1, GK-IT1 and DLEU7-AS1), which may be used as novel diagnostic biomarkers associated with GC, and may be functionally associated with GC development and progression. Thus, the present study provides novel insight into the mechanistic basis and biological functions of lncRNAs in GC.