Development of an immune-related signature for predicting survival outcome and immunotherapy response in osteosarcoma

Construction of immune-related risk signature in OS

Considering that the TARGET cohort is a special database for OS and contains the largest samples size. We therefore selected it as the discovery cohort to explore the immune-relevant signature that is linked with the survival risk of OS patients. The univariate Cox regression analysis was conducted in all of 1793 immune genes for the TARGET cohort, and a total of 328 prognostic genes were determined (all P < 0.05, Supplementary Table 1). The 328 immune prognostic genes were subsequently subjected to the LASSO regression model with 10-fold cross-validation setting to obtain a better regression model. The LASSO coefficient profile was generated against the log(k) sequence, and the minimize k method led to 21 optimal coefficients (Figure 1A and 1B). Finally, a developed model with 21 immune-related genes reached the optimal regression efficiency to speculate the prognostic ability.

Figure 1. Development of the immune-related risk signature. (A) Regression coefficient profiles of the prognostic immune genes generated by the LASSO model in the TARGET cohort. (B) Partial likelihood deviance representation of distinct gene combination in LASSO regression. The red dots indicated the partial likelihood of deviance values, the gray lines indicated the standard error (SE), the two vertical lines on the left and right respectively indicated optimal values by minimum criteria and 1-SE criteria. (C) Heatmap of the 21 immune-related genes within the identified immune signature and the risk score curve based on the LASSO coefficients. OS patients were stratified into high-risk and low-risk subgroups with the median risk score.

The determined immune-relevant genes involved in antigen processing and presentation (PSMC4), cytokines (STC2, TNFRSF11B, CORT, GAL, CMTM1, IL7, GCG), cytokine receptors (SSTR1, TNFRSF21, FGFRL1, SDC3, IL13RA2), chemokines (SEMA3E), antimicrobials (IL22, TMPRSS6, PPARG), BCR signaling pathway (CD79A, IGLV1-51, IGKVID-33), and TCR signaling pathway (TRAV9-2). Moreover, we established an immune risk signature to assess the risk score of each patient by combing the 21 genes mRNA expression levels and their corresponding coefficients weighted by the LASSO method (Figure 1C, Supplementary Table 2). A heatmap of the determined 21 genes expression profile and the scatterplot of overall survival with corresponding risk scores for OS were shown in Figure 1C.

The prognostic ability of the identified immune signature

To explore the immune risk signature contribute to OS prognosis evaluation, we partitioned the samples of the TARGET cohort into a low-risk subgroup (n = 42) and a high-risk subgroup (n = 43) by applying the median risk score as the cut-off value. We observed that OS patients with the low-risk scores were significantly connected with a preferable prognosis as compared with those in the high-risk group (Log-rank test P < 0.001; Figure 2A). This connection remained still significant in the multivariate Cox model with age, sex, and tumor site taken into account (HR: 0.02, 95% CI: 0.00–0.14, P < 0.001; Figure 2B).

Figure 2. The constructed immune risk signature was linked with OS prognosis. Kaplan-Meier survival curves according to distinct immune signature subgroups in (A) the TARGET discovery cohort, (C) GSE21257 cohort, and (E) TCGA sarcoma cohort. (B, D, F) Forest plot representation of the multivariate Cox regression models delineated the connection between immune risk signature and prognosis in the three datasets. Clinical confounding factors were taken into consideration.

To confirm that the 21-gene-based immune signature classifier harbored a similar prognostic ability in distinct populations, we further validate this connection in the GSE21257 cohort with 53 OS samples. In the corroboration dataset, we also observed that patients of low-risk group exhibited a favorable survival outcome than those in the high-risk group (Log-rank test P = 0.013; Figure 2C). Multivariate Cox regression model further demonstrated that the immune-related signature could serve as an independent indicator for OS prognosis even adjusting for clinical confounding factors (HR: 0.32, 95% CI: 0.13–0.81, P = 0.023; Figure 2D).

