The flow chart of analysis in this study is shown in Figure 1. The raw data of GSE15197 (including 13 lung samples from patients without IPAH and 18 lung samples from patients with IPAH), GSE48149 (including nine lung samples from patients without IPAH and eight lung samples from patients with IPAH), GSE113439 (including 11 lung samples from patients without IPAH and six lung samples from patients with IPAH), and GSE117261 (including 25 lung samples from patients without IPAH and 32 lung samples from patients with IPAH) were downloaded from the Gene Expression Omnibus (GEO) database.¹ GSE15197 was based on GPL6480 (Agilent-014850 Whole Human Genome Microarray 4 × 44K G4112F) and the expression matrix was extracted from raw data and normalized using “limma” package in R software (version 4.1.2; R Foundation for Statistical Computing, Vienna, Austria). GSE48149 was based on GPL16221 (Illumina HumanRef-8 v3.0 expression beadchip) and the expression matrix was obtained from raw data and normalized using “lumi” package in R software (version 4.1.2). GSE113439 and GSE117261 were based on GPL6244 (Affymetrix Human Gene 1.0 ST Array) and the expression matrixes were extracted from raw data and normalized using “oligo” package in R software (version 4.1.2). Thereafter, these expression matrixes were annotated using the “dplyr” and “limma” packages in R software (version 4.1.2). After that, the batch effect between these expression matrixes was removed using “sva” package in R software (version 4.1.2). Finally, these expression matrixes were merged as lung dataset (including 58 lung samples from patients without IPAH and 64 lung samples from patients with IPAH) for further analysis.

The raw data of GSE22356 (including ten PBMCs samples from patients without IPAH and eight PBMCs samples from patients with IPAH) and the normalized expression matrixes of GSE33463 (including 41 PBMCs samples from patients without IPAH and 30 PBMCs samples from patients with IPAH) were downloaded from the GEO database. GSE22356 was based on GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array) and the expression matrix was obtained from raw data and normalized using “affy” package in R software (version 4.1.2). Then, two expression matrixes were annotated using the “dplyr” and “limma” packages in R software (version 4.1.2), and the batch effect between two expression matrixes was removed using “sva” package in R software (version 4.1.2). Finally, two expression matrixes were merged as PBMCs dataset (including 51 PBMCs samples from patients without IPAH and 38 PBMCs samples from patients with IPAH) for further analysis.

Gene co-expression networks based on lung dataset were constructed using the WGCNA package in R software (version 4.1.2) (16). Soft-thresholding power was used to construct a weighted adjacency matrix. Relationships between a single gene and others in the analysis were incorporated, and the adjacency matrix was transformed into the topological matrix (TOM). Then, a hierarchical clustering analysis of genes was performed using 1 – TOM as the distance measure. Thereafter, modules were detected using a dynamic tree cut algorithm with a minimum module size of 50 and a minimum cut height of 0.99. The correlation between each module and IPAH was calculated and shown in a heatmap. Finally, the most relevant gene module was selected for further analysis.

The gene module which was significantly related to IPAH with highest correlation coefficient was applied to construct PPI network using the Search Tool for the Retrieval of Interacting Genes (STRING) online tool² and visualized using Cytoscape software (version 3.9.1; Institute for Systems Biology, Seattle, WA, USA). Then, cytoHubba, a plugin of Cytoscape software was used to identify the top ten key genes via the MCC method, and the corresponding PPI networks were constructed.

The Database for Annotation, Visualization and Integrated Discovery (DAVID, 2021 Update)³ were used to conduct GO and KEGG pathways enrichment analysis of genes in the most relevant gene module with the threshold of p value < 0.05.

The gene set files used in this study were downloaded from the Molecular Signatures Database (version: Human MSigDB v2022.1.Hs).⁴ The enrichment scores of GO and KEGG pathways terms in each group of lung dataset were calculated based on lung dataset using the GSEA software (version 4.2.3), and terms enriched in the IPAH group were identified. A nominal p-value (NOM p-val) of < 0.05 and false-discovery rate q value (FDR q-val) of <0.25 were considered as significantly enriched in the IPAH group.

