Bioinformatics and Machine-Learning Identification of Pulmonary Hypertension Biomarkers and Candidate Therapeutic Compounds

Summary

This article presents a reproducible bioinformatics workflow integrating public transcriptomic datasets, machine learning, external validation, quantitative reverse transcription PCR, Connectivity Map screening, and molecular docking to identify pulmonary hypertension biomarkers and candidate therapeutic compounds.

Abstract

This study aimed to identify pulmonary hypertension (PH)-associated molecular biomarkers and candidate small-molecule compounds using public transcriptomic data and independent validation resources. Three Gene Expression Omnibus datasets (GSE22356, GSE33463, and GSE48149) were integrated following normalization, probe annotation, and ComBat batch-effect correction. Differential expression analysis, weighted gene co-expression network analysis, functional enrichment analysis, protein-protein interaction network analysis, and three machine-learning algorithms were used to identify core feature genes. Diagnostic performance was evaluated using receiver operating characteristic curves. External validation included an independent lung-tissue cohort (GSE117261), a pulmonary artery single-cell RNA-sequencing dataset (GSE210248), and quantitative reverse transcription PCR validation in independent lung-tissue samples. Connectivity Map-based drug repositioning and molecular docking were used to screen candidate compounds. Seventy-eight differentially expressed genes were identified, and CXCL10, JUN, IFIH1, MX1, and TLR7 were selected as core feature genes. In the independent GSE117261 lung-tissue cohort, JUN showed the strongest external support, whereas replication of the other genes was variable. Quantitative reverse transcription PCR in 20 biologically independent pulmonary arterial hypertension samples and 20 control samples confirmed upregulation of all five genes. The apparent five-gene qRT-PCR model and 100 repeated stratified five-fold cross-validation analyses both yielded an area under the curve of 1.000, although the small cohort requires cautious interpretation and independent prospective validation. Single-cell analysis of GSE210248 supported altered communication between immune and structural cells and smooth muscle cell phenotypic switching. BRD-K91900765/VX-745 ranked highest in Connectivity Map screening. MAPK14/p38α, its established pharmacological target, was included as a positive-reference docking protein, whereas docking against the five biomarker-associated proteins was treated as exploratory. These findings support the five genes as candidate PH biomarkers and VX-745 as a computational drug-repositioning hypothesis requiring experimental validation.

Introduction

Pulmonary hypertension (PH) is a progressive cardiopulmonary syndrome characterized by persistently elevated pulmonary arterial pressure, increased pulmonary vascular resistance, and eventual right ventricular failure. Current hemodynamic criteria define PH as a mean pulmonary arterial pressure at rest of >20 mmHg, as measured by right heart catheterization¹. Among the different clinical subtypes, pulmonary arterial hypertension (PAH) is one of the most severe forms and is characterized by progressive pulmonary vascular remodeling. Its pathological features include endothelial dysfunction, abnormal proliferation and migration of pulmonary arterial smooth muscle cells, adventitial fibroblast activation, extracellular matrix deposition, inflammatory cell infiltration, and narrowing or obliteration of the distal pulmonary arteries². These changes indicate that PH/PAH is not only a disorder of vasoconstriction but also a complex vascular remodeling disease driven by coordinated molecular, cellular, and immune-inflammatory mechanisms.

Current PAH therapies mainly target the prostacyclin, endothelin, nitric oxide–soluble guanylate cyclase, and phosphodiesterase type 5 pathways³˒⁴. Although these treatments improve symptoms, exercise capacity, and hemodynamic parameters, their effects remain largely vasodilatory and hemodynamic. Their ability to reverse established pulmonary vascular remodeling is limited, and many patients continue to experience disease progression despite combination therapy. Therefore, identifying novel molecular biomarkers and therapeutic candidates that reflect the remodeling process represents an important unmet need. In particular, immune-inflammatory activation, interferon-related signaling, Toll-like receptor pathways, chemokine-mediated immune recruitment, and smooth muscle cell phenotypic switching have emerged as potential contributors to PH/PAH progression⁵˒⁶.

