The anoikis-related gene signature predicts survival and correlates with immune infiltration in osteosarcoma

Data processing

The gene symbols for 434 ARGs were obtained from the study by Chen S et al. [12]. The gene expression profiles from GSE16088 [13] and GSE28425 [14] datasets were used to identify DEGs using the “limma” R package [15]. The “sva” and “limma” R packages were used to integrate and standardize the GSE16088 and GSE28425 datasets. DEGs were defined as those with p <0.05 and log2|fold change| values >1. The “ggplot2” [16] and “pheatmap” [17] packages were used to generate volcano plots and heatmaps, respectively.

Data for the OS model were obtained from the TCGA and TARGET databases for training cohorts, and the validation cohorts’ clinical data and gene expression profiles were obtained from the GSE16091 [13] and GSE21257 [18] datasets. Integration and standardization of GSE16091 and GSE21257 datasets were implemented using the “sva” and “limma” R packages. For the RFS model, the gene expression profiles and clinical information for the training cohort were extracted from the GSE39058 [19] dataset, while the validation cohorts were downloaded from the TARGET and TCGA databases.

Identification of hub genes

The previously described method [20] was used to identify hub genes. The STRING database [21] was used to retrieve data on DEG protein-protein interactions (PPI). Using Cytoscape 3.7.2 software, the PPI network was developed and visualized. Candidate hub genes were identified using the three methods, MCODE, degree, and betweenness, and the hub genes were identified based on the overlap of the three methods’ results.

Prognostic model construction and validation

First, the OS prognostic model was constructed and validated. TCGA and TARGET data on mRNA expression and clinical information were collated and combined. The “survival” package was then used to perform Kaplan-Meier (KM) analysis and generate survival curves. Candidate model genes were those with p <0.05. The “glmnet” R package was then used to perform the least absolute shrinkage and selection operator (LASSO) Cox regression analysis to construct the OS model. Each patient’s risk score was also calculated. The risk score formula was as follows: RS= ∑Coef_genei x Exp_genei. Patients were classified into high-risk and low-risk groups based on their median risk scores. The survival curves of the two groups were plotted to compare OS. The “survivalROC” R package was used to plot the receiver operating characteristic (ROC) curve to evaluate the OS model’s predictive capability. Similarly, risk scores were calculated for patients in the validation (GSE16091 and GSE21257) cohort, and survival curves and ROC curves were visualized.

We used the same methods as above to construct the RFS model. The training cohort was GSE39058, and the validation cohort was TCGA. Finally, we developed nomograms using the “rms” package to predict OS and RFS more intuitively and conveniently.

GO enrichment, KEGG pathway, GSVA analysis and biological networks

We used the ‘clusterProfiler” R package to conduct GO and KEGG pathway analyses to further investigate the potential functions and mechanisms of hub genes and prognosis-related genes [22]. Then, we performed GSVA analyses to explore the reasons for the different prognosis of patients with high or low risk as previously described [23]. In addition, Transcription factor (TF)-microRNA (miRNA) coregulatory networks of prognosis-related and hub genes were established using NetworkAnalyst 3.0.

Immune scores and tumor-infiltrating immune cell analysis

The ESTIMATE algorithm calculated immune and stromal scores for the training (TCGA-TAEGET) cohort and validation (GSE16091 and GSE21257) cohort samples. Tumor-infiltrating immune cells (TIICs) were assessed in the training (TCGA-TAEGET) cohort and validation (GSE16091 and GSE21257) cohort using the CIBERSORT online tool [24].

Cell culture and transfection

From the Cell Bank of the Chinese Academy of Sciences (Shanghai, China), MG63 and U2OS cell lines were purchased. DMEM with 10% FBS (Biological Industries, Shanghai, China) was used to culture both cell lines. GeneChem (Shanghai, China) provided OGT-targeting siRNA and negative control siRNA (si-nc). Supplementary Table 1 provides comprehensive details about the si-1 and si-2 sequences. Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) was used to transfect the cells with the siRNAs.

Quantitative real-time PCR (RT-PCR) and Western blotting (WB)

