Metastasis-related gene signature associates with immunity and predicts prognosis accurately in patients with osteosarcoma

Identification of MRGs

Based on metastatic and non-metastatic patients, we identified 28 DGEs. We found 16 upregulated and 12 downregulated genes (Figure 1A). The expressions of the 28 MRGs are shown in Figure 1B. Next, the upregulated genes were enriched in this pathway. GO analysis showed that it was significantly enriched in bone development, bone morphogenesis, and regulation of Wnt and kinase signaling pathways. Enrichment analysis of the KEGG pathway showed that MRGs were mainly involved in the Wnt signaling pathway, the JAK-STAT signaling pathway, and cytokine-cytotoxic receptor interactions (Figure 1C), which are known to drive the metastasis of other cancers.

Figure 1. Identification of metastasis-related genes and functional enrichment analysis. (A) Volcano plot showing differentially expressed genes between primary and metastatic osteosarcoma. (B) Heatmaps of metastasis-related genes in the Target Database. (C) GO analysis of the metastasis-related genes. KEGG pathway enrichment analysis of the metastasis-related genes.

Expression pattern of OS patients based on MGRs

To explore the relationship between MGRs and the prognosis of OS, we performed a survival analysis of each MGR in patients with OS from the TARGET database. High levels of SEMA5A, HILPDA, WIF1, ALPL, CXCL8, LOX, and SEC61G were significantly associated with poor overall survival in OS, whereas ASPA, RBM24, GPR174, NLGN1, SLFN13, IGF2, and PDLIM3 showed the opposite trend (Figure 2A). The interaction network showed the interactions, regulatory relationships, and risk factors between MGRs in patients with OS (Figure 2B).

Figure 2. The prognostic value of MRGs in patients with OS. (A) Kaplan-Meier curves for the MRGs in OS patients from TARGET database. (B) A network of correlations including MRGs in the TARGET cohort. (p < 0.05 *; p <0.01 **; p < 0.001 ***).

Construction of prognosis model and Norman diagram based on MGRs and identification and functional analysis of MRS

First, based on the 28 MGRs, we conducted univariate Cox and LASSO region analyses to identify five risk genes (Figure 3A, 3B). The risk score is calculated as follows:

$$\begin{array}{l}Risk Score=(0.272989391508956\\ \times LOXExp)\\ +(0.18788859580782\\ \times WIF1Exp)\\ +(0.232621633568678\\ \times SEMA5AExp)\\ +(0.282143937841433\\ \times HILPDAExp)\\ +[(-0.343452464841895)\\ \times NLGN1Exp]\end{array}$$

According to Kaplan-Meier analysis, the OS of the low-risk group was considerably superior to that of the high-risk group (Figure 3C). The ROC curve also demonstrated that the model was thoroughly validated. The AUCs of the validation sets at three, five, and seven years were 0.711, 0.724, and 0.679, respectively. (Figure 3D). Additionally, a clinical correlation heat map and risk model revealed that the high-risk group had higher LOX, HILPDA, WIF1, and SEMA5A expression levels, whereas the low-risk group had higher expression levels of NLGN1.

Figure 3. The construction of MRSs and enrichment analysis. (A) Forest plot of the univariate Cox regression analysis for MRGs. (B) LASSO regression and stepwise multivariate Cox regression were used to construct the five-MRG signature. (C) Kaplan-Meier survival analysis of patients. (D) ROC curve analysis of the five-MRG signature. (E) Expression patterns of MRSs in high-and low-risk groups. (F, G) Forrest plot of the independent prognostic factors in OS. (H–J) Nomogram for predicting the 1-, 3-, and 5-years OS of osteosarcoma patients. (I) Calibration curves for validating the established nomogram. (J) The ROC curves of the nomograms compared for 1-, 3-, and 5-years OS in osteosarcoma patients, respectively. (K, L) Circle plot through Gene Ontology (GO) analysis visualizing the biological processes enriched by DEGs. (L) GO analysis of differentially expressed genes between the high- and low-risk groups. (M) Differences in risk scores between distinct survival status and the proportion of survival status of patients.

