Screening of disorders associated with osteosarcoma by integrated network analysis

Osteosarcoma is a common primary malignant bone tumor. It is also the most common primary bone tumor in children and adolescents under 20 years old. It is also the third most common cancer in childhood and adolescence [1–3]. As a malignant tumor, osteosarcoma mainly affects long bones, but may also involve other bones in the body. It has a bimodal distribution and will reach its peak in the second decade of life and adulthood [4]. On the one hand, due to the lack of lymphatic system in the skeleton, metastatic spread in bone and sarcoma is completely hematogenous. The most common sites are lung and bone [5]. On the other hand, osteosarcoma is considered to be a highly vascularized bone tumor, which mainly enters the lung through early metastatic diffusion of intratumoral blood vessels [6]. Osteosarcoma belongs to primary bone tumors. Its recurrence and metastasis potential may be due to cell subsets with stem cell-like characteristics, which maintain the regeneration ability of the whole tumor. Osteosarcoma is associated with the most common genetic abnormalities in human cancers [7]. Contrary to other malignant tumors, many genetic and environmental factors can lead to the development of osteosarcoma. Therefore, osteosarcoma is defined phenotypically rather than molecularly [8,9]. Osteosarcoma is a highly hereditary and unstable tumor with a high incidence. Because of local recurrence, metastasis, and chemotherapy resistance, it has a poor prognosis [10–12]. Despite significant improvements in treatment strategies in recent years, for most patients with metastatic or recurrent osteosarcoma, treatment outcomes are still poor [13]. Studies have shown that osteosarcoma is an invasive cancer in the skeletal system, in which altered oncogenes and tumor suppressor genes can participate in tumor cell migration, angiogenesis, cell apoptosis, and proliferation [14].

Although osteosarcoma is a rare malignant tumor, it is listed as one of the main causes of cancer-related death in the pediatric age group [15]. Therefore, in recent years, many biologists and pharmacists have carried out a series of in-depth studies on the effective treatment of osteosarcoma. The data show that the optimal treatment for osteosarcoma patients is neoadjuvant chemotherapy followed by complete surgical resection and adjuvant chemotherapy [2]. The treatment of osteosarcoma now includes systemic chemotherapy and local control surgery, multidrug chemotherapy, and active surgical techniques, which to some extent improves the survival rate of patients [16,17]. Hornicek et al. have shown that EZH2 is essential for the growth and metastasis of osteosarcoma, and epigenetic therapy targetting EZH2 through specific inhibitors pharmacology can constitute a new method for the treatment of osteosarcoma [12]. It is reported that miRNA-27a can promote the proliferation, migration, and invasion of human osteosarcoma cells. Therefore, miRNA-27a can be used as a non-invasive diagnostic and prognostic biomarker for osteosarcoma patients [18]. In addition, the down-regulation of miRNA-152 can also be considered as a predictor of diagnosis and prognosis in osteosarcoma patients [19]. As a tumor suppressor gene, inactivation of NPRL2 contributes to the development of tumors. The expression of NPRL2 is negatively correlated with the survival rate of osteosarcoma patients, and it has important value as a prognostic factor of osteosarcoma [20]. At present, the best treatment of osteosarcoma includes multidrug chemotherapy and active surgical removal of all the affected parts of the disease. This conclusion has also been confirmed by numerous experimental results [21].

Here, based on the gene-phenotype correlation network, we propose a comprehensive analysis to screen the disorders associated with osteosarcoma typing. We have predicted potential disorder molecules associated with disease, which not only provides a new insight into the typing and treatment of osteosarcoma but also provides rich resources and guidance for biologists to further design experiments

We collected the expression microarray dataset of osteosarcoma (GSE94805) from NCBI Gene Expression Omnibus database [22]. Then, differentially expressed genes (DEGs) of normal versus stationary osteosarcoma and normal versus senile osteosarcoma were identified with R language limma package [23] (P<0.05).

First, DEGs shared by the two DEG sets were screened and then coexpression network analysis was performed on these genes expression profiles with weighted gene coexpression network analysis [24] to explore coexpression modules in osteosarcoma. The weighted value of correlation coefficient is used to calculate the correlation coefficient between any two genes (Pearson coefficient) by taking the N power of the correlation coefficient. The connection between genes in the network obeys scale-free networks, which makes the algorithm more biologically meaningful. Finally, a hierarchical clustering tree is constructed by the correlation coefficient between genes. Different branches of the clustering tree represent different gene modules, and different colors represent different modules.

