Heat shock protein-related diagnostic signature and molecular subtypes in ankylosing spondylitis: new pathogenesis insights

Abstract

Heat shock proteins (HSP) have been associated with a range of persistent inflammatory disorders; however, little research has been conducted on the involvement of HSP in the development of ankylosing spondylitis (AS). The research aims to identify a diagnostic signature based on HSP-related genes and determine the molecular subtypes of AS. We gathered the transcriptional data of patients with AS from the GSE73754 dataset and conducted a literature search for HSP-related genes (HRGs). The logistic regression model was utilized for the identification of hub HRGs associated with AS. Subsequently, these HRGs were employed in the construction of a nomogram prediction model. We employed a consensus clustering approach to identify novel molecular subgroups. Subsequently, we conducted functional analyses, encompassing GO, KEGG, and GSEA, to elucidate the underlying mechanisms between these subgroups. To assess the immunological landscape, we employed the xCell algorithm. Through logistic regression analysis, the four core HRGs (CCT2, HSPA6, DNAJB14, and DNAJC5) were confirmed as potential biomarkers for AS. Subsequent stratification revealed two distinct molecular phenotypes, designated as Cluster 1 and Cluster 2. Notably, Cluster 2 was characterized by the upregulation of pathways pertinent to immune response and inflammation. Our research suggests that the CCT2, HSPA6, DNAJB14, and DNAJC5 exhibit potential as effective blood-based diagnostic biomarkers for AS. These findings contribute to a deeper comprehension of the underlying mechanisms involved in the development of AS and offer potential targets for personalized therapeutic interventions.

Keywords:

• Ankylosing spondylitis
• diagnostic biomarkers
• subtypes
• heat shock protein
• immune infiltration

Introduction

Ankylosing spondylitis (AS) is a chronic inflammatory rheumatic disease primarily affecting the axial skeleton, leading to pain and progressive spinal stiffness in affected individuals [1,2]. Progressive joint ankylosis can occur as a consequence of this, leading to a significant impact on the patient’s quality of life. An early diagnosis of ankylosing spondylitis is essential in order to alleviate the patient’s burden [3]. Despite advances in understanding the etiology of AS, the precise pathogenic mechanisms remain elusive. The role of genetic predisposition is well-established, with HLA-B27 positivity being a significant risk factor [4,5]. However, not all individuals with this genetic marker develop AS, indicating that other factors, including environmental triggers and additional genetic contributions, play crucial roles in disease development [6,7].

Heat shock proteins (HSP) are highly conserved molecules that primarily function as molecular chaperones to assist in the proper folding of other proteins [8]. Emerging evidence implicates HSP in the pathogenesis of various autoimmune and inflammatory diseases due to their involvement in immune response modulation and antigen presentation [9,10]. In rheumatoid arthritis, it has been observed that HSP exhibit immunoregulatory properties owing to their immunogenic nature [11]. Notably, HSP27, 60, 70, and 90 have been implicated in this process [12]. Several investigations have provided evidence of their involvement in mitigating inflammation. The ability of HSP to suppress arthritis can be attributed to their capacity to promote the release of the anti-inflammatory factor [13,14]. Their potential role in AS, however, has been relatively unexplored. Given the stress-induced nature of HSP expression and the stress-associated pathophysiology of inflammatory diseases, investigating the involvement of HSP-related genes (HRGs) in AS could provide novel insights into the disease’s pathogenesis.

Our study focuses on delineating the contribution of HRGs to AS development. By leveraging bioinformatics approaches and the extensive transcriptional data available from the GSE73754 dataset, we aimed to identify a set of HRGs that could serve as a diagnostic signature for AS. The application of logistic regression models facilitated the identification of hub HRGs with a strong association with AS, offering potential as diagnostic biomarkers. Furthermore, the heterogeneity of AS suggests the existence of distinct molecular subtypes with unique pathogenic pathways and clinical outcomes. We applied consensus clustering to the identified HRGs to unravel these subtypes, hypothesizing that such classification could lead to personalized therapeutic strategies. Functional analyses, including Gene Ontology (GO) enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA), provided the framework for understanding the biological differences between the molecular subtypes identified. The immunological landscape of AS was also investigated using the xCell algorithm. Through this approach, we aimed to enrich our understanding of the immune cell interactions and pathways that underlie the identified molecular subtypes.

