Identification of potential gene signatures associated with osteosarcoma by integrated bioinformatics analysis

Introduction

Osteosarcoma (OS) is a type of primary malignant bone cancer that causes public health concern, especially in children and adolescents (Isakoff, Meltzer & Gorlick, 2015; Lindsey, Markel & Kleinerman, 2017). Several treatment strategies of OS such as surgical resection, traditional adjuvant chemotherapy and radiotherapy have been favored by clinical oncologists in the past few decades (Nagarajan et al., 2011). Accordingly, the 5-year survival rate has been raised to approximately 70% (Bielack et al., 2002). However, it is estimated that 80% of OS patients may suffer from the micro-metastasis, which cannot be detected at early diagnosis (Messerschmitt et al., 2009). Although multiple methods for the diagnosis and treatment of OS have been developed, new methods for the prevention and treatment of OS are still needed. The pathogenesis of OS progression remains not fully understood. Therefore, the identification of effective diagnostic makers and exploring underlying molecular etiology of OS is a pressing need.

The emergence of high-throughput sequencing technology has become an effective way to illuminate the pathogenic genes in a variety of human diseases, which help to explore pathogenesis and develop biomarkers. Many research groups have screened multiple biomarkers associated with OS by analyzing gene expression data. For example, Sun et al. evaluated the difference of genes in the expression level between OS metastasis and OS non-metastasis and discovered that TGFB1, LFT3, KDM1A, and KRAS may participate in the occurrence of OS. Xiong et al. (2015) found that CCT3, COPS3 and WWP1 were involved in the OS development by integrating gene expression and genomic aberration data. Liu, Zhao & Chen (2017) constructed a co-expression network based on a Gene Expression Omnibus (GEO) dataset and identified many potential biomarkers such as CTLA4 and PBF for diagnosis and treatment of OS. However, the molecular mechanisms of OS initiation and development have not been fully explored.

In this study, we retrieved GEO database and obtained four OS datasets. Subsequently, the differentially expressed genes (DEGs) between OS samples and normal samples were obtained and subjected to functional analysis. A protein–protein interaction (PPI) network was constructed followed by the establishment of transcription factor (TF)-target regulatory network. Following this, we downloaded two published GEO datasets of OS as the validation set for assessing the expression levels of key candidate genes. Finally, the receiver operating characteristic (ROC) analysis was conducted to evaluate the diagnostic values of key candidate genes. This study will discover novel gene signatures associated with OS, providing new trains of thought for the diagnosis and treatment of OS.

Materials and Methods

Results

GO and KEGG enrichment analysis of DEGs

The GO enrichment analysis of DEGs showed that a total of 214 GO terms were enriched, including 194 GO-BP terms, 15 GO-CC terms and 5 GO-MF terms (Table S1). Specifically, for GO-BP analysis, these DEGs were strongly associated with positive regulation of cell death, myeloid cell activation involved in immune response and regulation of protein kinase activity. Meanwhile, protein domain specific binding and transcription factor binding were significantly enriched GO-MF terms. Many DEGs were primarily involved in multiple GO-CC terms, such as anchored component of membrane and tertiary granule. The top 20 clusters of significantly enriched GO terms were displayed in Fig. 1. In addition, these DEGs were markedly enriched in 20 KEGG pathways such as regulation of lipolysis in adipocytes, protein processing in endoplasmic reticulum and pathways in cancer (Table 2). Notably, the top 20 up-regulated genes did not enrich in any KEGG pathway. However, four of top 20 down-regulated genes played vital roles in multiple significantly enriched KEGG pathways, including FABP4 (fatty acid binding protein 4), CXCL12 (C-X-C motif chemokine ligand 12), CXCL12 (C-X-C motif chemokine ligand 12) and CAT (Catalase). More specifically, FABP4 was involved in regulation of lipolysis in adipocytes and CXCL12 participated in pathways in cancer and axon guidance (Table 2). CYP4F3 was closely correlated with arachidonic acid metabolism pathway and CAT was significantly enriched in biosynthesis of amino acids and AMPK signaling pathway (Table 2).

Table 2:

The top 20 significantly enriched KEGG pathways.