GSVA was applied to evaluate GO and KEGG pathway terms enriched in each sample by converting the lung dataset into a gene set expression matrix using the “GSVA” package in R (version 4.1.2). After that, the differentially enriched terms between two groups were identified using R (version 4.1.2) with the threshold of p < 0.05. The differentially enriched terms were visualized using the “pheatmap” package in R (version 4.1.2).

CIBERSORT⁵ is an algorithm that can analyze the relative abundance of 22 types of immune cells in each sample, including T-cells, B-cells, and macrophages. The parameters applied in this study were as follows: (I) 100 deconvolutions (Perm) and (II) p < 0.05. The analysis was based on lung dataset using “e1071” package in R software (version 4.1.2).

The DEGs between two groups of lung dataset were identified using “limma” package in R (version 4.1.2) with the thresholds of p value < 0.05 and absolute value of log2 fold change (FC) > 0.3. The immune-related genes (IRGs) list was downloaded from Immport⁶. The IRDEGs were identified by taking the intersection between DEGs of lung dataset and IRGs list. The transcription factors of upregulated and downregulated IRDEGs were predicted using DAVID (2021 Update) with threshold of p value < 0.05. The interaction between IRDEGs and transcription factors was visualized using Cytoscape (version 3.9.1).

The DEGs between two groups of PBMCs dataset were identified using “limma” package in R (version 4.1.2) with the thresholds of p-value < 0.05 and absolute value of log2 fold change (FC) > 0.3. The IRDEGs of PBMCs dataset were identified as mentioned above.

The intersection between IRDEGs of lung dataset and IRDEGs of PBMCs dataset was taken and genes in the intersection were regarded as candidate hub genes. The predictive models for IPAH were constructed using least absolute shrinkage and selection operator (LASSO) regression and extreme gradient boosting (XGBoost) algorithm, respectively.

To construct predictive model for IPAH using LASSO regression, samples in PBMCs dataset were randomly divided into a training set (70%) and a testing set (30%). The predictive model was based on samples in the training set and candidate hub genes were selected as potential features in the model. “Glmnet” package in R (version 4.1.2) was applied to screen features and construct the model. The accuracy of the model was further validated using samples in testing set.

Similarly, samples in PBMCs dataset were randomly divided into a training set (70%) and a testing set (30%) before constructing model using XGBoost. “Xgboost” package in R (version 4.1.2) was applied to construct the model which was based on samples in the training set with candidate hub genes as features. In addition, feature importance of was analyzed and ranked.

The intersection between features in model constructed by LASSO regression and top 5 features in model constructed by XGBoost algorithm was taken, and genes in the intersection were identified as hub genes.

Blood samples were collected from patients with IPAH (n = 10) and patients without IPAH (n = 10) in Xiangya Hospital. The study was approved by the Ethics Committee of Xiangya Hospital, Central South University. Informed consent was obtained from all patients. PBMCs were separated from blood samples using Lymphocyte Separation Medium (40504ES60, Yeasen Biotechnology Co., Ltd., China) according to the manufacturer’s recommendation. Total RNA of PBMCs was extracted using RNAliquid Blood RNA Kit (RN2302, Aidlab Biotechnologies Co., Ltd., China) and reverse transcribed using the Evo M-MLV RT Mix Kit (AG11728, Accurate Biotechnology (Hunan)Co., Ltd., ChangSha, China) according to the manufacturer’s recommendations. Quantitative reverse transcription-polymerase chain reaction (RT-qPCR) was performed on a ViiA 7 system (Applied Biosystems) using All-in-One qPCR Mix (QP001, GeneCopoeia Inc.) for 40 cycles. GAPDH was used as a control for normalization. The primers used in this study are shown in Supplementary Table 1.

