Most osteosarcomas develop from mesenchymal stem cells with osteogenic potential and are highly aggressive and distantly metastatic.[7]. The incidence in the general population is 2–3/million/year, but is higher in adolescents, with a peak annual incidence of 8–11/million/year at the age of 15–19 years[8]. It is histologically subdivided into conventional, low-grade central, periosteal, parosteal, capillary dilated, chondroblast and small cell types[9]. Current treatment of osteosarcoma is mostly based on chemotherapy regimens with adriamycin or surgery. Before multidrug chemotherapy, 90% of patients with osteosarcoma died from pulmonary metastases[10]. Data from next-generation sequencing studies have gradually revealed new patterns of genomic alterations in osteosarcoma, and may even identify new potential therapeutic targets. There is no doubt that epigenetic modifications also play an invaluable role in the development of osteosarcoma, with hypermethylation and hypomethylation observed in tissues. Genes that have been identified to have a role include TP53, RB1 and members of the RecQ DNA decapping enzyme family (BOX 1), which are mostly associated with drug resistance and metastasis in osteosarcoma. However, to date, research on common molecular therapeutic targets in osteosarcoma remains disappointing[11].
While current treatment for osteosarcoma is based on chemotherapy and surgery, a growing number of potential gene-targeted therapies are in development[12]. Therefore, it is crucial to observe and study the mechanisms of osteosarcoma development at the molecular level. Based on this, differentially expressed genes (DEGs) have been effectively used to predict OS phenotypes such as invasion, metastasis, angiogenesis, and drug resistance, which help us to be able to intervene in target genes from multiple aspects, and thus can play a role in inhibiting the development of osteosarcoma or improving the survival rate.
In this study, we performed a specific collation of gene expression profiles from the GEO database for GSE14359 and used R (version 4.1.1) to perform a specific code to analyze this dataset. We identified a total of 625 DEGs using the limma package, including 314 up-regulated genes and 311 down-regulated genes. Given that the number of screened DEGs was too large to facilitate our subsequent merit-based analysis, we immediately performed a WGCNA-weighted gene co-expression network analysis to identify 32 gene modules and selected the MEgrey60 module containing 594 genes based on absolute values as a way to find genes with the higher association. the DEGs and WGCNA-identified modules were taken to intersect to obtain candidate differential genes A total of 248 were obtained. A PPI network of differential candidate genes encoding proteins were constructed from the STRING database (minimum required interaction score = 0.990) and the 50 most important associated genes were screened. Among these genes, EGFR, FLT1, VAV2, and PXN had the highest modality. From this, the core gene network map was obtained by CYTOSCAPE software based on the PPI network interactions and the corresponding node files. Finally, the core ranking histogram of key genes was obtained from the above files, and the top 23 genes as shown in the figure, in order, were EGFR, FLT1, PXN, VAV2, MMP3, SPP1, ADAMTS2, AOX1, CHI3L1, CNDP2, COX7B, CTSB, CTSL, DDX18, EFNA1, ELL ENTPD1, FGF7, FGFR3, FHL2, GATM, HEY1, HEY2. GO enrichment analysis of candidate differential genes in OS using R software and correlation analysis confirmed that the candidate differential genes are mainly involved in cellular processes such as osteogenesis, epithelial cell proliferation, osteoblast differentiation and molecular functions such as growth factors and fibronectin binding. This finding is also consistent with biological processes that play a key role in OS development and progression. Li et al. found that the degree of ossification within soft tissue masses in common type osteosarcoma correlated with tumor surgical stage, chemotherapy sensitivity, and prognosis, and that regulation of PLK2 enriched for TAp73 affects osteogenic differentiation and prognosis in human osteosarcoma, also confirming the role of the ossification process in osteosarcoma[13]. Our data suggest that this biological process involves the largest number of genes. Luo et al. found that CDKN2B-AS1 can play an oncogenic role in osteosarcoma by promoting cell proliferation and epithelial-to-mesenchymal transition, suggesting that epithelial cell proliferation is an important part of osteosarcoma development[14]. Our enrichment analysis also confirmed exactly this, and Baumann's team demonstrated that collagen can be post-translationally modified by prolyl and lysyl hydroxylation followed by glycosylation of hydroxylysine, and the experimental conclusion illustrates that complete loss of collagen glycosylation decreases osteosarcoma cell proliferation and viability[15]. In other words, the invasion of bone tissue by tumor cells affects the dynamic balance between bone resorption and bone formation[16]. The bone remodelling process is significant in osteosarcoma development. The role of fibronectin binding in osteosarcoma inhibition is also significant, and some studies have also confirmed that it enhances the apoptotic response of U2-OS and contributes to the fight against osteosarcoma[17]. It is well known that vascular endothelial growth factor A (VEGF) is one of the most important growth factors that regulate vascular development and angiogenesis. VEGF is involved in different stages of bone repair, including the inflammatory phase, endochondral ossification, healing tissue formation and intramembranous ossification during bone remodelling.[18] These biological processes are all involved in our work. These biological processes were confirmed in our analysis, and it is reasonable to speculate that phenotypic studies on osteosarcoma could focus on these aspects.
In addition, candidate differential genes were enriched in the Kyoto Encyclopedia of Genes (KEGG) pathways mainly in adherent spots, Rap1, PI3K-Akt, MAPK, and insulin resistance, and only a few pathways with relatively more enriched genes are listed in this paper for subsequent studies. It has been shown that focal adhesion kinase (FAK) overexpression and phosphorylation may predict more aggressive biological behavior in osteosarcoma and may be an independent predictor of poor prognosis[19] Feng et al. also demonstrated the involvement of FAK in the migration of human osteosarcoma cells[20]. We know that cancer cachexia is easily driven by inflammation, and metabolic alterations (e.g. increased energy expenditure, elevated plasma glucose, insulin resistance and excess catabolism). Insulin insensitivity decreases glucose uptake in organs and leads to loss of skeletal muscle and adipose tissue.[21] Imerb also found in his study of insulin resistance that it leads to impaired skeletal regulation and imbalances in bone homeostasis and pathological bone.[22]. Insulin resistance, a common metabolic feature and a risk factor for many diseases are caused by the inability of insulin to perform its normal metabolic role, as well as by nutritional imbalances and abnormal lipid accumulation in skeletal muscle, liver and adipose tissue.[23] The PI3K/Akt pathway is considered to be one of the most important oncogenic pathways in human cancers. There is increasing evidence that this pathway is frequently over-activated in OS and contributes to disease onset and progression including tumorigenesis, proliferation, invasion, etc.[24]. Certainly, the role of MAPK signaling pathway in the development of osteosarcoma has also been long established[25–27] The role of MAPK signaling pathway in the development of osteosarcoma has also been well established. In addition, the expression levels of HIF-1, Rap1, PI3K and Akt were also found to be elevated in OS cells, suggesting that the Rap1/PI3K-Akt pathway could be a therapeutic target for OS.[28]. In conclusion, our KEGG enrichment results were consistent with previous studies, and most genes were enriched in two pathways, PI3K/Akt and MAPK, suggesting that the future direction of research on osteosarcoma may have to focus on these two pathways.
We constructed a PPI network of candidate differential genes encoding proteins, and then based on the core ranking histogram of key genes, we selected the top-ranked genes such as EGFR, FLT1, VAV2, PXN, and FGF7 for analysis. As shown by the survival curves, the survival curves of the key genes did not show statistical significance, except for FGF7. As seen by the FGF7 survival curve, the FGF7 high expression group had a higher survival rate than the FGF7 low expression group until 100 months of survival. The rest of the Hub genes failed to show differences. Thus, we speculate that FGF7 may serve as a marker to predict survival outcomes in OS patients, which can be focused on in later validation experiments. However, it has also been shown that FGF7 induces the proliferation of osteosarcoma cells and accelerates the secretion of EMT and inflammatory mediators[29]. This is contrary to our conclusion, so the follow-up study should make FGF7 a focus and increase the sample size appropriately for validation again. As for Epidermal growth factor receptor (EGFR), it can activate various signaling pathways such as PI3K/Akt, Jak/STAT, etc. Targeted regulation of EGFR expression can inhibit the development of osteosarcoma.[30]. Therefore, combined with survival analysis, we believe that although EGFR cannot predict the survival of patients, it can be a potential therapeutic target. Tsuchida, while studying the molecular mechanism of tumor progression after chemotherapy, demonstrated that cisplatin (CDDP) upregulates FLT1 expression, leading to the survival and proliferation of highly oncogenic side population (SP) cells in cell lines such as osteosarcoma[31] Yin et al. slowed the growth of osteosarcoma tumors in this mouse model by FLT1 gene modification.[32]. Therefore, FLT1 could also be an intervention target for osteosarcoma. It has been shown that VAV2 is associated with the migration and invasion of sarcoma cells[33]. PXN, a focal adhesion-associated protein, significantly promotes signaling between protein networks[34]. There is evidence that PXN may be associated with cell proliferation and migration[35]. However, PXN is currently poorly studied in osteosarcoma, suggesting that we should subsequently deepen our research on PXN as a therapeutic target in osteosarcoma. Also based on the analysis of gene expression and clinical correlation, we should also pay sufficient attention to the gender and age of tumor tissue origin. The comparison of promoter methylation level analysis suggests that we should also take TP53 mutation typing into full consideration.
Finally, we performed immune cell infiltration analysis of the tissues and found that infiltration of M2 macrophages, unactivated dendritic cells, and neutrophils accounted for a greater proportion of OS tissues, while infiltration of naive B cells and unactivated CD4 + memory T cells was statistically less. This suggests that we should focus on M2 macrophages in osteosarcoma research to inhibit tumor cell development by modulating immune infiltrating cells in the osteosarcoma microenvironment. In the correlation analysis of core genes with immune and infiltrating cells, key genes were found to be generally more correlated with naive B cells. On the correlation analysis between core genes and immune checkpoints, three immune checkpoints, ADAMTS2, AOX1, and CHI3L1, were found to be the most correlated with core genes. These two correlation analyses also provided us with the next direction in osteosarcoma research.
The identification of differential genes and further bioinformatics analysis have been carried out in the past few years, and these findings may provide new ideas for the study of the mechanisms of OS onset, development, metastasis, invasion, and drug resistance. However, a limitation of this study is that it was not validated with clinical specimens, which may need to be added in future studies.