All study participants were sourced from the Healthy Zhejiang One Million People (HOPE) Cohort. A total of 67 osteosarcoma patients were selected from the bone cancer patient subcohort of the HOPE cohort, initiated in June 2020 and ongoing. These patients were newly diagnosed, pathologically confirmed cases of osteosarcoma who had not undergone any surgical or chemotherapy treatment. Fifty healthy controls were chosen from the general population subcohort of the HOPE cohort. To reduce confounding effects, osteosarcoma patients and healthy controls were matched by age (± 5 years) and sex. This research received approval from the Institutional Review Board of the Second Affiliated Hospital of Zhejiang University. All participants provided signed informed consent. Clinical data were obtained from medical records, while follow-up information was gathered through regular telephone interviews. Patients with metastases were defined as those diagnosed with metastases during follow-up visits at our hospital or confirmed through telephone follow-up. The definitions of recurrence and death were established using similar criteria. Epidemiological data were collected via face-to-face interviews using a standardized questionnaire.

Each participant had 20 mL of blood collected into three vacutainer tubes (Thermo Fisher Scientific, USA): two lavender-colored tubes containing EDTA-Na and one red-colored tube without additives. Plasma was separated from the tubes containing EDTA-Na and stored at -80°C for future use. Tumor tissue samples were collected during surgery, immediately snap-frozen, and stored at -80°C for future use.

Additionally, mRNA expression data for 39 osteosarcoma tumor tissue samples were retrieved from the Gene Expression Omnibus (GEO) dataset GSE21257. scRNA-seq data for osteosarcoma tumor tissues were obtained from the GEO dataset GSE152048.

We used the ProcartaPlex human immune checkpoint panel (Thermo Fisher, USA) in a 96-well plate format to analyze plasma samples in order to quantify immune checkpoint proteins, including CD48, B7-H2, TIMD-4, B7-H6, CD134, B7-H5, CD47, S100A8/A9, and B7-H3. The assays were conducted using the FLEXMAP 3D system (Luminex, USA) and the xPONENT^® 4.3 software. To reduce batch effects in the Luminex assays, we included replicate standard samples across all batches for calibration and processed all samples within a single run. A double-blind methodology was maintained throughout the detection process. The quantification procedures have been detailed in our previous study (13). The lower limit of quantification for the proteins is provided in Supplementary Table 1 .

Total RNA was isolated using the AllPrep DNA/RNA/miRNA Universal Kit (Thermo Fisher, USA), and reverse transcription was performed with HiScript III All-in-one RT SuperMix Perfect for qPCR (Vazyme, CN). cDNA amplification was then carried out using the ChamQ Universal SYBR qPCR Master Mix (Vazyme, CN) system. For each gene, two pairs of highly specific primers were initially designed, and the primer with the highest amplification efficiency was chosen for further analysis. GAPDH served as an internal control to determine the relative mRNA levels of CD48, ICOSLG, TIMD4, NCR3LG1, TNFRSF4, VSIR, CD47, S100A8, and S100A9. All assays were performed in triplicate.

