Multiomics Profiling and Clustering of Low-Grade Gliomas Based on the Integrated Stress Status

Background Although the prognosis of low-grade glioma is better than that of glioblastoma, there are still some groups with poor prognosis. The integrated stress response contributes to the malignant progress of tumors. As there had limited research focused on the integrated stress status in LGG, it is urgent to profile and reclassify LGG based on the integrated stress response. Methods Information of glioma patients was obtained from the Chinese Glioma Genome Atlas, The Cancer Genome Atlas, and the GSE16011 cohorts. Statistical analyses were conducted using GraphPad Prism 8 and R language. Results We summarized and quantified four types of integrated stress responses. Relationships between these four types of stress states and the clinical characteristics were analyzed in low-grade glioma. We then reclassified the patients based on these four scores and found that cluster 2 had the worst prognosis, while cluster 1 had the best prognosis. We also established an accurate integrated stress response risk signature for predicting cluster 2. We found that immune response and suppressive immune cell components were more enriched in the high-risk group. We also profiled the genomic differences between the low- and high-risk groups, including the nonmissense mutation of driver genes and the copy number variations. Conclusion Low-grade glioma patients were divided into three clusters based on the integrated stress status, with cluster 2 exhibiting malignant transformation trends. The signature adequately reflected the traits of cluster 2, while a high risk score indicated a worse prognosis and an enriched inhibitory immune microenvironment.


Introduction
Glioma is the most common and deadly primary malignant tumor in the central nervous system [1]. Gliomas can be divided into high-grade glioma and low-grade glioma (LGG). We defined grade 4 glioma as glioblastoma(GBM) [2] and grade 2 and 3 glioma as LGG. GBM is highly aggressive, and even with standardized treatment, the median survival time of patients is only 14.6 months [3]. Compared with GBM, LGG is less aggressive and their progression is slightly slower [4,5]. However, increasing studies have reported that many LGG patients have a tendency toward malignant transformation, which can transform into higher grades, leading to adverse outcomes [6]. The outcomes of glioma patients are highly variable and heterogeneous [7], even in patients with the same tumor grade. Although isocitrate dehydrogenase 1 (IDH1) mutations and chromosome 1p/19q codele-tions have been identified as good prognostic markers for LGG, there is still an urgent need to develop a new classification system for LGG, which can identify LGG patients with a tendency toward malignant transformations.
Studies have shown that the integrated stress response (ISR) is related to the mechanisms and progressions of many complex diseases [8], including cancer, diabetes, and metabolic diseases. The ISR state refers to an evolutionarily conserved intracellular signal network, which helps cells, tissues, and organisms adapt to various stimuli from the microenvironment. The ISR includes protein homeostasis defects, nutritional deficiencies, viral infections, and oxidative stress. The ISR restores cell homeostasis by reediting gene expressions. Different types of stress responses are represented by four special kinases (PERK, GCN2, PKR, and HRI) [9]. The organs most affected in the ISR are the brain and pancreas [10]. Exosomes induced by endoplasmic reticulum stress promote the immune escape of breast cancer by regulating the expression of PD-L1 in macrophages [11]. The redox state affects T cell activation and subsequent supervision of redox-sensitive immune regulatory transcription factors such as NF-κB, NFAT, and AP-1, which are involved in the pathogenesis of inflammation-related diseases [12]. Multiple studies have indicated that the ISR is tightly embedded in the innate immune response of cells, and all four ISR kinases play a role in immunity and inflammation [13].
Most of the carcinogenic pathways found in human cancers cause various forms of protein synthesis disorders, while the ISR controls protein synthesis and protein stability. Three ISR kinases (PKR, PERK, and GCN2) have been correlated with cancer [14], but there has been minimal research describing the status of the ISR in LGG. There is therefore an urgent need to analyze its distribution status and clinical value in LGG, to reclassify LGG from the perspective of integrated stress and to identify subgroups of LGG patients with malignant tendencies.

Summary of Gene Sets.
Four gene sets were summarized to represent corresponding biological processes of the ISR. They were endoplasmic reticulum stress (unfolded protein response), viral infection, nutritional deficiency, and oxidative stress. The specific genes of each gene set are shown in Supplementary Table 1.

Classification of LGG Based on the Four Types of ISRs.
The four types of integrated stress states were quantified by the ssGSEA algorithm [15], and unsupervised clustering was performed to obtain the most reasonable classification according to these four factors. Survival curves of four clusters of samples were plotted using Prism 8 software (Graph-Pad, San Diego, CA, USA). Differential genes were screened out using the limma package between the best and the worst clusters in survival quality.

Development and Validation of a ISR-Related Signature
Using LASSO Regression Model. The "glmnet" R package was performed to filter the prognosis-related ISR genes by LASSO Cox regression analysis with a ten-fold crossvalidation. After identifying the significant genes and their regression coefficients (β), we calculated the risk score of each LGG patient by the formula as follows: risk score = ∑ N i=1 Exp i * β i .

