Gene signatures with predictive and prognostic survival values in human osteosarcoma

Data collection and preprocessing

Training set: The TARGET-OS RNA-sequencing dataset (presented as fragments per kilobase million, FPKM), corresponding clinical characteristics and prognosis information were downloaded from UCSC Xena (Goldman et al., 2019) (https://xena.ucsc.edu/). Patients with expression profiles but no prognostic information and clinical characteristics were excluded. Finally, 84 patients with OS were included in a training set. FPKM data were converted to TPM data and annotated using gencode.v22.annotation.gene.probeMap.

Validation set: The gene expression data GSE21257 (Buddingh et al., 2011) (GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array) for 53 patients with OS were downloaded from the GEO database, and corrected and annotated with R software.

Construction of gene signatures

A linear regression multiple regression model was developed for the underlying expression levels of genes for prognostic risk scores. The method chosen by lasso Cox was 10-fold cross-validation. According to the median cutoff value (the cutoff value refers to the content before the brackets of the HR value in each dataset) of the risk score, patients with OS were divided into high-risk and low-risk groups. Model prediction efficiency using the training set was evaluated by Kaplan–Meier log-rank test, time-dependent ROC curve analysis, Cox regression analysis, and risk factor score for validation and test sets. A nomogram was constructed using Iasso’s guidelines.

Weighted correlation network analysis of genes

Based on the variance of gene expression in TARGET-OS data, the top 5000 genes were selected for WGCNA (Langfelder & Horvath, 2008). This analysis proceeded as: check for outliers in all samples, construct sample tree with hclust, and remove outliers according to cut height. To explore the correlation between expression data and clinical phenotypes, the sample tree and characteristic heat map were visualized. Subsequently, the strength of associations between pairs of nodes of the adjacency matrix aij was calculated as: sij = |cor (xi, xj)| aij = Sijβ. xi and xj are genes i and j. The vector of expression value, sij, indicates the strength of Pearson’s correlation coefficient between genes i and j. The aij coding network connects genes i and j. β value is a soft threshold (power value). Further, the Scale-Free Topology Fit Index (scale-free R2) range from 0 to 1 is used to determine the scale-free topological model. Selecting a set of soft threshold powers (range: 1 to 20) assists in calculating scale-free topological model fitting. The soft threshold of β = 7 was used to define the adjacency matrix. The corresponding scale-free R2 value is 0.87, suggesting a satisfactory scale-free topology model.

In coexpression networks, the highest absolute association genes were clustered into the same module to generate a clustering dendrogram. Relationships between clinical traits and risk scores were analyzed by Pearson correlation and results were visualized by heat map analysis. Genes in the module were analyzed for gene ontology (GO) and KEGG pathway enrichment. Moreover, Cytoscape (Smoot et al., 2010) (version 3.7.2) was used to visualize the weighted coexpression network. gene ontology (GO) analysis, which includes annotation of biological processes (BPs), molecular functions (MFs), and cellular components (CCs), serves as a major bioinformatics tool to annotate genes and analyze the biological processes of these genes. Huimei Wang et al. analyzed the biological classification of DEGs, showing that changes in BPs of DEGs were significantly enriched in positive regulation of associated cell response, by GO analysis (Wang et al., 2018).

Functional enrichment analysis

Clusterprofiler (Yu et al., 2012) was used to study modules related to biological function for determining functional and pathway enrichment. A multiple testing correction was performed using hypergeometric test functions and the Benjamini–Hochberg method. The GOplot (Walter, Sánchez-Cabo & Ricote, 2015) package was used to visualize the enrichment analysis.

Statistical analysis

Statistical analysis was conducted in R software (version 3.6.1) with the following packages: “glmnet” (Friedman, Hastie & Tibshirani, 2009), “survivalROC” (Heagerty, Lumley & Pepe, 2000), “WGCNA,” and “clusterProfiler.” All the statistical tests were two sided, and P-values of < 0.05 were considered statistically significant.

Construction of genes classifier for OS

Prognostically significant genes for OS in TARGET-OS data were analyzed and a total of 1151 genes were significantly associated with poor prognosis. These 1151 candidate genes were included in the lasso prognostic classifier for further screening and model construction. A twenty-gene classifier for the OS was developed (Figs. 1A and 1B). The gene information in the model is shown in Table 1. Patients were divided into a high-risk group (n = 42) and low-risk group (n = 42) based on risk scores. The median risk score was set as the cutoff (the cutoff value refers to the content before the brackets of the HR value in each dataset). Table 2 shows clinical characteristics of patients with OS in the training set according to their high risk and low risk scores. The Kaplan–Meier log-rank test suggested a significant difference between high-risk and low-risk groups in the training set (P < 0.001; Fig. 1C). In the time-dependent analysis of the ROC curve, AUCs for OS in the first, third, and fifth years were 0.94, 0.98, and 0.97, respectively (Fig. 1D).

Identification of a prognostic risk score model based on the training set

Based on univariate regression analysis (Cox’s proportional hazard model), forest plots of selected genes with p-value and hazard ratios are shown in Fig. 2A. HR, Hazard Rate in Fig. 2A, shows that the p-value of selected genes was less than 0.05, showing the statistical difference, which cannot be ignored. Only when 95% CI conclude 1, it could be proved that these selected genes signatures were not significantly associated with the prognosis of patients with OS. Although the extremum of HR is very close to 1, they show that these selected genes might be associated with the prognosis of patients with OS. Risk scores and survival status of patients with OS are shown in Figs. 2B and 2C, respectively. These results indicate that patients with high risk scores have poor outcomes compared with patients with low risk scores (Fig. 2D). Using these data, a nomogram combining the classifier with clinicopathological features to predict the survival probability of patients with different risk scores was prepared (Fig. 3A). The calibration chart demonstrated that predicted three-year and five-year survival rates were very close to observed ratios (Fig. 3B).

Validation of twenty-gene signature for survival prediction in the validation set

Subsequently, the validation set was used to assess the power of the twenty-gene signature in predicting prognosis. OS of the low-risk group was superior to the high-risk group (P < 0.05; Fig. 4A). Also, time-dependent ROC analysis is used to assess the effectiveness of risk models in predicting outcomes. Areas under ROC was 0.81, 0.78, and 0.76 at one, three and five years, respectively, suggesting that the twenty-gene signature displays good predictive power (Fig. 4B). The C-index value is 0.944. The risk score (Fig. 4C) and survival status (Fig. 4D) of patients with OS, and distribution of risk scores of twenty-gene expression profiles are shown in a heat map of 53 patients in the validation set (Fig. 4E). These results suggest that the twenty-gene signature shows good prognostic value for patient survival. This finding is further validated since patients with high risk scores were associated with poorer prognoses compared with patients with low risk scores. Previous studies have investigated gene factors in an attempt to identify new prognostic OS markers (Goh et al., 2019; Guan, Guan & Song, 2020). Hence, I have compared the prognostic gene set for osteosarcoma, Cox univariate and multivariate analysis showed that BACE2, ING2 ALOX5AP, HLA-DMB, HLA-DRA, and SPINT2 were not the independent prognostic factors for osteosarcoma (Table 3). Table 3 shows the determinants of BACE2, ING2 ALOX5AP, HLA-DMB, HLA-DRA, and SPINT2 for all the patients. In the univariate analysis, These gene signatures exhibited a negative correlation with age, gender, specific and primary tumor site, but illustrate a positive correlation with recurrence and risk score. After multiple stepwise analysis, these biomarkers correlated independently with Recurrence (hazard ratio (HR) 5.374, 95% CI [1.183–24.404], p = 0.029) and risk score (HR 9.869, 95% CI [4.663–20.887], p < 0.001).

Gene co-expression network analysis

The relationship between risk scores and gene expression profiles was evaluated using the top 5000 genes in the variance filter for WGCNA analysis. The red line (cut height = 8000) was used to remove the abnormal samples in the sample tree. TARGET-40-PAUXOZ-01A was excluded after removing outliers (Fig. 5). The sample dendrogram and trait heat map placed selected samples into different sample clusters that provided clinical trait information (Fig. 6A). Independence and average connectivity of coexpression modules were determined by power (β) and scale R² value. A series of soft thresholds and corresponding performance power were plotted. The threshold for scale R² value was set at 0.85. The power value of seven is the threshold that first reaches the scale R² value, and was chosen as the soft threshold to construct and identify coexpression modules (Figs. 6B and 6C). The expression matrix was converted to an adjacency matrix, and subsequently into a topology matrix. Based on TOM, genes are clustered based on criteria of mixed dynamic shear trees using average connection–level distance. The minimum number of genes in each module was set at 7. Finally, a total of 11 modules were identified with the WGCNA package (Fig. 6D). Statistics of gene numbers in each module are shown in Table 4. Associations of these modules with clinical characteristics (including sex, age, relapse, OS, metastasis, major cancer sites, specific cancer sites, and risk scores) are shown in Fig. 6E. T module showed the highest correlation with the risk score. Further, the magenta module also was negatively correlated with sex and OS and positively correlated with major cancer sites. GSEA analysis between high-risk and low-risk groups can better illustrate the risk-score related to the biological process. But GSEA analysis between high-risk and low-risk groups only considers the expression level of gene sets. Analysis of WGCNA is based on the clinical relevant factors, which makes it sense, although its risk-score is low.

Functional enrichment analysis in the gene coexpression network of the magenta module

GO and KEGG enrichment analysis was used to characterize biological functions of genes in the magenta module as they related to risk scores. GO analysis showed that key genes related to OS were mainly enriched in peptide–antigen binding, clathrin-coated endocytic vesicle membrane, peptide binding, and MHC class II protein complex (Fig. 7A). KEGG analysis showed that the key genes related to OS were mainly enriched in antigen processing and presentation, phagosome, cell adhesion molecules, Staphylococcus aureus infection (Fig. 7B). Subsequently, the interaction network for enrichment pathways in the magenta module was visualized (Fig. 7C). Moreover, an interaction network was constructed to visualize genes in the coexpression magenta module (Fig. 7D).