| Objective:Osteosarcoma is the most common primary malignant bone tumor in children and adolescents[1].Despite its low incidence,it has a high degree of malignancy and a 5-year survival rate of 60-70%.The main clinical manifestations are bone and joint pain and localized mass.The pathological features were mainly neoplastic osteogenesis.Chemotherapy resistance and long-term lung metastasis were the leading causes of death.The 5-year survival rate for patients with recurrence or distal metastases is less than 20%[1].In particular,there has been no significant progress in the clinical treatment of osteosarcoma in the last 30 years.Currently,surgery in combination with neoadjuvant and adjuvant chemotherapy,with a focus on cytotoxic drugs,is the primary treatment[2].Immunotherapy,as a new tumor treatment method,has achieved remarkable curative effect in a variety of malignant tumors.However,osteosarcoma is a highly heterogeneous tumor that lacks the expression of immune-related neoantigens.Although there has been evidence suggesting the feasibility of immunotherapy for osteosarcoma,there have been only a few studies of immunotherapy for osteosarcoma to date,which is far from sufficient to meet the needs of clinical diagnosis and treatment.By using biomedical data,combined with cutting-edge molecular pathology and bioinformatics analysis tools,a comprehensive analysis of the genomic information of osteosarcoma,multi-dimensional mining of osteosarcoma-related multi-gene data,combined with molecular biology experiments for verification,thus building an immune risk assessment model of osteosarcoma,which has positive significance for predicting the survival of patients and the response to immunotherapy.Aimed at:1.Using bioinformatics combined with large-scale data analysis,a set of genes with abnormal expression in osteosarcoma compared to normal bone tissue were screened and two key genes were selected by combining PPI,KM curve,ROC curve,physical and chemical properties and other factors.A two-gene model was constructed to predict the prognosis using these two genes and the survival risk was assessed.2.The ICD,GSEA,ESTIMATE,TIDE,R language,IMvigor210 and other biogenetic methods were used to analyze the relationship between immuno-dual gene models and immune scores,immune checkpoints,immune therapy responses,and to screen the immune checkpoints most relevant to the immuno-dual model.3.The CIBERSORT,R language,and TISCH databases were used to explore immune cells associated with the dual gene model of immunity,and single cell analysis of key immune cells in osteosarcoma was conducted to explore the role of immune cells in the development of osteosarcoma.Methods:1.Screening of immune core genes in osteosarcoma and the construction of multi-gene models for immune prediction.①Expression data for osteosarcoma were downloaded from the Gene Expression Omnibus(GEO)database.(Screening criteria:Sample size>5 cases and Normal control tissue).②The pheatmap package in the R was used to analyze the differences in expression of the ICD signature genes between osteosarcoma and normal tissue.③The differentially expressed genes(DEGs)between osteosarcoma and normal specimens were identified.All differentially expressed genes were grouped according to their differentially expressed multiples.The results of the GO/KEGG analysis were used to select appropriate groups of differentially represented multiples.④All genes in the group were intermingled with genes in the IMMPORT immune database,and STRING was used to construct PPI network for intermingled genes,so as to explore the enrichment pathway of intermingled genes.⑤The core immune-related genes in osteosarcoma were obtained using LASSO.KM survival analysis was used to analyze the effect of individual genes on survival,and permutations and combinations of multiple genes were performed.The ROC curve was employed to pick the appropriate two genes.⑥A multivariate COX regression analysis was used to construct the immuno-dual gene model using the R of survival package.⑦Univariate and multivariate COX analysis,KM and ROC were used to verify the relevance of the model to survival outcomes.The nomogram package of R software was used to predict patient survival and explore the relevance of clinical features of the TARGET database to the immuno-dual gene model.⑧The solid tissues of osteosarcoma in this model were validated by the utilization of IHC to express immune checkpoints such as TIM-3,LAG-3 and PD-L1.2.Correlation analysis of immuno-dual gene model and microenvironmental immunometric indicators in osteosarcoma.①The pheatmap package in the R was used to analyze differences in expression of ICD marker genes between high and low risk groups.②ESTIMATE was used to predict the relationship between the model and immunity,including immune score,stromal score,and tumor purity.③GSEA was used to explore the pathway enrichment of the model.④SPSS was used to explore the correlation between the model and immune checkpoints including PD-L1,LAG3,TIM3,TIGIT and CTLA4 in different datasets.⑤IHC was used to validate the expression of the immune checkpoints such as TIM-3,LAG-3 and PD-L1 in the solid tissues of osteosarcoma in this model.⑥Methods such as IMvigor210 and TIDE have been used to predict a patient’s response to immunotherapy through the immuno-dual gene model.3.Correlation analysis between osteosarcoma immuno-dual gene model and immune cells and exploration of single cells.①CIBERSORT,Spearman and Wilcoxon tests were used to screen for model-related immune cells,and KM survival analysis was used to screen the immune cells.②IHC experiments were used to validate the expression of key immune cells in the solid tissues of the osteosarcoma model.③Appropriate single-cell datasets were selected from the GEO database and downloaded to obtain single-cell sequencing data.④After cell quality control,sample batch removal and other processing using the R software package seurat,the single-cell data was divided into multiple cell clusters,and the cell clusters were annotated by the expression of TOP 10 genes in the TISCH database.⑤The R software package pheatmap was used to analyze the correlation between the risk score of the immune model and the fraction of each cell cluster in the sample.⑥The R software package seurat was used to classify the myeloid cell population,pseudo-temporal trajectory analysis was used to explore the differentiation and maturation processes of immune cells by package monocle3,and the role of key immune cells in the development of osteosarcoma was combined with this analysis.Results:1.Screening of immune core genes in osteosarcoma and the construction of multi-gene models for immune prediction.①Data screening:In the GEO database,four osteosarcoma datasets with suitable conditions were screened:GSE12865,GSE42352,GSE16088,and GSE28424.Among them,GSE12865 dataset contained 14 samples,including 12 osteosarcoma samples and 2 normal tissue samples.The GSE42352 dataset contains 87 samples,including 84 osteosarcoma samples and 3 normal tissue samples.The GSE16088 dataset contains 20 samples,including 14 osteosarcoma samples and 6 normal tissue samples.The GSE28424 dataset contains 23 samples,including 19 osteosarcoma samples and 4 normal tissue samples.②The results showed that the expression of 33 ICD marker genes varied between osteosarcoma and normal tissue,with 21 of them being lower in osteosarcoma tissue.③According to the differential gene screening criteria:①p<0.05;(2)|log2FC |≤1.0(less than two differentially expressed multiples),1.0≤| log2FC |≤2.0(less than four difference ratio greater than two),2.0≤| log2FC |≤3.0(less than eight difference ratio greater than four),3≥|log2FC |(differentially expressed ratio greater than eight)will be differentially expressed genes into four groups.The GO/KEGG pathway analysis showed that genes in the fourth group with a difference of more than eight times were associated with some immune pathways.To interlink the 811 genes in the fourth group with IMMPORT immunology database revealed 111 immune related hub genes.PPI network interaction map shows that 108 genes are related to each other and are involved in signal reception,tumorigenesis and immune processes.④Six core genes(STC2,FPR1,VCAM1,PTN,IL13RA2 and HCST)were screened out by LASSO,which were significant in KM survival except for STC2 and IL13RA2.COX regression model was constructed by permutation and combination of the six genes.The ROC curve showed that the STC2 gene has a higher AUC when combined with other genes,and there was almost no difference between the three-gene combination model and the two-gene combination model.STC2 and FPR1 genes were selected by combining the physical and chemical properties of the six gene expression proteins.⑤The expression levels of STC2 and FPR1 genes were obtained in TAREGT database,and the immuno-dual gene model was constructed by R software survival package and multivariate COX regression:Risk Score=(0.5663)*STC2+(-1.3968)*FPR1,and the study population was further divided into high and low risk groups according to the median.⑥KM was used to analyze the prognosis of patients with osteosarcoma with high and low risk scores in the TARGET database,and the results showed that the higher the score,the worse the prognosis.AUC was used to judge the reference value of the model.The AUC values of 1,3 and 5-year survival were all greater than 0.7,and the 5-year survival value was 0.857,indicating that the model is more reliable.⑦Univariate and multivariate COX analyses showed that only immuno-dual gene models and metastasis could be considered as independent prognostic factors.These two metrics were used to construct a nomogram to predict patient survival at 1,3 and 5 years.There was no significant difference in the distribution of clinical indicators between the high and low risk groups.⑧The qPCR results showed that STC2 was more expressed in osteosarcoma patients,while FPR1 was more expressed in normal individuals.2.Correlation analysis of immuno-dual gene model and microenvironmental immunometric indicators in osteosarcoma.①The expression of the ICD marker gene was higher in the low-risk group than in the high-risk group.②The GSEA results showed that the enriched pathways in the high-risk group were related to glutamate and tooth growth,and did not include immune pathways.In the low-risk group,most of the enriched pathways were associated with immunity.③The ESTIMATE results showed that the risk score of the model was negatively correlated with the immune score and positively correlated with the tumor purity.That is,a population with a higher risk score was associated with a lower immune score and a higher tumor purity.④The correlation between the model and PD-L1,LAG-3,TIM-3,TIGIT,CTLA4 and other immune checkpoints indicated that the model was more likely to be negatively correlated with immune checkpoints,with TIM-3 having the highest correlation coefficient with the model and the strongest correlation.Samples with higher expression of TIM-3 had a better prognosis.⑤PD-L1 was barely expressed in osteosarcoma tissue,while TIM-3 and LAG-3 were more expressed in high-grade osteosarcoma than in low-grade osteosarcoma.⑥The immuno-dual gene model based on the IMvigor210 database were used to predict how patients would respond to immunotherapy,showing that the low-risk group had higher TMB values and were more likely to respond to immunotherapy.Two genes in the model were verified to be involved in immune response and tumor mutation burden.Similar results were obtained in the TIDE,namely that immunotherapy may be more effective in the low-risk group,and that CD8/PD-L1 expression was found to be higher in the low-risk group.3.Correlation analysis between osteosarcoma immuno-dual gene model and immune cells and exploration of single cells.①CIBERSORT was used to explore the expression of 22 immune cells in high and low-risk populations.The results showed that:there were significant differences in naive B cells,memory B cells,CD8+T cells,non-activated CD4+memory T cells,M1 macrophages,M2 macrophages,activated mast cells and neutrophils.The proportion of naive B cells,CD8+T cells,activated mast cells,and memory B cells was higher in the high-risk group.While the fraction of activated CD4+memory T cells,M1 macrophages and M2 macrophages was low.Naive B cells,memory B cells,M1 and M2 macrophages,activated mast cells,and neutrophils were selected by combining the results of Spearman’s and Wilcoxon test analyses of the immuno-dual gene model and the immune cell.Analysis of immune cell and KM survival outcomes showed that higher expression of activated mast cells and lower expression of M2 macrophages were associated with poorer survival.Considering the role of M2 macrophage specificity in osteosarcoma tissues,M2 macrophages are considered as key immune cells for further investigation.②CD 163 and CD68 were extensively expressed in osteosarcoma.CD68 was primarily expressed on osteoclast-like giant cells,and CD 163 was primarily expressed on monocyte-like macrophages.③The osteosarcoma dataset GSE162454 was selected from the GEO database,which contained six single-cell samples.④After the single cell data were processed by seurat package of R,the magnet genes were screened and obtained under the condition of 200<RNA<6000,and 20 cell clusters were obtained.Based on the expression of TOP10 genes such as IFNG,CTSL,NKG7 in TISCH,the 20 cell clusters were divided into eight cell populations:malignant cells,myeloid cells,CD4+T cells,CD8+T cells,osteoblast,fibroblast cells,plasma cells,and endothelial cells.⑤The proportion of cell populations in each sample was analyzed using the pheatmap package of R,and it was found that malignant and myeloid cells had relatively high infiltration in all samples.The risk score of the immuno-dual gene model was positively correlated with the infiltration fraction of the malignant cells and negatively correlated with that of the myeloid cells.⑥The seurat package of R reclassified and annotated the myeloid cells and obtained seven subsets of M1 macrophages,M2 macrophages,CD14+monocytes,CXCL8+monocytes,dendritic cells,neutrophils,and NKG7+NK/T cells.Pseudo-time trajectory analysis showed that when monocytes differentiated into M2 macrophages,the expression of STC2 was decreased,while the expression of FPR1 and TIM-3 was increased.Conclusion:1.The immune core genes of osteosarcoma were screened to construct a multi-gene immune prediction model.A immuno-dual gene model was constructed by screening for differential immune-related genes in osteosarcoma.Riskscore=(0.5663)*STC2+(-1.3968)*FPR1,this immune model can be used as an independent prognostic risk indicator.2.Correlation analysis between the immuno-dual gene model of osteosarcoma and immunometric indicators of the microenvironment and their clinical significance.The immuno-dual gene model was negatively correlated with the expression of ICD marker genes,immune score and immune checkpoint.The low-risk group is enriched in multiple immune pathways and may respond better to immunotherapy.The immuno-dual gene model can be used to assess the universal expression of the immune checkpoint TIM-3 in osteosarcoma.3.Correlation analysis between the immuno-dual gene model of osteosarcoma and immune cells and exploration of single-cell sequencing.The immune model was negatively correlated with M2 macrophage infiltration.The potential role of M2 macrophages in the osteosarcoma microenvironment was preliminarily analyzed by single-cell sequencing. |