Function Enrichment
Analyses. GSEA (http://www .broadinstitute.org/gsea/index.jsp) was conducted to determine whether the selected gene sets showed statistical differences between different groups. Gene ontology (GO) [16] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [17] pathway enrichment analyses were used to analyze the enriched biological processes. Enriched ontological terms and pathways with P < 0:05 were selected and presented in a bubble map using the R package.
2.6. Immune Response and Tumor Microenvironment-Related Analyses. The immune cell genes were downloaded from http://cibersort.stanford.edu/ [18], and the immune cell score of each sample was quantified by ssGSEA. Student's t -test was used to identify immune cells with significant differences between different groups. Stromal score, immune score, and glioma purity were calculated to evaluate nontumor cells within the microenvironment [19]. The correlations between ISR score and tumor purity and stromal score and immune score were calculated and exhibited. Principal component analysis (PCA) was used to profile the distribution patterns of different groups on the basis of the immune-related transcriptome expression matrix.

Statistical
Analysis. SPSS (SPSS, Chicago, IL, USA), Prism 8 (GraphPad), and R 3.6 (https://www. http://rproject.org/) software were used for statistical analyses. Student's t -test was used to identify differences in expressions, and Pearson's correlation was used to calculate correlations. A value of P < 0:05 was considered to indicate statistical significance ( * P < 0:05, * * P < 0:01, * * * P < 0:001, and * * * * P < 0:0001, as indicated in the figures and legends). Low and high-expression groups were classified according to the median expression. Survival analyses were conducted using Kaplan-Meier plotting, and the log-rank test was used to evaluate the differences between stratified groups.

General Description of the ISR State in
LGG patients. The ISR includes the ability to adapt to various stress states in the cell, including the endoplasmic reticulum stress state, nutrient deprivation, viral infection, and redox imbalance. Because few investigators have studied the relationship between ISR and the malignant progression of LGG, we systematically analyzed the clinical value of the ISR in LGGs. We first selected four types of representative gene sets related to the ISR from the official website of GSEA and used the ssGSEA algorithm to quantify the four responses of LGG patients in TCGA, to show the corresponding relationships between the four response scores and LGG clinical characteristics using heat map (Figure 1(a)). We found that endoplasmic reticulum stress, nutritional deprivation, and redox imbalance in the ISR correlated with age, IDH1, 1p19q, and other significant LGG survival prognostic markers. These three stress scores were significantly increased in LGG  3 BioMed Research International patients with relatively poor prognoses, but the difference in viral infection scores was not significant. Subsequently, we used the median score as a critical value to analyze the prognostic evaluation value of the four stress states using log-rank survival curves. The results showed that redox state imbalance and the nutritional deprivation state had the most significant impact on the survivals and prognoses of LGG patients, and the survival time of patients in the high scoring group was significantly shortened (Figures 1(f)-1(i)). Patients with the high risk score of the endoplasmic reticulum stress, nutritional deprivation, and redox imbalance were more enriched in the IDH wild (Figures 1(b)-1(d)). Patients with IDH mutations possessed significantly higher viral infection risk score than wild ones (Figure 1(e)).

Reclassification of
LGG Patients Based on the Four Types of Stress Response Scores. According to previous results, we performed unsupervised clustering of LGG patients based on the four stress response scores, and the k-means clustering results showed that the best distinguishing categories involved three categories (Figures 2(a) and 2(b)). The survival curves of the three types of LGG patients showed that the survival prognosis of cluster 2 was the worst, while the survival prognosis of cluster 1 was the best (Figure 2(g)). We further described different types of ISRs in the three cluster patients using an overlay diagram. Among the cluster 2 subgroups of patients with relatively poor survival prognoses, the endoplasmic reticulum stress state and redox imbalance state scores were relatively high (Figures 2(c) and 2(e)), while the virus infection status score was relatively low. The virus infection stress score was relatively high in the cluster 1 subgroup of patients with relatively good survival prognoses ( Figure 2(f)). To further characterize the underlying mechanism of the large differences in survival prognoses between the cluster 1 and cluster 2 subgroups, we screened the differential genes between the two groups of patients ( Figure 2(h)) and performed functional enrichment analysis. The enrichment results suggested that cluster 2 was enriched with several immune-related biological processes such as the inflammatory immune response and immune cell chemotaxis (Figures 2(i) and 2(j)). Analysis of the immune microenvironment indicated that inhibitory immune components such as macrophages and neutrophils were more enriched in cluster 2 ( Supplementary Fig. 1A).       Table 2). To obtain the best prognostic markers, the LASSO Cox regression analysis was conducted (Figure 3(a)). A total of 11 gene signatures were generated, and the risk scores of the signature based on the regression coefficients were finally calculated (Supplementary Table 4). The ISR risk score of cluster 2 was significantly higher than that of cluster 1 (Figure 3(h)). To estimate the effectiveness of the ISR risk score in predicting prognoses, we constructed a receiver operating characteristic curve (Figures 3(b) -3(d)) based on 1-year, 2-year, and 3-year survival times and then calculated the area under the curve. The results showed a good predicting efficiency in survival in the discovery cohort. Next, we used the median score as a critical value to analyze the prognostic evaluation value of the three cohorts using the log-rank survival curve. The results showed that among the three cohorts, TCGA-RNA sequencing cohort,    10 BioMed Research International CGGA-RNA sequencing cohort, and CGGA-microarray cohort, there were significant differences in survival (P < 0:0001). Univariate and multivariate Cox regression analyses were conducted to identify the independent prognostic factors. The results suggested that the ISR risk score could be used as an independent survival-related factor (Figures 3(i) and 3(j)). We obtained the same conclusion in the validation group of the CGGA cohort (Supplementary Table 3). We also established a nomogram to predict the 1-year, 2-year, and 3-year prognoses of LGG patients (Figure 3(k)) and used the calibration chart to evaluate the model (Figure 3(l)). It showed that the ISR risk score was effective in predicting the prognosis of LGG patients.

The High-Risk Group with an Enhanced Local Immune
Phenotype. Because cluster 2 was enriched with immunerelated biological processes and infiltrated with inhibitory immune cell components, we further characterized the immunophenotypic differences between the high-and low-ISR-risk groups. We first performed PCA based on immune-related genes, which was summarized from the terms of "Immune Response" and "Immune Process" [20,21]. The results showed that the low-and high-risk group

12
BioMed Research International patients were robustly distributed with different locations and patterns, suggesting that low-and high-ISR-risk groups exhibited a totally different immune status (Figures 4(a)  and 4(b)). Moreover, we performed GSEA based on a series of immune relevant terms, such as "leukocyte activation," "leukocyte-mediated immunity," and "leukocyte chemotaxis." All these terms were significantly enhanced within the high-risk group (Figures 4(c) -4(h)). Together, the results showed that the high-risk group was characterized by an enhanced local immune phenotype.

The High-Risk Group Was Infiltrated with an Inhibitory
Immune Microenvironment. To further characterize the differences of immune microenvironments between the low-and high-risk groups, we downloaded the gene sets data of immune cells from the CIBERSORT website and quantified the immune cell score of each sample using the ssGESA algorithm. The violin diagram shown in Figure 5(a) shows that M0 and M2 type macrophages were more abundant in the high-risk group. In addition, we found that the risk score was positively correlated with the immune score and stromal score, while negatively correlated with the tumor purity (Figures 5(c) and 5(g)). Moreover, the risk score was positively correlated with the expression of several classical inhibitory immune checkpoint-related genes (Figures 5(b) and 5(f)). These results suggested that the high risk score robustly indicated an enhanced local immune response as well as a suppressed immune microenvironment.       To investigate the differences between the high-and low-ISR-risk score groups at the genomic level, we analyzed somatic mutations and copy number variations (CNVs) from the TCGA RNA sequencing cohort. The overall pattern of classical somatic mutations in LGG stratified by the ISR risk score is shown in Figures 6(a) and 6(b)). We found 7 significant somatic mutation genes between the high-and low-risk groups (P < 0:05). In the high-risk score group, the somatic mutation profile showed that TTN and CIC had higher mutation frequencies in the high-risk score group.
To identify the specific CNV events of the high-and low-ISR-risk groups, we profiled the copy number variant events between these two groups. The results suggested that the high-risk group had amplification of partial fragments in chromosomes 4, 7, and 12 and partial deletion of fragments in chromosome 9 (Figures 6(c) and 6(d)). Although there was no significant difference in the copy number of important chromosomal fragments such as 1p and 19q between the two groups, amplified genes such as TTN in the high-ISR-risk group may be involved in promoting the malignant progression of LGG.

