LGG patients with clinical and gene expression files were downloaded from a public database The Cancer Genome Atlas (TCGA; http://cancergenome.nih.gov/) by using the R package TCGA-Assembler (15). LGG patients with survival information were procured from the UCSC Cancer Genomics Browser (https://xenabrowser.net/datapages/) (16). We used OS and PFS as endpoints and removed those without radiotherapy information (n = 29) and survival information (n = 3). Hypoxia-related genes were extracted from GeneCards (https://www.genecards.org/). Genes with a common symbol name in TCGA were selected. The flowchart is summarized in Figure 1 . Finally, we obtained 466 patients with 5,403 hypoxia-related genes for analysis. Gene expression and clinical profiles of 443 LGG patients (CGGA693 dataset) (17, 18) and 182 LGG patients (CGGA325 dataset) (19, 20) were downloaded as external validation datasets from the Chinese Glioma Genome Atlas (CGGA) dataset (http://www.cgga.org.cn/). The RNA-seq transcriptome data were estimated as log2(x + 1) transformed. The cleaned clinical data are summarized in Supplementary Tables 1 , 2 .

WGCNA can identify highly related genes in thousands of genes and cluster them into modules and then was used to establish the relationship between phenotypic traits and gene expression data. By calculating the correlation degree between the gene module and the external clinicopathological information, we can obtain the module genes highly related to the clinicopathological information and obtain the hub genes. WGCNA can be implemented by R package WGCNA (21).

In our study, WGCNA was performed to discover radiotherapy-related genes. We analyzed hypoxia-related genes and clinical data, including OS status, PFS status, age, grade, radiotherapy, and treatment response. First, hierarchical clustering analysis was utilized to exclude the outliers. Subsequently, the “pickSoftThreshold” function was performed to estimate the value of the powers. The R-squared criterion was set to 0.9. Pearson’s correlation matrices were used for all pairs of genes, and the weighted adjacency matrix was constructed using the power function. After the power was selected, the adjacency matrix was converted to a topological overlap matrix (TOM). Genes with similar expression profiles were classified into gene modules, and hierarchical clustering was performed by the class average method based on TOM. The minimum gene size in each module was set as 30. To further analyze the modules, the dissimilarity of module eigengenes was calculated, and some modules were merged. The merged cutoff threshold was set to 0.2, which meant that modules with a similarity higher than 0.8 were merged into one module. Then, the correlations between modules and clinical factors of LGG were investigated by using Pearson’s correlation test. Finally, the genes of the most significant radiotherapy-related module were chosen for subsequent analysis.

To obtain the function of genes in the radiotherapy-related module, we performed the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses by using the R clusterprofiler (22). The GO analysis included biological processes (BPs), cellular component (CC), and molecular function (MF).

In our study, radiosensitivity for the patients was defined in terms of survival benefit ( Supplementary Figure 1 ). 1) In patients who received radiotherapy, patients in group A had a better survival rate than the patients in group B. Then patients in group A could be defined as radiosensitive patients (RS group). 2) In patients who did not receive radiotherapy, the survival rate of group A (RS group) was not better (equal or worse) than that of the other group.

The radiosensitivity prediction model was constructed in the patients who had received radiotherapy. Patients who received radiotherapy provided more information related to radiosensitivity. We defined a radiosensitivity prediction model for the patients satisfying both of the following criteria: the constructed radiation sensitivity correlation model constructed can be used to divide the population into a high-risk group and a low-risk group. 1) In the radiotherapy patients, the survival rate of the low-risk group was higher than that of the high-risk group. 2) There was no significant difference in survival between the high- and low-risk groups in the group that did not receive radiation therapy. The low-risk group was defined as the RS group, and the high-risk group was defined as the radioresistant (RR) group. Therefore, the radiosensitivity model can select the patients with better benefit from radiotherapy, while the results in the population without radiotherapy can better show that the model is related to radiosensitivity.

A log-rank test was applied to assess the relationship between the expression of genes in radiotherapy-related modules and the OS of radiotherapy patients in TCGA. Whole LGG patients were divided into the high- and low-expression level groups using the median gene expression level as a cutoff point. The RNAs with log−rank p < 0.05 in the radiotherapy patients and log-rank p > 0.05 in the non-radiotherapy patients were identified as radiosensitivity−related RNAs. Then, the least absolute shrinkage and selection operator (Lasso) regression was performed to narrow the range of genes in patients who received radiotherapy. The radiosensitivity-related risk score was computed as follows:

R a d i o s e n s i t i v i t y − r e l a t e d   r i s k   s c o r e   =   Σ ( β R N A n   ×   e x p r R N A n )

