Construction of molecular subtype model of osteosarcoma based on endoplasmic reticulum stress and tumor metastasis-related genes

Abstract

Introduction: Osteosarcoma, the prevailing primary bone malignancy among children and adolescents, is frequently associated with treatment failure primarily due to its pronounced metastatic nature.

Methods: This study aimed to establish potential associations between hub genes and subtypes for the treatment of metastatic osteosarcoma. Differentially expressed genes were extracted from patients diagnosed with metastatic osteosarcoma and a control group of non-metastatic patients, using the publicly available gene expression profile (GSE21257). The intersection of these gene sets was determined by focusing on endoplasmic reticulum (ER) stress-related genes sourced from the GeneCards database. We conducted various analytical techniques, including functional and pathway enrichment analysis, WGCNA analysis, protein-protein interaction (PPI) network construction, and assessment of immune cell infiltration, using the intersecting genes. Through this analysis, we identified potential hub genes.

Results: Osteosarcoma subtype models were developed using molecular consensus clustering analysis, followed by an examination of the associations between each subtype and hub genes. A total of 138 potential differentially expressed genes related to endoplasmic reticulum (ER) stress were identified. These genes were further investigated using Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) pathways. Additionally, the PPI interaction network revealed 38 interaction relationships among the top ten hub genes. The findings of the analysis revealed a strong correlation between the extent of immune cell infiltration and both osteosarcoma metastasis and the expression of hub genes. Notably, the differential expression of the top ten hub genes was observed in osteosarcoma clusters 1 and 4, signifying their significant association with the disease.

Conclusion: The identification of ten key genes linked to osteosarcoma metastasis and endoplasmic reticulum stress bears potential clinical significance. Additionally, exploring the molecular subtype of osteosarcoma has the capacity to guide clinical treatment decisions, necessitating further investigations and subsequent clinical validations.

Keywords: Endoplasmic reticulum stress; Metastasis; Osteosarcoma; Subtype analysis.

Fig. 1

Flow chart.

Fig. 2

Difference analysis (Pink is the metastatic group and green is the non-metastatic group): (A) Boxplot before GSE21257 data processing. (B) Boxplots after GSE21257 data processing. (C) PCA graph. (D)GSE21257 dataset differential genes volcano plot. Red is up-regulated differential genes, blue is down-regulated differential genes, and grey is non-significant genes. (E) GSE21257 datasetdifferential genes heat map. Blue is the metastatic group, pink is the non-metastatic group, green is low expression, and red is high expression. (F) Venn diagram of the intersection of endoplasmic reticulum stress-related genes and differential genes. 138 potential ER stress-related differentially expressed genes. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

1. 2

   Molecular Subtype Construction

Fig. 3

Molecular typing of osteosarcoma: (A) The cumulative distribution function (CDF) curve. (B) The CDF Delta area curve. (C) Heat map of molecular subtype sample clustering. (D) The spatial distribution of different clusters.

1. 3

   Enrichment analysis

Fig. 4

Enrichment analysis: (A) GOBP enrichment bubble chart of intersection genes. The vertical axis is the BP name, the horizontal axis is the number of enriched genes, and the dot size is the ratio of the number of enriched genes to the total number of uploaded genes. The larger the ratio, the larger the dots. The redder the dot color, the more significant the P value. (B) GOCC enrichment bubble chart of intersection genes. The vertical axis is the CC name, the horizontal axis is the number of enriched genes, and the size of the dots is the ratio of the number of enriched genes to the total number of uploaded genes. The larger the ratio, the larger the dots. The redder the dot color, the more significant the P value. (C) GOMF enrichment bubble chart of intersection genes.The vertical axis is the MF name, the horizontal axis is the number of enriched genes, and the size of the dots is the ratio of the number of enriched genes to the total number of uploaded genes. The larger the ratio, the larger the dots. The redder the dot color, the more significant the P value. (D) PathwaysEnrichment bubble chart of intersection genes. The vertical axis is the pathways name, the horizontal axis is the number of enriched genes, and the size of the dots is the ratio of the number of enriched genes to the total number of uploaded genes. The larger the ratio, the larger the dots. The redder the dot color, the more significant the P value. (E) GO enrichment histogram of intersection genes. The horizontal axis is the number of genes, and the vertical axis is the GO item. (F) KEGG network diagram. Different colors represent different types of KEGG pathways, and bold fonts represent pathways.The more significant the P value, the larger the dots, and the two-point line represents the correlation between functions. (G) TOP5 KEGG GSEA diagram. Enrichment Score polyline section.The horizontal axis is the sorted gene, and the vertical axis is the corresponding ES.The peak in the line graph is the Enrichment score of this genes set, and the genes before the peak are the core genes under the genes set. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

1. 4

   WGCNA analysis

Fig. 5

WGCNA analysis: (A) Soft Threshold. (B) Module and Trait Data Heatmap. Orange means positive correlation, blue means negative correlation, the darker the color, the stronger the correlation. (C)Module Clustering Plot. (D)Blue Significant Difference Module Scatter Plot. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

1. 5

   PPI Interaction Network

Fig. 6

PPI interaction analysis of key genes: (A) Venn diagram of the intersection of ER stress-related differential genes and key genes in the blue module. (B) PPI interaction network diagramThe blue dots are key genes, and the connectionsare interactions. (C) PPI interaction NetworkAnalyzer visualization. The smaller the significance, the larger the dots are.The thickness of edge is the combine score, and the color is down-regulated to up-regulated from blue to red. (D) Hub gene interaction diagram. The color from red to light is the MCC score. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

1. 6

   Immune infiltration correlation analysis

Fig. 7

Analysis of immune infiltration: (A) Histogram of immune infiltration distribution. The horizontal axis is the cell type, and the vertical axis is the estimated proportion. (B) Histogram of the distribution of immune infiltration samples.The horizontal axis is the sample, the vertical axis is the estimated proportion, and different colors represent different immune cells. (C) Heat map related with immune infiltration. Correlation of cellular immune infiltration in each sample, positive correlation in blue, negative correlation in red. (D) Violin plot of the correlation of cellular immune infiltration between the metastatic group and the non-metastatic group. The horizontal axis is the type of immune cells, the vertical axis is the cell immune infiltration score, the blue is the metastatic group, the red is the non-metastatic group. (E–N) Bar graph related with hub genes immune infiltration. The horizontal axis is the correlation score, and the vertical axis is the immune cell type. The size of the dots is the correlation, the higher the correlation, the larger the dots. The redder the color, the more significant the P value. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

1. 7

   Correlation analysis between Hub genes and different molecular subtypes

Fig. 8

Correlation diagram between Hub genes and different molecular subtypes (A–J) Violin plots of correlation between different molecular subtypes and hub genes.