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 effective batch-effect reduction (Figure 1).

Using thresholds of an absolute log2 fold change >0.585 and an adjusted P value < 0.05, 78 differentially expressed genes were identified, including 44 upregulated and 34 downregulated genes (Figure 2A). The heatmap showed that interferon-related genes, including XAF1, MX1, IFI44L, EPSTI1, PARP9, IFIH1, CXCL10, GBP1, STAT1, SAMHD1, TNFSF10, and TLR7, were generally upregulated in PH-related samples. In contrast, erythroid-related genes, including HBG1, HBD, ALAS2, CA1, and SLC4A1, tended to be downregulated (Figure 2B).

Figure 1: Batch-effect correction. (A) Boxplots of the merged expression matrix before and after ComBat correction. (B) Principal component analysis plots showing sample distributions before and after batch-effect correction. Please click here to view a larger version of this figure.

Figure 2: Differentially expressed genes. (A) Volcano plot showing upregulated and downregulated genes in PH-related samples. (B) Heatmap showing differentially expressed genes between the control and disease groups. Please click here to view a larger version of this figure.

Weighted gene co-expression network construction and functional enrichment

Sample clustering showed stable overall clustering with no obvious outliers (Figure 3A). The scale-free topology fit index approached 0.8 at a power of 9, and β = 9 was selected for network construction (Figure 3B). Gene clustering and dynamic module identification yielded multiple co-expression modules (Figure 3C). The blue module showed the strongest association with PH status (r = 0.55, P = 1 × 10⁻¹⁹), whereas the turquoise and grey modules also showed correlations with PH (Figure 3D).

Consensus genes obtained by intersecting the key weighted gene co-expression network analysis module genes with the differentially expressed genes were enriched in antiviral immune defense, NF-κB and JAK-STAT regulation, inflammatory-factor responses, cytokine and chemokine receptor binding, and transcriptional regulation (Figure 4A). Kyoto Encyclopedia of Genes and Genomes enrichment analysis identified cytokine-cytokine receptor interaction, chemokine signaling, NOD-like receptor signaling, Toll-like receptor signaling, tumor necrosis factor signaling, and interleukin-17 signaling (Figure 4B), supporting immune-inflammatory dysregulation as a molecular basis of pulmonary vascular remodeling.

Figure 3: Weighted gene co-expression network analysis. (A) Sample clustering tree and trait heatmap. (B) Soft-threshold selection plot. (C) Gene dendrogram and module colors. (D) Module-trait relationship heatmap. Please click here to view a larger version of this figure.

Figure 4: Functional enrichment analysis. (A) Gene Ontology enrichment results for the consensus genes. (B) Kyoto Encyclopedia of Genes and Genomes pathway enrichment results for the consensus genes. Please click here to view a larger version of this figure.

Protein-protein interaction network and hub-gene screening

The STRING protein-protein interaction network constructed from the consensus genes revealed an interconnected immune-inflammatory network (Figure 5A). Degree-based CytoHubba ranking showed that FN1, CD44, JUN, TGFB1, CXCL8, and BCL2 had high connectivity (Figure 5B). These hub genes may participate in inflammatory signaling, cell adhesion, extracellular matrix remodeling, and pulmonary vascular structural remodeling.

Figure 5: Protein-protein interaction network and hub genes. (A) STRING protein-protein interaction network of the consensus genes. (B) Degree-ranked hub genes identified using CytoHubba. Please click here to view a larger version of this figure.

Machine-learning feature selection and diagnostic performance

Least absolute shrinkage and selection operator regression identified six candidate genes with non-zero coefficients after cross-validation (Figure 6A, B). Support vector machine-recursive feature elimination retained eight genes and yielded a cross-validation accuracy of 0.883 and an error of 0.117 (Figure 7A). The random forest out-of-bag error stabilized when the number of trees was ≥100, and IFIH1, JUN, and TLR7 ranked among the top genes according to the importance score (Figure 7B).

