Single-cell RNA-seq reveals T cell exhaustion and immune response landscape in osteosarcoma

Background: T cell exhaustion in the tumor microenvironment has been demonstrated as a substantial contributor to tumor immunosuppression and progression. However, the correlation between T cell exhaustion and osteosarcoma (OS) remains unclear.

Methods: In our present study, single-cell RNA-seq data for OS from the GEO database was analysed to identify CD8+ T cells and discern CD8+ T cell subsets objectively. Subgroup differentiation trajectory was then used to pinpoint genes altered in response to T cell exhaustion. Subsequently, six machine learning algorithms were applied to develop a prognostic model linked with T cell exhaustion. This model was subsequently validated in the TARGETs and Meta cohorts. Finally, we examined disparities in immune cell infiltration, immune checkpoints, immune-related pathways, and the efficacy of immunotherapy between high and low TEX score groups.

Results: The findings unveiled differential exhaustion in CD8+ T cells within the OS microenvironment. Three genes related to T cell exhaustion (RAD23A, SAC3D1, PSIP1) were identified and employed to formulate a T cell exhaustion model. This model exhibited robust predictive capabilities for OS prognosis, with patients in the low TEX score group demonstrating a more favorable prognosis, increased immune cell infiltration, and heightened responsiveness to treatment compared to those in the high TEX score group.

Conclusion: In summary, our research elucidates the role of T cell exhaustion in the immunotherapy and progression of OS, the prognostic model constructed based on T cell exhaustion-related genes holds promise as a potential method for prognostication in the management and treatment of OS patients.

1 Introduction

Osteosarcoma (OS) is the most prevalent aggressive bone tumor occurring in children and adolescents, constituting the majority of all bone tumor cases (1). The conventional approach to OS entails a blend of surgery and rigorous multi-agent chemotherapy. However, the prognosis for OS patients remains exceedingly grim, primarily attributed to delayed diagnosis and early-stage distant metastasis (2). Therefore, it is crucial to explore innovative and efficacious therapies aimed at enhancing the prognosis of OS patients.

Immunotherapy holds significant promise in the treatment of malignant tumors in humans, numerous recent studies have highlighted its considerable potential in tumor therapy, with preclinical trials providing robust support (3). Immune checkpoint inhibitors (ICIs) exhibit considerable potential for immunotherapy in OS. They can navigate the genomic complexity of OS, leading to enhanced overall outcomes (4). Although immunotherapy for OS has demonstrated promising therapeutic effects in some studies, it has yet to substantially improve patient prognosis (5, 6). In clinical trials of OS, the response to ICIs has not been favorable and trial results are not yet satisfactory (7). This may be attributed to the immune microenvironment in OS, which suppresses T cell function (8). Immune cells constitute the cellular foundation of immunotherapy, of which CD8⁺ T cells serving as a pivotal component of cancer immunotherapy (9). Activated CD8⁺ T cells mature into cytotoxic T lymphocytes (CTLs) and represent a key component of the immune system’s antitumor response, CTLs are associated with increased survival rates in various types of cancer and play a crucial role in immune surveillance, targeting and eliminating cancer cells (10). The optimal approach for achieving tumor eradication will likely entail a combination of therapies that promote immune activation and T cell initiation, counteract immunosuppressive signals in the tumor microenvironment, and sustain the presence of T cells in cancerous tissue.

T cell exhaustion entails a progressive, hierarchical, and negatively regulatory process affecting T cells within the tumor microenvironment (11). Classical inhibitors targeting PD-1 and CTLA-4 largely exert their anti-tumor effects by mitigating functional exhaustion (12). However, the precise underlying mechanisms of these inhibitors necessitate further investigation. The recent advancement of biomarkers unveiled potential molecular regulatory targets for CD8⁺ T cells in the intricate tumor heterogeneity of OS. Moreover, the potential correlation between alterations in exhaustion expression profiles and immune checkpoints has presented avenues for research (13).

In this study, we aim to delve into potential molecular regulatory targets and core regulatory genes associated with T cell exhaustion in the intricate tumor heterogeneity of OS. We developed a multi-biomarker model based on genes linked to T cell exhaustion, which functions in evaluating the tumor microenvironment, predicting immunotherapy response, and forecasting the prognosis of diverse OS patients. It has great potential to play a vital role in guiding clinical practice in the future.

2 Methods

2.1 Obtaining the raw data

