Necroptosis-Related lncRNAs: Predicting Prognosis and the Distinction between the Cold and Hot Tumors in Gastric Cancer

Background In the face of poor prognosis and immunotherapy failure of gastric cancer (GC), this project tried to find new potential biomarkers for predicting prognosis and precision medication to ameliorate the situation. Methods To form synthetic matrices, we retrieved stomach adenocarcinoma transcriptome data from Genotype-Tissue Expression Project (GTEx) and The Cancer Genome Atlas (TCGA). Necroptosis-related prognostic lncRNA was identified by coexpression analysis and univariate Cox regression. Then we performed the least absolute shrinkage and selection operator (LASSO) to construct the necroptosis-related lncRNA model. Next, the Kaplan–Meier analysis, time-dependent receiver operating characteristics (ROC), univariate Cox (uni-Cox) regression, multivariate Cox (multi-Cox) regression, nomogram, and calibration curves were made to verify and evaluate the model. Gene set enrichment analyses (GSEA), principal component analysis (PCA), immune analysis, and prediction of the half-maximal inhibitory concentration (IC50) in risk groups were also analyzed. For further discussing immunotherapy between the cold and hot tumors, we divided the entire set into two clusters based on necroptosis-related lncRNAs. Results We constructed a model with 16 necroptosis-related lncRNAs. In the model, we found the calibration plots showed a good concordance with the prognosis prediction. The area's 1-, 2-, and 3-year OS under the ROC curve (AUC) were 0.726, 0.763, and 0.770, respectively. Risk groups could be a guide of systemic treatment because of significantly different IC50 between risk groups. Above all, clusters could help distinguish between the cold and hot tumors effectively and contribute to precise mediation. Cluster 2 was identified as the hot tumor and more susceptible to immunotherapeutic drugs. Conclusion The results of this project supported that necroptosis-related lncRNAs could predict prognosis and help make a distinction between the cold and hot tumors for improving individual therapy in GC.


Introduction
Regrettably, more often than not, most gastric cancer (GC) patients are diagnosed at an advanced stage with a poor prognosis end [1]. As a result, GC is the third most common cause of cancer death globally (8.2% of 9.6 million cancer deaths in 2018 worldwide) [2]. Systemic treatments are the only choice of patients who cannot be surgically treated. Chemotherapeutics and target therapeutics are common systemic treatments and are frequently reported with treatment failure and side reactions [3]. It urges us to find new exploitable therapeutic strategies. Immunotherapy has had transformed the treatment landscape for malignancies and achieved a lot. But just one-third of patients respond to checkpoint inhibitors in most cancers [4]. Besides failing in the induction of cell death, the cold tumor, lacking preexisting immunity, is also the reason for immunotherapy resistance [5]. erefore, it is imperative to study how to augment immunotherapy in GC.
As most tumors have innate apoptosis resistance, the induction of other cell death mechanisms, such as necroptosis, has gradually been recognized as promising therapy strategies [4]. Necroptosis, as a novel programmed form of necrotic cell death different from apoptosis, can enhance CD8 + leukocyte-mediated antitumor immunity by activating RIPK1 and RIPK3 within the tumor microenvironment (TME) [6]. In another report, a necroptotic cancer cellmimicry nanovaccine can potentiate antitumor immunity in mice by inducing expansion of NKG2D + natural killer cells and CD8 + T cells [4]. At the same time, necroptosis works in generating an immunosuppressive TME to promote malignancies via CXCL1 and Mincle, which also hints at necroptosis as a potential immunotherapy target in GC [7].
Long noncoding RNA (lncRNA) can control genes by influencing their translation or interacting directly with proteins and other RNA species [8]. Linc00176 releases miRNAs, such as miR-9 and miR-185, to downregulate target mRNAs leading to necroptosis of hepatocellular carcinoma cells. TRINGS, the p53-inducible lncRNA, can protect cancer cells from necroptosis by inhibiting TRAP-GSK3β-NF-κB necrotic signaling. Cardiomyocyte necroptosis can be regulated by lncRNA through RIPK1/RIPK3 [9]. Besides, lncRNAs have been reported to promote tumors inflammation and help malignancies evade immune destruction [10]. e study of necroptosis-related lncRNA has not been widely mentioned as a potential therapeutic target in GC. erefore, acquiring more necroptosis-related lncRNAs knowledge can help us understand the roles of necroptosis and lncRNAs in immunotherapy clearly. e distinction between the cold and hot tumors and turning a cold tumor into a hot tumor will improve the antitumor effects of immunotherapy. It will bring a breakthrough in immunotherapy, while the mechanisms of other cell death remained to be fully elucidated in GC at this stage. But we still lack a simple and effective method for distinguishing tumors [5]. As lncRNAs are spoken highly of acting as new cancer biomarkers in bodily fluids, we tried to regroup patients based on necroptosis-related lncRNAs and identify the hot tumor effectively for improving prognosis and augmenting precise mediation in clinical practice [10,11].