The intersection of the least absolute shrinkage and selection operator, support vector machine-recursive feature elimination, and random forest results identified five core genes: CXCL10, JUN, IFIH1, MX1, and TLR7 (Figure 8A). Single-gene receiver operating characteristic analysis showed moderate-to-good diagnostic discrimination, with area under the curve values of 0.842 for IFIH1, 0.833 for JUN, 0.827 for TLR7, 0.814 for CXCL10, and 0.759 for MX1 (Figure 8B).

Figure 6: Least absolute shrinkage and selection operator regression analysis. (A) Coefficient path generated by the least absolute shrinkage and selection operator regression. (B) Cross-validation error plot. Please click here to view a larger version of this figure.

Figure 7: Support vector machine-recursive feature elimination and random forest analyses. (A) Support vector machine-recursive feature elimination feature-selection plot. (B) Random forest model and gene-importance ranking. Please click here to view a larger version of this figure.

Figure 8: Machine-learning summary. (A) Venn diagram showing the intersection of the three feature-selection algorithms. (B) Receiver operating characteristic curves for the five core genes. Please click here to view a larger version of this figure.

External validation in GSE117261

GSE117261 was used as an independent lung-tissue validation cohort and was not included in differential-expression screening, weighted gene co-expression network construction, or machine-learning feature selection (Table 3). The completed validation results showed heterogeneous replication across the five genes (Table 4). CXCL10 was increased (log₂ fold change = 0.677; P = 0.0410; FDR = 0.144) and yielded an AUC of 0.639 (95% CI, 0.491–0.786), with a cutoff of 6.245, sensitivity of 0.724, and specificity of 0.600. JUN was increased (log₂ fold change = 0.463; P = 0.00248; FDR = 0.0194) and yielded an AUC of 0.714 (95% CI, 0.593–0.835), with a cutoff of 8.708, sensitivity of 0.707, and specificity of 0.720.

IFIH1 (log2 fold change = 0.107; P = 0.371; FDR = 0.591; AUC = 0.543, 95% CI, 0.398–0.687), MX1 (log2 fold change = 0.109; P = 0.488; FDR = 0.690; AUC = 0.475, 95% CI, 0.337–0.614), and TLR7 (log2 fold change = -0.050; P = 0.543; FDR = 0.733; AUC = 0.546, 95% CI, 0.414–0.679) did not meet the prespecified external-support criteria. TLR7 also showed a direction opposite to the qRT-PCR result. The exploratory five-gene model fitted and evaluated within GSE117261 yielded an apparent AUC of 0.740, whereas repeated nested cross-validation yielded an AUC of 0.656. Thus, JUN received the strongest independent support, CXCL10 showed limited directionally consistent evidence, and replication of IFIH1, MX1, and TLR7 was weak or discordant.

Item	Description
Dataset accession	GSE117261
Data source	Gene Expression Omnibus (GEO)
Sample type	Human lung-tissue transcriptomic microarray data
Sample size	58 PAH samples and 25 failed-donor control samples
Platform	GPL6244 / Affymetrix Human Gene 1.0 ST Array
Validation objective	Expression differences, single-gene ROC analysis, and exploratory five-gene combined ROC modeling for CXCL10, JUN, IFIH1, MX1, and TLR7
Role in this study	Independent external-validation dataset; not included in the original training, WGCNA, or feature-selection analyses

Table 3: Basic information for the GSE117261 independent external-validation dataset. The table summarizes the dataset source, sample type, sample size, platform, validation objectives, and role of GSE117261 in the study.

