Data Processing
The expression data and clinicopathological information of samples in training cohort and validation cohort were download respectively. The clinicopathological characteristics of the two cohorts are detailed in Table 1. All 153 apoptosis-related pathways or biological processes were downloaded and corresponding 1807 nonredundant genes were used to take intersection with training and validation cohorts. Ultimately, 1476 mutual genes were identified.
The Construction of ARS
Univariate Cox regression analysis was performed for 1476 mutual genes, and 180 genes with P-value < 0.05 in univariate analysis were included in the LASSO regression analysis (Table S2) to construct the ARS. The ten-fold cross-validation was performed to determine the penalty parameter (λ) of the model. A total of 22 genes (RPS6, IFITM3, GRN, ATF4, MYC, DYNLL2, G6PD, BNIP3, PTGIS, BCL10, TRIM32, MAGEA3, PDK2, NMNAT1, EN1, TNFRSF11B, PPARG, SYNGAP1, GAL, GRIK2, MCF2, TERT) were included in the LASSO model (Figure 2). The coefficient of each gene was shown in Table S3.
The Power of ARS in Different Cohorts
In training cohort, the high-score group exhibited worse overall survival data as compared to the low-score group, with P-value < 0.01 (Figure 3E). Meanwhile, the distribution patterns of risk scores and survival status were plotted (Figure 3A and 3C). With the increasing ARS score, the overall survival time decreased and mortality increased. The 3- and 5-years AUC (AUC = 0.937, 0.947) of ARS was shown in Figure 3G, calculated by the “survivalROC” package.
These results indicate that the higher ARS represent the worse prognosis, and ARS has a good ability to predict the prognosis of osteosarcoma patients. The same result was also shown in the validation cohort. we calculate the ARS in validation cohort according to the LASSO formula. The high-score group exhibited worse overall survival as compared with the low-score group (P-value < 0.05; Figure 3F). The trend of survival time and mortality was the same as in the training cohort (Figure 3B and 3D). We found that the 3- and 5-years AUC (AUC = 0.771, 0.737) of ARS could still accurately predict survival state of osteosarcoma patients in validation cohort (Figure 3H).
Relationship Between ARS and Clinicopathological Features
To explore the relationship between ARS and clinicopathological features including age, sex, location and metastasis, we plotted the Kaplan-Meier curve in a different subgroup of clinicopathological features. As shown in Figure 4, the high-score group exhibited worse overall survival compared with the low-score group in each subgroup except in location of hand.
ESTIMATE, Immune Profile and Immune Checkpoint Molecules
The ESTIMATE algorithm revealed the immune (-1439.104 to 2560.632) and stromal scores (-689.9769 to 1927.611) of the training cohort (Figure 5A). It can be seen from the heatmap that there are significant differences in tumor purity, immune score and stromal score between the two groups. To further elucidate the relationship between ARS and immune/stromal score, both of the scores in ARS subgroups were compared. The Figure 5B-C shows that both immune scores and stromal scores in the low-score group were higher than that in the high-score group. This verified that immune score and ratio of tumor to stroma were significantly associated with ARS.
CIBERSORT package was performed to calculate the expression of 22 immune cells in different ARS groups. The result showed that NK cells resting is significantly more prevalent in the high-score group, while Macrophages M1 is significantly more prevalent in the low-score group (Figure 6).
To analyze the expression levels of immune checkpoint proteins in ARS subgroup, the expression of PD1, PDL1, and CTLA4 was calculated. The samples with low-score had higher expression of CTLA4 and PDL1 (Figure 6B-C), suggesting that patients with low-score ARS may respond better to immune checkpoint inhibitors targeting CTLA4 and PDL1.
Construction of a Prognostic Nomogram
The clinicopathological features and ARS of training cohort were used to construct a prognostic nomogram. Ultimately, the metastasis state and ARS were identified as the independent factors for prognosis of osteosarcoma according to univariate and multivariate cox regression analysis by the criterion of P-value < 0.05. A prognosis-related nomogram was constructed by the “rms” package with independent factors (Figure 7C). The tAUC of clinicopathological features, ARS and nomogram were shown in Figure 7B. The figure shown that the prediction ability of ARS (AUC = 0.887, 0.976, 0.937, 0.937, 0.947) and nomogram (AUC = 0.932, 0.984, 0.939, 0.939, 0.948) was higher than all other clinicopathological features. Then, the calibration curves of 1-, 3- and 5-year also shown that the prediction ability of the nomogram was very stable in different time points (Figure 7D-F).