We also curated 262 sarcoma samples from the TCGA cohort to further corroborate the prognostic value of the identified immune signature. Consistently, the significantly better survival outcome was observed in patients with low-risk scores under univariate Kaplan-Meier analysis (Log-rank test P = 0.007; Figure 2E) and multivariate Cox regression model (HR: 0.61, 95% CI: 0.40–0.91, P = 0.007; Figure 2F). ROC curves of the identified immune signature risk scores of the 3 cohorts used in this study (TARGET cohort area under curve [AUC]: 0.951; GSE21257 cohort AUC: 0.623; TCGA sarcoma cohort: 0.701) were illustrated in Supplementary Figure 1.

The identified immune signature predictive of lymphocyte infiltration and immune response

Since the risk signature was identified based on the immune-related genes, we hypothesized that its existence may be associated with lymphocyte infiltration, signaling pathways enrichment, and immune response. Therefore, we achieved a multiple-group boxplot by using ssGSEA to infer and visualize the relative infiltration levels of 31 immune cell subpopulations in low-risk versus high-risk subgroups from the TARGET cohort (Figure 3A). Anti-tumor lymphocyte subtypes, like central memory CD4⁺ and CD8⁺ T cells, effector memory CD8⁺ T cells, cytotoxic cells, and monocytes were markedly enriched in patients with the low-risk scores (all P < 0.05). However, the regulatory T cells (Treg), which belong to pro-tumor leukocytes, were exhibited decreased infiltration in the low-risk subgroup (P < 0.01). Furthermore, GSEA analysis against the OS gene expression data and Hallmark reference datasets demonstrated the biological signaling pathways associated with this immune signature. Genes involved in interferon γ response, interferon α response, and inflammatory response signaling pathways were remarkably enriched in patients from the low-risk group (Figure 3B–3D).

Figure 3. The risk signature predictive of immune infiltration and immunotherapy response indicators. (A) A total of 31 lymphocytes infiltration in low-risk versus high-risk subgroups was evaluated by the Angelova et al. algorithm. Pathways involved in immune response, such as (B) interferon γ response, (C) interferon α response, and (D) inflammatory response were markedly enriched in OS patients with low-risk signature scores. The elevated enrichment of immunotherapy efficacy predictors of (E) T cell-inflamed signature, (F) IFNγ signature, and (G) cytolytic activity signature were also connected with low-risk scores.

The three well-known signatures, including IFNγ response, T cell-inflamed GEP, and cytolytic activity signatures were identified to be associated with immunogenicity and immunotherapy response. Herein, we utilized these three signatures to investigate the association between the novel determined signature and immune response. In the TARGET cohort, the significantly elevated enrichment of T cell-inflamed signature (Figure 3E), IFNγ signature (Figure 3F), and cytolytic activity signature (Figure 3G) was noticed in the low-risk subgroup (Wilcoxon rank-sum test all P < 0.001). These findings suggested that OS patients with the reduced enrichment of this immune signature may be more responsive to immunotherapy.

The immune signature for the prediction roles of immunotherapy efficacy

Immune checkpoint inhibitors (ICI) treatments represented by anti-PD-1/PD-L1 or anti-CTLA-4 agents have dramatically improved the survival outcome of advanced cancer patients. We therefore retrospectively acquired the transcriptomic expression profiles and clinicopathologic characteristics from an anti-PD-L1-treated cohort (Imvigor210) of advanced urothelial cancer (UC), so as to evaluate the link between the identified immune risk signature and immunotherapy efficacy. In this UC ICI cohort, patients with the low-risk scores harbored a significantly prolonged survival outcome as compared with those in the high-risk group (Log-rank test P = 0.015; Figure 4A). The link between low-risk signature scores and better ICI prognosis still existed in the multivariate Cox model with multiple clinical factors adjusted (HR: 0.66, 95% CI: 0.51–0.87, P = 0.003; Figure 4B). The significantly elevated response rate of ICI treatment was also detected in patients of low-risk subgroup (30.2% vs. 15.4%, Fisher exact test, P = 0.004, Figure 4C; Kruskal-Wallis H test, P < 0.001, Figure 4D). Further analysis demonstrated that the markedly enhanced tumor mutation burden (TMB) and neoantigen burden (NB) were both observed in OS patients with the low-risk scores, which are closely predictive of immunotherapy outcome (Wilcoxon rank-sum test both P < 0.001; Figure 4E and 4F).