In summary, our research presents a comprehensive analysis of HRGs in the context of AS, offering new perspectives on its pathogenesis and highlighting potential blood-based biomarkers for diagnosis. The discovery of molecular subtypes within AS patients provides a stepping stone toward personalized medicine, guiding future research toward targeted therapeutic interventions that cater to the specific molecular characteristics of AS subgroups.

Methods

Dataset collection and data processing

The GSE73754 dataset, downloaded from the NCBI GEO public database, consisted of a cohort of 20 individuals in the healthy control (HC) group and 52 individuals in the AS group. Prior to the analytical process, the raw matrix files underwent extraction and standardization utilizing the “affy” software package. Subsequently, the gene expression matrix was obtained from the probe expression matrix by employing the annotations from the GPL10558 microarray platform. In addition, we initially amassed a compilation of 97 genes associated with HSP from the relevant literature [15,16].

Differential expression analysis of HRGs

The differential expression analysis between the AS and HC groups was conducted using the limma package in R software, with a significance threshold of p.adj < 0.05. The heatmap and box plot visualizations of the HRGs associated with differential expression were generated using the pheatmap package and ggplot2 package. Furthermore, the significance of each gene’s importance was determined through Friends analysis performed using the R package GOSemSim.

Functional enrichment analysis

We utilized the clusterprofler package to conduct Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses for the identified HRGs or DEGs. This analysis allowed us to gain valuable insights into the functional annotations and pathways associated with these genes.

Construction and evaluation of the predictive nomogram

HRGs were further refined through a logistic regression analysis utilizing the rms package. In addition, the nomogram was constructed using the rms package. To assess the diagnostic efficacy of both the nomogram and the gene in predicting AS, a Receiver Operating Characteristic (ROC) curve was constructed.

Consensus cluster analysis

To thoroughly investigate the association between HSP and AS, we utilized the ConsensusClusterPlus tool to execute a cluster analysis on AS patients, focusing on their HRGs expression profiles. Classification was accomplished by utilizing the most suitable “k” value obtained from the Cumulative Distribution Function (CDF) curve. Moreover, we conducted a comparative analysis of these subtypes by employing the limma package. Subsequently, the ggplot2 and ComplexHeatmap packages were utilized to generate a visually volcano plot and heatmap, respectively.

Gene set enrichment analysis (GSEA)

To explore the molecular mechanism of AS subgroups, we utilized GSEA to compare differential signaling pathways between the two subgroups. To perform this analysis, we utilized the clusterprofler package. To obtain the gene set information, we downloaded the subset “h.all.v7.4.symbols.gmt” from the Molecular Signatures Database. A statistical significance threshold of p < 0.05 was considered.

Immune microenvironment analysis

xCell was utilized in this study to quantify the infiltration level of immune cells based on gene expression profiles [17]. To evaluate the discrepancy in the proportion of infiltrating immune cells between the two subgroups, ComplexHeatmap and ggplot2 packages were employed to generate a visually heatmap and box plot. Additionally, the relationship between the expression of core HRGs and the infiltration of immune cells was investigated using the ggplot2 package.

Collection of clinical blood samples to validate core gene expression

A total of 14 individuals, including 7 healthy subjects and 7 individuals diagnosed with AS, were recruited for blood sample collection at the First Affiliated Hospital of Heilongjiang University of Traditional Chinese Medicine. Informed consent was obtained from all participants prior to sample collection, ensuring compliance with ethical guidelines. Furthermore, the study procedures were approved by the Ethics Committee at the First Affiliated Hospital of Heilongjiang University of Traditional Chinese Medicine.