Exploring the role of gene is often an effective means to study the molecular mechanism of disease, and the function and pathway module genes involving in are helpful for characterizing the effect of gene module on the occurrence and development of disease. Therefore, we used R language clusterProfiler package [25] to carry out GO function (P value cutoff = 0.01; qvalue cutoff = 0.01) and KEGG pathway (P value cutoff = 0.05; qvalue cutoff = 0.2) enrichment analysis, respectively.

First, all human transcription factor (TF)-target data and human ncRNA-protein data (score > 0.5) were downloaded from TRRUST V2 database [26] and RAID 2.0 database [27], respectively. Furthermore, pivot analysis based on these interaction data was performed to predict the regulatory relationships between TFs, ncRNAs and modules. Pivot analysis refers to searching regulators with at least two interactors in module, and the number of interactors was verified to be significant using the hypergeometric test (P value < 0.01). Such regulators were thought be key regulators significantly regulating modules.

All blood samples were confirmed by experienced pathologists and informed consents were obtained from all patients. Human tissue samples were collected according to the International Ethical Guidelines for Biomedical Research involving human and subjects. This research was approved by the Orthopaedic Department of Binzhou Medical University Hospital and carried out in line with the regulations of the Binzhou Medical University Hospital.

Specifically, total RNA in the blood was extracted and transcribed into cDNA with a reverse transcription kit and qPCR reaction was conducted with the SYBR qPCR Detection Kit. The qPCR program began the initial 3-min denaturation step at 95°C to stimulate the hot-start iTaqTM DNA polymerase, followed by 45 cycles of denaturation at 95°C for 10 s, and annealing and extension at 60°C for 45 s. The internal reference genes were β-actin and U6.

To screen out potential dysregulated molecules that are closely related to the occurrence and development of osteosarcoma, we identified DEGs of normal versus stationary osteosarcoma (5669 DEGs) and normal versus senile osteosarcoma (8346 DEGs) based on the expression microarray dataset of osteosarcoma. Finally, a total of 3767 DEGs shared by the two DEG sets were obtained for further analysis (Supplementary Table S1).

To further investigate the role of DEGs in osteosarcoma, we first performed coexpression analysis based on the expression of 3767 DEGs. A total of 20 coexpression modules were excavated as dysfunction modules of osteosarcoma, involving 3757 DEGs (Figure 1). Genes can influence the occurrence and development of diseases through their own molecular functions (MF) and involved biological processes (BPs). Module as a set of genes will play more significant role in the pathogenesis of disease. Further, the GO function and KEGG pathway enrichment were performed on all module genes. The analysis identified 7433 MF entries, 5274 cell component entries, 4267 BP entries, and 2477 KEGG pathway entries (Figure 2, Supplementary Table S2). In addition, we analyzed network connectivity based on dysfunction module and identified 38 key endogenous genes, including ARF1, HSP90AB1, and TUBA1B.

Transcriptional regulation of genes is key regulation in the occurrence and development of diseases, and TFs are an important factor in this process. In the present study, we performed pivot analysis to predict key TFs significantly regulating modules based on TF-gene relationships. The results identified 53 TFs and 63 TF-module interaction pairs (Figure 3, Supplementary Table S3). Statistical analysis showed that E2F1 had significant regulatory relationships with the three dysfunction modules, both NFKB1 and EGR1 regulate the two dysfunction modules.

Furthermore, ncRNAs are important regulators of post-transcriptional regulation. Similarly, pivot analysis was performed for key ncRNAs based on relationships between ncRNAs and genes. We have obtained 858 ncRNAs that have significant regulatory effects on modules, involving 1160 ncRNA-module interaction pairs (Figure 4, Supplementary Table S4). Statistical analysis of the results showed that MALAT1 had a regulatory effect on the seven dysfunction modules, miR-590-3p and TUG1 interacted with six dysfunctional modules. These key TFs/ncRNAs may participate in the process of osteosarcoma through regulating disease-related dysfunction modules. Furthermore, the expression level of key genes was verified by qPCR (Figure 5). We found that the expression trend of key genes was consistent with the previous results.

Osteosarcoma is the most common primary bone tumor in children and adolescents, and the prognosis of advanced osteosarcoma patients with metastatic evidence is poor [28]. Pathological interpretation of malignant bone tumors is one of the most challenging areas in surgical pathology, so even if primary osteosarcoma is uncommon, it also shows significant morphological heterogeneity and has a wide range of biological characteristics [29]. Although in recent years biomedical scientists have conducted in-depth studies on the pathogenesis and treatment mechanism of osteosarcoma, there is still a lack of systematic and in-depth exploration of key disorders in the classification of osteosarcoma. In order to further understand the pathogenesis of osteosarcoma, we have integrated abundant resources to conduct a variety of analysis of osteosarcoma, such as gene differential analysis, gene coexpression analysis, transcriptional, and post-transcriptional regulation analysis.