ID	Term	P_value	Count	Gene Symbols
hsa04923	Regulation of lipolysis in adipocytes	0.0001	10	ADRB2, FABP4, PDE3B, PIK3CD, PLIN1, PTGER3, IRS2, ABHD5, PNPLA2, ADCY4
hsa04141	Protein processing in endoplasmic reticulum	0.0001	19	BAG1, EIF2S1, HSPA2, LMAN1, MAN1A1, MAP3K5, P4HB, EIF2AK2, RPN2, RRBP1, SSR2, PREB, SEC24A, SEC61G, UBQLN2, DNAJB11, UBQLN4, DNAJC1, UBE2J2
hsa05200	Pathways in cancer	0.0002	34	BCL2L1, CASP8, CDKN2A, CEBPA, CKS1B, COL4A1, DVL1, E2F1, EPAS1, FGF7, FGF13, FLT3, MTOR, GLI2, GLI3, GNAQ, PIK3CD, PTEN, PTGER2, PTGER3, RARB, RXRA, CXCL12, SKP2, STAT5A, TGFB3, TGFBR2, VEGFA, NCOA4, FGF16, PIAS2, ARHGEF11, LEF1, ADCY4, SHC1, CALM1, CHEK1, ELK1, FDPS, GPS2, MSX2, RANBP1, TLN1, VCAM1, TLN2
hsa01230	Biosynthesis of amino acids	0.0004	11	CTH, ENO1, GAPDH, PC, SHMT1, TALDO1, TKT, TPI1, SDS, RPIA, AADAT, CAT, MDH2, ME1, MCEE
hsa04974	Protein digestion and absorption	0.0019	11	COL2A1, COL4A1, COL5A2, COL10A1, COL11A1, COL17A1, CPA3, SLC7A8, COL18A1, COL27A1, SLC16A10
hsa04722	Neurotrophin signaling pathway	0.0022	13	CALM1, MAPK14, MAP3K5, NFKBIE, NTF3, NTRK2, PIK3CD, PRKCD, RPS6KA2, SHC1, MAGED1, PRDM4, IRAK3
hsa04152	AMPK signaling pathway	0.0023	13	CPT1A, ELAVL1, MTOR, LEP, PIK3CD, PPP2R5A, PPP2R5D, PRKAB2, IRS2, CREB5, CAMKK2, STRADB, CREB3L1, PRKCD, PRKCE, PTEN, RPS6KA2, TBC1D4, JAK2, NFKBIE, RXRA, SOCS2, CAT, ADCY4, CALM1, ELK1, FLOT2, PDE3B, PRKAR2B, SHC1, INPP5K, HSPA2
M00007	Pentose phosphate pathway, non-oxidative phase, fructose 6P =>ribose 5P	0.0026	3	TALDO1, TKT, RPIA
hsa03060	Protein export	0.0028	5	OXA1L, SRP54, SRP72, SEC61G, SRPRB
hsa04933	AGE-RAGE signaling pathway in diabetic complications	0.0041	11	COL4A1, MAPK14, JAK2, PIK3CD, PRKCD, PRKCE, STAT5A, TGFB3, TGFBR2, VCAM1, VEGFA, CD247, MTOR, NFKBIE, RXRA, STAT6, IL27RA, LHB, SHC1, SOCS2, AOX1, BCL2L1, IL5RA, LEP, PIAS2
hsa04360	Axon guidance	0.0045	16	EFNA1, EFNA3, EFNA5, EFNB1, EPHA3, EPHB2, EPHB6, FES, PIK3CD, CXCL12, SEMA7A, SEMA3A, RHOD, SEMA4C, NTNG2, PLXNA4
hsa04015	Rap1 signaling pathway	0.0053	18	CALM1, MAPK14, EFNA1, EFNA3, EFNA5, FGF7, FGF13, GNAQ, ITGB3, PFN2, PIK3CD, TLN1, VEGFA, FGF16, RAPGEF3, APBB1IP, TLN2, ADCY4, BCL2L1, COL2A1, COL4A1, MTOR, ITGA7, JAK2, PPP2R5A, PPP2R5D, PTEN, RXRA, TLR4, CREB5, CREB3L1, THEM4
hsa00590	Arachidonic acid metabolism	0.0055	8	ALOX15, GGT5, GPX3, GPX7, LTA4H, CYP4F3, PLA2G5, PTGES, GSTA4, GSTM3, SRM, GGCT
hsa04932	Non-alcoholic fatty liver disease (NAFLD)	0.0060	14	CASP8, CEBPA, EIF2S1, LEP, MAP3K5, NDUFB9, NDUFS5, PIK3CD, PRKAB2, RXRA, UQCRH, IRS2, NDUFB11, NDUFA4L2
hsa00620	Pyruvate metabolism	0.0066	6	GLO1, HAGH, MDH2, ME1, PC, LDHD
hsa04924	Renin secretion	0.0073	8	ADRB2, CALM1, GNAQ, PDE1A, PDE1C, PDE3B, PTGER2, CLCA4
hsa04750	Inflammatory mediator regulation of TRP channels	0.0101	10	CALM1, MAPK14, GNAQ, PIK3CD, PRKCD, PRKCE, PTGER2, TRPA1, TRPV4, ADCY4, ADRB2, ATP2B4, PPP2R5A, PPP2R5D, CREB5, RAPGEF3, CACNG7, CREB3L1, HSPA2, SHC1, PLA2G5, PTGIR, ARHGEF11, PPP1R14A, MYL6B, GPX3, GPX7, SLC26A4, ELK1, GNRH2, LHB
hsa01524	Platinum drug resistance	0.0144	8	BCL2L1, CASP8, CDKN2A, GSTA4, GSTM3, MAP3K5, PIK3CD, PMAIP1, BCL2A1, DFFA, DFFB, EIF2S1, PARP2, DIABLO, DAB2IP, MAPK14, VCAM1, CREB5, CREB3L1, MLKL, CALM1, ITGB3, PLAT, VEGFA, TRPV4
hsa04071	Sphingolipid signaling pathway	0.0146	11	MAPK14, MS4A2, GNAQ, MAP3K5, PIK3CD, PPP2R5A, PPP2R5D, PRKCE, PTEN, SPTLC2, SGMS1
hsa00983	Drug metabolism-other enzymes	0.0147	6	CDA, CES1, DPYD, UCK2, CES2, UCKL1