The LGG patients were divided into the high- and low-risk groups with the median radiosensitivity-related risk score as the cutoff. The Kaplan–Meier method was used to plot survival curves. Time-dependent receiver operating characteristic (ROC) curve analysis was used to evaluate the prognostic value. The radiosensitivity-related risk score was validated in TCGA and two CGGA datasets.

Univariate and multivariate Cox proportional hazard regression analyses were used to examine whether the radiosensitivity-related risk score was an independent prognostic factor. The forest plot was plotted to show the hazard ratio (HR) and 95% CIs.

To evaluate the 1-, 3-, and 5-year survival probability for patients with LGG, a nomogram model including all independent prognostic factors was built for LGG patients in TCGA. The nomogram model was validated with the PFS of TCGA and two CGGA datasets.

All statistical analyses were performed using R software (4.0.2). WGCNA was performed by using the “WGCNA” R package. Lasso analysis was conducted by using the “glmnet” R package. A bilateral p-value <0.05 was considered statistically significant.

WGCNA was performed in TCGA-LGG dataset to determine the coexpression network most highly associated with the radiotherapy modules. The hclust function was used to determine if there were any outliers ( Supplementary Figure 2 ). A total of 466 samples were in the clusters after removing 8 outliers in the samples based on the average linkage method. When the soft threshold power value was β = 7 and the scale R² = 0.84, the average connectivity of the RNA group was high, and the connectivity between genes conformed to the scale-free network distribution ( Figures 2A, B ). The scale-free topological fitting index R-square was calculated to reach 0.84 ( Figures 2C, D ). Next, the TOM was constructed ( Figure 3A ), and a topological overlapping heatmap was depicted of the TOM including the top 400 genes ( Figure 3B ). A total of 13 modules were identified from the RNA coexpression network after merging modules with a similarity higher than 0.8. The relationships between gene modules and clinical traits are shown in a heatmap ( Figure 4A ). Thus, the black module was considered to have the highest correlation with radiotherapy (r = 0.69, p < 0.001) and was considered a “radiotherapy-related module” ( Figure 4B ).

GO analysis was performed to analyze the function of the radiotherapy-related module ( Figure 4C ). We discovered that the radiotherapy-related module was functionally associated with responding to oxygen levels, responding to decreased oxygen levels, and responding to hypoxia, cardiac chamber morphogenesis, cardiac chamber development, regulation of mRNA stability, and regulation of RNA stability. CCs include the transcription repressor complex, nuclear speck, transcription regulation complex, and methyltransferase complex.

We selected modules related to radiotherapy for further analysis. A total of 231 genes were subjected to the log-rank test in radiotherapy patients and non-radiotherapy patients. Thirty-six radiosensitivity−related genes were identified in the univariate analysis. Subsequently, the Lasso Cox regression model was used to identify the most robust markers for prognosis ( Figures 5A, B ). Finally, 12 genes (AGK, ETV4, PARD6A, PTP4A2, RIOK3, SIGMAR1, SLC34A2, SMURF1, STK33, TCEAL1, TFPI, and UROS) were included in the model. The radiosensitivity-related risk scores were calculated based on the linear combination of the expression levels of genes multiplied by the corresponding Lasso coefficients. The radiosensitivity-related risk score was as follows: radiosensitivity-related risk score = 0.16864 * AGK + 0.14242 * ETV4 − 0.14386 * PARD6A + 0.00584 * PTP4A2 − 0.09746 * RIOK3 − 0.23212 * SIGMAR1 + 0.07849 * SLC34A2 + 0.18813 * SMURF1 + 0.04491 * STK33 − 0.4922 * TCEAL1 + 0.01536 * TFPI − 0.4694 * UROS.

Then, the patients with LGG in TCGA dataset were divided into the high-risk (n = 233) or low-risk groups (n = 233) according to the median risk score. The Kaplan–Meier analysis revealed that OS time was significantly increased in the low-risk group compared with the high-risk group in patients who received radiotherapy (p < 0.001, Figure 5C ). There was no difference in OS between the high-risk group and low-risk group in patients who did not receive radiotherapy (p = 0.54, Figure 5D ). The low-risk group was defined as an RS group, and the high-risk group was defined as an RR group. The risk score distribution of each patient in TCGA is shown in Figure 5E .

Then, ROC analysis was used to evaluate the predictive efficiency of the radiosensitivity-related risk score model in the 1-, 3-, and 5-year survival rates (1-year area under the curve (AUC): 0.935 (0.904–0.967); 3-year AUC: 0.856 (0.778–0.933); 5-year AUC: 0.787 (0.704–0.87), Figure 5F ).