High-throughput transcriptomic datasets provide valuable resources for identifying disease-associated molecular signatures in PH/PAH. However, studies based on a single dataset are often limited by small sample sizes, batch effects, platform heterogeneity, and insufficient validation. Differential expression analysis can identify genes with altered expression but may not fully capture disease-related co-expression modules or network-level interactions. Weighted gene co-expression network analysis (WGCNA) can identify gene modules associated with disease traits, whereas protein-protein interaction (PPI) network analysis can reveal highly connected genes within biological networks. Machine-learning methods can also prioritize genes with diagnostic or classification value. However, reliance on a single algorithm may introduce model-specific bias. Integrating differential expression analysis, WGCNA, PPI network analysis, and multiple machine-learning algorithms may therefore improve the robustness of biomarker discovery.

Another major challenge in transcriptomic biomarker studies is biological interpretation. Bulk-tissue signals may reflect changes in gene expression within resident vascular cells, immune cell infiltration, or altered proportions of multiple cell populations. Single-cell RNA sequencing provides an opportunity to place bulk-derived candidate genes within a cellular context. In PH/PAH, pulmonary vascular remodeling involves endothelial cells, smooth muscle cells, fibroblasts, monocytes/macrophages, lymphocytes, and other immune or structural cells. Disease progression is also associated with altered cell-cell communication and smooth muscle cell phenotypic switching. Thus, combining bulk transcriptomic screening with single-cell validation may help determine whether candidate biomarkers are associated with immune activation, vascular structural remodeling, or an imbalance in multicellular communication.

In addition to biomarker discovery, transcriptomic signatures can be used for computational drug repositioning. The Connectivity Map (CMap) links disease-associated gene-expression profiles with small molecules that may reverse or modulate those signatures⁷. When combined with compound curation and molecular docking, this strategy can generate experimentally testable therapeutic hypotheses. Although CMap prediction and molecular docking cannot establish drug efficacy, they can prioritize candidate compounds for future target-binding assays, cell-based experiments, and animal-model validation.

An integrated and reproducible workflow was developed to identify PH/PAH biomarkers and candidate therapeutic compounds. Three public Gene Expression Omnibus transcriptomic datasets were integrated following normalization and batch-effect correction. Differential expression analysis, WGCNA, functional enrichment analysis, PPI network analysis, and three machine-learning algorithms were used to screen robust feature genes. Receiver operating characteristic analysis, an independent lung-tissue validation cohort, pulmonary artery single-cell RNA-sequencing evidence, and quantitative reverse transcription PCR validation in independent samples were used to further evaluate the selected genes. Finally, CMap-based drug repositioning and molecular docking were applied to identify candidate compounds. The novelty of the study lies in its multilayered validation framework, which connects bulk transcriptomic discovery, machine-learning prioritization, independent validation, experimental quantitative reverse transcription PCR confirmation, single-cell mechanistic interpretation, and computational compound screening. The study hypothesis was that PH/PAH is driven by a coordinated immune-inflammatory and vascular remodeling program and that robust genes within this program may serve as candidate biomarkers and provide drug-repositioning opportunities.

Protocol

The public Gene Expression Omnibus (GEO) datasets analyzed in this study contained de-identified transcriptomic data from previously published studies and did not require additional ethical approval. The Ethics Committee of Huaihua University approved the independent human lung-tissue quantitative reverse transcription PCR (qRT-PCR) validation study (approval no. 2024(A05112)). Written informed consent was obtained from all participants or their legally authorized representatives before sample collection. The approval and consent procedures applied to all 20 pulmonary arterial hypertension (PAH) and 20 control lung-tissue samples included in the qRT-PCR validation. The research tools used for this protocol are listed in the Table of Materials.

1. Collection and preprocessing of public transcriptomic datasets

Pulmonary hypertension (PH)-related microarray datasets GSE22356, GSE33463, and GSE48149 were obtained from the GEO database. PH/pulmonary arterial hypertension (PAH) and control samples were extracted according to the original phenotype annotations. Expression matrices and platform annotation files were downloaded using reproducible R scripts and the GEOquery package.