Discussion
LGG comprise a group of primary brain tumors caused by glial cells. The LGG patients have a better prognosis than GBM patients ,but they are still incurable, and most of them will relapse or even transform into high-grade gliomas. Previous studies have confirmed that the ISR is related to the mechanisms of a variety of cancers [22,23]. In addition, there are few studies on the relationship between the ISR and LGG, so there is an urgent need to quantitate and characterize the ISR state in LGG.
By quantitating the four types of stress responses of LGG patients, we found that endoplasmic reticulum stress, nutritional deprivation, and redox imbalance in the ISR were related to age, IDH1, 1p19q, and other LGG significant survival prognostic markers. The above three stress scores were significantly increased in the relatively poor prognoses cluster 2 LGG patients, although the viral infection stress score reflected a different situation. In the majority of human cancers, including breast cancer, pancreatic cancer, lung cancer, skin cancer, prostate cancer, brain cancer, and even liquid malignant tumors, enhanced endoplasmic reticulum stress responses have been found [24]. Studies have shown that continuous activation of the endoplasmic reticulum stress sensor could make malignant cells more tumorigenic, metastatic, and drug-resistant [25,26]. In addition, the endoplasmic reticulum stress response further hinders the development of protective anti-cancer immunity by regulating the function of myeloid cells in the tumor microenvironment [27]. The endoplasmic reticulum stress state can activate the unfolded protein response, which is necessary for the development of T cells, but the continuous increase of endoplasmic reticulum stress can lead to damages of T cell function, as well as decreased survival [28,29]. Oxidative stress is the main intracellular signal transducer that maintains autophagy [30]. Once tumors are formed, autophagy allows tumor cells to survive under stress conditions, thereby promoting tumor progression [31,32]. Regarding nutritional deprivation, there is evidence suggesting that cancer cells can adjust signal pathways to adapt to the new environment and continue to survive [33,34]. This characteristic of cancer cells is considered to be one of the prerequisites for tumor progression and chemotherapy resistance [35].
Based on these considerations, we know that the occurrence and development of LGG are affected by multiple factors, so a single aspect of stress factors cannot represent its occurrence or development. More angles and more dimensions should be considered in making an overall analysis. We therefore considered grouping the LGG patients in TCGA according to the four kinds of stress factors, with LGG patients divided into three clusters with the best characteristics. Using integrated bioinformatics analyses, we found that in cluster 2, which also exhibited the worst survival condition in these three clusters, the scores of oxidative stress and endoplasmic reticulum stress were significantly higher than those of the other two groups, while the viral infection stress score was significantly lower than those of the other two groups. In cluster 1 with the best prognosis, we found that the viral infection stress score was significantly higher than those of the other two groups, while the oxidative stress risk score was significantly lower than those of the other two groups. Importantly, the nutritional deprivation stress score in the cluster 3 group was significantly lower than that of cluster 1 groups, and the viral stress score was also very high (but with no significant difference compared with cluster 2). However, the prognosis of cluster 3 was significantly worse than that of cluster 1. This showed that redox imbalance stress may play an absolute role in cluster 3, suggesting that the occurrence and development of LGG might be regulated by multiple factors.
The classification that we established based on four kinds of ISR is relatively difficult to implement in specific clinical cases. To make our research results easier for clinical use, we established an integrated stress response-related risk score to represent the traits of cluster 2. Based on this score, the LGG patients in the other two databases were divided into the high-and low-risk groups, resulting in a significant difference (P < 0:0001), which was also a good independent factor for survival predictions. Using PCA of TCGA database, we found significant differences in immune functions between the high-risk and low-risk groups. Immunosuppressive cells such as M2, M0, regulatory T cells, and neutrophils were more enriched in the high-ISR-risk score group, which indicated that the reorganization of the immune microenvironment may have contributed to the poor survival of the high-risk group. Similarly, in the validation group of CGGA cohort, we found that immune infiltration and cell chemotaxis were also more enriched in high-risk patients (Supplementary Fig. 3A-B). We also found that the ISR risk score was positively correlated with the expressions of immune checkpoint genes, which further validated our hypothesis that high ISR risk was involved in immune suppression.
We analyzed the gene mutations and copy number variations between the high-risk and low-risk groups and found that the high-risk group had a higher frequency of CIC and TTN mutations. Some report that TTN was at the top ranking of mutated genes in multiple solid tumors, including the gastric adenocarcinoma, small cell lung cancer, and colorectal adenocarcinoma. TTN mutation could also be used to represent increased tumor mutation burden (TMB) and correlated with objective response to immune checkpoint blockade (ICB). [36,37]. Bettegowda et al.'s study indicated that there was a unique relationship between recurrent chromosomal aberrations and mutations in CIC and suggested that mutations in CIC were key events in the development of oligodendrogliomas [38].
While we profiled copy number variations, we found that the high-risk score group had an amplification of chromosomes 4, 7, and 12. We performed functional analyses of the corresponding gene set and found that it was enriched in angiogenesis and response to external stimuli (Supplementary Fig. 2H).

Conclusions
Our study provided new insights into the relationships between the integrated stress responses and clinical traits of LGG patients. Furthermore, we identified cluster 1 with the most potential for malignant transformation, which exhibited the worst prognosis and inhibitory immune microenvironment. More importantly, we designed an ISR signature to predict cluster 2, which may facilitate the classification and prediction of individual LGG patients, although the model needs to be further validated in prospective studies and multicenter clinical trials.