The single-cell sequencing data (GSE162454), along with microarray data (GSE16091 and GSE21257) pertaining to OS, were acquired from the GEO database (http://www.ncbi.nlm.nih.gov/geo). Additionally, data from 84 distinct OS patients’ samples were retrieved from the TARGETs database. All datasets were accompanied by clinical and prognostic information, which was employed for subsequent analyses.

2.2 Data processing of single-cell RNA sequence

Data analysis and quality assessment were conducted using the R package “Seurat” (version 4.3.0; http://satijalab.org/seurat/). Cells with expression of fewer than 250 genes or with a percentage of mitochondrial genes exceeding 20% of the total expressed genes were excluded from the analysis. Additionally, cells with unique molecular identifiers (UMI) resulting in log10(UMI) > 0.80 were also removed. Subsequently, potential doublets were identified and eliminated using the R package “DoubletDecon” (version 1.1.6; http://EDePasquale/DoubletDecon).

2.3 Data integration and dimensionality reduction

The feature counts for each cell underwent a transformation, involving division by the total counts for that cell, followed by multiplication by 10,000. Subsequently, the results were logarithmically transformed and then normalized by adding 1, thus preventing the computation of the logarithm of 0. Before proceeding with the normalization of the expression matrix, the top 2000 highly variable genes (HVGs) were identified, centered, and scaled. Subsequent to this, a principal component analysis (PCA) was conducted based on these HVGs. Following that, the R package “Harmony” (https://github.com/immunogenomics/Harmony) was employed to integrate the cellular data from six samples and mitigate any potential batch effects.

2.4 Cell-clustering and annotation

The clustering analysis relied on the embedding of the Harmony algorithm, executed through the “FindNeighbors” and “FindClusters” functions within the “Seurat” package. The resulting clusters were visualized on a two-dimensional plot generated via the UMAP method. For subcluster analysis, akin procedures were applied, encompassing variable gene identification, dimensionality reduction, Harmony for cell integration, and cluster identification for the distinct clusters derived from the overall analysis. The annotation of clusters was performed using established cellular markers drawn from the literature. Detailed information regarding the cellular markers can be found in the Supplementary Tables.

2.5 Identification and analysis of CD8⁺ T cell subtypes

CD8⁺ T cells were isolated and subsequently re-clustered using the “Seurat” package in R. Single-cell pseudotime trajectories were constructed employing the “Monocle2” package in R. Following this, a weighted correlation network analysis (WGCNA) was conducted to identify the core gene sets within CD8⁺ T cellular clusters using the “hdWGCNA” package in R. To explore intercellular communication between all cell clusters, the R package “Cellchat” was utilized. The levels of immune checkpoint molecules between clusters were assessed based on the immune checkpoint expression profile. Differential functional status regarding Gene Ontology (GO) and KEGG pathways for each CD8⁺ T cell cluster were analyzed using the “ClusterProfiler” R package. Additionally, the GSEA pathways obtained from MSigDB (gsea-msigdb.org) were evaluated using the “fgsea” R package. Furthermore, differences in HALLMARK pathways between the clusters were determined through gene set variation analysis (GSVA) using the “GSVA” R package.

2.6 Construction and validation of the T cell exhaustion signature

The CD8⁺ T cell exhaustion genes with prognostic potential in the TARGETs dataset were identified through Univariate Cox regression analysis (P<0.05). Subsequently, a combination of six machine learning algorithms was employed, which included the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm (14), Boruta feature selection algorithm (15), survival support vector machine (survival-SVM) based on 10-fold cross-validation (16), Boosting in Cox regression (Cox-boost) (17), Extreme Gradient Boosting(XG-boost) (18), and generalized boosted regression modeling (GBM) (19), to further refine the valuable T cell exhaustion signature. In constructing the model, the output of the biomarkers from the machine learning models was intersected, followed by the utilization of multiple Cox regression to calculate the weight of each gene. The TEX-score formula is as follows:

$$\text{TEX}-\text{Score}=\text{X}1(\text{coefficient\hspace{0.17em}of\hspace{0.17em}multi}-\text{COX\hspace{0.17em}of\hspace{0.17em}gene}1)\ast \text{Y}1(\text{expression}-\text{level\hspace{0.17em}of\hspace{0.17em}gene}1)+\text{X}2\ast \text{Y}2+\text{X}3\ast \text{Y}3$$

Based on the median value of the TEX-score, patients in the OS TARGETs cohort and the meta-cohort (formed by combining data from GSE21257 and GSE16091 using the R package “Combat”) were stratified into high and low TEX-score groups. Subsequently, Kaplan-Meier survival analysis and receiver operator characteristic curves (ROC) between these two groups were conducted using the “survminer”, “survival”, “rms”, and “timeROC” R packages.

2.7 Clinical characteristic and nomogram establishment

Uni-Cox and multi-Cox regression analyses were employed to assess the correlation and independence of the TEX-score in conjunction with clinical parameters in the meta-cohort. In order to delineate disparities between patient subgroups, a nomogram was developed. This nomogram is capable of accurately forecasting an individual’s probability of encountering an event in a clinical setting, incorporating independent clinical prognostic factors like age, gender, metastasis, and TEX-score. The performance of the nomogram in prognostic prediction was subsequently evaluated using calibration and ROC curves, validating its predictive capability for prognosis (20).

2.8 Evaluation of immune-related characteristics for the TEX-signature

The immune cell components in each sample were computed using the Tumor Immune Estimation Resource (TIMER), single sample gene set enrichment analysis (ssGSEA), and Microenvironment Cell Populations-counter (MCP-counter) algorithm (21). Additionally, the “ESTIMATE” package was utilized to estimate both stromal and immune scores, enabling the quantification of the Tumor Microenvironment (TME) in malignant tumors (22). The cancer immune cycle, encompassing seven distinct steps (TIP, hrbmu.edu.cn), as well as various immune indicators calculated by the “easier” package, were used to gauge the immune capacity of the TME (23). Furthermore, an examination was conducted into the expression levels of co-stimulatory, co-inhibitory, and HLA molecules. Parameters including T cell-inflamed gene expression profile (GEP), cytotoxic activity (CYT), and IFN-γ were computed in accordance with previously established methodologies (21, 24, 25). TME signatures, independently developed by Kobayashi, were gathered and computed utilizing Gene Set Variation Analysis (GSVA) (26).

2.9 Prediction of immunotherapy

The immunotherapy data was sourced from several datasets, namely GSE91061 (melanoma), GSE126044 (lung adenocarcinoma), Nathon (melanoma), and Mel-ucla (2016, metastatic-melanoma), which were utilized to forecast the response to immunotherapy (27, 28). Additionally, GSE79671 (glioblastoma) and GSE61676 (non-small cell lung cancer) were employed to assess the effectiveness of antivascular drugs within high and low TEX-score groups (29, 30). The TEX-score was calculated independently for each dataset. Subsequently, drug screening was conducted for patients with differing TEX-scores using the “oncopredict” package.

2.10 Quantitative real-time polymerase chain reaction

Ethical approval was obtained from the Medical Ethics Committee for tissue specimens acquired from the Third Xiangya Hospital of Central South University (Approval No. fast-23816). These specimens were stored at a temperature of -80°C. A total of three pairs of samples were collected from OS patients who underwent tumor resection, including tumor tissue and paratumor tissue. Total RNA from tissues was isolated using the TRIzol reagent by Thermo Fisher Scientific, based in Waltham, MA, USA. The cDNA was synthesized from 2μg of RNA utilizing the RevertAid™ First Strand cDNA Synthesis Kit (Thermo Fisher Scientific). Quantitative real-time polymerase chain reaction (qRT-PCR) was performed using SYBR Green Master Mix (Q111-02, Vazyme). The quantification of relative gene expression levels was conducted using the 2^-△△CT method. The primer sequences are shown in Supplementary Table 1.

2.11 Western blot analysis and Immunofluorescence

Protein samples were collected using RIPA buffer (Beyotime, China) and the protein concentration was determined using a bicinchoninic acid (BCA) assay kit (Chinese Biotechnology Company). A total of 20 μg of protein was separated by 12% SDS-PAGE and transferred onto PVDF membranes (Bio-Rad). After blocking with 5% non-fat milk at room temperature for one hour, the membranes were incubated overnight with antibodies diluted in antibody solutions against RAD23A (Immunoway), SAC3D1 (Immunoway), GAPDH (Immunoway), and PSIP1 (Proteintech). Following washing, the membranes were then incubated with ananti-rabbit IgG solution at room temperature for one hour, followed by additional washing and visualization. For histological analysis, the specimens were fixed in 4% paraformaldehyde after removal, and the fixed tissues were embedded in paraffin for sectioning and subsequent staining. The antibodies used for immunofluorescence staining were as follows: anti-human CD8 (Abcam); anti-human RAD23A (Immunoway); anti-human SAC3D1(Immunoway), anti-human PSIP1(Proteintech), Alexa 546-conjugated anti-rabbit IgG (Invitrogen). Cell nuclei were counterstained with DAPI (Sigma Aldrich). Imaging was performed using a Zeiss Axio Observer Z1 LSM 710 BiG confocal microscope (Carl Zeiss), and fluorescence images were captured using Zen 2012 software (Carl Zeiss). Images were pseudocolored for overlay, cropping, resizing, and enhancing contrast and brightness using Photoshop and Illustrator (Adobe Systems) or ImageJ (NIH).

2.12 Statistical analysis

The statistical analyses were performed using R (version 4.2.2) and RStudio. A prognostic model for OS was developed employing Combined LASSO regression, Boruta, survival-SVM, Cox-boost, XG-boost, and GBM. For survival analysis and assessing the diagnostic value of the TEX-signature, Kaplan-Meier curves and the Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) were employed, respectively. In cases of normally distributed variables, significant quantitative differences between and among groups were determined using a two-tailed t-test or one-way ANOVA, as applicable. Conversely, for non-normally distributed variables, significant quantitative differences were assessed using a Wilcoxon test. A statistical P-value<0.05 was considered to be statistically significant.

3 Results

3.1 Single-cell analysis explored cell subtypes in OS

After controlling data quality and curating single-cell sequencing data from 6 OS patients, a total of 31,398 cells were screened and visualized through uniform manifold approximation and projection (Supplementary Figures S1 and 2A). The optimal number of cell populations was determined using the Seurat package, resulting in the differentiation of all cells into 13 distinct clusters (Figure 1A). Using the differential expression of genes between these 13 major clusters, combined with corrections for cell-specific cell markers for all subpopulations, an annotated classification of each cellular subpopulation within the osteosarcoma tumor microenvironment was performed. This included both immune (such as myeloid cells, NK/T cells, and B cells) and non-immune cells (such as osteoblastic OS cells, endothelial cells, OCs, and CAFs) (Figure 1B, C and Supplementary Figures 2B-H). After a comprehensive examination of the landscape and dynamics of immune cells, all groups of immune cells were re-grouped and annotated (Figures 1D, E). The NK/T cells underwent a similar process, while CD8⁺ T cells were specifically singled out for subsequent research (Figures 1F-H). Additionally, we observed a low expression level of CD4 within the T cell subpopulation (Figure 1H).

Figure 1

Figure 1 Different cell clustering in single cell sequencing data of osteosarcoma. (A) Identification 13 types of cells in single cell sequencing data. (B) Co-heatmap of marker genes for different cell types. (C) The 7 cell types were identified by marker genes. (D, E) Extraction, recombination and annotation of immune cells. (F-H) Screening of CD8+ T cell and NK/NKT cell cluster. OS: Osteoblastic-OS cell, EC: Endothelial cell, OC: Osteoclast, CAF: Fibroblasts, MC: Myeloid cell, NK/T: NK/NKT/T cell.

3.2 Analysis of CD8⁺ T cell site differentiation, cell population communication, and functional enrichment

The pseudo-time series analysis revealed the differentiation status among distinct clusters of CD8⁺ T cells and the rearrangement of cell types within the Tumor Microenvironment (TME) of OS (Figures 2A, B). Subsequently, we delved into the communication network among cell populations and found that CD8⁺ T cell cluster one exhibited a more active interaction status and weight compared to CD8⁺ T cell cluster two (Figure 2C). Classical inflammatory activation pathways such as TNF, OSM, and IFN-II signaling pathways displayed heightened activity in CD8⁺ T cell cluster one and Myeloid cells. Similarly, we investigated the signaling pattern and weight of cytokine families like IL-1, IL-2, IL-4, and IL-6 pathways, along with signaling pathways including TGF-beta, CCL, CD40, complement, and TRAIL in clusters (Figure 2D, E and Supplementary Figure 3A). In summary, CD8⁺ T cell cluster one demonstrated more pronounced advantages than cluster two across most inflammatory and immune activation pathways. Further visualization of the ligand interaction signal intensity revealed that the ligand interaction between CD8⁺ T cell cluster two and others was relatively attenuated in comparison to CD8⁺ T cell cluster one (Figure 2F). Furthermore, gene enrichment analysis was conducted between the two subtypes of CD8⁺ T cells to validate our hypothesis. The GSVA results of KEGG terms demonstrated a strong association of cluster one with cytokine, JAK-STAT, and T cell receptor signaling pathways (Figure 2G). Additionally, KEGG analysis was carried out by evaluating the up- and down-regulated differentially expressed genes in subgroup one of CD8⁺ T cells. It revealed that major pathways in cytotoxicity mediated by NK cells, necroptosis, TNF, and NOD signaling were positively enriched (Figure 2H). As for GO terms, a plethora of immune processes exhibited significant enrichment in cluster one, including inflammatory response, lymphocyte migration, proliferation and activation, as well as T cell differentiation and activation regulation (Figure 2I). Furthermore, the hallmark pathways of GSEA in cluster one indicated that molecules and pathways associated with immune function were highly activated (Figure 2J). In contrast, cluster two exhibited greater enrichment in metabolism-related pathways and lacked immune activation in GO, KEGG, and GSEA analyses (Figure 2H, Supplementary Figures 3B, C). Finally, the overall expression status of co-stimulatory and co-inhibitory molecules was compared between the two clusters (Figure 2K). Combining these findings with our previous results, we propose a process of functional exhaustion in the differentiation of CD8⁺ T cells between the two subsets.

Figure 2

Figure 2 Further analysis of CD8+ T cell. (A) Trajectory plots showing different clusters in CD8+ T cells. (B) Rearrangement of cell types in TIME (C) The number and weights/strength of interactions in the cell-cell communication network. (D) Weight of classical inflammatory activation pathways. (E) Signal pattern and weight of the cytokine family (F) Visualization of the signal intensity of ligand interaction. (G) GSVA analysis of CD8+ T cells. (H-J) The outstanding enrichment of GO, KEGG, and GSEA terms in cluster 1. (K) Box plots comparing co-stimulatory and co-inhibitory molecules between two clusters. CD8+T-C1: CD8+T cell cluster 1, CD8+T-C2: CD8+t cell cluster 2.

3.3 Exploring the genetic changes associated with exhaustion phenotype

We delved into the core gene-level alterations within the differentiation trajectory of CD8⁺ T cell subpopulations and found that there were noteworthy disparities in core genes between different clusters (Figures 3A, B). We utilized the hd-WGCNA algorithm to compute the gene expression profiles of the two CD8⁺ T cell subsets, and then categorized the core genes between these subsets into distinct gene modules to identify the core gene sets. Finally, we verified the correlation between genes and modules in the network. By setting β to 14, we achieved an R-squared value of 0.85, which established a scale-free network (Figure 3C). The genes were segregated into respective modules via hierarchical clustering, and a gene similarity heatmap was generated based on the topological overlap matrix (Figures 3D-F). The core genes were predominantly concentrated in the turquoise module, with a remarkably high correlation of 96% (Figure 3G). Further analysis revealed a strong correlation between genes within the block and the block (Figure 3H).

Figure 3

Figure 3 Exploring the genetic changes associated with exhaustion phenotype. (A) Volcano map of differentially expressed genes in CD8+ T cell clusters. (B) Heat map of differentially expressed genes in CD8+ T cell clusters. (C) Scale independence and average connectivity. (D) Network heatmap plot of all genes. (E) Cluster dendrogram. (F) Eigengene adjacency heatmap. (G) Heatmap of module–trait correlations. (H) Correlation between gene significance and module membership.

3.4 Machine learning to build TEX-signature

In our initial investigation into the biomarkers of T cell exhaustion, which have prognostic significance for patients with OS, we identified 668 genes that were commonly present in both cohorts (Figure 4A). Following univariate analysis of these exhausted core genes, feature selection was performed using six machine learning algorithms, including LASSO, XGboost, GBM, Boruta, CoxBoost, and survival-SVM (Figure 4B, Supplementary Figure 4A). The C-index values for all the algorithms exceeded 0.8, indicating the strong performance of each model (Figure 4C). Subsequently, we selected the intersection of the biomarkers obtained from the machine learning model to construct a refined model (Figure 4D). Three target genes, RAD23A, SAC3D1, and PSIP1, were screened along with their corresponding coefficients calculated through multivariate analysis (Figure 4E). Expression of target genes between the two cell clusters was also visualized (Figure 4F). Moreover, we defined the TEX-score as the sum of the product of the expression values and the correlation coefficients of these three genes separately. Comparisons of the prognostic status of patients with high TEX-scores against those with low TEX-scores revealed that in both the TARGETs cohort and the meta-cohort, patients with high TEX-scores exhibited worse clinical outcomes, while those with low TEX-scores demonstrated better outcomes (Figures 4G, H). Additionally, the area under the ROC curve demonstrated the excellent diagnostic efficacy and predictive ability of the model at 1, 2, 3, 5, and 10 years in both the TARGETs cohort and the meta-cohort (Figures 4I, J).

Figure 4

Figure 4 Construction and validation of the TEX-Signature. (A) Venn diagram of exhausted genes and core genes of the turquoise module. (B) machine learning for valuable models. (C) The C-index of various algorithms. (D) Venn plot showing the intersection of valuable TEX-genes based on six machine learning algorithms. (E) coefficient of three target genes calculated by multivariate analysis (F) Density plots of gene expression intensity for target genes. (G) Kaplan-Meier survival curve of OS between patients with a relative high score of TEX-Signature and a low score of TEX-Signature in the TARGETs cohort. (H) Kaplan-Meier survival curve of OS between patients with a high score of TEX-Signature and a low score of TEX-Signature in the meta-cohort. (I) Time-dependent ROC curve at 1, 2, 3, 5, and 10 years in the TARGETs cohort. (J) Time-dependent ROC curve at 1, 2, 3, 5, and 10 years in the meta-cohort.

3.5 Construction of a nomogram and clinical characteristic subgroup analysis

We developed a nomogram that incorporated age, gender, metastasis, and TEX-scores for clinical prediction (Figure 5A, Supplementary Figures 4B–E). The AUC values for 1-, 2-, 3-, 5-, and 10-year Overall Survival (OS) for the nomogram were 0.96, 0.89, 0.82, 0.79, and 0.79, respectively, indicating that our model exhibited strong and consistent predictive capability (Figure 5B). Furthermore, the calibration plots illustrated the level of agreement between the predicted OS and the actual OS (Figure 5C). To further underscore the predictive potential of the TEX-signature, we conducted subgroup analyses based on available clinical features in the TARGETs and GSE21257 databases. The signature demonstrated accurate and robust performance across these subgroups. According to Kaplan-Meier survival analysis, the low TEX-score group consistently exhibited a superior prognosis compared to the high TEX-score group within subgroups stratified by OS type, gender, age, or metastasis (Figures 5D–L and Supplementary Figures 4F–J). In addition, there is a tendency for the TEX-score to decrease with age, suggesting to some extent that they may be generalizable (Supplementary Figures 4K-L).

Figure 5

Figure 5 Construction of a nomogram and subgroup analysis. (A) Nomogram based on gender, age, metastasis and the TEX-score. (B) ROC curves for predicting 1-, 2-, 3-, 5- and 10-year survival in the TARGETs database. (C) The calibration plot for the probability of 1-, 2-, 3- and 5- year overall survival of OS patients. (D-H) Kaplan–Meier survival analysis for OS patients with diverse clinical characteristics of osteosarcoma type (D), gender (E, F), age (G, H) in TARGETs cohort. (I) Boxplot of TEX scores between metastatic and non-metastatic patients in the GSE21257 dataset (J-L) Kaplan-Meier survival analysis of metastatic and non-metastatic patients in the TARGETs (J) and the GSE21257 dataset (K-L).

3.6 Immune characteristics related to the TEX-signature

We investigated the association between TEX-signature and immune cell infiltration as well as immunomodulators in both the TARGETs cohort and the meta-cohort to evaluate the impact of TEX-signature in OS. Patients with high TEX-scores displayed a strong correlation with tumor purity, whereas the low TEX-score group exhibited a more favorable immune microenvironment and matrix score (Figure 6A and Supplementary Figure 5A). Moreover, we discovered that the low TEX-score group had positive associations with immune cells such as CD8⁺ T cells, macrophages, natural killer cells (NK cells), NK T cells, B cells, and central and effector memory T cells, all of which play important roles in positive immune regulation and immune-mediated killing utilizing multiple algorithms including ESTIMATE, TIMER, MCP-counter, and ssGSEA (Figure 6A, Supplementary Figures 5B–D). However, myeloid-derived suppressor cells (MDSC) were also more enriched in the low TEX-score group (Figure 6A and Supplementary Figures 5B-D). Additionally, the TEX-score exhibited negative correlations with most immune modulators, classified as antigen presentation, co-stimulatory, co-inhibitory, receptor, and others in the TARGETs cohort. The expression status of all immune checkpoint molecules was also depicted (Figures 6B, C). Furthermore, we confirmed the expression levels of co-stimulatory, co-inhibitory, and HLA molecules in the meta-cohort (Supplementary Figures 5E, F). Finally, we explored several immunotherapy indices in both the TARGETs cohort and meta-cohort. High levels of GEP, CYT, and IFN-γ were significantly associated with a low TEX-score, all of which are determinants of a potentially improved immunotherapy response (Figures 6D-F and Supplementary Figures 5G-I). Interestingly, there was no statistically significant difference in IFN in the TARGETs cohort, although our validation in the meta-cohort indicated that the results were still meaningful. Our results showed a clear intrinsic correlation between the immune microenvironment and TEX-scores, with the low TEX-score patients having an “immune-heat response phenotype” reflecting a better immunotherapeutic potential, whereas the high TEX-score group showed an “immune-poor state”.

Figure 6

Figure 6 Immune-related characteristics of TEX-signature in TARGETs and Meta-cohort. (A) Heatmap showing the correlation between TEX-score and immune infiltrating. (B) Heatmap showing the correlation between TEX-score and immune modulators. (C) Box plots comparing the expression status of all immune checkpoint molecules between low and high TEX-score groups. (D-F) Boxplot and scatter plot displaying the levels of GEP, CYT, and IFN-γ between low and high TEX-score groups. *p<0.05, **p<0.01, ***p<0.001.

3.7 Potential biological process related to TEX-signature

We examined the immunity cycle of cancer to elucidate the relationship between immune processes and TEX-score across the entire dataset, several steps of the immune cycle were found to be more activated in the low TEX-score group in our study. They included cancer antigen presentation, recruitment of T cells, CD8⁺ T cells, Th1 cells, NK cells, macrophages, B cells, infiltration of immune cells into tumors, and killing of cancer cells (Figure 7A). Moreover, we gathered various indices, including the T cell inflammatory microenvironment signature (T-cell-inflamed), which is based on the combined potential of IFN-γ and T-cell-associated inflammatory genes in predicting the response to PD-1 blockade, as well as the immunological characteristics of Roh (Roh-IS) associated with immune activation related to tumor rejection, and the immunological characteristics of Davoli (Davoli-IS), defined by the expression of cytotoxic CD8⁺ T cell, NK cell markers, and immuno-expanded label (Ayers-expIS), which is produced by genes highly associated with IFN-γ signature genes. All of these scores were highly significant in the low TEX-score group (Figure 7B). The immune resistance program (resF-down, resF-up, and resF) represents the efficacy of immune resistance in the tumor microenvironment, with patients in the high TEX-score group exhibiting stronger immune resistance (resF, resF-up), while lower levels of immune resistance (resF-down) were present in the low TEX-score group (Figure 7B). Furthermore, we examined the signatures developed by Kobayashi in the TARGETs cohort, where a low TEX-score was associated with higher levels of recognition of tumor cells, innate immunity, T cells, IFN-γ response, Tregs, and MDSCs, while proliferation levels were positively correlated with TEX-scores (Figure 7C). Additionally, we found that transcription factors associated with inflammation and tumor suppression, such as USF1, USF2, RFX5, TP53, ETS1, SPI1, GATA, and STAT1, were highly expressed in the low TEX-score group. Factors that play a bidirectional role in proliferation and immunity, including NF-κB, STAT5B, and STAT6, were also highly expressed in the low TEX-score group. Other major potential tumor growth factors, such as POU2F2, RUNX1, ERG, REL, and JUN, were also relatively increased in the low TEX-score group. Except for FOXO1, KLF4, and SMAD4, the highly expressed transcription factors promoted OS proliferation, metastasis, and drug resistance in the high TEX-score group, including the E2F family, MYC, TFDP1, ZEB1, TFAP2C, LEF1, FOSL1, TCF7L2, TWIST1, GLI2, and FOXO3 (Figure 7D).

Figure 7

Figure 7 Potential biological process of the TEX-signature. (A) Boxplot showing the differences in the cancer immunity cycle between low and high TEX-score groups. (B) The differences in immune-related indexes collected from the “easier” package in groups. (C) Heatmap showing the correlation between TEX-score and immune level developed by Kobayashi. (D) Heatmap showing the correlation between the TEX-score and transcription factors. *p<0.05, **p<0.01, ***p<0.001.

3.8 Function enrichment and metabolism of TEX-signature

We examined the similarities and differences between TEX-score subgroups at the level of specific biological functions and pathways. The GO analysis of biological processes primarily encompassed positive reactions of leukocytes, such as immune migration, adhesion, activation, and phagocytosis in the low TEX-score group. This specifically included the activation and differentiation of CD8⁺ T cells, B cells, and myeloid cell-mediated immunity (Figure 8A). The cellular components identified in the GO analysis were associated with membrane and filopodium components. Molecular functions included IgG binding, immunoglobulin binding, immune receptor activity, and serine-type peptidase activity (Figure 8B). The GO analysis results for the high TEX-score group were also notable. There was enrichment related to ion channels, both voltage-dependent and independent in biological processes. Terms like transporter complex and channel complex were enriched for cellular components. Bone morphogenesis and ossification were significantly enriched in terms of molecular functions (Supplementary Figure 6A). Additionally, we conducted a KEGG analysis that showed a significant enrichment in the low TEX-score group. We visually compared the enrichment status of the corresponding pathways in the two subgroups. Pathways including the Toll-like receptor, T cell receptor, NOD-like receptor, leukocyte trans-endothelial migration, Fc-γ receptor-mediated phagocytosis, cytokine-cytokine receptor interaction, chemokine signaling pathway, and B cell receptor, all of which were involved in immune processes, were significantly associated with the low TEX-score group (Figure 8C). Furthermore, numerous Hallmark signaling pathways of GSVA correlated with the low TEX-score, included complement, IL6-JAK-STAT3, and IL2-STAT5 signaling pathway, inflammatory response, IFN-α and IFN-γ response, and TNFα-NF-kB signaling pathway. As for the high TEX-score group, the results were consistent with what we obtained previously (Figure 8D).

Figure 8

Figure 8 Function enrichment and metabolism of the TEX-signature. (A) Biological process of GO analysis in the low TEX-score group. (B) Cell component and molecular function of GO analysis in the low TEX-score group. (C) KEGG analysis in the low TEX-score group. (D) GSVA analysis of low and high TEX-score groups.

3.9 Predictive efficacy of TEX-signature for therapy

To further investigate the potential value of the TEX-signature in therapy response, we proceeded to validate its efficacy in multiple published therapy datasets. The predictive capacity of the TEX-signature was well-evidenced by Disease Control Rate (DCR) in the context of immunotherapy. Patients with low TEX-scores in the GSE91061 cohort exhibited significantly improved DCR compared to patients with high TEX-scores, and the ROC curve confirmed the robustness of the TEX-score in predicting therapy response (Figure 9A). Similarly, patients with low TEX-scores in the GSE126044 dataset demonstrated a higher likelihood of responding positively to immunotherapy (Figure 9B). Patients with low TEX-scores in the Nanthon dataset experienced extended survival times and were more inclined to respond to immunotherapy (Figure 9C). Patients with low TEX-scores in the Mel-ucla dataset exhibited a superior DCR (Figure 9D). Turning to anti-angiogenic therapy, patients with low TEX-scores in the GSE79671 dataset were more prone to positive responses to anti-angiogenic therapy (Figure 9E). Likewise, in the GSE61676 dataset, patients with low TEX-scores demonstrated prolonged survival times and were more likely to respond favorably to anti-angiogenic therapy (Figure 9F). In addition, we verified the predictive value of the TEX-signature for chemotherapies. As shown, the TEX-score showed a significant correlation with major chemotherapeutic agents, including docetaxel-tanespimycin, regorafenib, sorafenib, topotecan, pazopanib, and paclitaxel. Patients with high TEX-scores appeared to be more likely to respond positively to chemotherapies (Figures 9G, H and Supplementary Figures 6B-E). Taken together, our research revealed that patients with low TEX-scores could potentially benefit more from certain treatment options.

Figure 9

Figure 9 Predictive value of TEX-signature in therapy response. (A-D) Distribution of immunotherapy responses and ROC curves in the GSE91061 cohort, GSE126044 cohort, Nanthon dataset and Mel-ucla dataset based on the TEX-signature (E, F) Distribution of anti-angiogenic therapy responses and the ROC curve between low and high TEX-score groups in the GSE79617 and GSE61676 dataset. (G, H) The sensitivity of docetaxel-tanespimycin and regorafenib in two groups.

3.10 The validation of TEX-related gene expression

To validate the expression patterns of TEX-related genes in osteosarcoma (OS) patients, we performed RT-PCR and Western Blot analyses on tumor tissues and adjacent non-tumor tissues from three patients. The results showed that, compared to adjacent non-tumor tissues, the expression of RAD23A, SAC3D1, and PSIP1 was significantly upregulated in OS tissues (Figures 10A, B). Additionally, we characterized the localization expression of the three target genes in CD8 T cells using immunofluorescence in OS tissue (Figure 10C). Therefore, we propose that dysregulated expression of these genes may lead to T-cell exhaustion and promote OS progression.

Figure 10

Figure 10 Validation of TEX-related gene expression. (A) QRT-PCR analysis and (B) Western blot analysis of RAD23A, PSIP1, and SAC3D1. (C) Immunofluorescence co-localization of CD8 and TEX-related genes was performed under 40x magnification. *P < 0.05, **P < 0.01, ***P < 0.001.

4 Discussion

Chemotherapy and surgical resection have long been the mainstay of treatment for OS (31). Unfortunately, there has been limited advancement in the treatment of OS over the last three decades, particularly in contrast to the notable progress made in developing novel therapies for other types of cancer (32). This stagnation in innovation has regrettably not translated into improved survival rates for patients dealing with OS. The need for further research and breakthroughs in treatment options for OS remains paramount.

Immunotherapy, an emerging therapeutic approach, has made significant strides in treating various cancer types. However, its impact on OS has been relatively limited (33). Effectively reshaping the immunosuppressive tumor microenvironment is crucial for the success of immunotherapy (34). Nevertheless, the intricate interplay of factors including the complex immune microenvironment, tumor heterogeneity, and individual variations in OS poses formidable challenges to harnessing the full potential of immune-based treatments (35). The core of immunotherapy is the activated T cells, particularly CD8⁺ T cells, whose functional state closely correlates with immune response effectiveness. However, CD8⁺ T cells are often altered or exhausted due to prolonged exposure to high levels of persistent antigen and inflammatory stimuli during tumor progression. These exhausted T cells lose their ability to eliminate tumor cells (11, 36). Immune checkpoint inhibitors such as anti-PD-1 antibodies and anti-CTLA-4 antibodies, regulatory cytokines, and metabolic reprogramming targeting the tumor microenvironment work to reverse the exhausted T cell states, restore their functionality, and reactivate immune responses. Previous studies have indicated that the number of tumor-infiltrating lymphocytes (TIL) is significantly higher in OS compared to other sarcomas (37), which suggests that immune checkpoint inhibitors may be able to leverage the abundance of TIL in OS, offering hope for immunotherapy in this context. However, the relationship between T cell exhaustion and OS remains inadequately understood.

CD8⁺ T cells, originating from CD34 hematopoietic stem cells located in the bone marrow, can be activated by endogenous antigenic peptides presented in MHC class I molecules, thereby exerting anti-tumor immunity (38). When the functionality of CD8⁺ T cells is compromised, the body’s anti-tumor immune capacity diminishes, elevating the risk of tumor growth and cancer metastasis (39). Through an exploration of the molecular and functional attributes of distinct CD8⁺ T cell subgroups in OS, we observed indications of functional exhaustion within the tumor immune microenvironment. In contrast to relatively exhausted CD8⁺ T cells, their more active counterparts demonstrated heightened engagement in cellular interactions, with most immune-related pathways exhibiting elevated activity. These pathways encompassed inflammatory activation pathways, TNF family members, complement C3, cytokine family, immune response, cell-cell adhesion, necroptosis, and T cell activation. The relatively exhausted cell subgroup exhibited heightened expression of markers like LAG-3, TOX, CTLA-4, aligning with prior studies elucidating mechanisms associated with CD8⁺ T cell exhaustion (40). This expression profile may potentially impact the immune response and prognosis of OS patients. These findings inspire us to further refine and scrutinize exhaustion models, seeking additional insights to advance immunotherapy for OS.

Following the application of six machine learning algorithms, we identified RAD23A, SAC3D1, and PSIP1 as genes associated with T cell exhaustion, forming the basis for an OS prognostic model. RAD23A, also known as RAD23 or HR23A, is involved in nucleotide excision repair and the regulation of intracellular protein degradation (41). Previous pan-carcinoma analyses have indicated a significant positive correlation of RAD23A in various cancers (42). It participate in processes such as nuclear translocation of AIF during cell death induction and enhances resistance to chemical agents by modulating autophagic response (43). RAD23A may mediate T cell exhaustion through diverse pathways and is recognized as an immune function biomarker, substantiating its inclusion in the prognostic model (44). SAC3D1, or SHD1, is implicated in centrosome duplication and mitotic progression, potentially mediating cell cycle regulation via centrosome amplification (45). SAC3D1 is involved in immune response, as well as its association with metabolism-related signaling pathways, positions it as a key player in T cell exhaustion and provides valuable insights for prognosis and immunotherapy effectiveness in various cancers (46, 47). PSIP1, also known as LEDGF/p75, participates in various biological processes and plays a significant role in lens epithelium differentiation into fiber cell terminals (48). The precise influence of these genes on the occurrence and progression of T cell exhaustion, particularly in relation to CD8⁺ T cell exhaustion in OS, deserves further exploration. Indeed, understanding the intricate interplay between the tumor immune microenvironment and T cell exhaustion is crucial for unraveling the complexities of cancer progression and devising effective therapeutic strategies (49, 50). The observations made in this study regarding immune cell infiltration and immune checkpoints between individuals with high and low TEX-scores shed light on potential avenues for enhancing immune efficacy. The heightened presence of CD8+ cells, macrophages, NK cells, NK T cells, B cells, and monocyte cells in the low TEX-score group signifies a more active immune response, which aligns with the notion of reduced T cell exhaustion. Moreover, the association of immune pathways with TEX-score subgroups provides valuable insights into the potential effectiveness of immunotherapy in OS. The activation of the WNT/β-catenin signaling pathway in the high TEX-score group is particularly noteworthy, as this pathway is known to play a crucial role in T cell differentiation and effector function (51). Previous studies have demonstrated that the activation of WNT/β-catenin signaling can suppress the effector functions of CD8+ T cells, further emphasizing its relevance in the context of T cell exhaustion (52). The expression of immune checkpoint molecules can originate from various cell types, including tumor cells, regulatory T cells (Treg), fibroblasts, or their extracellular vesicles. In the low-TEX group, co-stimulatory molecules such as CD40/CD40LG and CD96 were relatively upregulated, accompanied by elevated levels of TNF-α, GZMA, and IFN-α/IFN-γ. Meanwhile, we observed a relative upregulation of some inhibitory immune checkpoints, this may be attributed to heightened antigen presentation stimulation due to high HLA expression and a tumor’s self-protective effect induced by sustained inflammatory responses. The upregulation of CD274, IDO1, etc. on the tumor surface by T cell activation and IFN-α/IFN-γ stimulation has been demonstrated (53, 54). Furthermore, the Meta-dataset’s high-TEX group showed increased expression of TOX and VCTN1, and the relationship between immune checkpoint regulation and immune infiltration was more intricate than we had first thought. Patients with elevated immune checkpoint levels may also exhibit higher levels of immune activation, and this group of patients may experience better clinical benefit from combination immunotherapy (55). This highlights the potential for interventions aimed at reversing the state of immune exhaustion in the tumor microenvironment, a development that could have far-reaching implications for cancer therapy. Reassuringly, recent clinical successes in reversing T cell exhaustion underscore the promising potential of such approaches (56).

As we know, T cell exhaustion is a prolonged and persistent process characterized by the upregulation of various immune inhibitory factors and impaired functionality, such as compromised release of IFN-γ and granzymes, within the tumor immune microenvironment (TIME) under inflammatory stimuli. Despite being the mainstay of immunotherapy, classical immune inhibitors like Anti-PD1 and Anti-CTLA4, represented by immune checkpoint blockade (ICB), unfortunately, fail to provide long-term benefits for a significant proportion of patients (57). The restoration of exhausted T cell functions is often limited, and they can rapidly revert to their pre-treatment state. Current research has identified CD8+ and Th1-type T cell markers, including IFN-γ, PRF1, and TAP1, to be significantly correlated with patients’ responses to immunotherapy (58). Additionally, scholars have found that early PD-1 blockade combined with CAR-T therapy can achieve better prognosis improvement (59). Therefore, for patients with relatively low tumor heterogeneity, high immune infiltration, and limited exhaustion, immunotherapy may attain better long-term efficacy (60). In our research cohort, besides the significant correlation of important indicators such as GEP, IFN-γ, and CYT with low exhaustion levels, the consistent performance of scores like Roh-is, Davoli-is, and RIR further supports our hypothesis, affirming the favorable prognosis of low TEX and providing support for our hypothesis. The validation of this model in multiple therapy datasets across different tumor types further strengthens its predictive efficiency. The findings regarding the sensitivity of patients to anti-angiogenic drugs and conventional chemotherapy drugs for OS highlight the potential clinical utility of the TEX-signature in guiding treatment decisions. However, it’s important to acknowledge the need for further validation and clinical implementation. This study sets a promising foundation for future research and potential advancements in the treatment of OS.

It should be mentioned that there are a few of restrictions. Firstly, due to tumor heterogeneity and limited sample size, the study’s findings are based on single-cell sequencing data from a relatively small sample number, which may not fully capture the heterogeneity present in osteosarcoma. Further validation in larger cohorts would provide more robust and generalizable results. Second, while the study identifies core genes in the TEX-signature, further molecular experiments are necessary to elucidate the functional roles of these genes and understand the underlying molecular mechanisms of CD8⁺ T cell exhaustion. Finally, the study did not specifically address the prediction of metastasis in osteosarcoma. It’s important to acknowledge that the model’s performance in this regard remains uncertain, because we were unable to get reliable metastasis-related data.

In conclusion, our study aims to analyze the immune microenvironment and tumor heterogeneity in OS using single-cell sequencing data, identifying distinct differentiation trajectories of CD8⁺ T cells in different individuals, and conducting a thorough evaluation of CD8⁺ T cells, which holds promise in shedding light on new avenues for OS immunotherapy.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The animal study was approved by the Medical Ethics Committee of the Third Xiangya Hospital of Central South University. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

QF: Conceptualization, Funding acquisition, Investigation, Writing – original draft, Data curation, Validation, Visualization. YW: Funding acquisition, Software, Conceptualization, Data curation, Investigation, Project administration, Visualization, Writing – original draft. JC: Data curation, Formal analysis, Methodology, Supervision, Writing – review & editing. BP:Methodology, Funding acquisition, Investigation, Validation, Writing – original draft. XZ: Validation, Funding acquisition, Writing – review & editing. RL: Investigation, Writing – original draft, Conceptualization, Data curation, Formal analysis. YD: Funding acquisition, Resources, Software, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by Natural Science Foundation of Hunan Province (2023JJ30855), Changsha Natural Science Foundation (kq2208364) and the Fundamental Research Funds for the Central Universities of Central South University (CX20220313).

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.2024.1362970/full#supplementary-material

Supplementary Figure 1 | (A) Density plot of transcript counts detected in each sample as a proportion of cell number. (B) Total number of genes detected in each sample. (C) Distribution of the percentage of mitochondrial genes in each sample. (D) Correlation analysis of total transcript number with mitochondrial genes, total gene number and PCA clustering analysis of all cells.

Supplementary Figure 2 | (A) Two-dimensional visualization of cell distribution across samples. (B-H) The “umap” visualization of cell subpopulation-specific marker expression levels.

Supplementary Figure 3 | (A) Communication networks among cell populations of multiple signaling pathways in cluster 2. (B) The enrichment of GO and GSEA terms in cluster 2.

Supplementary Figure 4 | (A) Risk coefficients for potential TEX-related genes calculated by univariate analysis. (B-E) The risk correlation analysis for each clinical characteristic and TEX-score in (B-C) The TARGETs and (D-E) the GSE21257 dataset. (F) Comparison of TEX scores of metastatic patients in the TARGETs cohort. (G-J) Kaplan–Meier survival analysis for OS patients with diverse clinical characteristics of age (G, H), and gender (I, J) in the GSE21257 dataset. (K-L) Distribution of TEX scores for all patients by age.

Supplementary Figure 5 | (A-C) Immune characteristics related to the TEX-signature. (D) Heatmap showing the correlation between TEX-score and immune infiltrating. (E-F) The expression status of immune checkpoint and antigen-presenting molecules between low and high TEX-score groups. (G-I) The levels of CYT, IFN-γ and GEP between low and high TEX-score groups.

Supplementary Figure 6 | (A) Function enrichment and metabolism of the high TEX-score. (B-E) The sensitivity of Sorafenib, Topotecan, Pazopanib, and Paclitaxel between low and high TEX-score.

References

1. Matsuoka K, Bakiri L, Wolff LI, Linder M, Mikels-Vigdal A, Patiño-García A, et al. Wnt signaling and Loxl2 promote aggressive osteosarcoma. Cell Res. (2020) 30:885–901. doi: 10.1038/s41422-020-0370-1

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Isakoff MS, Bielack SS, Meltzer P, Gorlick R. Osteosarcoma: current treatment and a collaborative pathway to success. J Clin Oncol. (2015) 33:3029–35. doi: 10.1200/JCO.2014.59.4895

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Kciuk M, Yahya EB, Mohamed Ibrahim Mohamed M, Rashid S, Iqbal MO, Kontek R, et al. Recent advances in molecular mechanisms of cancer immunotherapy. Cancers (Basel). (2023) 15. doi: 10.3390/cancers15102721

CrossRef Full Text | Google Scholar

4. Wedekind MF, Wagner LM, Cripe TP. Immunotherapy for osteosarcoma: Where do we go from here? Pediatr Blood Cancer. (2018) 65:e27227. doi: 10.1002/pbc.27227

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Bielack SS, Smeland S, Whelan JS, Marina N, Jovic G, Hook JM, et al. Methotrexate, doxorubicin, and cisplatin (MAP) plus maintenance pegylated interferon alfa-2b versus MAP alone in patients with resectable high-grade osteosarcoma and good histologic response to preoperative MAP: first results of the EURAMOS-1 good response randomized controlled trial. J Clin Oncol. (2015) 33:2279–87. doi: 10.1200/JCO.2014.60.0734

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Tawbi HA, Burgess M, Bolejack V, Van Tine BA, Schuetze SM, Hu J, et al. Pembrolizumab in advanced soft-tissue sarcoma and bone sarcoma (SARC028): a multicentre, two-cohort, single-arm, open-label, phase 2 trial. Lancet Oncol. (2017) 18:1493–501. doi: 10.1016/S1470-2045(17)30624-1

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Somaiah N, Conley AP, Parra ER, Lin H, Amini B, Solis Soto L, et al. Durvalumab plus tremelimumab in advanced or metastatic soft tissue and bone sarcomas: a single-centre phase 2 trial. Lancet Oncol. (2022) 23:1156–66. doi: 10.1016/S1470-2045(22)00392-8

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Kawano M, Itonaga I, Iwasaki T, Tsumura H. Enhancement of antitumor immunity by combining anti-cytotoxic T lymphocyte antigen-4 antibodies and cryotreated tumor lysate-pulsed dendritic cells in murine osteosarcoma. Oncol Rep. (2013) 29:1001–6. doi: 10.3892/or.2013.2224

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Sun L, Su Y, Jiao A, Wang X, Zhang B. T cells in health and disease. Signal Transduct Target Ther. (2023) 8:235. doi: 10.1038/s41392-023-01471-y

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Russell JH, Ley TJ. Lymphocyte-mediated cytotoxicity. Annu Rev Immunol. (2002) 20:323–70. doi: 10.1146/annurev.immunol.20.100201.131730

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Wherry EJ. T cell exhaustion. Nat Immunol. (2011) 12:492–9. doi: 10.1038/ni.2035

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Nakamoto N, Cho H, Shaked A, Olthoff K, Valiga ME, Kaminski M, et al. Synergistic reversal of intrahepatic HCV-specific CD8 T cell exhaustion by combined PD-1/CTLA-4 blockade. PloS Pathog. (2009) 5:e1000313. doi: 10.1371/journal.ppat.1000313

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Bi J, Jin X, Zheng C, Huang C, Zhong C, Zheng X, et al. Checkpoint TIPE2 limits the helper functions of NK cells in supporting antitumor CD8+ T cells. Adv Sci (Weinh). (2023) 10:e2207499. doi: 10.1002/advs.202207499

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Van Calster B, Van Smeden M, De Cock B, Steyerberg EW. Regression shrinkage methods for clinical prediction models do not guarantee improved performance: Simulation study. Stat Methods Med Res. (2020) 29:3166–78. doi: 10.1177/0962280220921415

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Kursa MB, Jankowski A, Rudnicki WR. Boruta - A system for feature selection. Fundamenta Informaticae. (2010) 101:271–86. doi: 10.3233/FI-2010-288

CrossRef Full Text | Google Scholar

16. Van Belle V, Van Calster B, Van Huffel S, Suykens JAK, Lisboa P. Explaining support vector machines: A color based nomogram. PloS One. (2016) 11:e0164568. doi: 10.1371/journal.pone.0164568

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Tutz G, Binder H. Generalized additive modeling with implicit variable selection by likelihood-based boosting. Biometrics. (2006) 62:961–71. doi: 10.1111/j.1541-0420.2006.00578.x

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Ester M, Kriegel HP, Xu X. XGBoost: A scalable tree boosting system. In Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. Geographical Anal. (2022), 785. doi: 10.1111/gean.12315

CrossRef Full Text | Google Scholar

19. Cham H, West SG. Propensity score analysis with missing data. Psychol Methods. (2016) 21:427–45. doi: 10.1037/met0000076

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wang W, Xiang M, Liu H, Chu X, Sun Z, Feng L, et al. A prognostic risk model based on DNA methylation levels of genes and lncRNAs in lung squamous cell carcinoma. PeerJ. (2022) 10:e13057. doi: 10.7717/peerj.13057

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. (2017) 18:248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Yoshihara K, Shahmoradgoli M, Martínez 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

23. Chen DS, Mellman I. Oncology meets immunology: the cancer-immunity cycle. Immunity. (2013) 39. doi: 10.1016/j.immuni.2013.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Roh W, Chen P-L, Reuben A, Spencer CN, Prieto PA, Miller JP, et al. Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci Transl Med. (2017) 9. doi: 10.1126/scitranslmed.aah3560

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. (2017) 127:2930–40. doi: 10.1172/JCI91190

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kobayashi Y, Kushihara Y, Saito N, Yamaguchi S, Kakimi K. A novel scoring method based on RNA-Seq immunograms describing individual cancer-immunity interactions. Cancer Sci. (2020) 111:4031–40. doi: 10.1111/cas.14621

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Yang L, Wei S, Zhang J, Hu Q, Hu W, Cao M, et al. Construction of a predictive model for immunotherapy efficacy in lung squamous cell carcinoma based on the degree of tumor-infiltrating immune cells and molecular typing. J Transl Med. (2022) 20:364. doi: 10.1186/s12967-022-03565-7

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Nathanson T, Ahuja A, Rubinsteyn A, Aksoy BA, Hellmann MD, Miao D, et al. Somatic mutations and neoepitope homology in melanomas treated with CTLA-4 blockade. Cancer Immunol Res. (2017) 5:84–91. doi: 10.1158/2326-6066.CIR-16-0019

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Urup T, Staunstrup LM, Michaelsen SR, Vitting-Seerup K, Bennedbæk M, Toft A, et al. Transcriptional changes induced by bevacizumab combination therapy in responding and non-responding recurrent glioblastoma patients. BMC Cancer. (2017) 17:278. doi: 10.1186/s12885-017-3251-3

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Baty F, Joerger M, Früh M, Klingbiel D, Zappa F, Brutsche M. 24h-gene variation effect of combined bevacizumab/erlotinib in advanced non-squamous non-small cell lung cancer using exon array blood profiling. J Transl Med. (2017) 15:66. doi: 10.1186/s12967-017-1174-z

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Strauss SJ, Frezza AM, Abecassis N, Bajpai J, Bauer S, Biagini R. Bone sarcomas: ESMO-EURACAN-GENTURIS-ERN PaedCan Clinical Practice Guideline for diagnosis, treatment and follow-up. Ann Oncol. (2021) 32:1520–36. doi: 10.1016/j.annonc.2021.08.1995

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Cole S, Gianferante DM, Zhu B, Mirabello L. Osteosarcoma: A Surveillance, Epidemiology, and End Results program-based analysis from 1975 to 2017. Cancer. (2022) 128:2107–18. doi: 10.1002/cncr.34163

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Evdokimova V, Gassmann H, Radvanyi L, Burdach SEG. Current state of immunotherapy and mechanisms of immune evasion in ewing sarcoma and osteosarcoma. Cancers (Basel). (2022) 15. doi: 10.3390/cancers15010272

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Bedognetti D, Ceccarelli M, Galluzzi L, Lu R, Palucka K, Samayoa J, et al. Toward a comprehensive view of cancer immune responsiveness: a synopsis from the SITC workshop. J Immunother Cancer. (2019) 7:131. doi: 10.1186/s40425-019-0602-4

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Zhang N, Bevan MJ. CD8(+) T cells: foot soldiers of the immune system. Immunity. (2011) 35:161–8. doi: 10.1016/j.immuni.2011.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Franco F, Jaccard A, Romero P, Yu YR, Ho PC. Metabolic and epigenetic regulation of T-cell exhaustion. Nat Metab. (2020) 2:1001–12. doi: 10.1038/s42255-020-00280-9

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Van Erp AEM, Versleijen-Jonkers YMH, Hillebrandt-Roeffen MHS, Houdt L, Gorris Mark A J, Van Dam LS. Expression and clinical association of programmed cell death-1, programmed death-ligand-1 and CD8+ lymphocytes in primary sarcomas is subtype dependent. Oncotarget. (2017) 8:71371–84. doi: 10.18632/oncotarget.v8i41

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Gascoigne NRJ. Do T cells need endogenous peptides for activation? Nat Rev Immunol. (2008) 8:895–900. doi: 10.1038/nri2431

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Gonzalez SM, Taborda NA, Rugeles MT. Role of different subpopulations of CD8+ T cells during HIV exposure and infection. Front Immunol. (2017) 8:936. doi: 10.3389/fimmu.2017.00936

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Scott AC, Dündar F, Zumbo P, Chandran SS, Klebanoff CA, Shakiba M. TOX is a critical regulator of tumour-specific T cell differentiation. Nature. (2019) 571:270–4. doi: 10.1038/s41586-019-1324-y

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Schauber C, Chen L, Tongaonkar P, Madura K. Rad23 links DNA repair to the ubiquitin/proteasome pathway. Nature. (1998) 391:715–8. doi: 10.1038/35661

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Zhao S, Wu Y, Wei Y, Wei Y, Xu X, Zheng J. Identification of biomarkers associated with CD8+ T cells in coronary artery disease and their pan-cancer analysis. Front Immunol. (2022) 13:876616. doi: 10.3389/fimmu.2022.876616

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Tan X, Wang H-C, Liang R-Y, Chuang SM. Human Rad23A plays a regulatory role in autophagy. Biochem Biophys Res Commun. (2016) 478:1772–9. doi: 10.1016/j.bbrc.2016.09.025

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Park S, So M-K, Huh J. Methylation of DNA repair genes as a prognostic biomarker in AML of a TCGA-LAML cohort. Clin Lab. (2022) 68. doi: 10.7754/Clin.Lab.2021.211025

CrossRef Full Text | Google Scholar

45. Liu Y, He M, Ke X, Chen Y, Zhu J, Tan Z, et al. Centrosome amplification-related signature correlated with immune microenvironment and treatment response predicts prognosis and improves diagnosis of hepatocellular carcinoma by integrating machine learning and single-cell analyses. Hepatol Int. (2023). doi: 10.1007/s12072-023-10538-5

CrossRef Full Text | Google Scholar

46. Nakajima H, Tamura T, Ito M, Shibata F, Kuroda K, Fukuchi Y, et al. SHD1 is a novel cytokine-inducible, negative feedback regulator of STAT5-dependent transcription. Blood. (2009) 113:1027–36. doi: 10.1182/blood-2008-01-133405

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Lu W, Huang J, Shen Q, Li J. Identification of diagnostic biomarkers for idiopathic pulmonary hypertension with metabolic syndrome by bioinformatics and machine learning. Sci Rep. (2023) 13:615. doi: 10.1038/s41598-023-27435-4

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Singh DP, Ohguro N, Kikuchi T, Sueno T, Reddy VN, Yuge K, et al. Lens epithelium-derived growth factor: effects on growth and survival of lens epithelial cells, keratinocytes, and fibroblasts. Biochem Biophys Res Commun. (2000) 267:373–81. doi: 10.1006/bbrc.1999.1979

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Van Der Heide V, Humblin E, Vaidya A, Kamphorst AO. Advancing beyond the twists and turns of T cell exhaustion in cancer. Sci Transl Med. (2022) 14:eabo4997. doi: 10.1126/scitranslmed.abo4997

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Babar Q, Saeed A, Tabish TA, Sarwar M, Thorat ND. Targeting the tumor microenvironment: Potential strategy for cancer therapeutics. Biochim Biophys Acta Mol Basis Dis. (2023) 1869:166746. doi: 10.1016/j.bbadis.2023.166746

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Gattinoni L, Ji Y, Restifo NP. Wnt/beta-catenin signaling in T-cell immunity and cancer immunotherapy. Clin Cancer Res. (2010) 16:4695–701. doi: 10.1158/1078-0432.CCR-10-0356

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Li X, Xiang Y, Li F, Yin C, Li B, Ke X. WNT/β-catenin signaling pathway regulating T cell-inflammation in the tumor microenvironment. Front Immunol. (2019) 10:2293. doi: 10.3389/fimmu.2019.02293

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Dong H, Strome SE, Salomao DR, Tamura H, Hirano F, Flies DB, et al. Tumor-associated B7-H1 promotes T-cell apoptosis: a potential mechanism of immune evasion. Nat Med. (2002) 8:793–800. doi: 10.1038/nm730

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Mclane LM, Abdel-Hakeem MS, Wherry EJ. CD8 T cell exhaustion during chronic viral infection and cancer. Annu Rev Immunol. (2019) 37:457–95. doi: 10.1146/annurev-immunol-041015-055318

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Liu L, Liu Z, Gao J, Liu X, Weng S, Guo C, et al. CD8+ T cell trajectory subtypes decode tumor heterogeneity and provide treatment recommendations for hepatocellular carcinoma. Front Immunol. (2022) 13:964190. doi: 10.3389/fimmu.2022.964190

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Dall'olio FG, Marabelle A, Caramella C, Garcia C, Aldea M, Chapu N, et al. Tumour burden and efficacy of immune-checkpoint inhibitors. Nat Rev Clin Oncol. (2022) 19:75–90. doi: 10.1038/s41571-021-00564-3

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Davoli T, Uno H, Wooten EC, Elledge SJ. Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science. (2017) 355. doi: 10.1126/science.aaf8399

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Papatriantafyllou M. Regulatory T cells: distilling regulatory T cell inducers. Nat Rev Immunol. (2013) 13:546. doi: 10.1038/nri3506

CrossRef Full Text | Google Scholar

59. Sznol M, Chen L. Antagonist antibodies to PD-1 and B7-H1 (PD-L1) in the treatment of advanced human cancer. Clin Cancer Res. (2013) 19:1021–34. doi: 10.1158/1078-0432.CCR-12-2063

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Vitale I, Shema E, Loi S, Galluzzi S. Intratumoral heterogeneity in cancer progression and response to immunotherapy. Nat Med. (2021) 27:212–24. doi: 10.1038/s41591-021-01233-9

PubMed Abstract | CrossRef Full Text | Google Scholar