In order to explore the mechanism of osteosarcoma potential pathogenic genes, we first analyzed the genes expression profiles of osteosarcoma, obtained 3767 DEGs, and clustered them into 20 coexpression modules. Subsequently, in view of the results of enrichment analysis, we found that genes in the 20 modules obtained from osteosarcoma mainly participate in nuclear chromosome segregation, microtubule organizing center part and G2/M transition of mitotic cell cycle, and other functional pathways. In addition, we identified 53 TFs that significantly regulating these 20 modules, involving 63 TF-module interaction pairs. It was found that E2F1 had a significant regulatory relationship with the three dysfunction modules and was considered as a key regulation molecule in the process of osteosarcoma. This conclusion has also been confirmed in Zhang and Ding’s research. E2F1 is identified as a new regulator of osteogenic differentiation of osteosarcoma cells induced by ATRA. Meanwhile, E2F1 level has also been found to be associated with poor prognosis of osteosarcoma patients [30,31]. On the other hand, both NFKB1 and EGR1 have been found to regulate the two dysfunction modules, which may also be key regulators involved in the pathogenesis of osteosarcoma. The NF-κB TF, which activates the transcription of genes that regulate a variety of fundamental BPs, including immune response, cell survival, and development [32,33]. NF-κB is significantly involved in the development of breast cancer and other cancers [34]. Nuclear factor-κB1 (NFκB1) gene plays an important role in the pathogenesis of osteosarcoma, which indicates that osteosarcoma is closely related to the gene polymorphism of NFKB1 [35]. According to Matsunoshita et al., EGR1 can be used as a tumor promoter or inhibitor. Experiments show that the strong expression of EGR1 can prevent osteosarcoma cells from migrating into blood vessels, thus inhibiting osteosarcoma to a certain extent [36].

In addition, ncRNA has been considered to be an important regulator of disease occurrence and development, especially inflammatory-related diseases [37]. We conduct pivot analysis based on the targetting relationship between ncRNA and genes. The predicted results show that 858 ncRNAs have significant regulatory effects on modules, including 1160 ncRNA-module interaction pairs. On the one hand, MALAT1 plays an important role in regulating seven dysfunction modules. A number of studies have shown that MALAT1 has carcinogenic effect in osteosarcoma tumorigenesis and metastasis, so it can be considered that MALAT1 may be a promising therapeutic target for osteosarcoma patients [38,39]. On the other hand, miRNAs-590-3p and TUG1 also play an important role in the regulation of six dysfunctional modules, and also play an important role in module dysfunction. MiR-590-3p is an anticancer miRNA, which can inhibit the proliferation and metastasis of osteosarcoma cells [40]. TUG1 knockdown can also inhibit the proliferation, migration, and invasion of osteosarcoma cells and promote cell apoptosis. It can be found that TUG1 knockdown plays an important role in the development of osteosarcoma [41]. At the same time, we also screened a series of genes with the greatest connectivity as the internal gene of the module. A total of 38 endogenous genes were obtained, which may represent potential key regulators of osteosarcoma typing. Statistical analysis of these endogenous genes showed that some of them could be identified as having a regulatory role in the pathogenesis of osteosarcoma. Amongst them, the discovery of ARF1 is of great significance for improving the diagnosis of osteosarcoma [42].

A series of regulatory factors predicted in the present study have certain regulatory effects on the potential pathogenesis of osteosarcoma in varying degrees. However, in addition to the above key factors, other unmentioned ncRNA and transcription factors may play a role in the pathogenesis of osteosarcoma, which need further exploration. Overall, our study is based on the results of phenotypic correlation network analysis, and screening guides the identification of key disorders in osteosarcoma typing. This will not only provide a new way for biologists and pharmacists to study the pathogenesis of osteosarcoma but also provide valuable reference for their follow-up treatment.

We are grateful to the institution Binzhou Medical University Hospital who supported our research.

Y.D. wrote the manuscript and performed part of this research. K.Z. conceived and designed the study. Z.S., X.G., and Q.F. performed the majority of experiments.

The authors declare that there are no sources of funding to be acknowledged.

The authors declare that there are no competing interests associated with the manuscript.

ATRA

all-trans retinoic acid

DEG

differentially expressed gene

GO

GO Ontology

KEGG

Kyoto Encyclopedia of Genes and Genomes

MF

molecular function

ncRNA

non-coding RNA

qPCR

quantitative real-time polymerase chain reaction

TF

transcription factor