Qualification of Necroptosis-Related lncRNA to Forecast the Treatment Outcome, Immune Response, and Therapeutic Effect of Kidney Renal Clear Cell Carcinoma

Background Kidney renal clear cell carcinoma (KIRC) is considered as a highly immune infiltrative tumor. Necroptosis is an inflammatory programmed cell death associated with a wide range of diseases. Long noncoding RNAs (lncRNAs) play important roles in gene regulation and immune function. lncRNA associated with necroptosis could systematically explore the prognostic value, regulate tumor microenvironment (TME), etc. Method The patients' data was collected from TCGA datasets. We used the univariate Cox regression (UCR) to select prediction lncRNAs that are related to necroptosis. Meanwhile, risk models were constructed using LASSO Cox regression (LCR). Kaplan–Meier (KM) analysis, accompanied with receiver operating characteristic (ROC) curves, was performed to assess the independent risk factors of different clinical characteristics. The evaluated factors are age, gender, disease staging, grade, and their related risk score. Databases such as Gene Ontology (GO), Kyoto encyclopedia of genes and genomes (KEGG), and Gene set enrichment analysis (GSEA) were used to search the probable biological characteristics that could influence the risk groups, containing signaling pathway and immue-related pathways. The single-sample gene set enrichment analysis (ssGSEA) was chosen to perform gene set variation analysis (GSVA), and the GSEABase package was selected to detect the immune and inflammatory infiltration profiles. The TIDE and IC50 evaluation were used to estimate the effectiveness of clinical treatment on KIRC. Results Based on the above analysis, we have got a conclusion that patients who show high risk had higher immune infiltration, immune checkpoint expression, and poorer prognosis. We identified 19 novel prognostic necroptosis-related lncRNAs, which could offer opinions for a deeper study of KIRC. Conclusion The risk model we constructed makes it possible to predict the prognosis of KIRC patients and offers directions for further research on the prognostication and treatment strategies for KIRC.


Introduction
Renal cell carcinoma (RCC) is a branch of urologic tumors that extensively occurred in the world. Studies have reported that RCC is the third occurred tumor among the urinary system, the incidence of RCC is next to prostate cancer and bladder cancer [1], and almost 30% of the patients were present with distant metastases when they were diagnosed [2]. The five-year survival rate of metastatic RCC (mRCC) is only 10%, worser than nonmetastatic RCC [3]. KIRC is the most frequent pathological subtype in adults and is responsible for 80%-90% of the RCC cases [1]. In most patients with KIRC, proper surgical method remains the preferred treatment. However, the tumor that is not sensitive to chemotherapy and radiotherapy is more likely to metastasis or recur compared with other pathological RCCs [4]. Fortunately, recent studies have demonstrated that KIRC is sensitive to immunotherapy and some big data on clinical trials have proven its worth in KIRC [5]. KIRC has been reported to be linked to significant infiltration of immune cells, and the clinical outcomes differ based on the type of cell involved [6]. Thus, identifying the cells related to immune factors would prove helpful.
lncRNAs is one type of the transcribed noncoding RNAs (ncRNAs). The length of it was longer than 200 nucleotides and they are distributed widely in cells [7]. However, lncRNAs exhibited vital roles in multifarious functions of gene expression, such as transcription, chromatin organisation, and translation [8]. Recently, some studies reported that tumor-related lncRNAs could regulate the progression of cancer by affecting the tumor microenvironment (TME) [9], cell differentiation [10], and apoptosis [11]. In addition, lncRNAs could significantly influence the immune system, including immune cell infiltration and immune activation [12]. Some studies reported that the lncRNA LINK-A could downregulate antigen presentation by inactivating the PKA pathway [13]. Moreover, the clinical progression and prognosis of some tumors, including lung cancer, prostate cancer, and BC, are associated with dysfunction of lncRNAs [7]. For instance, the lncRNA HOXD-AS1 was found to have high expression in castration-resistant prostate cancer (CRPC) cells. The multiplication could be inhibited by knockdown of HOXD-AS1 and it could be sensitive to chemotherapy after knockdown. Therefore, HOXD-AS1 showed important role in cancer development [14].
With the development of bioinformatics, several studies have been published recently regarding signature construction based on lncRNAs to handle the therapeutic effect on patients with KIRC. For instance, Sun et al. constructed a 5-immunerelated-lncRNA signature to distinguish whether KIRC patients prognosis is good or not. In addition, the study analyzed the relationship between lncRNA and mRNA to find out the behavior and relationship of these RNAs [15]. Cui et al. reported that seventeen autophagy-associated lncRNAs were successfully identified and a risk profile associated with KIRC prognosis was constructed. This feature is a valid prognostic indicator and not dependent on other features for patients with KIRC [16]. However, the prediction of lncRNAs associated with necroptosis in KIRC and their relationship with immune status has not been clearly described.
In recent years, there are many ways of cell death that have been discovered and via a number of different pathways, including apoptosis, necrosis, programmed necrosis, pyroptosis, iron death, and autophagy. Apoptosis, which is the well-known programmed cell death, had characteristic morphological change with a number of specific biochemical processes. Necrosis is the uncontrolled cellular death, which is often followed by spillage of the cellular contents into surrounding tissues. For the other forms of cell death, for example, pyroptosis is also a kind of programmed cell death with collateral damage (nuclear integrity is maintained) and autophagy, which is a mechanism for both killing stressed cells and to recycle cellular components. Necroptosis is a freshly detected mechanism of cell programmed death mediated by RIP1, MLKL, and RIP3 [17,18]. More and more studies are available suggesting that necroptosis is caught up in various diseases, such as cardiovascular disease, cancers, and neuroinflammation [18][19][20]. Additionally, a recent study showed that necroptosis may boost the cancer metastasis and T cells death in tumors [21]. Necroptosis serves as one of the programmed cell death in the cell, it contains the features of necrosis combined with apoptosis, suggesting it might cause and enhance antitumor immunity of tumors [17]. Park et al. have found that the key regulatory genes in necroptosis could influence the therapeutic effect in non-small-cell lung cancer [22]. In alcoholic cirrhosis, RIPK3-mediated necroptosis was always associated with poor prognosis [23]. Nonetheless, how necroptosis affect the prognosis and inflammation mechanism in KIRC is not yet clear.
Here we established risk signatures to explore the connection between necroptosis-related lncRNAs (NRLs) and the prognosis of KIRC. In addition, we studied how NRLs influenced the tumor microenvironment (TME) and their drug sensitivity in KIRC. We have provided novel prognostic predictors and data for a clearer understanding of the immune infiltrates of necroptosis in patients with KIRC.

