A novel molecular signature for predicting prognosis and immunotherapy response in osteosarcoma based on tumor-infiltrating cell marker genes

Background: Tumor infiltrating lymphocytes (TILs), the main component in the tumor microenvironment, play a critical role in the antitumor immune response. Few studies have developed a prognostic model based on TILs in osteosarcoma.

Methods: ScRNA-seq data was obtained from our previous research and bulk RNA transcriptome data was from TARGET database. WGCNA was used to obtain the immune-related gene modules. Subsequently, we applied LASSO regression analysis and SVM algorithm to construct a prognostic model based on TILs marker genes. What’s more, the prognostic model was verified by external datasets and experiment in vitro.

Results: Eleven cell clusters and 2044 TILs marker genes were identified. WGCNA results showed that 545 TILs marker genes were the most strongly related with immune. Subsequently, a risk model including 5 genes was developed. We found that the survival rate was higher in the low-risk group and the risk model could be used as an independent prognostic factor. Meanwhile, high-risk patients had a lower abundance of immune cell infiltration and many immune checkpoint genes were highly expressed in the low-risk group. The prognostic model was also demonstrated to be a good predictive capacity in external datasets. The result of RT-qPCR indicated that these 5 genes have differential expression which accorded with the predicting outcomes.

Conclusions: This study developed a new molecular signature based on TILs marker genes, which is very effective in predicting OS prognosis and immunotherapy response.

1 Introduction

Osteosarcoma (OS), the most common primary bone tumor in children and adolescents (1, 2), is characterized by high aggressiveness and poor prognosis (3). It is a relatively rare cancer with global incidence rate of approximately 4.8 per million per year (4, 5). Although treatments including surgery, chemotherapy and immunotherapy have made some progress in recent years, the overall survival rate is still unsatisfactory and fluctuates between 60% and 70% (6, 7). Therefore, it is particularly important to identify reliable prognostic markers and construct new molecular models to accurately predict the survival trend of patients with osteosarcoma.

The tumor microenvironment is a harbor for tumor cells and influences tumor development (8). Tumor infiltrating lymphocytes (TILs), one of the most common cell composition in tumor microenvironment, is composed of innate immune cells, adaptive immune cells, immunoreactive cells (e.g. cytotoxic T lymphocytes) and immunosuppressive cells (e.g. regulatory T cells) (9–11). Secrete cytokines derived from TILs can induce migration and aggregation of CXCR5-expressing B and T cells to suppress tumor progress. Chemotactic B cells and T cells form well-organized structures can prevent tumor metastasis to other sites (12, 13). In the context of antitumor immunotherapy, TILs are gradually gaining attention. In laryngeal squamous cell carcinoma, Sara et al. concluded that the number of TILs was positively correlated with PD-L1 expression and a good prognosis (14). In colorectal cancer, Yu-jie Liang et al. concluded that patients with high levels of Foxp3+ T cells had a better prognosis (15, 16). Thus, TILs could be used as biomarkers with good prognostic predictive value in a variety of tumors (12, 17). Although some prognosis gene models of OS have been established (18), few studies have reported and developed a prognostic risk model based on TILs in OS.

We firstly found TILs marker genes based on single cell RNA sequencing (scRNA-seq) data. Subsequently, in order to screen the prognostic genes and established risk model, LASSO regression analysis and machine learning SVM algorithm were utilized. Ultimately, we developed a prognostic risk model successfully, which may offer a novel reference for predicting the prognosis of OS patients.

2 Materials and methods

2.1 Preliminary experiment and raw data acquisition

