$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
Overview of analysis workflow
The overall design of this study’s transcriptomic and machine learning-based analysis is illustrated in Figure 1, encompassing key steps: collection of Renin-Angiotensin System-related genes (RASRGs); screening of RAS-related differentially expressed genes (RASRDEGs) from hypertension datasets; functional enrichment analysis (GO/KEGG/GSEA); immune infiltration analysis (CIBERSORT); construction of protein-protein interaction (PPI) and regulatory networks; machine learning-based key gene selection (logistic regression, random forest [RF]); and evaluation of the hypertension diagnostic model. A complete list of software, databases, and online tools used in this study is provided in the Table of Materials.
Data download
Hypertension datasets GSE753608 and GSE74144 (Homo sapiens) were obtained via the R package GEOquery9 from the GEO database10. GSE75360 derived from peripheral blood mononuclear cells (platform: GPL10558) included 10 hypertension and 11 control samples; GSE74144 derived from white blood cells (platform: GPL13497) included 14 hypertension and 8 control samples (Table 1). Protein-coding RASRGs (1,264) were initially identified via GeneCards11 (keyword: "Renin-Angiotensin System") and PubMed (keyword: "Renin-Angiotensin System")12,13. The intersection of these RASRGs with genes in GSE75360/GSE74144 yielded 1,159 final RASRGs14. The two datasets were processed separately because they were generated on different microarray platforms. Probe annotation was performed according to the corresponding GPL platform annotation files, and normalized gene expression matrices were used for downstream analyses. Box plots were used to compare expression distributions before and after normalization.
Hypertension-related renin-angiotensin-related differentially expressed genes
Samples in the GSE75360 dataset were categorized into the hypertension group and control group. The limma software was employed to conduct differential gene expression analysis between the two groups14, with differentially expressed genes (DEGs) identified by the threshold of |logFC| > 0.45 and p-value < 0.05. The results of this differential analysis were visualized via volcano plots (generated using the R package ggplot2).
To obtain RASRDEGs, DEGs meeting the above threshold (|logFC| > 0.45, p-value < 0.05) were cross-referenced with RAS-related genes (RASRGs), and the intersection result was presented via a Venn diagram. Subsequently, the expression patterns of the identified RASRDEGs were visualized as a heatmap using the R package pheatmap, and the chromosomal localization of RASRDEGs was displayed via chromosome maps generated using the R package RCircos15.
Differentially expressed gene validation and ROC curve analysis
An intergroup plot was constructed to analyze RASRDEG expression differences between hypertension/control in GSE7536016, R package pROC was used to plot ROC curves and calculate AUC (0.5–0.7: low accuracy; 0.7–0.9: moderate; >0.9: high) for RASRDEG diagnostic efficacy.
Correlation analysis
Spearman’s correlation analysis was performed on RASRDEG expression in GSE75360; results were visualized via heatmap (R package ggplot2) (|r| < 0.3: no/weak correlation; 0.3–0.5: weak; 0.5–0.8: moderate; >0.8: strong).
Enrichment analysis of GO and KEGG
GO (Gene Ontology, 2024 release, http://geneontology.org/) is a widely used resource for large-scale functional enrichment, covering three domains: biological processes (BP), cellular components (CC), and molecular functions (MF)17. KEGG (Kyoto Encyclopedia of Genes and Genomes, Release 109.0, 2024, https://www.genome.jp/kegg/) stores data on genomes, biopathways, diseases, and drugs18.
RASRDEGs were subjected to GO annotation and KEGG pathway enrichment analysis using the R package clusterProfiler19. Enrichment test method: hypergeometric test; multiple test correction method: Benjamini-Hochberg (BH) method. Screening criterion: adjusted p-value < 0.05.
Gene set enrichment analysis (GSEA)
For cohort-level GSEA, all genes tested in the differential expression analysis of GSE75360 were ranked in descending order by logFC and used as the input gene list for clusterProfiler19. No DEG prefiltering was applied before GSEA. The c2 gene set collection from MSigDB20. Parameters: seed = 2022, 10–500 genes per set; screening criteria: corrected p < 0.05 (Benjamini-Hochberg, BH method), FDR < 0.2521.
Construction of hypertension diagnostic model
To identify key genes associated with hypertension, we employed two types of machine learning algorithms: logistic regression and random forests (RF). Logistic regression (binary dependent variable: hypertension/control) screened RASRDEGs with p < 0.05. Random Forest (RF, R package randomForest): parameters set.seed(520), ntree = 1000; MeanDecreaseGini (variable importance indicator) was extracted, and top15 RASRDEGs were selected. RASRDEGs were screened with a p value < 0.05 as the standard.
The RF (Random Forest) algorithm, an ensemble learning method under the Bagging category (integrating multiple decision trees), was applied via the R package randomForest22 (parameters: set.seed(520), ntree = 1000). MeanDecreaseGini (reflecting variable importance by average purity decrease during node splitting) of feature genes was extracted, and the top 15 RASRDEGs were selected. Finally, a Venn diagram of genes screened by logistic regression and RF was plotted to identify hypertension-related key genes.
Validation of hypertension diagnostic model
A logistic regression model was built based on key genes; linear predicted value (η) was calculated as:
The R package pROC16 was used to plot ROC curves and evaluate the model’s efficacy in predicting hypertension risk. A nomogram was constructed via the R package rms23 to visualize the contribution of each key gene to the logistic regression model (reflecting the association between key genes and hypertension risk). Calibration curves were generated to assess the consistency between predicted and actual hypertension probabilities; decision curve analysis (DCA, R package ggDCA24) was performed to evaluate the model’s clinical utility (net benefit) in GSE75360 and GSE74144.
Single-gene GSEA
GSEA explores the role of genes associated with a specific gene in biological processes/pathways/diseases by analyzing its expression, aiding in understanding the gene’s functional role. For each focal gene in GSE75360, samples were split at the median into high- and low-expression groups. Differential expression analysis was then performed across all tested genes, and genome-wide logFC values were ranked from highest to lowest before GSEA with clusterProfiler19. No DEG prefiltering was applied before GSEA. Parameters: seed = 2020, 10–500 genes per set (c2 gene set collection from MSigDB21). Screening criteria: p < 0.05 (adj. p corrected via BH method).
Immune infiltration analysis (CIBERSORT)
The CIBERSORT algorithm25 (based on linear support vector regression) deconvoluted the transcriptome matrix to estimate immune cell composition in mixed samples (data with immune cell enrichment score > 0 were selected). The final immune cell infiltration matrix of GSE75360 was visualized via a proportion bar chart. Spearman’s correlation was used to analyze immune cell-immune cell and key gene-immune cell associations, with results presented as a correlation heatmap (R package pheatmap) and correlation bubble plot (R package ggplot2), respectively.
Protein-protein interaction (PPI) network
PPI networks are systems of interconnected proteins regulating biological processes via interactions. Using the STRING database26, a PPI network for key genes was constructed (minimum interaction score: 0.150, low confidence). Renin-angiotensin-related hub genes were selected by screening interacting genes. The GeneMANIA database27, which identifies functionally similar genes using genomic and proteomic datasets, was used to predict functionally similar genes of key RAS genes and to construct a protein interaction network.
Construction of regulatory network
mRNA-TF network: Transcription factors (TFs) regulate gene expression via post-transcriptional interaction with target genes. TFs targeting hub genes and their regulatory relationships were retrieved from the ChIPBase database28, and the mRNA-TF network was visualized using Cytoscape29.
mRNA-miRNA network: miRNAs modulate multiple target genes (single targets may be co-regulated by multiple miRNAs). StarBase v3.030 was used to identify miRNAs associated with RASRDEGs, and the mRNA-miRNA network was visualized via Cytoscape.
mRNA-drug network: Toxicogenomic databases31 were used to predict direct/indirect drug targets of hub genes. The mRNA-drug network (showing gene-drug interactions) was visualized with Cytoscape to complete network construction.
Ang II-induced HUVEC model
Human umbilical vein endothelial cells (HUVECs) were maintained at 37°C in a humidified incubator with 5% CO2. Cells were maintained in complete endothelial cell culture medium supplemented with fetal bovine serum and antibiotics according to the supplier’s instructions. To establish an in vitro hypertension-related endothelial injury model, HUVECs were treated with angiotensin II (Ang II; 100 nM) for 48 h. Vehicle-treated cells were used as the control group.
For gene intervention experiments, small interfering RNAs targeting CST3 or FURIN (si-CST3 and si-FURIN), corresponding negative control siRNA (si-NC), CST3 or FURIN overexpression plasmids (oe-CST3 and oe-FURIN), and the corresponding empty-vector control (oe-NC) were transfected into HUVECs using a commercial transfection reagent according to the manufacturer’s protocol. After transfection, cells were exposed to Ang II and then harvested for expression validation and functional assays. Knockdown and overexpression efficiencies were confirmed by qRT-PCR and western blotting.
qRT-PCR
Total RNA was isolated from HUVECs with a standard RNA extraction reagent, and complementary DNA was generated using a reverse transcription kit. SYBR Green chemistry was used for qRT-PCR. Expression levels of LRP1, CTSD, MTHFR, AUTS2, FURIN, CST3, FCER1G, TBXAS1, IL-6, TNF-α, VCAM1, ICAM1, and eNOS were normalized to GAPDH. and calculated by the 2−ΔΔCt method.
Western blotting
For western blot analysis, proteins were extracted with RIPA lysis buffer and quantified using a BCA assay. Equal protein amounts were resolved by SDS-PAGE and transferred to PVDF membranes. After blocking, membranes were incubated with primary antibodies against CST3, FURIN, TBXAS1, or GAPDH and then with suitable secondary antibodies. Bands were detected by chemiluminescence, and densitometry was normalized to GAPDH. Secreted CST3 in culture supernatants was quantified with an ELISA kit following the manufacturer’s protocol.
Cell viability
Cell viability was assessed using the Cell Counting Kit-8 (CCK-8) assay. Briefly, transfected and Ang II-treated HUVECs were seeded into 96-well plates, and absorbance at 450 nm was measured at 0, 24, 48, and 72 h after addition of the CCK-8 reagent. Cell migration was evaluated using Transwell chambers. After the indicated interventions, cells were seeded into the upper chambers, and migrated cells on the lower membrane surface were fixed, stained, and counted under a microscope in randomly selected fields.
Inflammatory test
To evaluate inflammatory activation, oxidative stress, and endothelial function, IL-6, TNF-α, VCAM1, ICAM1, and eNOS .mRNA levels were detected by qRT-PCR. Nitric oxide (NO) levels in the culture supernatant were measured using a commercial NO assay kit, and intracellular reactive oxygen species (ROS) levels were detected using DCF fluorescence according to the manufacturer’s instructions.
Statistical analysis
Transcriptomic processing and modeling were performed in R. Continuous variables were assessed for normality with the Shapiro-Wilk test. For two-group comparisons, independent-samples t-tests were used for normally distributed variables, whereas Wilcoxon rank-sum tests were used for non-normal variables. For three or more groups, one-way analysis of variance with appropriate post hoc testing was used when normality and homogeneity of variance assumptions were met; otherwise, the Kruskal-Wallis test was applied. CCK-8 time-course data were analyzed using two-way analysis of variance. Spearman correlation coefficients were calculated for association analyses. Unless otherwise stated, experimental results are shown as mean ± SD, and two-tailed p < 0.05 was considered significant.