Gene/model	PAH n	Control n	log2 fold change	P value	FDR	AUC	AUC 95% CI	Youden cutoff	Sensitivity	Specificity
CXCL10	58	25	0.677	0.041	0.144	0.639	0.491–0.786	6.245	0.724	0.6
JUN	58	25	0.463	0.00248	0.019	0.714	0.593–0.835	8.708	0.707	0.72
IFIH1	58	25	0.107	0.371	0.591	0.543	0.398–0.687	5.3	0.759	0.44
MX1	58	25	0.109	0.488	0.69	0.475	0.337–0.614	6.833	0.724	0.4
TLR7	58	25	-0.05	0.543	0.733	0.546	0.414–0.679	4.065	0.138	1
Five-gene model (apparent/in-sample)	58	25	/	/	/	0.74	0.621–0.859	0.613	0.845	0.6
Five-gene model (repeated nested CV)	58	25	/	/	/	0.656	0.517–0.795	0.655	0.845	0.52

Table 4: Actual expression-difference and ROC validation results for GSE117261. The table reports the PAH and control sample sizes, log₂ fold changes, P values, false-discovery-rate-adjusted values, AUCs, 95% confidence intervals, Youden-index cutoffs, sensitivities, specificities, and interpretations for the five genes and exploratory combined models.

Quantitative reverse transcription PCR validation

Quantitative reverse transcription PCR validation included 20 biologically independent PAH lung-tissue samples and 20 biologically independent control samples, with three technical replicates averaged for each biological sample. CXCL10, JUN, IFIH1, MX1, and TLR7 were significantly upregulated in PAH (Table 5; Figure 9A). The mean relative-expression values were approximately 3.470 for CXCL10, 2.560 for JUN, 2.760 for IFIH1, 2.650 for MX1, and 2.580 for TLR7. The corresponding P values/FDR values were 1.43 × 10⁻⁷/7.15 × 10⁻⁷, 4.17 × 10⁻⁵/4.17 × 10⁻⁵, 1.10 × 10⁻⁵/1.38 × 10⁻⁵, 1.58 × 10⁻⁶/2.63 × 10⁻⁶, and 1.37 × 10⁻⁶/2.63 × 10⁻⁶, respectively.

Single-gene ROC analysis based on qRT-PCR expression showed AUCs of 0.988 for CXCL10 (95% CI, 0.961–1.000), 0.880 for JUN (95% CI, 0.758–1.000), 0.908 for IFIH1 (95% CI, 0.820–0.995), 0.945 for MX1 (95% CI, 0.874–1.000), and 0.948 for TLR7 (95% CI, 0.886–1.000) (Figure 9B; Table 5). The five-gene logistic regression model achieved an apparent AUC of 1.000 (DeLong 95% CI: 1.000–1.000), with sensitivity and specificity of 1.000 (Figure 9C). Expression-direction consistency between quantitative reverse transcription PCR validation and GSE117261 was visualized using a heatmap (Figure 9D). Because the same 40 samples were used for both fitting and evaluation, this reflected the apparent in-sample performance. In 100 repeats of stratified five-fold cross-validation using an L2-regularized logistic-regression model, the pooled out-of-fold AUC also remained 1.000 (95% CI, 1.000–1.000), and every repeat yielded an AUC of 1.000. Despite this internal stability, the cohort was small and independent prospective validation remains necessary. Directional comparison with GSE117261 showed concordant increases for CXCL10, JUN, IFIH1, and MX1, but discordant direction for TLR7 (Figure 9D).

Figure 9: Quantitative reverse transcription PCR validation. (A) Boxplots showing the relative expression of CXCL10, JUN, IFIH1, MX1, and TLR7 in 20 biologically independent PAH lung-tissue samples and 20 biologically independent control samples. Each biological sample was measured in three technical replicates, and the mean Ct value was used for analysis. For each boxplot, the center line represents the median, the box represents the interquartile range, the whiskers extend to 1.5 times the interquartile range, and individual points beyond the whiskers represent outliers. (B) Single-gene ROC curves based on qRT-PCR expression values; AUCs and DeLong 95% confidence intervals are shown. (C) ROC curves for the five-gene logistic-regression model, showing both apparent/in-sample performance and pooled out-of-fold performance from 100 repeated stratified five-fold cross-validation analyses. (D) Heatmap showing expression-direction consistency between qRT-PCR and GSE117261; CXCL10, JUN, IFIH1, and MX1 were concordantly increased, whereas TLR7 was discordant. Please click here to view a larger version of this figure.