A radiosensitivity model was validated in TCGA with PFS as the endpoint. The Kaplan–Meier plots indicated that patients in the RR group exhibited worse PFS than patients in the RS group in patients who received radiotherapy (p < 0.001, Figure 6B ). Time-dependent ROC analysis results showed that the AUCs of the radiosensitivity model were 0.74, 0.676, and 0.732 at survival times of 1, 3, and 5 years, respectively ( Figure 6C ). Plots of risk score distribution are shown in Figures 6A, D, G .

Patients in the CGGA693 and CGGA325 datasets were divided into the RS group and RR group based on the median risk score in each dataset. The Kaplan–Meier analysis showed that patients in the RS group had a more favorable outcome than patients in the RR group in patients who received radiotherapy (CGGA693, p < 0.001; CGGA325, p < 0.001; Figures 6E, H ). These results indicated the accuracy of the radiosensitivity-related signature in predicting the outcomes of LGG patients. ROC curves were used to evaluate the predictive accuracy for 1-, 3-, and 5-year survival. AUC values revealed the high predictive value of the radiosensitivity-related risk score for LGG patients (CGGA693: 1-year AUC: 0.636 (0.527–0.746); 3-year AUC: 0.655 (0.587–0.732); 5-year AUC: 0.643 (0.572–0.713); CGGA325: 1-year AUC: 0.696 (0.567–0.862); 3-year AUC: 0.745 (0.644–0.845); 5-year AUC: 0.731 (0.639–0.823), Figures 6F, I ).

Then, univariate and multivariable Cox regression analyses were conducted to evaluate whether the radiosensitivity-related risk score is an independent prognostic factor for LGG. The results indicated that factors such as risk score and grade were significantly correlated with patient survival in both TCGA dataset and two CGGA datasets. Age (HR: 1.054, 95% CI: 1.038–1.071, p < 0.001), tumor grade (HR: 2.715, 95% CI: 1.736–4.247, p < 0.001), and risk score (HR: 2.712, 95% CI: 1.763–4.171, p < 0.001) were significantly associated with OS. The univariate analysis indicated that a high-risk score was significantly correlated with poor OS. The multivariate Cox regression results showed that the radiosensitivity-related risk score was an independent prognostic factor for LGG patients after adjusting for clinical factors such as age, sex, tumor grade, race, and IDH1. When OS was used as an endpoint, the HR was 2.029 (95% CI: 1.407–3.448, p < 0.001, Figure 7A ). When PFS was used as an endpoint, the HR was 2.170 (HR: 2.170, 95% CI: 1.526–3.086, p < 0.001, Figure 7B ). In CGGA datasets, we adjusted for clinical factors such as age, sex, tumor grade, race, IDH2, and X1p19q2, and the multivariate Cox regression results also demonstrated that the radiosensitivity-related risk score was an independent prognostic factor for LGG (CGGA693: HR: 1.730, 95% CI: 1.215–2.463, p = 0.003; Figure 7C ). Unfortunately, the multivariate Cox regression result was not significant in the CGGA325 dataset (CGGA325: HR: 1.609, 95% CI: 0.902–2.871, p = 0.112; Figure 7D ). We consider that there are too few patients in the CGGA325 database.

Age, tumor grade, and risk score were listed as candidate indicators for nomogram construction. Then, an optimal nomogram was established combining age, tumor grade, and risk score to predict a certain clinical outcome ( Figure 8A ). Figure 8B shows that AUCs of the nomogram for 1-, 3-, and 5-year OS were 0.947 (0.915–0.978), 0.888 (0.83–0.946), and 0.850 (0.779–0.922), respectively, which were better than those of the models with a single risk score model. Figure 8C demonstrates that AUCs of the nomogram for 1-, 3-, and 5-year PFS were 0.74 (0.665–0.815), 0.676 (0.601–0.750), and 0.732 (0.638–0.826), respectively. We also used two CGGA datasets to verify a nomogram model. Figure 8D demonstrates that the AUCs of nomogram at 1, 3, and 5 years were 0.64 (95% CI: 0.542–0.739), 0.669 (95% CI: 0.602–0.736), and 0.64 (95% CI: 0.569–0.711), respectively, for CGGA693. Figure 8E demonstrates that the AUCs of nomogram at 1, 3, and 5 years were 0.740 (95% CI: 0.608–0.872), 0.787 (95% CI: 0.694–0.88), and 0.79 (95% CI: 0.707–0.872), respectively, for CGGA325. The nomogram model for the prognostic model displayed superior predictive performance as compared with the risk score in the CGGA325.