Quantitative RT-PCR and WB were carried out as previously described [20]. RT-PCR reagents were purchased from Accurate Biology (Changsha, China), which included the SYBR Green Premix Pro Taq HS qPCR Kit, AG RNAex Pro Reagent, and Evo M-MLV RT Premix Kit. Supplementary Table 1 shows the primer sequences. Beyotime Biotechnology (Shanghai, China) provided WB reagents, such as BeyoECL Plus and RIPA buffer Kit. Proteintech (Wuhan, China) supplied anti-OGT (11576-2-AP) and anti-GAPDH (10494-1-AP) antibodies.

Cell migration and invasion assays

Transwell assays were performed as previously described [25] to assess the effects of OGT expression on MG63 and U2OS cell migration and invasion. Finally, cells that migrated into the lower chamber were counted using a microscope in five different fields.

Statistical analysis

Most bioinformatics and statistical analyses in this study were performed using R software (version 4.0.0), including integration and normalization of mRNA expression data, LASSO Cox regression analysis, survival analysis, ROC analysis, CIBERSORT, and ESTIMATE. The mean ± standard deviation of three independent experiments was used to represent all quantitative data obtained from in vitro experiments. GraphPad Prism 8.0 (GraphPad, La Jolla, CA, United States) was used to analyze differences across three groups using one-way ANOVA. p <0.05 denoted a statistical significance level.

Data availability statement

All datasets presented in this study are included in the article Supplementary Material.

DEGs and hub genes in osteosarcoma

The GSE16088 and GSE28425 datasets of normal bone and tumor tissues were used to identify differentially expressed ARGs. Figure 1A, 1B show the Volcano plot and heatmap of ARGs in GSE16088 and GSE28425, respectively. Finally, 71 ARGs out of 434 ARGs were identified as DEGs (Supplementary Table 2).

Figure 1. Differentially expressed genes and hub genes. (A) Volcano plot of 434 ARGs in GSE16088 and GSE28425 datasets. (B) Heatmap of DEGs in GSE16088 and GSE28425 datasets. (C) Protein-protein interaction (PPI) network of 71 DEGs and the critical module. Identifying the first 10 ARGs and constructing the corresponding PPI network using the degree (D) and betweenness (E) topological method. (F) Venn calculation is applied to identify seven hub ARGs.

After mapping the 71-DEG PPI network with STRING, the Cytoscape MCODE plug-in was used to reconstruct the network (Figure 1C). The first 10 genes were chosen, and the degree and betweenness topological approaches were used to construct the associated PPI network (Figure 1D, 1E). Seven genes—FN1, CD44, HRAS, TP53, PPARG, CTNNB1, and VEGFA—were shown to be hub genes by crossing the three methods’ results (Figure 1F).

Prognostic model construction and validation

First, we conducted KM analysis to investigate how DEGs affected the training (TCGA-TARGET) cohort’s OS before analyzing ARGs’ effects on osteosarcoma patients’ prognoses. The results revealed that OS was substantially correlated with nine candidate ARGs (Supplementary Figure 1). Subsequently, the OS model was constructed using LASSO regression analysis. Figure 2A shows the corresponding confidence interval for each lambda and each gene’s lambda-valued coefficient trajectory. Seven signature genes were chosen for the OS model. The following formula was used to calculate the risk score: (0.106*Exp CD24) + (0.135*Exp FASN) + (-0.206*Exp MMP2) + (0.272*Exp EIF2AK3) + (0.029*Exp ID2) + (-0.569*Exp PPARG) + (-0.243*Exp PIK3R3). The training cohort’s median risk score classified patients into low- and high-risk groups. According to the KM analysis results, patients in the high-risk group had considerably worse (p <0.001, Figure 2B). We next performed ROC analyses to evaluate aforementioned ARGs’ predictive capabilities. The outcomes showed that the OS predictive model’s area under the ROC curve (AUC) values were 0.838, 0.832, and 0.83 over 3, 5, and 7 years, respectively (Figure 2B). The GSE16091 and GSE21257 datasets were then used for verification analysis. Patients at high risk had shorter OS duration, in line with the findings of the training cohort (p <0.05, Figure 2C). Over 3, 5, and 7 years, the OS model’s AUCs were 0.653, 0.647, and 0.561, respectively (Figure 2C).