ScRNA-seq sequencing data was obtained from our previous research, in which we collected tumor tissue from 6 OS patients from our hospital and conducted single-cell RNA sequencing (19). The basic information of the patients was showed in Table S1. Tumor tissues collecting from surgery were cut into pieces of approximately 1 mm³ in size and converted into cell suspensions for use. The sequencing work was performed in a double-end sequencing mode and i.e. 150 bases were measured at both Read 1 and Read 2 ends. We followed the 10X Genomics’ official process to perform upstream analysis using Cell Ranger software (version 4.0.0). We aligned the single-cell data with the human genome sequencing reference library GRCh38. The barcode.tsv file, gene.tsv file and matrix.mtx file were then derived by pairing read lengths, generating feature barcode matrices, performing clustering and other secondary analyses (19, 20). Furthermore, the data used for prognostic signature model development were from TARGET databank (https://ocg.cancer.gov/programs/target). Two external verification databases named GSE21257 and GSE16091 were downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo).

2.2 Identification of TILs marker genes

In this study, single-cell RNA sequencing data was analyzed by “Seurat” and “SingleR” R packages (21). In order to seek osteosarcoma cells, we rigorously filtered the raw matrix data of each cell based on three filtering criteria: nFeature_RNA > 300, nFeature_RNA < 4500 and percent. mt < 15. First, we set the mode to “LogNormalize” when we normalized the data through “NormalizeData” function of the “Seurat” R package. The normalized data was then transformed into a Seurat object, which was processed by “FindVariableFeatures” function to identify the genes of ideal cells. Next, “RunPCA” function of the “Seurat” R package analyzes the top 2000 highly variable genes of the target cells. The result of PCA were presented as PCA scatter plots. The top 30 principal components (PCs) were identified by JackStraw analysis. The cell clustering analysis was performed by using “FindNeighbors” and “FindClusters” functions of the “Seurat” R package. The clustering results were visualized as t-distributed random neighborhood embeddings (t-SNE) by “RunTSNE” function. We used “FindAllMarkers” function of the “Seurat” R package to identify differentially expressed genes (DEGs) for each cluster according to the criteria of adjusted p < 0.05 and |log2(FC)| >0.25. The “SingleR” R package was leveraged to pair the DEGs with marker genes from various cell types in the human primary cell atlas, thus enabling the annotation of cell clusters. Besides, we implement metabolomic analysis of TILs with “scMetabolism” R package (22).

2.3 Immune-related co-expression analysis

We first extracted the expression profile data of TILs marker genes from TARGET cohorts which includes 85 RNA expression matrix. On the basis of expression profile data of TILs marker genes, we calculated tumor-associated stromal scores, immune scores, ESTIMATE scores and tumor purity using “estimate” R package (23). We then used “WGCNA” R package to find the modules that were most significantly positively correlated with stromal scores, immune scores and ESTIMATE scores but most significantly negatively correlated with tumor purity. These genes included in the selected modules were defined as immune-associated TILs marker genes and used in the construction of a prognostic signature model.

2.4 Construction of a prognostic signature based on TILs marker genes

We first integrated the expression profile data of immune associated TILs marker genes with the clinical phenotype data of 85 TARGET cohorts, and then used univariate Cox regression analysis to screen for genes with prognostic value. In the LASSO regression analysis implemented by “glmnet” R package, we used “cv.glmnet” function to perform a 10-fold cross-validation of prognosis genes and then obtained genes with non-zero β coefficients. Simultaneously, according to the survival status, patients in the TARGET cohorts were assigned into survival or death group. In the machine learning SVM regression analysis, “svmRadial” function of the “e1071” R package was used to cross-verify the prognostic genes and aggregate the genes with the lowest error. The genes obtained from both analysis methods were intersected and common genes were extracted. The multivariate Cox regression analysis was the final step in creating the TILscore and assessing risk. The risk score of TILScore was calculated based on the equation “risk score = Σexpgenei*βi”. The “expgene” is the expression value of the model gene and “β” is the risk coefficient of the model gene. Median risk score splitted OS patients of the TARGET cohorts into two risk groups.

2.5 Validation of prognostic signature based on immune-related TILs marker genes

To evaluate the prognostic performance of TIL score, we used “ survivalROC” R package to obtain time-dependent ROC curves (24), and area under curve (AUC) values reflect the predictive ability of the model for patients’ overall survival at 1, 3, and 5 years. We invoked “survminer” R package in the Kaplan-Meier survival analysis to investigate survival differences between the high-risk and low-risk groups. In addition, we examined the possibility of TILScore as an independent prognostic factor by independent prognostic analysis. We combined TILScore with other clinical phenotypes to form a clinical nomogram. The clinical nomogram initially predicted the prognosis of different patients at 1,3,5 years.

In order to further verify the reliability of the model, we collected the expression matrix and clinical information from two GEO databases (GSE21257 and GSE16091). After removing the batch effect, the two databases were merged by Combat function of “sva” package. Subsequently, the merged database will be used as external verification database.

2.6 Tumor immune landscape assessment and immunotherapy response prediction

First, we extracted the expression profiles of 2044 TILs genes with cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT) algorithm and gained the infiltration abundance of 22 immune cell types. Subsequently, we analyzed the proportion of immune cell infiltrates in the two risk groups. In addition, we performed differential analysis of the expression matrix data of the immune checkpoint genes and observed the expression levels between the two risk groups.

2.7 Function and pathway enrichment analysis