Probe annotation and gene-symbol mapping were performed consistently across the datasets. When multiple probes mapped to the same gene, the mean expression value was calculated. Quantile normalization was applied, and genes with low expression or low variance were removed. The datasets were merged, and batch effects were corrected using the ComBat algorithm in the sva package⁸. The correction was evaluated using boxplots and principal component analysis.

2. Identification of differentially expressed genes

The limma package was used to compare expression levels between PH and control samples in the batch-corrected expression matrix⁹. A linear model was fitted, and empirical Bayes statistics were applied. Differentially expressed genes were defined using an adjusted P value <0.05 and an absolute log2 fold change > 0.585. The results were visualized using volcano plots and heatmaps.

3. Weighted gene co-expression network construction

A weighted gene co-expression network was constructed using the WGCNA package¹⁰. Sample clustering was performed to detect outliers. The soft-thresholding power was selected based on the scale-free topology fit index. Gene modules were identified using the dynamic tree-cutting algorithm. Module eigengenes were correlated with the PH phenotype, and the disease-associated module with the strongest correlation was selected. Genes in the key module were intersected with the differentially expressed genes to obtain consensus genes.

4. Functional enrichment analysis

Gene Ontology biological process, cellular component, and molecular function categories were analyzed using clusterProfiler¹¹. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis was performed to identify signaling pathways¹². A P value < 0.05 and a q value < 0.2 were used as the enrichment thresholds, and enriched terms were visualized using bubble plots¹¹.

5. Protein-protein interaction network construction and hub-gene identification

The consensus gene list was submitted to the STRING database, with Homo sapiens selected as the species and an interaction confidence threshold > 0.4¹³. The interaction file was imported into Cytoscape, and the CytoHubba plug-in was used to rank genes by node degree. Highly connected genes were defined as hub genes.

6. Selection of diagnostic feature genes using machine learning

Three independent feature-selection algorithms were applied. First, least absolute shrinkage and selection operator logistic regression was performed using the glmnet package and 10-fold cross-validation to identify genes with non-zero coefficients¹⁴. Second, recursive feature elimination with support vector machines was applied to remove redundant features and select the feature subset that achieved the highest cross-validation accuracy¹⁵. Third, a random forest model was constructed, and the features were ranked by mean decrease in Gini impurity¹⁶. The intersection of the gene sets derived from the three algorithms was used to define the final set of core feature genes. The pROC package was used to generate receiver operating characteristic curves and calculate area under the curve values¹⁷.

7. Validation of core genes using independent bulk and single-cell datasets

GSE117261 was used as an independent external lung-tissue validation cohort containing 58 PAH samples and 25 failed-donor control samples¹⁸. This dataset was not used in the discovery differential-expression analysis, weighted gene co-expression network construction, or machine-learning feature selection. The expression matrix was normalized and annotated, and differential expression was analyzed using limma v3.68.0. The Benjamini-Hochberg false-discovery-rate correction was applied across the entire annotated transcriptome. Single-gene receiver operating characteristic (ROC) curves were calculated using pROC v1.19.0.1, DeLong 95% confidence intervals, and Youden-index cutoffs. An exploratory five-gene logistic regression model was fitted within GSE117261, and its internal performance was additionally evaluated using repeated nested cross-validation.

GSE210248 (Table 1) was used as a single-cell pulmonary artery validation dataset containing samples from three patients with PAH and three healthy donors¹⁹. The data were processed using Seurat v5.5.1 for quality control, normalization, dimensionality reduction, clustering, and cell annotation²⁰. Major cell populations, including endothelial cells, smooth muscle cells, fibroblasts, monocytes/macrophages, and T/natural killer cells, were identified. Cell-cell communication was analyzed using CellChat v2.1.2 and the CellChatDB.human ligand-receptor database²¹. A CellChat object was created from the normalized Seurat expression matrix and cell-type metadata. Overexpressed genes and ligand-receptor interactions were identified; communication probabilities were calculated; interactions involving cell groups with fewer than 10 cells were removed; and pathway-level communication networks were inferred and aggregated. This dataset was used only for external mechanistic validation and not for model training.

