A novel stratification framework based on anoikis-related genes for predicting the prognosis in patients with osteosarcoma

Background Anoikis resistance is a prerequisite for the successful development of osteosarcoma (OS) metastases, whether the expression of anoikis-related genes (ARGs) correlates with OS prognosis remains unclear. This study aimed to investigate the feasibility of using ARGs as prognostic tools for the risk stratification of OS. Methods The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases provided transcriptome information relevant to OS. The GeneCards database was used to identify ARGs. Differentially expressed ARGs (DEARGs) were identified by overlapping ARGs with common differentially expressed genes (DEGs) between OS and normal samples from the GSE16088, GSE19276, and GSE99671 datasets. Anoikis-related clusters of patients were obtained by consistent clustering, and gene set variation analysis (GSVA) of the different clusters was completed. Next, a risk model was created using Cox regression analyses. Risk scores and clinical features were assessed for independent prognostic values, and a nomogram model was constructed. Subsequently, a functional enrichment analysis of the high- and low-risk groups was performed. In addition, the immunological characteristics of OS samples were compared between the high- and low-risk groups, and their sensitivity to therapeutic agents was explored. Results Seven DEARGs between OS and normal samples were obtained by intersecting 501 ARGs with 68 common DEGs. BNIP3 and CXCL12 were significantly differentially expressed between both clusters (P<0.05) and were identified as prognosis-related genes. The risk model showed that the risk score and tumor metastasis were independent prognostic factors of patients with OS. A nomogram combining risk score and tumor metastasis effectively predicted the prognosis. In addition, patients in the high-risk group had low immune scores and high tumor purity. The levels of immune cell infiltration, expression of human leukocyte antigen (HLA) genes, immune response gene sets, and immune checkpoints were lower in the high-risk group than those in the low-risk group. The low-risk group was sensitive to the immune checkpoint PD-1 inhibitor, and the high-risk group exhibited lower inhibitory concentration values by 50% for 24 drugs, including AG.014699, AMG.706, and AZD6482. Conclusion The prognostic stratification framework of patients with OS based on ARGs, such as BNIP3 and CXCL12, may lead to more efficient clinical management.


Introduction
Osteosarcoma (OS) is the most common primary bone malignancy that seriously threatens the health of children and adolescence worldwide (1). With the widespread use of multiagent chemotherapy in the 1980s, the 5-year survival rate of patients with OS without metastasis has substantially increased from 20% to 65% (2, 3). However, the prognosis of patients with recurrence and metastasis remains dismal, with only a 12% 4month event-free survival rate in single-arm phase II clinical trials (4). Owing to the heterogeneity of the tumor itself and the complex mechanism of metastasis, metastatic cells vary phenotypically from primary tumor cells and may more easily escape immune surveillance and survive chemotherapy (5,6). This may explain why patients with the same clinical or pathological conditions, who receive the same treatment regimen, may have different clinical outcomes. Additionally, the lack of reliable markers and therapeutic targets makes it difficult to identify patients who are at a high risk of metastasis and may maximally benefit from specific therapies. Therefore, for a considerable period, the prediction and risk stratification of distant metastases before treatment remain the cornerstone for therapeutic decisions.
Anoikis is a type of programmed cell death caused by cell detachment from the extracellular matrix (ECM) (7). As a protective mechanism, anoikis can prevent the ectopic growth of somatic cells at inappropriate body sites. However, tumor cells with malignant potential usually show resistance to anoikis, allowing them to escape from the primary tumor site and colonize secondary sites (8). Specific factors, such as the integrin family, growth factor family, and metabolic intermediates have been reported as drivers of anoikis resistance, thus enhancing tumor recurrence and metastasis (7,9). With the rapid development of omics technology and bioinformatic approaches, accumulating evidence suggests that tumor metastasis is tightly regulated at multiple levels by anoikis-related genes (ARGs). These genes remodel tumor cells by reprogramming lipid and amino acid metabolism to promote anoikis resistance. For example, the GDH1mediated metabolic reprogramming of glutaminolysis mediates lung cancer metastasis (10). Upon TGF-b2 stimulation, PKC-zeta-mediated translocation of CD36 facilitates the uptake of fatty acids by tumor cells and supports their invasiveness (11). Advances in research on ARGs have provided new opportunities to clarify the mechanisms of tumor metastasis. However, to the best of our knowledge, few studies are available on the role of ARGs in OS metastasis. However, whether abnormal expression of ARGs is associated with poor OS prognosis has not been fully explored (5). The identification and characterization of ARGs in OS may provide new biomarkers and therapeutic targets for diseases.
Considering the potential correlations among ARGs, OS metastasis, and clinical outcomes, this study focused on the identification of ARGs that may influence the prognosis of OS. Following prognosis-associated clinical parameter screening, a prognostic risk model was established, and its ability to identify high-risk populations with OS was evaluated. Finally, we explored the specific immunological features of the high-and low-risk patients. The results of this study provide valuable references for clinicians to stratify at-risk patients with OS, which may lead to more efficient clinical management.

