Information from 33 datasets was downloaded from the University of California, Santa Cruz (UCSC) Xena The Cancer Genome Atlas (TCGA) hub (https://xena.ucsc.edu) (Goldman et al., 2020). The dataset included 407 gastric cancer tissues, 32 normal tissues, and 9951 other tumors from different types of cancers, and the RNA expression matrix was transformed to log₂ (FPKM+1). Ensemble transcript IDs were converted to their corresponding GENCODE v19 using the Gene Transfer Format (GTF) annotation files from humans. The enhancer RNAs and target gene information were obtained from the putative literature, which was previously identified by the PreSTIGE method (Corradin et al., 2014).

To avoid bias, patients with a survival time of less than 1 month were excluded. A total of 375 patients passed the quality control and were used in the following analysis as shown in Table 1. The survival-associated eRNAs were screened using the Cox regression model, with age, sex, and tumor stage adjusted as covariates. We set p < 0.05 as the significance cut-off value.

The R package maftools were used to compare the mutant frequencies of significantly mutated genes between CDK6-AS1 as high- and low-expression groups. Mutant types including frame shift, deletion, splice site, frameshift insertion, missense mutation, nonsense mutation, multiple hist, and in-frame deletion were considered in the analysis.

Gene Ontology (GO) functional analysis was performed using the clusterProfiler package in R software, and a Reactome pathway analysis of eRNA-related coding genes was performed based on co-expression analysis. Specifically, the GO analysis revealed the function of the biology process (BP), cell component (CC), and molecular function (MF). To avoid accumulation of type-I errors, enrichment items meeting the false discovery rate (FDR) < 0.05 were considered significant.

The expression data of CDK6-AS1 and its target gene CDK6 were obtained at the pan-cancer level as previously described, and patients were classified into low- and high-expression groups according to the median value of CDK6-AS1 expression. The Cox regression method was used to compare the overall survival difference between the two groups. Covariates of sex, age, and tumor stage were adjusted in the Cox model, and Spearman’s coefficient was applied to the correlation analysis.

Univariate Cox regression analysis and Kaplan-Meier analysis were used to screen 6 genes co-expressed with CDK6-AS1. Gastric cancer patients in the TCGA data set were randomly divided into the training set and internal test set. The aforementioned six genes were used for the LASSOCox regression analysis. By using the cross-verification error curve, the best tuning parameter λ is selected through the minimum 10-fold cross-verification in the training set. Based on these six genes, a risk-scoring model is established. Risk score = 0.2637*CTHRC1 + 0.0132*PFN2 + 0.1384*PRSS35 + 0.0355*RTN4 - 0.072*SMPD3 + 0.5459*SYCP2L.

GC Patients in the internal test set and an external cohort were divided into a high-risk group and a low-risk group by the optimal cut-off value of the risk score. The overall survival (OS) rates between the high-risk and the low-risk groups were analyzed by the Kaplan–Meier OS analysis. A two-sided log-rank p < 0.05 was considered significant. The time-dependent prognostic value of the prognostic signature was evaluated using the R package“survival ROC.” Area under the curve (AUC) values were used to evaluate the time-dependent prognostic values of the prognostic signature. An AUC >0.60 was considered to be acceptable.

To evaluate the relationship between tumor-infiltrating lymphocytes (TIL) and the expression of CDK6-AS1 in gastric cancer, we estimated the expressed fraction of TIL cells using the ssGSEA algorithm by comparing the gastric cancer gene expression matrix with those of the signatures from nine reported TIL cell types (Li et al., 2016). The relationship of the proportion matrix for the nine TIL cells with CDK6-AS1 was calculated by Spearman’s correlation analysis.

The R package pRRophetic (Geeleher et al., 2017), based on the pharmacogenomics database of the Cancer Genome Project (CGP) cell line data and the Cancer Cell Line Encyclopedia (CCLE), was used to predict chemotherapeutic sensitivity for gastric cancer by an estimation of IC50 (the maximal inhibitory concentration). Default settings were used for the prediction model, including “stomach cancer” for reference tissue type and “cvFold = 10” for ridge regression model training.

R software (Version 3.6.2) was used to perform analyses in this study. The statistical results are expressed as mean ± standard deviation (M ± SD), and the data comparison of the two groups was analyzed with the Wilcoxon rank-sum test. A value of p < 0.05 was used to determine the statistical significance.

Twenty-three eRNAs were identified, eight of which met the criteria (Spearman r ≥ 0.3 and FDR < 0.05) and were included (Supplementary Table S1). Of these, CDK6-AS1 exhibited the lowest Cox model p-value and was therefore considered a candidate marker. Patients in the CDK6-AS1 high-expression group had a shorter survival than those in the low-expression group (3-year OS:HR = 1.68, p = 3.84 × 10⁻³; 5-year OS:HR = 1.62, p = 5.64 × 10⁻³, Figures 1A,B). In addition, CDK6-AS1 shows a higher expression in unpaired and paired tumor tissues compared to normal tissues (unpaired: p = 8.00 × 10⁻³, paired: p = 0.046, Figures 1C,D). A positive correlation between CDK6-AS1 and its target gene CDK6 was observed (Spearman r = 0.38, p = 1.68 × 10⁻¹⁴). The connections between the clinical features of gastric cancer patients and the CDK6-AS1 expression were further investigated. It was found that CDK6-AS1 had a higher expression in patients aged <60 years (p = 0.022, Figure 2A). CDK6-AS1 was significantly linked to the clinical stage (III vs. II, p = 0.048, Figure 2C). Other clinical characteristics were not clearly correlated with the CDK6-AS1 expression (p > 0.05, Figures 2B, D–H). As driver gene mutations are crucial to tumor growth, the frequencies of significantly mutated genes were compared between the high- and low-CDK6-AS1 expression groups. It was noted that several classic gastric cancer driver genes were more frequently mutated in the low-CDK6-AS1 expression group than in the high-expression group, such as ARID1A and PIK3CA (Figures 3A,B, Supplementary Table S2).

To further explore the function and related pathways CDK6-AS1 involved in gastric cancer, a co-expression analysis between CDK6-AS1 and other protein-coding genes in 375 TCGA gastric cancer cases was performed. It was found that 595 transcripts presented a significant correlation with CDK6-AS1 (Spearman r > 0.30 & FDR <0.05). GO and Reactome enrichment analyses were performed. Results from TOP10GO pathways in BP, MF, and CC are shown in Figure 4A, Supplementary Table S3. In BP, the terms were mainly related to RNA transportation. In MF, the terms were related to exoribonuclease activity. In CC, the term is involved in the nuclear chromosomal region. Enrichment from the Reactome pathway database indicated that CDK6-AS1 related co-expressed genes were mainly involved in the cell cycle and mitosis (Figure 4B, Supplementary Table S5), and key signals, in tumor cell proliferation, as shown by Spearman’s correlation (r > 0.3 are showed in Supplementary Table S5). Taken together, CDK6-AS1 and its related genes may be involved in gene transcription and cell cycle processes that are essential for malignant progression.

To determine CDK6-AS1 expression, prognosis, and correlation with its target gene at the pan-cancer level, we analyzed 33 tumor cohorts from the TCGA database. CDK6-AS1 displayed higher expression in tumor tissues than in adjacent normal tissues in 16 tumor types: COAD, DLBC, ESCA, GBM, HNSC, KIRP, LAML, LIHC, LUSC, OV, PAAD, PCPG, READ, SARC, THYM, and UCS. Meanwhile, six types showed higher expression in normal tissues than in malignant ones: ACC, BRCA, KICH, PCPG, and TGCT (Figure 5A). In terms of survival analysis, high expression of CDK6-AS1 was related to poor prognosis in BLCA, HNSC, KIRC, LGG, LUAD, MESO, and THCA (Figures 5B–H), while showing better prognosis in UVM (Figure 5I). Further analysis was performed to determine the correlation between CDK6-AS1 expression and its target gene, CDK6. We found that CDK6-AS1 was correlated with its target in 29 tumor types (Supplementary Table S6). Taken together, when CDK6-AS1 is dysregulated it could influence the prognosis in a range of different cancer types. Representative immunohistochemical staining was used for gastric cancer tumor-infiltrating target gene CDK6. Scale bar, 50 mm (Figure 5J). Western blot results showed that CDK6 was significantly expressed in metastatic gastric cancer cells MKN-45 and in the in situ gastric cancer cell line HGC-27, while GES-1 was not significantly expressed in normal gastric mucosa epithelial cells (Figure 5K).

In order to identify predictive genes and construct a prognostic model, six genes related to DEG and OS co-expressed with CDK6-AS1 were crossed to obtain 6 differentially expressed genes related to CDK6-AS1 and OS. Then the six genes were analyzed by LASSOCox regression and the tuning parameter lambda (λ) was selected by using the cross-validation error curve. The prognostic models were constructed when the λ value was minimum (Figure 6A) (Risk score = 0.2637*CTHRC1 + 0.0132*PFN2 + 0.1384*PRSS35 + 0.0355*RTN4 - 0.072*SMPD3 + 0.5459*SYCP2L, and their LASSO coefficient curves are shown in Figure 6B. The relationship between survival status/risk score, mRNA expression heat map of 6 genes, and survival time (days)/risk score showed that the prognostic model had a good prognostic effect, and the OS of gastric cancer patients in the high-risk group was worse than that in the low-risk group (p < 0.001) (Figures 6C–F).

TIL are important players in the TME and have been reported to influence gastric cancer survival rates. Therefore, the association between CDK6-AS1 expression and the proportion of infiltrating lymphocytes was analyzed. Higher CDK6-AS1 expression showed negative tendencies with almost all inferred immune cell enrichment scores (Figures 7A–I) and was among Th_elper (Spearman r = −0.20, e < 0.001, Figure 7C), T_reg (Spearman r = −0.10, e = 0.044, Figure 7D), and neutrophil cell clusters (Spearman’s r = −0.14, p = 0.008, Figure 7I). Taken together, it appears that high-CDK6-AS1 expression may hamper immune TME cell infiltration, both in innate and adaptive cell clusters, which may have an effect on anti-tumor immunity.

Since chemo-sensitivity or -resistance is related to the clinical prognosis of gastric cancer, we explored the chemosensitivity of the high- and low-CDK6-AS1 expression groups. The ridge regression model was used to predict individual drug sensitivities. A commonly used chemotherapy drug in gastric cancer therapy, cisplatin, showed greater sensitivity in the high-CDK6-AS1 expression group than in the low-expression group (p = 0.024, Figure 8A). Conversely, paclitaxel showed higher sensitivity in the low-expression group (p = 0.040). We also explored the expression of immunotherapy and targeted therapeutic markers. PD-L1 expression was higher in the low-CDK6-AS1 expression group, which may indicate anti-PD-L1 therapy (Figure 8B). We can also prove that CDK6-AS1 is related to immunotherapy by TMB analysis, and the results show that the immunotherapy effect is better in the group with low expression of CDK6-AS1 (Figure 8C).