Item	Description
Dataset	GSE210248
Data type	10x Genomics/droplet-based single-cell RNA sequencing; high-throughput transcriptomic profiling
Human samples	Three PAH pulmonary artery samples and three healthy donor pulmonary artery samples
Tissue source	Ex vivo pulmonary artery tissue, primarily reflecting the cellular ecology of the pulmonary vascular wall and vascular-remodeling process
Main analytical purpose	Cell-type localization, smooth-muscle-cell phenotypic switching, immune-structural cell communication, and mechanistic consistency validation of candidate genes

Table 1: Basic information for the GSE210248 single-cell validation dataset. The table summarizes the dataset accession, sequencing platform, tissue source, sample composition, and analytical purpose of the single-cell pulmonary artery validation analysis.

8. Validation of gene expression by qRT-PCR

The qRT-PCR validation included 20 biologically independent PAH lung-tissue samples from patients with PH/PAH and 20 biologically independent control lung-tissue samples. Total RNA was extracted using the Total RNA Extraction Kit. RNA concentration and purity were assessed using a spectrophotometer, and RNA integrity was evaluated by agarose gel electrophoresis. Only RNA samples with A260/280 values between 1.8 and 2.1 and no visible degradation were included.

Equal amounts of RNA were reverse-transcribed into complementary DNA using the Solarbio Universal RT-PCR Kit (AMV; catalog no. RP1200). Quantitative PCR for CXCL10, JUN, IFIH1, MX1, and TLR7 was performed using SYBR Green PCR Master Mix on a Real-Time PCR System. Each biological sample was analyzed in three technical replicates, together with no-template and no-reverse-transcription controls. The mean Ct value of the three technical replicates was used for subsequent analysis; technical replicates were not treated as independent observations. Primers spanning exon-exon junctions and producing 80–200 bp amplicons were used (Table 2). Primer specificity was verified using NCBI Primer-BLAST and melting-curve analysis²².

β-actin (ACTB) was used as the internal reference gene to normalize the expression levels of the target genes. Relative expression was calculated using the 2^-ΔΔCt method²³. Two-sided Mann-Whitney U tests were used for between-group comparisons based on the data distribution, and Benjamini-Hochberg false-discovery-rate correction was applied across the five genes. Single-gene ROC curves were generated with DeLong 95% confidence intervals, and optimal cutoffs were selected using the Youden index. The five-gene logistic regression model was initially fitted and evaluated on the same 40 biological samples; this estimate was therefore defined as the apparent in-sample performance. To assess potential overfitting, 100 stratified five-fold cross-validation repeats were performed using an L2-regularized logistic regression model, and the pooled out-of-fold ROC performance was calculated.

Gene	RefSeq accession	Forward primer (5′–3′)	Reverse primer (5′–3′)	Product size (bp)	Tm (°C)	Exon-spanning
CXCL10	NM_001565.4	GTCAAGCCAT AATTGTTC	ATAGTGCCAG GGTAGAGT	141	46.1	Yes
JUN	NM_002228.4	ACAAGTGGCA GAGTCCCG	CGCCCAAGTT CAACAACC	152	54.5	Yes
IFIH1	NM_022168	GCACAGAGCG GTAGACCCT	GCCCTGAAGC ACGAGATG	182	54.7	Yes
MX1	NM_002462.5	TTAGCCGTGG TGATTTAGC	CAAGGTGGAG CGATTCTG	156	52.3	Yes
TLR7	NM_016562.4	ATTGCCCTCGT TGTTATA	TTCCTGGAGTT TGTTGAT	179	48.1	Yes
ACTB	NM_001101.3	CTCACCATGGAT GATGATATCGC	AGGAATCCTTCT GACCCATGC	194	56.2	Yes

Table 2: Primer sequences used for quantitative reverse transcription PCR. The table lists the target genes, RefSeq accession numbers, forward and reverse primer sequences, product sizes, melting temperatures, and exon-spanning status of the primers used for qRT-PCR.

9. Candidate compound screening and molecular docking

The upregulated and downregulated core-gene signatures were submitted to the Connectivity Map database to identify small molecules predicted to reverse the PH-associated expression profile⁷. Candidates were ranked by Logit score and prediction probability.

