Differentially expression analysis
A total of 4989 up-regulated genes and 2349 down-regulated genes were revealed between tumor group and normal group. Among them, there were 4412 up-regulated and 1805 down-regulated mRNAs, as well as 155 up-regulated and 316 down-regulated lncRNAs between two groups (Supplementary Figure 2). The verification results based on the GEPIA2 database show that there were totally 2454 down-regulated genes and 5213 up-regulated genes. Moreover, there were 2229 intersected down-regulated genes (94.89% of GEPIA2 analysis results) and 4778 intersected up-regulated genes (95.77% of GEPIA2 analysis results), which indicating a reliable results of current differentially expression analysis.
Tumor immune score and cytolytic activity analysis
The Box plots were used to reveal the tumor immune score and cytolytic activity between tumor samples and normal samples. The result showed that the stromal score (Supplementary Figure 3A), immune score (Supplementary Figure 3B) and ESTIMATE score (Supplementary Figure 3C) in GBM group were significantly higher (all P < 0.05) than those in normal group. Moreover, the cytolytic activity score in GBM group was significantly higher (P < 0.05) than that in normal group (Supplementary Figure 3D).
Immune cell infiltration analysis based on CIBERSORT and TIMER
In the expression matrix of CIBERSORT, a total of 335 genes were included in the 547 genes of LM22 signal matrix (accounting for 61.4%). A total of 22 immune cell infiltration abundance matrices were estimated by CIBERSORT algorithm. In all 276 samples, only 66 samples were valid (P < 0.05), including 65 GBM samples and 1 normal sample (Figure 1A). Among these 66 samples, the infiltration ratio of M2 macrophages cells in GBM samples reached an average of 0.458 and a maximum of 0.759. The finding was confirmed in TCIA database that macrophages M2 cells had the highest average infiltration based on 604 samples (Figure 1B). However, in the unique normal sample (GTEX-ZYY3-3126-SM-5SI9), the infiltration ratio of M2 Macrophages was only 0.109, which was much lower than the GBM sample. Due to the small sample size, the result could not be verified in this data. Furthermore, the clustering heat map of 22 immune cell infiltration abundances was showed in Supplementary Figure 4. This result showed that the infiltration degree of macrophages M2 cells in GBM samples was significantly higher than that of other immune cells.
Based on TIMER, the immunocyte infiltration histograms of GBM vs. normal showed that the infiltration degree of Macrophage cells in GBM samples were significantly higher (P < 0.05) than those in normal samples (Figure 1C). Meanwhile, the Violin chart of infiltration abundance showed that except for B cells (P = 0.489) and CD4 T cells (P = 0.006), the infiltration abundance between two groups in all other cells were significantly different (all P < 0.01). The heatmap for GBM immune cell subsets was showed in Supplementary Figure 5. Furthermore, the ROC analysis of infiltration abundance of 6 immune cells and OS showed that Macrophage cell infiltration abundance was significantly related to sample survival time (Supplementary Figure 6). The lower the infiltration abundance, the longer the overall survival of the patient. Based on the results of the two algorithms, Macrophages, which had different abundance in GBM vs. normal and associated with the prognosis of GBM, were selected as the focus for the further investigation.
Immune cell associated DEMs/DElncRNAs investigation and enrichment analysis
Correlation analysis between Macrophages M2 infiltration abundance and gene expression was performed to explore the immune cell associated DEMs and DElncRNAs. The result showed that there were 168 immune-related DEMs and 5 immune-related DElncRNAs (including RP11-161H23.5, RP11-713C5.1, AC007620.3, RP11-74E22.8 and RP11-126K1.6) were revealed as immune cell associated genes by screening with P < 0.015.
The GO function analysis showed that the immune-related DEMs were mainly assembled in functions like Cognition (GO:0050890; Gene: BRSK1), Signal release (GO:0099643; Gene: FOX4 etc.) and Regulation of neurotransmitter transport (GO:0051590; Gene: BRSK1) (Figure 2A). Furthermore, the KEGG pathway enrichment analysis showed that the immune cell associated DEMs were mainly enriched in pathways like Endocrine resistance (hsa01522), Viral protein interaction with cytokine and cytokine receptor (hsa04061), as well as Chemokine signaling pathway (hsa04062) (Figure 2B).
CeRNA network construction
Person correlation analysis of lncRNA-mRNA interaction revealed that totally 137 mRNAs were positively correlated with 5 lncRNAs with r > 0.5 and P < 0.01. The top-10 mRNA-lncRNA interactions according to the correlation coefficient were showed in Supplementary Table 1. The result showed that all the 10 mRNAs including RUNDC3A, BRSK1 and SPTBN2 were interacted with the same lncRNA RP11-74E22.8. Moreover, there were totally 41 mRNAs, 126 miRNAs and 140 interactions in current mRNA-miRNA relations. Furthermore, the lncRNA-miRNA interaction analysis revealed a total of 580 interactions, 5 lncRNAs and 501 miRNAs.
A total of 37 miRNAs, 25 mRNAs and 5 lncRNAs (Supplementary Table 2 and Supplementary Table 3) were included in the current ceRNA interaction network (Figure 3A). The results showed that several lncRNA-miRNA-mRNA interactions such as RP11-161H23.5-miR6849-3p-KIF18B and RP11-74E22.8-miR-339-5p-BSK1 were involved in current ceRNA network. Moreover, the PPI network constructed by 25 mRNAs in current ceRNA network was showed in Figure 3B.
Risk prediction model construction based on prognostic analysis
Based on the COX univariate regression analysis, the mRNAs including SOX4, KIFC1, ELFN2, PLP2 and MEX3A, as well as lncRNA PR11-161H23.5 were related to the prognosis of PFS. Meanwhile, the K-M survival curve analysis showed that SOX4, KIFC1, ELFN2, MEX3A, and RP11-161H23.5 were significantly correlated with prognosis. However, the result of multivariate COX regression analysis showed that only RP11-161H23.5 was associated with prognosis (Supplementary Figure 7). Then, the risk prediction model based on the expression of SOX4, KIFC1, ELFN2, PLP2 and MEX3A and PR11-161H23.5was used to predict the clinical end point of PFS. The result showed that the AUC of the prediction model was 0.828 with a specificity of 0.758 and a sensitivity of 0.909, which could be used to predict the clinical end point of PFS (Figure 4).
According to the median value of risk score (0.08325018), 165 samples were divided into high risk (H) groups and low risk (L) groups. The result of samples sorted by risk score was showed in the Figure 5A. Meanwhile, the scatter plot of sample grouping and PFS survival time showed a shorter survival time in H group (Figure 5B). Moreover, the K-M curve analysis showed that the PFS survival time of H group was significantly lower than that of L group (Figure 5C). Furthermore, the clustering heat map of 6 gene expression levels among the sample groups showed that the expression levels of SOX4, KIFC1, ELFN2, and MEX3A were lower in the L group when compared to the H group. Meanwhile, the expression of PLP2 and RP11-161H23.5 were higher in samples of L group when compared to the H group (Figure 6). The box plot analysis for 6 gene expression showed that there was significant difference in the expression of each key gene between two groups (Figure 7).
Prognosis verification based on external data sets
Since there was no PFS survival information in the external data sets (GEO and CGGA datasets), the effect of the risk model constructed in the previous step could not be verified. Thus, only the differences of 6 key genes and OS prognosis were enrolled for current verification analysis. A total of 693 samples were downloaded from the CGGA database. After screening, there were 249 GBM and rGBM samples. The analysis results showed that KIFC1, ELFN2 and PLP2 were significantly related to OS (Supplementary Figure 8). Furthermore, COX survival analysis was performed based on 4 GEO validation datasets. The results showed that except for RP11-161H23.5, the COX survival analysis results of the other four genes were contribute to the survival (Supplementary Figure 9).