Figure 4. The roles of the immune signature in evaluating ICI treatment efficacy for advanced urothelial cancer (UC). (A) Kaplan-Meier survival plots of low-risk versus high-risk subgroups in the anti-PD-L1-treated UC cohort. (B) Forest plot of the multivariate Cox regression models with sex, ECOG score, smoking status, and immune phenotype taken into account to exhibit the association of the identified immune signature with ICI survival. (C) The proportion of therapeutic advantages to anti-PD-L1 therapy in low versus high-risk OS subgroups. (D) Distribution of immune signature risk scores in patients with distinct ICI treatment effect. (E) Tumor mutation burden and (F) neoantigen burden distribution in low-risk versus high-risk OS subgroups with the genomic data from UC immunotherapy cohort.

Another immunotherapeutic cohort of advanced/ metastatic melanoma patients (treated with anti-PD-1/ PD-L1 or anti-CTLA-4 agents) was also curated for validating the ICI predictive roles of the immune risk signature. In this melanoma population, we observed that patients in the low-risk group exhibited both significantly extended overall survival (OS) and progression-free survival (PFS) outcomes (Log-rank test P = 0.019 and 0.021, respectively; Figure 5A and 5B). Under multivariate-adjusted analysis, these connections remained still significant even controlling for sex, stage, and therapy type (OS: HR, 0.51, 95% CI, 0.30–0.87, P = 0.016, Figure 5C; PFS: HR, 0.59, 95% CI, 0.38–0.93, P = 0.028, Figure 5D). Conformably, the higher proportion of ICI-responsive melanoma patients was also noticed in the low-risk subgroup (49.2% vs. 30.0%, Fisher exact test P = 0.039; Figure 5E). The above results indicated that the immunotherapy efficacy may be evaluated by using the newly identified immune-related risk signature.

Figure 5. The roles of the immune signature in evaluating ICI treatment efficacy for melanoma. (A) Overall survival (OS) rate and (B) progression-free survival (PFS) rate differences of low- and high-risk groups were illustrated by Kaplan-Meier survival curves. Multivariate Cox regression models were conducted with sex, stage, and treatment type taken into account to elucidate the connections of the immune signature with (C) OS and (D) PFS outcomes. (E) The proportion of melanoma ICI responders in patients with low-risk versus high-risk signature scores was illustrated.

Collection of genomic data and clinical information

Gene expression profiles and clinical characteristics of 85 OS samples were retrospectively acquired from the TARGET project, which was defined as the discovery cohort. For the validations, 53 OS samples and 262 sarcoma samples respectively derived from the GSE21257 cohort and TCGA project were employed. The detailed clinical features and sequencing platforms were shown in Supplementary Table 3. Two immunotherapeutic cohorts with gene expression data were utilized to confirm the ICI roles of the identified immune signature. A total of 348 patients with advanced urothelial cancer in the IMvigor210 cohort [28], which was treated with atezolizumab (anti-PD-L1 agents), were collected; and their transcriptome data and clinical indexes were acquired from http://research-pub.gene.com/IMvigor210CoreBiologies. Another immunotherapy cohort contained 121 melanoma samples [29], which received anti-PD-1/PD-L1 or combination treatments, was also curated for further validating the risk signature response to ICI treatments. The detailed clinical characteristics of the two ICI cohorts were illustrated in Supplementary Table 4 for urothelial cancer and Supplementary Table 5 for melanoma. The list of comprehensive immune-related genes, containing a total of 1793 genes, was obtained from the ImmPort database [14] and shown in Supplementary Table 6.

