Overview of mutation profiles.
A total of 508 mutation profiles were downloaded from TCGA. The 'maftools' package was used to visualize the results based on the mutation data. The waterfall plot was used to display the mutation information including the name and frequency of mutant genes, different colors were used to represent different mutation types (Figure 2).
The most common mutation type was 'Missense Mutation' and SNP was the most frequent variant. C>T, C>A, and T>C ranked in the top three of the SNV class, the median variants per sample was 25. Moreover, the top ten mutated genes in LGG ranked with percentage included IDH1(77%), TP53(46%), ATRX (33%), CIC (20%), TTN (12%), FUBP1(9%), MUC16(6%), NOTCH1(7%), PIK3CA (7%), EGFR (6%) (Figure 3A-F).
And their gene symbols were shown in the gene cloud (Figure S1). Besides, a significant interaction of mutant genes (co-occurrence and mutual exclusion) was exhibited in the interaction plot, among which EGFR-IDH1, CIC-TP53 showed the most mutually exclusive relationships and ATRX-TP53 showed most co-occurrence relationships. (Figure 3G)
High-TMB is associated with unfavorable prognosis.
Clinical data of 513 LGG patients were included in this study. The clinicopathological parameters included age, sex, tumor grade, TMB level, and TMBPI (Table 1). Tables
Table 1. Parameters of 513 LGG patients from TCGA
Parameters
|
Number (%)
|
Vital status
|
|
Alive
Dead
|
388(75.6)
125(24.4)
|
Age
≤45
>45
Gender
Female
Male
Grade
Grade Ⅱ
Grade Ⅲ
Unknown
TMB
High level
Low level
Undetected
TMBPI
High level
Low level
Unknown
|
310(62.2)
203(37.8)
289(56.3)
284(43.7)
248(48.3)
264(51.5)
1(0.2)
241(47.0)
260(50.7)
2(2.3)
249(48.5)
250(48.7)
14(2.8)
|
Patients were divided into two groups, high-TMB (n=241) and low-TMB (n=260) levels using the median TMB as the cutoff value. The survival analysis showed that the high-TMB group was associated with a worse prognosis compared with the low-TMB group (P<0.001, Figure 4A).
Furthermore, box plots displayed that the TMB level was higher in older patients and grade III (P<0.001) patients (Figure 4B, C). However, no significant differences were observed in associations of TMB with gender (Figure 4D).
Comparison of gene expression profiles between low-TMB and high-TMB groups
A total of 394 DEGs were revealed after differential analysis (See supplementary material). Heatmap indicated that the expression level of genes in the high-TMB group (Figure 5A) was lower than in the low TMB group. GO enrichment analysis revealed that the DEGs were mainly associated with cell division and composition of the extracellular matrix (Figure 5C). Besides, KEGG analysis showed several enriched pathways (ECM-receptor interaction, human T-cell leukemia virus 1 infection, cell cycle, transcriptional misregulation in cancer, and P53 signaling pathway, etc.) (Figure 5D).
GSEA analysis revealed 26 pathways with the criterion of FDR<0.05. Some of the pathways were involved in mismatch repair, p53 signaling pathway, cell cycle, glutathione metabolism (Figure 5E). For details see supplementary materials for GO, KEGG, and GSEA analysis results.
Differences in infiltration level of immune cells in TME
The level of 6 types of immune cells in TME (B cells, CD4+ T cells, CD8+ T cells, Neutrophils, Macrophages, and Dendritic cells) was related to survival prognosis of LGG patients (Figure S2). Multivariable Cox proportional hazard model of TIMER model, which was constructed as Surv(LGG) ~ B_ cell + CD8_Tcell + CD4_Tcell + Macrophage + Neutrophil + Dendritic, showed that high macrophage infiltration was an independent risk factor (HR=1845.804, p = 0.000, Table 2).
Table 2. Multivariate Cox regression analysis of lower-grade glioma by TIMER database.
Cell types
|
coef
|
HR
|
95%CI_l
|
95%CI_u
|
P-value
|
sig
|
B cell
|
0.304
|
1.355
|
0.005
|
362.209
|
0.915
|
-
|
CD8+ Tcell
|
6.076
|
435.270
|
0.523
|
362042.429
|
0.077
|
·
|
CD4+ Tcell
|
-3.956
|
0.019
|
0.000
|
33.325
|
0.299
|
-
|
Macrophage
|
7.521
|
1845.804
|
44.941
|
75811.118
|
0.000
|
***
|
Neutrophil
|
-5.201
|
0.006
|
0.000
|
8.569
|
0.165
|
-
|
Dendritic
|
2.136
|
8.469
|
0.202
|
355.685
|
0.263
|
-
|
R square= 0.103 (max possible= 9.09e-01), Likelihood ratio test p= 4.94e-10, Wald test p= 1.19e-12 Score (log rank) test p= 1.4e-13. P-value Significant Codes: 0 ≤ *** < 0.001 ≤ ** < 0.01 ≤ * < 0.05 ≤. < 0.1
Besides, based on the CIBERSORT algorithm, we obtained an estimation of the abundances of 22 immune cells in LGG, and results were visualized with bar plot. (Figure 6A).
Moreover, a comparison of the proportion of different immune cells and TMB levels (high-TMB vs. low-TMB) showed an increase in the proportion of macrophages M0 and M2 (P=0.033, P=0.009), and a decrease in the proportion of monocytes (P=0.006, Figure 6B).
Identification of 7 hub immune genes.
The intersection was obtained between DEGs from high- and low-TMB groups and immune-related genes. A total of 24 genes were screened (|logFC|>1.5, FDR <0.05) for further analysis (Figure 5B). To identify hub immune genes, the 24 genes were further screened by Cox regression and Kaplan-Meier analysis. Out of the 24 genes, 7 were considered hub TMB-related immune genes. They included HOXB4, IGF2BP3, IL13RA2, NPNT, PLA2G2A, TNFAIP6, and WISP1. All p-values computed by Cox regression and Kaplan-Meier analysis were less than 0.05 (Figure 7A).
Hub immune genes influence the infiltration of immune cells in TME.
The SCNA module on the TIMER webserver was used to explore the relationship between the 7 hub immune genes and the infiltration of immune cells in TME. The results indicated that the two most common alterations were Arm-level Gain and Deep Deletion, and these had a significant effect on the level of immune cells in TME, especially PLA2G2A, IGF2BP3 and NPNT (Figure 7B). The rest were shown in Figure S3.
Establishment and assessment of TMBPI model for LGG.
To predict LGG patients’ prognosis, the TMBPI model was established based on the 7 hub immune genes. Their expression data and coefficients were calculated by multivariate Cox regression analysis and the model was constructed as follows: PI= (0.248935×IGF2BP3+0.045491×NPNT+0.073771× WISP1-0.245454 × HOXB4 + 0.109997×TNFAIP6+0.017986×IL13RA2 +0.025691×PLA2G2A). Further, 499 out of 513 cases from the TCGA database were divided into high-risk (n=249) and low-risk (n=250) groups based on the median of TMBPI as the cutoff value. Besides, 413 cases from the CGGA database were divided into high-risk (n=215) and low-risk (n=216) groups. Compared with the low-risk group, the high-risk group was associated with poor prognosis (P<0.001, Figure 8A-B). The ROC curves showed that the AUC at 1-, 3- and 5-years overall survival (OS) prediction were 0.84, 0.83, and 0.73, respectively for the TCGA database. Consistently, the AUC at 1-, 3- and 5-years were 0.75, 0.75, and 0.79, respectively for the CGGA database (Figure 8C-D).