DOI: 10.7717/peerj.11496/table-2

Notes:

KEGG; Kyoto Encyclopedia of Genes and Genomes.

Evaluation the expression level and diagnostic values of key genes

An external dataset (GSE39262) was obtained from GEO database, which contained 10 human osteosarcoma cell lines and five untransformed cell lines samples. The platform for this dataset was GPL96 [HG-U133A] Affymetrix Human Genome U133A Array. GSE126209 was downloaded from the GEO database, which inclued 12 osteosarcomas tumors and 11 adjacent normal tissues samples. The platform for this dataset was GPL20301 Illumina HiSeq 4000. Seven key genes (CAMP, METTL7A, TCN1, LTF, CXCL12, CYP4F3 and NETO2) were selected to verify in GSE39262. Among them, CAMP, METTL7A, TCN1, LTF and CXCL12 acted as hub genes in PPI network. METTL7A, CYP4F3, TCN1, LTF and NETO2 were key genes in TF-gene network. Moreover, Pax-6 regulated four key genes (TCN1, CYP4F3, NETO2 and CXCL12). The gene differential expression analysis of the GSE39262 dataset revealed that NETO2 was significantly up-regulated while CXCL12, METTL7A and TCN1 were significantly down-regulated, which were consistent with our integration analysis (Fig. 4; Table S2). The gene differential expression analysis of the GSE126209 dataset displayed that NETO2 was significantly up-regulated while CXCL12 was significantly down-regulated, which were consistent with our integration analysis (Fig. 5; Table S3). ROC analysis is a commonly used method to evaluate the value of genetic diagnosis and has been used in previously biomedical works (Le et al., 2019; Le, Yapp & Yeh, 2019; Thi, Trang & Khanh, 2020). Additionally, the results of the GSE39262 dataset showed that five genes had good diagnostic values for OS (CXCL12, CYP4F3, METTL7A, NETO2 and TCN1; Fig. 6). The AUC of CXCL12 was 1.000 and the specificity and sensitivity of this model were 100.0% and 100%, respectively. The AUC of CYP4F3 was 0.840 and the specificity and sensitivity of this model were 80.0% and 80.0%, respectively. The AUC of METTL7A was 0.900 and the specificity and sensitivity of this model were 100.0% and 70.0%, respectively. The AUC of NETO2 was 0.860 and the specificity and sensitivity of this model were 100.0% and 80.0%, respectively. The AUC of TCN1 was 0.820 and the specificity and sensitivity of this model were 80.0% and 90.0%, respectively. Ultimately, the results of the GSE126209 dataset showed that four genes had good diagnostic values for OS (CXCL12, METTL7A, NETO2 and TCN1; Fig. 7). CXCL12 was 1.000 and the specificity and sensitivity of this model were 100.0% and 100%, respectively. The AUC of METTL7A was 0.856 and the specificity and sensitivity of this model were 100.0% and 83.3%, respectively. The AUC of NETO2 was 0.833 and the specificity and sensitivity of this model were 100.0% and 75.0%, respectively. The AUC of TCN1 was 0.879 and the specificity and sensitivity of this model were 81.8% and 83.3%, respectively.