Figure 2. Development and validation of the OS-prediction model for osteosarcoma. (A) The least absolute shrinkage and selection operator (LASSO) method of ARGs is associated with prognosis. Survival curve and ROC curve for low- and high-risk subgroups in the training cohort (TCGA-TAGET) (B) and the validation cohort (GSE16091 and GSE21257) (C).

We used the same method to construct the RFS predictive model. GSE39058 sets were used as the training cohort. Eight candidates’ ARGs were strongly related to RFS, according to KM analysis (Supplementary Figure 2). The RFS model was constructed using LASSO regression analysis. Figure 3A shows the corresponding confidence interval for each lambda and each gene’s lambda-valued coefficient trajectory. Four signature genes were chosen for developing the RFS model. The following formula calculated the risk score: (0.416*Exp BRMS1) + (0.295*Exp COL4A2) + (-0.032*Exp FGF2) + (0.372*Exp OGT). The training cohort’s median risk score divided patients into low- and high-risk groups. According to the KM analysis results, patients in the high-risk group had considerably worse RFS (p <0.001, Figure 3B). Over 3, 5, and 7 years, the RFS model’s AUCs were 0.875, 0.805, and 0.867, respectively (Figure 3B). The validation cohort was then created using the TCGA-TARGET sets. Patients in the high-risk group had shorter RFS duration, which was consistent with the training cohort results (p =0.009, Figure 3C). Over 3, 5, and 7 years, the OS model’s AUCs were 0.741, 0.803, and 0.678, respectively (Figure 3C). The risk score and clinical parameters’ prognostic value were also assessed using univariate and multivariate Cox analyses. The findings revealed that patients in the TCGA-TARGET dataset had an independent risk score for OS (Supplementary Figure 3A) and RFS (Supplementary Figure 3B). In summary, osteosarcoma patients’ outcomes were accurately predicted by the OS and RFS models.

Figure 3. Development and validation of the RFS-prediction model for osteosarcoma. (A) The least absolute shrinkage and selection operator (LASSO) method of ARGs is associated with prognosis. Survival curve and ROC curve for low- and high-risk subgroups in the training cohort (GSE39058) (B) and the validation cohort (TCGA-TAGET) (C).

Building predictive nomograms

Nomograms (Figure 4A, 4B) were constructed to enable us to better predict osteosarcoma patients’ OS and RFS using prognosis-related ARGs. Each gene was assigned a corresponding point on the point scale. According to each patient’s predicted 3-, 5-, and 7-year survival rates, the points were then added up to determine the total score for each patient.

Figure 4. The nomograms for predicting 3-, 5-, and 7-year OS and RFS of the training cohort. (A) Nomogram for predicting 3- and 5-year OS in the training cohort (TCGA-TAGET) and calibration plots. (B) Nomogram for predicting 3- and 5-year RFS in the training cohort (GSE39058) and calibration plots.

Potential functions and mechanisms

GO and KEGG analyses of the hub and prognosis-related ARGs were carried out using the “clusterProfiler” R package to investigate their functions and possible mechanisms in osteosarcoma. As shown in Figure 5A, for GO enrichment analysis, the seven hub ARGs were enriched in cell growth, fibroblast proliferation, transcription regulator complex, protein-C terminus binding, and phosphatase binding. The prognosis-related ARGs were enriched in lipid metabolic process and smooth muscle cell proliferation (Figure 5B). For KEGG pathway analysis, the seven hub ARGs were enriched in virus infection and some cancers (Figure 5A). The prognosis-related ARGs were enriched in AGE-RAGE, AMPK, and Relaxin signaling pathways. (Figure 5B). Then, GSVA analyses were performed to explore the reasons for the different prognosis of patients with high or low risk. For OS model, the Ribosome pathway was activated primarily in high-risk patients, whereas low-risk patients mainly activated other pathways, such as leukocyte transendothelial migration, natural killer cell mediated cytotoxicity, and complement and coagulation cascades (Figure 5C). For RFS model, high-risk subgroups mainly activated cytokine- cytokine receptor interaction pathway, chemokine signaling pathway, NOD like receptor signaling pathway and TOLL like receptor signaling pathway, whereas low-risk subgroups mainly activated lysine degradation pathway (Figure 5D). Furthermore, the TF-miRNA coregulatory network of the seven hub ARGs (Figure 5E) and prognosis-related ARGs (Figure 5F) was constructed using NetworkAnalyst.