Identification of differentially expressed ARGs among OS and normal samples
Using the limma R package (Version 3.50.1), differentially expressed genes (DEGs) in samples from the OS and normal sample groups were obtained from the GSE16088 and GSE19276 datasets with the criteria of |logFC| > 1 and adj.P.Value < 0.05. In the GSE99671 dataset, DEGs among the OS and normal sample groups were screened using the DESeq2 R package (Version 1.34.0), and the identification thresholds were |logFC| > 1 and adj.P.Value < 0.05. The upregulated DEGs from the GSE16088, GSE19276, and GSE99671 datasets were crossed to obtain common upregulated DEGs. Similarly, common down-regulated DEGs in the above three datasets were obtained. Common upregulated and downregulated DEGs were combined to obtain DEGs between OS and normal samples. The combined DEGs were intersected with ARGs to obtain DEARGs.

Identification of anoikis-related clusters in patients with OS
Consistent clustering analysis was performed on samples from the TARGET-OS dataset based on DEARGs using the R package ConsensusClusterPlus (version 1.58.0) (12). The consistency clustering effect was assessed using principal component analysis (PCA). Next, the clinical characteristics of the different clusters of patients with OS were analyzed using the Chi-square test. In addition, differences in functional enrichment between different clusters were compared using the gene set variation analysis (GSVA) R (version 1.42.0) (13).

Risk model of patients with OS
Univariate Cox analysis of DEARGs in the TARGET-OS dataset was used to identify prognosis-related genes of patients with OS. Prognosis-related genes were tested for the proportional hazards (PH) hypothesis, and a multifactorial Cox model was constructed (14). Based on the median risk score, the patients with OS were classified into high-and low-risk groups. Kaplan-Meier (K-M) survival curves were used to compare the survival probability of patients with OS between the low-and high-expression groups based on the median expression value of each prognostic-related gene (15).
In the TARGET-OS dataset, differences in survival between patients in the two risk groups were compared using K-M survival curves. The receiver operating characteristic (ROC) curve of the risk score was plotted using the survivalROC R package (Version 3.2-13) (16). In addition, differences in the clinical traits of patients with OS among the risk groups were explored using the TARGET-OS dataset, and the association between the clinical characteristics and risk scores of patients with OS was assessed.

Nomogram model and GSVA between the high-and low-risk groups
In the TARGET-OS dataset, the independent prognostic value of clinical characteristics and risk scores was explored using univariate and multifactorial Cox analyses. Then, a nomogram model was obtained by combining factors with independent prognostic values, and the reliability of the nomogram model was assessed using decision curve analysis (DCA), ROC curves, and calibration curves (16). The GSVA was performed on samples from both risk groups.

Immune infiltration landscape
In the TARGET-OS dataset, the estimate R package (version 1.0.13) was used to compare stromal scores, immune scores, ESTIMATE scores, and tumor purity in samples from the highand low-risk groups. The infiltration level of 30 tumor microenvironment (TME) cells in samples from both groups was assessed using single-sample gene set enrichment analysis (ssGSEA) (17). Correlations between differential immune cells and risk scores were calculated using Spearman analysis. In addition, the survival probability of patients with OS in the two groups of each TME cell type, according to the median ssGSEA score, was evaluated using K-M survival analysis. TME cells with different infiltrations between the two groups, TME cells that were related to the survivability of patients with OS, and TME cells related to the risk score were intersected to obtain important TME cells for the risk scores of patients with OS, as described above.