Gene/model

PAH n

Control n

Control 2^-ΔΔCt, mean ± SD

PAH 2^-ΔΔCt, mean ± SD

Direction

P value

FDR

AUC

AUC 95% CI

Youden cutoff

Sensi- tivity

Specifi- city

Validation type

CXCL10

20

20

1.099 ± 0.502

3.470 ± 1.043

Upregulated

1.43E-07

7.15E-07

0.987

0.961–1.000

2.028

1

0.95

Single-gene qRT-PCR analysis

JUN

20

20

1.158 ± 0.738

2.560 ± 1.609

Upregulated

4.17E-05

4.17E-05

0.88

0.758–1.000

1.423

0.85

0.9

Single-gene qRT-PCR analysis

IFIH1

20

20

1.091 ± 0.440

2.760 ± 1.398

Upregulated

1.10E-05

1.38E-05

0.908

0.820–0.995

1.579

0.8

0.85

Single-gene qRT-PCR analysis

MX1

20

20

1.132 ± 0.582

2.650 ± 1.085

Upregulated

1.58E-06

2.63E-06

0.945

0.874–1.000

1.779

0.9

0.9

Single-gene qRT-PCR analysis

TLR7

20

20

1.080 ± 0.444

2.580 ± 1.284

Upregulated

1.37E-06

2.63E-06

0.947

0.886–1.000

1.764

0.85

0.9

Single-gene qRT-PCR analysis

Five-gene model (apparent/in-sample)

20

20

Not appli- cable

Not appli- cable

Not applicable

/

/

1

1.000–1.000

0.998

1

1

Same 40 biological samples used for model fitting and evaluation

Five-gene model (100× repeated 5-fold CV)

20

20

Not appli- cable

Not appli- cable

Not applicable

/

/

1

1.000–1.000

0.716

1

1

Internal cross-validation using L2-regularized logistic regression

Table 5: Complete qRT-PCR expression and ROC results for CXCL10, JUN, IFIH1, MX1, and TLR7, including the five-gene combined-model analyses. The table reports the PAH and control sample sizes, relative-expression values, expression directions, P values, false-discovery-rate-adjusted values, AUCs, 95% confidence intervals, Youden-index cutoffs, sensitivities, specificities, and validation types for the individual genes and combined models.

Single-cell transcriptomic validation in GSE210248

GSE210248 provided cell-level mechanistic support by showing that PAH pulmonary arterial remodeling was accompanied by altered communication between immune cells and vascular structural cells. This observation was consistent with the bulk transcriptomic enrichment of inflammatory responses, chemokine signaling, Toll-like receptor signaling, and tumor necrosis factor signaling.

Single-cell evidence suggested that the PAH pulmonary artery signaling network shifted toward structural cells, including smooth muscle cells and fibroblasts. Smooth muscle cells exhibited multiple states, including oxygen-sensing/pericyte-like, contractile, synthetic, and fibroblast-like. These results support a disease model in which immune-inflammatory activation and vascular structural-cell remodeling jointly drive PH/PAH progression.

Candidate compound screening and molecular docking

Connectivity Map screening identified BRD-K91900765 as the highest-ranked candidate compound among the top 10 hits, with a Logit score of 10.13 and a predicted probability of 0.085 (Figure 10). Compound curation showed that BRD-K91900765 corresponds to VX-745/neflamapimod, a selective p38α/MAPK14 inhibitor with a PubChem compound identifier of 3038525 and a molecular mass of 436.27 g/mol (Table 6).

