We searched the Gene Expression Omnibus (GEO) database for studies on patients with AS and expression profiles from human blood samples were obtained. Three studies were identified and from them a dataset created of 76 cases and 45 healthy controls. Details are in Table 1.

For the two sets of gene expression data from the Illumina platform (http://www.ncbi.nlm.nih.gov/geo/), we downloaded the expression data after quantile normalization. Since there is no quantile normalization available for the GSE11886 dataset (http://www.ncbi.nlm.nih.gov/geo/), we downloaded the raw data and processed the data using the R package affy (version 4.2.0) [8] to acquire the probe expression profiles, which were further normalized using the RMA approach [9]. The result was an expression matrix of 54,675 gene probes for the 17 samples. The probes were annotated with R package hgu133plus2.db; probes matching multiple genes were deleted, and where multiple probes were matched to one gene, the median of the probes was used as the gene expression value. Subsequently, 23,313 gene expression profiles were acquired.

The R software package limma [10] was used to log${}_2$-transform and normalize the matrix data in each GEO dataset. DEGs between normal and AS samples were then detected using the limma package. The DEGs recognized from the three datasets were genetically integrated with the R package “RobustRankAggreg” [11], which is based on the Robust Rank Aggregation (RRA) approach. Those genes that consistently ranked superior to expectations based on the null hypothesis of input irrelevance could be screened with the RRA approach. Therefore, an integration of gene expression values for samples from various datasets was not performed. In common with many who have utilized the RobustRankAggreg package [12, 13], we did not correct for batch effects. p 0.05 together with 1.2-fold change in DEGs were deemed to be statistically significant.

R package clusterProfiler [14] was exploited for analyzing KEGG and Gene Ontology (GO) pathway enrichment to determine GO terms that appeared disproportionately in three categories (cellular component, molecular function, and biological processes). KEGG pathways were also analyzed by visualization with R package DOSE [15]. In both analyses, q-values (adjusted p-values for multiple tests) less than 0.05 indicate statistical significance.

We downloaded the XML file of each KEGG pathway that was enriched in the DEGs from the KEGG website [16]. The R package XML was employed to extract entries, relationships, and group relationships from these XML files. The information on gene interactions was subsequently extracted with a combination of scripts to build the KEGG gene interaction network to visualize the network on Cytoscape [17], and examine the topological characteristics of the cytohubba network [18].

Visualization was performed using Cytoscape [17], with the degree and closeness of the network computed through the use of the Cytoscape [17] plug-in cytohubba [18]. Genes of degree $>$20 and closeness $>$5 were selected and intersected with the differential genes to identify the final characteristic genes.

A support vector machine (SVM) [19] was used to build a diagnostic predictive model using genes to predict healthy versus AS individuals. As a supervised learning model, SVM uses machine learning algorithms for data analysis and pattern recognition. SVM builds a hyperplane that can be applied to categorization and regression in infinite-dimensional or high-dimensional spaces. A group of training samples, each labeled as belonging to one of two classes is used by the SVM training algorithm to create a model, and new instances are then assigned to a category according to the model, thereby making it a non-probabilistic binary linear categorization. The three datasets were categorized into two validation datasets and one training dataset. The model was built with the training dataset, and its classification capability was confirmed using the 10-fold cross-validation approach. The constructed model was, in turn, used to predict the category of samples in the validation dataset. The model’s predictive power was evaluated using the area under the receiver operating characteristic curve (ROC) curve (AUC), and the model was analyzed for its specificity and sensitivity in predicting AS.

Subjects were AS patients aged between 16 and 61 who were hospitalized in Jishuitan Hospital for the first time. The inclusion criteria were: (1) the diagnosis was in accordance with the New York classification revised in 1984; and (2) patients had not yet received medical treatment, including non-steroidal anti-inflammatory drugs or disease-modifying anti-rheumatic drugs. The exclusion criteria were: (1) individuals with a history of severe infection, surgery, or tumor; (2) Individuals with severe heart, lung, liver, and kidney dysfunction; and (3) individuals with a history of autoimmune diseases such as rheumatoid arthritis. Included in the study was a control group consisting of healthy volunteers between 24 and 60 without acute or chronic diseases, such as rheumatoid arthritis, and no long-term use of any drugs. Table 2 shows the characteristics of the AS patients and healthy control subjects.

Our experiments were in line with the principles of the Declaration of Helsinki and were approved by the ethics committee of Beijing Jishuitan Hospital (ethics code: 201611-03). From all participants, fasting whole blood was collected in sodium heparin-coated tubes (Greiner Bio-one, Chonburi, Thailand) and stored at –80 °C until further analysis. All RNA samples had an A260/A280 ratio between 1.66 and 2.58 and concentrations ranging from 10.9–143 ng/mL and thus met the criteria for gene expression analysis.

Extraction of total RNA from blood samples was implemented with Trizol reagent (Tiangen Biochemical, Beijing, China) and later extracted RNA was reverse transcribed to cDNA using FastQuant RT Kit containing gDNase (Tiangen Biochemical). Amplification was performed on an ABI7500 real-time fluorescence quantitative PCR instrument utilizing the Power SYBR® Green PCR Kit (Applied Biosystems, Foster City, CA, USA). Table 3 gives the primer sequences. The procedure followed was: pre-denaturation 95 °C for 10 min; 95 °C for 15 s; and 60 °C for 60 s for 40 cycles. Three replicates were employed. Relative expression levels were determined via the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{CT}}$ method.

