RNA-seq data of 161 GBM patients were obtained from The Cancer Genome Atlas (TCGA) datasets. All patients had accurate follow-up data. For further analysis, RNA transcriptome profiling was then performed using log2-based transformation of FPKM values by R Studio software (version 1.2.5001). Figure 1 shows the paths of the data analysis.

PreSTIGE, an algorithm that could predict tissue-specific enhancers based on the H3K4me1 mark and tissue-specific expression of mRNAs, was applied to predict lncRNAs expressed from active tissue-specific enhancers and their target genes (Corradin et al., 2014; Vucicevic et al., 2015). Ensembl transcript ID and gene symbol were annotated by “biomaRt” R package (human GRCh38.p13) (Durinck et al., 2009; Smedley et al., 2009). The log-rank test was used to calculate the differences between eRNA expression (cutoff value = 50%) and survival time in GBM patients using a threshold of p < 0.05. Based on the result of log-rank test, Kaplan–Meier curves were plotted to indicate the survival time differences. According to the list of eRNAs and target genes, Spearman correlation test was used to evaluate correlations and Spearman rank correlation coefficient r > 0.4, p < 0.001 as the threshold.

All transcripts with significant correlation (Spearman rank correlation coefficient r > 0.4, p < 0.001) with prognosis-related eRNAs were identified as eRNA-related genes by Spearman correlation test. To investigate gene ontology and signaling pathway enrichment, eRNA-related genes were uploaded to online tool Metascape¹ website for gene annotation, functional enrichment, and interactome analysis (Ashburner et al., 2000; Lebrec et al., 2009; Gene Ontology Consortium, 2015; Zhou et al., 2019). The terms were considered as significant and grouped into clusters when p < 0.01 and the numbers of enriched genes ≥3 based on their membership similarities. According to results of functional enrichment analysis, we confirmed the immune-related eRNA for further research.

One hundred sixty-one GBM patients were divided into two groups, 80 AC003092.1_H patients and 81 AC003092.1_L patients, according to the expression of AC003092.1 (cutoff value = 50%). In every GBM patient, the single-sample gene-set enrichment analysis (ssGSEA) score was used to measure the enrichment levels of the 29 immune signatures (Supplementary Table 3) (Barbie et al., 2009; Hanzelmann et al., 2013). The infiltration of immune cell (immune score), stromal content (stromal score), and tumor purity were quantified by Estimation of STromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) (Yoshihara et al., 2013). The immunogenomic analysis was executed by R Studio software (version 1.2.5001).

The CIBERSORTx and the LM22 gene signature were used to estimate the fraction of immune cell in the GBM patients (Newman et al., 2019). CIBERSORTx² is an analytical tool from the Alizadeh Lab and Newman Lab to impute gene expression profiles and provide an estimation of the abundances of member cell types in a mixed cell population, using gene expression data. We uploaded RNA-seq data as the Mixture file and choose LM22 (22 immune cell types) as Signature matrix file. The successful deconvolution of a sample should be satisfied with the criteria of 1,000 permutations and p < 0.05. Mann–Whitney U test was used to compare the proportions of the immune cell phenotypes between AC003092.1_H and AC003092.1_L.

To confirm the RNA-sequencing data, we validated AC003092.1 and TFPI2 transcript levels in glioma samples and GBM cell lines by reverse transcription quantitative polymerase chain reaction (RT-qPCR). The isolated total RNA was sourced in 11 HGG patients and 5 LGG patients who were diagnosed at our institution in Sun Yat-sen University Cancer Center (SYSUCC) from February 1, 2019, to July 1, 2019 (Supplementary Table 4), and 7 glioblastoma cell lines (U251, U138, T98G, SKMG4, U87, U118, and LN229). U251, U138, T98G, SKMG4, U87, U118, and LN229 were obtained from the State Key Laboratory of Oncology in South China and were maintained in Dulbecco modified eagle medium supplemented with 10% fetal bovine serum. All cells were maintained in a humidified incubator at 37°C and 5% CO₂. The Ethics Committee of Sun Yat-sen University Cancer Center approved the procedures (no. GZR2018-244). Total RNA was isolated from tissue samples or cultured cells using RNA-Quick Purification Kit (ESscience, Shanghai, China). The cDNA was synthesized from 1,000 ng total RNA by Fast All-in-One RT Kit (ESscience, Shanghai, China). RT-qPCR was performed by Bio-Rad CFX96 Real-Time PCR System (Bio-Rad Laboratories, Inc., Hercules, CA, United States) and LightCycler 480 II Real-Time PCR System (Roche, Basel, Switzerland) with PowerUp SYBR Green Master Mix (Thermo Fisher Scientific, Waltham, MA, United States). Custom primers for AC003092.1 (forward: TTAGCAGCAAACCCAGAAC, reve rse: TGCTGAGGATACATGACGAA) and TFPI2 (forward: GCCTGAGAACTTTGAATGATGCTG, reverse: GGCCCTGT GTTTCTTATGTATCCTG) were obtained from Tianyi Huiyuan Bioscience and Technology Inc., Wuhan, China. Each sample was measured in triplicate, and the results were normalized to the housekeeping gene U6 (forward: CTCGCTTCGGCAGCACA, reverse: AACGCTTCACGAATTTGCGT) and GAPDH (for ward: GGAGCGAGATCCCTCCAAAAT, reverse: GGCTGTT GTCATACTTCTCATGG).

