A Novel Necroptosis-Related lncRNA Signature for Predicting Prognosis and Immune Response of Glioma

Glioma is one of the most common intracranial malignancies that plagues people around the world. Despite current improvements in treatment, the prognosis of glioma is often unsatisfactory. Necroptosis is a form of programmed cell death. As research progresses, the role of necroptosis in tumors has gradually attracted the attention of researchers. And lncRNA is regarded as a critical role in the development of cancer. Therefore, this study is aimed at establishing a prognostic model based on necroptosis-associated lncRNAs to accurately assess the prognosis and immune response of patients with glioma. The RNA sequences of glioma patients and normal brain samples were downloaded from The Cancer Genome Atlas (TCGA) and GTEx databases, respectively. The coexpression analysis was performed to identify the necroptosis-related lncRNAs. Then, we utilized LASSO analysis following univariate Cox analysis to construct a prognostic model. Subsequently, we applied the Kaplan-Meier curve, time-dependent receiver operating characteristics (ROC), and univariate and multivariate Cox regression analyses to assess the effectiveness of this model. And the functional enrichment analyses and immune-related analyses were employed to investigate the potential biological functions. A validation set was obtained from the Chinese Glioma Genome Atlas (CGGA) database. And qRT-PCR was employed to further validate the expression levels of selected necroptosis-associated lncRNAs. Seven necroptosis-related lncRNAs (FAM13A-AS1, JMJD1C-AS1, LBX2-AS1, ZBTB20-AS4, HAR1A, SNHG14, and LINC00900) were determined to construct a prognostic model. The area under the ROC curve (AUC) was 0.871, 0.901, and 0.911 at 1, 2, and 3 years, respectively. The risk score was shown to be an important independent predictor in both univariate and multivariate Cox regression analyses. Through functional enrichment analyses, we found that the differentially expressed genes (DEGs) were mainly enriched in protein binding and signaling-related biological functions and immune-associated pathways. In conclusion, we established and validated a novel necroptosis-related lncRNA signature, which could accurately predict the overall survival of glioma patients and serve as potential therapeutic targets.


Introduction
Glioma is the most common type of primary intracranial malignancy, accounting for approximately 84% of all malignant brain tumors (1). Traditional treatments of glioma are mainly based on surgery, supplemented by radiotherapy and chemotherapy, but the outcomes of patients are often poor (2). The five-year survival rate of high-grade glioma (HGG) patients is less than 5% (3). Although the prognosis of patients with low-grade glioma (LGG) is better than that of glioblastoma (GBM), they often face the risk of recurrence and transition to HGG (4,5). Therefore, new treatment strategies to handle this enormous burden are urgently needed.
Necroptosis is a novel type of programmed cell death mediated mainly through the activation of the RIPK1-RIPK3-MLKL signaling pathway by the TNFR superfamily, which ultimately leads to plasma membrane rupture, lysis, and the release of inflammatory factors (6,7). Among the mechanisms of drug resistance in tumors, the evasion of and resistance to apoptosis are considered to be the most important causes and often lead to the failure of traditional chemotherapy regimens (8). In contrast, related studies have shown that some conventional proapoptotic chemotherapeutic agents can bypass the resistance mechanism and exert antitumor effects by inducing necroptosis (9,10). This plays an extremely significant role in the death of tumor cells. Moreover, some reports have suggested that necroptosis plays a key role in limiting cancer metastasis, for example, shikonin could reduce osteosarcoma metastasis through induction of necroptosis (11,12).
Long noncoding RNAs (lncRNAs) have extremely important biological functions and play critical roles in cell cycling, proliferation, immune response, transcription, and translation (13)(14)(15). There are increasing evidences that lncRNAs are closely related to tumor development. For instance, HOTAIR is overexpressed in various cancers and correlates with the metastasis and invasion of these cancers (16). MALAT1 has been determined to be an important potential therapeutic target in lung adenocarcinoma and, in addition, it has a role in tumor growth inhibition in patients with glioma (17,18). Therefore, developing new lncRNA-related therapeutic strategies and relevant tumor subtypes seem to be of great practical importance.
Our study is aimed at creating a necroptosis-associated lncRNA model to predict the prognosis of patients with glioma and to explore the immune status and tumor microenvironment of patients in different subtypes to provide new insights for individualized treatment of glioma.

