The peripheral blood mononuclear cell (PBMC) samples utilized in this retrospective study were obtained from residual or discarded clinical specimens collected during routine diagnostic and therapeutic procedures from the General Hospital of Xinjiang Military Command with approval from the Institutional Review Board (2020RR0618). Written informed consent was obtained from all participants. Comprehensive clinical characteristics were presented in Supplementary Table S1. The single-cell RNA sequencing (scRNA-seq) cohort comprised 5 patients with high-altitude pulmonary hypertension (HAPH) and 5 matched healthy controls. We employed the SCOPIT (V1.1.4) to perform power analysis to determine the number of cells required per sample (Supplementary Figure S1; Supplementary Table S2), following established methodologies (13, 14). For bulk RNA sequencing, we assembled a composite cohort of 56 cases and 70 controls by integrating in-house samples with other group available dataset, which served as independent validation cohorts. The training and validation datasets were sourced from distinct institutions, with no overlap between them. The clinical baseline characteristics of the train and validation cohorts were summarized in Supplementary Table S3. Proteomic profiling was performed on PBMC samples from 18 HAPH patients and 24 healthy controls.

The PBMC samples were processed following an optimized protocol derived from established methodology (15). Cell viability was quantified (>90%) before loading suspensions onto a Chromium Single-Cell Controller (10 × Genomics) for library preparation. Sequencing reads were aligned to the GRCh38 reference genome using CellRanger (v8.0.0). Downstream analyses incorporated SCTransform normalization and Harmony integration via the Seurat (v4.0.2) (16) in R (v4.3.3) (17). Quality control retained cells expressing 500–4,000 genes with mitochondrial gene content <15%. Cell-type identification utilized canonical marker genes. Pseudobulk differential expression analysis was conducted using glmGamPoi (v1.12.2) (18), with significance criteria of |log2 (fold change)| > 0.5 and adjusted p < 0.05.

We employed the observed to predicted cell number (Ro/e) ratio and contribution scoring to identify hub immune cell subsets in HAPH (19). To quantify cluster-specific group preferences, we calculated the Ro/e ratio for each cell cluster across experimental groups using a validated methodology. Predicted cell numbers for cluster-group combinations were derived through chi-square testing. Contribution scores of cellular subgroups were computed following an established computational framework (20). Specifically, signature genes for each subset were defined as the top 100 differentially expressed genes (DEGs) between HAPH and control groups. These scores integrate both quantitative changes in cell populations and expression-level alterations of signature genes during disease progression. Final cluster contributions were determined by averaging the fold-change scores of all signature genes within each cellular subgroup.

To elucidate molecular mechanisms underlying HAPH progression, we conducted pseudotime trajectory reconstruction of myeloid cell differentiation using Monocle2 (v2.18.0) (21). The differentGeneTest function identified differentially expressed genes (DEGs) correlated with pathological transition from normal vasculature to HAPH. We applied the discriminative dimensionality reduction via learning a tree (DDRTree) algorithm for trajectory inference, followed by cellular ordering in reduced dimensional space. Pseudotemporal analysis revealed temporally upregulated HAPH-progression genes (adjusted p < 0.05), implicating their potential mechanistic contributions to disease advancement.

To identify genes associated with HAPH, we conducted bulk RNA-seq analysis PBMC samples from 56 HAPH patients and 70 healthy controls. Total RNA extraction was performed using the TRNzol Total RNA Extraction Reagent Kit. Samples meeting quality control thresholds [RNA integrity number (RIN) > 7.0 and 28S/18S ratio ≥ 0.7] were selected for library preparation. Sequencing was conducted on the Illumina HiSeq PE150 platform with 150-bp paired-end reads. Raw sequencing reads were aligned to the GRCh38 human reference genome using STAR software (v2.7.2a) (22), with subsequent gene-level quantification performed through htseq-count (v2.05) (23). The ComBat algorithm was used to eliminate batch effects. Differential gene expression analysis was conducted using DESeq2 (v1.40.2) (24), with statistically significant genes identified using thresholds of |log2 (fold change) | > 0.5 and adjusted p < 0.05 (25, 26).

To identify potential biomarkers associated with HAPH, we performed comparative proteomic analysis of PBMC samples obtained from 18 HAPH patients and 24 healthy controls. Cellular proteins were extracted using a standardized lysis buffer according to established protocols. Subsequent proteomic profiling was conducted using liquid chromatography–tandem mass spectrometry (LC–MS/MS). Raw data from data-dependent acquisition (DDA) experiments were analyzed with MaxQuant software (v1.5.3.30) (27) against the Human reference proteome (UniProtKB; 46,570 sequences). We generated a spectral library in Spectronaut with the following parameters: trypsin digestion; minimum peptide length of 7 amino acids; variable modifications including methionine oxidation and N-terminal acetylation; fixed carbamidomethylation of cysteine residues; and a peptide-spectrum match (PSM) false discovery rate (FDR) ≤ 1%. Remaining parameters retained default configurations. For data-independent acquisition (DIA) analysis, Spectronaut executed spectral deconvolution using the preconstructed library and implemented quality control through the mProphet algorithm, yielding high-confidence quantitative profiles. MSstats software (v4.8.7) was employed to reduce batch effects (28). Differential protein expression between HAPH and control groups was statistically evaluated using MSstats (v4.14.2) (29), with significance thresholds established at p < 0.05 and |log2 (fold change)| ≥ 0.5. The overlapping candidate genes identified from scRNA-seq, bulk RNA-seq, and proteomic data were designated as HAPH signature genes.