The three-dimensional structure of BRD-K91900765/VX-745 was obtained from PubChem under CID 3038525²⁴. Compound-related pharmacological information was curated from public drug databases, and structural descriptors were calculated using DrugBank and SwissADME^25,26. Protein structures were obtained from the RCSB Protein Data Bank using the following PDB identifiers: CXCL10, 1LV9; JUN, 1JUN; IFIH1, 3B6E; MX1, 5GTM; TLR7, 7CYN; and MAPK14/p38α, 1OUK²⁷. Blind cavity detection and molecular docking were performed using CB-Dock2 v2.0 with the AutoDock Vina v1.2.0 scoring engine^28,29. Protein and ligand files were uploaded to CB-Dock2, candidate cavities were automatically detected, and docking was performed within the cavity-specific boxes generated by the server. For each protein, the cavity identifier, Vina score, cavity volume, docking-box center, docking-box dimensions, and protein-ligand complex file were recorded. The pose with the most negative Vina score was selected as the top-ranked predicted conformation. MAPK14/p38α was included as the established pharmacological target and positive-reference docking protein for VX-745. Docking against CXCL10, JUN, IFIH1, MX1, and TLR7 was exploratory and was not interpreted as evidence of direct pharmacological targeting, binding, inhibition, or efficacy.

10. Statistical analysis and reproducibility control

All statistical analyses were performed in R unless otherwise specified. Two-sided P values < 0.05 were considered statistically significant. Multiple-testing correction was applied to the differential-expression, enrichment, external-validation, and qRT-PCR analyses as specified above. Cross-validation was used to assess the stability of the machine-learning and combined qRT-PCR models.

Results

Public transcriptomic data preprocessing and differentially expressed gene identification

Integration and ComBat correction of GSE22356, GSE33463, and GSE48149 reduced systematic differences among the datasets. Boxplots showed that sample expression distributions became more consistent after correction. Principal component analysis indicated that samples clustered mainly according to dataset source before correction but became more intermixed after correction, indicating effec...

Discussion