According to previous study, there are 2,695 lncRNA transcripts, and 2,303 predicted target genes have been identified by the PreSTIGE algorithm (Vucicevic et al., 2015). The eRNA-target genes pairs were confirmed by this eRNA transcripts database. For further analysis, we conversed transcript ID to gene ID by “biomaRt” R package (Durinck et al., 2009). After conversion, the 2,695 eRNA transcripts were mapped to 1,288 genes. Then log-rank test was used to analyze the relationship of eRNAs and over survival in RNA-sequence data from 161 GBM samples of TCGA database. Next, we identified that 70 eRNAs were significantly associated with overall survival (Supplementary Table 1, log-rank test, p < 0.05).

Finally, the correlation between eRNAs and their predicted target genes was evaluated, and only 28 eRNAs were included in the next study (Spearman rank correlation coefficient r > 0.4, p < 0.001; Table 1).

To further evaluate biological process that eRNAs may regulate, we uploaded the eRNA-related genes to the website Metascape to identify Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. GO terms and KEGG pathways functional enrichment analysis showed that 420 genes that were significantly correlated to lncRNA AC003092.1 enriched notably in immune-related GO terms (Supplementary Table 2 and Figure 2A). The terms were considered as significant and grouped into clusters when p < 0.01 and the numbers of enriched genes were ≥ 3 based on their membership similarities. Each node represents an enriched term and is first colored by its cluster ID (Figure 2B) and then by its p value (Figure 2C). In glioblastoma patients, the AC003092.1 showed a positive correlation with the gene expression of its predicted target gene TFPI2 (tissue factor pathway inhibitor 2) (Figure 2D, Spearman r = 0.6, p < 0.001). A previous study had shown that AC003092.1 could upregulate the expression of TFPI2 by inhibiting the function of miR-195 (Xu et al., 2018). According to the expression of AC003092.1, the patients were divided into two groups: AC003092.1 High (AC003092.1_H) and AC003092.1 Low (AC003092.1_L) (cutoff value = 50%). Survival analyses showed that the AC003092.1_L had a better prognosis than the AC003092.1_H (Figure 2E, log-rank test, p = 0.006).

To confirm the immune-related feature between AC003092.1_H and AC003092.1_L, we analyzed 29 sets of immune-associated genes (Supplementary Table 3) representing different immune cell types, functions, and pathways in each GBM patient by ssGSEA analysis (Barbie et al., 2009; Hanzelmann et al., 2013). The results of ssGSEA indicated that immune cell types, functions, and pathways were enriched in AC003092.1_H (Figure 3A). Although the AC003092.1_H had a stronger immune function, the overall survival is better in AC003092.1_L. So, we further confirm the expression of immune checkpoints. We found that PD-L1, ICOS, and TIM-3 expressed highly in AC003092.1_H [analysis of variance (ANOVA) test, p < 0.05] (Figure 3B). Therefore, CD8⁺ T cells in AC003092.1_H were exhausted, and AC003092.1_H patients might get more benefit from immune checkpoint inhibitors.

Moreover, the immune scores, stromal scores, ESTIMATE scores, and tumor purity were compared between AC003092.1_H and AC003092.1_L by ESTIMATE. The results showed that the immune scores (Kruskal–Wallis test, p < 0.001), ESTIMATE scores (Kruskal–Wallis test, p < 0.001), and stromal scores (Kruskal–Wallis test, p < 0.001) were higher in the AC003092.1_H group. However, the tumor purity was higher in the AC003092.1_L group (Kruskal–Wallis test, p < 0.001) (Figure 3C). From the above, there were more immune and stromal cells in AC003092.1_H, whereas there were more tumor cells in AC003092.1_L.

In order to analyze the difference of immune cell proportion between AC003092.1_H and AC003092.1_L, we applied CIBERSORTx, which is a deconvolution algorithm based on gene expression (Newman et al., 2019). CIBERSORTx was combined with LM22, which have 547 gene expression characteristics to distinguish 22 immune cell subsets and define p value of the proportions of immune cell subsets. One of the GBM patients in AC003092.1_L was excluded because of p > 0.05. According to the result of CIBERSORTx, the proportion of immune cells between AC003092.1_H and AC003092.1_L varies widely (Figure 4A). Monocytes and neutrophils had a higher proportion in the AC003092.1_H, whereas the percentages of naive B cell, T cells, follicular helper cells, and macrophages M1 were higher in AC003092.1_L (Mann–Whitney U test, p < 0.05) (Figure 4B).

The correlation between AC003092.1 and TFPI2 in 11 patients with high-grade gliomas, 5 patients with low-grade gliomas, and 7 glioblastoma cell lines was investigated using RT-qPCR. In high-grade gliomas and glioblastoma cell lines, there was a significantly positive correlation between AC003092.1 and TFPI2 (high-grade gliomas: Spearman rank correlation coefficient r = 0.8614, p < 0.001; glioblastoma cell lines: Spearman rank correlation coefficient r = 0.8200, p < 0.05) (Figures 5A,B and Supplementary Figure 1). There was no significant correlation between the two RNAs in low-grade gliomas, the same as the result in 529 low-grade gliomas of TCGA (Spearman rank correlation coefficient r = 0.8646, p > 0.05) (Figures 5C,D).