Files of scRNA-seq data (GSE152048) were accessed via the GEO database. Eleven tumor samples from 11 OS patients were included for analysis. The microarray datasets were based on the 10X Genomics platform. Among them, eight lesions were osteoblastic OS, including six primary, one recurrent, and one lung metastatic lesions; three were chondroblastic OS with each derived from primary, recurrent, and lung metastasis sites.

The scRNA-seq data were processed by using the Seurat R package according to the standard protocol (17). We excluded cells with less than 200 detected genes and genes that were detected in less than three cells and limited the mitochondrial gene proportion to less than 20%. The data normalization was performed by using the LogNormalize method. T-distributed stochastic neighbor embedding (t-SNE), a non-linear dimensionality reduction method, was applied after principal component analysis (PCA) for unsupervised clustering and unbiasedly visualizing cell populations on a two-dimensional map. Marker genes of each cluster were detected using the “FindAllMakers” function, and the criteria for identifying marker genes were set as follows: absolute log2 fold change (FC) > 1, and the minimum cell population fraction in either of the two populations was 0.25. The expression pattern of each marker gene among clusters was visualized by applying the “DotPlot” function in Seurat. Marker-based cell-type annotation was performed by using the SingleR package.

The clusterProfiler R package was used to perform the Gene Ontology (GO) functional enrichment analysis. The marker genes were assigned to various biological processes (BPs), cellular components (CCs), molecular functions (MFs), and pathways. Significant enrichment was set as P < 0.05.

The distribution of NR gene expression levels was evaluated in osteosarcoma tissues by calculating the quantile expression (Q1–Q4) using log10-transformed RNA-seq data. Expression levels were then classified as not expressed (Q1 (0%–25%: -Inf to 0.98) were defined as absent and Q2 (25%–50%: 0.98 to 2.32) as low expression) or expressed (Q3 (50%–75%: 2.32 to 2.93) as moderate and Q4 (75%–100%: 2.93 to 5.13) as high expression). The expression of NR in different types of osteosarcoma tissues was plotted using the “ggpubr” package, and significant difference was analyzed using t-test.

The gene expression profile GSE32981 including a panel of 23 osteosarcoma samples of primary and metastatic origin was downloaded from the GEO database (GPL3307 platform: ABI Human Genome Survey Microarray v2.0). The gene expression profile GSE154540 with CDK4 overexpression in the two-condition experiment (poor responder vs. good responder) was obtained from the GEO database (GPL5477 platform: Agilent-014950 Human Genome CGH Microarray 4x44K). The processed gene expression matrix files were extracted by the R package GEOquery (version 2.54.1). With the expression of NR4A1 in GSE32981 and GSE154540, the R package ggplot2 (version 3.3.5) was used to perform the comparison result in different groups, respectively. Wilcoxon test was used for the test method in the comparisons. P < 0.05 was considered statistically significant.

Survival analysis was used to identify biomarkers from significant feature genes. The log2(fpkm+1) expression data and clinical information of OS in the TARGET database were downloaded from the University of California Santa Cruz (UCSC) Genome Browser database (18). The results were visualized using Kaplan–Meier plots. The samples were divided into two groups (high and low) based on the median expression levels of NR4A1 or NR4A2, followed by the overall survival computed between OS and normal by the K–M survival curve.

In the analysis of GSE152048, we detected 11 main clusters in parallel using unbiased clustering of the cells. Based on the canonical markers and gene profiles in the osteosarcoma tissues, we plotted the t-SNE graph of different cell types ( Figure 1A ) and 11 cell types including “osteoblastic”, “myeloid”, “osteoblastic_proli”, “Osteoclast”, “TIL”, “Chondroblastic”, “Endothelial”, “MSC”, “Pericyte”, “Fibroblast”, and “Myoblast” ( Figure 1A ). In the further analysis, we examined the expression profiles of NRs in the 11 cell types. The distribution of the NRs in different cell types was visualized by the t-SNE plot, and a high expression of NRs was found in cell types including “myoblast,” “Osteoblastic,” “MSC,” “Endothelial,” and “Myeloid” ( Figure 1B ). The heatmap also visualized the expression of individual NR in 11 cell types ( Figure 1C ). Based on the results, we found that several NRs including PPARG, NR1D1, NR2F1, NR2F2, NR4A1, and NR4A2 were markedly expressed in the 11 cell types ( Figure 1C ).

For a further analysis, we plotted the t-SNE graphs based on the metastasis, primary, and recurrent tissues, and the t-SNE graphs showed that 11 cell types were detected in the metastasis, primary, and recurrent tissues ( Figure 2A ). Furthermore, we also analyzed the distribution of NRs in the 11 cell types from metastasis, primary, and recurrent tissues. In the metastasis osteosarcoma tissues, the NRs were mainly detected in cell types including “TIL” and “Osteoblastic”; in the primary osteosarcoma tissues, the NRs were mainly enriched in cell types including “TIL”, “Myoblast,” “Endothelial”, and “Myeloid”; in the recurrent osteosarcoma tissues, the NRs were mainly detected in “Chondroblastic”, “Osteoblast”, and “Pericyte” ( Figure 2B ).