The Encyclopedia of RNA Interactomes (ENCORI)⁷ database was used to construct ceRNA networks of hub genes (17). The micro-RNA (miRNA) which could interact with the hub genes was predicted using ENCORI with the threshold of number of supporting crosslinking-immunoprecipitation and high-throughput sequencing (CLIP-seq) experiments > = 2 and number of target-predicting programs > = 3. Then, the interaction pairs of miRNA-circRNA were predicted using ENCORI with the threshold of number of supporting CLIP-seq experiments > = 5 and supporting degradome-seq experiments > = 3. Finally, the circRNA-miRNA-mRNA interaction network was visualized using Cytoscape software (version 3.9.1).

The relative expression levels of mRNA were presented as mean ± standard deviation (SD). The differences in proportions of immune cells and ratio of macrophage M1/M2 between different groups were analyzed by the Mann–Whitney U test. Pearson’s correlation analysis was used to analyze the relationship between different immune cells and the relationship between IRDEGs and proportions of immune cells. The diagnostic value of the predictive model was determined using receiver operating characteristic (ROC) curve and the area under ROC curve (AUC). Expression levels of mRNA in different group were compared using Student t-test. A value of p < 0.05 was considered to be statistically significant. Statistical analyses were performed using SPSS version 19 (IBM Corporation, Armonk, NY, USA) and R (version 4.1.2).

WGCNA was applied to construct IPAH-related gene co-expression networks using the lung dataset. After analysis, four gene modules were obtained (Figure 2A). Then, the correlation between IPAH and each gene module was analyzed, and the results were displayed in Figure 2B. According to the absolute value of coefficient and p-value, the gray module was the most relevant module to IPAH. PPI network of the gray module was then constructed and the top ten key genes were identified using cytoHubba via MCC method. Top ten key genes included IL1B, IFNG, CD8A, CXCL10, ICAM1, CXCL12, CCL19, ITGAM, MMP9, and CXCL13 (Figure 2C), most of which were inflammatory cytokines or chemokines. To further elucidate the functions of genes in the gray module, enrichment analysis was conducted. The results suggested that these genes were significantly enriched in terms of the GO biological process (BP) like inflammatory response, neutrophil chemotaxis, cell adhesion, chemokine-mediated signaling pathway, and immune response (Figure 2D). And the KEGG pathways in which these genes enriched were cytokine-cytokine receptor interaction, complement and coagulation cascades, hematopoietic cell lineage, chemokine signaling pathway, and extracellular matrix (ECM)-receptor interaction (Figure 2E).

Gene set enrichment analysis was conducted to further explore the terms of GO BP and KEGG pathways which enriched in IPAH group. As shown in Table 1, terms of GO BP included positive regulation of release of sequestered calcium ion into cytosol, positive regulation of glycolytic process, positive regulation of myoblast differentiation, positive regulation of cell adhesion mediated by integrin, and regulation of T cell chemotaxis were significantly enriched in IPAH group. Besides, KEGG pathways like ABC transporters, WNT signaling pathway, cell adhesion molecules, T cell receptor signaling pathway, and TGF beta signaling pathway were significantly enriched in IPAH group.

GSVA was applied to analyze the differences in pathways between different groups and the differentially activated pathways were identified. Top 10 upregulated and downregulated terms of GO BP were shown in Figure 3A. The results indicated that negative regulation of CAMP mediated signaling, regulation of T cell chemotaxis, regulation of integrin activation, transforming growth factor beta activation, and positive regulation of endothelial cell apoptotic process were significantly activated in IPAH group, while oligosaccharide lipid intermediate biosynthetic process, tricarboxylic acid metabolic process, negative regulation of complement activation, neutrophil activation involved in immune response, and NADH regeneration were significantly downregulated in IPAH group. The heatmap of differentially activated KEGG pathways between CON and IPAH groups was displayed in Figure 3B. The results demonstrated that pathways like WNT signaling pathway, TGF beta signaling pathway, regulation of autophagy, cell adhesion molecules, and ECM receptor interaction were significantly upregulated in IPAH group, while terpenoid backbone biosynthesis, glycosphingolipid biosynthesis, amino sugar and nucleotide sugar metabolism, lysosome, and citrate cycle TCA cycle were significantly downregulated in IPAH group.