To establish a consensus HAPH prognostic signature exhibiting robust predictive accuracy and stability, we systematically integrated 12 machine learning algorithms encompassing 113 methodological combinations (30, 31). The signature development protocol consisted of four key phases: (a) systematic application of 113 algorithm combinations to construct predictive models using a leave-one-out cross-validation (LOOCV) framework; (b) comprehensive validation of all models in the independent dataset; and (c) quantitative evaluation using Harrell’s concordance index (C-index), with model selection based on maximal average C-index performance across validation cohorts; (d) The Hosmer–Lemeshow test and Brier score were computed to assess the calibration of the model.

To validate the HAPH signature, we acquired distinct tissue samples from the original sequencing cohort for Quantitative real-time PCR (qPCR) validation. Total RNA was isolated using the RNA Extraction Kit, followed by reverse transcription into cDNA using the cDNA Synthesis Kit. qPCR analysis was performed in six technical replicates using SYBR Green Master Mix on a CFX Connect Real-Time PCR Detection System (Bio-Rad Laboratories, United States). Gene-specific primers (Supplementary Table S4) were used in 20 μL reaction volumes with standardized cycling parameters: initial denaturation at 95°C for 3 min, followed by 40 cycles of 95°C for 10 s and 60°C for 30 s. The GAPDH was employed as an endogenous control for normalization of mRNA expression levels. Relative quantification was calculated using the comparative threshold cycle (2^{^−ΔΔCt}) method (26).

All statistical analyses were conducted using R (v4.3.3). Prior to parametric testing, data were assessed for normality using the Shapiro–Wilk test and homogeneity of variance using Levene’s test. For group comparisons meeting these assumptions, unpaired two-tailed Student’s t-tests were employed. Statistical significance was defined as p < 0.05 for all analyses.

Our single-cell transcriptomic atlas comprising 56,058 single cells from 5 HAPH patients and 5 healthy controls (Cntrl) (Figure 1A) systematically mapped disease-associated immune microenvironment. Following SCTransform normalization and Harmony integration to mitigate batch effects, unsupervised clustering of harmonized data identified 15 distinct clusters (Figure 1B). Canonical marker expression analysis (Figure 1C) resolved five hematopoietic lineages: T lymphocytes (CD3D+), myeloid cells (LYZ+), natural killer cells (NKG7+), B lymphocytes (MS4A1+), and platelets (PPBP+). The observed to predicted cell number (Ro/e) ratio revealed significant myeloid compartment expansion in HAPH (Ro/e = 1.26) with concomitant NK cell depletion (Ro/e = 0.89) compared to Cntrl (Figure 1D). Applying a computational framework to quantify cellular pathogenic potential, we identified myeloid subpopulations as dominant contributors (contribution score = 2.86), implicating their functional centrality in HAPH pathophysiology.

Pseudotime trajectory analysis revealed three transcriptionally distinct states within myeloid subpopulations. The HAPH and Cntrl groups displayed divergent distribution patterns along this trajectory, with HAPH samples predominantly localized to states 1–2 and control samples concentrated in state 3 (Figure 2A). This compartmentalization correlated with state-specific transcriptional signatures. Notably, pseudotime-dependent expression profiling identified 2,615 upregulated genes associated with HAPH progression (Figure 2A; Supplementary Table S5). KEGG pathway enrichment analysis demonstrated revealed these genes were significantly enriched in immune response pathways and erythrocyte differentiation processes (Figure 2B).

Differential expression analysis applied stringent thresholds (|log2 (fold change)| > 0.5, adjusted p < 0.05) to identify 144 significant differentially expressed genes (DEGs) using the bulk RNA-seq data (Figure 3A; Supplementary Table S6). Parallel data-independent acquisition (DIA) proteomic profiling quantified 420 and 456 plasma proteins in the HAPH and control groups, respectively. Comparative proteomic analysis revealed 77 differentially expressed proteins with |log2 (fold change)| > 0.5 and adjusted p < 0.05 (Figure 3B; Supplementary Table S7). Integrative multi-omics analysis demonstrated 22 consensus molecules consistently identified across scRNA-seq, bulk RNA-seq, and proteomic datasets. These overlapping biomolecules were subsequently designated as HAPH-associated signature genes (Figure 3C).

The 22 candidate genes identified through preliminary screening were subjected to integrated machine learning analysis to establish a diagnostic model for HAPH. Our in-house bulk RNA-seq cohort (n = 55) was utilized as the training set, while the external dataset (n = 71) served as independent validation set. Using a leave-one-out cross-validation (LOOCV) framework, we systematically evaluated 113 prediction models derived from 10 machine learning algorithms. Model performance was assessed through concordance index (C-index) evaluations across validation cohorts (Figure 4A). The random forest (RF) algorithm exhibited superior diagnostic performance, attaining the highest mean C-index (0.884). This optimized model demonstrated AUC values of 0.995 (95%CI 0.980–1.000) in the training cohort (Figure 4B) and 0.773 (95%CI 0.643–0.877) in the combined validation cohorts (Figure 4C). The Brier scores for both the training and validation sets are below 0.25, and the Hosmer-Lemeshow test p-value exceeds 0.05 (Supplementary Table S8), Moreover, the confusion matrix reveals accuracy rates of 0.964 (training set) and 0.704 (validation set), with both values exceeding 0.7 (Supplementary Figure S3). The final RF model incorporated six biomarker genes: HEMGN, HBG2, MYL9, ANK1, UBE2O, and RBPMS2. qPCR validation experiments confirmed significant upregulation of these biomarkers in HAPH patients relative to healthy controls (p < 0.05, Student’s t-test) (Supplementary Figure S1).