Exploratory docking of BRD-K91900765/VX-745 with the five biomarker-associated proteins yielded top Vina scores of -7.5 kcal/mol for CXCL10, -7.6 kcal/mol for JUN, -7.5 kcal/mol for IFIH1, -8.7 kcal/mol for MX1, and -8.1 kcal/mol for TLR7 (Tables 7–11; Figure 11A–E). These results indicated predicted structural compatibility only and did not establish the five proteins as direct pharmacological targets. Preliminary absorption, distribution, metabolism, excretion, and toxicity predictions suggested that the compound had several drug-like properties, although the relatively high calculated cLogP requires further evaluation (Table 12). Docking against the established VX-745 target, MAPK14/p38α (PDB ID: 1OUK), was included as a positive reference analysis. The top MAPK14 cavity, C1, yielded a Vina score of -7.9 kcal/mol, a cavity volume of 3560 ų, a docking-box center of (2, 22, 34), and dimensions of (22, 31, 31) (Table 13; Figure 11F).

Figure 10: Connectivity Map candidate compound ranking. Ranking of the candidate compounds identified through Connectivity Map screening. BRD-K91900765 was the highest-ranked compound, with a Logit score of 10.13 and a predicted probability of 0.085. Please click here to view a larger version of this figure.

Figure 11: Three-dimensional molecular docking diagrams for BRD-K91900765/VX-745. (A) Exploratory docking with CXCL10. (B) Exploratory docking with JUN. (C) Exploratory docking with IFIH1. (D) Exploratory docking with MX1. (E) Exploratory docking with TLR7. (F) Positive-reference docking with the established pharmacological target MAPK14/p38α (PDB ID: 1OUK). Panels A–E indicate predicted structural compatibility and do not establish direct pharmacological targeting. Please click here to view a larger version of this figure.

Item	Description
CMap/Broad ID	BRD-K91900765 (common batch format: BRD-K91900765-001-xx-x)
Common name/aliases	VX-745; neflamapimod; VRT-031745; VD-31745
Chemical name	5-(2,6-dichlorophenyl)-2-(2,4-difluorophenyl)sulfanylpyrimido[1,6-b]pyridazin-6-one
PubChem CID	3038525
CAS number	209410-46-8
Molecular formula / relative molecular mass	C19H9Cl2F2N3OS; 436.27 g/mol
Canonical SMILES	C1=CC(=C(C(=C1)Cl)C2=C3C=CC(=NN3C=NC2=O)SC4=C(C=C(C=C4)F)F)Cl
InChIKey	VEPKQEUBKLEPRA-UHFFFAOYSA-N
Established major pharmacological target	MAPK14/p38α; p38β inhibition has also been reported with lower selectivity than p38α

Table 6: Chemical and pharmacological information for BRD-K91900765/VX-745. The table summarizes the compound identifiers, aliases, chemical name, molecular formula, molecular mass, structural descriptors, and established pharmacological target of BRD-K91900765/VX-745.

CurPocket ID	Vina score (kcal/mol)	Cavity volume (ų)	Center (x, y, z)	Docking size (x, y, z)
C1	-7.5	7568	49, 15, 4	34, 33, 35
C4	-7	160	46, 0, 4	22, 22, 22
C3	-6	208	34, 15, 9	22, 22, 22
C2	-5.8	465	45, -6, 19	22, 22, 22
C5	-5.1	150	63, 0, 20	22, 22, 22

Table 7: Predicted docking pockets for BRD-K91900765/VX-745 with CXCL10 (PDB ID: 1LV9). The table reports the ranked cavity identifiers, Vina scores, cavity volumes, docking-box centers, and docking-box dimensions generated by CB-Dock2.

CurPocket ID	Vina score (kcal/mol)	Cavity volume (ų)	Center (x, y, z)	Docking size (x, y, z)
C3	-7.6	1593	26, 35, 70	22, 22, 22
C2	-7.5	1927	25, 21, 60	35, 22, 22
C1	-7.4	5420	32, 37, 48	35, 22, 31
C5	-6	464	26, 5, 48	22, 22, 22
C4	-5.3	683	59, 32, 39	22, 22, 22

Table 8: Predicted docking pockets for BRD-K91900765/VX-745 with JUN (PDB ID: 1JUN). The table reports the ranked cavity identifiers, Vina scores, cavity volumes, docking-box centers, and docking-box dimensions generated by CB-Dock2.