Table 3:

The top 6 TF that has the most downstream genes.

TF	Number∗	Gene Symbol
Pax-4	31	FLAD1, ZZZ3, GNL2, EDARADD, RTKN, LTF, CAMP, CHL1, ADH1C, TMEM154, IRX2, FABP4, LAPTM4B, SLC28A3, FGD1, CXCL12, HBD, TCN1, CAT, SLC35F2, DYRK4, METTL7A, OLFM4, AHSP, RAB37, TMEM97, SERPINB2, CYP4F3, PDCD5, BPI, KCNG1
1-Oct	22	FLAD1, GNL2, RPAP2, EDARADD, ZC3H8, CHL1, CAMP, IRX2, FABP4, TMEM65, SLC28A3, HBD, TCN1, CAT, TMEM132C, DYRK4, METTL7A, TMEM97, HOXB6, SERPINB2, BPI, KCNG1
Nkx2-5	19	FLAD1, MLLT11, RPAP2, EDARADD, CAMP, CHL1, TMEM154, LAPTM4B, TMEM65, SLC28A3, FGD1, TCN1, SLC35F2, DYRK4, AHSP, SERPINB2, CYP4F3, KCNG1, BPI
HNF-4	16	ZZZ3, ZNF593, GNL2, LTF, CHL1, CAMP, LAPTM4B, SLC28A3, FGD1, CXCL12, HBD, SLC35F2, METTL7A, RAB37, CYP4F3, KCNG1
FOXD3	14	GNL2, RPAP2, EDARADD, ZC3H8, CHL1, LTF, FABP4, FGD1, HBD, DYRK4, METTL7A, OLFM4, NETO2, RAB37
Pax-6	13	EDARADD, GNL2, CHL1, TMEM65, SLC28A3, CXCL12, HBD, TCN1, DYRK4, NETO2, SERPINB2, CYP4F3, KCNG1

DOI: 10.7717/peerj.11496/table-3

Notes:

Number^∗

indicates the number of genes regulated by the TF

TF

transcription factor

Discussion

OS is a common malignant bone tumor and originates from mesenchymal stromal cells (MSCs) (Xiao, Hogendoorn & Cleton-Jansen, 2013). The heterogeneous histopathological characteristics and complex genomic landscape of OS have been major challenges for elaborating underlying the molecular pathogenesis of OS. In this study, we included four OS datasets and identified 1,059 DEGs (569 up-regulated DEGs and 490 down-regulated DEGs) between OS and normal samples. These genes were significantly enriched in 54 KEGG pathways such as pathways in cancer. Moreover, CAMP, METTL7A, TCN1, LTF and CXCL12 served as hub genes in PPI network. METTL7A, CYP4F3, TCN1, LTF and NETO2 were key players in TF-target gene regulatory network. Interestingly, TCN1, CYP4F3, NETO2 and CXCL12 were all regulated by Pax-6. Additionally, the expression patterns of key genes (CAMP, METTL7A, TCN1, LTF, CXCL12, CYP4F3 and NETO2) were selected to verify in two published OS datasets (GSE39262 and GSE126209).