An integrated and reproducible workflow was developed to identify PH/PAH-associated molecular biomarkers and candidate therapeutic compounds by combining public transcriptomics, weighted gene co-expression network analysis, functional enrichment, protein-protein interaction network analysis, three machine-learning algorithms, external validation, single-cell transcriptomic interpretation, quantitative reverse transcription PCR confirmation, Connectivity Map screening, and molecular docking. CXCL10, JUN, IFIH1, MX1, and T...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
2× SYBR Green PCR Mastermix	Beijing Solarbio Science & Technology Co., Ltd.	Catalog no. SR1110	Dye-based quantitative real-time PCR amplification and fluorescence detection
AI21.msvmRFE.R and e1071	Custom R script with the CRAN e1071 package	AI21.msvmRFE.R; e1071 v1.7-17	Support vector machine-recursive feature elimination
Agarose	Beijing Solarbio Science & Technology Co., Ltd.	Catalog no. A8201; CAS 9012-36-6	Assessment of total RNA integrity by agarose gel electrophoresis
Agarose gel electrophoresis apparatus	Beijing Liuyi Biotechnology Co., Ltd.	Model DYCZ-24DN	Electrophoretic assessment of RNA integrity
AutoDock Vina scoring engine	Center for Computational Structural Biology, Scripps Research	v1.2.0; RRID: SCR_011958	Protein-ligand pose scoring within the CB-Dock2 workflow
CB-Dock2	Cao Laboratory, CB-Dock2 web server	v2.0; accessed July 2026	Blind cavity detection and molecular docking of VX-745 with the selected protein structures
CellChat	CellChat R package	v2.1.2	Inference and visualization of cell-cell communication from the single-cell expression matrix
CellChatDB.human	Distributed with the CellChat R package	CellChatDB.human; Secreted Signaling subset; minimum cell threshold = 10	Human ligand-receptor interaction database for CellChat
clusterProfiler	Bioconductor R package	v4.20.0; Bioconductor release 3.23	Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses
Connectivity Map (CMap/CLUE)	Broad Institute	L1000/CLUE resource; RRID: SCR_016204; accessed July 2026	Computational drug-repositioning analysis
Custom oligonucleotide primers	Beijing Solarbio Science & Technology Co., Ltd.	Custom synthesized; primer sequences provided in Table 2	Amplification of ACTB, CXCL10, JUN, IFIH1, MX1, and TLR7
cytoHubba	Cytoscape App Store	v0.1	Degree-based hub-gene ranking in the protein-protein interaction network
Cytoscape	Cytoscape Consortium	v3.10.4; RRID: SCR_003032	Protein-protein interaction network visualization and analysis
DrugBank	DrugBank Knowledgebase	v6.0; RRID: SCR_002700	Compound identity and pharmacological-information curation
Gel documentation system	Beijing Liuyi Biotechnology Co., Ltd.	Model WO-9413B	Visualization and recording of agarose-gel RNA integrity results
Gene Expression Omnibus (GEO)	National Center for Biotechnology Information	GSE22356, GSE33463, GSE48149, GSE117261, and GSE210248; RRID: SCR_005012	Retrieval of bulk and single-cell transcriptomic datasets
GEOquery	Bioconductor R package	v2.80.0; Bioconductor release 3.23	Programmatic downloading and import of GEO expression and phenotype data
glmnet	CRAN R package	v5.0	Least absolute shrinkage and selection operator logistic regression and regularized logistic modeling
limma	Bioconductor R package	v3.68.0; Bioconductor release 3.23; RRID: SCR_010943	Differential expression analysis and empirical Bayes statistics
NanoDrop spectrophotometer	Thermo Fisher Scientific	NanoDrop ND-1000; software v3.8	Measurement of RNA concentration and A260/280 and A260/230 purity ratios
NCBI Primer-BLAST	National Center for Biotechnology Information	Web tool; RRID: SCR_003095; accessed July 2026	Verification of primer specificity
pROC	CRAN R package	v1.19.0.1; RRID: SCR_024286	Receiver operating characteristic analysis, DeLong confidence intervals, and Youden-index cutoffs
Protein Data Bank (PDB)	RCSB Protein Data Bank	CXCL10: 1LV9; JUN: 1JUN; IFIH1: 3B6E; MX1: 5GTM; TLR7: 7CYN; MAPK14/p38α: 1OUK; RRID: SCR_012820	Retrieval of experimentally determined protein structures for molecular docking
PubChem	National Center for Biotechnology Information	PubChem CID 3038525; RRID: SCR_004284	Retrieval of the three-dimensional structure and chemical identifiers of BRD-K91900765/VX-745
R	R Foundation for Statistical Computing	v4.6.1; RRID: SCR_001905	Statistical computing, data processing, machine learning, and visualization
randomForest	CRAN R package	v4.7-1.2	Random forest feature selection and variable-importance ranking
Real-time PCR system	Stratagene, now Agilent Technologies	Mx3000P Real-Time PCR System	qRT-PCR amplification, fluorescence acquisition, melting-curve analysis, and Ct export
Seurat	CRAN R package; Satija Laboratory	v5.5.1; RRID: SCR_016341	Single-cell RNA-sequencing quality control, normalization, dimensionality reduction, clustering, and annotation
STRING	STRING Consortium	v12.0; RRID: SCR_005223	Protein-protein interaction network construction
sva (ComBat)	Bioconductor R package	v3.60.0; Bioconductor release 3.23	Correction of inter-dataset batch effects
SwissADME	Swiss Institute of Bioinformatics	Web server; accessed July 2026	Drug-likeness, physicochemical-property, and ADME prescreening
Total RNA Extraction Kit	Beijing Solarbio Science & Technology Co., Ltd.	Catalog no. R1200	Extraction and purification of total RNA from lung-tissue samples
Universal RT-PCR Kit (AMV)	Beijing Solarbio Science & Technology Co., Ltd.	Catalog no. RP1200	Reverse transcription of total RNA into complementary DNA
WGCNA	CRAN R package	v1.74	Weighted gene co-expression network construction and module-trait analysis