In this study, we performed Genes Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses with “clusterProfiler” R package (25). We used “clusterProfiler” R package to visualize the functions and pathways of TILs marker genes. The GO annotations are based on the genome-wide annotation package released by the Bioconductor project (org.Hs.eg.db). The KEGG annotations were queried through the web API in the “clusterProfiler” R package for the latest online KEGG database. P value <0.05 was considered as a significant enrichment criteria.

2.8 Cell cultures

In this study, human osteoblasts cell (OB) were used as a control group, while human OS cells including 143B and HOS cells as experimental cell. All cells were purchased from Procell Life Science&Technology Co.,Ltd. (Wuhan, China). 143B cells and HOS cells were cultured in 1640 medium (Gibco, USA) and MEM medium, respectively. The two medium was supplemented with 1% penicillin/streptomycin (Solarbio, Beijing, China) and 10% fetal bovine serum (FBS; Gibco). OB cells cultured in DMEM/F-12 medium. OS cells were cultured in a humidified 5% CO2 incubator at 37 °C while OB cells was in a humidified 5% CO2 incubator at 35 °C.

2.9 RT-qPCR assays

In accordance with manufacturer’s instructions, total RNA was extracted using RNA fast 200 Kit (Fastagen Biotech, Shanghai, China). RNA was reverse-transcribed into complementary DNA (cDNA) by using a cDNA synthesis kit (Takara, Japan). RT-qPCR was performed using SYBR Green (FastStart Universal SYBR Green Master ROX, Germany) on a StepOnePlus™ Real-Time PCR System (ABI7500). The PCR procedure was as follow: 95 °C for 10 minutes, followed by 40 cycles at 95 °C for 10 seconds and 60 °C for 1 minute. The gene expression level in cell lines was expressed as relative expression and calculated using the ^2-ΔΔCt method. The primer sequences can be found in Table S4.

2.10 Statistical analysis