Blood samples were subjected to RNA extraction using TRIzol reagent (ThermoFisher, USA). To perform qRT-PCR analysis, the SYBR qPCR Master Mix (Bio-Rad) was utilized. The expression of the primary core HRGs was assessed using the Roche LC480 Real-Time PCR System (Roche). β-actin was used as a reference gene for mRNA normalization. Therelative mRNA expression levels were quantified using the 2^–ΔΔCt method.

Results

Identification of HRGs in as

A total of 11 HRGs were identified. To elucidate the interrelationship among these 11 HRGs, a correlation analysis was conducted (Figure 1A). The findings demonstrated a strong synergistic effect between CCT2 and HSPE1 (R = 0.71), while DNAJC8 and HSPA6 exhibited a clear antagonistic action (R = −0.71). The Friends analysis revealed a descending arrangement of genes based on their mean similarity values with other genes. The gene at the topmost position, HSPA6, exhibited the highest level of similarity with other genes, thus indicating its significance as a key gene (Figure 1B). Furthermore, the expression levels of HSPA6, DNAJB14, DNAJC5, and DNAJC30 were found to be significantly elevated in AS blood samples compared to HC blood samples, while the expression levels of DNAJC8, CCT2, DNAJB6, HSPA5, HSPE1, DNAJC9, and CCT3 were significantly reduced in AS blood samples compared to HC blood samples (Figures 1C–D).

Figure 1. Identification of HRGs in AS. (A) The correlation heatmap of HRGs. Significant positive correlations are indicated by red circles, while significant negative correlations are marked by green circles. The size of the circle and the numeric value within correspond to the strength of the correlation coefficient, with statistical significance annotated by asterisks (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001). (B) Friends analysis of HRGs. To demonstrate the differential expression of HRGs between the HC and AS groups, we utilized a heatmap (C) and box plots (D). (**p < 0.01, ***p < 0.001).

Enrichment analysis of HRGs

The molecular biological functions of HRGs in AS were determined using GO terms and KEGG pathways (Figures 2A–B). In the biological process (BP) category, enrichment was observed in chaperone-mediated protein folding and protein folding. The cellular component (CC) category showed enrichment in the chaperone complex and chaperonin-containing T-complex. In the molecular function (MF) category, enrichment was observed in heat shock protein binding and unfolded protein binding. KEGG pathway analysis revealed significant enrichment in pathways such as protein processing in endoplasmic reticulum, prion disease, and lipid and atherosclerosis.

Figure 2. Enrichment analysis of HRGs. (A) GO term enrichment for BP, CC, and MF categories. Enrichment significance is represented as -Log10(P.adj), where P.adj denotes the adjusted p-value for multiple testing corrections. (B) KEGG pathway enrichment analysis. The size of each bubble indicates the count of genes associated with the pathway, while the color corresponds to the p-value. The x-axis, labeled as GeneRatio, represents the ratio of the number of differentially expressed genes involved in a given pathway to the total number of genes that make up that pathway.

Construction and evaluation of the nomogram

By conducting logistic regression analysis on the 11 HRGs, we identified four key genes, namely CCT2, HSPA6, DNAJB14, and DNAJC5, with p values < 0.05. These findings indicate that these genes possess the potential to serve as biomarkers for patients with AS (Table 1). In order to enhance the accuracy of predicting the progression of AS, we have constructed a nomogram that integrates four signature genes (CCT2, HSPA6, DNAJB14, and DNAJC5) (Figure 3A). The evaluation of this nomogram using ROC curve analysis yielded an impressive AUC value of 0.87, indicating a good predictive capacity (Figure 3B). The AUC values obtained for the signature genes were the following: CCT2-AUC: 0.809, HSPA6-AUC: 0.788, DNAJB14-AUC: 0.81, and DNAJC5-AUC: 0.718 (Figure 3C). These promising AUC values indicate the good predictive capabilities of the four candidate biomarkers for AS.

