Necroptosis-related lncRNAs and Hepatocellular Carcinoma Undoubtedly Secret

Background: The current study demonstrates that necroptosis is an important mechanism of carcinogenesis. However, the predictive value of necroptosis-associated long non-coding RNAs (LncRNAs) in hepatocellular carcinoma has not been demonstrated. The purpose of this study is to apply necroptosis-related lncRNAs to construct a predictive signature to predict the prognosis of hepatocellular carcinoma patients. Methods: The clinical and RNA-seq data were downloaded using The Cancer Genome Atlas (TCGA) database. Univariate and multivariate COX analysis was used to screen out suitable necroptosis-related lncRNAs, and then predictive signature was constructed. The TCGA data were randomly divided into high- and low-risk groups, and the Kaplan-Meier method was used to analyze the overall survival (OS) of the two groups to verify the predictive signature. Finally, a prognostic correlation model for predicting disease-freesurvival (cid:0)DFS(cid:0)in hepatocellular carcinoma(cid:0)HCC(cid:0) was constructed and validated. Results: We had screened out 9 necroptosis-related lncRNAs (BACE1-AS, LINC01188, LUCAT1, PI3KCD-AS2, Z83851.1, AC009283.1, AC012360.2, AC015908.3, AC103760.1), which were used to construct a predictive signature, and draw the receiver operating characteristic (ROC) curve by the high- and low-risk group. It was found that the area under the curve (AUC) of risk score was 0.874, and the AUC values of 1-,3-,and 5-years were 0.8, 0.759, 0.787. The TCGA data were randomly divided into two cohorts. In two cohorts, The OS of the high-risk groups were signicantly lower than that of the low-risk groups, and the AUC of the ROC curves of the two cohorts were 0.851, 0.804, 0.802 and 0.735, 0.716, 0.76 at 1-, 3-, and 5-years. After quantifying immune cell subsets and related functions, it was found that the inltration of active dendritic cells (aDCs), macrophages, mast cells, natural killer cells (NK), T regulatory cells (Tregs) were signicantly different, and the expressions of immune checkpoints CD86, LAIR1, CTLA4, VTCN1, TNFRSF18, CD80, CD276, HHLA2, TNFSF4, TNFRSF8, TNFRSF4, TNFRSF9, LGALS9, HAVCR2 and TNFSF15 were also signicantly different. Conclusion: This predictive signature can accurately predict the prognosis of hepatocellular carcinoma patients and provide guidance for the clinical treatment of hepatocellular carcinoma patients.


Introduction
Liver cancer is one of the most common malignant tumor in the world and the sixth most common cancer type that seriously threatens human health [1]. Traditional chemotherapy does not work well for patients with advanced liver cancer [2] Although immunotherapy related studies CHECKMATE-040 [3] KEYNOTE-224 [4] have shown good e cacy, not all patients can bene t from immunotherapy. Therefore, better diagnostic and therapeutic methods for HCC need to be discovered. In 2014, Yuan et al. [5] found that the expression of long non-coding RNA (lncRNA) is related to the occurrence and development of hepatocellular carcinoma. Subsequently, more and more studies have found that lncRNAs are closely related to tumor angiogenesis, tumor metastasis, and tumor treatment effects of hepatocellular carcinoma [6][7][8][9][10].
The main function of lncRNA is to regulate gene expression [11]. At present, people are mainly classi ed according to their positional relationship with adjacent coding genes [12]. According to whether lncRNAs overlap with the sequence of coding genes, they are divided into two categories: gene-related lncRNAs and intergenic lncRNAs. According to their transcription sites and molecular characteristics, they are divided into intergenic long non-coding RNAs, sense or reverse that overlap with protein-coding genes.
Sense transcripts, enhancer RNAs, and promoter upstream transcripts [13]. These lncRNAs can bind to chromatin regulators [14]. It also can bind to DNA [15]. These functions are closely related to the occurrence and development of tumors.
Many lncRNAs have been reported to be involved in inducing necroptosis in HCC cells [16,17]. At the same time, some studies have found that necroptosis is closely related to the development and treatment of liver cancer [18]. However, the predictive value of necroptosis-related lncRNAs in hepatocellular carcinoma has not been demonstrated. The purpose of this study was to apply necroptosis-related long noncoding LncRNAs to construct a predictive model to predict the prognosis of patients with hepatocellular carcinoma and validate it.