Immune microenvironment
First, immune-related genes were acquired from the ImmPort database (https://www.immport.org), and their expression in samples from the high-and low-risk groups was assessed using the TARGET-OS dataset. Correlations between prognostic-and immune-related genes were assessed using Pearson correlation analysis. Differences in the expression of human leukocyte antigen (HLA) genes between samples from both groups were compared, and correlations between prognostic-related genes and HLA genes were assessed. In addition, differences in the expression of immune checkpoints between both groups were evaluated, and the correlations between immune checkpoints and risk scores were further explored.

Sensitivity of patients with OS to therapeutic drugs
Treatment sensitivity to PD-1 and CTLA4 inhibitors was predicted in patients with OS from high-and low-risk groups using the SubMap algorithm in the TARGET-OS dataset (18). Differences in the 50% inhibitory concentration (IC50) of 138 drugs between both groups were compared using the pRRophetic algorithm (19).

Screening for DEARGs among OS and normal samples
Through differential analysis, 5345 DEGs among OS and normal samples were obtained from the GSE16088 dataset, comprising 2905 up-regulated genes (URGs) and 2440 downregulated genes (DRGs). A total of 902 DEGs between OS and normal samples were obtained from the GSE99671 dataset, comprising 265 URGs and 637 DRGs. A total of 1248 DEGs between OS and normal samples were screened in the GSE19276 dataset, comprising 501 URGs and 747 DRGs ( Figure 1A).
Genes screened from the three datasets with upregulated expression were crossed to acquire 11 commonly upregulated DEGs ( Figure 1B). Similarly, genes that were downregulated in the three datasets were crossed to acquire 57 common downregulated DEGs ( Figure 1C). Combining the common URGs and DRGs yielded 68 DEGs between the OS and normal samples.

Acquisition of anoikis-related clusters of patients with OS
According to the seven DEARGs, consistent cluster analysis of the TARGET-OS dataset samples showed that clustering was the best when K=2 and all samples were classified into cluster 1 and cluster 2 (Figures 2A-C). PCA showed good results for consistent clustering, with all samples clustered into two clusters ( Figure 2D). The race of patients with OS in the different clusters was significantly different (P<0.05); however, no significant difference was detected in the other clinical characteristics between both clusters ( Figure 2E). In addition, 27 pathways were enriched between both clusters by GSVA, including progesterone-mediated oocyte maturation, oocyte meiosis, cell cycle, folate biosynthesis, and steroid biosynthesis ( Figure 2F).

Prognosis-related gene expression levels were significantly different in OS and normal samples
In the TARGET-OS dataset, two prognosis-related genes (CXCL12 and BNIP3) were identified by univariate Cox analysis of the seven DEARGs (P<0.05) ( Figure 3A). Both CXCL12 and BNIP3 satisfied the PH hypothesis (P<0.05) and were used to construct a multifactorial Cox model ( Figure 3B). Patients in the high-risk group with relatively short survival times accounted for most patients with OS ( Figures 3C, D). CXCL12 was strongly expressed in the low-risk group, whereas BNIP3 was substantially expressed in the other group ( Figure 3E). In the GSE16088, GSE19276, and GSE99671 datasets, the expression level of BNIP3 was higher and that of CXCL12 was lower in the OS samples than those in the normal samples (adj.P < 0.05) (Figures 3F-H).
K-M survival curves demonstrated that, compared with that of the other groups, the probability of survival in patients from the BNIP3 high-expression group in the TARGET-OS dataset was greatly reduced (P<0.05) ( Figure 3I). In contrast, the probability of survival in patients from the CXCL12 high expression group was significantly higher (P<0.05) ( Figure 3J).

Survival probabilities for patients in the high-risk group were substantially lower compared to those in the low-risk group
In the TARGET-OS dataset, the probability of survival was greatly reduced for patients in the high-risk group compared to that in the other groups (P<0.05) ( Figure 4A). The area under the ROC curve (AUC) of the risk model was 0.7, indicating that the risk model had a strong predictive ability for TARGET-OS dataset patients ( Figure 4B).
In the validation dataset (GSE16091), high-risk patients with relatively short survival times accounted for most patients with OS ( Figures 4C, D). The expression patterns of CXCL12 and BNIP3 in the high-and low-risk groups were similar to those observed in the TARGET-OS dataset ( Figure 4E). The probability of survival was greatly reduced in patients in the high-risk group compared to those in the other groups (P<0.05) ( Figure 4F). Meanwhile, the risk model correctly estimated the prognosis of patients with OS because the ROC curves of the risk model in the GSE16091 dataset were all greater than 0.6 ( Figure 4G).

OS patients of Black or African American and metastatic OS patients had a high-risk score
In the TARGET-OS dataset, the race of the patients with OS differed between the high-and low-risk groups, whereas other clinical characteristics did not significantly differ (Table 1). Age and sex did not strongly correlate with the risk scores of patients with OS ( Figures 5A, B). Black or African Americans had significantly higher risk scores than those Asian patients ( Figure 5C). Patients with metastatic OS had substantially higher risk scores than those with nonmetastatic OS ( Figure 5D).

Nomogram model could reliably forecast the prognosis of patients with OS
Univariate Cox analysis confirmed the association between the risk score and metastasis and OS in the TARGET-OS dataset (P<0.05) ( Figure 6A). Next, the independent prognostic value of the risk score and metastasis was identified using multifactorial Cox analysis (P<0.05) ( Figure 6B).
To create a nomogram model, it was possible to reliably estimate OS prognosis by fusing the risk score with metastasis ( Figures 6C-E). In addition, the DCA curves showed that the nomogram model had greater reliability in predicting the survival of patients with OS ( Figure 6F).

3.7
The differentially expressed signaling pathways among high-and low-risk groups were mainly associated with tumorigenesis and metastasis The GSVA enriched 31 pathways that significantly differed between both risk groups, including IGA production facilitated by the intestinal immune network, thyroid disease of autoimmune origin, rejection of allografts, leishmaniasis, and cellular adhesion molecules (Figure 7).

High-risk score was characterized by low immune infiltration and elevated tumor purity
In the TARGET-OS dataset, we detected that samples in the highrisk group had significantly lower stromal, immune, and ESTIMATE scores and significantly higher tumor purity than those of the samples in the other groups (P<0.05) ( Figure 8A). The level of infiltration of 23 TME cells was significantly different from that of the other groups (P<0.05) ( Figure 8B), and 24 TME cells were significantly correlated with risk scores (P<0.05) ( Figure 8C). Individuals in the low-scoring groups of central memory CD8 T cells, activated B cells, macrophages, monocytes, effector memory CD8 T cells, CD56 bright natural killer cells, and natural killer T cells exhibited dramatically lower survival probabilities than those in the other groups (P<0.05) (Figures 8D-J).
The 23 TME cells that were differentially infiltrated between the two risk groups, the 24 TME cells that were related to the risk score, and the seven cells that were related to survival in patients with OS were intersected to obtain seven TME cells that were significant in the OS risk score ( Figure 8K).

3.9
The expression of the immune response gene set and HLA genes immune checkpoints were significantly lower in the high-risk group of patients with OS Eight immune response gene sets differed between both risk groups in the TARGET-OS dataset ( Figure 9A). Immune response gene sets were negatively correlated with BNIP3 expression and positively correlated with CXCL12 expression ( Figure 9B). The immune response gene sets most associated with CXCL12 and BNIP3 were both highly expressed in the low-risk group compared to those in the other groups ( Figures 9C-F).
In the TARGET-OS dataset, the expression of 17 HLAs was highly variable between both groups (P<0.05) ( Figure 10A). The HLAs were negatively associated with BNIP3 and actively associated with CXCL12 ( Figure 10B). The HLAs were most strongly associated with CXCL12 and BNIP3, and both were highly expressed in the low-risk group compared to those in the other groups ( Figures 10C-F). In addition, the expression of CD274, CTLA4, and HAVCR2 was noticeably lower in the high-risk group than that in the other groups (P<0.05) ( Figure 11A). The risk score was negatively correlated with HAVCR2, LAG3, PDCD1LG2, CD274, and CTLA4 ( Figure 11B).

Patients with OS in the high-risk group had reduced sensitivity to therapeutic agents
Patients in the low-risk group were susceptible to PD-1 inhibitor treatment (P<0.05) ( Figure 12A). A total of 81 drugs showed significantly different IC50 values between the two groups (P<0.05). Of these, 24 drugs, including AG.014699, AMG.706, AZD6482, and BI.D1870, may be candidates for treating patients in the high-risk group ( Figure 12B). In contrast, 57 drugs, including A.443654, A.770041, AKT inhibitor VIII, and AP.24534, may not be ideal for patients in the other group ( Figure 12C).

Discussion
Although the prognosis of localized OS has markedly improved owing to new therapeutic developments, long-term survival has stagnated over the past several decades (2-4). Metastasis, particularly lung colonization, is the most common cause of death in high-risk patients with OS (20). Anoikis is a physiological process that plays an important role in tissue homeostasis and development. Under pathological conditions, it is the main factor in tumor metastasis and therapy failure (7,8). Alterations in ARGs that lead to anoikis resistance are hallmarks of the malignant transformation of tumors (10,11,21,22). Regrettably, the relationship between ARGs and OS progression is much less recognized, which has limited the improvement in patient prognosis. The present study identified prognosis-related genes involved in OS progression. Based on BNIP3 and CXCL12 expression, this stratification framework can be used to effectively stratify the survival of patients with OS. Moreover, patients with high-risk OS presented unique patterns of immune characteristics and sensitivity to different chemotherapeutic agents. These findings provide a scientific basis for the discovery of new immunotherapeutic targets and efficient selection of existing drugs in clinical practice.
In this study, CEACAM1, CXCL12, and LTF were found to be downregulated in OS samples. CEACAM1 is a transmembrane cell adhesion molecule belonging to the CEA superfamily (23). Similar to the inhibitory signaling mode of PD-1, CEACAM1 represses the anti-tumor activity of T cells by dephosphorylating the downstream kinases of T cell receptor signaling (24)(25)(26). Monoclonal antibodies targeting CEACAM1 have been approved for the treatment of advanced and recurrent cancers in clinical trials (27). Downregulation of CXCL12 in OS facilitates the release of tumor cells from the bone and metastasis to other tissues (28). However, the protective coating formed by CXCL12 allowed malignant cells to escape immune attacks by T cells (29). After specific binding to CXCR4, the MAPK, PI3K, and phospholipase C pathways are activated, and the antitumor immune response is suppressed (30). LTF is widely considered a tumor suppressor. Consistent with the results of this study, patients with OS with low LTF levels showed lower survival rates, which may be attributed to the inhibitory effect of LTF on tumor cell proliferation (31). BNIP3, CDKN2A, PHLDA2, and UCHL1 are genes that we found highly expressed in OS samples. Upregulation of BNIP3 has been reporte d to enhance anoikis re sistance in hepatocarcinoma cells (32). A strong correlation between high BNIP3 levels and lower progression-free survival has also been observed in some platinum-resistant tumors (33). CDKN2A encodes a protein called p14ARF, which binds MDM2 in the nucleus and binds it to the nucleolus, thereby attenuating the ubiquitination degradation of p53 caused by MDM2. At  (36). The GA and AA genotypes of rs3217992 in CDKN2A may indicate higher malignancy, higher risk of lung metastasis, and poorer prognosis (37). PHLDA2 encodes a pleckstrin homology domain-containing protein that inhibits cell proliferation by suppressing AKT activation (38). Decreased PHLDA2 expression increases cell proliferation and reduces sensitivity to targeted agents in EGFR/ErbB2-driven cancer (39). How changes in PHLDA2 affect the development of OS requires further investigation. UCHL1 is a de-ubiquitinating enzyme that has been found to be over-expressed in some cancers and is considered a cancer promoter. UCHL1 downregulation decreases the proliferation, migration, and invasion of lung adenocarcinoma cells (40). It has also been proven to induce Correlations between risk scores and clinical features. metastasis of breast cancer cells by acting on TGF-b signaling (41). The high expression of UCHL1 in OS observed in this study suggests that it may play a role in tumor progression. Changes in the relevant pathways in the differential genes of the low-and high-risk groups were further obtained by GSVA analysis, which showed that pathways, such as JAK-STAT, tolllike receptor, and Hedgehog signaling were substantially enriched. STAT5 downregulation inhibits the proliferation, clonogenicity, and growth of OS cells (42). The toll-like receptor signaling pathway is remarkably differentially expressed in OS and is involved in the regulation of apoptosis, inflammation, and immunity (43). Similar to other studies, this study also confirmed the abnormality of the Hedgehog signaling pathway in OS, which may be related to tumor metastasis (44,45). Further investigation of the regulatory mechanisms of signaling pathways in different risk groups may provide additional information to better understand the heterogeneity within tumors.
The TME broadly consists of tumor, immune, and stromal elements and has been proven to determine the biological behavior of tumor cells. As expected, we found a pattern of low scores of immune, stromal, and ESTIMATE, and a high score of tumor purity in high-risk groups, which is commonly observed in malignant solid tumors; this was confirmed by the significantly lower survival rate we found in high-risk populations (46)(47)(48). To further characterize immune infiltration between the groups, we identified seven TME cells that were substantially different in the OS samples. The number of these cells was low in the high-risk group and negatively correlated with the risk score. BNIP3 overexpression accelerates the death of macrophages and T cells and promotes tumor proliferation and early metastasis (49,50). CXCL12 is an important chemokine in T and NK cells that helps macrophages polarize into tumor-associated macrophages (51, 52). In addition, CXCL12 mediates the progression of rectal cancer by promoting the retention of neutrophils in tumors and increasing their interactions with CD8+ T cells (53). These findings not only suggest that lowlevel infiltration of immune cells in OS high-risk samples may be associated with poor prognosis but also highlight the non-negligible role of BNIP3 and CXCL12 in the regulation of immune cell biological behavior.
Besides immune cell infiltration, the immune response gene set and HLA family genes were also found to be differentially expressed between the two risk groups in our study. Through the correlation analysis, we found that antigen processing and presentation and HLA-DMA were the two entries that most negatively correlated with BNIP3, whereas antimicrobials and HLA-DOA were most positively correlated with CXCL12. HLA families participate in tumor immunity (54). The downregulation of HLA genes may reduce antigen presentation and facilitate immune evasion (55). HLA-DMA variants have been reported to be associated with a higher risk of local recurrence in head and neck squamous cell carcinoma (56). HLA-DOA is also a key molecule in the antigen processing and presentation pathway and has been implicated in OS progression through the downregulation of HLA-DOA expression (57, 58). Another study has described the inhibitory role of HLA-DOA in B cell-mediated antigen presentation (59). Based on the present evidence, we hypothesized that the immune dysregulation associated with BNIP3 and CXCL12 may drive the poor prognosis of OS.
In the last few decades, immune checkpoint inhibitor-based immunotherapies have provided a huge boost to research on immune surveillance and have transformed the therapeutic landscape of cancer (60). However, immune checkpoint blocking therapy is less effective in treating OS, with only 5% of patients with OS achieving objective remission in a 2017 clinical trial of the PD-1 antibody (61,62). To clarify the differences in the response to immunotherapy between the high-and low-risk groups, we examined the expression levels of common immune checkpoints. We found that CD274, CTLA4, and HAVCR2 were substantially under-expressed in the high-risk groups and were negatively correlated with risk scores. Similar to previous studies, patients in the high-risk group showed insensitivity to PD-1 and CTLA4 FIGURE 7 GSVA enrichment analysis of the biological pathways between the two risk groups. Green indicates significant down-regulation, blue indicates significant up-regulation, and gray indicates no significance. GSVA, gene set variation analysis. Immune response gene sets in the two risk groups. inhibitors, which may partially explain the poor response of some patients with OS to immunotherapy in clinical practice (62,63). In addition, the 24 agents in our study showed sensitivity in high-risk groups and may help improve OS outcomes. In the future, combinations of immune checkpoint inhibitors with chemotherapy, targeted therapies, or novel therapies could potentially lead to new treatment strategies (2, 64, 65).

Conclusion
To the best of our knowledge, this is the first study to identify OS-related DEARGs and explore their predictive power for disease prognosis. The stratification framework based on BNIP3 and CXCL12 can effectively screen individuals at high risk for OS. Compared with those of the low-risk group, the high-risk group Prediction of response to immunotherapy. (A) The expression of three immune checkpoints was significantly different between the groups. (B) Correlation between the risk score and expression of immune checkpoints. had unique immune characteristics and sensitivity to drug therapy, which may provide a scientific reference for clinicians to develop efficient treatment strategies. Further experimental and clinical studies based on our results are promising to consistently improve the prognosis of patients with high-risk OS.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found within the article/supplementary material.

Author contributions
SZ and LR led the study design and prepared the manuscript. XZ, ZW, and QW performed the research. XZ wrote the manuscript. SZ revised the manuscript. All authors contributed to the article and approved the submitted version.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher. B C A FIGURE 12 Prediction of response to therapeutic drugs. (A) Sensitivity of different groups to PD-1 and CTLA4 inhibitors. (B) Agents with lower IC50 in the highrisk group than those in the low-risk group. (C) Agents with higher IC50 in the high-risk group than those in the low-risk group. IC50, 50% inhibitory concentration.