All data are expressed as mean $\pm$ standard deviation (SD) calculated using the GraphPad Prism 9.0 software (GraphPad Software Inc., San Diego, CA, USA). An unpaired two-tailed Student’s t-test was performed to compare the differences between healthy and AS groups. p 0.05 was considered statistically significant.

Gene chip data from three datasets, GSE2510, GSE73754, and GSE1188655, were acquired from the GEO database. Following standardization, the gene expression distribution of each sample appeared similar (Fig. 2A). Genes with differential expression between AS and healthy samples from all three data sets were then further integrated and analyzed using limma. Altogether, 73 DEGs were identified (Supplementary Table 1), with 31 and 42 genes increased in the healthy and AS group, respectively, as indicated by the differential multiple heatmaps in each dataset (Fig. 2B). The enrichment results of the differential multiple GSEA [20] in all data sets were analyzed (Fig. 2C). In each data set, up- and downregulated genes were predominantly enriched in groups with high and low differential multiples, respectively. The only exception was observed for the upregulated genes in data set GSE2510 where they were enriched in the group with lower differential multiplicity, which may be caused by sampling error.

Fig. 2.

Determination of differentially expressed genes (DEGs) between Healthy Control and AS Samples. (A) Box plot of the overall expression level of the normalized samples (green bars: normal samples, red bars: AS samples). (B) Heatmap of the down- and up-regulated DEGs from the integrated microarray analysis. Each column and row denote a dataset and a gene, separately. In each rectangle, the color denotes the log${}_2$-FC value. The gradient from green to red goes from descending to ascending expression. (C) The enrichment results of log${}_2$-FC in each dataset of differentially expressed genes. FC, fold change.

To gain insight into the functional implications of these 73 DEGs, GO and KEGG functional enrichment analyses were implemented on the 31 upregulated and 42 downregulated DEGs. While 42 increased genes were not enriched in any KEGG pathways and GO terms, 31 decreased genes were identified as enriched in multiple KEGG pathways and GO terms. In the bioprocess group, five enriched GO terms were noted, including “cellular defense response”, “cytolysis”, and “T cell costimulation” (Fig. 3A), suggesting that these genes are associated with the cellular defense response. Classification by cellular component displayed seven enriched GO terms, for example, “lumen side of endoplasmic reticulum membrane” and “MHC protein complex” (Fig. 3B). Seven molecular functions were found to be mediated by these genes and were primarily related to “heparin-binding” and “antigen-binding” (Fig. 3C). Furthermore, 21 KEGG pathways were identified including “Type I diabetes mellitus”, “Th17 cell differentiation”, and “Th1 and Th2 cell differentiation” (Fig. 3D). Previous studies have also found a close relationship between “Th17 cell differentiation” and AS [21, 22, 23]. Together these data suggest that the downregulated DEGs are closely related to immunity.

Fig. 3.

Functional enrichment analysis of 31 downregulated DEGs. (A) Enriched GO terms in the “biological process” category. (B) Enriched GO terms in the “cellular component” category. (C) Enriched GO terms in the “molecular function” category. (D) Enriched KEGG biological pathways. The x-axis indicates the ratio of DEGs, the y-axis shows the various categories, the different colors denote levels of significance, and the size represents the number of differentially expressed genes. All enrichment results are provided in the Supplementary Table 2.

