1. Identification Of Inflammatory Response-associated Demrnas Related To Osteosarcoma Metastasis
The mRNA expression matrix of GSE21257 which included 34 metastatic and 19 non-metastatic osteosarcoma patients were downloaded from GEO database. 247 up-regulated and 217 down-regulated DEmRNAs related to osteosarcoma metastasis were identified (Fig. 2A). Two hundred inflammatory response-associated genes were obtained from hallmark gene sets in MSigDB and then overlapped with previously obtained DEmRNAs (Fig. 2B). In the end, 24 inflammatory response-associated DEmRNAs were screened out (Table 1).
Table 1
The inflammatory response-associated DEmRNAs in GSE21257
Gene
|
Log2FC
|
P Value
|
PDPN
|
-0.72
|
0.041
|
AXL
|
-0.75
|
0.0079
|
BST2
|
-0.69
|
0.022
|
C5AR1
|
-0.62
|
0.00014
|
STAB1
|
-0.78
|
0.0021
|
CXCL9
|
-0.94
|
0.015
|
TNFRSF1B
|
-0.79
|
0.018
|
CYBB
|
-0.88
|
0.031
|
CD48
|
-0.82
|
0.00028
|
CXCL10
|
-1.26
|
0.0034
|
LCP2
|
-0.77
|
0.00060
|
IL10RA
|
-0.71
|
0.000084
|
SERPINE1
|
-0.77
|
0.023
|
CLEC5A
|
-0.90
|
0.00000013
|
CCL5
|
-0.79
|
0.0046
|
OLR1
|
-0.91
|
0.0012
|
GNA15
|
-0.77
|
0.0020
|
IL1B
|
-0.85
|
0.000077
|
MYC
|
0.63
|
0.012
|
PLAUR
|
-0.80
|
0.00032
|
CCL2
|
-0.71
|
0.016
|
RGS1
|
-0.81
|
0.0090
|
FPR1
|
-0.82
|
0.00025
|
CD14
|
-1.29
|
0.000057
|
DEmRNA:differentially expressed messenger RNA |
2. Construction and validation of the inflammatory response-associated prognostic signature in GEO and Target databases
Based on the 24 DEmRNAs acquired previously, univariate and multivariate cox regression analyses were utilized to establish the inflammatory response-associated prognostic signature. In univariate cox regression analysis, three genes were associated with osteosarcoma prognosis (Fig. 3A). Then two of these three DEmRNAs were screened out to construct the inflammatory response-associated prognostic signature according to the multivariate cox regression result (Table 2). The risk score formula was established through mRNAs’ regression coefficients multiplied by expression levels:
Table 2
The multivariate cox regression analysis of GSE21257
Gene
|
Coefficient
|
HR
|
HR 95Low
|
HR 95High
|
P Value
|
CLEC5A
|
-0.96
|
0.38
|
0.18
|
0.83
|
0.015
|
MYC
|
0.41
|
1.51
|
0.92
|
2.49
|
0.11
|
HR:hazard ratio |
Risk score= (0.41*MYC) + (-0.96*CLEC5A)
Based on the median value of risk score, 26 patients were assigned into high-risk group and the rest 27 were regarded as low-risk (Fig. 3B). At the same time, patients in high-risk group showed higher mortality rates and shorter survival time than low-risk patients (Fig. 3c). The KM survival curve also indicated worse survival rates existed in high-risk patients (Fig. 3D). Finally, the area under the time-dependent ROC curves (AUC) were 0.65 at 1-year, 0.80 at 3- and 5-year, which demonstrated the robustness of this inflammatory response-associated prognostic signature (Fig. 3E).
The reliability of the inflammatory response-associated gene signature was validated in Target database. Ninety-five patients with integrated follow-up information were included in the test. Based on the median value of risk scores, 48 patients were separated into low-risk group and the rest were considered as high risk (Fig. 3F). The risk plot (Fig. 3G) and KM survival analysis (Fig. 3H) showed that low-risk patients had better prognoses. In addition, the time-dependent ROC curve revealed the signature was robust (Fig. 3I).
3. Exploration of the relationship between clinical characters and inflammatory response-associated prognostic signature
The prognostic value of gene signature was investigated through univariate and multivariate cox regression analyses. Figure 4A and 4B showed that the risk score was the independent prognostic factor in osteosarcoma. The same result was validated in Target database. (Supplementary Fig. 1A and 1B).
After that, the correlations between clinical features and genes in signature were investigated. In GSE21257, CLEC5A and MYC were related to the metastasis status and inflammatory response-associated risk (Fig. 3E,3F,3I,3J). In Target database, CLEC5A and MYC were only correlated to risk scores (Supplementary Fig. 1F and 1J). The KM survival analyses indicated the expression of CLEC5A positively related to the survival rates of osteosarcoma patients in GSE21257, Target and online R2 databases (Supplementary Fig. 2A,2C,2E). Instead, the higher expression of MYC related to the worse prognosis of osteosarcoma patients (Supplementary Fig. 2B,2D,2F).
4. Establishment of the nomogram
The nomogram was a powerful tool that integrated clinical variables to evaluate the prognosis of patients. In this research, the nomogram was constructed in GSE21257 and Target databases to predict the 3- and 5-year survival rates of osteosarcoma patients. (Figure 5A and Supplementary Figure 3A). The 3-year and 5-year calibration curves were close to the standard curve in two datasets (Figure 5B,5C and Supplementary Figure 3B,3C). The C-index was 0.99 in GSE21257 and 0.77 in Target dataset, which indicated the high accuracy of the predictive model.
5. Functional analysis of DEmRNAs related to inflammatory response-associated prognostic signature
The DEmRNAs among high- and low-risk patients were explored in GSE21257 and Target databases. A total of 736 DEmRNAs were screened out in GSE21257 and 1619 DEmRNAs were identified in Target database (Figure 6A and 6B). Then 216 DEmRNAs overlapped in two databases were selected for further analysis (Figure 6C).
The biological functions of these DEmRNAs were investigated through GO functional and KEGG pathway analyses (Figure 6D and 6E). A number of inflammatory response-associated biological processes such as neutrophil activation, regulation of leukocyte activation and antigen processing and presentation were significantly enriched. Furthermore, some immune-related functions (i.e., regulation of immune effector process, Th17 cell differentiation) were also enriched, which implied the potential relationships between inflammatory response and immune alterations in osteosarcoma.
6.Evaluation of the correlation between inflammatory response and immune infiltration
The relationship between tumor microenvironment and inflammatory response-associated risks was evaluated through ESTIMATE software. The results showed that the risk score was negatively related to the immune and stromal scores in GSE21257 and Target datasets (Figure 7A,7B and Supplementary Figure 4A,4B). Then the immune infiltration of 22 immune cells was investigated by CIBERSORT algorithm. The statistically significant abundance ratios of immune cells in osteosarcoma were shown in Figure 7C and supplementary Figure 4C. Then the specific relationships between immune functions and inflammatory response-associated risks were evaluated through ssGSEA analysis. The outcome revealed that the vast majority of immune cells (14) and functions (11) were down-regulated in high-risk group in GSE21257 (Figure 7D,7E). At the same time, eight immune cells and twelve immune functions were dysregulated in Target dataset (Supplementary Figure 4D,4E). The overlapped items including eight immune cells and eleven immune functions were considered as the alterations of inflammatory response-associated immune infiltration (Figure 7F).
7. Validation of the inflammatory response-associated genes in cell lines
The expressive levels of inflammatory response-associated genes in human osteosarcoma and osteoclast cell lines were detected through qRT-PCR analyses. The results showed that MYC was upregulated and CLEC5A was downregulated in all osteosarcoma cell lines (Figure 8A,8B). Then osteosarcoma cell lines were divided into high- (143B, WELL5, MNNG/HOS) and low-metastatic (U20S, SAOS2, MG63) groups according to their metastasis ability. The MYC was up-regulated in high-metastatic osteosarcoma cell lines and CLEC5A showed the opposite trend but with no statistically significant (Figure 8C,8D).