Additionally, more individuals had metastatic disease in the high-risk group than in the low-risk group (Figure 3E). Univariate and multivariate COX analyses revealed that metastasis was a standalone prognostic risk factor for OS (P=0.01; hazard ratio [HR] =3.168) (Figure 3F, 3G). The Norman diagram includes independent prognostic criteria (Figure 3H). The calibration curve demonstrated some congruence between the observed data and the projected survival probability of the overall survival times of 1, 3, and 5 years (Figure 3I). Additionally, a comparison of the ROC curves demonstrated that the nomogram was more accurate than any clinical factor in forecasting patients’ long-term survival (Figure 3J). There were 135 MRS-related differential genes, and 39 genes were enriched in multiple pathways. These included extracellular matrix organization, regulation of cell-subtract adhesion, extracellular structure organization, response to amyloid-beta, ossification, extracellular matrix assembly, cell-subtract adhesion, and multicompolysaccharide metabolism processing (Figure 3K). In addition, BP, CC, and MF reflected the extensive heterogeneity of OS metastatic capacity (Figure 3L). We found that the group with a high-risk score corresponded to more patients with metastasis, death, and a poor prognosis (Figure 3M).

MRS based comprehensive immunoassay and drug screening

The association between the risk score and the number of immune cells was assessed using the CIBERSORT method. The low-risk group exhibited greater immune cell infiltration than the high-risk group (Figure 4A). Additionally, the low-risk group had higher Matrix and ESTIMATE scores than the high-risk group (Figure 4B). We hypothesized that the lower risk of OS in the low-risk group may be due to greater immune cell infiltration compared to the high-risk group, which showed lower levels of immune cell infiltration. Therefore, it is necessary to explore the potential use of immune checkpoint therapies. The results of the correlation between the risk score and immune checkpoints showed that the risk score was negatively correlated with CD86 and CD27 (Figure 4C). Next, the results of drug correlation and IC50 difference showed that low-risk patients were more sensitive to Axitinib, Docetaxel, and Nilotinib, whereas high-risk groups were more sensitive to Rapamycin, Imatinib, and DMOG (Figure 4D).

Figure 4. Distinct TME characteristics and mutation of OS patients according to the risk score. (A) Abundance of 23 infiltrating immune cell types. (B) Correlations between risk score and both immune and stromal scores. (C) Correlations between risk score and immune checkpoints. (D) Predicts the responsiveness of OS to chemotherapy.

Building unsupervised consensus clustering and comprehensive immune analysis of MRGs

Based on the mRNA expression levels of the 28 MGRs, OS samples in the TARGET database were clustered using unsupervised non-negative matrix decomposition. The best k value was determined using molecular clustering and comprehensive correlation analysis. We examined the k values between 2 and 9. For clustering and overall survival analysis, k = 3 was the most appropriate of all possible values (Figure 5A, 5B). Clusters A, B, and C of the MRGs were considered subgroups. Significant variations in MGRs expression and clinical traits were observed among the three groups when clinical factors and MGRs expression were compared. Among them, MRG cluster A had a good prognosis, whereas MRG cluster C had a poor prognosis. In addition, we will combine clustering with clinical features for the analysis. The findings demonstrated a substantial difference between MRG clusters A and C in OS metastasis. In addition, the expression patterns of LOX, HILPDA, WIF1, SEMA5A, and NLGN1 were consistent with the outcomes of the risk model, and the number of patients in MRG cluster A was much lower than that in MRG cluster C (Figure 5C). This is also consistent with the two survival curves that we analyzed earlier. In addition, we examined the impact of the MRG clusters on the TME and found that the infiltration of CD56dim natural killer cells increased in MRG cluster Group C (Figure 5D).

Figure 5. Identification of metastasis subtypes in OS. (A) t-SNE of the mRNA expression profiles of MRGs from the OS samples in the TARGET dataset confirmed the three clusters: A, B and C. (B) Kaplan-Meier curves for the three molecular patterns of OS patients. (C) Heatmap depicted the correlation between the subtypes and different clinicopathological characteristics. (D) Boxplots showed abundance of 23 infiltrating immune cell types.

Generation of gene cluster and comprehensive analysis of immunity

Gene clusters A and B were created based on the DEGs between the risk categories (Figure 6A). According to the survival study, patients in group B had a poorer prognosis than those in group A (Figure 6B). This previous study was consistent with the gene cluster expression pattern (Figure 6C). These genes directly affect the development and incidence of OS. Similarly, the TME characteristics of each gene cluster showed that almost all statistically significant immune cells differed significantly between the two gene clusters. Gene cluster A showed more abundant immune cell infiltration than cluster B (Figure 6D). This result shows that the TME characteristics of the gene cluster are consistent with those of the MRG clusters and the risk model.

