Bioinformatics analysis of the immune cell infiltration characteristics and correlation with crucial diagnostic markers in pulmonary arterial hypertension

Abstract

Background: Pulmonary arterial hypertension (PAH) is a pathophysiological syndrome, characterized by pulmonary vascular remodeling. Immunity and inflammation are progressively recognized properties of PAH, which are crucial for the initiation and maintenance of pulmonary vascular remodeling. This study explored immune cell infiltration characteristics and potential biomarkers of PAH using comprehensive bioinformatics analysis.

Methods: Microarray data of GSE117261, GSE113439 and GSE53408 datasets were downloaded from Gene Expression Omnibus database. The differentially expressed genes (DEGs) were identified in GSE117261 dataset. The proportions of infiltrated immune cells were evaluated by CIBERSORT algorithm. Feature genes of PAH were selected by least absolute shrinkage and selection operator (LASSO) regression analysis and validated by fivefold cross-validation, random forest and logistic regression. The GSE113439 and GSE53408 datasets were used as validation sets and logistic regression and receiver operating characteristic (ROC) curve analysis were performed to evaluate the prediction value of PAH. The PAH-associated module was identified by weighted gene association network analysis (WGCNA). The intersection of genes in the modules screened and DEGs was used to construct protein-protein interaction (PPI) network and the core genes were selected. After the intersection of feature genes and core genes, the hub genes were identified. The correlation between hub genes and immune cell infiltration was analyzed by Pearson correlation analysis. The expression level of LTBP1 in the lungs of monocrotaline-induced PAH rats was determined by Western blotting. The localization of LTBP1 and CD4 in lungs of PAH was assayed by immunofluorescence.

Results: A total of 419 DEGs were identified, including 223 upregulated genes and 196 downregulated genes. Functional enrichment analysis revealed that a significant enrichment in inflammation, immune response, and transforming growth factor β (TGFβ) signaling pathway. CIBERSORT analysis showed that ten significantly different types of immune cells were identified between PAH and control. Resting memory CD4⁺ T cells, CD8⁺ T cells, γδ T cells, M1 macrophages, and resting mast cells in the lungs of PAH patients were significantly higher than control. Seventeen feature genes were identified by LASSO regression for PAH prediction. WGCNA identified 15 co-expression modules. PPI network was constructed and 100 core genes were obtained. Complement C3b/C4b receptor 1 (CR1), thioredoxin reductase 1 (TXNRD1), latent TGFβ binding protein 1 (LTBP1), and toll-like receptor 1 (TLR1) were identified as hub genes and LTBP1 has the highest diagnostic efficacy for PAH (AUC = 0.968). Pearson correlation analysis showed that LTBP1 was positively correlated with resting memory CD4⁺ T cells, but negatively correlated with monocytes and neutrophils. Western blotting showed that the protein level of LTBP1 was increased in the lungs of monocrotaline-induced PAH rats. Immunofluorescence of lung tissues from rats with PAH showed increased expression of LTBP1 in pulmonary arteries as compared to control and LTBP1 was partly colocalized with CD4⁺ cells in the lungs.

Conclusion: LTBP1 was correlated with immune cell infiltration and identified as the critical diagnostic maker for PAH.

Keywords: Bioinformatics analysis; Hub gene; Immune infiltration; LTBP1; Pulmonary arterial hypertension.

Fig. 1

The flowchart of the bioinformatics analysis. GSE117261, GSE113439 and GSE53408 datasets were downloaded from the GEO database. After pre-processing and normalization of the data, the differentially expressed genes (DEGs) were identified in GSE117261 and the functional enrichment analyses of Gene Ontology and KEGG were performed. GSEA was conducted to investigate the potential biological pathways using the entire gene set. The immune landscape in the dataset was determined by the CIBERSORT algorithm. Lasso regression analysis was performed to identify 17 feature genes, and fivefold cross-validation was performed using RF and LR in GSE117261. ROC curve of 17 feature genes was performed to construct PAH prediction model in GSE53408 and GSE113439. WGCNA was performed to identify the modules associated with PAH. The intersection of genes in the modules screened and DEGs were used to construct PPI network and identification of the core genes. After the intersection of 17 feature genes and 100 core genes, four hub genes were identified. Pearson correlation analysis was performed to analyze the correlation between the hub genes and immune cell infiltration

Fig. 2

