Machine learning-based selection of immune cell markers in osteosarcoma: prognostic determination and validation of CLK1 in disease progression

Abstract

Introduction: Osteosarcoma (OS) is a malignancy of the bone that mainly afflicts younger individuals. Despite existing treatment approaches, patients with metastatic or recurrent disease generally face poor prognoses. A greater understanding of the tumor microenvironment (TME) is critical for enhancing outcomes in OS patients.

Methods: The clinical and RNA expression data of OS patients were extracted from the TARGET database. The single-cell RNA sequencing (scRNA-seq) data of 11 OS samples was retrieved from the GEO database, and analyzed using the Seurat package of R software. Copy number variation (CNV) was analyzed using the InferCNV software. The potential interactions between the different cells in the TME was analyzed with the CellChat package. A multi-algorithm-based computing framework was used to calculate the tumor-infiltrating immune cell (TIIC) scores. A prognostic model was constructed using 20 machine learning algorithms. Maftools R package was used to characterize the genomic variation landscapes in the patient groups stratified by TIIC score. The human OS cell lines MG63 and U2OS were used for the functional assays. Cell proliferation and migration were analyzed by the EdU assay and Transwell assay respectively. CLK1 protein expression was measured by immunoblotting.

Results: We observed higher CNV in the OS cells compared to endothelial cells. In addition, there was distinct transcriptional heterogeneity across the OS cells, and cluster 1 was identified as the terminal differentiation state. S100A1, TMSB4X, and SLPI were the three most significantly altered genes along with the pseudo-time trajectory. Cell communication analysis revealed an intricate network between S100A1+ tumor cells and other TME cells. Cluster 1 exhibited significantly higher aggressiveness features, which correlated with worse clinical outcomes. A prognostic model was developed based on TIIC-related genes that were screened using machine learning algorithms, and validated in multiple datasets. Higher TIIC signature score was associated with lower cytotoxic immune cell infiltration and generally inferior immune response and survival rate. Moreover, TIIC signature score was further validated in the datasets of other cancers. CLK1 was identified as a potential oncogene that promotes the proliferation and migration OS cells.

Conclusion: A TIIC-based gene signature was developed that effectively predicted the prognosis of OS patients, and was significantly associated with immune infiltration and immune response. Moreover, CLK1 was identified as an oncogene and potential therapeutic target for OS.

Keywords: immune cell markers; immunotherapy; intratumor heterogeneity; osteosarcoma; prognosis.

Figure 1

Detailed classification of OS cells. (A) t-SNE plot displaying the origin of the sorted OS cells. The different samples are color-coded. (B) t-SNE plot displaying the OS populations. The different clusters are color-coded. (C) Correlation matrix of marker gene expression in each cell cluster. (D) t-SNE plot showing unique expression of different marker genes in each annotated cluster. (E) Bar plot indicating the composition of different cell types from individual patients. (F) t-SNE plot illustrating the annotation results of cell types. (G) CNV heatmap of OS cells (endothelial cells as the internal reference). (H) CNV scores of different OS cell clusters. (I) t-SNE plot displaying the classification of identified OS cells based on CNV scores.

Figure 2

Trajectory and intercellular communication analysis of OS cells. (A) Differentiation trajectories, pseudo-time distribution, cell cluster distribution along pseudo-time, and the proportion of each cluster for all OS cells. (B) Relative alteration in the expression of S100A1, SLP1, and TMSB4X along pseudo-time. (C) Quantity and strength of intercellular communication between S100A1+ OS cells, TMSB4X+ OS cells, SLPI+ OS cells, and other cell types. (D) Bubble plot illustrating the interaction between S100A1+ OS cells, TMSB4X+ OS cells, SLPI+ OS cells, and the different cell ligands and receptors. (E) Bubble plot illustrating the interaction between different cell types and S100A1+ OS cells, TMSB4X+ OS cells, SLPI+ OS cells. (F) Enrichment analysis of S100A1+ OS cells, TMSB4X+ OS cells, SLPI+ OS cells.

Figure 3

TF analysis of OS cells. (A) Volcano plot showing the top 5 highly expressed genes in each cluster. (B) Violin plots and UMAP plots of the top 5 upregulated genes. (C, D) Heatmaps displaying the distribution of TFs in the different clusters.

Figure 4

Functional analysis of Aggressive and EMT phenotypes. (A) Cell trajectory analysis showing the expression pattern of different identified TFs in various differentiation states. (B, C) Invasion levels of the three clusters shown in t-SNE plot (B) and violin plot (C). (D, E) EMT levels of the three clusters displayed in t-SNE plot (D) and violin plot (E).

Figure 5