Figure 3. Evaluation of predictive biomarkers and model performance for AS risk assessment. (A) Nomogram for predicting the probability of disease risk based on the expression levels of four biomarkers: CCT2, HSPA6, DNAJB14, and DNAJC5. (B) ROC curve for the predictive model. (C) ROC curves for individual biomarkers CCT2, HSPA6, DNAJB14, and DNAJC5.

Table 1. The results of logistic regression analysis.

Characteristics	Total (N)	Univariate analysis		Multivariate analysis
		Odds Ratio (95% CI)	p value	Odds Ratio (95% CI)	p value
CCT2	72	0.003 (0.000 − 0.082)	<0.001	0.132 (0.001 − 12.655)	0.384
HSPA6	72	16.395 (3.006 − 89.407)	0.001	5.396 (0.626 − 46.483)	0.125
DNAJB14	72	1312.394 (18.474 − 93232.2041)	<0.001	455.032 (3.837 − 53964.8391)	0.012
DNAJC5	72	4908.927 (14.346 − 1679719.0099)	0.004	48.519 (0.052 − 45456.0683)	0.266

Furthermore, we conducted an analysis of marker gene expression levels in the GSE73754 dataset. As depicted in Figure 4A, the AS group exhibited decreased expression of CCT2 in comparison to the HC group, whereas the AR group demonstrated increased expression of HSPA6, DNAJB14, and DNAJC5 (p < 0.01 or p < 0.001). These findings were further corroborated in the clinical blood samples collected (Figure 4B).

Figure 4. Differential expression of candidate genes between HC and AS. Boxplot representation of gene expression levels for CCT2, HSPA6, DNAJB14, and DNAJC5 in GSE73754 dataset (A) and clinical samples (B), comparing HC and AS (**p < 0.01, ***p < 0.001).

Clustering analysis and differential gene expression profiling of as samples

As shown in Figure 5A, a heatmap generated from a hierarchical clustering analysis of HRGs expression data, revealing two main clusters of AS samples, Cluster 1 and Cluster 2. As shown in Figure 5B, a volcano plot that highlights the differential expression between Cluster 1 and Cluster 2. 958 genes that are upregulated in Cluster 1 relative to Cluster 2 are shown in red, whereas 887 genes downregulated are in blue. Genes not exhibiting statistically significant changes are colored grey. Figure 5C presents an additional heatmap focused on the most significantly differentially expressed genes between the two clusters. Overall, this analysis summarizes the clustering of AS samples based on HRGs expression profiles and identifies key genes with significant differential expression that may be critical for understanding the biological differences between the clusters.

Figure 5. Clustering analysis and differential gene expression profiling of AS samples. (A) Heatmap with hierarchical clustering based on gene expression data of HRGs, illustrating two distinct clusters (Cluster 1 and Cluster 2) of AS samples. (B) Volcano plot depicting the differential expression between Cluster 1 and Cluster 2. (C) Heatmap of the top 50 differentially expressed genes between the two clusters, arranged by hierarchical clustering.

Functional and pathway enrichment analysis of differentially expressed genes in different clusters of as

The molecular biological functions of differentially expressed genes in different clusters of AS were determined using GO terms and KEGG pathways. As shown in Figure 6A, key enriched terms include ribosome biogenesis, ribonucleoprotein complex biogenesis, and rRNA processing for BP; ribosomal subunit and ficolin-1-rich granule for CC; and structural constituent of ribosome, ribonucleoprotein complex binding, and rRNA binding for MF. KEGG pathway analysis revealed significant enrichment in pathways such as osteoclast differentiation, autophagy, natural killer cell mediated cytotoxicity, necroptosis, TNF signaling pathway, etc (Figure 6B). Furthermore, the GSEA analysis findings indicated significant disparities in signaling pathways, specifically IL6 JAK STAT3 signaling, inflammatory response, TNFA signaling via NFKB, PI3K AKT MTOR signaling, apoptosis, and hypoxia (Figure 6C). These findings provide a deeper understanding of the biological processes differentiating the AS clusters.