In “Chondroblastic,” no significant difference in the expression of NRs was detected between metastasis and primary osteosarcoma tissues; the expression of NRs in the recurrent osteosarcoma group was significantly higher than that in the primary osteosarcoma group and was lower than that in the metastasis group ( Figure 3 ). In “Endothelial”, no significant difference in the NR expression was detected among the three groups ( Figure 3 ). In “Fibroblast”, the expression of NRs was significantly higher in the metastasis and recurrent groups than that in the primary osteosarcoma tissues, and the expression of NRs in the recurrent group was even higher than that in the metastasis group ( Figure 3 ). In “MSC”, no significant difference in the expression of NRs was detected between metastasis and primary osteosarcoma tissues; the expression of NRs in the recurrent osteosarcoma group was significantly higher than that in the primary and metastasis osteosarcoma groups ( Figure 3 ). In “Myeloid”, the expression of NRs was significantly higher in the metastasis and recurrent groups than that in the primary osteosarcoma tissues, and the expression of NRs in the recurrent group was even higher than that in the metastasis group ( Figure 3 ). In “Myoblast”, no significant difference in the NR expression was detected among the three groups ( Figure 3 ). In “Osteoblastic”, “Osteoblastic_proli”, and “Osteoclast”, the expression of NRs was significantly higher in the metastasis and recurrent groups than that in the primary osteosarcoma tissues, and the expression of NRs in the recurrent group was even higher than that in the metastasis group ( Figure 3 ). In “Pericyte”, the expression of NRs was significantly higher in the primary and recurrent groups than that in the metastasis group, and the expression of NRs in the recurrent group was even higher than that in the primary group ( Figure 3 ). In “TIL”, no significant difference in the expression of NRs was detected between metastasis and primary osteosarcoma tissues; the expression of NRs in the recurrent osteosarcoma group was significantly higher than that in the primary and metastasis osteosarcoma groups ( Figure 3 ).

Furthermore, the individual gene of NRs was further analyzed in different cell types among metastasis, primary, and recurrent osteosarcoma tissues. In “Pericyte”, differentially expressed NRs such as NR4A1, NR3C1, NR2F2, NR2F1, NR1H2, and ROR were detected among metastasis, primary, and recurrent osteosarcoma tissues ( Figure 4 ). In “Fibroblast”, differentially expressed NRs such as NR4A1, NR4A2, NR3C1, NRF2F6, PARA, and THRA were detected among these three groups ( Figure 4 ). In “Myeloid”, differentially expressed NRs such as NR4A2, NR4A1, NR3C1, RXRA, NR1D2, PPARD, and RARA were identified among metastasis, primary, and recurrent osteosarcoma tissues ( Figure 4 ). In “Osteoblastic_proli,” differentially expressed NRs such as NR4A2, NR4A1, NR3C1, NR2F6, NR2F2, NR2F1, RXRA, RORA, and THRA were detected among these three groups ( Figure 4 ). In “Osteoblastic”, differentially expressed NRs such as NR4A2, NR4A1, NR3C1, NR2F6, NR2F2, NR2F1, VDR, NR1H2, RORA, and THRA among metastasis, primary, and recurrent osteosarcoma tissues ( Figure 4 ). In “Osteoclast”, differentially expressed RNAs such as NR4A2, NR4A1, NR3C1, ESRRA, NR2F6, RXRA, NR1H2, RORA, RARA, and THRA were detected among the three groups ( Figure 4 ).

Based on the above heatmap analysis, we found that NR4A1 and NR4A2 were highly differentially expressed among primary, metastasis, and recurrent osteosarcoma tissues. Thus, NR4A1 and NR4A2 were chosen for further survival analysis using the TARGET database. As shown in Figure 5A , a high expression of NR4A1 in the osteosarcoma tissues was significantly correlated with a shorter 5-year overall survival of patients with osteosarcoma ( Figure 5A ). On the other hand, there was no significant association between NR4A2 expression and the 5-year overall survival of patients with osteosarcoma ( Figure 5B ).

In order to further validate the expression of NR4A1 in the osteosarcoma tissues, we further explored two datasets, namely, GSE32981 and GSE154540. As shown in Figure 6A , the expression of NR4A1 was significantly higher in the metastasis osteosarcoma tissues than that in the primary osteosarcoma tissues ( Figure 6A ). In addition, we also examined if NR4A1 was differentially expressed in the osteosarcoma tissues from patients with different chemo-sensitivities. As shown in Figure 6B , the expression of NR4A1 was significantly higher in osteosarcoma tissues from patients with poor chemosensitivity than that from patients with good chemosensitivity ( Figure 6B ).

The percentage of osteoblasts with a high NR4A1 expression and low NR4A1 expression in the primary, metastasis, and recurrent osteosarcoma tissues was further analyzed in GSE152048. As shown in Figure 7A , the percentage of osteoblasts with a high NR4A1 expression was higher in the recurrent osteosarcoma tissues than that with a low NR4A1 expression ( Figure 7A ). The volcano plot illustrates the DEGs in the osteoblasts between the high NR4A1 expression group and low NR4A1 expression group ( Figure 7B ). The GO enrichment analysis showed that the up-regulated DEGs were mainly enriched in GO terms such as “response to unfold protein,” “response to topologically incorrect protein,” and “fat cell differentiation” ( Figure 7C ), while the down-regulated DEGs were mainly enriched in GO terms such as “extracellular structure organization,” “extracellular matrix organization,” and “ossification” ( Figure 7D ).