Construction of the immune signature by the LASSO regression model

Univariate Cox regression analysis was used to determine the prognostic immune genes based on the expression profiles of 1793 genes in the TARGET cohort. After obtaining the above prognosis-related genes, the LASSO regression model with 10-fold cross-validation (implemented by R glmnet package) was performed to identify the final genes mostly contributive to the survival model. By combing the identified specific gene expression values and corresponding coefficient, the immune risk signature was obtained. In the subsequent analyses, patients were stratified into the high-risk and low-risk subgroups with the median risk score of the immune signature as cut-off value.

Quantify lymphocyte infiltration and immunotherapy response predictor enrichment

A study published by Angelova et al. proposed an 812-immune-metagene to conclude distinct tumor infiltration of 31 immune cell subtypes and depict a comprehensive immune landscape [24]. In this work, by applying the specific feature genes (Supplementary Table 7) that belong to distinct immune cells, we calculated the infiltration scores of each cell type for each sample. CIBERSORT is the most frequently used method to evaluate lymphocyte infiltration in cancer research [23]; nevertheless, in the discovery cohort (i.e., TARGET dataset), a majority of samples could not obtain a spendable estimation value. We therefore employed Angelova et al. algorithm to evaluate immune infiltration.

Recent several studies have demonstrated the crucial roles of microenvironment-based signals for assessing cancer immunotherapy efficacy. The IFNγ-related signature, which is a signal locates in the central site of the antitumor immune and predicts the clinical response to PD-1 blockade [16]. T cell-inflamed gene expression profiles (GEP) signature, which is consisted of 18 inflammatory regulators linked with immune response [16]. Cytolytic activity signature [17] is defined based on the transcriptomic expression levels of two vital cytolytic modulators, GZMA and PRF1, which are markedly upregulated upon CD8⁺ T cell activation [30] and during productive clinical responses to anti-CTLA-4 and anti-PD-L1 treatment [31, 32]. According to the representative genes of the above signatures (Supplementary Table 8), we inferred the enrichment scores for each sample and compared them in low-risk versus high-risk OS subgroups.

GSEA and GSVA

Differential expression analysis of low- versus high-risk groups was conducted with the limma package [33]. The calculated t statistics of all genes were extracted and subsequently as input to perform gene set enrichment analysis (GSEA), which is embedded in the fgsea R package. Signaling pathways acquired from the Hallmark database were utilized as the annotation and reference circuits. The false discovery rate (FDR) and normalized enrichment score (NES) were calculated based on 1 million permutations. Single sample gene set enrichment analysis (ssGSEA) method within the GSVA package [34] was used to evaluate the enrichment scores of collected immune cells and signatures based on the specific feature gene panels.

Definition of mutation burden and ICI responders in the immunotherapy cohorts

Tumor mutation burden (TMB) and neoantigen burden (NB) were defined as the log2 transformation of total mutation counts and neoantigen counts, respectively. For defining patients who responded to ICI treatments, patients with complete response (CR) or partial response (PR) were regarded to be efficacious to immunotherapy, and the rest were non-responders.

Statistical analyses

R software (version 4.0.2) was utilized to perform all related calculations and plots. Gene expression profiles of identified 21 immune genes illustrated by heatmap were achieved with pheatmap R package. Kaplan-Meier method was applied to produce survival curves and the Log-rank test to compare the differences. Multivariate Cox regression model with clinical confounding factors taken into account was conducted under R forestmodel package. The relationship of two distinct groups with categorical variables was evaluated with Fisher exact test, while the continuous variable was Wilcoxon rank-sum test. Association between multiple subgroups and continuous variables was calculated with the Kruskal-Wallis H test. Two-sided P values less than 0.05 were considered to be statistically significant unless a particular specification.

Data availability

All data used in this study are acquired from publicly available OS cohorts.