Frequency and percentage [n (%)] were used to present categorical variables, while medians [interquartile range (IQR)] or mean ± standard deviation was used for continuous variables. Comparisons of categorical variables used Pearson’s χ² test, while continuous variables were compared using the t-test or Wilcoxon rank-sum test. False discovery rate (FDR) correction was applied to the p-values to reduce the risk of Type I errors. Pathway enrichment analysis of differential immune checkpoint proteins was conducted using Metascape (https://metascape.org) (14).

After adjusting for age, sex, alcohol consumption, smoking, hypertension, BMI, and diabetes, multivariable unconditional logistic regression was applied to estimate the relationship between each plasma immune regulatory protein and the risks of osteosarcoma incidence and metastasis. Protein levels were considered continuous variables and log-transformed to minimize skewness. Multivariable Cox proportional hazards analysis, adjusted for the same covariates, was utilized to identify independent prognostic factors for osteosarcoma.

Unsupervised consensus clustering of osteosarcoma patient samples was conducted using the R package ConsensusClusterPlus (version 1.60.0) (15). The K-means algorithm was iterated 1000 times to ensure the stability of the clustering. The optimal number of clusters (k) was determined based on the cumulative distribution function and the proportion of ambiguous clustering. Osteosarcoma patients were classified into two categories to define plasma immune regulatory protein subtypes. The distribution differences among the clusters were visualized using principal component analysis (PCA).

Survival differences were compared using the log-rank test and Kaplan-Meier curves. A multivariable Cox proportional hazards regression model incorporating plasma immune regulatory protein subtypes and clinical variables was constructed to predict metastasis at six months (180 days), one year (365 days), and two years (730 days) post-diagnosis. Tumor staging was excluded from the Cox model, as all patients had similar stages at enrollment. Model discrimination was evaluated using Harrell’s concordance index (C-index) and time-dependent receiver operating characteristic (ROC) curves. Calibration was assessed using calibration curves generated from 1,000 bootstrap resamples. Overall model performance was measured using the integrated Brier score (IBS) (16).

For scRNA-seq data processing, Seurat (version 5.1.0) was used. Cells were filtered by restricting the number of detected genes per cell, total molecule count, and mitochondrial gene expression percentage, removing empty droplets, doublets, multiplets, and low-quality cells. Batch correction was performed using the Harmony package (version 1.2.3). Data normalization was conducted using the NormalizeData function in Seurat, and the identification of the top 2,000 highly variable genes was performed with the FindVariableFeatures function, followed by centering and scaling the data. PCA was performed to determine the number of principal components needed for further analysis. Unsupervised cell clustering was carried out through the FindClusters function in Seurat (default Louvain algorithm), and two-dimensional visualization was generated using Uniform Manifold Approximation and Projection (UMAP) and t-Distributed Stochastic Neighbor Embedding (tSNE). Cell type annotation was conducted by identifying highly differentially expressed genes for each cell cluster using the FindAllMarkers function in Seurat (default non-parametric Wilcoxon rank-sum test and Bonferroni correction). Annotation was based on the top ten differential genes for each cluster and known marker genes from the literature.

All statistical analyses, except for pathway enrichment analysis, were performed in R (version 4.3.2), with all tests being two-sided and a significance threshold of 0.05.

This study included 117 participants, consisting of 67 osteosarcoma patients and 50 healthy controls. Clinical and demographic characteristics are presented in Table 1 . The median age difference between the groups was 3.5 years. Among all participants, approximately 80% were male, with almost no history of smoking or alcohol consumption, and most had normal BMI, with minimal cases of diabetes or hypertension. Within the patient group, 14 cases (20.9%) experienced metastasis, 4 cases (6.0%) resulted in mortality, and 3 cases (4.5%) showed recurrence. The median follow-up time was 394 days (range: 11–1018).

Luminex multiplex assays were performed to measure 10 soluble immune checkpoint proteins in all participants. Since S100A8 and S100A9 frequently exist as the S100A8/A9 complex in plasma, we measured this complex (17). The B7-H3 assay was excluded due to an invalid standard curve. Protein levels are detailed in Table 2 .

With the exception of B7-H2, all measured immune checkpoint proteins—CD48, TIMD-4, B7-H6, CD134, B7-H5, CD47, and S100A8/A9—were significantly elevated in osteosarcoma patients compared to healthy controls ( Figure 1A , p < 0.05). Multivariate unconditional logistic regression analysis showed that higher levels of these protein markers were significantly linked to an increased risk of osteosarcoma ( Table 3 ).

Levels of eight immune checkpoint proteins (CD48, B7-H2, TIMD-4, B7-H6, CD134, B7-H5, CD47, and S100A8/A9) were significantly elevated in metastatic patients compared to those without metastasis ( Figure 1B , p < 0.05). Further analysis indicated that elevated levels of these markers, with the exception of B7-H2, were significantly associated with a higher risk of metastasis in osteosarcoma patients ( Table 3 ).

Pathway enrichment analysis of the differentially expressed immune checkpoint proteins identified four significantly enriched pathways: GO:0051251 (positive regulation of lymphocyte activation), GO:0046649 (lymphocyte activation), GO:0050729 (positive regulation of inflammatory response), and R-HSA-1280218 (Adaptive Immune System). All pathways had -log10(P) > 1.3 (p < 0.05) ( Supplementary Figure 1 ), suggesting that these proteins are closely associated with lymphocyte activation and inflammatory processes.

Unsupervised consensus clustering of osteosarcoma patients based on the levels of eight identified soluble immune checkpoint proteins (CD48, B7-H2, TIMD-4, B7-H6, CD134, B7-H5, CD47, S100A8/A9) classified patients into two clusters ( Figure 2A , Supplementary Figure 2 , k = 2): 31 in cluster 1 and 36 in cluster 2. Principal component analysis demonstrated a clear separation between these clusters ( Figure 2B ). Immune checkpoint protein levels were lower in cluster 1 compared to cluster 2, leading to their classification as Osteosarcoma Immunity Type I (cluster 1) and Osteosarcoma Immunity Type II (cluster 2). Clinical features did not significantly differ between these immune subtypes ( Figure 2C , Supplementary Table 2 ).

Multivariate Cox regression models indicated that soluble immune checkpoint protein subtypes were strong predictors of metastasis (HR = 10.87, 95% CI: 1.50–78.48) and progression-free survival (PFS) (HR = 8.35, 95% CI: 1.66–42.13) ( Figures 2D, E ). Kaplan-Meier analysis showed that Type II had significantly shorter PFS compared to Type I (log-rank p = 0.032) ( Figure 2F ).

Using soluble immune checkpoint protein subtypes and clinical variables (age, sex, smoking, alcohol consumption, BMI, hypertension, and diabetes), we constructed a multivariable Cox proportional hazards regression model to predict metastasis in osteosarcoma patients ( Figure 3A ). Tumor staging was excluded as all patients were at the same stage at enrollment. The C-index of the model at six months (180 days), one year (365 days), and two years (730 days) was 0.902, 0.876, and 0.648, respectively. The model demonstrated strong discriminative ability throughout the follow-up period, with superior performance at six months and one year. The area under the ROC curve (AUC) at six months, one year, and two years was 0.833, 0.889, and 0.739, respectively, while the IBS was 0.039, 0.071, and 0.160. Calibration curves for predicting metastasis at six months, one year, and two years showed that the predicted outcomes closely aligned with the observed results ( Figures 3B-D ). Decision curve analysis ( Supplementary Figure 3 ) showed that, within the threshold probability range of 10% to 70%, the 1-year and 2-year models yielded greater net benefit than the ‘treat-all’ and ‘treat-none’ strategies, suggesting their clinical utility in identifying patients at high risk of metastasis. The 1-year model demonstrated the widest range of clinical benefit, indicating broader applicability. These findings suggest that our model, incorporating soluble immune checkpoint protein subtypes, is a valuable tool for predicting osteosarcoma metastasis.

To further investigate the relationship between soluble immune checkpoint proteins and osteosarcoma and to preliminarily explore their mechanisms of action, we collected tumor tissues from 24 osteosarcoma patients (9 with metastasis and 15 without metastasis) out of the 67 who underwent plasma testing. We analyzed the expression levels of genes encoding soluble immune checkpoint proteins, including CD48, B7-H2, TIMD-4, B7-H6, CD134, B7-H5, CD47, S100A8, and S100A9 (genes: CD48, ICOSLG, TIMD4, NCR3LG1, TNFRSF4, VSIR, CD47, S100A8, S100A9) in these tumor tissues. Expression levels of all nine genes were observed to be lower in tumor tissues from metastatic patients compared to non-metastatic patients ( Figure 4A ). Although only the difference in S100A8 reached statistical significance, a clear trend of reduced expression was evident for the other genes. This lack of statistical significance may result from the limited sample size, given the rarity of osteosarcoma.

To validate these findings, we analyzed the expression levels in tumor tissues from 39 osteosarcoma patients (20 with metastasis and 19 without metastasis) available in the GEO database. Since the database did not include data on B7-H6 and B7-H5 expression, these genes were excluded from the analysis. We found that the expression levels of CD48, B7-H2, S100A8, and S100A9 were significantly lower in metastatic patients, with differences reaching statistical significance ( Figure 4B ). The expression levels of other soluble immune checkpoint protein genes also showed a noticeable decrease in metastatic patients. Interestingly, this contrasts with the findings observed in peripheral blood.

Based on the expression levels of CD48, B7-H2, TIMD-4, CD134, CD47, S100A8, and S100A9 (genes: CD48, ICOSLG, TIMD4, TNFRSF4, CD47, S100A8, and S100A9), we performed unsupervised consensus clustering of osteosarcoma patients in the GSE21257 dataset and identified two distinct immune subtypes: osteosarcoma immune subtype I (cluster 1, n = 18) and osteosarcoma immune subtype II (cluster 2, n = 21) ( Figures 5A–C ).

Interestingly, patients classified as osteosarcoma immune subtype I exhibited a higher risk of metastasis and shorter metastasis-free survival (log-rank p = 0.002) ( Figures 5D, E ), which is in contrast to the findings observed in peripheral blood.

To further explore the potential mechanisms underlying our findings, we analyzed scRNA-seq data from tumor tissues of 11 osteosarcoma patients in the GEO database. Batch effect correction was performed on the data ( Supplementary Figure 4 ). Unbiased clustering of the cells identified 15 major cell clusters ( Supplementary Figure 5A , Supplementary Table 3 ). All immune cell populations were extracted and subsequently re-clustered, resulting in 12 distinct immune cell types ( Supplementary Figure 5B , Supplementary Table 4 ). Expression differences of representative marker genes within these cell populations were quantitatively assessed and presented, with biological annotations applied ( Supplementary Figure 5C ). We calculated the proportions of various immune cell types in different prognostic groups ( Supplementary Figure 5D ).

We quantified the expression of these factor genes in immune cells ( Figure 6A ) and analyzed their differential expression across various immune cell types among different prognostic groups ( Figures 6B, C , Supplementary Table 5 , 6 ). Comparing metastatic patients with primary patients, we observed a significant elevation in the expression of the S100A8 gene in M1 macrophages in metastatic patients (p=0.003). In M2 macrophages, the expression of B7-H5, CD47, S100A8, and S100A9 was significantly upregulated in metastatic patients (p=0.037, p=0.007, p=0.003, p<0.001, respectively). For mDCs, the expression of CD48, TIMD-4, and CD47 was markedly increased in metastatic patients (p=0.022, p<0.001, p=0.012, respectively), while the expression of S100A8 and S100A9 was significantly reduced (p=0.005, p=0.012, respectively). In memory T cells, metastatic patients showed a significant downregulation of S100A9 expression (p<0.001). Additionally, in monocytes, the expression of CD48, TIMD-4, S100A8, and S100A9 was upregulated in metastatic patients (p=0.041, p=0.009, p<0.001, p=0.008, respectively). In neutrophils, the expression of CD48 and S100A8 was also elevated in metastatic patients (p<0.001, p=0.019, respectively). However, CD48 expression was significantly reduced in Tregs (p<0.001), which may impact their immunosuppressive function.

In conclusion, we reveal the relationship between the expression levels of immune checkpoint factor genes across various immune cells and osteosarcoma metastasis.