3.1 Consistent clustering of glioma and the survival rate of clusters
The study was conducted as described in the flow chart (Fig. 1). Consistent clustering analysis of glioma patients was performed based on the expression of glycolysis genes in the TCGA database. Finally, patients were decided to be divided into 2 sample clusters, with cluster 1 containing 132 samples and cluster 2 containing 512 samples (Fig. 2A). The K-M suggested a significant difference in cumulative mortality at 20 years between the two subtypes (P < 0.0001; Fig. 2B). In the glioma cohort, the best outcome was found in patients classified as cluster 2. The next question was to determine the relationship between sample clustering results and clinicopathological data (treatment, age, gender, IDH, grade, and MGMT). Matching these results to all clinicopathological characteristics, significant correlations had been noted among treatment, age, IDH, grade, MGMT and the two clusters (Table 2). In summary, these results appeared to imply that glioma could be stratified according to the expression pattern of glycolysis genes.
Table 2
The relationship between sample clusters and clinical characteristics of glioma patients.
|
|
cluster1
|
cluster2
|
P
|
treatment
|
no PRT
|
16
|
108
|
3.04E-05
|
|
RT
|
7
|
67
|
|
|
PT
|
6
|
48
|
|
|
PT + RT
|
95
|
247
|
|
age
|
> 60
|
56
|
80
|
7.33E-11
|
|
< 60
|
68
|
390
|
|
gender
|
male
|
78
|
262
|
0.1831
|
|
female
|
46
|
208
|
|
IDH
|
IDH-mut
|
15
|
321
|
< 2.2e-16
|
|
IDH-wt
|
78
|
89
|
|
Grade
|
G2
|
4
|
186
|
< 2.2e-16
|
|
G3
|
19
|
195
|
|
|
G4
|
70
|
29
|
|
MGMT
|
Methylated
|
43
|
334
|
3.76E-12
|
|
Unmethylated
|
50
|
76
|
|
no PRT: no Pharmacotherapy and radiotherapy; PT: Pharmacotherapy; RT: radiotherapy.
Furthermore, we also assessed differences in pathway enrichment scores between the two clusters based on the KEGG pathway using the ssGSEA algorithm (Supplementary Table 6). Unexpectedly, glycolytic (‘Glycolysis / Gluconeogenesis’) and glycan synthesis and metabolic (‘Amino sugar and nucleotide sugar metabolism’, ‘Glycosaminoglycan biosynthesis-keratan sulfate’, ‘Mannose type O-glycan biosynthesis’, ‘Galactose metabolism’, etc.) pathways were significantly different between the two clusters. Also, immune response- (‘cAMP signaling pathway’, ‘IL-17 signaling pathway’, ‘Cell adhesion molecules’, ‘Intestinal immune network for IgA production’, etc.) and immune cell (‘Th17 cell differentiation’, ‘Th1 and Th2 cell differentiation’, etc.)-related pathways were notable between Cluster 1 and Cluster 2. Further, PD-L1 expression and PD-1 checkpoint pathway in cancer also showed substantial variations.
3.2 Establishment of a prognostic signature based on the DE-glycolysis-related lncRNAs
We identified a total of 2 DE-lncRNAs from both subtypes, relative to cluster 2, lncRNA FLJ16779 was the up-regulated gene and lncRNA AL390755.1 was the down-regulated gene (Fig. 3A; Supplementary Table 7). These genes were defined as DE-glycolysis-related lncRNAs for the subsequent analysis. To verify whether these 2 DE-lncRNAs were the risk factors for glioma, we performed a univariate Cox regression analysis in the training set. Coincidentally, all DE-lncRNAs met P < 0.05 (Fig. 3B). Subsequently, a multifactorial Cox regression analysis with a STEP function indicated these 2 lncRNAs as the optimal variables for constructing a prognostic signature (Fig. 3C). Specifically, AL390755.1 (P = 8.59E-14) was possible a risk factor for glioma (HR = 1.38), whereas FLJ16779 (P = 1.82E-6) was inferred to be the oncogene (HR = 0.77).
3.3 Assessment and validation of the effectiveness of the prognostic feature
Risk scores were calculated separately for each sample in the training and testing sets according to the aforementioned formula, and patients were categorized into high- and low-risk groups based on the median risk score of each set. In both the training and testing sets, with increasing risk scores, the number of patient deaths climbed sharply (Figs. 4A and D). Subsequently, K-M-survival analysis was performed on this prognostic feature, while AUC values were calculated by time-dependent ROC analysis for the accuracy of which in predicting survival, using the outcome variable. In the training set, K-M curves could effectively distinguish between high-risk and low-risk groups (Fig. 4B), and the risk scores all reached an AUC of 0.8 or more at 1 to 5 years (Fig. 4C). Similar results were also reproduced in the testing set (Figs. 4E and F). Based on the expression heatmap of the prognostic lncRNAs, which was observed in the high-risk group associated with poorer OS, AL390755.1 was overexpressed. In contrast, FLJ16779 was overexpressed in the low-risk group possessing a longer OS characteristic. Moreover, we found that the vast majority of Cluster 2 patients were in the low-risk group. These results suggested that the prognostic signature was highly sensitive and specific and that these signature lncRNAs could be utilized as prognostic biomarkers in clinically. Furthermore, the statistical tables of clinical information in the training and testing sets were displayed in Table 3.
Table 3. clinical information of glioma patients in the training and testing sets.
|
Training set
|
Testing sets
|
|
|
Expression
|
p
|
|
Expression
|
p
|
Total
(N=417)
|
High
(N=203)
|
Low
(N=214)
|
Total
(N=172)
|
High
(N=88)
|
Low
(N=84)
|
gender
|
|
|
|
0.832
|
|
|
|
0.687
|
female
|
184
(44.1%)
|
88
(43.3%)
|
96
(44.9%)
|
|
68
(39.5%)
|
33
(37.5%)
|
35
(41.7%)
|
|
male
|
233
(55.9%)
|
115
(56.7%)
|
118
(55.1%)
|
104
(60.5%)
|
55
(62.5%)
|
49
(58.3%)
|
|
age(years)
|
|
|
|
<0.001
|
|
|
|
<0.001
|
≥60
|
93
(22.3%)
|
74
(36.5%)
|
19
(8.9%)
|
|
42
(24.4%)
|
37
(42.0%)
|
5
(6.0%)
|
|
<60
|
324
(77.7%)
|
129
(63.5%)
|
195
(91.1%)
|
130
(75.6%)
|
51
(58.0%)
|
79
(94.0%)
|
|
cluster
|
|
|
|
<0.001
|
|
|
|
<0.001
|
1
|
86
(20.6%)
|
82
(40.4%)
|
4
(1.9%)
|
|
36
(20.9%)
|
34
(38.6%)
|
2
(2.4%)
|
|
2
|
331
(79.4%)
|
121
(59.6%)
|
210
(98.1%)
|
136
(79.1%)
|
54
(61.4%)
|
82
(97.6%)
|
|
type
|
|
|
|
<0.001
|
|
|
|
0.604
|
no PRT
|
85
(20.4%)
|
27
(13.3%)
|
58
(27.1%)
|
|
37
(21.5%)
|
16
(18.2%)
|
21
(25.0%)
|
|
RT
|
52
(12.5%)
|
14
(6.9%)
|
38
(17.8%)
|
22
(12.8%)
|
12
(13.6%)
|
10
(11.9%)
|
|
PT
|
38
(9.1%)
|
6
(3.0%)
|
32
(15.0%)
|
16
(9.3%)
|
7
(8.0%)
|
9
(10.7%)
|
|
PT+RT
|
242
(58.0%)
|
156
(76.8%)
|
86
(40.2%)
|
97
(56.4%)
|
53
(60.2%)
|
44
(52. 4%)
|
|
IDH
|
|
|
|
<0.001
|
|
|
|
<0.001
|
IDHmut-codel
|
114
(27.3%)
|
21
(10.3%)
|
93
(43.5%)
|
|
38
(22.1%)
|
7
(8.0%)
|
31
(36.9%)
|
|
IDHmut-non-codel
|
158
(37.9%)
|
43
(21.2%)
|
115
(53.7%)
|
|
68
(39.5%)
|
16
(18.2%)
|
52
(61.9%)
|
|
IDHwt
|
139
(33.3%)
|
134
(66.0%)
|
5
(2.3%)
|
|
63
(36.6%)
|
63
(71.6%)
|
0
(0%)
|
|
Missing
|
6
(1.4%)
|
5
(2.5%)
|
1
(0.5%)
|
|
3
(1.7%)
|
2
(2.3%)
|
1
(1.2%)
|
|
MGMT
|
|
|
|
<0.001
|
|
|
|
0.001
|
+
|
300
(71.9%)
|
101
(49.8%)
|
199
(93.0%)
|
|
123
(71.5%)
|
46
(52.3%)
|
77
(91.7%)
|
|
-
|
96
(23.0%)
|
81
(39.9%)
|
15
(7.0%)
|
|
40
(23.3%)
|
33
(37.5%)
|
7
(8.3%)
|
|
Missing
|
21
(5.0%)
|
21
(10.3%)
|
0
(0%)
|
|
9
(5.2%)
|
9
(10.2%)
|
0
(0%)
|
|
grade
|
|
|
|
<0.001
|
|
|
|
<0.001
|
G2
|
136
(32.6%)
|
25
(12.3%)
|
111
(51.9%)
|
|
55
(32.0%)
|
16
(18.2%)
|
39
(46.4%)
|
|
G3
|
148
(35.5%)
|
75
(36.9%)
|
73
(34.1%)
|
|
64
(37.2%)
|
29
(33.0%)
|
35
(41.7%)
|
|
G4
|
92
(22.1%)
|
92
(45.3%)
|
0
(0%)
|
|
41
(23.8%)
|
39
(44.3%)
|
2
(2.4%)
|
|
Missing
|
41
(9.8%)
|
11
(5.4%)
|
30
(14.0%)
|
|
12
(7%)
|
4
(4.5%)
|
8
(9.5%)
|
|
no PRT: no Pharmacotherapy and radiotherapy; PT: Pharmacotherapy; RT: radiotherapy; +: methylated, -: unmethylated;
3.4 Independent prognostic analysis
Here, risk scores and clinicopathological characteristics (age, gender, treatment type, cluster, MGMT, grade, and IDH) of the full Gliomas sample in the TCGA database were included in the Cox analysis to explore their potential for independent prognosis. Ultimately, age, IDH, MGMT, grade, and risk score were identified as independent prognostic factors for glioma through univariate and multivariate Cox analyses (Figs. 5A and B). These independent prognostic factors were then incorporated into the Nomogram to explore the prediction of the Nomogram model for patient survival at 1, 3, and 5 years, which with a c-index of 0.8705766 (Fig. 5C). The status of each variable corresponded to a score, with a higher total score for a patient indicating poorer survival for that patient. Besides, the calibration curves also demonstrated that the Nomogram model was effective in predicting the survival of patients (Fig. 5D).
3.5 Functional enrichment analysis of risk score-related DEGs
To further reveal the potential function of prognostic lncRNAs, we first screened 680 risk score-related DEGs between high- and low-risk groups (Supplementary Table 8) and performed GSEA on them to reveal potential functions. Figures 6A and B illustrated the top10 GO and KEGG results, respectively. GO analysis revealed that these DEGs were mainly enriched in terms related to immune responses (‘ADAPTIVE IMMUNE RESPONSE’, ‘ADAPTIVE IMMUNE RESPONSE BASED ON SOMATIC RECOMBINATION OF IMMUNE RECEPTORS BUILT FROM IMMUNOGLOBULIN SUPERFAMILY DOMAINS’, ‘HUMORAL IMMUNE RESPONSE’, ‘NEGATIVE REGULATION OF IMMUNE SYSTEM PROCESS’, etc.) and immune cell regulation (‘LEUKOCYTE PROLIFERATION’, ‘REGULATION OF LYMPHOCYTE ACTIVATION’, ‘REGULATION OF T CELL ACTIVATION’, etc.). Notably, a multitude of biological processes associated with neurons (‘NEURON PROJECTION TERMINUS’, ‘NEURON TO NEURON SYNAPSE’, ‘REGULATION OF NEURON PROJECTION DEVELOPMENT’, etc.), neurotransmitters (‘NEUROTRANSMITTER RECEPTOR ACTIVITY’, ‘NEUROTRANSMITTER SECRETION’, ‘NEUROTRANSMITTER TRANSPORT’, etc.), and synapses (‘POSTSYNAPTIC MEMBRANE’, ‘POSTSYNAPTIC NEUROTRANSMITTER RECEPTOR ACTIVITY’, ‘POSTSYNAPTIC SPECIALIZATION MEMBRANE’, ‘PRESYNAPTIC MEMBRANE’, etc.) were significantly enriched. KEGG results indicated that risk score-related DEGs were significantly associated with immune disorders such as ‘AUTOIMMUNE THYROID DISEASE’ and ‘SYSTEMIC LUPUS ERYTHEMATOSUS’. Consistently, immune response-related pathways (‘CYTOKINE CYTOKINE RECEPTOR INTERACTION’, ‘PRIMARY IMMUNODEFICIENCY’, ‘NATURAL KILLER CELL MEDIATED CYTOTOXICITY’, etc.) were also significantly enriched. Moreover, the ‘AMINO SUGAR AND NUCLEOTIDE SUGAR METABOLISM’ pathway was inevitably featured in the results. More detailed results of GO and KEGG analysis were reported in Supplementary Tables 9 and 10.
Moreover, we illustrated the network of interactions of 150 risk score-related DEGs with prognostic lncRNAs for |cor| ≥ 0.6 and P < 0.05 in Fig. 6C. Among them, 35 genes were common DEGs for both prognostic lncRNAs (Table 4). Elaborately, ARPP21, JPH3, NTNG2, KCNIP3, and ELFN2 were positively correlated with lncRNAs FLJ16779, but negatively associated with AL390755.1. The remaining 30 genes such as EMP3, TNFRSF12A, MCUB, and ANXA2 were positively associated with lncRNA AL390755.1, but negatively related to FLJ16779.
Table 4
List of common risk score-related DEGs for prognostic lncRNAs (|cor| ≥ 0.6 and P < 0.05).
Gene name
|
Correlation
|
FLJ16779
|
AL390755.1
|
ELFN2
|
0.696
|
-0.623
|
KCNIP3
|
0.677
|
-0.626
|
NTNG2
|
0.637
|
-0.645
|
JPH3
|
0.625
|
-0.614
|
ARPP21
|
0.614
|
-0.614
|
EMP3
|
-0.669
|
0.700
|
TNFRSF12A
|
-0.615
|
0.684
|
MCUB
|
-0.697
|
0.679
|
ANXA2
|
-0.630
|
0.670
|
ANXA1
|
-0.616
|
0.666
|
CD58
|
-0.636
|
0.663
|
RAB42
|
-0.626
|
0.600
|
FABP5
|
-0.633
|
0.601
|
MEOX2
|
-0.620
|
0.602
|
XKR8
|
-0.600
|
0.605
|
EFEMP2
|
-0.626
|
0.606
|
DPYD
|
-0.676
|
0.606
|
NAMPT
|
-0.632
|
0.610
|
PPIC
|
-0.630
|
0.610
|
TEAD3
|
-0.608
|
0.616
|
FAM114A1
|
-0.622
|
0.620
|
PINLYP
|
-0.636
|
0.620
|
SERPINH1
|
-0.609
|
0.621
|
RAB34
|
-0.669
|
0.623
|
OSMR
|
-0.620
|
0.628
|
SRPX2
|
-0.652
|
0.628
|
GDF15
|
-0.671
|
0.630
|
TNFAIP6
|
-0.650
|
0.631
|
S100A4
|
-0.634
|
0.631
|
STEAP3
|
-0.624
|
0.643
|
TUBA1C
|
-0.693
|
0.647
|
TIMP1
|
-0.694
|
0.651
|
PDPN
|
-0.638
|
0.655
|
SMIM3
|
-0.668
|
0.656
|
IGFBP2
|
-0.665
|
0.657
|
3.6 The impact of risk scores on the immune landscape of glioma patients
Inspired by the results of the ssGSEA, we hypothesized that prognostic lncRNAs may operate in the patient's immune microenvironment to influence patient outcomes. ESTIMATE analysis indicated that the immune, stromal, and ESTIMATE scores were significantly higher in the high-risk group than in the low-risk group, suggesting that the high-risk group had more components of the immune microenvironment (Fig. 7A). Subsequent immune cell enrichment analysis demonstrated that more immune cells, such as Natural killer T cells, Myeloid-derived suppressor cell (MDSC), and Type 1 T helper cell, were presented in the high-risk group (Fig. 7B). This evidence suggested that prognostic lncRNAs were involved in the altered immune microenvironment of patients.
3.7 The relationship between risk scores and patient response to ICI therapy
Over the past decade, ICIs have proven to be promising agents for many solid tumor malignancies, adapting this treatment strategy to an increasingly important role in Neuro-oncology (Zhao et al., 2019; Almquist et al., 2020). Therefore, we investigated the potential relationship between risk score and the expression of nine ICIs. Seven ICIs, excluded LAG3 and TIGIT were detected to have higher expression levels in the high-risk group, which may be related to the higher tumor grade of patients in which group (Fig. 7C). The TIDE algorithm then predicted the sensitivity of patients in the high- and low-risk groups to ICI treatment. TIDE scores were significantly higher in the high-risk group than in the low-risk group, suggesting that patients in the low-risk group may benefit from ICI treatment (Fig. 7D). Fig. 7. (A). Comparison of immune, stromal, and ESTIMATE scores between the high-risk group and low-risk group. (B). Volcano plots for the enrichment of immune cell types based on the normalized NES score between low-risk and high-risk patients from the GSEA. (C). The relationship between risk score and the expression of nine ICIs. (D). Comparison of TIDE scores between the high-risk group and the low-risk group.
3.8 Expression of AL390755.1 and FLJ16779 mRNA in NHA, U87 and A172 Cells.
By qRT-PCR validation of the expression of the two lncRNAs, we found that the expression of AL390755.1 was significantly lower in glioma cells than that in NHA cells, while the expression of FLJ16779 showed the opposite results (Fig. 8).