In this study, data was analyzed and generated using R software version 4.1.0 (http://www.R-project.org). P < 0.05 was the significance threshold and the “p.ucdilg” R function was used to adjust the P values for multiple analyses.

3 Results

3.1 Identification and annotation of TILs clusters

First, quality control criterias (nFeature_RNA > 300, nFeature_RNA < 4500, %. mt < 15) were used to filter single-cell RNA sequencing data of 6 patients with osteosarcoma. We obtained 24,611 genes from high quality cells and selected the top 2,000 high variance genes (Figures 1A, B). Next, we incorporated these 2000 high-variance genes into PCA analysis to reduce the dimensionality of the sequencing data. The cells were then subjected to co-expression clustering analysis to obtain 11 cell clusters (Figure 1C). Subsequently, each cluster were annotated by automatically coordinating the “SingleR” R package and manually based on marker genes. We defined the cells of cluster 1 and cluster 9 as TILs with the marker genes NKG7, CD3D, GZMK and GZMB (Figures 1D–F). After integrating the data of cluster 1 and cluster 9, a total of 2078 marker genes were obtained. In addition, we found that the single cell metabolic features of TILs are related to glycolysis or gluconeogenesis (Figure 1G).

FIGURE 1

Figure 1 Quality control of cells and the single-cell RNA sequencing to identify TILs markers. (A) Violin plot of 3 cell quality control standards. (B) The top 2,000 high variance genes. (C) t-SNE plots sorted by cell clusters. (D) t-SNE plots colored by the same cell type. (E) Different cell clusters identified by marker genes. (F) TILs identified by 6 marker genes. (G) single cell metabolic features of TILs are related to glycolysis or gluconeogenesis.

3.2 Functions and pathways of TILs marker genes

Enrichment analysis was used to understand the role of TILs in the anti-tumor immune response. First, we used the previous 2078 TILs marker genes in combination with genes expression matrix data from 85 TARGET cohorts to obtain expression profile data for 2044 marker genes (Figure 2A). The results of GO analysis showed that these marker genes are mainly involved in the biological processes of immune cell adhesion, cell migration and protein translation, including regulation of leukocyte cell-cell adhesion, focal adhesion, cell-substrate junction and cadherin binding. The KEGG analysis also confirmed the close association of these genes with immune cell adhesion and ribosomes (Figures 2B, C)

FIGURE 2

Figure 2 Biological functions and pathways of TILs marker genes. (A) Venn diagram of 2044 TILs marker genes from TARGET OS cohort and scRNA-seq data. (B) Bubble plot of the functions and pathways of TILs marker genes. (C) Network plot of the functions and pathways of TILs marker genes.

3.3 Establishment of a prognostic signature on the basis of five TILs marker genes

The optimal soft threshold was set to 7 and the results of WGCNA analysis showed that a total of 7 gene modules were identified. The 545 genes in the blue-green module were most significantly positively correlated with stromal scores, immune scores, and ESTIMATE scores, but most significantly negatively correlated with tumor purity (Figures 3A, B). Thus, the genes in the blue-green module are marker genes associated with immunity. Next, we performed a univariate Cox regression analysis on 545 marker genes and found that 140 genes were significantly associated with overall survival (Figure S1). Subsequently, we screened 11 genes and 34 genes from the 140 prognostic genes using LASSO Cox regression analysis and machine learning SVM regression analysis, respectively (Figures 3C-E and Table S2). The genes obtained from the two analysis methods were taken to intersect to obtain 6 common genes (Figure 3F and Table S2). After multivariate Cox regression analysis, 5 genes were remained and used to construct a prognostic signature model. We also performed risk assessment based on the model. In addition, the risk score of the prognostic model = (0.705 × EPHX2 expression value) + (0.478 × FDPS expression value) + (-0.35 × GBP1 expression value) + (-0.726 × MMD expression value) + (-0.815 × ZYX expression value) (Table S3).

FIGURE 3

Figure 3 Construction of the TILscore. (A) Module thresholds for WGCNA analysis. (B) Coefficients of different modules showing the correlation of TILs marker genes with 4 microenvironment scores. (C, D) Coefficient and parameter plots of LASSO analysis showing the 11 filtered candidate genes and pathways of TILs marker genes. (E) Machine learning SVM analysis curve showing 34 filtered candidate genes. (F) Venn diagram of common candidate genes from LASSO analysis and machine learning SVM analysis.

3.4 Survival analysis of TILScore

Their median value of the risk score was 0.874 and categorized the patients in TARGET cohorts into a low risk group (n = 43) and a high risk group (n = 42). Figure 4A showed the relationship among risk scores, clinical phenotype, and modeled genes. Risk model are all had prognostic value in different sub-groups except the group of age (Figure 4B and Figure S2). Kaplan-Meier analysis also showed a more advantageous overall survival rate in the low-risk group (Figure 4C). The AUC values of the ROC analysis were 0.852, 0.859, 0.872 at 1, 3, 5 years, which indicated the high predictive accuracy of the prognostic signature model (Figure 4D). Figures S3A, B showed the survival status of patients with different risk scores. Independent prognostic analysis using “age, sex, tumor metastasis, TILScore” revealed that TILScore (HR: 1.244,95% CI: 1.149-1.347, P<0.001) and tumor metastasis (HR: 4.691,95% CI: 2.160-10.187, P<0.001) were independent prognostic factors (Figures S4A, B). Finally, nomograms was constructed to predict the survival state at the 1st, 3rd, and 5th years (Figures S4C, D). To verify the reliability of the model, we combined two GEO datasets and ultimately obtained 82 examples. According to the quartile method, patients were divided into two groups by the risk score calculated by our formula. The result showed that patients in high risk group had low overall survival rate (p = 0.085) (Figure 4E). What’s more, the result of RT-qPCR released that, compared to OB cells, EPHX2 and FDPS were significantly overexpressed in 143B and HOS cells, while the expression of ZYX and GBP1 was significantly lower in 143B and HOS cells. The expression trend of MMD was not definite (Figure 4F). In conclusion, we successfully established a 5-gene TILScore based on TILs marker genes.

FIGURE 4

Figure 4 Survival analysis and predictive performance evaluation of TILscore. (A) Heat map of the relationship between different clinical phenotypes and two risk groups. (B) 17 years of age as the threshold for survival differences between the two risk groups. (C) Kaplan-Meier curves for survival analysis compared overall survival of OS patients in the high-risk and low-risk groups. (D) TILscore ROC curves predict the risk of death at 1, 3, and 5 years. (E) Kaplan-Meier curves for survival analysis based on external verification database. (F) The result of RT-qPCR. *p < 0.05; **p < 0.01; ***p < 0.001; ns is for no statistical significance.

3.5 Immune landscapes associated with TILscore

Because TILs are critical in the antitumor immune response, we explored the relationship between TILScore and immune cell infiltration. Differential analysis of four tumor-associated microenvironment score revealed that high-risk patients had lower immune scores, stromal scores and ESTIMATE scores (Figures S5A, B). Thus, the high-risk score was negatively correlated with the level of immune cell infiltration. The CIBERSORT algorithm calculated associations between five modeled genes and 22 immune cells (Figure 5A). The differences in the infiltration abundance of memory B cells, naive B cells, Treg cells, and γδ T cells were statistically significant in two risk groups (Figures 5B, C). Details of the infiltration of the four immune cells mentioned above were shown in Figures 5D-G. Genes of the prognostic signature may further affect the development of osteosarcoma by influencing the antitumor response through B and T cells.

FIGURE 5

Figure 5 TILscore-associated immune landscape. (A) Heat map of the correlation between the 5 genes of TILscore and the abundance of 22 immune cell infiltrates. (B) Heat map showing the infiltration abundance of 22 immune cell species in different risk groups. (C) Box plots showing differences in immune cell infiltration by risk group. Differences among naive B cells (D), memory B cells (E), Treg cells (F) and gamma delta T cells (G) in the high-risk and low-risk groups. *p < 0.05; **p < 0.01; ***p < 0.001.

3.6 The TILScore predicts immunotherapy benefits in OS patients

TILs are involved in antitumor immune responses, so we analyzed the relationship between TILScore and immune checkpoints. The results revealed that some critical immune checkpoints, such as CTLA4, LAIR1, HAVCR2, CD48, CD44, CD27, LGALS9 and LAG3, possess higher expression in low-risk group (Figures 6A-H).

FIGURE 6

Figure 6 Differential analysis of immune checkpoints to assess the predictive power of TILscore immunotherapy. Differences in expression levels of immune checkpoints CTLA4 (A), LAIR1 (B), HAVCR2 (C), CD48 (D), CD44 (E), CD27 (F), LGALS9 (G), and LAG3 (H) in the low- and high-risk groups. *p < 0.05; **p < 0.01; ***p < 0.001.

4 Discussion

Previous studies have demonstrated that tumor-infiltrating lymphocytes are strongly associated with the prognosis of different tumor types (26–28), whereas few prognosis models have been conducted in osteosarcoma. In recent years, by using scRNA-seq analysis, Peng Song et al. developed a prognostic signature based on NK cell marker genes in patients with lung adenocarcinoma (29), suggesting that scRNA-seq analysis is an excellent tools to the construction of prognosis model. Therefore, we searched for marker genes of TILs by scRNA-seq analysis in present study. Our results showed that TILScore proved to be a powerful predictive model for OS patients’ prognosis. GO and KEGG analysis turned out that most of these TILs genes were enriched in biological processes such as immune cell adhesion, cell migration and protein synthesis pathways. Therefore, these genes may influence the proliferation and developmental processes of osteosarcoma through immune cell migration and cellular protein synthesis (30, 31).

Our risk model consisted of five TILs marker genes (EPHX2, FDPS, GBP1, MMD and ZYX). FDPS, an osteoclast farnesyl pyrophosphate synthase, is closely associated with osteosarcoma formation. Its main role is to promote osteoclast bone resorption activity and to inhibit osteoclast apoptosis (32). Upregulation of FDPS expression is also significantly associated with patients’ poor prognosis (33). MMD, one of the proteins associated with monocyte-to-macrophage differentiation, can enhance ERK1/2 and Akt phosphorylation in macrophages after LPS stimulation. It is involved in the antitumor immune response (34). Consistent with our findings, this gene is closely associated with phagosomes of tumor-infiltrating lymphocytes. Higher expression levels of MMD were found in the low risk group of this study and MMD was associated with good prognosis. ZYX is a LIM structural domain protein, which is involved in cytoskeletal organization and tumorigenesis (35, 36). ZYX is involved in apoptosis and osteoblast differentiation processes (37, 38). Therefore, ZYX may be involved in the development of osteosarcoma. High expression levels of ZYX were associated with good prognosis. EPHX2, the coding gene of soluble epoxide hydrolase (sEH) protein, is related to activity of epoxide hydrolase and phosphatase (39). EPHX2 was demonstrated to be strongly associated with cancer prognosis and macrophage phagocytosis (40). GBP1, encoding Guanylate Binding Protein 1, is related to tumor progression and chemotherapy drug resistance (41). What’s more, GBP1 was also regarded as microbe-specific gatekeeper of macrophage apoptosis and pyroptosis (42). Combining literature reports and our results of bioinformatics and RT-qPCR, we hold that EPHX2 and FDPS are the oncogenes, while GBP1, MMD and ZYX are the anti-oncogene in osteosarcoma. Risk model based on these five genes is very reliable to predict the prognosis of patients.

The differences in the immune landscape between risk groups allowed us to see the potential value of TILScore in predicting immunotherapy response. In the present study, naive B cells, Treg cells and γδ T cells were more distributed in the high-risk group, but memory B cells were more distributed in the low-risk group. Naïve B cells are undifferentiated B cells without antigen stimulation and infiltrate a high percentage in tumor tissue (43). Memory B cells are generated in the germinal center response during T cell-dependent immune response. Unlike naïve B cells, memory B cells are involved in a faster and stronger immune response (44). Therefore, the antitumor immune responses of memory B cells are faster, thus the prognosis of the low-risk group is better. Previous studies confirmed that Treg cells are involved in the tumor development process by suppressing antitumor immunity and its high expression level represents a poor prognosis (45). γδ T cells, a congenital cytotoxic T cells, were highly infiltrated in the low-risk group, which may be related to the fact that osteosarcoma is an immunocompromised tumor (46). Moreover, LAIR1, HAVCR2, CD27, CTLA4, CD48, CD44, LAG3 and LGALS9 were highly expressed in the low-risk group, indicating that the low-risk group was more likely to benefit from more types of immunotherapy. In summary, TILScore, a reliable biomarker for predicting response to immunotherapy, predicting that low-risk patients were more likely to benefit from immunotherapy.

There are some limitations of our study. First, a full functional experiment to elucidate the specific mechanisms of TILs marker genes in osteosarcoma is still very important. Second, the data involved in our study were derived from single-cell RNA sequencing data and the TARGET database, but the sample size remains insufficient. Therefore, the predictive power of this prognostic feature has some limitations.

5 Conclusion

Overall, this study developed a 5-gene prognostic signature based on TILs marker genes, which performed well in predicting prognosis and immunotherapy response in patients with osteosarcoma. TILScore can be considered to be an independent prognostic factor to predict patient prognosis and to guide patients to benefit from different immunotherapy methods.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics statement

The studies involving human participants were reviewed and approved by Ethics Committee and Institutional Reviewer Board of the First Affiliated Hospital of the Guangxi Medical University. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

Written informed consent was obtained from the individual(s), and minor(s)’ legal guardian/next of kin, for the publication of any potentially identifiable images or data included in this article.

Author contributions

All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the National Natural Science Foundation of China (82160536) and Youth Science Foundation of Guangxi Medical University (GXMUYSF202128).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1150588/full#supplementary-material

Supplementary Figure 1 | Forest plot showing 140 prognostic genes screened by univariate Cox regression analysis.

Supplementary Figure 2 | Survival differences between the two risk groups in subgroups of other clinical indicators. (A, B): gender group; (C, D): metastatic group.

Supplementary Figure 3 | The distribution of risk score (A) and survival status (B).

Supplementary Figure 4 | (A, B) TILscore and tumor metastasis are independent prognostic factors. (C, D) Nomogram for predicting survival status of OS patients at the 1st, 3rd and 5th year.

Supplementary Figure 5 | Visualization of differences in microenvironment scores across risk groups. (A) Heat map showing the difference between the 4 microenvironment scores in the high and low risk groups. (B) Box plots of differences in basal scores, immune scores, and ESTIMATE scores across risk groups.

References

1. Ferguson JL, Turner SP. Bone cancer: Diagnosis and treatment principles. Am Fam Physician (2018) 98:205–13.

PubMed Abstract | Google Scholar

2. Brown HK, Tellez-Gabriel M, Heymann D. Cancer stem cells in osteosarcoma. Cancer Lett (2017) 386:189–95. doi: 10.1016/j.canlet.2016.11.019

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Lei T, Qian H, Lei P, Hu Y. Ferroptosis-related gene signature associates with immunity and predicts prognosis accurately in patients with osteosarcoma. Cancer Sci (2021) 112:4785–98. doi: 10.1111/cas.15131

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Gamerdinger D, Bodem F, Brussatis F, Otte P. [A computer-assisted electromechanical measuring system for the study of flexion and extension forces in the upper ankle joint]. Z Orthop Ihre Grenzgeb (1987) 125:641–3. doi: 10.1055/s-2008-1039702

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Simpson E, Brown HL. Understanding osteosarcomas. JAAPA (2018) 31:15–9. doi: 10.1097/01.JAA.0000541477.24116.8d

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Wu CC, Livingston JA. Genomics and the immune landscape of osteosarcoma. Adv Exp Med Biol (2020) 1258:21–36. doi: 10.1007/978-3-030-43085-6_2

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Jafari F, Javdansirat S, Sanaie S, Naseri A, Shamekh A, Rostamzadeh D, et al. Osteosarcoma: A comprehensive review of management and treatment strategies. Ann Diagn Pathol (2020) 49:151654. doi: 10.1016/j.anndiagpath.2020.151654

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Arneth B. Tumor microenvironment. Medicina (Kaunas) (2019) 56. doi: 10.3390/medicina56010015

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hanahan D, Coussens LM. Accessories to the crime: functions of cells recruited to the tumor microenvironment. Cancer Cell (2012) 21:309–22. doi: 10.1016/j.ccr.2012.02.022

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell (2011) 144:646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol (2018) 1711:243–59. doi: 10.1007/978-1-4939-7493-1_12

PubMed Abstract | CrossRef Full Text | Google Scholar

12. De Meulenaere A, Vermassen T, Aspeslagh S, Vandecasteele K, Rottey S, Ferdinande L. TILs in head and neck cancer: Ready for clinical implementation and why (Not)? Head Neck Pathol (2017) 11:354–63. doi: 10.1007/s12105-016-0776-8

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Paijens ST, Vledder A, De Bruyn M, Nijman HW. Tumor-infiltrating lymphocytes in the immunotherapy era. Cell Mol Immunol (2021) 18:842–59. doi: 10.1038/s41423-020-00565-9

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Pai SI, Westra WH. Molecular pathology of head and neck cancer: implications for diagnosis, prognosis, and treatment. Annu Rev Pathol (2009) 4:49–70. doi: 10.1146/annurev.pathol.4.110807.092158

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Bron L, Jandus C, Andrejevic-Blant S, Speiser DE, Monnier P, Romero P, et al. Prognostic value of arginase-II expression and regulatory T-cell infiltration in head and neck squamous cell carcinoma. Int J Cancer (2013) 132:E85–93. doi: 10.1002/ijc.27728

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Liang YJ, Liu HC, Su YX, Zhang TH, Chu M, Liang LZ, et al. Foxp3 expressed by tongue squamous cell carcinoma cells correlates with clinicopathologic features and overall survival in tongue squamous cell carcinoma patients. Oral Oncol (2011) 47:566–70. doi: 10.1016/j.oraloncology.2011.04.017

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Zhou Y, Tian Q, Wang BY, Yang J, Zhao SD, Yang J. The prognostic significance of TILs as a biomarker in triple-negative breast cancer: what is the role of TILs in TME of TNBC? Eur Rev Med Pharmacol Sci (2021) 25:2885–97. doi: 10.26355/eurrev_202104_25542

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Yang W, Wu H, Tong L, Wang Y, Guo Q, Xu L, et al. A cuproptosis-related genes signature associated with prognosis and immune cell infiltration in osteosarcoma. Front Oncol (2022) 12:1015094. doi: 10.3389/fonc.2022.1015094

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Liu Y, Feng W, Dai Y, Bao M, Yuan Z, He M, et al. Single-cell transcriptomics reveals the complexity of the tumor microenvironment of treatment-naive osteosarcoma. Front Oncol (2021) 11:709210. doi: 10.3389/fonc.2021.709210

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Rothzerg E, Feng W, Song D, Li H, Wei Q, Fox A, et al. Single-cell transcriptome analysis reveals paraspeckles expression in osteosarcoma tissues. Cancer Inform (2022) 21:11769351221140101. doi: 10.1177/11769351221140101

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol (2019) 20:163–72. doi: 10.1038/s41590-018-0276-y

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Wu Y, Yang S, Ma J, Chen Z, Song G, Rao D, et al. Spatiotemporal immune landscape of colorectal cancer liver metastasis at single-cell level. Cancer Discov (2022) 12:134–53. doi: 10.1158/2159-8290.CD-21-0316

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Heagerty PJ, Zheng Y. Survival model predictive accuracy and ROC curves. Biometrics (2005) 61:92–105. doi: 10.1111/j.0006-341X.2005.030814.x

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an r package for comparing biological themes among gene clusters. OMICS (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Santoiemma PP, Powell DJ Jr. Tumor infiltrating lymphocytes in ovarian cancer. Cancer Biol Ther (2015) 16:807–20. doi: 10.1080/15384047.2015.1040960

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Stanton SE, Disis ML. Clinical significance of tumor-infiltrating lymphocytes in breast cancer. J Immunother Cancer (2016) 4:59. doi: 10.1186/s40425-016-0165-6

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Schatton T, Scolyer RA, Thompson JF, Mihm MC Jr. Tumor-infiltrating lymphocytes and their significance in melanoma prognosis. Methods Mol Biol (2014) 1102:287–324. doi: 10.1007/978-1-62703-727-3_16

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Song P, Li W, Guo L, Ying J, Gao S, He J. Identification and validation of a novel signature based on NK cell marker genes to predict prognosis and immunotherapy response in lung adenocarcinoma by integrated analysis of single-cell and bulk RNA-sequencing. Front Immunol (2022) 13:850745. doi: 10.3389/fimmu.2022.850745

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Ozga AJ, Chow MT, Luster AD. Chemokines and the immune response to cancer. Immunity (2021) 54:859–74. doi: 10.1016/j.immuni.2021.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Challa S, Khulpateea BR, Nandu T, Camacho CV, Ryu KW, Chen H, et al. Ribosome ADP-ribosylation inhibits translation and maintains proteostasis in cancers. Cell (2021) 184:4531–46.e4526. doi: 10.1016/j.cell.2021.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Body JJ. Rationale for the use of bisphosphonates in osteoblastic and osteolytic bone lesions. Breast (2003) 12 Suppl 2:S37–44. doi: 10.1016/s0960-9776(03)80162-5

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Jiang H, Du H, Liu Y, Tian X, Xia J, Yang S. Identification of novel prognostic biomarkers for osteosarcoma: a bioinformatics analysis of differentially expressed genes in the mesenchymal stem cells from single-cell sequencing data set. Transl Cancer Res (2022) 11:3841–52. doi: 10.21037/tcr-22-2370

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Liu Q, Zheng J, Yin DD, Xiang J, He F, Wang YC, et al. Monocyte to macrophage differentiation-associated (MMD) positively regulates ERK and akt activation and TNF-alpha and NO production in macrophages. Mol Biol Rep (2012) 39:5643–50. doi: 10.1007/s11033-011-1370-5

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Partynska A, Gomulkiewicz A, Dziegiel P, Podhorska-Okolow M. The role of zyxin in carcinogenesis. Anticancer Res (2020) 40:5981–8. doi: 10.21873/anticanres.14618

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Bach I. The LIM domain: regulation by association. Mech Dev (2000) 91:5–17. doi: 10.1016/s0925-4773(99)00314-7

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Kotb A, Hyndman ME, Patel TR. The role of zyxin in regulation of malignancies. Heliyon (2018) 4:e00695. doi: 10.1016/j.heliyon.2018.e00695

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Hervy M, Hoffman L, Beckerle MC. From the membrane to the nucleus and back again: bifunctional focal adhesion proteins. Curr Opin Cell Biol (2006) 18:524–32. doi: 10.1016/j.ceb.2006.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Li H, Bradbury JA, Edin ML, Graves JP, Gruzdev A, Cheng J, et al. sEH promotes macrophage phagocytosis and lung clearance of streptococcus pneumoniae. J Clin Invest (2021) 131. doi: 10.1172/JCI129679

CrossRef Full Text | Google Scholar

40. Zhou Y, Li X, Guan A, Zhou H, Zhu Y, Wang R, et al. EPHX2 inhibits colon cancer progression by promoting fatty acid degradation. Front Oncol (2022) 12:870721. doi: 10.3389/fonc.2022.870721

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Cheng L, Gou L, Wei T, Zhang J. GBP1 promotes erlotinib resistance via PGK1−activated EMT signaling in non−small cell lung cancer. Int J Oncol (2020) 57:858–70. doi: 10.3892/ijo.2020.5086

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Fisch D, Bando H, Clough B, Hornung V, Yamamoto M, Shenoy AR, et al. Human GBP1 is a microbe-specific gatekeeper of macrophage apoptosis and pyroptosis. EMBO J (2019) 38:e100926. doi: 10.15252/embj.2018100926

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Zhong R, Chen D, Cao S, Li J, Han B, Zhong H. Immune cell infiltration features and related marker genes in lung cancer based on single-cell RNA-seq. Clin Transl Oncol (2021) 23:405–17. doi: 10.1007/s12094-020-02435-2

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Seifert M, Kuppers R. Human memory b cells. Leukemia (2016) 30:2283–92. doi: 10.1038/leu.2016.226

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Ohue Y, Nishikawa H. Regulatory T (Treg) cells in cancer: Can treg cells be a new therapeutic target? Cancer Sci (2019) 110:2080–9. doi: 10.1111/cas.14069

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Liu Q, Xu R, Xu X, Huang Y, Ma Z. Characteristics and significance of T lymphocyte subsets in peripheral blood of osteosarcoma mice. Transl Cancer Res (2022) 11:1503–9. doi: 10.21037/tcr-22-264

PubMed Abstract | CrossRef Full Text | Google Scholar