Landscape of CSRG
We first analyzed the differential expression levels of 278 CSRGs between IVDD and non-IVDD patients by using the "limma" package in R (Supplementary Table 1). Then we plotted histograms and heat maps to visualize 53 differentially expressed CSRGs (Figure 1A and B and Supplementary Table 2). We found that VEGFA, SP1, DUSP3, and CEBPB were highly expressed in IVDD patients and MAPKAPK5, CDKN2AIP, and PDCD10 were lowly expressed in IVDD patients compared to non-IVDD patients.
Functional and pathway enrichment of 53 DEGs
We performed enrichment analysis to further investigate the functional and molecular pathways involved in the 53 DEGs in IVDD. Biological process (BP) terms are related to cell aging (GO: 0007569), peptidyl-serine phosphorylation (GO: 0018105), and peptidyl-serine modification (GO: 0018209); cellular component (CC) terms are related to secretory granule lumen (GO: 0034774), vesicle lumen (GO: 0031983), and cytoplasmic vesicle lumen (GO: 0060205); molecular function (MF) terms are related to protein serine/ threonine kinase activity (GO: 0004674) and protein serine kinase activity (GO: 0106310) (Figure 1C and Supplementary Table 3). KEGG enrichment analysis showed that DEGs were highly enriched mainly in Chemical carcinogenesis - reactive oxygen species (hsa05208) and IL-17 signaling pathway (hsa04657) (Figure 1D and Supplementary Table 4).
Recognition of hub CSRGs
By importing the DEGs into the string database, we obtained a PPI network (Supplementary Figure 1). We next used Cytoscape for further visualization and applied the MCODE algorithm to explore the central network and genes. We finally identified the first sub-networks containing 21 candidate genes, namely CXXCL1, MMP9, VENTX, DUSP3, IRF5, MAPKAPK5, SOD1, WNT2, RAD21, TOP1, BCL6, SYK, PRKCD, SP1, CEBPB, VEGFA, BRCA1, ETS1, ATM, and MAPK14 (Figure 2A). Based on 21 candidate genes, we identified 6 key genes by the LASSO algorithm (Figure 2B). We identified 9 key genes by constructing RF models (Figure 2C). We determined that SVM-RFE had the highest accuracy when 11 key genes were selected by tenfold crossover curves (Figure 2D). Finally, we combined the results of PPI network analysis, LASSO analysis, SVM-RFE algorithm, and RF algorithm to obtain 4 hub CSRGs significantly associated with IVDD (Figure 2E). Moreover, we constructed a nomogram model to predict the incidence of IVDD based on the four selected hub genes (Figure 2F). The ROC curve analysis showed that the AUC values of DUSP3, MAPKAPK5, SP1 and VEGFA were 0.755, 0.676, 0.734 and 0.758, respectively, suggesting that the three biological markers had a high accuracy of predictive value (Supplementary Figure 2).
Identification of two different hub gene clusters
Based on these four hub genes, we used the "ConsensusClusterPlus" package in R software and the consistency clustering method to identify different hub gene patterns and finally obtained two different hub gene clusters A and B (Figure 3A and Supplementary Figure 3). Histograms were then plotted to demonstrate the differential expression levels of the four hub genes between the two clusters (Figure 3B). The histogram showed that the expression levels of DUSP3 and SP1 were higher in hub gene cluster A than in cluster B. In contrast, MAPKAPK5 expression levels were higher in cluster B. There was no significant difference in the expression levels of VEGFA in Cluster A and Cluster B. PCA suggested that the two hub gene clusters based on four Hub genes could be fully distinguished (Figure 3C). In addition, we used ssGSEA to determine the abundance of immune cells in IDD samples. This allowed us to determine the immune cell infiltration of the two hub gene clusters and the correlation between the four hub genes and immune cell levels, and histograms and heat maps were drawn to show the results (Figure 3D and E). We found that immune cells such as Activated dendritric cells, CD56dim natural killer cells, Immature dendritric cell, Macrophage, Monocyte, Type 1 T helper cell, and Plasmacytoid dendritric cells infiltrated at higher levels in hub gene clusterA (p<0.05). A heat map of the correlation between the four hub genes and the level of immune cell infiltration showed that DUSP3 was positively correlated with the level of immune cell infiltration, while VEGFA was negatively correlated with the level of immune cell infiltration. Based on two hub gene clusters, we selected 297 DEGs for GO functional enrichment analysis to explore the role of these DEGs in IDD (Figure 3F and Supplementary Table 5). The results showed that these DEGs were mainly enriched in the regulation of ion transmembrane transport (GO: 0034765), the cell cortex region (GO: 0099738) and GTP binding (GO: 0005525).
Identification of two different gene patterns
To further validate the hub gene clusters, we distinguished the different gene clusters by using a consistent clustering approach based on the 297 DEGs obtained from the two hub gene clusters. We obtained two different gene clusters and them found consistent with the hub gene pattern (Figure 4A and Supplementary Figure 4). Figure 4B shows the expression levels of 297 DEGs in two different gene clusters. We then similarly plotted histograms to show the expression levels of the four hub genes between the two different gene clusters (Figure 4C). We found that the results were similar to the hub gene cluster, where the expression levels of DSUP3 and SP1 were higher in the gene cluster. In addition, we also plotted histograms to show the difference in immune cell infiltration levels between the two gene clusters (Figure 4D). We found that the immune cell infiltration levels between the two gene clusters were similar to those of the two hub gene clusters, thus confirming once again the accuracy of the consensus clustering approach to grouping. To quantify the hub gene pattern, we used the PCA algorithm to calculate the Hub score for each sample and compared the hub scores between two different hub gene clusters or gene clusters (Figure 5A). The final result shows that hub gene cluster A and gene cluster A have a higher Hub score than hub gene cluster B and gene cluster B. In addition to this, the relationship between hub gene cluster, gene cluster, and Hub score is shown in the Sankey plot (Figure 5B).
Hub gene patterns and the role of gene patterns in IVDD
To further explore the relationship between hub gene clusters and IVDD, we examined the relationship between the expression levels of different clusters of genes that are closely associated with the development of IVDD. To further investigate the relationship between hub gene clusters and IVDD, we examined the relationship between the expression levels of different clusters and genes closely related to IVDD development, such as ATG7, CXCL8, MMP2, NLRP3, TGFB1, TIMP1. The results exhibited that the expression levels of AG7T,NLRP3,TGFB1,TIMP1 in hub gene cluster A and Genen cluster A were appreciably higher than those in hub gene cluster B and gene cluster B, while CXCL8 and MMP2 were the opposite (Figure 5C and D).
Validation of 4 hub CSRGs
To further validate our results, we used qPCR analysis and western blot to verify the expression levels of mRNA and protein in normal and IVDD cells of 4 hub CSRGs. The results exhibited that the mRNA and protein expression levels of DUSP3, SP1 and VEGFA were upregulated in IVDD cells compared to normal cells, and MAPKAPK5 was the opposite, which is consistent with the results of our analysis (Figure 6A and B). The above results confirm the reliability and stability of our analytical method.