Materials and Methods
2.1. Data Availability. The patients' material and their related RNA sequencing data were downloaded from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/) database. Transcribed RNA data were obtained from the fragments per kilobase million (FPKM) for our study. The lncRNAs genes were analyzed using the GENCODE project (https://www.gencodegenes.org/) [24]. Patients with unavailable survival information and incomplete data were excluded.

2.2.
Identification of Genes Associated with Necroptosis. 67 mRNAs related to necroptosis were extracted for identification [25]. We performed the Pearson correlation coefficient analysis in R software (version 4.0.4) to determine the lncRNAs that has a relationship with the pyroptosis-related genes. A correlation coefficient (jRj) value larger than 0.5 defined that a strong correlation exist, and p value less than 0.01 was regarded as the difference was different. The protein-protein interaction (PPI) network of the necroptosisrelated genes was analyzed using the Search Tool for the Retrieval of Interacting Genes (STRING) (https://string-db .org/). Thereafter, PPI network was observed by the Cytoscape software (version 3.7.1).

Qualification of the Necroptosis-Associated lncRNA
Prognostic Signature. The association between NRLs expression and survival data was assessed by using UCR analysis to identify necrosis-associated lncRNAs. NRLs which has been found to have a significant relationship (p < 0:05) were chosen as the necroptosis-related lncRNAs for KIRC. Subsequently, LCR analysis with the 'glmnet' package was applied to establish a prediction model of possible genes. The following formula could be utilized to calculate risk score:     Journal of Oncology      Journal of Oncology An individual risk score was assigned to each patient. Next, KIRC patients were separated as different risk groups using the median cut-off of risk score according to the risk model. Protective and risk prognostic factors were determined using the hazard ratios (HR) by the UCR and the multivariate Cox regression (MCR). The factor was considered risky when HR was >1 and protective when HR was <1.

Survival and ROC Analysis.
We analyzed survival with the R packages survival and survminer, and the differences were distinguished via the KM analysis. We calculated whether our model for different overall survival (OS) is sensitive and specific using the package timeROC of R (version 4.0.4). In addition, we used timeROC to evaluate the independent risk factors of different clinical factors, including age, gender, stage, grade, and risk score.

Construction of Alignment Diagram and PCA of the Risk
Genes. An alignment diagram was created on the basis of the NRLs with 'rms' package to evaluate the various years OS of KIRC patients. Plotting calibration curves was performed to estimate the accuracy of alignment charts. PCA was utilized to categorize the patients into groups according to the NRLs.
2.6. GO, KEGG, GSEA, and ssGSEA Analysis. For bioinformatics analysis, GO and KEGG were used to search possible biological characteristics that may influence the risk groups, including the changed signaling pathway. GSEA was used to explore the immune-related pathways. The ssGSEA was accompanied with the GSEABase package to explore the immune and inflammatory infiltration profiles.

Effectiveness of the Necroptosis-Related lncRNA
Trademark in Clinical Trial. The effectiveness of immunotherapy on KIRC was estimated using TCIA. Relationship between the risk score and immunotherapy sensitive genes including PDL1, PD1, CTLA4, and TIGIT was also checked. The IC 50 value of chemotherapeutic agents was selected to explore the response of KIRC to first-line targeted therapy based on the R package 'pRRophetic'.

Identification of Genes Related to Necroptosis.
We downloaded a total of 15,142 lncRNA expression profiles using the R package. We screened 67 necroptosis-related mRNAs analyzed using the Pearson correlation coefficient based on jRj larger than 0.5 and p value smaller than 0.01 to identify NRLs. At last, 2,180 NRLs were exported. Following this, we performed "limma" R package to get 428 DE necroptosis-related lncRNAs difference from the tumor tissues and normal tissue samples (Figures 1(a) and 1(b)). Then we used GO and KEGG to search the potential

3.2.
Construction of the Prognostic Signature. Following this, UCR analysis was utilized to separate the lncRNAs that have prognosis functions from the DE NRLs, and 348 lncRNAs have significant association with OS in KIRC patients (Supplementary Table S1). From a total of 348 lncRNAs, we identified 19 lncRNAs by LCR to build up the prediction model ( Figure 2 indicate UBE2Q1-AS1 could be a risk factor of prediction in KIRC (Figure 2(c)). We separated 258 patients equally in the high and low-risk group according to the median. We performed the PCA to know the risk patterns in KIRC to estimate the effectiveness of the risk model (Figures 2(d)-2(f)). We can see that a risk model containing 19 lncRNAs had great efficiency to separate patients into different risk groups.

Survival Analysis and Proof of the NRLs Trademark.
We took the KM survival analysis to figure out the OS of the risk signature. We found that the high-risk group was more likely to die than the other group (p < 0:001, Figure 3(a)). We also approved the accuracy of the risk model with the ROC curve. The AUC value of one-year OS was 0.763, while the value for three-and five-year OS was 0.758 and 0.804, respectively. These results indicated that the prediction risk model can precisely forecast the OS (Figure 3(b)). We also found that with an increase in the risk score, the high-risk group has more possibilities to die (Figures 3(c)-3(d)). The NRLs expression in the risk signature was also visualized ( Figure 3(e)).

Affirmation of the NRLs Trademark.
To confirm the truthfulness of our risk trademark in predicting the prognosis in KIRC, the KM survival analysis was further executed. Better OS was found in the low-risk groups (Figure 4(a)).
The AUC values of ROC curve suggested the well predictivity of the risk trademark ( Figure 4(b)). As the risk scores were increasing, the more patients were dead (Figures 4(c)-4(d)). NRLs expression was observed in the testing set ( Figure 4(e)), we found different set expressed more in various groups.       Figure 5(c)). The Chi-squared analysis suggested that higher risk seems to have higher levels of grade, AJCC stage, T stage, and M stage ( Figure 5(d)). We built a nomogram model using risk scores to predict OS in KIRC patients ( Figure 5(e)). The predictions of OS were effective as presented in the calibration plot ( Figure 5(f)).

Autonomous Prognostic Factors and
The results suggested that both the risk and nomogram models were accurate. The prognostic value of various clinical features was demonstrated in the DCA plot ( Figure 5(g)). To detect the prognosis in diverse clinical elements, we evaluated the survival differences of KIRC patients in various risk groups. Except for N stage, patients in the low-risk groups have longer OS than the other group. (Figures 6(a)-6(n)).
3.6. Functional Assessment of the Risk Feature. We performed GSEA to assess the action of the risk model. Processes significantly influenced the development of cancer, including MYC targets V2, DNA repair, IL-6/JAK/STAT3 signaling, and immune response. These processes existed more in the high-risk group, whereas metabolic processes were embellished in the low-risk group (Figures 7(a)-7(t)). The high-risk group can upregulate several pathways and processes linked with tumor progression and immune response, suggesting that necroptosis might influence the treatment outcomes of immunotherapy according to analysis using GSEA.
Then "CIBERSORT" was applied to determine how immune cell is expressed, indicating immune cells such as plasma cells (p < 0:001), T cells CD8 (p = 0:02), T cells follicular helper (p = 0:004), T cells regulatory (Tregs) (p < 0:001), NK cells resting (p = 0:044), and macrophages M0 (p = 0:04) expressed more in the high-risk groups, which further confirmed our conclusions (Figure 8(e)). Additionally, as the risk scores were increasing, there were higher levels of the aforementioned cells (Figure 8(f)).      Immunotherapy scores data was collected from TCIA database to differentiate the immune responses of the two groups. We found that patients without CTLA4 and PD-1 expressed had no differ-ences in immunotherapy scores (Figure 9(a)). But either one or two of them positive would lead to greater immunotherapy scores (Figures 9(b)-9(d)). Subsequently, we examined whether there exist a relationship in the risk groups and  Journal of Oncology       (Figures 9(e) and 9(f)). But more low-risk patients are sensitive to sorafenib (p = 0:048), sunitinib (p < 0:001), and temsirolimus (p < 0:001) (Figures 9(g)-9(i)). In conclusion, our prognostic model can be a potential indicator of the effectiveness of clinical treatment.
3.9. The Correlation between our Risk NRLs and the Related Genes. We analyzed how our prognostic NRLs could influence each other, and what interested us was that the levels of the prognostic NRLs including IGFL2-AS1, LINC01943, U62317.1, LASTR, LINC01126, AC026401.3, MYOSLID, APCDD1L-DT, AL162586.1, and NARF-IT1 were positive in increasing risk scores (Figure 10(a)). Both lncRNA-    27 Journal of Oncology mRNA expressed network was built up according to our risk signature to find the related necroptosis genes (Figure 10(b)). The Sankey diagram demonstrated the protective and risk factors of NRLs and the related mRNAs ( Figure 10(c)). Finally, we explored the biological function of the related mRNAs, which were significantly associated with the procession of the cell death and progression of the cancer (Figure 10(d)).

Discussion
lncRNAs have been identified as crucial regulators of various kinds of cellular processes since they could function as tumor suppressors. The upregulation of lncRNA will promote the proliferation and invasion of tumor, while knockdown of its expression suppresses this process. It has been reported in many studies that the lncRNAs are altered in many types of cancers, and therefore the aberrant lncRNAs expression levels can be applied as effective diagnostic markers, and deregulated lncRNAs can be used as targets in cancer treatment. In our study, a 19-NRL risk model was built up by us to estimate the prognosis of KIRC patients. The model presented unique advantages. We used the UCR to select prognostic lncRNAs that are related to necroptosis. In the meantime, risk models were constructed using LCR. Finally, 258 patients were equally separated in the high-or low-risk group to know the 19 lncRNAs. The KM and ROC curve analyses were performed to know the treatment effect of KIRC patients, which revealed that the model was a powerful prediction tool. In addition, this model could assess the different clinical characteristics; the evaluated factors are age, gender, disease staging, grade, and their related risk score. The estimated risk score was independent with excellent sensitivity and specificity. Furthermore, via GSEA, the high-risk group was found to enhance tumor development and progression, which confirmed the differences in prognostic property based on classification by risk.
Some previous studies have contributed to the construction of the lncRNA-related model to predict the immune infiltration landscape in KIRC. Chen et al. identified four lncRNA predictive risk scoring models and found that higher risk scores were associated with higher levels of immune infiltration in the KIRC microenvironment. Higher risk score will increase activation of six immune cells, the cell types were mentioned before. [26]. Based on the extensive