Figure 6. Functional and pathway enrichment analysis of differentially expressed genes in different clusters of AS. (A) GO term enrichment for BP, CC, and MF categories. (B) KEGG pathway enrichment analysis. (C) GSEA results visualized as enrichment plots for selected pathways.

Immune infiltration analysis

To further explore the immune profiles in different subgroups of AS, we employed the xCell algorithm to calculate the proportions of distinct immune cell subsets in AS. Figure 7A displays a heatmap with hierarchical clustering representing the relative abundance of various immune cell types across two distinct groups, Cluster 1 and Cluster 2. Moreover, the Cluster 2 group demonstrated significantly elevated levels of neutrophils and NKT cells in comparison to the Cluster 1 group. Conversely, the Cluster 2 group exhibited diminished levels of B cells, CD4+ memory T cells, CD4+ naive T cells, CD4+ T cells, CD4+ Tcm cells, CD8+ T cells, CD8+ Tcm cells, NK cells, and Th1 cells (Figure 7B). The results provide insights into the immune landscape of the two HSP-related clusters, highlighting specific differences in immune cell composition that may be relevant for understanding the immune responses or pathological states associated with AS.

Figure 7. The immunological landscape of HSP-related clusters. (A) The distribution of immune cell subpopulations between Cluster 1 and Cluster 2 groups is presented in the heat map. (B) The levels of immune cell subpopulations between Cluster 1 and Cluster 2 groups are displayed in the box plot. *p < 0.05, **p < 0.01, ***p < 0.001.

In addition, the Figure 8 displays a correlation matrix showing the relationship between the expression levels of specific genes (CCT2, HSPA6, DNAJB14, and DNAJC5) and various immune cell types’ abundance. For example, HSPA6 and DNAJC5 expression shows a significant positive correlation with the abundance of neutrophils, as indicated by the large red circle with four asterisks. In contrast, CCT2 shows a significant negative correlation with the abundance of neutrophils and NKT, as indicated by the large green circle with four asterisks.

Figure 8. The correlations between key genes and immunological characteristics in AS. Red circles indicate positive correlation, while green circles indicate negative correlation. The intensity of the color and the size of the circle correspond to the strength of the correlation, with a color gradient from light to dark representing increasing correlation strength from low to high. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

Discussion

AS is a chronic inflammatory disorder characterized by inflammation of the axial skeleton, leading to progressive fusion of the spine [18]. In recent years, there has been growing interest in the role of HSP in chronic and aging inflammatory diseases [19]. However, there is limited research investigating the involvement of HSP in the pathogenesis of AS. In this study, we aimed to identify a diagnostic signature based on HSP-related genes (HRGs) and determine the molecular subtypes of AS, providing new insights into the underlying mechanisms of the disease.

Through the analysis of transcriptional data from AS patients using the GSE73754 dataset, we identified four core HRGs (CCT2, HSPA6, DNAJB14, and DNAJC5) that were significantly associated with AS. These HRGs have shown potential as effective blood-based diagnostic biomarkers for AS. CCT2, functioning as a specialized autophagy receptor for aggregated proteins, may represent a critical therapeutic target in the amelioration of pathologies associated with proteostasis dysfunction [20]. CCT2 may be regarded as an independent prognostic biomarker across a spectrum of neoplasms, wherein the expression profile of CCT2 varies considerably among different tumor types and is indicative of divergent prognostic implications [21]. HSPA6, a member of the HSP70 chaperone family, exhibits partial evolutionary conservation across mammalian species. Recent research has elucidated that HSPA6 is implicated not only in oncogenesis and the progression of neoplastic conditions but also in the pathogenesis of various non-neoplastic disorders [22]. HSPA6 enhances the growth of gastric cancer through the Hippo signaling pathway, serving as a new predictive marker for prognosis and a possible target for treatment [23]. DNAJB14 is a constituent of the type II transmembrane Hsp40 protein family, exerting a pivotal role in modulating the activity of Hsp70 ATPase on the intracellular side of the endoplasmic reticulum membrane [24]. DNAJC5 plays a critical role in maintaining the proper functioning of nerve tissues. When DNAJC5 is mutated or removed in organisms such as mice, drosophila, and worms, it results in impaired synaptic activity, degenerative changes in the nervous system, and a decrease in lifespan [25,26]. Furthermore, there is a strong correlation between DNAJC5 mutations and the development of various human diseases [27,28]. Our findings suggest that these HRGs may play a crucial role in the development and progression of AS.