CurPocket ID	Vina score (kcal/mol)	Cavity volume (ų)	Center (x, y, z)	Docking size (x, y, z)
C1	-7.5	583	15, 4, 17	22, 22, 22
C3	-7.2	146	19, 21, 24	22, 22, 22
C4	-6.6	141	31, 19, 15	22, 22, 22
C2	-6.2	273	38, 9, 24	22, 22, 22
C5	-6.2	125	27, -7, 10	22, 22, 22

Table 9: Predicted docking pockets for BRD-K91900765/VX-745 with IFIH1 (PDB ID: 3B6E). The table reports the ranked cavity identifiers, Vina scores, cavity volumes, docking-box centers, and docking-box dimensions generated by CB-Dock2.

CurPocket ID	Vina score (kcal/mol)	Cavity volume (ų)	Center (x, y, z)	Docking size (x, y, z)
C1	-8.7	523	-18, -14, -5	22, 22, 22
C4	-7.3	262	-13, -6, -9	22, 22, 22
C3	-7	306	12, 12, -7	22, 22, 22
C2	-6.9	408	-3, 1, -12	22, 22, 22
C5	-6.2	179	24, 25, 13	22, 22, 22

Table 10: Predicted docking pockets for BRD-K91900765/VX-745 with MX1 (PDB ID: 5GTM). The table reports the ranked cavity identifiers, Vina scores, cavity volumes, docking-box centers, and docking-box dimensions generated by CB-Dock2.

CurPocket ID	Vina score (kcal/mol)	Cavity volume (ų)	Center (x, y, z)	Docking size (x, y, z)
C2	-8.1	7347	112, 137, 148	33, 28, 35
C3	-8.1	2604	124, 124, 176	30, 22, 22
C1	-8	7676	137, 112, 148	34, 29, 35
C5	-6.9	1976	117, 157, 88	22, 22, 22
C4	-6.6	2040	132, 92, 91	22, 22, 22

Table 11: Predicted docking pockets for BRD-K91900765/VX-745 with TLR7 (PDB ID: 7CYN). The table reports the ranked cavity identifiers, Vina scores, cavity volumes, docking-box centers, and docking-box dimensions generated by CB-Dock2.

Category	Parameter	Result	Interpretation
Physicochemical property	Molecular weight	436.27 g/mol	Below 500 Da, meeting the Lipinski molecular-weight threshold
Physicochemical property	cLogP	Approximately 5.49	Slightly higher than 5, suggesting high lipophilicity and the need to consider solubility and nonspecific binding
Physicochemical property	TPSA	Approximately 47.26 Ų	Low polar surface area, consistent with potentially favorable membrane permeability
Drug-likeness	HBA/HBD	May-00	Meets the Lipinski thresholds for hydrogen-bond acceptors and donors
Drug-likeness	Rotatable bonds	3	Low conformational flexibility, favorable for stable binding conformations
Structural alerts	PAINS/Brenk alerts	Not detected	No common pan-assay interference or reactive structural alerts detected
Toxicity prediction	Ames mutagenicity	Predicted non-Ames toxic	Suggests a low predicted mutagenic risk; experimental validation is still required
Toxicity prediction	Carcinogenicity	Predicted non-carcinogenic	Suggests a relatively low predicted long-term carcinogenic risk; experimental validation is still required
Pharmacokinetic note	Oral availability/brain penetrance	Literature and databases indicate an orally available, brain-penetrant small molecule	Consistent with its development background as a p38α inhibitor; reevaluation is still needed for PH indications

Table 12: Preliminary physicochemical, drug-likeness, ADMET, and toxicity predictions for BRD-K91900765/VX-745. The table summarizes predicted physicochemical properties, drug-likeness measures, structural alerts, toxicity endpoints, and pharmacokinetic characteristics. These computational predictions are preliminary and do not replace experimental pharmacokinetic or toxicological validation.