Materials and Methods
2.1. Acquisition and Collation of Raw Data. The RNA sequences (FPKM value) and matching clinical data of 698 glioma patients (project: TCGA-GBM, TCGA-LGG; Disease Type: gliomas) were obtained from TCGA database (https:// portal.gdc.cancer.gov). The RNA expression data of 1152 normal brain samples were downloaded from the GTEx database (https://gtexportal.org/home/datasets). We combined TCGA and GTEx data and then performed the "nor-malizeBetweenArrays" algorithm of the "limma" package to eliminate the batch effects. And the genetic mutation files of GBM and LGG patients were downloaded from TCGA database. Clinical and transcriptome data of 1018 glioma samples were achieved from the CGGA database (http:// www.cgga.org.cn/) for external validation. We then excluded patients with unknown clinical features. Finally, the training and testing cohorts were constructed after integration and standardization.
2.2. Screening of Differentially Expressed Necroptosis-Related lncRNAs. We captured 67 necroptosis-related genes from previous studies and listed them in Table S1 (19)(20)(21). Firstly, we applied the "maftools" package to generate mutation waterfall plots for LGG and GBM separately. Then, we downloaded the human gene annotation file from the GENCODE website (https://www.gencodegenes .org). Through the "gene biotype" in the file, we can know which gene belongs to protein-coding RNA and which gene belongs to long noncoding RNA. We used a Perl script to extract the lncRNA matrix from the gene expression matrix of the training set based on the "gene biotype" in the annotation file. And finally, a total of 732 lncRNAs were screened out. Next, we performed correlation analysis to determine the necroptosis-related lncRNAs with the criterion of correlation coefficients > 0:6 and p < 0:001 (Table S2). And the "limma" package was utilized to identify the differentially expressed lncRNAs between normal and tumor tissues (FDR < 0:05 and jlog 2 FCj ≥ 1). Finally, we took the intersection of the two results to obtain the differentially expressed necroptosisassociated lncRNAs and then listed them in Table S3.

Construction and Validation of a Prognostic lncRNA
Signature. The univariate Cox analysis was conducted to identify the differentially expressed necroptosis-associated lncRNAs with prognostic significance, and the results were placed in Table S4. Then, we utilized the "glmnet" package to perform LASSO regression of these lncRNAs. For minimizing the deviations, we used the "cv.glmnet" function to choose the best lambda. And eventually, the prognostic signature was established. Each glioma patient was assigned an individual risk score, and the formula for the risk score was defined as follows (where n, exp ðlncRN A i Þ, and coef ðlncRNA i Þ represented the quantity of those screened lncRNAs, the relevant expression level of each candidate lncRNA, and the regression coefficient of each lncRNA, respectively): For predicting the accuracy of the prognostic model, the "timeROC" package was employed to do the ROC analysis. And we utilized the "ggplot2" and "Rtsne" packages to perform the principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis. Then, univariate and multivariate Cox regression analyses were applied to assess whether the risk score could be used as an independent predictor of prognosis. Finally, we used the "rms" package to integrate the risk score, gender, age, and tumor grade to draw a nomogram and applied calibration curves to assess whether the prediction results matched the actual ones.

GO and KEGG Enrichment
Analyses. The Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses were conducted using the Sangerbox tools, a free online platform for data analysis (http://www.sangerbox.com/tool). And the criterion was p < 0:05 and FDR < 0:1.

Genetic Mutation and Necroptosis-Associated lncRNA
Screening. The genetic mutation diagrams of necroptosisrelated genes in GBM and LGG were shown in Figures 1(a) and 1(b). We found that 36.92% (144/390) of GBM patients had genetic mutations in necroptosisassociated genes. The highest mutation frequency among these genes was in EGFR, followed by ATRX and IDH1. While, in LGG patients, mutations were found in 89.13% (451/506) of the samples, the highest mutation rate was in IDH1, followed by ATRX and EGFR. Missense mutation was the most common type of variant in both GBM and LGG samples, and C>T was the most common in the SNV class. And as shown in Figure 1(c), we identified 12  7 BioMed Research International differentially expressed necroptosis-associated lncRNAs between normal and tumor samples. We then presented a heatmap to demonstrate the differential expression between the two categories in Figure 1(d). As the heatmap showed, 6 lncRNAs were upregulated and 6 lncRNAs were downregulated in the tumor category.

Construction and Validation of the Prognostic lncRNA
Signature. Based on the univariate regression analysis, we found that all 12 differentially expressed necroptosisassociated lncRNAs were significantly related to prognostic value (p < 0:001) (Figure 2(a)). We then did LASSO regression analysis on these lncRNAs and finally determined 7 lncRNAs to construct the prognostic model (Figures 2(b) and 2(c)). According to the model, we calculated the risk score for each glioma patient and classified the patients into two subgroups by median value. And we plotted heatmaps of the model-related lncRNAs and clinical characteristics between the two risk categories in the training cohort and the CGGA cohort, respectively (Figures 2(e) and 2(f)). The distributions of age, grade, IDH mutation status, 1p19q codeletion status, and MGMT methylation status were found to diverge between the low-and high-risk categories. MGMT methylation, IDH mutation, and 1p19q codeletion were more frequently seen in patients in the low-risk group compared to the high-risk group. And patients of advanced age and high pathological grade were more often clustered in the high-risk group. In addition, we also created a Sankey    (Figure 2(d)). The next step was to verify the effectiveness of the model. As shown in Figures 3(a) and 3(b), patients were separated into two subgroups based on the median cutoff value. And we found that patients in the highand low-risk groups had dramatically different survival sta-tus (Figures 3(c) and 3(d)). The t-SNE and PCA analyses revealed that patients in both risk categories were distributed in different directions in the training cohort (Figures 3(e) and 3(f)) and validation cohort (Figures 3(g) and 3(h)). Depending on the Kaplan-Meier analysis, we discovered that      (Figure 3(l)). And we observed superior predictive ability in our risk model compared to clinical characteristics (Figures 3(m) and 3(n)).

Creating and Verifying a Nomogram.
Based on the Cox regression analysis, we constructed a nomogram to predict the prognosis of glioma patients (Figure 4(e)). The nomogram allowed us to score each patient and finally derive a total score to project the 1-, 2-, and 3-year OS occurrences. And as presented in Figure 4(f), we also plotted 1-, 2-, and 3-year calibration curves to assess the predictive value of the nomogram.

Functional Enrichment Analyses.
The key biological activities of DEGs linked with the risk model were deter-mined using GO and KEGG functional analyses. As shown in Figures 5(a)-5(c), the GO analysis suggested that the DEGs were significantly enriched in biological functions such as cell signaling and protein binding. Moreover, the KEGG pathway analysis indicated that the DEGs were remarkably concentrated in immune-related pathways, including phagosome, antigen processing, and presentation (p < 0:05, Figure 5(d)).

Exploration of the Immune Responses in Different Risk
Subgroups. We performed the ssGSEA to investigate the immune responses of different risk categories. As presented in Figure 6(a), we found that most of the immune cells had higher infiltration levels in the high-risk group, except for the Th1 cells and neutrophils which infiltrated more in the low-risk category. Furthermore, in the high-risk subgroup, all 13 immune pathways were shown to be more active (Figure 6(b)). We observed comparable findings in the validation set (Figures 6(c) and 6(d)). As shown in Figures 6(e) and 6(f), we explored the differences in immune checkpoints expression between the high-and low-risk categories, and most immune checkpoints were more positive in the high risk class. All of these reflected that the immune response was more active in the high-risk group of glioma patients, and we could select appropriate therapeutic targets to improve the efficacy of immunotherapy. Figure 7, patients in the high-risk subgroup had higher stromal and immune scores. This may imply that the higher potential for patients in the high-risk group to benefit from immunotherapy. More details can be found in Table S6. 3.8. Validating the Expression Levels of Selected Necroptosis-Associated lncRNAs. As shown in Figure 8, we found that these selected necroptosis-associated lncRNAs were -log10(pvalue) 8       14 BioMed Research International differentially expressed in tumor tissues and normal tissues, and the expression levels were also different between tumor tissues of different grades. These results further confirmed our previous conclusions.

Discussion
Immunotherapy first originated in the nineteenth century. In 1893, William Coley first treated cancer by activating the immune system. However, due to the immune evasion of tumor cells, which led to unsatisfactory treatment effect, the cancer immunotherapy did not attract much attention (22). In contrast, over the past few decades, immunotherapy has gradually become the first-line treatment option for many tumors as the tumor microenvironment and immune checkpoints have been studied in depth. Some studies have shown that PD-1/PD-L1 blockade has unprecedented efficacy in B cell lymphomas and can be used as an alternative  TNFRSF14  NRP1  LAIR1  CD244  LAG3  CD200R1  ICOS  CD40LG  CTLA4  CD48  CD28  HAVCR2  ADORA2A  TNFRSF4  CD276  CD80  LGALS9  TNFSF14  ICOSLG  TMIGD2  VTCN1  IDO1  PDCD1LG2  TNFSF18  PDCD1  TNFSF9  TNFRSF8  TNFRSF25  CD40  TNFRSF18  TNFSF15  CD274  TNFSF4  CD86  HHLA2  CD70 CD44 TNFRSF9  TNFRSF14  NRP1  LAIR1  CD244  LAG3  CD200R1  ICOS  CD40LG  CTLA4  CD48  CD28  HAVCR2  ADORA2A  TNFRSF4  CD276  CD80  LGALS9  TNFSF14  CD160  ICOSLG  TMIGD2  VTCN1  IDO1  PDCD1LG2  PDCD1  TNFSF9  TNFRSF8  CD27  CD40  TNFRSF18  TNFRSF15  TIGIT  CD274  TNFSF4  CD86  HHLA2  CD70  CD44  to the failure of conventional treatment regimens (23,24). According to Forde et al., PD-1 blockade is remarkably effective in the preoperative treatment of patients with early-stage lung cancer and has few side effects, which can significantly reduce the possibility of tumor recurrence (25). And some related studies found that CTLA-4 blockade can dramatically improve overall survival in melanoma (26,27). Nevertheless, while the promising results of immunotherapy have been achieved, not all patients benefit from it, so we need to explore thoroughly which patients are suitable for immunotherapy in order to maximize the benefits of treatment.
Since most tumors are innately resistant to apoptosis, inducing new types of cell death has become a new strategy for cancer therapy. Necroptosis has gradually aroused attention as a form of proinflammatory programmed cell death. And recent research has shown that immunogenic substances released by necroptosis can exert powerful antitumor immune effects in concert with immune checkpoint blockade (28). Furthermore, a study by Yatim et al. found that necroptosis induces strong cross-priming of the immune system by coupling RIPK1 and NF-κB to activate CD8(+) T cells (29). Therefore, combining necroptosis and immune checkpoint blockade seems to be a direction that holds great promise.
Traditional treatments for glioma include surgery, chemotherapy, and radiotherapy, but these do not significantly improve the survival of patients. The presence of the blood-brain barrier and the unique immune microenvironment of the brain make the treatment of glioma extremely difficult (30). Despite the remarkable results of many novel treatment options in other tumors, few drugs have been approved by the Food and Drug Administration (FDA) for glioma. However, a growing number of studies suggest that immunotherapy may become an important tool for the treatment of glioma in the future. A study by Wainwright et al. revealed that combined targeted inhibition of IDO, CTLA-4, and PD-L1 significantly prolonged the survival of glioma mice (31). In addition, Reardon et al. applied the combination of anti-CTLA-4 and anti-PD-1 therapy in glioma mouse models and cured 75% of the mice; moreover, these mice developed robust immune memory response (32).
In our study, we developed a prognostic signature based on necroptosis-associated lncRNAs. With this model, we assigned a unique risk score to each patient and classified the patients into the high-and low-risk groups. We found that the prognosis of patients in the high-risk group was significantly worse than that of the low-risk group. Validation analysis by ROC suggested that the model has strong predictive power and can assess the prognosis of glioma patients well. In addition, we created a nomogram to improve the clinical applicability of the model. By performing functional enrichment analyses of differential genes between the highand low-risk categories, we found that these differential genes were mainly enriched in immune function-related pathways. Subsequently, we performed immune-related analysis, and most of the immune cells had higher infiltration levels in the high-risk group, such as CD8(+) T cells, macrophages, and Treg cells, suggesting the complexity of the immune microenvironment and potentially high response reactivity in high-risk group patients. The analysis of immune checkpoints revealed that the vast majority of immune checkpoints were highly expressed in patients of the high-risk group, especially the obvious targets like PD-1, IDO1, IDO2, and CTLA-4. And through the exploration of tumor microenvironment, we found that the immune score and stromal score of patients in the high-risk group were significantly higher than those in the low-risk group. All these implied that modulating the immune microenvironment of patients by targeting immune checkpoints has great potential value in the high-risk group patients.
Our model included 7 lncRNAs (FAM13A-AS1, JMJD1C-AS1, LBX2-AS1, ZBTB20-AS4, HAR1A, SNHG14, and LINC00900). Among them, LBX2-AS1 was found to be closely associated with the progression of several cancers. According to Yang et al., LBX2-AS1 is highly expressed in gastric cancer tissues and mainly affects the proliferation of gastric cancer cells (33). According to Cao et al., LBX2-AS1 is considered to be an oncogenic lncRNA in ovarian cancer, and targeted knockdown of LBX2-AS1 significantly reduces the ability of ovarian cancer cells to grow and invade (34). In addition, it was found that high expression of LBX2-AS1 in glioma patients was associated with poor prognosis (35). HAR1A is a tumor suppressor, and in oral cancer patients, knockdown of HAR1A promotes the expression of ALPK1 and leads to oral cancer progression (36). According to Dong et al., SNHG14 plays a critical role in the mechanism of drug resistance in breast cancer, and knockdown of SNHG14 significantly improves trastuzumab efficacy in breast cancer patients (37). According to Zhao et al., SNHG14 is highly expressed in diffuse large B cell lymphoma, and the SNHG14/miR-5590-3p/ZEB1 pathway would enable tumor immune evasion by regulating the immune checkpoint PD-1; targeted inhibition of SNHG14 expression is likely to improve the efficacy of immunotherapy (38).  Figure 7: Investigation of the tumor microenvironment. Different stromal score and immune score between the high-and low-risk groups. 16 BioMed Research International At the same time, there are some limitations in this work. First, most of the data we analyzed were obtained from publicly available databases, so it is necessary for us to conduct more profound in vivo and in vitro trials to validate our results. Second, some lncRNAs in our model have not been studied in depth and await further exploration. This motivates us to keep exploring further in order to contribute to the treatment of glioma.

Conclusions
In conclusion, the necroptosis-associated lncRNA signature we created can accurately predict patient's prognosis, and the combination of necroptosis-related lncRNAs and immune checkpoints blockade has the potential to be a promising new therapeutic option. Our study provides new ideas for individualized treatment of glioma patients.