Figure 5. Function enrichment analysis. GO function enrichment and KEGG pathway analysis of hub genes (A) and prognosis-related ARGs (B). GSVA analysis for OS (C) and RFS (D) in low- and high-risk subgroups of TCGA cohort. Transcription factor-miRNA coregulatory network of hub genes (E) and prognosis-related ARGs (F).

Immune infiltration differences between risk groups

Because the tumor microenvironment (TME) influences malignant tumor progression [26], we used the ESTIMATE and CIBERSORT tools to assess the differences between risk groups in immune cells and immune infiltration. Both the validation (GSE16091 and GSE21257) and training (TCGA-TARGET) cohorts showed that the high-risk group had lower immune and stromal scores than the low-risk group, according to the ESTIMATE results (Figure 6A). Low immune scores correspond with low survival rates (Figure 6B). We then compared TIIC differences between risk groups. Results from the TCGA-TARGET cohort revealed that the high-risk group had considerably fewer infiltrating T cells CD4 memory activated than the low-risk group (Figure 6C). GSE16091 and GSE21257 cohort results showed that the high-risk group had significantly fewer infiltrating T cells gamma delta and M2 than the low-risk group, whereas the high-risk group had substantially higher T cells follicular helper, mast cells activated, plasma cells, and infiltrating B cells naïve (Figure 6D). Differences in immune checkpoint gene expression between risk groups were also investigated. In the TCGA-TARGET cohort (Figure 6E), the high-risk group had significantly lower mRNA levels of RGMB, CTLA4, CD96, CD274, HAVCR2, LAG3, TIGIT, CD28, and PDCD1LG2 than the low-risk group. Moreover, the seven OS-related ARGs were associated with different TIICs infiltrating (Figure 6F).

Figure 6. Immune infiltration between risk groups. (A) Comparisons between the high- and low-risk groups in terms of stromal score and immune score in the TCGA-TARGET cohort and GSE16091 + GSE21257 cohort. (B) The Sankey diagram of risk score, immune score, and survival status. Violin plot comparing the proportions of TIICs associated with low- and high-risk scores designated by the OS-prediction model for the training cohort (TCGA-TARGET) (C) and the validation cohort (GSE16091 and GSE21257) (D). (E) Differences in expression of immune checkpoint genes between low- and high-risk groups in TCGA-TARGET cohort. (F) Correlation heatmap of prognostic ARGs and TIICs.

OGT knockdown inhibited the invasion and migration capability of osteosarcoma cells

Although O-linked N-acetylglucosamine transferase (OGT) has been identified as an ARG associated with patient prognosis, its functions in osteosarcoma remain unclear. Therefore, OGT’s effects on osteosarcoma cells were investigated in vitro. Compared to normal bone tissues, tumor tissues had higher OGT mRNA expression levels (Figure 7A and Supplementary Table 2). RT-PCR was used to measure OGT mRNA expression in various osteosarcoma cell lines (Figure 7B). The mRNA level of OGT was high in multiple osteosarcoma cell lines. siRNAs (si-1, si-2, and si-nc) were then used to transfect MG63 and U2OS cells. (Figure 7C). Subsequently, transwell assays revealed that OGT knockdown decreased MG63 (Figure 7D) and U2OS (Figure 7E) cell migration and invasion. OGT-related pro-tumor mechanisms were also investigated using KEGG pathway analysis (Figure 7F). In summary, OGT promoted osteosarcoma metastasis and might be a prognostic biomarker.

Figure 7. OGT knockdown inhibits the migration and invasion capability of osteosarcoma cells. (A) OGT expression in normal bone tissues and osteosarcoma tissues (GSE16088 and GSE28425). (B) OGT mRNA expression level in osteosarcoma cells. (C) OGT knockdown in MG63 and U2OS cells was confirmed using RT-PCR and WB. (D, E) Migration and invasion capability of MG63 and U2OS cells were significantly weakened by the downregulation of OGT expression (**p <0.01, ***p <0.001). (F) KEGG pathway analysis of OGT in osteosarcoma.