The GEO (http://www.ncbi.nlm.nih.gov/geo) database is a free public gene-expression data repository containing microarray and high-throughput sequencing data. Gene-expression datasets from GSE42352 (12, 13), including 103 osteosarcoma cases and 15 normal controls, were collected from GEO as the test cohort. Using the median value of GNG12 expression as the cutoff, the 103 patients were divided into high (n = 52) and low (n = 51) GNG12 mRNA expression groups. Significant differentially expressed genes (DEGs) between high and low GNG12-expression groups were screened with the limma package (http://www.bioconductor.org/packages/release/bioc/html/limma.html) in R version 3.6.3 (http://www.R-project.org/) (14). Genes were considered differentially expressed with an adjusted p < 0.5 and |log fold change| (|logFC|) ≥ 1. Heat maps were constructed using R package “pheatmap”, version 1.0.12 (https://cran.rproject.org/web/packages/pheatmap/index.html), and volcano plots were generated with “ggplot2” version 3.3.3 (https://cran.r-project.org/web/packages/ggplot2/index.html). In addition, gene expression profiles and clinical data of osteosarcoma cases were collected from TARGET (https://ocg.cancer.gov/programs/target). The TARGET cohort (n = 101) were divided into high/low groups using the same cutoff as the GEO cohort (high, 51 cases; low, 50 cases). The same method as described above was used to determine DEGs between TARGET high/low groups. Common DEGs between the GSE42352 and TARGET datasets were represented with a Venn diagram made in ggplot2.

The Second Affiliated Hospital of Nanchang University provided 78 osteosarcoma specimens from January 2012 to December 2016 as an external validation cohort. Samples were provided as formalin fixed paraffin-embedded blocks.

This research was approved by the Ethics Committee of The Second Affiliated Hospital of Nanchang University [Review (2020) No. (086)]. All participants provided informed consent.

Functional enrichment analyses (Gene Ontology, GO; https://www.geneontology.org and Kyoto Encyclopedia of Genes and Genomes, KEGG; https://www.genome.jp/kegg/) were performed on DEGs between high and low GNG12-expression groups in the test cohort, using clusterProfiler version 3.14.3 (http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html). The results were visualized with the “ggplot2” package (https://cran.r-project.org/web/packages/ggplot2/index.html). Xiantao Academic online (https://www.xiantao.love) was used to generate a diagram of the functional enrichment network analysis.

Significant differences in function and pathways between high/low GNG12-expression groups were determined using GSEA in clusterProfiler and Xiantao Academic online. For each analysis, gene set permutations were performed 1,000 times to obtain a normalized enrichment score (NES). Enrichment was considered significant with adjusted p < 0.05, false discovery rate (FDR) q < 0.25, and |NES| > 1. The selected reference gene set was c2.cp.v7.2.symbols.gmt (Curated), and results were visualized in ggplot2.

A protein–protein interaction (PPI) network of DEGs was constructed using Metascape online tools (https://metascape.org) with the following parameters: min network size = 3 and max network size = 500. Crucial proteins in this PPI network were screened using the Molecular Complex Detection (MCODE, http://apps.cytoscape.org/apps/mcode), a plug-in of Cytoscape version 3.7.2 (https://cytoscape.org/).

CIBERSORT (https://cibersort.stanford.edu/) is an algorithm developed to characterize the cellular composition of complex tissues through gene expression profiles compared with a signature matrix (LM22). LM22 comprises 547 genes that define 22 immunization cell subtypes. This analysis was used to calculate infiltration abundance of 22 immune cells (T cells, B cells, plasma cells, natural killer cells, myeloid subgroups) in 103 osteosarcoma samples across high and low GNG12-expression groups. We also analyzed any correlations between the two GNG12 groups in terms of expression distribution among the 22 infiltrating immune cells; the correlation heat map was then plotted. Red represents a positive correlation, and blue represents a negative correlation, with darker colors and values closer to 1 indicating stronger correlations.

The UCSC XENA platform (https://xenabrowser.net/datapages/) is a pan-cancer HTSeq-TPM database of TCGA and GTEx, processed uniformly by Toil. Based on XENA data, differential GNG12 mRNA expression between tumor and normal tissues was determined in R and visualized with ggplot2. Spearman’s correlations were calculated between GNG12 expression in the TARGET cohort with the expression of nine DEGs (GREM1, CAV1, PTGDS, PRB2, PRB1, NDUFB9, AMBN, FAT3, and IBSP) common across GSE42352 and TARGET datasets. Again, the analysis was performed in R and the results visualized in ggplot2.

Kaplan–Meier (KM) analyses of the TARGET database were performed using R packages survival (https://www.rdocumentation.org/packages/survival/versions/2.42-3) and survminer (https://cran.rstudio.com/web/packages/survminer/index.html), to determine whether GNG12 expression accurately predicted overall survival time (OS) and progression-free survival (PFS). The timeROC package (https://cran.r-project.org/web/packages/timeROC/index.html) was used to analyze 1-, 3-, and 5-year survival prognosis of GNG12-related lncRNAs, miRNAs, and genes. All outcomes were visualized in ggplot2.

An ROC curve was generated using the pROC package (https://cran.r-project.org/web/packages/pROC/). The ordinate was drawn with the true positive rate (sensitivity), and the abscissa was the false positive rate (1 − specificity). This analysis determined whether GNG12 expression can distinguish between 15 standard samples and 103 osteosarcoma samples. Also determined was the cutoff value that produces the highest likelihood ratio for assessing the recognition threshold of GNG12 as a biomarker for distinguishing osteosarcoma from normal tissue.

To further verify the suitability of GNG12 expression for predicting survival, immunohistochemistry was performed on paraffin sections following the standard protocol (GNG12, Ab204757, 1:100). All slides were observed and photographed under a XSP-C204 microscope (CIC). GNG12 expression was assessed using the H-score.

Known survival outcomes of osteosarcoma patients with high and low H-scores were subjected to KM analysis using the R package “survival” and visualized in “surviminer”. The accuracy of GNG12 expression in predicting survival prognosis (1, 3, and 5 years) was verified with timeROC and visualized in ggplot2. An AUC >0.5 directly reflects the prognostic value of the biomarker. Additionally, an AUC closer to 1 and to the (0, 1) point indicates greater authenticity of the prognostic prediction.

Statistics were performed in Microsoft Excel and R version 3.6.3. The log-rank test was used to perform the KM survival analysis. Continuous variables that violated the normality assumption according to the Shapiro–Wilk normality test were analyzed with the Wilcoxon rank sum test. Categorical variables were analyzed with chi-square tests. Significance was set at p < 0.05 (ns, p > 0.05; *p < 0.05; **p < 0.01; ***p < 0.001).

GNG12 expression was significantly different between the standard samples (n = 15) and osteosarcoma samples (n = 103) ( Figure 1A ). After screening for DEGs between high and low GNG12-expression groups ( Figure 1B ), we found 210 coexpressed genes (165 upregulated and 45 downregulated, Supplementary Table 1 ) that we have visualized in a volcano plot ( Figure 1D ) and a heat map ( Figure 1E ). Using the same method to divide the TARGET cohort (n = 101) into high/low GNG12-expression groups ( Figure 1C ), we obtained 347 DEGs (178 upregulated and 169 downregulated, Supplementary Table 2 and Figures 1F, G ). We then screened out nine common DEGs (six upregulated and three downregulated) across the two DEG sets ( Figures 1H, I ).

The GO and KEGG enrichment results for the GEO cohort ( Figures 2A, B ) indicated that the following terms were enriched: “protein targeting to endoplasmic reticulum (ER)”, “establishment of protein localization to ER”, “focal adhesion”, and “cell adhesion molecule binding” ( Figure 2C ). For the TARGET cohort, enriched terms were “ossification”, “cartilage development”, “integrin binding”, “PI3K-Akt signaling pathway”, and “Focal adhesion” ( Figure 3 ).

We selected msigdb.v7.0.entrez.gmt as the reference gene set for our GSEA. For the GEO cohort, DEGs were involved in “extracellular matrix organization” from the Reactome database, “core matrisome” from NABA, “cytoplasmic ribosomal proteins” from WP, and “cell adhesion molecules cams” from KEGG ( Figure 4A ). For the TARGET cohort, DEGs were involved in “class A 1 rhodopsin like receptors”, “GPCR ligand binding”, and “diseases of metabolism” from the Reactome database, and “Matrisome” from NABA ( Figure 4B ).

Based on the Metascape online tool, we constructed two PPI networks of DEGs from the GEO and TARGET cohorts ( Figures 5A, C ). We then screened hub-gene clusters using the MCODE clustering algorithm. In the GEO cohort, screened hub genes were linked in the following networks: MCODE_1 (cytoplasmic translation, cytosolic ribosome, and ribosome; EIF3A, RPL23, RPL14, RPL7, RPL27A, RPS28, RPLP1, RPS3A, RPL10A, and RPN2), MCODE_2 (NLS-bearing protein import into the nucleus, protein import into the nucleus, and import into the nucleus; KPNA2, RGPD6, RGPD3, RGPD8, SUMO2, ACTR3, WIPF1, ARPC5, and GMFG), and MCODE_3 (ER, growth factor activity, and extracellular matrix structural constituent; VCAN, CCN1, AMBN, CDH2, VGF, and IL6; Figure 5B ). Hub genes of the TARGET cohort were involved in more networks. These were MCODE_1 (extracellular matrix structural constituent conferring tensile strength, collagen trimer, and protein digestion and absorption; COL8A2, COL9A1, COL11A1, COL11A2, COL15A1, COL22A1, COL24A1, P4HA3, COL14A1, COL27A1, COL2A1, COL3A1, COL6A1, and COL6A2), MCODE_2 (in utero embryonic development, peptidyl-tyrosine phosphorylation, and peptidyl-tyrosine modification; FGFR4, ISLR, FGF2, IGF2, TGFB3, CAV1, TP53, ACTA2, ALDH1A1, PYGM, GPI, GFPT2, and H2AW), MCODE_3 (G-protein beta-subunit binding, neuropeptide hormone activity, and heterotrimeric G-protein complex; GNG12, GAL, CORT, NPB, GNG4, GNGT2, HCRT, HRH1, PTGFR, ADRA1D, and TBXA2R), MCODE_4 (ER lumen, biomineral tissue development, and biomineralization; MELTF, PENK, MEPE, SPP2, MOTUM, DMP1, FBN1, AMBN, AMELX, and GPC3), MCODE_5 (anchored component of membrane; NTM, ALPL, NTNG1, CD109, and LY6K), MCODE_6 (lamellar body, mutivesicular body and respiratory gaseous exchange by respiratory system; SFTPB, SFTPC, SFTPA1, and SFTPA2), and MCODE_7 (Wnt signaling pathway; CCND2, TCF7L2, and FOSL1; Figure 5D ).

For further immune infiltration analysis, we used CIBERSORT package and LM22 algorithm to calculate the infiltration abundance of 22 kinds of immune cells between high/low expression level of GNG12 groups in test and validation cohort, including different B cells, T cells, natural killer cells, plasma cells, and different myeloid subsets. The results were visualized as two violin maps ( Figures 6A, B ) and two correlation heat maps ( Figures 6C, D ) to reveal the differences in the expression of immune cells between groups.

Evaluation of differential GNG12 expression in pan-cancer revealed that GNG12 was downregulated in 12 tumor types (i.e., ACC, BLCA, BRCA, CESC, KIRC, LAML, OV, PCPG, PRAD, TGCT, UCEC, and UCS) and upregulated in 11 (CHOL, COAD, DLBC, ESCA, GBM, LGG, LIHC, PAAD, READ, STAD, and THYM; Figure 7A ). When we used a gene-expression heat map to examine the nine common DEGs across GEO and TARGET cohorts, we found that GREM1 and CAV1 expression were positively correlated with GNG12 expression in the TARGET cohort. Additionally, NDUFB9, AMBN, FAT3, and IBSP expression were negatively correlated ( Figure 7B ). We then generated six scatterplots of correlations between the expression of these six genes and of GNG12 ( Figures 7C–H ).

Using limma, we obtained 22 differentially expressed lncRNAs (DELs) from the TARGET database and visualized them on a heat map ( Figure 8A ). We then used the miRcode website to predict miRNAs that were highly conserved with DELs. Furthermore, we identified mRNAs targeted by the miRNAs through a series of databases and then constructed a ceRNA network with DEGs ( Figure 8B ).

The results of our KM analysis indicated that OS and PFS differed significantly between low and high GNG-expression groups ( Figures 8C, F ). In addition, the ceRNA network revealed that low HOTTIP and GNG12 expression led to high risk. Elevated hsa-miR-27a-3p expression also resulted in high risk ( Figures 8D, E ). Through time-dependent ROC curves, we found that GNG12 was a reliable predictor of patient prognosis ( Figure 8G ).

We used an ROC curve to analyze whether GNG12 expression could distinguish between 15 normal samples and 103 tumor samples from the GEO dataset. We found an AUC of 0.920, indicating that GNG12 expression was an excellent biomarker for distinguishing osteosarcoma from normal tissue. ( Figure 8H ).

We used 78 patients with osteosarcoma as a validation cohort to confirm the prognostic reliability of GNG12 ( Table 1 and Supplementary Table S3 ). GNG12 expression levels (high vs. low) of these patients were confirmed using immuno-histochemistry (IHC) ( Figures 9A, B ). The results of KM analysis showed that patients with low GNG12 expression tend to have poor survival prognosis ( Figure 9C ). Analysis of time-dependent ROC curves revealed that 1-/3-/5-year AUCs were 0.961, 0.826, and 0.808 ( Figure 9F ). In addition, GNG12 expression was lower in metastatic than in non-metastatic osteosarcoma ( Figure 9D ); patients with metastatic osteosarcoma have a poor survival prognosis ( Figure 9E ).