Figure 6. Prognosis and TME characteristics in two metastasis gene clusters for OS patients. (A) Consensus matrix heatmap defining two gene clusters according to the prognostic DEGs. (B) Kaplan-Meier survival analysis for patients in the two gene clusters. (C) Clinical features of the two gene clusters. (D) Boxplots showed abundance of 23 infiltrating immune cell types.

Construction of MRGs score and comprehensive immune analysis in patients with OS

The MRGs score and the distribution and survival of the gene and MRGs clusters are represented in the Sanggi diagram (Figure 7A). In addition, we have made further verifications. The expression of MRG scores in the gene and MRG clusters (Figure 7B, 7C). Subsequently, the survival rate was the same as before. A high MRGs score indicated a poor prognosis (Figure 7D). Interestingly, the expression patterns of LOX, HILPDA, WIF1, SEMA5A, and NLGN1 were consistent with those of previous analyses (Figure 7E). Additionally, immune cell infiltration in the low MRGs score group was more noticeable according to TME characteristic data (Figure 7F). Additionally, the matrix score of the low-risk group was higher than that of the high-risk group, and the immune score of the high-risk group was higher than that of the low-risk group (Figure 7G). The group with a high metastasis score corresponded to a greater number of patients with metastasis, death, and a poor prognosis (Figure 7H). This conclusion is consistent with that of the risk-scoring group.

Figure 7. Development and validation of the metastasis scoring system and distinct TME characteristics for OS. (A) Sankey Diagram of metastasis clusters, gene clusters, metastasis score, and clinical outcomes. Differences in metastasis score between (B) the three metastasis subtypes and (C) the two gene clusters. (D) Kaplan-Meier analysis of the OS between the two metastasis score groups. (E) Heatmap depicted the correlation between the metastasis score and different clinicopathological characteristics group. (F) Boxplots showed abundance of 23 infiltrating immune cell types. (G) Correlations between risk score and both immune and stromal scores. (H) Differences in metastasis scores between distinct survival status and the proportion of survival status of patients.

Verification of risk prediction model in external queue

We used two external queues to evaluate the reliability and stability of five MRS predictions. The risk score for each patient on these lines was determined using the same method. The patients were then divided into high- and low-risk groups according to the previously set median of the training queue. The survival rates of the two patient groups (Figure 8A, 8D) and the amount of MRS expression were mostly compatible with the expression mode of the training set risk model (Figure 8B, 8E). In the GSE21257 cohort, the area under the curve (AUC) for 1-, 3-, and 5-year survival rates were 0.673, 0.693, and 0.629, respectively. In the GSE39055 cohort, the (AUC) for 1-, 3-, and 5-year survival rates were 0.729, 0.521, and 0.518, respectively (Figure 8C, 8F).

Figure 8. External verification of models OS and MRSs expression pattern. (A–C) Validation of external validation set GSE21257. (D–F) Validation of external validation set GSE39055.

Differential expression of HILPDA, LOX, NLGN1, SEMA5A and WIF1 in hFOB1.19 and MG63 cells

We used qRT-PCR and western blotting to study the expression of HILPDA, LOX, NLGN1, SEMA5A, and WIF1 in MG63 cells. Compared to hFOB1.19, HILPDA, LOX, SEMA5A, and WIF1 showed increased mRNA and protein expression in MG63 cells, whereas NLGN1 showed the opposite trend (Figure 9A, 9B).

Figure 9. Expression difference of HILPDA, LOX, SEMA5A, WIF1, NLGN1 in hFOB1.19 cells and MG63 cells. (A) Compared with hFOB1.19 cells, the mRNA levels of HILPDA, LOX, SEMA5A and WIF1 in MG63 cells increased significantly, while NLGN1 showed the opposite. (B) Compared with hFOB1.19 cells, the protein levels of HILPDA, LOX, SEMA5A and WIF1 in MG63 cells increased significantly, while NLGN1 showed the opposite.(p < 0.05 *; p <0.01 **; p < 0.001 ***).