The infiltrated immune cells in different samples were analyzed using CIBERSORT and the overall relative abundances of 22 types of immune cells were shown in Figure 4A. Then, the difference in proportions of immune cells between CON and IPAH groups was analyzed and the results revealed that the proportions of T cells CD4 memory resting and macrophage M1 were significantly greater in IPAH group, while the proportions of monocytes and neutrophils were significantly lower in IPAH group (Figure 4B). The correlation between 22 kinds of immune cells was further analyzed and shown in Figure 4C. Finally, the macrophage M1/macrophage M2 ratio of different samples was calculated and the result indicated that the ratio was significantly higher in IPAH group (Figure 4D), which suggested severe inflammatory state in patients with IPAH.

With the thresholds of p-value < 0.05 and absolute value of log2 FC > 0.3, there were 103 IRDEGs in lung dataset, including 56 upregulated IRDEGs and 47 downregulated IRDEGs (Figure 5A). Transcription factors of IRDEGs were predicted using DAVID and the results demonstrated that CDP, IRF2, and ISRE were the most relevant transcription factors of upregulated IRDEGs (Figure 5B). But for downregulated IRDEGs, there was no transcription factor that met the threshold of p value < 0.05. Top three upregulated IRDEGs were CXCL9, EDN1, and CXCL10, and expression levels of these genes were shown in Figure 5C. The correlation between top three upregulated IRDEGs and significantly changed immune cells was analyzed and results indicated that expression levels of CXCL9 and CXCL10 were negatively associated with the proportion of monocytes while positively associated with the proportion of macrophages M1 (Figure 5D).

Using the thresholds mentioned above, there were 100 IRDEGs in PBMCs dataset, including 56 upregulated IRDEGs and 44 downregulated IRDEGs (Figure 6A). Then, the intersection between IRDEGs of the lung dataset and PBMCs dataset was taken (Figures 6B, C). The intersection contained eight genes, five of which were upregulated genes including CXCL10, PPBP, IFIH1, STAT1, and CMTM5, and three of which were downregulated genes including IL6, VIPR1, and HSPA5. These genes were identified as candidate hub genes and expression levels of candidate hub genes in different datasets were shown in Figures 6D, E.

Predictive models were constructed based on data from the PBMCs dataset with candidate hub genes as potential features. First, LASSO regression was applied to construct the predictive model. According to LASSO regression, CXCL10, IFIH1, CMTM5, IL6, and VIPR1 were selected as features in the model, and the values of AUC of the model in training and testing sets were 0.8956 and 0.9424, respectively (Figures 7A–C). Then, XGBoost algorithm was utilized to construct the predictive model and all eight candidate hub genes were included in the model as features. The importance of features in the model was calculated and ranked (Figure 7D). The values of AUC of the model in training and testing sets were 0.9979 and 0.8303, respectively (Figures 7E, F). The intersection between features in model constructed by LASSO regression and top five features in model constructed by XGBoost algorithm was taken (Figure 7G). Genes in the intersection included CXCL10, IFIH1, and VIPR1. To further screen hub genes, mRNA expression levels of these genes in PBMCs samples from patients with or without IPAH were determined using RT-qPCR (Figures 7H-J). And the results suggested that the mRNA expression level of CXCL10 was significantly increased in PBMCs samples from patients with IPAH, while the mRNA expression level of VIPR1 was significantly decreased in PBMCs samples from patients with IPAH. In addition, the mRNA expression level of IFIH1 was significantly decreased in PBMCs samples from patients with IPAH, which was not consistent with the result of IRDEGs. Therefore, CXCL10 and VIPR1 were identified as hub genes. Considering the advantages of miRNA and circRNA as biomarkers, ceRNA networks of hub genes were then constructed. The ceRNA network of CXCL10 was constructed successfully and shown in Figure 8. But for VIPR1, there was no interaction pair of miRNA-mRNA according to the threshold of prediction.