CAMP, also known as hCAP18 or LL37, is an antimicrobial peptide gene in human (Larrick et al., 1995). The C-terminal of the protein product of CAMP contains a 37-amino acid-long peptide with broad spectrum-antibacterial activity (Vandamme, Luyten & Schoofs, 2012). There are positive expressions of CAMP in the multiple cell systems, such as epithelial cells, neutrophils and macrophages (Dhawan et al., 2015; Frew et al., 2014; Li et al., 2018). Wu et al. (2012) suggested that bone marrow stroma could express CAMP, which may be a potential ex vivo priming factor for hematopoietic stem progenitor cells to promote hematopoietic reconstitution after transplantation. Later, Coffelt et al. (2019) discovered that CAMP expression level was elevated in MSCs compared to that in ovarian cancer cells. Herein, our analysis showed that CAMP was the most down-regulated gene in patients suffering from OS. Besides, CAMP acted as a hub gene in PPI network, suggesting that this gene may be involved in the pathologic mechanism of OS. Although the underlying role of CAMP on the initiation and progression of OS has not been investigated, available evidence showes that CAMP plays significant roles in several cancers, including breast cancer, lung cancer and pancreatic cancer (García-Quiroz et al., 2016; Sainz Jr et al., 2015; Von Haussen et al., 2008). More notably, existing data indicated that CAMP had either carcinogenic or anti-cancer effects (Chen et al., 2018; Wu et al., 2010). Therefore, the influence of CAMP on OS occurrence and development needs to be further clarified in future.

Our gene differential expression revealed that CXCL12 and TCN1 were down-regulated in OS patients, which were verified in a validation dataset. Moreover, these two genes also acted as hub genes in PPI network. In addition, up-regulated NETO2 and down-regulated CYP4F3 had high degree in TF-gene regulatory network. Interestingly, CXCL12, TCN1, NETO2 and CYP4F3, regulated by Pax-6, exhibited important diagnostic values for OS. CXCL12 is also called stromal cell-derived factor-1 (SDF-1) and can bind to G-protein-coupled chemokine receptor CXCR4 (Nagasawa, 2014). Increasing studies suggested that CXCL12/CXCR4 axis played pivotal roles in tumor growth and development (Balkwill, 2004; Lu et al., 2015; Perissinotto et al., 2005). Li et al. (2018). highlighted that epigenetic regulation of CXCL12 by DNA methyltransferase 1 was associated with the metastasis and immune response in OS. Previous reports also indicated that down-regulation of CXCR4 induced OS cell apoptosis via suppressing PI3K/Akt/NF-κβ pathway (Pollino et al., 2019). However, there is no directive evidence to support the involvement of TCN1, NETO2 and CYP4F3 in OS. Notably, Pax-6 is a highly conserved evolutionarily TF and belongs to paired box TF family (Mansouri & Gruss, 1996). Several studies have pointed out that Pax6 participated in the regulation of cancer cell proliferation and progression (Shyr et al., 2010; Zong et al., 2011). Yang et al. (2019) established a TF-top 20 DEGs regulatory network by integrating and analyzing three GEO datasets (GSE66673, GSE49003 and GSE37552), and found that Pax-6 down-regulated BMP6 expression in non-metastatic OS samples. Taken together, we inferred that CXCL12/TCN1/NETO2/CYP4F3-Pax-6 axis may be implicated in the pathogenesis of OS, and four genes (CXCL12, TCN1, NETO2 and CYP4F3) were novel diagnostic biomakers for OS.

METTL7A and LTF are reported to act as tumor suppressor genes (Qi et al., 2017; Zhang et al., 2011; Zhang et al., 2015). Similarly, our findings showed that the expressions of METTL7A and LTF were decreased in OS samples. Moreover, these two genes were both hub genes in PPI analysis and key gene nodes in TF-gene regulatory analysis. These results implied that METTL7A and LTF may be correlated with underlying mechanisms of OS. However, the potential effects of METTL7A and LTF down-regulation on OS progression needs to be further investigated.

Although we have identified multiple novel gene signatures associated with OS, there are still limitations in this work. Our conclusion was drawn based on an integrated bioinformatic analysis. Therefore, additional experiments are required to confirm our findings. In addition, a larger sample size verification will also improve the reliability of our conclusion. Moreover, the clinical information should be collected to evaluate the diagnostic value of biomarkers for OS patients. Finally, the biological significances of key biomarkers will be investigated in model systems or cell lines.

In summary, a total of 1,059 DEGs were identified between OS and normal samples. Among them, up-regulation of NETO2 and down-regulation of METTL7A, TCN1, and CXCL12 may be potential gene signatures related to OS. Pax-6 was also probably associated with the pathological process of OS. However, a comprehensive bioinformatics analysis with larger sample size and in vivo or in vitro assays should be performed to confirm our results.

Supplemental Information