Knockdown of HILPDA, LOX, SEMA5A, WIF1 or overexpression of NLGN1 can effectively inhibit migration and invasion of MG63 cells

The scratch test was used to analyze the migration of MG63 cells. Compared to the NC group, the migration area of the si-HILPDA, si-LOX, si-SEMA5A, si-WIF1, and oe-NLGN1 groups decreased significantly, indicating that the low expression of HILPDA, LOX, SEMA5A, WIF1, and other genes and the overexpression of NLGN1 effectively inhibited the migration ability of MG63 cells (Figure 10A).

Figure 10. Knock out HILPDA, LOX, SEMA5A, WIF1 or overexpress NLGN1 to observe the migration and invasion ability of MG63 cells. (A, B) Knockdown of HILPDA, LOX, SEMA5A, WIF1 and overexpression of NLGN1 can reduce the migration and invasion ability of MG63 cells. (p < 0.05 *; p <0.01 **; p < 0.001 ***).

We used a Transwell assay to evaluate the invasive ability of MG63 cells. Compared to the NC group, the number of MG63 cells passing through the chamber in the si-HELPDA, si-LOX, si-SEMA5A, si-WIF1, and oe-NLGN1 groups decreased significantly, indicating that the low expression of HILPDA, LOX, SEMA5A, WIF1, and other genes and the overexpression of NLGN1 can effectively inhibit the invasiveness of MG63 cells (Figure 10B).

Collection of OS datasets

The OS database was obtained from the TARGET and GEO databases. The TARGET database contained 88 tumor samples, GSE21257 contained 53 samples, and GSE39055 contained 37 samples. The TARGET database was used for the training set, and GSE21257 and GSE39055 datasets were used for the verification set.

Screening of transfer related differential genes

The differential expression of mRNA was examined using the “limma” software package of R (version: 3.4.1), which was based on the initial count of the TARGET dataset and the associated clinical data. The screening of threshold mRNA differential expression was specified as “adjusted p = 0.05 and Log₂ (fold change) >1 or Log₂ (fold change) < −1“. We investigated the upregulated genes and conducted a functional enrichment analysis of the data to further validate the possible functions of the potential targets. The “ClusterProfiler” package in R was used to examine the GO function of prospective mRNA and enhance the KEGG pathway in order to better comprehend the carcinogenesis of the target gene.

Consensus clustering analysis

Consensus unsupervised clustering analysis was carried out using the R package “CONSENSUSClusterPlus” based on the expression profiles of 28 MRGs. OS patients in the TARGET database were divided into three molecular subgroups of MRGs: clusters A, B, and C.

Relationship between MRGs: clinical features and prognosis of OS

We evaluated the correlation among OS clustering, clinicopathological features, and survival outcomes to assess the clinical importance of various clustering approaches. The Kaplan-Meier technique was used to compare the overall survival rates.

Construction of risk label and Norton

The weighting coefficient was determined using least absolute contraction and selection operator (LASSO) Cox regression analysis, and the risk model was created in R using the glmnet package. Based on the best critical value of the risk score, patients were separated into high- and low-risk groups, and the Kaplan-Meier technique was used to determine the survival rate of the patients. To define the ROC curve and obtain the AUC for survival analysis, use the “Survival ROC” tool in R. Finally, a Norman Chart was updated to include clinical and independent prognostic variables.

Identification and enrichment analysis of differentially expressed genes based on risk model

To find DEGs between various risk score groups, R’s “limma” package is utilized. Genes were deemed to be expressed differently when their adjusted p-values were 0.05, and | log2 (FC) | was greater than 1.0. Finally, three DEGs were predicted to be prognostic genes and analyzed. The “clusterProfiler”, “enrichplot”, and “ggplot2” packages of R were used for the enrichment analysis of MRS-related differential genes in OS. The screening criteria were p<0.01 and q<0.05.

Construction of MRGs score

Based on the three prognostic DGEs, we calculated the MRG scores for each patient with OS using PCA. Subsequently, the MRGs gene score was established using PCA, and principal components 1 and 2 were extracted as feature scores. The score formula was as follows: The MRGs score = Σ (PC1i+PC2i).

TIME cell infiltration evaluation

