Abstract Introduction Osteosarcoma is a bone-derived malignancy that often leads to lung metastasis and death. Material and methods The RNA-seq data of TARGET-osteosarcoma were collected from TARGET database. [43]GSE16088 and [44]GSE12865 datasets of osteosarcoma x from Gene Expression Database (GEO) were donwloaded. ConsensusClusterPlus was used for molecular subtype classification. Univariate Cox and Lasso regression was employed to develop a risk model. To analyze the regulatory effects of model feature genes on the malignant phenotype of osteosarcoma cell lines, qRT-PCR, Transwell and wound healing assays were performed. The abundance of immune cell infiltration was assessed using MCP-Counter, Gene Set Enrichment Analysis (GSEA), and ESTIMATE. The Tumor Immune Dysfunction and Exclusion (TIDE) software was employed to evaluate immunotherapy and response to conventional chemotherapy drugs. Results Three clusters (C1, C2 and C3) were classified using 39 necroptosis score-associated genes. In general, C1 and C2 showed better prognosis outcome and lower death rate than C3. Specifically, C2 could benefit more from immunotherapy, while C3 was more sensitive to traditional medicines, and C1 had higher immune cell infiltration. Next, an 8-gene signature and a risk score model were developed, with a low risk score indicating better survival and immune cell infiltration. ROC analysis showed that 1-, 3-, and 5-year overall survival of osteosarcoma could be correctly predicted by the risk score model. Cellular experiments revealed that the model feature gene IFITM3 promoted the osteosarcoma cell migration and invasion. Furthermore, the overall survival of osteosarcoma patients from TARGET and validation datasets can be accurately evaluated using the nomogram model. Conclusions Our prognostic model developed using necroptosis genes could facilitate the prognostic prediction for patients suffering from osteosarcoma, offering potential osteosarcoma targets. Keywords: Necroptosis, Osteosarcoma, Prognosis, Risk model, Tumor microenvironment (TME) 1. Introduction Osteosarcoma is an aggressive bone tumor that frequently occurs to children and adolescents [[45]1]. Lung metastasis is detected in around 15 %–20 % of patients with osteosarcoma at their first diagnosis [[46]2,[47]3]. According to histological findings, osteosarcoma can be divided into high grade, medium grade and low grade [[48]4]. Neoadjuvant chemotherapy has provided more treatment choices than before, but the overall survival of patients with osteosarcoma has not been significantly improved [[49]5]. Hence, finding new prognostic factors and therapeutic targets is necessary for managing osteosarcoma. Necroptosis [[50]6] is a non-caspase-dependent cell death characterized by both necrosis and apoptosis and is regulated by receptor-interacting protein (RIP). RIP1 fulfills a crucial function in the signal transduction of necrotizing apoptosis. Under abnormal conditions of the body, significantly increased production of caspase-8 will initiate the apoptotic signaling pathway [[51]7], which can completely inhibit RIP1-mediated proteolytic activation of procaspase-8, block apoptosis signaling pathway, and then transform apoptosis into necrotizing apoptosis. Different stimuli are associated with necrotizing apoptosis, and TNF-α combined with Z-VAD-FMK (T/Z) is the classic model of inducing necroptosis [[52]8,[53]9]. In 2016, study reported that pancreatic ductal adenocarcinoma has high-expressed RIPK1 and RIPK3, which are the key components of necroptosis, but the suppression of RIPK3-mediated necrotizing apoptosis could delay the cancer progression [[54]10]. Necroptosis is also present in the necrotic areas in mice with breast cancer, and the inflammatory response caused by necroptosis can cause lung metastasis [[55]11]. Apart from these findings, numerous studies have been conducted on necroptosis in osteosarcoma [[56][12], [57][13], [58][14]]. In recent years, tumor diagnostic and prognostic biomarkers are increasingly discovered using public databases. Microarray technology is becoming increasingly useful in gene analysis. Gene expression profile database is a useful tool with significant clinical value in medical oncology [[59][15], [60][16], [61][17]]. Developing risk models for the adjuvant detection and management of multiple cancer types can improve the precision and effectiveness of treatment protocols [[62]18,[63]19]. Thus, in this study, TARGET-osteosarcoma, related osteosarcoma datasets of [64]GSE12865 and [65]GSE16088 were retrieved from the TARGET database and GEO database, respectively, and the potential key genes were identified to develop a prognostic risk model for the cancer. The patterns of immune infiltration were analyzed to provide a novel possibility for the pathogenesis and treatment strategy of osteosarcoma. 2. Materials and methods 2.1. Raw data We downloaded the RNA sequencing (RNA-seq) data of “TARGET-Osteosarcoma” from the TARGET ([66]https://ocg.cancer.gov/programs/target) database. A sum of 79 osteosarcoma samples with 19,533 genes were retained by removing samples without survival time, survival status or clinical follow-up information while retaining mRNAs encoding proteins ([67]Supplementary Table 1). Additionally, the dataset [68]GSE21257 (45 osteosarcoma samples) and [69]GSE39058 (36 osteosarcoma samples) were downloaded from GEO database by [70]GPL10295 and [71]GPL14951 [[72]20].The data from the TARGET-osteosarcoma dataset were used as a training cohort, while those from the [73]GSE21257 and [74]GSE39058 datasets were used as a validation cohort. In total, 74 necroptosis-correlated genes were collected following a previous study [[75]21]. 2.2. Molecular subtypes Firstly, necroptosis score was calculated for TARGET-osteosarcoma using ssGSEA. Genes (cor >0.5 and P < 0.05) related to the necroptosis score were selected by Hmisc package. Univariable Cox analysis was performed to screen the prognostic genes from TARGET-osteosarcoma dataset using the Coxph function of R package survival, and the liminal value was defined when p < 0.05. Using these genes, the samples in the TARGET-osteosarcoma were subjected to molecular subtyping using the R package Consensus Cluster Plus 1.52.0 [[76]22]. Applying Pam arithmetic and “spearman” distance, 500 bootstraps were conducted with each bootstrap containing at least 80 % of the samples from the TARGET-osteosarcoma dataset. The Cluster number k was between the range of 2 and 10, and according to cumulative distribution function (CDF) and AUC, the optimum k was determined. 2.3. Analysis of immune infiltration R software ESTIMATE [[77]23] was employed to calculate overall immunocyte infiltration (Immune Score), stroma level (Stromal Score), and the combination (ESTIMATE Score) for the samples in the TARGET-osteosarcoma cohort using Wilcox.test analysis. 2.4. SsGSEA The GSEA was used to estimate a total of 28 subpopulations of TILs including the major innate immunity-correlated cell types including natural killer T (NKT) cells, macrophages, monocytes, etc. 2.5. TIDE TIDE [[78]24,[79]25] algorithm ([80]http://tide.dfci.harvard.edu) was run to evaluate the myeloid suppressor cells (MDSC), dysfunction of tumor infiltration cytotoxic T lymphocytes (CTL) (Dysfunction), IFNG, Exclusion of CTL by immunosuppressive factors (Exclusion), and M2 subtypes of M2 (TAM), which are the three types of cells that restrict T cell invasion into tumors. 2.6. Drug sensitivity analysis Sensitivities of Erlotinib, MG-132, Z-LLNle-CHO, Dasatinib, CGP-60474 and WH-4-023 to IC50 were predicted using pRRophetic [[81]26]. 2.7. Gene set enrichment analysis (GSEA) GSEA was conducted on the basis of the H.all.v7.5.1.Entrez.GMT gene set in the MSigDB database [[82]27], with FDR<0.05 indicating significant enrichment. 2.8. Development of a prognosis model for osteosarcoma To filter DEGs and prognostic genes from the molecular subtypes classified, the limma package and univariable Cox analysis were employed. Lasso regression was conducted using the glmnet package [[83]28] to reduce the range of genes. The following formula was calculated to assess the prognosis of osteosarcoma: [MATH: RiskSco re=∑k=0n βi×Expi :MATH] where βi and Expi refer to the Cox regression coefficient and expression of the i gene, respectively. The median RiskScore value was used to divide low- and high-risk groups of the samples in the training and independent validation datasets. The prognostic prediction was further assessed using KM survival curve and ROC. 2.9. MCP-counter and CIBERSORT algorithm The 8 immune populations and 2 stromal populations were quantified using MCP-counter for all the samples [[84]29]. The CIRERSORT algorithm was run to assess the landscape of immune infiltration of 22 immune cells in different risk groups [[85]30]. 2.10. Nomogram development The model was further evaluated by developing a nomogram with other clinicopathological features of patients with osteosarcoma (e.g. age, gender, metastasis) using RMS software R package. The decision curve was plotted to test the nomogram prediction. 2.11. Cell culture and transfection The osteosarcoma cell lines hFOB1.19 (BNCC338626) and Saos-2 (BNCC338485) were purchased from BNCC (Beijing). Dulbecco's Modified Eagle Medium (Gibco, 11,965–092) added with 10 % fetal bovine serum (Gibco, 26,140–095) and 1 % antibiotics (Gibco, 15,070–063) were used for cell culture at 37 °C in 5 % CO[2]. Next, cell transfection with negative control (NC) and IFITM3 siRNA (Sagon Co., China) was conducted with the use of Lipofectamine 2000 (Invitrogen, USA). 2.12. QRT-PCR TRIzol reagent was used for total RNA extraction (Invitrogen, Carlsbad, California, USA) from cells. RNA was reverse-transcribed into cDNA using the Qiagen One-Step RT-PCR kit (Qiagen Gmbh, Germany) and subjected to qRT-PCR experiments. Amplification was performed using SYBR Green (Bimake, Houston, Texas, USA) on the ABI 7500 system (Thermo Fisher Scientific, USA). The primers for qRT-PCR are shown in [86]Supplementary Table 2. 2.13. Wound-healing assay Osteosarcoma Saos-2 cellss were seeded into 6-well plate to grow until they covered the entire bottom. A scratch was made vertically using a 200 μL pipette tip. After washing the cells twice with phosphate-buffered saline, the cells were photographed under an inverted microscope at 0 h and 48 h after the scratching. The wound healing rate was calculated as [(Width at 0 h - Width at 48 h)/Width at 0 h] × 100 %. The experiment was performed in triplicate. 2.14. Transwell analysis Saos-2 cells were seeded into the upper chamber of the Transwell apparatus (Corning, USA) containing serum-free medium. Subsequently, the lower chamber was supplemented with PRMI-1640 medium with 10 % FBS. After 1-day incubation at 37 °C, cells remaining on the upper chamber were removed, whereas the rest cells were fixed by 4 % paraformaldehyde and subsequently dyed by crystal violet for 30 min. A microscope was employed to observe the cells and count the cell number. 3. Results 3.1. Three molecular subtypes were defined using the necroptosis genes There were a total of 116 genes associated with necroptosis score in TARGET-osteosarcoma dataset. Next, 39 prognostic genes (34 protective genes and 5 risk genes) for osteosarcoma were screened by univariate cox survival analysis. Using the 39 genes, ConsensusClusterPlus was employed to group patients in TARGET-osteosarcoma dataset. As shown in the relative change in area under CDF curve ([87]Fig. 1A and B), C1, C2 and C3 were clustered when K = 3 ([88]Fig. 1C). Specifically, the prognosis were the worst in C3 in TARGET-osteosarcoma ([89]Fig. 1D) and [90]GSE39058 dataset ([91]Fig. 1E), whereas the samples in C1 and C2 had better survival. Risk genes had high expressions in C3, while protective genes were high-expressed in C1 and C2 ([92]Fig. 1F). Fig. 1. [93]Fig. 1 [94]Open in a new tab Classification of 3 molecular subtypes. A: Cumulative distribution function. B: Delta area of Cumulative distribution function. C: when k = 3, 3 molecular subtypes were identified. D: KM survival curve of 3 molecular subtypes in TARGET-osteosarcoma. E: KM survival curve of 3 molecular subtypes in [95]GSE39058 dataset. F: Heatmap of necroptosis genes expressions in 3 molecular subtypes. The distribution of clinical characteristics in C1, C2 and C3 was shown in [96]Fig. 2, which demonstrated that patients’ status in TARGET-osteosarcoma cohort was noticeably different. Fig. 2. [97]Fig. 2 [98]Open in a new tab Clinical features distribution in 3 molecular subtypes. 3.2. Immune microenvironment and immunotherapy analysis for the three molecular subtypes Immune cell infiltration in the TARGET-osteosarcoma cohort was assessed based on the expressions of immune cells. Compared to the other two subtypes, C1 had the highest StromalScore, ImmuneScore and ESTIMATEScore ([99]Fig. 3A) and the highest infiltration of 22 out of 28 immune cells ([100]Fig. 3B). We also observed that most of the immune checkpoint genes were high-expressed in C1 than in C2 and C3 ([101]Fig. 3C). The results of TIDE showed higher Exclusion, TIDE, and MDSC in C2 group and higher IFNG and Dysfunctions in C1. In addition, C3 had more tumor-associated macrophage M2 subtypes than in C1 and C2 ([102]Fig. 3D). Moreover, the IC50 of Erlotinib, MG-132, Z-LLNle-CHO, Dasatinib, CGP-60474 and WH-4-023 was higher in C3, suggesting a lower sensitivity of C3 patients to these drugs ([103]Fig. 3E). Fig. 3. [104]Fig. 3 [105]Open in a new tab Immune microenvironment and immunotherapy analysis. A: Differences of the 3 molecular subtypes in ESTIMATEScore, StromalScore, and ImmuneScore. B: Differences of 28 kinds of immunes score in the 3 clusters. C: Differences of the expression of immune checkpoint genes in the 3 clusters. D: TIDE analysis among C1, C2, C3. E: The estimated IC50 for drug in TARGET-osteosarcoma was shown in the box plots. (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001). 3.3. Pathways analysis for the three molecular subtypes The GSEA showed that compared with genes C3, genes in C1 in the TARGET cohort were significantly enriched to 21 pathways, including epithelial mesenchymal transition (EMT), IFN-α response, interferon-γ response, etc. ([106]Fig. 4A). Subsequently, we also compared the differential pathways between the three subtypes (including C1vsC2, C1vsC3, C2vsC3). Among the differential pathways between C2 and C3, the enrichment scores in EMT and UV response DN pathways were the highest ([107]Fig. 4B). As shown in [108]Fig. 4C, the immune-related pathways (including interferon response and IL6-JAK-STAT3 signaling pathway) in patients in C1 subtype were activated but cell cycle-related pathways were inhibited. Fig. 4. [109]Fig. 4 [110]Open in a new tab Pathways enrichment analysis. A: GSEA analysis of C1vsC3 in TARGET-osteosarcoma. B: GSEA analysis of C1vsC3, C1vs C2, C2vs C3 in TARGET-osteosarcoma. C: Radar plot of activated pathways in C1vsC2, C1vsC3, and C2vsC3 in the TARGET-osteosarcoma cohort. 3.4. Establishment of a prognosis risk model A total of 114 prognosis genes (48 risk genes and 66 protective genes) were selected out of 983 DEGs from C1, C2, and C3 ([111]Fig. 5A). The trajectory and confidence interval of lambda in Lasso analysis was presented ([112]Fig. 5B and C). When lambda = 0.1547, a risk model was developed using 8 genes: RiskScore=(−0.516*IFITM3)+0.068*ACTA2+(−0.108*GBP1)+(−0.355*APBB1IP)+(− 0.26*GJA5)+0.116*CGREF1+(−0.309*CDK6)+0.26*TAC4. Fig. 5. [113]Fig. 5 [114]Open in a new tab Identification of hub genes. A: 114 genes associated with prognosis were identified. B: As lambda changes, the trajectory of the independent variable was shown. C: Confidence interval under lambda. D: The distribution of Lasso coefficients of genes related to in the necroptosis-based signature. 3.5. Validation of the prognosis risk model The risk score for all the patients in TARGET-osteosarcoma dataset was calculated with the formula ([115]Fig. 6A). Subsequently, osteosarcoma patients were categorized into two risk groups (high and low risks) by the risk score based on the threshold of “0". In the TARGET-osteosarcoma cohort, the AUC for 1-year, 3-year and 5-year survival was 0.8, 0.9 and 0.88, respectively ([116]Fig. 6B). The KM survival curve showed a better OS in low-risk group in TARGET-osteosarcoma dataset ([117]Fig. 6C). In [118]GSE39058 cohort, the AUC for 1-, 3- and 5-year survival was 0.8, 0.86 and 0.9, respectively ([119]Fig. 6D), with the low-risk samples having better survival ([120]Fig. 6E). In [121]GSE21257 cohort, the AUC for 1-, 3- and 5-year survival was 0.73, 0.7 and 0.73, respectively ([122]Fig. 6F), with the low-risk samples having longer survival ([123]Fig. 6G). Fig. 6. [124]Fig. 6 [125]Open in a new tab Validation of the risk score. A: The risk score, survival status, and gene expression in TARGET-osteosarcoma dataset. B: ROC analysis of the risk score in TARGET-osteosarcoma. C: KM survival curve between low group and high group in TARGET-osteosarcoma. D: ROC analysis of the risk score in [126]GSE39058 dataset. E: KM survival curve between the two groups in [127]GSE39058 dataset. F: ROC analysis of the risk score in [128]GSE21257 dataset. G: KM survival curve between the two groups in [129]GSE21257 dataset. 3.6. Effects of the model feature genes on the malignant phenotype of osteosarcoma cell lines Cellular experiments were used to analyze the influences of the model genes on osteosarcoma progression of and to detect key genes involved in the malignant phenotype of osteosarcoma. Molecular detection results showed that IFITM3, APBB1IP, CDK6, GBP1, TAC4, and GJA5 expressions were relatively upregulated, and GJA5, CGREF1 mRNA levels were downregulated in osteosarcoma cell lines, with IFITM3 having the highest expression level ([130]Fig. 7A–H). Wound healing and Transwell experiment demonstrated that silencing IFITM3 gene suppressed the migration ([131]Fig. 7I–J) and invasion ([132]Fig. 7K–L) of osteosarcoma cells, respectively. Overall, IFITM3 is a gene significantly affecting the prognostic outcomes of patients suffering from osteosarcoma and had a tangible regulatory effect on osteosarcoma cell lines, which verified the reliability of the prognostic model. Fig. 7. [133]Fig. 7 [134]Open in a new tab Effects of model feature genes on osteosarcoma cell. A–H: qRT-PCR results of model feature genes IFITM3, APBB1IP, CDK6, GBP1, TAC4, and GJA5. I–J: Scratch healing experiments and the quantified results of osteosarcoma cell line of si-IFITM3. K–L: Transwell experiment and the quantified results of osteosarcoma cell line of si-IFITM3. 3.7. Distribution of the risk score in clinical characteristics The risk score for samples in different clinical characteristics was computed, including metaststic, gender, age, cluster and status. The samples in metaststic, C3 and Dead had a higher risk score ([135]Fig. 8A), while those in gender, age, metaststic had low risk score and better survival ([136]Fig. 8B). Fig. 8. [137]Fig. 8 [138]Open in a new tab Clinical features analysis for both the high and low groups. A: The differences in the risk score between different clinicopathological groups in the TARGET-osteosarcoma cohort. B: KM curves for both the high-risk and low-risk groups among different clinicopathological groups in the TARGET-osteosarcoma cohort. 3.8. Analysis on immune microenvironment between the two risk score groups Immune microenvironment differences between the two groups were observed, with the group of low risk score having higher ESTIMATEScore, StromalScore, and ImmuneScore ([139]Fig. 9A). MCP-Counter analysis demonstrated that the scores of 6 out of 10 types of immune cells in the low group were higher ([140]Fig. 9B). The ssGSEA analysis showed that the scores of 15 of 28 types of immune cell scores were higher in the low group ([141]Fig. 9C). Additionally, we calculated the infiltration of 22 immune cells in different risk groups using the CIBERSORT and obtained similar results to the MCP-Counter algorithm. When compared to the high-risk group, low-risk patients had higher CD8+T cells infiltration ([142]Supplementary Fig. 1). Further study revealed that most immune cells from 28 immune cells showed negative correlation with the risk score ([143]Fig. 9D). Fig. 9. [144]Fig. 9 [145]Open in a new tab Immune microenvironment analysis. A: The differences of the high group and low group in ESTIMATEScore, StromalScore, ImmuneScore. B: The differences of the two groups in 10 kinds of immunes score. C: The differences of the two groups in 28 kinds of immunes score. D: The correlation between 28 immune cells scores and the risk score. (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001). 3.9. Pathways analysis for the two risk score groups The ssGSEA analysis demonstrated a negative correlation between the immune pathways and risk score ([146]Fig. 10A). GSEA analysis revealed that 12 pathways and 30 pathways were activated in the group with a high risk score in TARGET-osteosarcoma and [147]GSE39058 dataset, respectively, and that most pathways were associated with cell cycle ([148]Fig. 10B). Fig. 10. [149]Fig. 10 [150]Open in a new tab Pathways enrichment analysis. A: Results of correlation analysis between the risk score and KEGG pathways. NPRS showed a correlation higher than 0.4. B: By comparing high and low group, normalized enrichment scores of Hallmark pathways were shown in the heatmap. 3.10. Clinical pathology features and the risk score to predict osteosarcoma survival The decision tree was produced based on gender, age, metastatic and the risk score in the TARGET-osteosarcoma cohort. Five different risk subgroups ([151]Fig. 11A) with significant differences in OS ([152]Fig. 11B) were classified. Patients in C1, C2 and C3 were in low risk score group, while those in C4 and C5 were in high risk score group ([153]Fig. 11C). The 5 subgroups had significance differences in survival status ([154]Fig. 11D). Both metastatic and the risk score were validated as independent prognosis factors by univariate and multivariate Cox survival analysis ([155]Fig. 11E and F). We also constructed a nomogram together with the risk score and metastatic to predict 1-, 3-, 5-year OS of osteosarcoma patients ([156]Fig. 11G). The risk score had higher AUC ([157]Fig. 10H) and the reliability of the nomogram was proven by the calibration curve ([158]Fig. 11I). DCA further validated the risk score as an indicator for effectively predicting osteosarcoma prognosis in comparison to other clinical variables ([159]Fig. 11J.). Fig. 11. [160]Fig. 11 [161]Open in a new tab Construction of nomogram. A: To optimize risk stratification, a survival decision tree was established according to the full-scale annotations of patients including age, metastatic, NPRS. B: The overall survival outcomes of the 5 risk subgroups were significantly different. C: The distribution of low and high group samples in the 5 risk subgroups. D: The status distribution in 5 risk subgroups. E: Univariate Cox analysis. F: Multivariate Cox analysis. G: The development of a nomogram. H: The survival prediction of the nomogram was the most accurate in comparison to other clinicopathological features. I: 1-, 3- and 5- years calibration curves for the nomogram. J: Decision curve of the nomogram. 4. Discussion Previous studies have confirmed crucial role of necroptosis in cancer development and metastasis [[162]31,[163]32]. In this work, 39 necroptosis-related genes in TARGET-osteosarcoma were screened to classify C1, C2, and C3 as three molecular subtypes of osteosarcoma. Patients in C1 and C2 had enhanced immune status and better survival but limited immunotherapy benefit. The prognostic outcomes of patients with osteosarcoma could be correctly evaluated by the 8 necroptosis-correlated prognostic signature genes. The prognostic outcomes were better for patients who had a low risk score, and the nomogram combining both the risk score and risk factors could produce an accurate prognostic prediction. A total of 8 DEGs (ACTA2, APBB1IP, CDK6, CGREF1, GBP1, GJA5, IFITM3, and TAC4) were contained in the present model. The ROC curve and calibration curve showed that the model had an accurate prediction for osteosarcoma patients in the training and validation cohorts. Among all the 8 genes related to osteosarcoma prognosis, 4 genes (ACTA2, CDK6, IFITM3) were either involved in necroptosis-related pathogenesis of osteosarcoma or were important indicators of overall survival or relapse-free survival [[164][33], [165][34], [166][35]]. As an interferon-induced transmembrane protein, IFITM3 is high-expressed in the tumor region of gastric cancer. Depletion of IFITM3 can cause HGF-triggered inhibition of cell growth and migration by suppressing AKT/c-MYC signaling in gastric cancer [[167]36]. Gan et al. found that knockdown of IFITM3 could impair the growth of the oral squamous cell carcinoma cell lines by inhibiting the proliferation of the cancer cells and inducing cell cycle arrest, senescence and apoptosis [[168]37]. Consistently, we observed that silencing IFITM3 expression greatly affected the osteosarcoma cell migration and invasion. In addition, Tang et al. observed that patients with high-expressed ACTA2 have a noticeably more favorable prognosis than those with low-expressed ACTA2, suggesting that high-expressed ACTA2 in osteosarcoma patients with lung metastasis indicates a better prognosis [[169]38]. Recent discoveries showed that CDK6 has a crucial function in the advancement of various human cancers, and that a downregulation in its expression is associated with a negative prognostic outcome. Based on an analysis of the GEPIA and STARBASE databases, Zhao et al. found that CDK6 is an oncogene in pancreatic cancer, and that CDK6 is significantly associated with PD-L1, PD-L2, and HAVCR2 (three immune checkpoints), immune cells infiltration, and immune biomarkers [[170]39]. However, previous study was less concerned with the remaining 5 model genes in the osteosarcoma prognosis or their roles as new potential biomarkers for the cancer. Therefore, these key gene targets should be analyzed more extensively. Past study has analyzed the roles of TME, in particular immune microenvironment [[171]40,[172]41]. Immune infiltration analysis demonstrated a negative relationship between the immune cell score and the risk score. Compared to the high-risk patients, those with a low risk had more immune cells and greater activation of immune pathways. Studies found that tumors could induce the generation of type II NKT cells and secrete IL-13 in turn, which causes MDSC aggregation in the tumor microenvironment, activates STAT6 signaling pathway and suppresses the function of CD8^+ T cell [[173]42,[174]43]. Based on scRNA-seq transcriptomic data and immunohistochemical staining analysis, Xie et al. confirmed that EVI2B is present on CD8^+ T cells derived from osteosarcoma patients. EVI2B enhances the expressions of granzyme A and K in CD8^+ T cells, resulting in a robust cytotoxic effect against tumor cells [[175]44]. NK cells within the TME of osteosarcoma are innate immune cells that display cytokine-secreting and cytotoxic properties along with cytotoxic and helper T cells [[176][45], [177][46], [178][47]]. Both CD8^+ T cells and NK cells possess the ability to eliminate cancer cells through their cytotoxic activities, yet the synergistic relationship between them is intricate. MHC-I expression downregulated by specific cancer cells escape the detection of CD8^+ T cells while enhancing NK cell activation by eliminating a significant inhibitory cue [[179]48]. At present, the application of T cells and NK cells in tumor immunotherapy has attracted extensive research attention. Apart from hematological tumors, memory-like NK cells and modified NK cells have great potential in treating breast cancer, colorectal cancer, ovarian cancer, non-small cell lung cancer, liver cancer [[180]49]. This suggested that the tumor microenvironment of low-risk patients may be more conducive to immune cell survival and function, which was strongly associated with a better prognosis for patients. In conclusion, differentially expressed necroptosis-related genes were identified in this study to develop a prediction model for osteosarcoma prognosis, and the association between immune activities and the risk score was systematically assessed. However, the mechanisms of necroptosis in tumor immunology should be further investigated by future studies. Ethical approval and consent to participate Not applicable. Consent for publication Not applicable. Ethical guidelines Not applicable. Availability of data and materials The datasets generated and/or analyzed during the current study are available in the [[181]GSE16088] repository, [[182]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE16088]; in [[183]GSE12865] repository, [[184]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12865]; in [[185]GSE21257] repository, [[186]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE21257]; in [[187]GSE39058] repository, [[188]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39058]. Funding National Natural Science Foundation of China (No.81772864, No.81772861) supported this project. CRediT authorship contribution statement Yiying Bian: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Investigation, Data curation, Conceptualization. Jixiang Shi: Validation, Software, Resources, Methodology, Investigation. Ziyun Chen: Supervision, Software, Methodology, Investigation, Formal analysis, Data curation. Ji Fang: Visualization, Supervision, Resources, Formal analysis. Weidong Chen: Visualization, Supervision, Software, Project administration, Methodology. Yutong Zou: Visualization, Resources, Investigation, Formal analysis, Conceptualization. Hao Yao: Visualization, Software, Project administration, Methodology, Investigation. Jian Tu: Visualization, Validation, Supervision, Formal analysis, Data curation, Conceptualization. Yan Liao: Visualization, Supervision, Resources, Project administration, Investigation. Xianbiao Xie: Writing – original draft, Visualization, Validation, Software, Funding acquisition, Formal analysis. Jingnan Shen: Writing – review & editing, Visualization, Supervision, Project administration, Funding acquisition, Formal analysis. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements