We downloaded the datasets from the GEO database, and those that met the following criteria were included in our study: 1) Studies including both peripheral blood mononuclear cell (PBMC) samples from SSc-PH patients and PBMC samples from SSc patients without pulmonary hypertension. 2) Studies whose data and platform information were complete. Three datasets (GSE19617, GSE3346, and GSE22356) were included in our study. Specifically, GSE19617 contains 17 PBMC samples from SSc-PH patients and 25 PBMC samples from SSc patients without pulmonary hypertension, GSE33463 contains 42 PBMC samples from SSc-PH patients and 19 PBMC samples from SSc patients without pulmonary hypertension, and GSE22356 contains 10 PBMC samples from SSc-PH patients and 10 PBMC samples from SSc patients without pulmonary hypertension. Table 1 presents information about the datasets utilized in this study.

First, the array probes in the three datasets were transformed into matched gene symbols based on the platform annotation information. Then, to decrease the sample selection bias caused by the different distributions in the training and test sets, it was necessary to make the ratio of the sample size of the treatment group to the sample size of the control group in the training set close to the ratio of the sample size of the treatment group to the sample size of the control group in the test set, so we merged the mRNA expression data in GSE19617 and GSE33463 as the training set and selected GSE22356 as the test set (Bickel et al., 2007). GSE19617 was based on the GPL6480 platform, in which the mRNA expression data had been normalized by the researchers; GSE33463 and GSE22356 were based on the GPL6947 and GPL570 platforms, respectively, in which the mRNA expression data were not normalized (Pendergrass et al., 2010). We used the R package, “limma,” to normalize the mRNA expression data in the GSE33463 and GSE22356 datasets. Subsequently, to remove the batch effect caused by different platforms and different normalization methods, after studying the literature, we found that “ComBat” in the R package, “sva,” can efficiently remove the batch effect among data generated by different laboratories on account of different platforms (Johnson et al., 2007; Thillaiyampalam et al., 2017; Tang et al., 2021). Therefore, we merged the normalized mRNA expression data from GSE19617 and GSE33463 and used “ComBat” in the R package, “sva,” to remove the batch effect (Leek et al., 2012). Through our literature review, we found that among the dimensionality reduction algorithms, both t-distributed stochastic neighbor embedding (tSNE) and uniform manifold approximation and projection (UMAP) can effectively analyze sample-to-sample heterogeneity and detect batch effects (Yang Y. et al., 2021; Xiang et al., 2021). Therefore, we evaluated the efficacy of removing the batch effect by tSNE and UMAP. In addition, we also used these two methods to analyze the difficulty in distinguishing SSc-PH patients from SSc patients without pulmonary hypertension. SangerBox was used to visualize the results (Shen et al., 2022).

We obtained 113 SRP-related genes from the “REACTOME_SRP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE.v7.5.1” gene set in the MSigDB database. We used SSc-PH patients and SSc patients without pulmonary hypertension as the treatment group and control group, respectively, and used the R package, “limma,” with a p-value < 0.05 as the criterion to filter out the differentially expressed SRP-related genes (DE-SRPGs) between the treatment and control groups in the training set (Ritchie et al., 2015). The p values were calculated using the Wilcoxon rank sum test. In addition, we verified the expression patterns of the DE-SRPGs in the test set.

We performed an unsupervised clustering analysis of the SSc-PH patients in the training set based on the DE-SRPGs using the R package, “ConsensusClusterPlus,” (Wilkerson and Hayes, 2010). According to the clustering effect, the clustering stability was higher when k = 2. Therefore, we categorized the SSc-PH patients from the training set into two SRP clusters (SRPcluster A and SRPcluster B) based on the unsupervised clustering results. To further evaluate the relationships among SRPcluster A, SRPcluster B, and the control group, we performed dimensionality reduction of the training set using tSNE and UMAP based on the expression of DE-SRPGs.

To explore the differences in signal transduction between SRPcluster A and SRPcluster B, we downloaded the file, “c2.cp.kegg.v2022.1.Hs.symbols.gmt,” from the MSigDB database for gene set enrichment analysis (GSEA). We performed GSEA using the R package, “clusterProfiler,” and the statistical significance was set to an adjusted p-value of <0.05 (Yu et al., 2012; Wu et al., 2021). Then, GO annotation and KEGG pathway enrichment analysis of the DE-SRPGs were performed using the R package “clusterProfiler.” Significantly enriched signaling pathways were identified using a p-value <0.05 as the criterion. The results were visualized using the R packages, “ggplot2″ and “ComplexHeatmap” (Gu et al., 2016). A single-sample gene set enrichment analysis (ssGSEA) of 29 immune gene sets was performed using the R package, “GSVA” (Hänzelmann et al., 2013). The enrichment scores of 29 immune gene sets in each sample were calculated. Similarly, we obtained 14 SSc-PH-related pathway gene sets from the MSigDB database and performed ssGSEA on 14 SSc-PH-related pathway gene sets. Then, we compared the normalized ssGSEA scores of the treatment and control groups and the normalized ssGSEA scores of SRPcluster A and SRPcluster B. The metagenes of 14 SSc-PH-related pathways are shown in Supplementary Table S1.

For DE-SRPGs, we performed LASSO (least absolute shrinkage and selection operator) regression and SVM-RFE (support vector machine-recursive feature elimination) to identify the optimal signal recognition particle-related diagnostic genes (SRP-DGs) for SSc-PH. For both LASSO regression and SVM-RFE, the seed setting was 123. LASSO regression analysis was performed using the R package, “glmnet,” and SVM-RFE using the R package, “e1071” (Friedman et al., 2010). The SRP-related markers that were identified by the two algorithms were intersected, the intersecting genes were identified as the SRP-DGs, and the accuracy of the SRP-DGs for diagnosis in the training and test sets was evaluated using the receiver operating characteristic curve (ROC). We also compared the expression levels of SRP-DGs in SRPcluster A, SRPcluster B, and the control group.

To further analyze the diagnostic efficacy of the SRP-DGs, we constructed an SRP scoring system based on the SRP-DGs. We referred to the method of previous studies and performed a principal component analysis based on the expression levels of SRP-DGs and used principal component 1 and principal component 2 as feature scores (Sotiriou et al., 2006; Zhang et al., 2020; Zhang et al., 2022). The formula for calculating the SRPscore is: S R P s c o r e = ∑ P C 1 i + P C 2 i

In the formula, “i” represents the expressions of SRP-DGs. Then, we compared the SRPscore values of the control and treatment groups. Subsequently, we categorized the samples with SRPscore >0 as the high SRPscore group and those with SRPscore ≤0 as the low SRPscore group and analyzed the correlation between SRPscore and SRPcluster. Finally, we compared the normalized ssGSEA scores of 29 immune gene sets in the high SRPscore group with the low SRPscore group. We used ROC to evaluate the accuracy of the SRPscore values for diagnosis in the training and test sets. We then compared SRPscore values in SRPcluster A, SRPcluster B, and the control group.

To predict the risk of SSc-PH, we constructed a nomogram based on the expression levels of the SRP-DGs using the R package, “rms.” We then plotted a calibration curve to determine the extent to which the predicted values corresponded to reality. We carried out a decision curve analysis (DCA) and plotted a clinical impact curve to determine whether clinical decisions based on the nomogram model were beneficial to patients.

We constructed an ANN model using the SRP-DGs. After the gene expression data were normalized using the min-max normalization method, the seed was set to 123. An ANN model was constructed using the R package, “neuralnet.” The ANN consists of three layers: 1) Input layer, which includes the gene expressions of the seven SRP-DGs normalized by the min-max method; 2) hidden layer, which includes the gene expressions of the seven SRP-DGs normalized by the min-max method and the weights of the seven SRP-DGs; 3) output layer, which represents the results of determining whether the samples belong to the control group or treatment group. The number of neurons in the hidden layer should be two-thirds of the number of neurons in the input layer plus two-thirds of the number of neurons in the output layer, and should be in the range between the number of neurons in the input layer and the number of neurons in the output layer (Sheela and Deepa, 2013). Therefore, we set the number of neurons in the hidden layer to six and used the ROC to evaluate the predictive performance of the ANN in the training and test sets.

To assess the correlations between SRP-DGs with immune features and SSc-PH-related pathways, we calculated Spearman’s rank correlation coefficients and p values of the SRP-DGs with normalized ssGSEA scores of 29 immune gene sets and 14 SSc-PH-related pathways, which were visualized using the R package, “ggplot2.”

To further explore SRP-DGs’ functions, we performed GSEA based on single gene batch correlation analysis for each SRP-DG. The idea is to calculate Spearman’s rank correlation coefficients and p values for all genes in the training set with a single gene and to perform GSEA for genes that are significantly positively and negatively correlated with a single gene, respectively, thus simulating the possible involvement of a single gene in activation and suppression of signaling pathways. See Supplementary Table S2 for the code.

Using the Enrichr platform (https://maayanlab.cloud/Enrichr/), we entered the gene names of the SRP-DGs and screened for drugs associated with the SRP-DGs based on the DSigDB database in the “Diseases/Drugs” module with a criterion of p-value < 0.05 (Kuleshov et al., 2016).

Data heterogeneity and batch effects exist between datasets from different studies, which will adversely affect subsequent analyses if not correctly handled. Figures 1A,C show the tSNE plot and UMAP plot of the samples from the GSE19617 and GSE33463 datasets. As shown in Figures 1A,C, there was a clear difference between GSE19617 and GSE33463. Therefore, we must remove the batch effect before proceeding with the analysis. We used the “ComBat” function from the R package, “sva” to remove the batch effect. In the “sva” package, the “sva” function can be used for variable estimation, and the “ComBat” function removes batch effects, thereby reducing dependencies, stabilizing error rate estimates, and improving the reproducibility of the analysis (Leek et al., 2012). Figures 1B,D show the tSNE and UMAP plots for the samples from GSE19617 and GSE33463 after we removed the batch effect using the ComBat function. The results showed that the batch effect between GSE19617 and GSE33463 was removed and the data could then be used for subsequent analyses. In addition, tSNE (Figure 1E) and UMAP (Figure 1F) for the control group (SSc patients without pulmonary hypertension) and treatment group (SSc-PH patients) revealed no significant differences between the control and treatment groups, suggesting diagnostic difficulties.

We performed a differential analysis of 113 SRP-related genes present in the PBMC samples from SSc-PH patients versus SSc patients without pulmonary hypertension in the training set. The results showed that 30 differentially expressed SRP-related genes (DE-SRPGs) were identified using p < 0.05 as the criterion (Supplementary Table S3). Figure 2A is a box plot of the 30 DE-SRPGs. Notably, all 30 DE-SRPGs were downregulated in SSc-PH. Subsequently, we verified the expression patterns of the DE-SRPGs in the test set. Due to platform differences, the expressions of RPL10, RPL13A, RPL21, RPL23, RPL4, and RPSA were missing in the test set (GSE22356) among the 30 DE-SRPGs. In the test set, all 24 DE-SRPGs were also significantly downregulated in SSc-PH, except for 6 DE-SRPGs that were missing due to platform differences (Figure 2B, Supplementary Table S3). This suggests that SRP-related dysfunctions and defects may occur in the pathogenesis of SSc-PH. For further investigation, we clustered SSc-PH patients based on the DE-SRPGs and performed gene enrichment analysis.

To further analyze the role of SRP-related genes in SSc-PH, we performed unsupervised clustering of the PBMC samples from SSc-PH patients in the training set using the expression values of 30 DE-SRPGs with the R package, “ConsensusClusterPlus.” The consensus matrix indicates that at k = 2, the number of patients in each cluster was equally distributed, none of the clusters contained abnormally high or abnormally low numbers of patients, and the correlation between the two clusters was low (Figure 3A). When k = 2, the CDF curve was flat (Figure 3B). Figure 3C shows the variations in the area under the CDF curve for k = 2–9. Finally, SSc-PH patients from the training set were categorized into two clusters: SRPcluster A and SRPcluster B.

We performed dimensionality reduction of the training set using tSNE and UMAP based on the expression of DE-SRPGs. The tSNE plot (Supplementary Figure S1A) and the UMAP plot (Supplementary Figure S1B) indicated that SRPcluster A was closer to the control group than SRPcluster B, suggesting that SRPcluster A may be an SSc-PH subtype closer to the control group in the two clusters.

We performed gene enrichment analysis to explore the potential signaling pathways involved in the SRP gene signature. GSEA indicated that in SRPcluster A, “primary immunodeficiency,” “ribosome,” “RNA degradation,” “spliceosome,” and “T cell receptor signaling pathway” were the significant processes (Figure 3D). In SRPcluster B, the major processes included the “chemokine signaling pathway,” “complement and coagulation cascades,” “Leishmania infection,” “lysosome,” and “regulation of actin cytoskeleton” (Figure 3E).

To further explore the functions of the DE-SRPGs, we performed GO annotation (Figure 4A) and KEGG enrichment analysis (Figure 4B) on 30 DE-SRPGs. Supplementary Table S4 shows the complete GO annotation results, and Supplementary Table S5 shows the complete KEGG enrichment analysis results. The GO annotation results suggested that the 30 DE-SRPGs were predominantly enriched in functions and pathways associated with SRP. The most abundant GO biological process (BP) was protein targeting to ER (Figure 4A). Meanwhile, the most abundant GO cellular component (CC) was cytosolic ribosome (Figure 4A), and the most abundant GO molecular function (MF) was structural constituent of ribosome (Figure 4A). The results of the KEGG enrichment analysis indicated that 30 DE-SRPGs were significantly enriched in three signaling pathways: ribosome, coronavirus disease—COVID-19, and protein export (Figure 4B).

Subsequently, to determine the relationships among the DE-SRPG expression patterns and immune characteristics, we quantified the normalized ssGSEA scores of typical immune cells and pathways (Figure 5). Between SRPcluster A and SRPcluster B, there were a series of immune cells and pathways that were significantly upregulated in SRPcluster A, including APC costimulation, B cells, T-cell costimulation, T follicular helper cells (Tfh), and tumor-infiltrating lymphocytes (TIL) (Figure 5A). B cells, dendritic cells (DCs), NK cells, T helper cells, and the type II IFN response showed significant differences between the control group (SSc patients without pulmonary hypertension) and treatment group (SSc-PH patients) (Figure 5B). Among the 14 SSc-PH-related pathways, complement and coagulation cascades, complement system, endothelin, interleukin-1, interleukin-8, and osteopontin were significantly upregulated while immunoglobulin was significantly downregulated in SRPcluster B compared with SRPcluster A (Figure 6A). Complement activation, complement and coagulation cascades, complement system, interleukin-12, and troponin were significantly upregulated, while immunoglobulin and interleukin-5 were significantly downregulated in the treatment group compared to the control group (Figure 6B).

We further screened the SRP-related diagnostic genes (SRP-DGs) for SSc-PH in the DE-SRPGs with two machine learning algorithms. The results showed that we identified nine SRP-related potential diagnostic markers with the LASSO regression algorithm (Figure 7A). Meanwhile, 16 SRP-related potential diagnostic markers were identified by the SVM-RFE algorithm (Figure 7B). Finally, they were intersected to obtain eight SRP-related diagnostic genes (SRP-DGs), namely, RPL10, RPL32, RPS12, RPS14, RPS23, RPS3, RPS7, and SRP9 (Figure 7C).

From the previous results, we learned that the expression of RPL10 was missing in the test set (GSE22356) among the eight SRP-DGs due to platform differences. This leads to the fact that if the diagnostic model is constructed using all eight SRP-DGs, it will cause inconsistencies between the model in the training set and test set and cause difficulties in validation. Therefore, we used seven SRP-DGs, namely, RPL32, RPS12, RPS14, RPS23, RPS3, RPS7, and SRP9, to construct the diagnostic model for SSc-PH.

Subsequently, we plotted the ROC curves for the seven SRP-DGs. The results showed that RPL32, RPS12, RPS14, RPS23, RPS3, RPS7, and SRP9 had good diagnostic efficacy in the training set, with area under the ROC curve (AUC) values of 0.682, 0.639, 0.634, 0.615, 0.707, 0.676, and 0.644, respectively (Figure 7D). When the seven SRP-DGs were combined into one signature, the AUC value was 0.769 (Figure 7E). We also validated the diagnostic efficacy of the seven SRP-DGs in the test set. The results showed that the AUC values of RPL32, RPS12, RPS14, RPS23, RPS3, RPS7, and SRP9 were 0.875, 0.765, 0.965, 0.875, 0.940, 0.790, and 0.785 in the test set, respectively (Figure 7F). When the seven SRP-DGs were combined into one signature, the AUC value was 1.000 (Figure 7G).

As a result, these seven SRP-DGs can effectively distinguish SSc-PH patients from SSc patients without pulmonary hypertension and have better diagnostic efficacy when combined.

In addition, we compared the expression levels of seven SRP-DGs in SRPcluster A, SRPcluster B, and the control group (Supplementary Figure S2A–G). The results indicated that the expression levels of seven SRP-DGs in SRPcluster A were closer to those in the control group than in SRPcluster B.

To more accurately quantify the personalized SRP-related gene expression pattern of each patient, we constructed a scoring system, SRPscore, based on the seven SRP-DGs. Supplementary Table S6 shows the SRPscore values of the samples in the training set, and Supplementary Table S7 shows the SRPscore values of the samples in the test set. We visualized the attributes of each SSc-PH patient using alluvial plots (Figure 8A). The results showed that most SSc-PH patients with high SRPscore values belonged to SRPcluster B, whereas the majority of SSc-PH patients with low SRPscore values belonged to SRPcluster A (Figure 8A). The SRPscore values in the treatment group were significantly higher than those in the control group in both the training and test sets (Figures 8B,C).

Subsequently, we investigated whether patients in the high SRPscore group had a different type of immune cell infiltration than those in the low SRPscore group. The ssGSEA results showed that the patients in the high SRPscore group had significantly lower normalized ssGSEA scores for “B cells,” “check-point,” “T cell co-inhibition,” “T cell co-stimulation,” “Tfh” (T follicular helper cells), “Th2” (T helper 2 cells) and “TIL” (tumor-infiltrating lymphocytes) but significantly higher normalized ssGSEA scores for “parainflammation” and “type I IFN response” than those in the low SRPscore group (Figure 8D).

Then, we tested whether the SRPscore values could be used as an independent diagnostic biomarker to distinguish SSc-PH patients from SSc patients without pulmonary hypertension. The results showed that the AUC values of the SRPscore were 0.723 and 0.910 in the training set (Figure 8E) and test set (Figure 8F), respectively, thus validating the diagnostic efficacy of the SRPscore.

We compared the SRPscore values in SRPcluster A, SRPcluster B, and the control group (Supplementary Figure S2H). The results showed that the SRPscore values in SRPcluster A were closer to the control group than SRPcluster B, further demonstrating that the SRP-related gene expression patterns in SRPcluster A were closer to that of the control group.

To further investigate the relationships among the SRP-DGs and risk of SSc-PH, we constructed a nomogram model using seven SRP-DGs (RPL32, RPS12, RPS14, RPS23, RPS3, RPS7, and SRP9) to predict the risk of pulmonary hypertension complications in patients with SSc (Figure 9A). The calibration curve indicated that the nomogram model was relatively accurate in predicting SSc-PH (Figure 9B). The decision curve demonstrated that the predictions made using the nomogram model could be beneficial to patients (Figure 9C). Moreover, the clinical impact curve indicated the good predictive capacity of the nomogram model (Figure 9D).

After normalizing the expressions of the seven SRP-DGs using the min-max method, we constructed an ANN model to predict whether the samples belonged to the control group or treatment group (Figure 10A). The output results of the artificial neural network are shown in Supplementary Table S8. Then, we compared the prediction results of the ANN model with the actual grouping information and evaluated the model prediction accuracy. Subsequently, we performed ROC to evaluate the prediction performances of the ANN in the training and test sets. The results showed that the AUC values for the training and test sets were 0.999 and 0.860, respectively (Figures 10B,C). Table 2 shows the complete results of the prediction accuracies and AUC values of the ANN for the training and test sets. Overall, the ANN model was credible and has potential as an independent diagnostic predictor of SSc-PH. The results also confirmed that SRP-related genes are likely to play an essential role in the pathogenesis of SSc-PH.

We calculated Spearman correlation coefficients and p values for the expressions of seven SRP-DGs with normalized ssGSEA scores for 29 immune gene sets and 14 SSc-PH-related pathways. The results showed that the seven SRP-DGs were related to a series of immune cells, functions, and SSc-PH-related pathways (Figure 11).

For example, RPL32 was significantly positively correlated with “aDCs” (activated dendritic cells), “APC co-inhibition” (Figure 11A), and “interleukin-5” (Figure 11B) while significantly negatively correlated with “vitamin D” (Figure 11B) (p < 0.001). RPS14 was significantly positively correlated with “check-point,” “T cell co-stimulation,” and “TIL” (Figure 11A) while significantly negatively correlated with “complement activation,” “complement and coagulation cascades,” and “complement system” (Figure 11B) (p < 0.001). RPS3 was significantly positively correlated with “TIL” (Figure 11A) while significantly negatively correlated with “complement activation,” “complement and coagulation cascades,” and “complement system” (Figure 11B) (p < 0.001). RPS12 was significantly negatively correlated with “complement and coagulation cascades,” “complement system,” “endothelin,” and “troponin” (Figure 11B) (p < 0.001). RPS7 was significantly positively correlated with “interleukin-5” (Figure 11B) (p < 0.001). The results suggest that SRP-related genes may influence the immune microenvironments of SSc and SSc-PH patients and disease progression by regulating dendritic cells, T cells, and B cells. Meanwhile, SRP-related genes may regulate functions and pathways that play essential roles in pulmonary vascular remodeling, such as complement activation, the complement system, complement and coagulation cascades, endothelin, troponin, interleukin, and vitamin D.

Supplementary Figure S3 shows the results of single gene batch correlation analysis-based GSEA for RPL32, RPS3, RPS7, and RPS12. Supplementary Figure S4 shows the results of single gene batch correlation analysis-based GSEA for RPS14, RPS23, and SRP9.

The results indicated that all seven SRP-DGs might inhibit heme metabolism. RPS3, RPS7, RPS12, RPS14, RPS23, and SRP9 might inhibit coagulation. RPL32, RPS3, RPS12, RPS14, RPS23, and SRP9 might be involved in the inhibition of interferon-alpha response. RPS3, RPS12, RPS14, RPS23, and SRP9 are likely to be involved in suppressing the interferon-gamma response. RPL32, RPS3, RPS12, and RPS14 may be involved in inhibiting the complement system. RPS3, RPS7, RPS12, and RPS14 potentially activate DNA repair. In addition, the results demonstrated that SRP-DGs might be linked to a series of functions and pathways such as epithelial-mesenchymal transition, IL-6/JAK/STAT3 signaling, TNF alpha signaling, mTORC1, oxidative phosphorylation, inflammatory response, and apoptosis.

Based on the DSigDB database, we used the Enrichr platform to identify drug molecules associated with the seven SRP-DGs with a p-value < 0.05. The combined scores reflect the correlations between drugs and genes, and higher combined scores indicate stronger correlations between drugs and genes. Table 3 lists the drugs with the top ten rankings in their combined score and p values < 0.05. The results indicate that 2,6-DICHLORO-4-NITROPHENOL CTD 00000815 has a strong affinity for RPS3, while Fenbuconazole CTD 00004512 is likely to have a regulatory effect on RPS7.