AS is a type of rheumatic disease characterized by immune-mediated inflammation [29]. AS disease activity can be effectively monitored by the systemic immune-inflammation index, which has shown a significant increase in AS patients [30]. Further functional analysis revealed that Cluster 2, one of the identified molecular subtypes, exhibited upregulation of pathways related to immune response and inflammation. This finding indicates that Cluster 2 may represent a more inflammatory subtype of AS, potentially providing insight into the pathogenesis of the AS. The identification of these molecular subtypes in AS may have important clinical implications. By stratifying AS patients into different subgroups based on their molecular profiles, it may be possible to develop personalized therapeutic interventions. For example, patients in Cluster 2, characterized by heightened immune response and inflammation, may benefit from targeted anti-inflammatory therapies. On the other hand, patients in Cluster 1, which displayed a distinct molecular phenotype, may require different treatment approaches. This personalized approach could lead to improved outcomes and better management of AS patients.

AS is a chronic inflammatory autoimmune disease that affects the spine and other joints [31]. In contrast to other systemic autoimmune disorders, AS displays a prominent involvement of the innate immune system. This is evident through the abnormal functioning of various innate and innate-like immune cells, including mucosal-associated invariant T cells, neutrophils, T cells, and mast cells, within the susceptible areas associated with the disease [2]. The mutual interaction between CD4 + T cells and human leukocyte antigen B27 results in a series of cytokines and chemokines, triggering inflammatory reactions and causing erosion of the bones in AS [32]. According to a recent research, individuals in the active stage of AS exhibit an immune system imbalance in terms of frequency, encompassing a variety of immune cell types such as B cells, CD4 + T cells, Tregs, CD8 + T cells, and Th cells [33]. In the present study, we employed the xCell algorithm to evaluate the immunological landscape associated with AS. The Cluster 2 group demonstrated significantly elevated levels of neutrophils and NKT cells in comparison to the Cluster 1 group. Conversely, the Cluster 2 group exhibited diminished levels of B cells, CD4+ memory T cells, CD4+ naive T cells, CD4+ T cells, CD4+ Tcm cells, CD8+ T cells, CD8+ Tcm cells, and Th1 cells. This analysis provided further insights into the immune cell composition and functional differences between the identified molecular subtypes. Understanding the immunological landscape of AS can aid in the development of novel therapeutic strategies and the identification of potential targets for intervention.

Conclusions

In conclusion, our study identified four core HSP-related genes (CCT2, HSPA6, DNAJB14, and DNAJC5) as potential blood-based diagnostic biomarkers for AS. These genes can be used in a nomogram prediction model for improved diagnosis. We also discovered two distinct molecular subgroups in AS, with Cluster 2 showing upregulation of immune response and inflammation pathways. Our findings provide valuable insights into the development of AS and offer potential targets for personalized therapeutic interventions.

Ethics approval and consent to participate

The study was approved by the Ethics Committee of the First Affiliated Hospital of Heilongjiang University of Traditional Chinese Medicine.

Consent for publication

Not applicable.

Authors’ contributions

Geqiang Wang wrote the manuscript. Yongji Li, Jiaxing Liu, Qian Zhang, and Weixin Cai analyzed the data and produced the figures. Xiaodong Li reviewed and edited the manuscript.

Acknowledgment

Not applicable.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

All data used in the present study were available from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The accession number is as follows: GSE73754).

Additional information

Funding

Not applicable.