Correlation of cell cluster with prognosis in OS. (A–C) Impact of abundance of clusters 0, 1, and 2 on survival. (D–F) Time-dependent ROC curves of clusters 0, 1, and 2. (G) Forest plot showing the results of multifactor Cox analysis.

Figure 6

Identification of TIIC-RNA at single-cell level. (A) t-SNE plot of identified TME cells and OS cells. (B) t-SNE plot of identified OS cells and 5 types of immune cells. (C) Violin plot showing differentially expressed genes in the identified immune cells. (D) t-SNE plot of identified immune cells and OS cells. (E) Volcano plot displaying differentially expressed genes between immune cells and OS cells. (F) Venn diagram classifying intersecting genes identified by six ML algorithms.

Figure 7

Construction of TIIC prognosis model. (A) Univariate Cox regression analysis of TIIC-related genes. (B–F) Dimension reduction of 22 prognostic genes using (B) CoxBoost algorithm, (C, D) LassoCox algorithm, and (E, F) random survival forest algorithm. (G) Venn diagram showing prognostic genes identified by all three ML algorithms. (H) Kaplan-Meier survival curves of OS patients of different TIIC feature scores in TARGET-OS and other validation datasets. (I) ROC curves of TIIC scores for predicting 1- to 5-year overall survival in the TARGET-OS and other validation datasets.

Figure 8

Comparison of the prognostic value of TIIC score and other prognostic models. (A) Circos plot showing different clinical factors in the TIIC-low and TICC-high groups. (B) C-index bar plot of TIIC score and various clinical factors in TARGET-OS and other validation datasets. (C–F) C-index plots of TIIC score and 42 prognostic models in TARGET-OS and other validation datasets.

Figure 9

Biological characteristics of the TIIC signature in TARGET dataset. (A) Results of GSVA based on MsigDB showing the biological properties associated with TIIC score. (B) t-SNE plots illustrating the differences in GO terms and KEGG pathways between TIIC-low and TIIC-high groups. (C) Enrichment analysis of differentially expressed genes between the TIIC-low and TIIC-high groups based on Metascape. (D) GSEA results depicting the enrichment of GO and KEGG terms between the TIIC-high and TIIC-low groups.

Figure 10

Immunological features of the TIIC signature in TARGET dataset. (A, B) Relationship between the TIIC score, immune infiltrating cells and immune regulatory genes. (C) Abundance of associated pathways in high TIIC group and low TIIC group.

Figure 11

Prediction of the immunotherapeutic response based on TIIC signature scores. (A) Survival analysis of the IMvigor cohort based on TIIC scores. (B–D) Correlation between TIIC score and the immunotherapeutic response in the (B) IMvigor, (C) GSE179351, and (D) GSE35640 datasets. (E) Survival analysis of the Braun dataset based on TIIC scores. (F–H) Correlation between TIIC score and the immunotherapeutic response in the (F) Braun, (G) GSE91061, and (H) GSE103668 datasets. (I) Survival analysis of the Nathanson dataset based on TIIC scores. (J) Correlation between TIIC score and immune therapeutic response in the Nathanson dataset. (K) Survival analysis of the GSE78220 dataset based on TIIC scores. (L–O). Correlation between TIIC score and immune therapeutic response in the (L) GSE78220, (M) GSE126044, (N) GSE165252, and (O) TARGET datasets.

Figure 12

Metabolic characteristics associated with TIIC scores in the TARGET dataset. (A) Results of GSVA based on KEGG pathways for 11 metabolic categories in the TIIC score groups. (B) Differences in metabolic pathways between the TIIC-high and TIIC-low groups. (C) Correlation between TIIC feature scores and KEGG pathways in GSVA.

Figure 13

Mutation landscape in the TARGET dataset. (A) Waterfall plot of the top 50 mutated genes in the TARGET dataset. (B) Mutation landscapes of OS patients grouped by TIIC score. (C) Exclusive and co-occurring mutations in the OS patients with different TIIC scores. (D) Distribution of CNVs in the OS patients stratified by TIIC score, with FGA, FGG, and FGL as features.

Figure 14

CLK1 promotes OS proliferation and migration. (A) Immunoblot showing CLK1 protein expression in the OS tissues. (B) Colonies formed by the CLK1-overexpressing and CLK1-knockdown MG63 and U2OS cell lines. (C) The up-regulation and down-regulation of EDU staining of the MG63 and U2OS cell lines under the condition of CLK1 overexpression and knockout, reflecting the proliferating capability of the OS cells. (D) Flow cytometry assessing the MG63 and U2OS cell lines in the phase of G2/M state and corresponding quantification results. (E) Transwell analysis assessing the migration of the MG63 and U2OS cell lines under the condition of CLK1 overexpression and knockout and the corresponding quantification results. (*, denotes a significance level of p < 0.05; **, p < 0.01; ***, p < 0.001).