3.1 A holistic landscape of immunological features
Based on the tumor stroma and immune characteristics of osteosarcoma, the acquired osteosarcoma samples were analyzed using the ESTIMATE algorithm to reveal the level of immune infiltration of tumor samples in osteosarcoma. StromalScore, ImmuneScore, ESTIMATEScore, and TumorPurity were calculated and analyzed in the tumor samples. The scores were grouped according to their median scores, which were named as high and low groups, and Kaplan-Meier curves were performed according to the survival time and status, to explore the relationship with tumor prognosis(Fig. 1A-D). In StromalScore, ImmuneScore, and ESTIMATEScore, high grouping was significantly associated with high survival, while low TumorPurity was associated with high survival. Of these, ImmuneScore was most significantly associated with survival. The data obtained from the above suggests that osteosarcoma can be analyzed prognostically based on ImmuneScore, followed by ssGSEA analysis.
【Figure 1】
The ESTIMATE algorithm was utilized to conduct survival analysis on patient groups categorized based on high and low immune scores, and the results were visually represented through Kaplan-Meier (K-M) curves. (A) ESTIMATE score, (B) Tumor purity score, (C) Stromal score, and (D) Immune score.
The impact of immunity on osteosarcoma was explored using ssGSEA. Subsequently, unsupervised clustering was performed on the ssGSEA results, and the number of clusters with the highest average within-group concordance was 2(Fig. 2B). Therefore, based on their immune characteristics, the samples were categorized into two groups, and the 30 immune cells were subsequently displayed along with clinical characteristics including gender, age, survival time, survival status, and whether the tumor has metastasized or not(Fig. 2A).There are noticeable variations in immunological features observed between the clustered groups.Groups with rich immune profiles were named high immune groups, therefore, groups with lower immune profiles were named low immune groups.To further investigate the relationship between immune groupings and the levels of TME, immune activation, and tumor cell infiltration, several analyses of osteosarcoma TME were performed. As expected, StromalScore, ImmuneScore, and ESTIMATEScore were higher and TumorPurity was lower in the hyperimmune group(Fig. 2C-F), suggesting that the hyperimmune group was associated with higher survival. Immune checkpoints and immune cell abundance were also significantly higher in the hyperimmune group(Fig. 2G).Subsequent Principal Component Analysis (PCA) was performed and the immune characteristics differed significantly between the two groups(Fig. 2H), therefore, we hypothesized that the immune groupings found above could well distinguish the immune and genetic characteristics of the samples.
【Figure 2】
Immune subtype identification and comparative analysis. The symbols *represent p-values less than 0.05,**represent p-values less than 0.01,*** represent p-values less than 0.001,**** represent p-values less than 0.0001. (A). ssGSEA analyses of 84 samples, divided into two groups based on 30 ssGSEA scores per sample. (B). Unsupervised clustering of the samples based on their immunological characteristics, where the number of clusters with the highest average within-group agreement is 2. (C-F).C-F shows, in order, the differences in ESTIMATEScore, ImmuneScore, StromalScore, and TumorPurity between the high and low immunity groups.(G). A box plot is utilized to display the levels of immune cell infiltration in groups categorized as either having high or low immunity. In this visualization, the red boxes correspond to the high immunity group, while the blue boxes correspond to the low immunity group.(H). PCA analysis of the two immune subtypes, with purple and yellow dots representing the immunity − high and − low groups, respectively.
3.2 Mutation analysis
The analysis of gene mutations revealed a notably elevated mutation rate in the low immunity group, with TP53 being identified as the gene exhibiting the highest frequency of mutations(Fig. 3A). Tumor Mutation Burden (TMB) analysis was carried out and visualized using scatter plots(Fig. 3B).Following this, an analysis of TMB in the high and low-immunity cohorts demonstrated a variation in TMB levels between the two groups, with TMB showing an elevation in the low-immunity cohort(Fig. 3C). We hypothesized that this may be related to the poorer prognosis of the low-immunity group.
【Figure 3】
Mutations between high and low immunity groups. The symbols**represent p-values less than 0.01. (A). Mutation status of genes in the high and low immunity groups. (B). TMB distribution of all samples. (C). Bar graph showing the difference in TMB between the high and low immunity groups.
3.3 Screening and evaluation of differentially expressed genes
The examination of gene expression variations between immune subtypes identified a pool of 836 differentially expressed genes for subsequent scrutiny, with 697 genes showing up-regulation (83.37%) and 139 genes demonstrating down-regulation (16.63%)(Fig. 4A).Enrichment analyses of DEGs by GO annotation and KEGG pathway enrichment analyses identified 497 biological processes (BPs), 94 cell components (CCs), 118 molecular functions (MFs) and 69 KEGG pathways(Fig. 4C,D).The DEGs were also analyzed for GSEA enrichment(Fig. 4B).
【Figure 4】
Enrichment analysis was conducted on the DEGs identified in the two immune subtypes. (A).Volcano diagram showing the regulation of DEGs expression, with green, grey, and red dots representing down-regulation, unregulation, and up-regulation, respectively. (B).Bubble diagram showing the top ten pathways according to Gene Set Enrichment Analysis (GSEA). (C).Bubble plots showing the top 10 enriched GO BP, CC,and MF. (D).Bubble plots showing the top 20 enriched terms of the KEGG pathway, with the size of the dots representing the number of enrichments.
As shown in the figure, the top ten biological processes were filtered according to P-value, most of which were related to immunity, including immune response, inflammatory response, innate immune response, neutrophil chemotaxis, and adaptive immune response. in addition, KEGG pathways are also related to immunity, including Phagosome, Antigen processing and presentation, B cell receptor signaling pathway and Th17 cell differentiation, etc.In summary, DEGs and immunity are closely related, subsequently, GSEA enrichment analysis of DEGs was performed to further investigate the pathway differences between immune subtypes, based on P < 0.05, a total of 47 enriched pathways were obtained, and the top ten were Tuberculosis (NES = 3.11,p.value = 6.73E-10) in order, Leishmaniasis (NES = 3.01,p.value = 3.03E-09), Staphylococcus aureus infection (NES = 3.01,p.value = 7.05E-09), Th17 cell differentiation (NES = 2.90,p.value = 4.69E-08), Th1 and Th2 cell differentiation (NES = 2.82,p.value = 1.13E-07), Hematopoietic cell lineage (NES = 2.72,p.value = 5.33E-07), Influenza A (NES = 2.71,p.value = 4.87E-07), Systemic lupus erythematosus (NES = 2.60,p.value = 1.83E-06), Phagosome (NES = 2.66,p.value = 2.21E-06) and Neutrophil extracellular trap formation (NES = 2.46,p.value = 9.22E-06). These findings imply that the activation of the immune system in the tumor microenvironment is implicated in the development of osteosarcoma.
3.4 Screening and Evaluation of IR-DEGs
From the TIMER database, 1811 immune-related genes were obtained, and these genes were intersected with DEGs to obtain 221 IR-DEGs(Fig. 5A).The IR-DEGs were entered into the STRING online database to obtain the protein- protein interaction network(Fig. 5B). The IR-DEGs were subjected to GO annotation and KEGG pathway enrichment analyses(Fig. 5C,D), resulting in 418 BP, 57 CC, 57 MF and 69 KEGG pathways.These IR-DEGs were mainly enriched in biological pathways such as immune response, inflammatory response, and adaptive immune response. Among KEGG-enriched pathways, the top five were Cytokine-cytokine receptor interaction, Viral protein interaction with cytokine and cytokine receptor, Rheumatoid arthritis, Chemokine signaling pathway and Graft-versus-host disease.GSEA enrichment analyses were also performed for IR-DEGs(Fig. 5E). Their top five were Th17 cell differentiation (NES = 2.56, p.value = 1.35E-05), Th1 and Th2 cell differentiation (NES = 2.46,p.value = 3.33E-05), Staphylococcus aureus infection (NES = 2.30,p.value = 0.0002), Tuberculosis (NES = 2.09,p.value = 0.0005), Leishmaniasis (NES = 2.13,p.value = 0.001). Similar to the above results. According to the results, it is clear that immune activation, especially T cells, is important for the development of osteosarcoma.
【Figure 5】
Identification and enrichment analysis of differentially expressed genes (DEGs) associated with the immune system. (A) Venn diagram showing 221 immune-associated DEGs overlapping 836 DEGs and 1811 IRGs.(B) Protein interaction network diagram of DEGs (C).Bubble plots showing the top 10 enriched GO BP, CC, and MF. (D). .Bubble plots showing the top 20 enriched terms of the KEGG pathway, with the size of the dots representing the number of enrichments.(E) Bubble plots showing the top 20 enriched terms of the KEGG pathway, based on GSEA analyses of the top 10 pathways with the highest gene enrichment.
Following that, an examination of the connections and relationships among these IR-DEGs was carried out through the analysis of the protein-protein interaction network,and associations were found for a variety of IR-DEGs, with the most significant correlations between IL6, IL10, CD4, CD8A, IL1B, TNF, and CCL5 and the other immune IR-DEGs.
The interaction network is characterized by the presence of Interleukin (IL) family genes (including IL10, IL1B and IL6), T-Cell Surface Glycoprotein genes (including CD8A, CD86 and CD4), C-C Motif Chemokine Receptor genes (including CCR7 and CCR5) and C-C Motif Chemokine Ligand genes (including CCL2 and CCL5), which are among the hub nodes. Subsequently, the MCC algorithm was applied to obtain the top ten genes, and prognostic survival analysis and immune infiltration analysis were performed on these ten hub genes(Fig. 6,7).
The results showed that ten genes were closely associated with immune cells, and among all hub genes, Tumor Infiltrating Lymphocytes (TIL) infiltration levels were the highest, and mast cells and dendritic cell infiltration levels were lower. In the K-M survival analysis of the hub genes, CD4, CD8A, CCR5, and CCL5 were prognostically significant.
【Figure 6】
Analysis of immune cell infiltration was conducted on the ten hub genes, visualized using lollipop charts.From left to right, the order is CCR5, TNF, IL10, IL6, CD8A, CD4, IL1B, CCR7, CCL5, and CCL2.
【Figure 7】
Ten hub genes were subjected to survival analysis.
3.5 Risk modelling
Twenty-six genes were screened from 221 IR-DEGs using one-way COX analysis based on p < 0.05, and then the 26 genes were further analyzed by applying LASSO algorithm using ten-fold cross-validation when lambda.1se = 0.1745826(Fig. 8A,B), a risk model was developed using two candidate genes, namely PDK1 and PPARG.Utilizing the expression levels and coefficient values of the candidate genes, a risk score model was created through the application of the subsequent formula: risk score = (0.6839441 × PDK1 expression) + (-0.6420120 × PPARG expression). The risk score formula was utilized on all samples to calculate the risk score for each sample. Subsequently, The samples were divided into two cohorts at random, with no statistical distinction observed between the pair.K-M curve analysis and ROC curve analysis were then performed on the two cohorts and the original combined cohort(Fig. 8C-H).ROC curves with area under the curve (AUC) values greater than 0.5 were considered statistically different. The results showed that the p-value of the K-M curve was less than 0.05 for the training set, the test set, and the merged set, so this feature was considered to have prognostic value. Whereas, in all three cohorts, the ROC curves indicated that the AUC values for 1-year, 3-year, and 5-year were above 0.5,and the feature had a higher predictive sensitivity for 3 and 5 years (AUC values were greater than 0.7 for both 3 and 5 years). In all three cohorts, a higher survival advantage was demonstrated for low-risk scores, and low-risk scores were strongly correlated with high survival.PDK1 and PPARG were analyzed according to the risk score grouping, and it was found that the expression of PDK1 was positively correlated with the risk score, while the expression of PPARG was inversely correlated with the risk score(Fig. 8I). Then, the differences in TMB between the high and low-risk score groups were compared based on the previously calculated TMB values, and as expected, there was a difference between the two groups (P < 0.05)(Fig. 8J), with higher TMB values in the high-risk score group, which may be related to the poorer prognosis of the high-risk score group.
【Figure 8】
Construction of the total risk profile.The symbols *represent p-values less than 0.05,*** represent p-values less than 0.001. (A, B) Employing the LASSO method for the identification of key candidate genes.(C-E) Survival analysis K-M curves for the training cohort, validation cohort, and initial combined cohort, respectively. (F-H) ROC curves were generated to assess the prognostic significance of risk features in the training, validation, and original merged groups.(I) Differences in PDK1 and PPARG gene expression between high and low-risk score groups. (J) Differences in TMB between high and low-risk score groups.
3.6 Performance of risk models
The samples were categorized based on the median risk score, resulting in the formation of a high-risk group and a low-risk group. Subsequently, the variations in immune checkpoint-associated genes between these two groups were analyzed(Fig. 9A). The results showed that most of the immune checkpoint-associated genes were differentially expressed, and of the genes that were differentially expressed, all had higher gene expression values in the low-risk score group. Subsequently, the two candidate genes were grouped according to expression(Fig. 9B,C), and most of the genes were differentially expressed. high expression of PDK1 was positively correlated with low expression of immune-related genes, while high expression of PPARG was positively correlated with high expression of immune-related genes. The two candidate genes were shown to have opposite roles. Subsequently, CIBERSORT immune infiltration analysis was performed(Fig. 9D), which showed differences in five immune cell subpopulations. Subsequently, the two candidate genes were grouped according to expression, and differential analysis of the five immune cell subpopulations obtained by CIBERSORT immune infiltration showed that high expression of PDK1 was positively correlated with high expression in Macrophages M0 and inversely correlated with high expression in Macrophages M2, and the opposite was true in PPARG(Fig. 9E,F). Once again, the two candidate genes were shown to have opposite roles. Finally, the samples were analyzed using the xCell online website to obtain the MicroenvironmentScore, which was significantly different between the high-risk scoring and low-risk scoring group and was higher in the low-risk scoring group(Fig. 9G).
【Figure 9】
Role of risk models for immune checkpoints.The symbols *represent p-values less than 0.05,**represent p-values less than 0.01,*** represent p-values less than 0.001. (A)Comparison of the expression levels of immune checkpoint-related genes in the group with high-risk scores and the group with low-risk scores. (B, C) The expression levels of immune checkpoint-related genes between the two groups by dividing the gene expression of PDK1 and PPARG into high and low groups according to the median. (D) CIBERSORT immune infiltration analysis of all samples. (E, F) Infiltration abundance levels of some immune cells by dividing gene expression of PDK1 and PPARG into high and low groups by median. (G) MicroenvironmentScore in the high and low-risk scoring groups.