Acquisition of Information of Patients with GC.
For getting synthetic data matrices about stomach adenomas and adenocarcinomas and normal stomach tissue, the RNA transcriptome datasets (HTSeq-Counts and HTSeq-FPKM) and the relevant clinical information were downloaded from Genotype-Tissue Expression Project (GTEx) (https://www.gtexportal.org/) and e Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). en we converted the FPKM value to the TPM value of the synthetic matrix by data.table, tibble, dplyr, and tidyr R packages. As a result, we got two synthetic data matrices. e Counts value matrix was just for identifying differentially expressed lncRNAs, while the TPM value matrix was for the other analyses. To reduce statistical bias in this analysis, stomach adenomas and adenocarcinomas patients with missing overall survival (OS) values or short OS values (<30 days) were excluded. With relevant clinical information, we retrieved 306 patients and divided them into the train risk group and test risk group randomly by Strawberry Perl and caret R package. e ratio was 1:1.

Selection of Necroptosis-Related Genes and lncRNAs.
e necroptosis gene set M24779.gmt contains eight necroptosis genes, and it was downloaded from the Gene Set Enrichment Analysis (GSEA) (http://www.gsea-msigdb.org/ gsea/index.jsp). In addition, with previous reports about necroptosis, we finally obtained the profile of 67 necroptosisrelated genes (Appendix T1). en we found 5,022 differentially expressed lncRNAs (Log 2 fold change (FC) > 1, false discovery rate (FDR) < 0.05, and p < 0.05) after screening the synthetic data matrix by Strawberry Perl and limma R package [12]. Correlation analysis was performed between 67 necroptosis-related genes and differentially expressed lncRNAs in the combined matrices.

Establishment and Validation of the Risk Signature.
According to the clinical data of GC cases in the TCGA and GTEx, univariate Cox proportional hazard regression analysis was used to screen lncRNAs related to survival from necroptosis-related lncRNA (p < 0.05). en, we made the Lasso regression performed with 10-fold cross-validation and a p value of 0.05 as well as run for 1,000 cycles. For each cycle, a random stimulation was set up 1,000 times in order to prevent overfitting. en a model was established. e 1-, 2-, and 3-year time-dependent receiver operating characteristics (ROC) curves of the model were plotted by the calculation procedure. We calculated the risk score with the following formula: risk score � n k�1 coef lncRNA k * expr lncRNA k , (1) where the coef (lncRNAn) was the short form of the coefficient of lncRNAs correlated with survival and expr (lncRNAn) was the expression of lncRNAs. According to the median risk score, subgroups including low-and high-risk groups were established [12,13]. We used the chi-square test to analyze the relationship between the model and clinical factors in order to evaluate the prognostic value of the constructed model.

Independence Factors and ROC.
We developed univariate Cox (uni-Cox) and multivariate Cox (multi-Cox) regression analyses to evaluate whether the risk score and clinical characteristics were independent variable factors and made ROC to compare different factors in predicting outcome.

Nomogram and Calibration.
With rms R package, the risk score, age, and tumor stage were used to set up a nomogram for the 1-, 2-, and 3-year OS and correction curves based on the Hosmer-Lemeshow test to illustrate whether the prediction outcome showed good consistence with the practical.

e Investigation of the TME and Immune Checkpoints.
According to the result of GSEA, we decided to analyze the immune-cell factors in risk groups. We could calculate the immune infiltration statuses among the GC patients from the TCGA including TIMER, CIBERSORT, XCELL, QUANTISEQ, MCPcounter, EPIC, and CIBERSORT on TIMER2.0 (http://timer.cistrome.org/). In another way, we could download the profile of infiltration estimation for all TCGA tumors on the same website. Wilcoxon signed-rank test, limma, scales, ggplot2, and ggtext R packages were performed in analyzing the differences in immune infiltrating cell content explored, and the results were shown in a bubble chart [13]. Besides, we also made comparisons about TME scores and immune checkpoints activation between low-and high-risk groups by ggpubr R package.

Exploration of the Model in the Clinical Treatment.
en we used the R package pRRophetic to evaluate their therapy response determined by the half-maximal inhibitory concentration (IC50) of each GC patient on Genomics of Drug Sensitivity in Cancer (GDSC) (https://www. cancerrxgene.org/) [14]. 16 Prognostic lncRNAs. For exploring GC response to immunotherapy, we decided to explore potential molecular subgroups by ConsensusClusterPlus (CC) R package based on the prognostic lncRNAs expression [15]. Principal component analysis (PCA), T-distributed stochastic neighbor embedding (t-SNE), and Kaplan-Meier survival were made by Rtsne R package.

Clusters Based on
Besides, we made immunity analysis and drug sensitivity comparison by GSVA Base and pRRophetic R package.

Necroptosis-Related lncRNAs in GC Patients.
e flow of the study was exhibited in Figure 1. From e Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression Project (GTEx) matrix, we obtained 204 normal samples (174 samples from GTEx) and 343 tumor samples. According to the expression of 67 necroptosis-related genes and differentially expressed lncRNAs (|Log 2 FC| > 1 and p < 0.05) between normal and tumor samples, we finally got 387 necroptosis-related lncRNAs (correlation coefficients > 0.4 and p < 0.001) [12,16]. Of them, 194 were upregulated, and the others were downregulated (Figure 2(a)). e network figure and data between necroptosis-related genes, such as AXL and BCL2, and lncRNAs were shown in Figure 2(b) and Appendix D1.

Construction and Verification of the Model.
According to univariate Cox (uni-Cox) regression analysis, we found 16 necroptosis-related lncRNAs significantly correlated with overall survival (OS) (all p < 0.05) and made a heat map (Figures 3(a) and 3(b)). To avoid overfitting the prognostic signature, we performed the Lasso regression on these lncRNAs and extracted 16 lncRNAs related to necroptosis in GC when the first-rank value of Log(λ) was the minimum likelihood of deviance (Figures 3(c) and 3(d)). Besides, we could find 10 lncRNAs were upregulated and the others were downregulated in the Sankey diagram ( Figure 3(e)).
With the risk score formula, the distribution of risk score, the survival status, survival time, and the relevant expression standards of these lncRNAs of patients were compared between low-and high-risk groups in the train, test, and entire sets. ese all indicated the high-risk group had worse prognoses (Figures 4(a)-4(l)). Besides, the conventional clinicopathologic characteristics, age, gender, grade, stage, T, M, and N also performed the same results ( Figure 4(m)).
According to three independent prognostic factors, risk score, age, and TNM stage (all p < 0.05 in multi-Cox), we built a nomogram for predicting the 1-, 2-, and 3-year OS incidences of GC patients ( Figure 5(c)). We also utilized the 1-, 2-, and 3-year calibration plots to attest that the nomogram had a good concordance with the prediction of 1-, 2-, and 3-year OS ( Figure 5

GSEA.
To investigate differences in biological functions between risk groups, we utilized GSEA software to explore the high-risk group in the KEGG pathway in the entire set ( Figure S1A). Seven of the top ten pathways with enrichment in the high-risk group were highly correlated with tumor invasion, and the others were correlated with immunity such as "leukocyte transendothelial migration" (all p < 0.05; FDR < 0.25; |NES| > 1.5) (Figure 6(a)) [17]. erefore, we tried to make an immunity analysis in the model.

3.6.
e Investigation of Immunity Factors and Clinical Treatment in Risk Groups. More immune cells were associated with the high-risk group on different platforms exhibited at the immune cell bubble chart and in the document such as macrophage M1, T cell CD4+ naïve, immune score at XCELL, T cell CD8 + , T cell CD4 + at TIMER, T cell CD4 + at QUANTISEQ and macrophage at MCPcounter and EPIC (all p < 0.05) ( Figure 6(b)) (Appendix D2). Besides, we also found that the higher risk score had more association with immune cells such as dendritic cells resting, which had been reported as a part of immunotherapy in GC ( Figure 6(c)) [18]. All of these showed the high-risk group had a higher immune infiltration status. e high-risk group had a higher immune score and a higher ESTIMAT (microenvironment) score, signifying a different TME from the low-risk group ( Figure 6(d)). Most immune checkpoints also showed better activation in the high-risk group ( Figure 6(e)). It implied that we could choose appropriate checkpoint inhibitors for GC patients regrouped by the risk mode [19]. Consistent with reports, the high-risk group, with a higher immune score, had a lower IC50 of 12 immunotherapeutic drugs such as bryostatin 1 ( Figure 6(f )) [20]. What's more, we could also find that 16 chemical or targeted drugs, which applied to GC therapy, showed lower IC50 in the high-risk group ( Figure S2).

Distinguishing between the Cold and Hot Tumors and Precision Medicine in Clusters.
Referring to previous research, different clusters, known as subtypes, usually showed different immune microenvironments leading to different immunotherapeutic responses [21,22]

necroptosis-related lncRNAs expression (Figures 7(a) and
S3A) [15]. T-distributed stochastic neighbor embedding (t-SNE) indicated two clusters could be distinguished clearly ( Figure 7(b)). In addition, we employed principal component analyses (PCA) to verify that both risk groups and clusters have different PCA (Figure 7(c)). Moreover, cluster 1 had better OS (p � 0.006) in the Kaplan-Meier analysis ( Figure 7(d)). GSEA was also employed to investigate clusters' biological functions. Seven of the top ten pathways with enrichment in cluster 2 were related to immunity (p < 0.05; FDR < 0.25; |NES| > 1.9), such as "natural killer (NK) cell-mediated cytotoxicity" (Figures 7(e) and Figure S1B). NK cells and their cytotoxicity acted as an important role in immunity and cancer [23]. To verify its relations with risk, a chart was also made. Cluster 1 was significantly associated with the low-risk group, and cluster 2 was associated with the high-risk group ( * * means p < 0.01) (Figure 7(f )). e results below about cluster 2 might contribute to patients' immunotherapy in risk groups. Concerning the comparison of the single sample GSEA (ssGSEA) scores for immune cells and immune functions, 15 immune cells, such as CD8 + T cells, and 12 immune functions, such as inflammation-promoting, had more relations with cluster 2 (Figure 7(g)). Cluster 2 was more highly infiltrated by immune cells based on analyses of the different platforms (Figure 7(h)) (Appendix D3). Cluster 2 had a higher immune score and a higher ESTIMAT (microenvironment) score, signifying a different TME from cluster 1 (Figure 7(i)). Almost all the immune checkpoints expressed more activity in cluster 2, such as HAVCR2 (TIM3), LAG3, and CD274 (PD-L1) (Figure 7(j)). CD8 + T cells, the function of inflammation-promoting, high immune score, activation of TIM3, LAG3, and PD-L1 played vital roles in the hot tumor. erefore, we could consider cluster 2 as the hot tumor while cluster 1 as the cold tumor [5,24]. It might result in different immunotherapeutic responses [22,25]. With the notion of the cold and hot tumors, cluster 2 was more susceptible to immunotherapy. With drug sensitivity comparison, we found nine immunotherapeutic drugs, such as shikonin, showed different IC50 solely in clusters as well as 16 chemical or targeted drugs that applied to systemic treatments in GC (Figure 7(k) and Figure S3B) [26,27]. Because of clusters based on these lncRNAs, we might further study immunotherapy responses and potentiate precise medication in GC patients.

Discussion
Immunotherapy can ameliorate the situation of frequently reported treatment failure, but it is not a panacea for all diseases [4]. Because of immunosuppressive TME, some patients had poor immunotherapy responses. erefore, we introduced the notion that cold and hot tumors refer to an immune-based classification of tumors rather than conventional cancer-based for improving immunotherapy. e highly infiltrated tumor with a high immune score usually is regarded as the hot tumor while the noninfiltrated tumor with a low immune score as the cold tumor. e higher activity of checkpoints, such as TIM3 and LAG3, is also the one of characteristics of the hot tumor. In the face of the hot tumor, we can treat patients with T-cell-targeting immunotherapies, microbiome modulation, or other    Journal of Oncology      immunotherapeutic drugs. But it is not easy for the cold tumor because it fails to unleash preexisting immunity with low degree T cells. CD8 + T cells can kill cancer cells by releasing PRF1, GNLY, or GZM and break tolerance as a preexisting immune response, enhancing immunotherapy via the PD-1/PD-L1 immune inhibitory axis. So it is wise to turn a cold tumor into a hot tumor rather than just give other treatments [5,28].
In the study, we built 16 necroptosis-related lncRNAs mode and tried to identify the cold and hot tumors. Patients were regrouped into low-and high-risk groups by the model and made some analyses such as Kaplan-Meier analysis, GSEA, and IC50 prediction. Although we found risk groups could be a guide in predicting prognoses and systemic treatments, we could not identify the hot tumor by risk groups. Referring to reports, molecular subtypes, also known as clusters, are associated with tumor immune suppression and microenvironments [1,29]. Different subtypes have different immune and TME scores leading to different prognoses and immunotherapy responses [22,25]. us, we divided patients into two clusters based on the expression of these lncRNAs [15]. As expected, the two clusters had different immune microenvironments. Cluster 1 had an immunosuppressive TME. At the same time, there were more CD8 + T cells highly infiltrated, more active function of inflammation-promoting, higher immune score, and higher activity of TIM3, LAG3, and PD-L1 in cluster 2, which could be identified as the hot tumor definitely [5,24]. What's more, cluster 2 was more sensitive to immunotherapeutic drugs. Necroptosis-related lncRNAs could not only predict prognosis but also be a guide for individual therapy. Above all, these lncRNAs, as liquid biopsies, could distinguish between the cold and hot tumor briefly and effectively compared with the tumor biopsy by imaging mass cytometry or other experiments [28].
Besides, in the Sankey diagram, we found some of these lncRNAs were related to star genes such as BCL2. FRMD6-AS2 and LINC02657 were associated with AXL, which contributed to immunotherapy by reprogramming the immunological microenvironment and PD-1 inhibitors [30]. BCL2, as a hot point of PD-1 immunotherapy, correlated with AC012409.3 and AC069549.1. It also could reactive apoptosis to overcome immunotherapy failure and keep durable antitumor responses [31]. RNF31 was a potential immunotherapy target of GC and might be regulated by AC004596.1 [32]. TRAF2, associated with AL355574.1,   could augment immunotherapy by increasing the susceptibility of tumors [33]. FRMD6-AS2 suppressed tumor growth, migration, and invasion via the Hippo signaling pathway [34]. And Hippo signaling pathway is engaged in programmed cell death by regulating kinds of target such as YAP [35]. LINC02657 was reported that its overexpression would keep cancer cells from programmed cell death by regulating SART3 [36]. e other lncRNAs were firstly discovered. Newly acquired necroptosis-related lncRNAs knowledge could help us develop a better mechanistic understanding of GC, which would bring a breakthrough into clinical practice.
ere were still some shortcomings and deficiencies though we had utilized many methods to asset our model. As a retrospective study, it was susceptible to the inherent biases of this research paradigm [37]. Although the activation of checkpoints performed significantly between risk groups and clusters, we could not make a comparison of corresponding checkpoint inhibitors IC50, such as PD-1 inhibitors, as a result of insufficient data on GDSC. We had performed internal validation by the test and entire sets in the model, but it was difficult to do external validation for prognoses. Even if we had retrieved all information of GSE84437 series and GSE62254 series matrices from Gene Expression Omnibus (GEO) (https:// www.ncbi.nlm.nih.gov/geo/), we could not get appropriate information of lncRNAs because commercial microarray data had biases and limitations compared with GTEx and TCGA. However, the immune cell bubble and the immune cell heat map showed the results from multiple platforms, which might be recognized as external validation in a sense. Given the above analyses and previous reports, we felt that our model was reasonable and acceptable for future clinical tests [13,38]. Collecting more clinical datasets could help reassert the value of these necroptosis-related lncRNAs, which would be in our plans.

Conclusions
Necroptosis-related lncRNAs could predict prognosis and help propose an exploitable therapeutic strategy by identifying the cold and hot tumors, which would make great progress in individual therapy and improve patients' prognoses. Targeting necroptosis and lncRNAs will be a promising pathway for overcoming systemic treatments failure and expanding the field of immunotherapy. erefore, the mechanisms and relationships, among necroptosis, lncRNAs, immunity, and GC, were worth being fully elucidated and validated.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this study.

Authors' Contributions
Zirui Zhao and Haohan Liu contributed equally to this study. Zirui Zhao and Haohan Liu contributed to the conceptualization; Zirui Zhao contributed to the methodology, analysis, visualization, and original draft preparation; Zirui Zhao, Haohan Liu, Xingyu Zhou, Deliang Fang, Xinde Ou, Jinning Ye, Jianjun Peng, and Jianbo Xu reviewed and edited the manuscript; Jianbo Xu contributed to the project administration and funding acquisition. All authors have read and agreed to the published version of the manuscript.