PMT scores were positively associated with the risky behaviors of gliomas
To explore the general pattern of PMT status in glioma, we compared the distribution of PMT scores according to WHO grade, histology, and molecular subtype. We found that PMT scores were increased with increased risky progression of gliomas. Cases with higher grades from TCGA RNAseq dataset exhibited higher PMT scores (Fig. 1a). According to histopathologic classification, cases with higher PMT scores were enriched in order of GBM, astrocytoma, oligoastrocytoma, and oligodendroglioma. Regarding the five molecular subtypes described in our previous study [17], PMT scores were higher in the early-GBM (LGG-IDH-wt) and GBM with wild-type IDH (GBM-IDHwt, Fig. 1a). Similar PMT score distributions were verified in TCGA microarray, CGGA RNAseq, and microarray, Rembrandt, and GSE16011 datasets (Fig. 1b, c, Additional file 2: Fig. S1A–D).
According to TCGA microarray gene expression profiles, the gliomas were classified into four subtypes [3]. The MES subtype had the riskiest progression and most unfavorable survival prognosis. Here, we found that PMT score was significantly higher in the MES subtype than in the PN subtype in multiple public glioma datasets (Fig. 1a–c, Additional file 2: Fig. S1B–D). Meanwhile, glioma cases, even GBM, with higher PMT scores had poorer survival prognoses (Fig. 1d, e, Additional file 2: Fig. S1E–M). Furthermore, there was no difference in survival between the cases with low and high PMT scores without radiotherapy or chemotherapy. However, cases with higher PMT scores that received radiotherapy or chemotherapy still suffered reduced survival time (Fig. 1f–i, Additional file 2: Fig. S2). These results suggested that the PMT balance plays a critical role in the risky progression of glioma.
Identification of PMT related lncRNAs by WGCNA analysis
To determine the PMT-related lncRNAs with potential biological functions, we conducted WGCNA analysis to generate lncRNA co-expression modules. After removing outlier cases (Additional file 2: Fig. S3A and B), 3308 DELs were divided into different modules by WGCNA cluster analysis in the training set (TCGA RNAseq set, Fig. 2a). Here, we chose β = 4 as the soft threshold to build a scale-free network (Additional file 2: Fig. S3C), and then obtained 12 lncRNA co-expression modules for further analyses (Fig. 2a). All non-co-expression lncRNAs were gathered into the grey module (Fig. 2a).
To assess correlations between each module and PMT score, overall lncRNA expression in the respective modules was estimated by module signature (MS). Then correlation analysis between MS and clinical characteristics was performed. Co-expression modules with an absolute R value of > 0.6 were selected for further study (Fig. 2b). We found that the turquoise module was best correlated with PMT score, with an R value close to 0.9. Additionally, we found that the brown, yellow, and black modules were highly negatively related to PMT scores. All four above-mentioned modules were tightly associated with high-risk clinical and molecular features, such as WHO grade, Histopathology, molecular classification, TCGA subtype, Chr 7 gain/Chr 10 loss status, IDH status, and 1p/19q status. Thus, the turquoise, brown, yellow, and black modules were identified as modules of interest for further analyses.
Next, we screened core lncRNAs associated with PMT by setting the GS threshold to > 0.5 and the MM threshold to > 0.7 (Fig. 2c). This analysis identified 149 lncRNAs that were closely related to PMT (Fig. 2d and Additional file 1: Table S2). As expected, a strong relationship between several of these lncRNAs and mesenchymal transition in solid tumors had been reported in the literature. For example, mesenchymal transition was reinforced by CRNDE via the Wnt/β-catenin pathway in osteosarcoma [18], and lncRNA PVT1 sponged miR-195 to enhance EMT and induce therapeutic resistance in cervical cancer [19]. Therefore, we identified 149 lncRNAs that might play important roles in regulating mesenchymal transition.
Establishing a PMT-related 10-lncRNA risk signature
To further explore the prognostic value of these PMT-related core lncRNAs in glioma, we performed univariate Cox analyses on these lncRNAs in all glioma, non-GBM, and GBM cohorts in TCGA RNAseq set. This yielded 17 lncRNAs with significant prognostic value (P < 0.05) in the three cohorts that were selected to generate the lncRNA risk signature (Fig. 3a, b). The lncRNAs were of two types: i.e., protective and risky; 14 lncRNAs with HR > 1 were defined as risky, whereas three lncRNAs with HR < 1 were defined as protective.
To predict the prognosis of glioma cases by PMT-related lncRNAs, we used TCGA RNAseq set as the training dataset to perform LASSO regression analysis on these 17 lncRNAs. Finally, 10 of the lncRNAs were selected to construct a risk signature (Fig. 3c). Meanwhile, we used coefficients obtained from LASSO analysis to calculate risk scores for each of the TCGA and CGGA glioma cases (Fig. 3d). To analyze the prognostic value of the PMT-related 10-lncRNA risk signature, glioma cases in TCGA and CGGA RNAseq sets were dichotomized based on the median risk score. There was a significant difference in survival between the low- and high-risk groups (Fig. 3e, j). Stratified analyses of TCGA and CGGA datasets showed that higher risk scores were associated with unfavorable prognosis in all WHO grades (Fig. 3f–i, k–n). In summary, we established a PMT-related lncRNA signature with robust prognostic value in glioma.
The risk signature was closely related to clinical and molecular features in glioma
The heatmaps in Fig. 4a and Additional file 2: Fig. S4A showed the distribution of clinical and molecular features according to risk score. There was a significantly different distribution in age, Karnofsky performance score (KPS), WHO grade, histopathology, molecular subtype, TCGA subtype, Chr 7 gain/Chr 10 loss status, 1p/19q status, and related key molecular events (IDH, MGMT promoter, TERT promoter, PTEN, ATRX, TP53, and EGFR) between the high- and low-risk groups (Fig. 4a, Additional file 1: Table S3, Additional file 2: Fig. S4A). We also analyzed the correlation between risk score and clinical, pathological, and molecular characteristics. Cases with risky clinical factors and genomic events (MGMT promoter unmethylated, wild-type IDH, and 1p/19q non-codeletion) had higher risk scores (Fig. 4b, Additional file 2: Fig. S4B–G). Receiver operating characteristic curve (ROC) analyses indicated that risk signature might effectively predict the MES subtype in glioma. The areas under the corresponding ROC curves in TCGA and CGGA were 96.67% and 94.32%, respectively, which were both more than the areas of LINC00152 and LOXL1-AS1, and which were previously reported to be associated with PMT in glioma (Fig. 4c, Additional file 2: Fig. S4H). Furthermore, GSEA confirmed that glioma cases with high risk scores were more likely to be the MES subtype, whereas cases with low risk scores were more often the PN subtype (Fig. 4d, Additional file 2: Fig. S4I and J), indicating that the lncRNA risk score could discriminate PMT status. Additionally, we found that three (TP73-AS1, LINC01503, and CRNDE) of ten PMT-related lncRNAs highly expressed both in TCGA (Additional file 2: Fig. S9A) and CGGA (Additional file 2: Fig. S9B). TP73-AS1 and CRNDE were reported to promote malignant progression of glioma [20,21,22,23]. However, researches about the impact of LINC01503 on glioma were rare. Here, we selected LINC01503 as an example to validate the impact of lncRNAs on regulating PMT in glioma cell. In TCGA and CGGA, we found the expression level of LINC01503 was positively correlated to the mesenchymal marker (CD44; Additional file 2: Fig. S9C), but negatively associated with the pro-neural marker (OLIG2; Additional file 2: Fig. S9C). In U87 glioma cell, suppression of LICN01503 could significantly inhibit the expression of CD44, but increase OLIG2 expression, suggestive of mesenchymal-to-proneural transition (Additional file 2: Fig. S9D–F). Together, these findings showed the lncRNA risk signature was positively correlated with risky behaviors in glioma.
Prognostic value of the risk signature in different stratified cohorts of high-grade glioma cases
The prognosis of high-grade (Grade III/IV) glioma cases is poor. To accurately evaluate the prognosis of these cases, we analyzed the prognostic value of the PMT-related 10-lncRNA risk signature in different stratified cohorts of high-grade glioma cases. In TCGA, we found that the risk score had significant predictive value in cohorts stratified by age and KPS (Fig. 5a–d). There was also a significant difference in survival time between the high- and low-risk groups regardless of whether the patients received radiotherapy or chemotherapy (Fig. 5e–h).
Next, the prognostic value of the risk score in different molecular stratified cohorts was assessed. Under different statuses, such as IDH, MGMT promoter, ATRX, and TP53, cases with high risk scores had decreased survival time compared with cases with low risk scores (Fig. 5i–p). High scores indicated reduced OS in cases with wild-type EGFR and PTEN, but not in cases with mutant EGFR and PTEN (Fig. 5q–t). Similar results were validated in the CGGA dataset (Additional file 2: Fig. S5). Additionally, multivariate Cox analysis suggested that the risk signature was an independent prognostic factor in high-grade glioma (Fig. 6). In summary, the risk signature could significantly predict the prognosis of high-grade glioma.
Biological functions of lncRNA in the risk signature
To explore the different functional features of GBM cases between the high- and low-risk groups, we first screened the coding genes that were positively correlated with the risk score (R > 0.3, p < 0.05) in GBM cases in the training dataset (TCGA, Additional file 1: Table S4). We then used the ClueGO plugin of Cytoscape to functionally annotate these coding genes. Then, we found that the most relevant biological processes in high-risk group were cell activation involved in immune response, cell adhesion and collagen-containing extracellular matrix, aminoglycan metabolic process, and positive regulation of tyrosine phosphorylation of STAT protein (Fig. 7a, b). These functional enrichments in the high-risk group were further verified by GSEA analyses (Fig. 7c). Analyses of the CGGA dataset revealed similar results (Additional file 1: Table S4; Additional file 2: Fig. S6). Additionally, previous studies have shown that STAT3 activation was positively associated with irradiation-induced PMT [24]. Taken together, we inferred that the 10 lncRNAs of the signature might facilitate biological processes that regulated PMT in GBM.
GBM cases with high risk scores exhibited enhanced immunosuppression
Recently, the MES subtype was reported to recruit more immunosuppressive cells than the PN and CL subtypes [25]. LncRNAs are important regulators of PMT, as described above; however, the correlation between PMT-related lncRNAs and immune responses in glioma are unknown. Therefore, we further analyzed immune responses related to the PMT-related lncRNA risk signature in GBM by calculating the scores of different immune cell subsets [26] by ssGSEA in the GSVA package, and immune score, stromal score, and tumor purity by the ESIMATE package [27] in GBM. Pearson’s correlation analysis was performed to evaluate the relationship between these tumor microenvironment indicators and risk score. We found that risk score was positively correlated with immune score and stromal score, but negatively correlated with tumor purity (Fig. 8a). Meanwhile, enrichment of eosinophils, macrophages, neutrophils, and regulatory T cells was positively associated with risk score, while central memory T cells and helper T cells were inversely related to risk score (Fig. 8b, c, Additional file 2: Fig. S7A and B). Further analyses showed that the PMT-related lncRNA risk score was positively related to the expression levels of the immunosuppressive genes described previously [28] (markers of Tregs, immunosuppressive signaling pathways, tumor-supportive macrophage chemotactic and skewing molecules, and immune suppressors; Fig. 8d, Additional file 2: Fig. S7C and D). Next, we intersected RSGs from TCGA and CGGA with immunosuppressive genes (Additional file 2: Fig. S7E). Finally, eight effectors were identified to be significantly positively correlated with the risk signature (Fig. 8e, Additional file 2: Fig. S7F). Together, these findings indicated that the ten lncRNAs could reinforce PMT to enrich certain immune cell subsets in the tumor microenvironment.
Constructing the PMT-related lncRNA/miRNA/mRNA network
Previous studies have shown that lncRNAs can act as ceRNAs to sponge miRNAs and release the coding mRNA to enhance its biological function. To explore whether any of the 10 PMT-related lncRNAs acted as ceRNAs, we predicted the potential miRNAs that could bind these lncRNAs using the DIANA website and intersected them with the downregulated miRNAs in GBM. Then, 15 miRNAs were selected for subsequent analyses. The potential target genes of these miRNAs were predicted by the starBase v3.0 database. Next, these target genes were intersected with the coding genes positively related to risk score in both TCGA and CGGA. In total, 82 potential target genes were identified (Additional file 2: Fig. S8A).
Subsequently, we used Cytoscape to build an initial ceRNA network based on the filtered lncRNAs, miRNAs, and mRNAs (Additional file 2: Fig. S8B). Meanwhile, we constructed a protein interaction network to analyze relationships among the 82 target genes. Then, nine core genes were screened by the MCODE method, with k = 2 and degree = 3 (Additional file 2: Fig. S8C). We further generated a core lncRNA/miRNA/hub gene network that included six lncRNAs, 10 miRNAs, and nine core coding mRNAs (Additional file 2: Fig. S8D). Further GO analyses by ClueGO showed that the functions of these nine core genes were involved in regulating leukocyte aggregation, and positive regulation of adhesion organization and regulation of apoptosis (Additional file 2: Fig. S8D). Moreover, some of the ceRNAs in the network had previously been reported to affect risky behaviors and immune responses. For example, miR-512-3p binded to CD44, resulting in suppression of invasion and cell adhesion in breast carcinoma [29]. Additionally, the mesenchymal marker CD44 was identified to significantly positively regulate PD-L1 expression in breast cancer and non-small cell lung cancer [30]. Overall, these results suggested that several PMT-related lncRNAs could promote mesenchymal transition and elicit immunosuppression via this ceRNA network in GBM.