The tumor immune estimation resource (TIMER) was utilized to evaluate the level of immune cell infiltration between different clusters using the Estimation of Stromal and Immune cells in malignant tumor tissues using Expression data (ESTIMATE) method. Single-sample gene set enrichment analysis (ssGSEA) was performed to determine the relative abundance and activity levels of each immune cell type in OS, followed by CIBERSORT analysis to assess the fraction of invading immune cells in OS patients with specific expression patterns, in order to explore differences between immune cell subtypes.

Immune checkpoint and drug sensitivity of risk score

The immune checkpoint inhibitor response to a single sample or subtype was predicted using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm [15]. The medications sensitive to high-risk populations are screened based on the risk score and the IC50 value using the R “pRRophic” package.

Cell transfection

Biomics Biotech created and manufactured si-HILPDA, si-LOX, si-SEMA5A, si-WIF, and their respective negative controls (NCs) (Nantong, China). Their sequences were listed in Supplementary Table 1. Genesyntech created and manufactured oe-NLGN1 (Homo sapiens neuroligin 1, transcript variant 1, mRNA; NCBI Reference Sequence: NM_001365923.2), and a negative control. Lipofectamine 3000 reagent was used to transfect MG63 cells (Invitrogen, Waltham, MA, USA) according to the manufacturer’s recommendations. The transfection efficiency is shown in Supplementary Figure 1.

RNA extraction and qRT PCR detection

Total RNA from hFOB1.19 and MG63 cells should be extracted using TRIzol (Invitrogen), and 2ug of total RNA should be reverse-transcribed into cDNA using the TransScript One Step gDNA Removal and cDNA Synthesis SuperMix kit (TransGen Biotech, Beijing, China). The ABI Prism 7500 rapid real-time PCR equipment and TransStart Green qPCR SuperMix kit (TransGen Biotech, Beijing, China) were used to measure the expression levels of HILPDA, LOX, NLGN1, SEMA5A, and WIF1 (Applied Biosystems, StepOnePlus, Thermo Fisher Scientific, Waltham, MA USA). Actin was used as the internal reference for mRNA at the same time. Shanghai Sangong Biotechnology Co., Ltd. created all the primers (Shanghai, China). The primer sequences used in the experiments are listed in Table 1.

Western blot analysis

Radioimmunoprecipitation (RIPA) lysis buffer was used to fully lyse the cells, and the resulting protein samples were collected, separated by electrophoresis, and transferred onto PVDF membrane (Millipore, Billerica, MA, USA). The PVDF membrane was then sealed with a fast-sealing solution and incubated with primary antibodies against HILPDA (diluted 1:1000; Affinity), LOX (diluted 1:1000; Proteintech), NLGN1(diluted 1:10000; Proteintech), SEMA5A (diluted 1:2000; Bioss), WIF1 (diluted 1:3000; Proteintech), and GAPDH(diluted 1:50000; Proteintech) at 4° C. The corresponding secondary antibodies were incubated at room temperature for 1.5 h. Finally, the protein bands were visualized using an ECL kit.

Cell migration and invasion assays

The scratch test was used to analyze the migratory ability of MG63 cells. After the cells were digested, the seeds were spread onto a well plate. After the cells were digested, the seeds were spread onto a well plate. When the degree of cell fusion reached 95%, a sterile plastic tip was used to scratch the cells. After PBS washing (1-2 times), the cells were cultured in serum-free medium for 24 h. The test results were analyzed by Image J, and the migration areas of the cells in each group were compared.

The invasive ability of the MG63 cells was analyzed using a cross-pore chamber. According to 1 × 106 cells/ml, 200ul ratio was plated in a Transwell chamber (8-um; Corning, NY, USA). The same cells or 500ul MEM medium containing 10% fetal bovine serum were used in the inferior chamber and incubated for 24 h. The cells were fixed with 4% paraformaldehyde for 30 min, and 0.5% crystal violet dye was used for 20–30 min. The cells in the upper chamber were gently wiped with a cotton ball and the migrated cells were photographed under a high-magnification field of view using an inverted microscope (×200).

Statistical analysis

All statistical calculations were carried out using R 3.4.1. One-way ANOVA was used to ascertain statistical significance in the experimental data, which are presented as the mean standard deviation. Statistical significance was set at a bilateral p value of 0.05.

Availability of data and materials

The original contributions presented in this study are included in the article/Supplementary Material. Further inquiries can be directed at the corresponding author.

Consent for publication

All authors have agreed on the contents of the manuscript.