For the 21 KEGG pathways identified, the XML files were downloaded from the KEGG website, and genetic interaction information was extracted using the R package XML. Specifically, gene ID was transformed into a gene symbol, a KEGG pathway gene interaction network was built and 622 genes were derived. Genetic analysis of the crosstalk between these pathways indicates that a close relationship exists among most pathways, and some genes are involved in multiple pathways (Fig. 4A). A total of 1663 genetic interactions were obtained based on the regulatory information of genes in the KEGG pathway (Fig. 4B). In addition, most of the genes involved in multiple pathways are closely related to other genes in the network.

Fig. 4.

KEGG pathway gene interaction network analysis. (A) Crosstalk among genes; color, and size between KEGG pathways indicate the number of pathways in which a particular gene is involved. (B) The interaction network; color and size of KEGG pathway genes indicate the number of pathways in which genes participate.

We analyzed the pathways involved in the KEGG gene interaction network and found that in total, 248 (39.8%) genes were involved in only one pathway (Fig. 5A). A small number of genes are involved in multiple pathways. We also observed a power-law degree distribution in our network (Fig. 5B), where the degree of most genes is less than 20. Additionally, closure of the network revealed that the majority of the nodes have lower closeness and only a few have high closeness (Fig. 5C). The network distribution of betweenness was also computed and the majority of the nodes exhibit a lower betweenness (Fig. 5D). Nodes with high degree, high betweenness or high closeness are viewed as significant in a network.

Fig. 5.

Network topological property analysis to screen AS key genes. (A) Distribution of the number of pathways involved in the pathway gene network. (B) Degree distribution of the network. (C) Closeness distribution of the network. (D) Betweenness distribution of the network.

We considered a node to be the hub gene of the access network when it met the criteria of degree $>$20, closeness $>$5, and betweenness $>$0. Six genes met this criteria (KLRD1, HLA-DRB3, HLA-DRB5, IL2R$\beta$, CD247, CXCL10). Among these genes, HLA-DRB3 and HLA-DRB5 are HLA class II histocompatibility antigens. Major histocompatibility complex (MHC) is an important risk factor for AS; IL2RB is an important regulator of cytokine signaling and RET signaling in the immune system and has been reported as a potential transcriptional marker for AS [24]; CXCL10 is a discriminant marker of autoimmune joint disease [25, 26]. In sum, these six genes are strongly linked to the development and progression of autoimmune joint diseases, and can potentially be used as markers for AS.

To construct a diagnostic model, GSE2510 was used as the training dataset (n = 32, AS = 16) while GSE73754 (n = 72, AS = 52) and GSE11886 (n = 17, AS = 8) were used as validation datasets. Six previously chosen hub genes were utilized as features in the training dataset to acquire their corresponding expression profiles. After the z-score, the support vector machine classification model was established. A 10-fold cross-validation indicated the classification accuracy was 93.8%, with 30 of the 32 samples categorized correctly (Fig. 6A). The model’s specificity and sensitivity for AS were both 93.8%. The AUC was 0.94 (Fig. 6D). Validation of the model was then conducted with the GSE73754 dataset; 62 out of 72 samples were categorized correctly for an 86.1% classification accuracy. The model’s specificity and sensitivity for AS in this dataset were 60% and 96.2%, respectively (Fig. 6B) and the AUC was 0.78 (Fig. 6E). Finally, for the GSE11886 dataset, 15 of 17 samples were correctly categorized with 88.2% classification accuracy. The model’s specificity and sensitivity for AS were 88.9% and 87.5%, separately (Fig. 6C) and the AUC was 0.88 (Fig. 6F).

Fig. 6.

Establishment and Testing of the Diagnostic Model. (A) Categorization results of the diagnostic model for the GSE2510 training dataset. (B) Categorization results of the diagnostic model for the GSE73754 validation dataset. (C) Categorization results of the diagnostic model for the GSE11886 validation dataset. (D) Receiver operating characteristic curve (ROC) curve of the diagnostic model for the training dataset (area under the curve (AUC) = 0.9375). (E) ROC curve of the diagnostic model for the GSE73754 validation dataset (AUC = 0.781). (F) ROC curve of the diagnostic model for the GSE11886 validation dataset (AUC = 0.882).

We validated these gene expressions in clinical samples. The findings demonstrated a significant difference in IL2R$\beta$ expression versus normal controls, which agrees with the findings of our bioinformatics analysis (Fig. 7).

Fig. 7.

Expressions of IL2R$\mathrm{\beta }$ in AS samples. The expression of IL2RB was lower in AS samples relative to controls. ****p 0.0001.

These findings suggest that the diagnostic prediction model established in this research can differentiate effectively between healthy individuals and those with AS and that IL2R$\beta$ can be used as a reliable biomarker for AS diagnosis.