CurPocket ID	Vina score (kcal/mol)	Cavity volume (ų)	Center (x, y, z)	Docking size (x, y, z)
C1	-7.9	3560	2, 22, 34	22, 31, 31
C5	-7.5	254	-9, 28, 60	22, 22, 22
C3	-7.4	351	12, 7, 37	22, 22, 22
C2	-6.3	827	18, 5, 28	22, 22, 22
C4	-6.3	334	-16, 17, 38	22, 22, 22

Table 13: Predicted docking pockets for BRD-K91900765/VX-745 with its established pharmacological target MAPK14/p38α (PDB ID: 1OUK), included as the positive-reference analysis. The table reports the ranked cavity identifiers, Vina scores, cavity volumes, docking-box centers, and docking-box dimensions generated using the same docking workflow applied to the five biomarker-associated proteins.

Collectively, the discovery analyses and qRT-PCR results support CXCL10, JUN, IFIH1, MX1, and TLR7 as candidate PH/PAH biomarkers associated with immune-inflammatory dysregulation and pulmonary vascular remodeling, although independent GSE117261 replication was strongest for JUN and variable for the other genes. BRD-K91900765/VX-745 is a computationally prioritized drug-repositioning candidate with a plausible mechanism of MAPK14/p38α inhibition; target-binding, cellular, pharmacokinetic, toxicity, and animal-model validation are required before therapeutic interpretation.

DATA AVAILABILITY:

All public transcriptomic datasets used in this study are available from the Gene Expression Omnibus database under accession numbers GSE22356, GSE33463, GSE48149, GSE117261, and GSE210248. All coding files, processed datasets, de-identified qRT-PCR raw and analyzed data, model outputs, and molecular docking input/output files have been consolidated into a structured Zenodo repository. The repository includes a README that describes each file, software and package versions, script execution order, and complete reproduction steps - https://zenodo.org/records/21682282

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 TLR7 were consistently prioritized as core feature genes and collectively mapped to an immune-inflammatory and interferon-related molecular axis. These findings support the concept that PH/PAH is not only a hemodynamic disorder but also a complex vascular remodeling disease involving immune activation, inflammatory signaling, innate nucleic acid sensing, and structural and cellular phenotypic changes^2,5,6.

The diagnostic potential of the five genes was supported by multi-algorithm feature selection and discovery-cohort ROC analysis. Independent validation in GSE117261 was heterogeneous rather than uniform: JUN met the prespecified FDR and AUC criteria, CXCL10 showed a directionally consistent nominal increase without transcriptome-wide FDR significance, IFIH1 and MX1 showed limited replication, and TLR7 showed a discordant direction. These results do not support the claim that all five genes were independently validated and indicate possible effects of cohort composition, tissue heterogeneity, platform differences, and disease severity. In contrast, qRT-PCR in 20 PAH and 20 control lung-tissue samples confirmed significant upregulation of all five genes and favorable single-gene ROC performance.

The five-gene qRT-PCR logistic model achieved an apparent AUC of 1.000 (95% CI, 1.000–1.000), and its pooled out-of-fold AUC remained 1.000 in 100 repeated stratified five-fold cross-validation analyses. Nevertheless, the model was developed in only 40 biological samples, and complete separation in a small retrospective cohort can produce optimistic and unstable performance estimates. The panel should therefore be considered an exploratory molecular signature rather than a clinically validated diagnostic tool. Larger multicenter cohorts, prespecified fixed model coefficients, protein-level validation, immunohistochemistry, and prospective testing are required before clinical translation.

Among the five core genes, CXCL10 may promote immune cell recruitment and local inflammatory amplification within the pulmonary vascular microenvironment. IFIH1 and TLR7 are involved in innate nucleic acid sensing and may reflect activation of antiviral-like inflammatory pathways. MX1 is a classical interferon-stimulated gene and may represent a downstream marker of type I interferon pathway activation. JUN is a stress-responsive transcription factor that links inflammatory stimulation to cell proliferation, apoptosis, and tissue remodeling. Together, these genes suggest a biologically coherent model in which innate immune activation and interferon-related signaling interact with vascular remodeling processes in PH/PAH. This interpretation is consistent with previous evidence that inflammation, immunity, and interferon-related pathways contribute to PAH pathobiology^2,5,6.