The box diagram of the gene expression matrix before and after normalization. A GSE113439 expression profile before and after normalization; B GSE113439 expression profile before and after normalization; C GSE53408 expression profile before and after normalization. The red color represents PAH lung tissue samples, and the blue color represents normal lung tissue samples

Fig. 3

Screening and functional enrichment analysis of DEGs. A The volcano plot shows DEGs between PAH and normal control. The blue dots represent downregulated DEGs, the red dots represent upregulated DEGs and the grey dots represent the genes that were not significantly changed. B The DEGs were visualized by the heatmap with red color for upregulation and blue color for downregulation. The enrichment analysis of DEGs includes GO functional analysis (C) and KEGG pathway enrichment (D); E GSEA plot showed that top five enriched KEGG pathways were positively correlated to PAH; F GSEA plot showed that TCA cycle was negatively correlated to PAH

Fig. 4

The difference in immune cell infiltration between PAH and control. A The proportions of 21 immune infiltrating cells showed in the histogram; B The clustering of 19 types of immune cells in a heatmap; C A correlation heatmap of the proportions of 19 immune cells; D The violin plot showing the differences in 21 types of immune cells between PAH and normal control

Fig. 5

Establishment of PAH prediction model. A, B LASSO coefficient spectrum for identification of feature genes of PAH and control samples from GSE117261. C The result of ROC analysis of internal validation datasets in GSE117261 using LR and RF; D Externally validated ROC analysis in the GSE113439 and GSE53408 datasets. LASSO: least absolute shrinkage and selection operator; RF: random forest; LR: logistic regression

Fig. 6

Construction of WGCNA network. A Analysis of the scale-free index and mean connectivity for various soft-threshold powers (1 ~ 20). B Dendrogram of 15 modules of genes with different colors. C Correlation heatmap showing 15 modules of different colors associated with PAH. The scatter plots showing gene distribution within the dark olive green module (D) and the dark green modules (E), respectively

Fig. 7

Construction of a PPI network. A The Venn diagram of the intersection of the upregulated DEGs and the genes in the dark olive module. B The protein interaction network of genes in the intersection. The red color represents high weight; the blue color represents low weight. C Screening of 50 upregulated hub genes. D The Venn diagram of the intersection of the downregulated DEGs and the genes in dark green modules. E Protein interaction network of genes in the intersection. The red color represents high weight and the blue color represents low weight. F Screening of 50 downregulated hub genes

Fig. 8

Identification of the hub genes. A The Venn diagram of the intersection of 17 feature genes and 100 core genes; B 4 hub genes marked in the volcano plot of all genes. C The difference of the four key genes between PAH and control. D ROC curves of 4 hub genes as independent diagnostic indicators

Fig. 9

Correlation analysis between hub genes and immune cells. A A heatmap of the correlation between the hub genes and immune cells. The red color represents a positive correlation, and the blue color represents a negative correlation, * P < 0.05; ** P < 0.01; ***P < 0.001. B The relationship between resting memory CD4.⁺ T cells with CR1, TXNRD1, LTBP1 and TLR1; C The correlation of monocytes to CR1, TXNRD1, LTBP1 and TLR1; D The correlation of neutrophils to CR1, TXNRD1, LTBP1 and TLR1

Fig. 10

LTBP1 was upregulated in monocrotaline (MCT)-induced PAH and (PDGF-BB)-induced PASMCs. A, B mPAP and RVHI were increased in rats of MCT-PAH. C HE staining of the lung in MCT-induced PAH rats and control (Scale bar = 25 μ m). D, E WT% and WA % of pulmonary arteries were calculated from HE staining. F Western blotting and quantification of LTBP1 and β-actin in the lungs of (MCT)-induced PAH (n = 9) and control (n = 8). G Immunohistochemistry analysis of LTBP1 in lungs, Scale bar = 50 μm. H Immunofluorescence of lung tissues with LTBP1 (green) and CD4 (red); Scale bar = 50 μm; I Western blot and quantification of LTBP1 and β-actin in the (PDGF-BB)-induced PASMCs and control, n = 4. Data represent the mean ± SD and Student t-test was used to compare the two groups. ****: P < 0.0001 vs Ctrl; *: P < 0.05 vs Ctrl. MCT: monocrotaline; Ctrl: control; MCT-PAH: MCT-induced pulmonary arterial hypertension