Identification of critical genes associated with human osteosarcoma metastasis based on integrated gene expression profiling

Abstract

Osteosarcoma is the most common type of malignant bone cancer, which often affects teenagers and young adults. The present study aimed to screen for critical genes and microRNAs (miRNAs/miRs) involved in osteosarcoma. A total of four microarray datasets (accession numbers GSE32981, GSE21257, GSE14827 and GSE14359) were downloaded from the Gene Expression Omnibus database. Following data preprocessing, module analysis was performed to identify the stable modules using the weighted gene co‑expression network analysis (WGCNA) package. The differentially expressed genes (DEGs) between metastatic samples and non‑metastatic samples were screened, followed by gene co‑expression network construction, and Gene Ontology function and Kyoto Encyclopedia of Genes and Genomes pathway analyses. Subsequently, prognosis‑associated genes were screened and a miRNA‑target gene regulatory network was constructed. Finally, the data for critical genes were validated. WGCNA analysis identified six modules; blue and yellow modules were significantly positively associated with osteosarcoma metastasis. A total of 1,613 DEGs were screened between primary tissue samples and metastatic samples. Following comparison of the genes in the two (blue and yellow) modules, a total of 166 DEGs were identified (metastatic samples vs. non‑metastatic samples). Functional enrichment analysis demonstrated that these DEGs were mainly involved in 'defense response', 'p53 signaling pathway' and 'lysosome'. By utilizing the clinical information in GSE21257, 10 critical genes associated with osteosarcoma prognosis were obtained, including CTP synthase 2 (CTPS2), tumor protein p53 inducible protein 3 (TP53I3) and solute carrier family 1 member 1 (SLC1A1). In addition, hsa‑miR‑422a and hsa‑miR‑194 were highlighted in the miRNA‑target gene network. Finally, matrix metallopeptidase 3 (MMP3) and vascular endothelial growth factor B (VEGFB) were predicted as critical genes in osteosarcoma metastasis. CTPS2, TP53I3 and SLC1A1 may serve major roles in osteosarcoma development, and hsa‑miR‑422a, hsa‑miR‑194, MMP3 and VEGFB may be associated with osteosarcoma metastasis.

Figure 1.

Analysis process for the four microarray datasets. DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; miRNA/miR, microRNA; WGCNA, weighted gene co-expression network analysis.

Figure 2.

Identification of stable gene modules associated with osteosarcoma as determined by weighted gene co-expression network analysis. (A) Correlation values between any two datasets from GSE21257, GSE32981, GSE14827 and GSE14359. The charts represent correlations between GSE21257-GSE32981, GSE21257-GSE14359, GSE21257-GSE14827, GSE32981-GSE14359, GSE32981-GSE14827 and GSE14359-GSE14827. (B) Cluster dendrogram based on the dynamic tree (GSE21257, GSE32981, GSE14827 and GSE14359). Different dendrogram colors represent various modules.

Figure 3.

Assessment of the stability of the modules. (A) Adjacency function definition for the genes. The left chart represents the power selection diagram of adjacency matrix weight parameter. The horizontal axis represents weight parameters of the power, while the vertical axis represents the square values of correlation coefficient between log (k) and log [p (k)]. A higher square value indicates the scale-free distribution of these data. The red line represents the standard line while square value reached 0.9. The right chart represents the mean connectivity of genes under different adjacency matrix weight parameters. (B) Multidimensional scaling plot of genes in each module. The X- and Y-axes represent the first and second principal components, respectively. (C) Cluster dendrogram of modules in the four datasets, GSE21257, GSE32981, GSE14359 and GSE14827. (D) Heat map for the correlation between each module and clinical factors. The horizontal axis represents clinical factors and the vertical axis represents different colored modules; the color changes from green to pink indicate changes from negative to positive, the numbers in the grid indicate the correlation coefficient and the numbers in parentheses indicate the significance of the correlation (P-value).

Figure 4.

DEG screening and gene co-expression network analysis. (A) Heat map for the significant DEGs. Black bars represent metastatic osteosarcoma samples and white bars represent non-metastatic osteosarcoma samples. (B) Venn diagram of key genes screened according to the WGCNA method and using the MetaDE package. (C) Co-expression network of overlapping genes. Blue and yellow represent the genes screened from blue and yellow modules, respectively. The equilateral and inverted triangles represent upregulated genes and downregulated genes; the green and gray lines represent negative and positive correlations, respectively. DEGs, differentially expressed genes; WGCNA, weighted gene co-expression network analysis.

Figure 5.

Functional annotation of the key overlapping genes in the co-expression network. (A) GO annotation. The horizontal axis represents the number of genes and the vertical axis represents the name of the GO terms. The size of the dot represents a significant P-value; larger dots and lower P-values indicate a higher significance. (B) Kyoto Encyclopedia of Genes and Genomes pathway analysis for genes in the network. The color changes from purple to light pink represent changes in significance from high to low. The numbers in each component represent the number of genes involved in a pathway. GO, Gene Ontology.

Figure 6.

Kaplan-Meier survival curve analysis for the top three genes, (A) SLC1A1, (B) CTPS2 and (C) TP53I3, associated with the prognosis of osteosarcoma. The black and red curves represent low expression and high expression sample groups, respectively. CTPS2, CTP synthase 2; SLC1A1, solute carrier family 1 member 1; TP53I3, tumor protein p53 inducible protein 3.

Figure 7.

miRNA-target gene network associated with osteosarcoma. Blue and yellow colors represent the genes screened from blue and yellow modules, respectively. The equilateral and inverted triangles represent upregulated and downregulated genes, respectively. The red squares refer to miRNAs associated with osteosarcoma. The green and gray lines represent the negative and positive connections. The red lines represent the interactions between miRNAs and target genes. miRNA/miR, microRNA.

Figure 8.

Expression and prognostic validation for critical genes, including (A) VEGFB, (B) MMP3 and (C) CTPS2. The black and red curves represent low expression and high expression osteosarcoma sample groups, respectively. CTPS2, CTP synthase 2; MMP3, matrix metallopeptidase 3; VEGFB, vascular endothelial growth factor B.

Figure 9.

Expression and prognostic validation for critical genes, including (A) MOSPD2, (B) FAP and (C) SLC1A1. The black and red curves represent low expression and high expression osteosarcoma sample groups, respectively. FAP, fibroblast activation protein α; MOSPD2, motile sperm domain containing 2; SLC1A1, solute carrier family 1 member 1.