Patients and datasets
Downloading the necroptosis-related passway in the Kyoto Encyclopedia of Genes and Genomes (KEGG) and consulting the latest literature to obtain 123 necroptosis-related genes. Downloading the corresponding clinical and prognostic data for The Cancer Genome Atlas liver cancer (TCGA-LIHC) dataset and the fragments per kilobase of transcript per million mapped reads (FPKM)-standardized RNAseq data from the TCGA website (https://portal .gdc.cancer.gov/); obtained the LncRNA expression of 374 hepatocellular carcinoma tissues and 50 normal liver tissues, and obtained the survival time and clinical data of hepatocellular carcinoma patients. We obtained disease-free survival (DFS) data for 310 hepatocellular carcinoma patients from the cBioPortal database (https://www.cbioportal.org/). Immune in ltration results from the TIMER database(https://cistrome.shinyapps.io/timer/). Immunohistochemical results from proteinatlas database(https://www.proteinatlas.org/). gene expression pro le results from GEPIA database(http://gepia.cancer-pku.cn/index.html). All experimental procedures are shown in Figure  1 2.2. Functional enrichment analysis of differentially expressed necroptosis-related genes Using |log 2 fold change(FC)>1|and a false discovery rate (FDR)<0.05 as screening criteria to obtain necroptosis-related differentially expressed genes (DEGs). We performed Gene Ontology (GO) and KEGG analyses in the"ggplot2"package.

The necroptosis-related lncRNA predictive signature
Using the "limma" package to calculate the correlation between necroptosis-related genes. Using the correlation coe cient p<0.001 and |R2|>0.4 as the screening criteria, a total of 765 necroptosis-related lncRNAs with expression values were obtained.Univariate COX analysis was used to obtain necroptosisrelated lncRNAs related to the prognosis of hepatocellular carcinoma patients, and then multivariate COX analysis was used to obtain necroptosis-related lncRNAs, which laid the foundation for the construction of necroptosis-related lncRNAs Predictive Signature. The computational formula used for this analysis was as follows Coef represents the coe cient value, and x represents the selected necroptosis-related lncRNAs. This formula was used to calculate the risk score for each hepatocellular carcinoma patient.

Nomogram
A Nomogram that can predict 1, 3, and 5-year survival of hepatocellular carcinoma patients was constructed using risk scores and clinicopathological features such as age, gender, stage, and TNM stage. A calibration curve was used to test whether the predicted survival rate was consistent with the actual survival rate.

Functional enrichment analysis
Hepatocellular carcinoma patients were divided into high-and low-risk groups according to the median risk score. Which pathway genes were predominantly enriched were analyzed using Gene set enrichment analysis (GSEA). Nominal p<0.05 and FDR<0.25 were the thresholds for statistical signi cance.

Statistical analysis
All statistical analyses were performed with R software (version 4.1.2). The expression of necroptosisrelated DEGs in cancer tissues and normal tissues were analyzed by Wilcoxon test. Univariate Cox regression was used to analyze the relationship between necroptosis-related lncRNAs and overall survival (OS), and multivariate Cox analysis was used to screen necroptosis-related lncRNAs to construct predictive signals. The Kaplan-Meier method and LOG-RANK test were used to analyze the OS of patients in high-and low-risk groups. The receiver operating characteristic (ROC) curve was drawn with the "Survival ROC" software package, and the area under the curve (AUC) value was determined. SsGSEA uses the "GSVA" software package.

Results
3.1. Functional enrichment analysis of differentially expressed necroptosis-related genes 67 necroptosis-related DEGs were obtained by differential analysis of the data screened from the TCGA database, including 63 up-regulated genes and 4 down-regulated genes ( gure 2A, S1, S2). We further analyzed the regulatory relationship between these differentially expressed genes( gure 2B). These genes were then subjected to KEGG pathway and GO enrichment analysis. The KEGG enrichment analysis showed that the DEGs were mainly enriched in Necroptosis, NOD-like receptor signaling pathway, Linoleic acid metabolism, In ammatory mediator regulation of TRP channels, alpha-Linolenic acid metabolism, Lipid and atherosclerosis, GnRH signaling pathway, Shigellosis, Arachidonic acid metabolism ( gure 2C, S3, S4). GO enrichment analysis showed that In the biological process category, the DEGs were mainly enriched in regulation of I-kappaB kinase/NF-kappaB signaling, I-kappaB kinase/NF-kappaB signaling, programmed necrotic cell death, In the cellular components category, the DEGs are mainly enriched in CD40 receptor complex, autolysosome, secondary lysosome, In the molecular function category, the DEGs are mainly enriched in phospholipase A1 activity, calcium-dependent phospholipase A2 activity, ubiquitin protein ligase binding gure 2D, S5, S6 .

Validation of the necroptosis-related lncRNA predictive signature
Calculateing the risk score of all patients, Then we divided patients into two groups by the median risk score. As risk scores increased, so did the number of patient deaths ( gure 3A). We then performed PCA principal component analysis, which showed that patients grouped according to risk lncRNA was more speci c than all gene, all lncRNA, and necroptosis-gene ( gure 3B). Kaplan-Meier analysis method was used to analyze the OS time of high-risk group and low-risk group to determine the value of risk score in predicting the prognosis of patients with hepatocellular carcinoma. The result showed that the OS of the high-risk group was signi cantly shortened ( gure 3C, p 0.001). Subsequently, Univariate COX regression analysis showed that tumor stage, T grade, N grade and risk score were the factors affecting the prognosis of patients ( gure 3D, p 0.05). Multivariate COX regression analysis showed that only risk score was an independent risk factor affecting the prognosis of patients ( gure 3E, p 0.05). According to the drawn ROC curve, it was found that the AUC of the risk score was 0.874, indicating that the risk score is better than other clinical features for the prognosis of hepatocellular carcinoma patients ( gure 3F). And the AUC of 1-year, 3-year and 5-year were 0.8, 0.759, 0.787, indicating that the risk score has good accuracy in predicting survival time( gure 3G). Differences in lncRNA expression between high-and lowrisk groups were analyzed according to clinical characteristics ( gure S7). A nomogram was constructed by clinical features and risk score, which predicts 1-, 3-, and 5-year prognosis in patients with hepatocellular carcinoma ( gure 3H). Calibration curves show good agreement between predicted 1-, 3-,and 5-year survival rates and patients' actual OS ( gure 3I).
3.4. The accuracy of the predictive signature in predicting the prognosis of hepatocellular carcinoma patients with different clinical characteristics In order to verify the accuracy of the Predictive Signature in predicting the prognosis of hepatocellular carcinoma patients with different clinical characteristics, hepatocellular carcinoma patients were grouped according to age, gender, stage, T stage and grade. The OS of low-risk group were signi cantly higher than the high-risk group( gure 4A). These results suggest that predictive features can predict the prognosis of patients with hepatocellular carcinoma regardless of clinical features.

Validation of predictive signature
To verify the validity of the predictive signature, the TCGA data were randomly divided into two cohorts, and the clinical characteristics of two cohorts were shown in Table 1. The results were the same as expected, and the OS of patients in the high-risk group in the rst cohort was signi cantly lower than that in the low-risk group ( gure 4B p 0.01). Another cohort also showed the same result( gure 4C p 0.01). And the AUC of the ROC curve for the rst cohort at 1-, 3-, 5-year were 0.851, 0.804, 0.802 ( gure 4D). The AUC of the ROC curve for another cohort at 1-, 3-, 5-year were 0.735, 0.716, 0.765 ( gure 4E). Box plots were then drawn according to different clinical characteristics and risk scores( gure 5A), This suggests that there are differences in risk scores for different clinical features. 3.6. Gene enrichment analysis Differences in prognosis of patients in high-and low-risk groups were investigated using GSEA. We found that the high-risk groups had cell cycle passway, complement and coagulation cascades passway, drug metabolism cytochrome p450 passway, metabolism of xenobiotics by cytochrome p450 passway, oocyte meiosis passway, peroxisome passway, fatty acid beta oxidation, azurophil granule, aromatase activity and purine metabolism passway signi cant enrichment ( gure 5B, 5C), It shows that high-risk patients are closely related to tumors and a variety of substances metabolism-related pathways.

3.7.The correlation between risk scores and immune cells and functions
To further explore the correlation between risk scores and immune cells and function, we quanti ed immune cell subsets and associated functions. The results showed that active dendritic cells (aDCs), macrophages, mast cells, natural killer cells (NK), T regulatory cells (Tregs) were different between highrisk and low-risk groups. salience ( gure 6A). Likewise, there were signi cant differences in immune function between the two groups ( gure 6B). Then we analyzed the differences of different immune checkpoints between the two groups, and found that the expressions of CD86 , LAIR1, CTLA4,  VTCN1,TNFRSF18,CD80,CD276,HHLA2, TNFSF4, TNFRSF8, TNFRSF4, TNFRSF9,LGALS9,HAVCR2 and TNFSF15 were signi cantly different( gure 6C). To further explore the relationship between 9 necroptosis-related lncRNAs (BACE1-AS, LINC01188, LUCAT1, PI3KCD-AS2, Z83851.1, AC009283.1, AC012360.2, AC015908.3, AC103760.1) and immune cell in ltration, we analyzed their co-expression Expression of necroptosis-related genes in tumor tissue, among which H2AFX, TRAF2,HSP90AB1, CDKN2A, PLK1, and FTL are highly expressed in tumor tissue. The relationship between them and immune cell in ltration was also analyzed online using the TIMER database, and the results showed that they were closely related to a variety of immune cell in ltration( gure S8, immune in ltration level from the TIMER database(https://cistrome.shinyapps.io/timer/).).

Analysis of Gene Mutation in Patients with Hepatocellular Carcinoma
In order to further explore the reasons for the occurrence and development of hepatocellular carcinoma, we conducted gene mutation analysis( gure 6D). First, we analyzed the genes with the highest mutation rates. The results show that these mutated genes were signi cantly associated with necroptosis-related genes( gure 6E). We also analyzed the bases of gene mutation and found that the most common T base changed to A base( gure 6F). These mutated genes are enriched in RTK-RAS, WNT, NOTCH, PI3K and other pathways( gure 6G). Subsequent drug-mutated gene correlation analysis showed that the therapeutic effect of DRUGGABLE GENOME and CLINICALLY ACTIONABLE may be better( gure 6H). Finally, we found that the heterogeneity of hepatocellular carcinoma is more obvious, which suggests that the drug resistance of liver cancer will be very strong( gure 6I).
3.9. the necroptosis-related lncRNA predictive signature for DFS Then we constructed the necroptosis-related lncRNA predictive signature for DFS.We collected DFS data on hepatocellular carcinoma patients from the cBioPortal database, including 310 patients. After univariate Cox regression analysis, 33 necroptosis-related lncRNAs were found to be signi cantly associated with DFS in patients with hepatocellular carcinoma. After multivariate Cox regression analysis, 7 necroptosis-related lncRNAs were obtained, which were used to construct predictive signatures. Risk score=(0.357×ZFPM2-AS1 expression)+ (-0.332×AC026401.3 expression)+ (0.596×SREBF2-AS1 expression)+ (0.293×AL391427.1 expression)+ (0.405×MIR4458HG expression)+ (0.523×AC005229.4 expression)+ (0.512×LYRM4-AS1 expression). The risk score for each patient was calculated according to the formula, and the patients in the entire dataset were divided into high-risk and low-risk groups based on the median. Kaplan-Meier survival curve analysis showed that the DFS in the high-risk group was signi cantly shorter than that in the low-risk group ( gure S9 p 0.001). The AUCs for 1-, 3-, and 5-year survival were 0.677, 0.728, and 0.726 ( gure S10). To validate the predictive signature, 274 patients were randomized into two cohorts. The validation results were the same as expected, and the DFS of patients in the high-risk group in the rst cohort was signi cantly lower than that in the lowrisk group ( gure S11 p 0.01). Another cohort also showed the same result ( gure S12 p 0.01). And the AUC of the ROC curve of the rst cohort at 1-, 3-, 5-year were 0.726, 0.749, 0.734 ( gure S13). The AUC of the ROC curve of another cohort could reach 0.602, 0.719, 0.766 at 1-, 3-, 5-year ( gure S14).

Discussion
Hepatocellular carcinoma is the main pathological type of primary liver cancer [19]. The occurrence of hepatocellular carcinoma is related to factors such as chronic hepatitis virus infection, long-term exposure to carcinogens (such as a atoxin, etc.), excessive drinking, non-alcoholic fatty liver disease, hemochromatosis, and α-1 antitrypsin de ciency [20]. In addition, studies have shown that necroptosis plays a key role in the development of hepatocellular carcinoma. Our study also demonstrated that necroptosis-related lncRNAs can accurately predict the prognosis of HCC patients.
There are various ways of cell death, including apoptosis, autophagy, ferroptosis and necroptosis.
Different studies have proved that these cell death ways are related to the occurrence of tumors [21,22]. However, unlike apoptosis, necroptosis is a necroptosis complex formed by the combination of RIPK1 and RIPK3 after apoptosis is inhibited, which induces high expression of MLKL and forces cells to lose cell membrane integrity [23]. Necroptosis can promote tumor cell extravasation and migration For example, tumor cells use the expressed amyloid precursor as a mediator and induce necroptosis of endothelial cells through DR6 on endothelial cells, leading to tumor cell extravasation and metastasis [24]. When RIPK1/RIPK3 mediates tumor necroptosis, the released soluble factor CXCL1 promotes tumor growth and invasion through an immunosuppressive tumor microenvironment induced by suppressive macrophages [25]. Genetic inactivation of RIPK3 or TNFR1 alleviates clinical symptoms of colitis and cancer, suggesting that RIPK3-mediated necroptosis promotes chronic in ammation and colorectal cancer development [26,27]. In addition, sorafenib is a multi-targeted kinase inhibitor, which is mainly used to treat kidney cancer, liver cancer, thyroid cancer and some leukemia patients. Sorafenib targets the necroptotic pathway in addition to the Braf/Mek/Erk and VEGFR pathways and alleviates RIPK1/3-mediated pathology in acute in ammation [28]. Our experiments also demonstrated a close relationship between necroptosis and hepatocellular carcinoma.
Many studies have demonstrated that lncRNAs can act as key regulators to affect tumor invasion and metastasis [29][30][31], tumor cell proliferation and apoptosis [32][33][34][35],tumor angiogenesis [34,[36][37][38] and regulate tumor drug resistance [39,40]. More and more studies have shown that abnormally expressed lncRNAs are key factors in cancer detection and can be used as non-invasive tumor markers to play an irreplaceable role in the early diagnosis of hepatocellular carcinoma [7,[41][42][43]. Further research has found that lncRNA can be used as a key breakthrough point in molecular targeted therapy. Silencing, blocking, destroying oncogenic lncRNA or restoring the function of tumor suppressor lncRNA can effectively inhibit the malignant biological behavior of hepatocellular carcinoma cells and regulate tumor resistance, which are all important for the treatment of hepatocellular carcinoma [10]. A large number of studies have con rmed that abnormally expressed lncRNA is closely related to the prognosis of liver cancer. It can be used as a potential molecular marker to predict the degree of tumor tissue differentiation, lymph node metastasis, TNM stage, and tumor size, which is of great signi cance for the clinical diagnosis and treatment of hepatocellular carcinoma [10]. Recent studies have reported that necroptosis-related lncRNAs play a key role in hepatocellular carcinoma [18]. Our study also demonstrated that necroptosis-related lncRNAs are closely related to the prognosis of patients with hepatocellular carcinoma, and can even construct a model that effectively predicts the prognosis of patients. This proves that there is enough space for research in this direction, and the author's team will continue to pay attention and report on this direction.
The results of immune cell and immune function analysis showed that the high-risk group had higher in ltration levels of Treg cells, higher in ltration levels of NK cells, and higher in ltration levels of macrophages. This is consistent with the results of previous studies. Previous studies have found that hepatocellular carcinoma tissue forms a speci c immune microenvironment at an early stage, and Treg cells can promote the proliferation and metastasis of tumor cells. Patients with higher levels of Treg cell in ltration have poorer prognosis [44][45][46]. Studies have shown that TAMs isolated from the liver cancer microenvironment are mainly M2 type, and more evidence supports that TAMs with M2 phenotype promote tumor progression through complex autocrine and paracrine pathways closely related to tumor malignant proliferation, invasion and metastasis [34,47]. There is also a clear relationship between immune cells and immune function and patient prognosis, and some studies have even established an immune prognosis model to identify high-risk patients with low survival rates [48][49][50]. Our study also found that there were signi cant differences in the expression of immune checkpoints CD86, LAIR1,  CTLA4, VTCN1, TNFRSF18, CD80, CD276, HHLA2, TNFSF4, TNFRSF8, TNFRSF4, TNFRSF9, LGALS9, HAVCR2 and TNFSF15 between the high-risk group and the low-risk group. Related studies can also be conducted to explore whether these immune checkpoints can be used as diagnosis and treatment's target.
However, our study also has certain limitations. We only use data from the TCGA database for validation, so we also need real-world clinical data for validation to test the accuracy of the predictive model. Second, the mechanism of action of necroptosis-related lncRNAs in hepatocellular carcinoma needs further experimental veri cation. In conclusion, the necroptosis-related lncRNA model can independently predict the prognosis of hepatocellular carcinoma patients, and provide a basis for exploring the possible mechanism and clinical e cacy of necroptosis-related lncRNA in hepatocellular carcinoma.

Conclusions
Our study has identi ed the prognostic value of necroptosis-related lncRNAs in hepatocellular carcinoma and established a prognostic prediction signature. The veri cation results are the same as expected. Our study has also found that there are signi cant differences in the expression of immune checkpoints IDO2, TNFRSF14, BTNL2, TNFRSF25, and CD276 between the high-risk group and the low-risk group. Related studies can also be conducted to explore whether these immune checkpoints can be used as diagnosis and treatment's target.  Figure 1 Experimental ow chart Figure 2 67 necroptosis-related DEGs were obtained by differential analysis of the data screened from the TCGA database, including 63 up-regulated genes and 4 down-regulated genes ( gure 2A). We further analyzed the regulatory relationship between these differentially expressed genes( gure 2B). The KEGG enrichment analysis showed that the DEGs were mainly enriched in Necroptosis, NOD-like receptor signaling pathway, Linoleic acid metabolism, In ammatory mediator regulation of TRP channels, alpha-Linolenic acid Divided patients into two groups by risk score . As the risk score increases, the number of patient deaths also increases( gure 3A). The OS time of the high-risk group was signi cantly shortened ( gure 3C, p 0.001). PCA principal component analysis showed that patients grouped according to risk lncRNA was more speci c than all gene, all lncRNA, and necroptosis-gene( gure 3B). Univariate COX regression analysis showed that tumor stage, T grade, N grade and risk score were the factors affecting the prognosis of patients ( gure 3D, p 0.05). Multivariate COX regression analysis showed that only risk score was an independent risk factor affecting the prognosis of patients ( gure 3E, p 0.05). According to the drawn ROC curve, it was found that the AUC value of the risk score was 0.748, indicating that the risk score is better than other clinical factors for the prognosis of hepatocellular carcinoma patients ( gure 3F). And the AUCs of 1 year, 3 years and 5 years were 0.8, 0.759, 0.787 respectively, indicating that the risk score has good accuracy in predicting survival time.( gure 3G). A nomogram was constructed using clinical features and risk score, which predicts 1, 3, and 5-year prognosis in patients with hepatocellular carcinoma ( gure 3H). Calibration curves show good agreement between predicted 1, 3, 5-year survival rates and patients' actual OS ( gure 3I).

Figure 4
The OS of low-risk groups were signi cantly higher than the high-risk groups( gure 4A). The OS of patients in the high-risk group in the rst cohort was signi cantly lower than that in the low-risk group ( gure 4B p 0.01). Another cohort also showed the same result( gure 4C p 0.01). And the AUC of the ROC curve for the rst cohort at 1, 3, 5-years was 0.851, 0.804, 0.802 ( gure 4D). The AUC of the ROC curve for another cohort at 1, 3, 5-years was 0.735, 0.716, 0.765 ( gure 4E).  were different between high-risk and low-risk groups. salience ( gure 6A). Likewise, there were signi cant