Single-cell validation provided a mechanistic context for the bulk-derived findings. GSE210248 suggested that PAH pulmonary arterial remodeling was accompanied by altered communication between immune cells and vascular structural cells, including smooth muscle cells, fibroblasts, endothelial cells, and monocytes/macrophages. The presence of multiple smooth muscle cell phenotypic states, including contractile, synthetic, oxygen-sensing/pericyte-like, and fibroblast-like, supports a disease model in which immune activation and structural remodeling of cells occur concurrently. This cellular evidence is important because bulk transcriptomic signals may arise from altered cell proportions, immune cell infiltration, or transcriptional changes in resident vascular cells. The single-cell analysis, therefore, places CXCL10, JUN, IFIH1, MX1, and TLR7 within a multicellular pulmonary vascular remodeling ecosystem rather than within a single-cell-type process^18,19,20,21.

The drug-repositioning analysis identified BRD-K91900765, corresponding to VX-745/neflamapimod, as the highest-ranked computational candidate. VX-745 is a selective p38α/MAPK14 inhibitor, and its relationship to inflammatory stress pathways makes it mechanistically plausible in the context of PH/PAH-associated inflammation³⁰. Docking against MAPK14/p38α was therefore included as a mechanistically relevant positive-reference analysis. By contrast, docking against CXCL10, JUN, IFIH1, MX1, and TLR7 was exploratory and indicated only predicted structural compatibility; it did not demonstrate that these biomarker-associated proteins are direct VX-745 targets or establish direct binding, target inhibition, or therapeutic efficacy. A more biologically plausible hypothesis is that VX-745 may indirectly modulate the identified immune-inflammatory and interferon-related transcriptional signature by inhibiting MAPK14. Connectivity Map predictions, docking scores, and ADMET estimates remain computational evidence. Further work should include biochemical target-binding assays, experiments in pulmonary arterial endothelial and smooth muscle cells, inflammatory stimulation models, pharmacokinetic and toxicological evaluations, and animal model validation.

Recent experimental studies in hypoxic pulmonary hypertension have also highlighted the importance of communication between neutrophils and pulmonary vascular cells. HCK-mediated interactions between neutrophils and pulmonary arterial smooth muscle cells and SERPINB3-mediated interactions between neutrophils and endothelial cells have been reported as contributors to pulmonary vascular remodeling³¹˒³². The SERPINB3–STAT1/3 axis is particularly relevant to the present findings because interferon-related, STAT1, and JAK–STAT signals were identified in the transcriptomic analyses. Together, these observations support the interpretation that immune-cell activation and communication with vascular structural cells may contribute to PH/PAH progression.

The study has several strengths. Multiple public datasets and batch-effect correction were used to reduce dataset-specific bias. Differential expression analysis, weighted gene co-expression network analysis, protein-protein interaction analysis, and three machine learning algorithms were combined to improve feature robustness. Independent bulk validation, qRT-PCR confirmation, and single-cell evidence provided complementary but nonidentical layers of evidence. The heterogeneous GSE117261 findings and the small qRT-PCR cohort also highlight important limitations, including incomplete external replication, potential tissue- and platform-specific effects, and the risk of overfitting. Biomarker discovery was also extended to candidate-compound screening, but the docking analyses remain hypothesis-generating. Future studies should validate the five genes in larger independent cohorts using spatial transcriptomics, proteomics, immunohistochemistry, and organoid or vascular-on-chip models. Overall, CXCL10, JUN, IFIH1, MX1, and TLR7 remain candidate PH/PAH biomarkers linked to immune-inflammatory and interferon-related vascular remodeling, while BRD-K91900765/VX-745 is a computational drug-repositioning candidate whose therapeutic relevance requires experimental validation.

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