The underlying molecular mechanisms and biomarkers of Hip fracture combined with deep vein thrombosis based on self sequencing bioinformatics analysis

Abstract

Background

Thrombus formation is a severe complication in orthopedic surgery, significantly increasing mortality in patients with fractures. Therefore, identifying feature genes to determine thrombus presence in fracture surgeries is critical.

Methods

Whole blood samples were collected from 18 patients with fractures with thrombosis (YES_thrombus) and 18 patients with fractures without thrombosis (NO_thrombus) from the Second Hospital of Shanxi Medical University, China, and used for transcriptome sequencing and quality control to generate the YES_thrombus dataset. Candidate genes were identified by overlapping differentially expressed genes (DEGs) with key module genes from weighted gene co-expression network analysis (WGCNA). Functional enrichment analysis was then performed to explore the roles of the candidate genes. Feature genes were further refined by intersecting results from three machine learning algorithms and constructing an artificial neural network (ANN). Diagnostic performance was assessed using receiver operating characteristic (ROC) curves. Additionally, single-gene gene set enrichment analysis (GSEA) was conducted, and correlations between feature genes and differential immune cells were analyzed. The competing endogenouse RNA (ceRNA) regulatory network for feature genes was also constructed. Finally, quantitative reverse transcriptase PCR (qRT-PCR) was used to validate gene expression.

Results

Seven candidate genes were selected, with functional enrichment analysis linking them to the autophagosome and PPAR signaling pathways. Five feature genes with excellent diagnostic performance were identified. Single-gene GSEA enrichment showed that the feature genes were primarily associated with the cytosolic ribosome and oxidative phosphorylation. The correlation analysis revealed that aDC exhibited the strongest negative correlation with WDR81 and the strongest positive correlation with RGS1. The ceRNA regulatory network encompassed three feature genes, five miRNAs, and 236 lncRNAs. Expression analysis indicated that, with the exception of WDR81, other genes were significantly upregulated in the NO_thrombus group. qRT-PCR validation confirmed that the expression of AAED1, ARL4A, and WDR81 matched sequencing results.

Conclusions

In conclusion, five feature genes (RGS1, HSF2, ARL4A, AAED1, and WDR81) were identified, and functional enrichment analyses were conducted, providing a foundation for predicting the diagnosis of fractures associated with thrombosis.

Background

Hip fractures (HF), including subtrochanteric, femoral neck, and intertrochanteric fractures, are among the most common fractures in older adults with osteoporosis. Sing et al. reported that over 70% of HF cases occur in females, with more than 40% affecting individuals aged 85 years and older. Although less common than in women, HF in men tends to be more severe, as one-third of men who experience an HF die within a year [1, 2]. Most patients with HF require prolonged bed rest, leading to reduced physical activity, muscle atrophy, impaired blood circulation, increased red blood cell aggregation, elevated blood viscosity, and heightened clotting risk, ultimately increasing the likelihood of deep vein thrombosis (DVT) [3]. Preoperative DVT is a frequent complication in patients with HF [2], with its incidence reaching up to 50%. Abnormal coagulation in blood vessels leads to DVT, and thrombus detachment can result in pulmonary embolism (PE), a major contributor to sudden clinical death.

Additionally, HF combined with DVT has significant negative consequences, including disability, depression, cardiovascular disease, and, in cases of venous thromboembolism (VTE), a mortality rate of 10–30% [4]. As the global population ages, the incidence of HF in the elderly continues to rise, along with the occurrence and recurrence of DVT [5]. This situation incurs substantial healthcare costs due to extended hospital stays and rehabilitation. While independent risk factors for DVT and predictive factors for its recurrence have been identified, and routine anticoagulation therapy has been implemented in cases of single lower limb fractures, it remains unclear whether novel genes play a role in the etiology of DVT. Consequently, elucidating the biological processes and pathways linked to DVT is essential, as it may uncover potential drug targets for preventing deep DVT following HF.

Recent advancements in sequencing technologies, coupled with bioinformatics analysis, have enabled the identification of differentially expressed genes (DEGs) and the exploration of associated biological processes [6]. With the continuous development of bioinformatics methodologies, novel approaches to studying DVT have emerged. Using transcriptome data from fracture individuals with thrombosis, this study identified feature genes through PPI networks and machine learning algorithms, explored the biological pathways involved, and examined the immune characteristics of these genes in fracture-associated thrombosis. In vivo validation was conducted to assess the expression of feature genes in DVT. Through these analyses, characteristic genes that may predict unstable DVT and serve as promising therapeutic targets were identified for managing this condition.

Methods

Source of data

Inclusion and exclusion criteria

1) Inclusion criteria were as follows: Patients admitted from April 2023 to April 2024, aged ≥ 60 years, with HF resulting from osteoporosis combined with low-energy trauma (typically caused by falls), including femoral neck fractures and intertrochanteric fractures. DVT was confirmed via vascular ultrasound at three time points during hospitalization: within 24 h of admission, one day before surgery, and one week post-surgery. The diagnostic criteria for DVT via color Doppler ultrasound included solid echogenic filling in the blood vessels, lack of compression in the vessel lumen with the probe, and either bypassing or absence of blood flow signals. Two orthopedic doctors independently reviewed the diagnostic information for all patients with HF, while two ultrasound specialists and two orthopedic doctors independently reviewed the diagnostic information for HF individuals with DVT.

2) Exclusion criteria were as follows: (1) Pathological fractures; (2) Previous HF; (3) Femoral shaft fractures; (4) Distal femoral fractures; (5) Subtrochanteric fractures; (6) Complications and sequelae of HF (non-union, delayed healing, osteomyelitis, osteoarthritis, stiffness); (7) Complications related to prosthetic devices; (8) Femoral head necrosis; (9) Hip dislocation; (10) Removal of internal fixation devices. Additionally, patients with HF caused by high-energy trauma, typically seen in younger individuals, were excluded. Patients with DVT resulting from kidney disease, rheumatic disease, hematological disorders, immune system disorders, infections, tumors, or those who developed DVT prior to injury were also excluded.

This study was approved by the institutional review boards of The Second Hospital of Shanxi Medical University (protocol code (2023) YX No. (262), approval date: September 13, 2023), and written informed consent was obtained from all participants. Whole blood samples were collected from 18 fracture patients with thrombosis (YES_thrombus) and 18 fracture patients without thrombosis (NO_thrombus). All samples were utilized for transcriptome sequencing and are referred to as the YES_thrombus dataset.

RNA extraction and library construction for transcriptome sequencing

Total RNA was extracted from whole blood samples using the TRIzol reagent (Invitrogen, CA, USA) following the manufacturer’s instructions. RNA quantity and purity were measured using a NanoDrop ND-1000 spectrophotometer. The integrity of the RNA was verified by agarose gel electrophoresis. mRNA was selectively captured using Oligo (dT) magnetic beads (Dynabeads Oligo [dT], Thermo Fisher, USA). The fragmented RNA was reverse transcribed into complementary DNA (cDNA) using SuperScript™ II Reverse Transcriptase (Invitrogen, CA, USA). E. coli DNA polymerase I (NEB, USA) and RNase H (NEB, USA) were employed for second-strand synthesis. The resulting DNA-RNA hybrids were converted into DNA duplexes. The second strand was digested using UDG enzyme (NEB, MA, USA), followed by PCR amplification to create a library with a fragment size of 300 bp ± 50 bp. Sequencing was performed on the Illumina Novaseq 6000 platform (LC Bio Technology Co., Ltd., Hangzhou, China) following standard procedures for paired-end sequencing.

Processing of sequencing data

High-throughput transcriptome sequencing was conducted on 36 whole blood samples using the Illumina platform. Raw sequencing data were processed with Trimmomatic (v0.39) [7]and FastQC (v0.11.9) software [8] to remove adapter sequences and low-quality reads. The filtered data were aligned to the reference genome (HG37) using Hisat2 software (v2.2.1) [9]. Gene count values were obtained using the featureCounts software.

Identification of differentially expressed genes (DEGs)

A gene expression matrix was generated for each sample, and normalization was performed using the DESeq2 package (v1.26.0) [10]. Expression data were visualized with a box plot. DEGs between the YES_thrombus and NO_thrombus groups were identified using DESeq2 in the YES_thrombus dataset (adjusted P value < 0.05 and|log2FC| > 0.5) [10]. DEGs were visualized through a volcano plot and heatmap, which were generated using the ggplot2 (v3.3.2) [11] and pheatmap (v4.1.0) [12] packages, respectively.

Selection of key module genes by weighted gene co-expression network analysis (WGCNA)

To identify gene sets with high covariance, samples were clustered, and outliers were removed using the goodSamplesGenes function from the WGCNA (v1.71) package [13]. The optimal soft threshold was determined based on R² and mean connectivity. Modules were subsequently identified, and the module most strongly associated with thrombus was determined using dynamic tree-cutting. Genes from the key module were selected for further analysis.

Screening and functional enrichment of candidate genes

Candidate genes were identified by overlapping DEGs and key module genes. The expression of these candidate genes was visualized through a volcano plot and box plot. To investigate the functional roles of the candidate genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the clusterProfiler package (v4.0.2) [14] (P value < 0.05).

Machine learning screening of feature genes, construction of diagnostic models, and protein-protein interaction (PPI) network

To evaluate feature genes, three machine learning models were developed: the Least Absolute Shrinkage and Selection Operator (LASSO), Support Vector Machine-Recursive Feature Elimination (SVM-RFE) and Random Forest (RF). LASSO logistic regression was performed using the R software glmnet package (v4.0-2) [15], in which the parameters were set to famil = binomial and type.measure = class. then, the lambda.min values determined by cross-validation were used and the significant feature genes were screened. Important gene features were screened by recursive feature elimination (RFE) using the rfecontrol function in the caret package (v0.4.6) [16] in conjunction with the SVM and RF algorithms.SVM-RFE utilised the svmRFE function with input gene expression data and grouping information to remove the least important features through 20-fold cross-validation iterations and recorded model performance metrics to select the subset of features with the lowest error rate.RF-RFE, on the other hand, uses the randomForest function to construct the model, assesses the variable importance of each gene, iteratively removes genes with lower importance, and records the model performance, and ultimately selects the number of features with optimal performance. Genes identified as significant by all three machine learning algorithms were classified as feature genes. The diagnostic performance of the feature genes and models was assessed by generating receiver operating characteristic (ROC) curves using the pROC software (v1.18.0) [17]. Additionally, confusion matrix analysis was performed to evaluate model performance.

Finally, to explore the interactions between feature genes, a protein-protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interaction Genes/Proteins (STRING) website.

Construction of artificial neural network (ANN)

To validate the effectiveness of the biomarkers in distinguishing between fracture samples with thrombus and simple fracture samples, an ANN diagnostic model was constructed based on logistic regression principles using the R software package neuralnet (v1.44.2) [18]. In addition, ROC curves were plotted using the pROC software (v1.18.0) to assess the validity of the model [17]. Finally, we also used a k-fold cross-validation approach to comprehensively and objectively assess the efficacy of the model. In this process, we divided the data into multiple subsets, each containing balanced distributions of positive and negative samples, which were used for model training and testing, respectively. In each validation, the model’s accuracy, AUC, and detailed categorical performance in the confusion matrix were recorded.

Correlation analysis of clinical indicator

After obtaining the feature genes, a correlation analysis was performed to examine their relationship with clinical indicators. The Spearman correlation method was used (P value < 0.05 and|cor| > 0.3), and the results were visualized using correlation heatmaps and scatter plots.

Single-gene gene set enrichment analysis (single-gene GSEA) of feature genes

In order to further understand the abundant regulatory pathways and biological functions of the biomarkers, we calculated the correlation coefficients of the expression of all genes with the biomarkers individually using pearson, which was used as a ranking criterion. Subsequently, single-gene GSEA analysis (p-adjust < 0.05) was performed using the R software clusterProfiler package (v4.0.2) [14]. The GO database was used with version 05/03/2023 and the KEGG pathway database was used with KEGG version 101.0+ (12/07/2023). Based on the significance values, we selected the top 10 results for visualisation.

Ingenuity pathway analysis (IPA)

To understand the classical biological pathways involved in the biomarkers, we performed IPA functional enrichment analysis of them using IPA native software. Firstly, the differential expression data were imported, core analysis was selected and P < 0.05 was set to recalculate and filter out the meaningful results. Subsequently, pathway enrichment analysis was performed based on the IPA results, and upstream and downstream regulatory networks were constructed, aiming to clarify the key biological regulatory mechanisms.

Immune-infiltration analysis

To explore immune-infiltration levels between the YES_thrombus and NO_thrombus groups, the proportion of 29 immune cell subtypes was computed for each sample using the single sample gene set enrichment analysis (ssGSEA) algorithm. The differential immune cell distributions between the two groups were compared, and box plots were generated. A correlation analysis was also performed between the differential immune cells and feature genes.

Construction of a competing endogenous RNA (ceRNA) regulatory network

To investigate the potential regulatory mechanisms of the feature genes, the miRWalk and Starbase databases were utilized to predict miRNAs targeting feature genes. The miRNAs predicted by both databases were considered as co-miRNAs. Additionally, the relationships between lncRNAs and miRNAs were predicted using the Starbase and miRnet databases. A lncRNA-miRNA-mRNA network was then constructed to further explore the regulatory interactions.

Expression level validation and quantitative reverse transcriptase PCR (qRT-PCR)

The expression of feature genes was extracted from the YES_thrombus dataset and compared between the YES_thrombus and NO_thrombus groups. A correlation analysis was also conducted between the feature genes.

At the Second Hospital of Shanxi Medical University, whole blood samples were collected from ten thrombotic patients and ten normal individuals. For each sample, 500 µL of whole blood was collected, and 1 mL of TRIzol reagent (Ambion, USA) was added for RNA extraction. The RNA concentration and purity were evaluated, and the sample volume was adjusted accordingly. Reverse transcription was performed following the manufacturer’s instructions, and the reaction system was programmed to run at 25 °C for 5 min, 50 °C for 15 min, 85 °C for 5 s, and stored at 4 °C. Finally, qRT-PCR was carried out using a CFX96 real-time quantitative fluorescence PCR apparatus. The 2^−ΔΔCT method was applied to normalize the expression of feature genes to GAPDH. Intergroup differences were then analyzed. Table 1 contains the primer information.

Table 1 qRT-PCR primers for RGS1, HSF2, ARL4A, AAED1, WDR81, and GAPDH

Statistical analysis

All bioinformatics analyses were conducted in R language. The correlation analysis was performed using the Spearman method, with p < 0.05 indicating statistical significance.

Results

Quality control and preprocessing of raw sequencing data

Transcriptome sequencing data were obtained for each sample (Additional file 1). After filtering and quality assessment, the sequencing error rate at each base position was found to be less than 0.5%, making the data suitable for subsequent analyses (Additional file 2). Comparison of the data revealed an average unique alignment rate of 90% across the 36 samples, confirming the quality of the data for further analysis.

Screening of DEGs and key module genes

A box plot of gene expression levels across the 36 samples showed uniform expression, supporting their suitability for further study (Fig. 1a). A total of 42 DEGs were identified between the YES_thrombus and NO_thrombus groups, including 21 upregulated and 21 downregulated genes (Fig. 1b-c, Additional file 3). WGCNA was used to explore genes associated with thrombus formation, with clustering results indicating no outlier samples (Fig. 1d). The optimal soft threshold was determined to be 7, where R² equaled 0.85 and average connectivity approached zero (Fig. 1e). Five modules were identified using the dynamic tree-cutting algorithm, and after merging similar modules with a MEDissThres of 0.2, five modules remained (Fig. 1f-g). Among these, MEturquoise (Cor = 0.34, P = 0.04) exhibited the highest correlation with thrombus and was thus designated as the key module (Fig. 1h). The 3,230 genes within the key module were selected for further analysis.

Fig. 1

Quality control and preprocessing of raw sequencing data. (a) Box plot of gene expression for all samples. The horizontal axis indicates the name of the sample and the vertical axis indicates the total expression of the gene for each sample. (b) Volcano plot depicting differentially expressed genes (DEGs) between YES_thrombus and NO_thrombus samples. Dots within the dotted line in the upper left corner indicate down-regulated genes obtained by screening, and dots within the dotted line in the upper right corner indicate down-regulated genes obtained by screening. (c) Heatmap illustrating DEGs between YES_thrombus and NO_thrombus samples. The top graph represents a plot of the density distribution of differential gene expression in patients. The bottom graph represents the heat map of differential gene expression. Each small square indicates the expression of a different gene in each sample after normalisation, and its colour indicates how much it is present, with the higher the expression the redder the colour and the lower the bluer the colour. (d) Cluster dendrogram of 36 sequenced samples. (e) Identification of the optimal soft thresholds in WGCNA. (f) WGCNA module height parameters for merging. (g) Gene clustering dendrogram derived from module gene clustering and comparable expression patterns. (h) Heatmap of module correlations with thrombus

Screening and functional enrichment of candidate genes

A total of 7 candidate genes (RGS1, HSF2, ARL4A, KLHL15, AAED1, MMP1, and WDR81) were selected by intersecting DEGs and key module genes (Fig. 2a). Expression analysis revealed that all genes, except for WDR81, were downregulated in the YES_thrombus group (Fig. 2b-c). Enrichment analysis indicated that these candidate genes were associated with 33 GO-biological process (BP) terms, 7 GO-cell component (CC) terms, 7 GO-molecular function (MF) terms, and 7 KEGG pathways. GO term annotations suggested involvement in processes such as autophagosome formation and cellular response to misfolded proteins (Fig. 2d, Additional file 4). KEGG pathway analysis identified associations with the PPAR signaling and relaxin signaling pathways, among others (Fig. 2e, Additional file 5).

Fig. 2

Acquisition and functional enrichment of candidate genes. (a) Differentially expressed genes (DEGs) intersected with key module genes to yield candidate genes. (b-c) Differential expression of candidate genes between YES_thrombus and NO_thrombus groups, visualized in volcano (b) and box-and-line (c) plots. (d-e) Functional exploration of candidate genes: (d) Gene Ontology (GO) enrichment analysis; (e) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The closer the colour to blue the higher the significance, the lighter the colour the lower the significance. * p-value < 0.05, ** p-value < 0.01, *** p-value < 0.001, **** p-value < 0.0001

RGS1, HSF2, ARL4A, AAED1, and WDR81 were identified as featured genes

From the feature selection process, six LASSO-feature genes (RGS1, HSF2, ARL4A, AAED1, MMP1, and WDR81) were identified (Fig. 3a-b). Similarly, six SVM-feature genes (WDR81, ARL4A, AAED1, MMP1, HSF2, and RGS1) and six RF-feature genes (WDR81, AAED1, ARL4A, HSF2, KLHL15, and RGS1) were obtained (Fig. 3c-e). Five feature genes (RGS1, HSF2, ARL4A, AAED1, and WDR81) were selected by overlapping the LASSO, RF, and SVM feature genes (Fig. 3f).

Fig. 3

Machine learning screening of feature genes. (a-b) Six feature genes selected when λ was set to 0.0347 in the LASSO algorithm. (c-d) The six feature genes corresponding to the smallest model error filtered by support vector machine recursive feature elimination (SVM-RFE). (e) Random forest (RF) model used to screen six feature genes. (f) Overlap of feature genes identified by the three machine learning methods, shown in a Venn diagram

ROC curve analysis demonstrated that the feature genes exhibited strong predictive performance (Fig. 4a), with area under the curve (AUC) values exceeding 0.75 in the YES_thrombus dataset (Fig. 4a). Additionally, the diagnostic model achieved an excellent performance, with an (area under the curve) AUC of 0.917 (Fig. 4b). Confusion matrix analysis confirmed the model’s capacity to accurately identify thrombus complications post-fracture (Fig. 4c). Furthermore, the PPI network indicated that the feature genes were predominantly associated with the extrinsic component of the membrane and purine ribonucleoside binding functions (Fig. 4d).

The constructed ANN had better diagnostic performance

In the constructed ANN model, the predictive weight values of each biomarker ranged from − 9.72 to 3.13. Specifically, the predicted weight value of RGS1 was − 5.02, HSF2 was − 0.68, ARL4A was − 0.12, AAED1 was 0.51, and WDR81 was 1.01 (Fig. 4e). Subsequently, we plotted ROC curves to evaluate the performance of the model. The results showed that in the training set, the model had an AUC value of 0.9907 and its optimal cut-off value was 0.549.At this cut-off value, the model had a specificity of 0.944 and a sensitivity of 0.944 (Fig. 4f). The k-fold cross-validation results showed that an average accuracy of 86.7%, an average AUC of 0.90, an average sensitivity of about 0.83, and an average specificity of about 0.90, indicating that the model has a good and stable classification ability in most data division situations (Table 2, Additional file 6). These results indicated that the ANN model was able to accurately utilise the five biomarkers for diagnosis.

Fig. 4

Feature genes and diagnostic models for predictive power assessment. (a) Receiver operating characteristic (ROC) curves for feature genes. (b) ROC curves for diagnostic models. (c) Confusion matrix assessing the diagnostic model’s ability to differentiate patients with thrombotic complications after a fracture. (d) Co-expression networks between feature genes and genes with reciprocal relationships. Each dot represents a gene, the colour within the dot represents the possible function of the corresponding gene, the line colour represents the type of inter-gene interaction and the line thickness represents the intensity. (e) Constructing artificial neural networks. (f) ROC curves for artificial neural network models

Table 2 Results of k-fold cross-validation performance evaluation of ANN models

Correlation analysis of feature genes with clinical indicators

Correlation analysis revealed that RGS1 exhibited significant positive associations with five clinical indicators, with the strongest correlation observed with hematocrit (Fig. 5a-b, Additional file 7). In contrast, WDR81 showed a notable negative relationship with four clinical markers, while ARL4A demonstrated the strongest negative correlation with platelet volume distribution width (Fig. 5a-c, Additional file 7).

Fig. 5

Correlation analysis of feature genes with clinical indicators. (a) Bubble plots of correlations between characterised genes and clinical indicators. (b) Scatter plot showing the correlation between RGS1 and hematocrit. (c) Association between ARL4A and platelet volume distribution width, shown in a scatter plot

Single-gene GSEA of feature genes

To further explore the pathways and biological processes associated with the feature genes, single-gene GSEA analyses were performed. Functional enrichment results indicated that the feature genes were primarily enriched in GO terms related to cytosolic ribosomes, large ribosomal subunits, and other processes (Fig. 6, Additional files 8–12). KEGG pathway analysis revealed that HSF2, ARL4A, AAED1, and WDR81 were mainly enriched in pathways involving ribosomes and oxidative phosphorylation (Fig. 7, Additional files 13–16), while RGS1 was primarily linked to the NF-kappa B signaling pathway and TNF signaling pathway (Fig. 7, Additional file 17).

Fig. 6

Single-gene gene set enrichment analysis (GSEA) of feature genes, conducted using Gene Ontology (GO) as the background gene set

Fig. 7

Single-gene gene set enrichment analysis (GSEA) of feature genes conducted with Kyoto Encyclopedia of Genes and Genomes (KEGG) as the background gene set

IPA of feature genes

To elucidate the molecular mechanisms underlying the function of feature genes, a classical pathway analysis was conducted. IPA identified 28 significantly enriched pathways, including endocrine system disorders and organismal injury and abnormalities (Fig. 8a). Notably, RGS1 was regulated by a larger number of factors, whereas HSF2 activated multiple factors (Fig. 8b).

Fig. 8

Ingenuity Pathway Analysis (IPA) of feature genes. (a) Enrichment results of feature genes in the classical pathways of IPA. (b) Diagram illustrating the regulatory network of feature genes and their interactions with multiple factors

Immune-related analyses of feature genes

Subsequently, differences in the immune microenvironment between the YES_thrombus and NO_thrombus groups were examined. The bar chart displayed the proportion of 29 immune cell types in each sample (Fig. 9a). Significant variations were observed in major histocompatibility complex (MHC) class I and aDCs between the two groups (Fig. 9b). Further correlation analysis between feature genes and differential immune cells revealed that WDR81 had the strongest negative association with aDCs, whereas RGS1 showed the strongest positive correlation with aDCs (Fig. 9c-e, Additional file 18).

Fig. 9

Relationship between immune cells and feature genes. (a) Proportion of immune cell infiltration in the sequencing dataset. (b) Analysis of immune cell variances between YES_thrombus and NO_thrombus samples. (c-e) Heatmap (c) and scatter plots (d-e) illustrating the correlation between feature genes and differential immune cells. Each grid colour and size table shows the size of the absolute value of the spearson correlation coefficient, the redder the colour and the larger the shape of the circle the larger the absolute value of the spearson correlation coefficient. * p-value < 0.05, ** p-value < 0.01, *** p-value < 0.001, **** p-value < 0.0001

The CeRNA regulatory network of feature genes

Prediction analyses led to the identification of 6 co-miRNAs and 237 lncRNAs (Fig. 10a-b). The network consisted of 3 feature genes, 5 miRNAs (e.g., hsa-miR-93-5p, hsa-miR-320b), 236 lncRNAs (e.g., THUMPD3-AS1, SH3BP5-AS1, THRB-IT1), and 314 edges (Fig. 10c). Specific mRNA-miRNA pairs included HSF2-hsa-miR-139-5p, while miRNA-lncRNA pairs such as hsa-miR-139-5p-SH3BP5-AS1 were identified (Additional files 19–20).

Fig. 10

Construction of the competitive endogenous RNA (ceRNA) network. (a) Co-miRNAs predicted in the miRWalk and Starbase databases based on characterized genes. (b) Intersecting lncRNAs predicted in the Starbase and miRNet databases based on feature genes. (c) Diagram of the ceRNA network. Blue diamonds indicate lncRNAs, green pentagons indicate miRNAs, and red circles indicate mRMA

Expression analysis and PPI network of feature genes

Box plot analysis showed that, except for WDR81, all other feature genes were more highly expressed in the NO_thrombus group than in the YES_thrombus group (Fig. 11a). qRT-PCR validation confirmed significant differences in the expression of all five feature genes between the two groups, with AAED1, ARL4A, and WDR81 showing expression patterns consistent with the RNA sequencing results (AAED1 and ARL4A were downregulated in the thrombus group, while WDR81 exhibited the opposite trend) (Fig. 11b). Furthermore, RGS1 demonstrated strong positive correlation with HSF2 and a strong negative correlation with WDR81 (Fig. 11c), with WDR81 showing a negative correlation with all other feature genes, consistent with the expression analysis findings (Fig. 11c).

Fig. 11

Expression validation and correlation analysis of characterized genes. (a) Expression visualization of feature genes in the YES_thrombus dataset. (b) Visualization of feature gene expression in clinical samples. (c) Heatmap of correlations between feature genes. Each grid colour and size table shows the size of the absolute value of the spearson correlation coefficient, the redder the colour and the larger the shape of the circle the larger the absolute value of the spearson correlation coefficient. * p-value < 0.05, ** p-value < 0.01, *** p-value < 0.001, **** p-value < 0.0001

Discussion

HF represents a major public health concern, marked by high morbidity, a disproportionate risk of disability and mortality, and significant healthcare expenditures as the population ages [19]. DVT, a common complication in HF patients, is primarily due to hypercoagulability and reduced mobility, with an increased prevalence observed in patients undergoing delayed surgery [20]. The mortality rate associated with DVT, particularly as VTE develops, ranges from 10 to 30% [21]. Inflammation has been shown to exacerbate venous thrombosis by inducing vascular endothelial cell damage, leading to the release of adhesion molecules [22, 23]. Thrombosis further activates endothelial cells, triggering the release of cytokines such as interleukin-6, tissue factor, and tumor necrosis factor (TNF)-α. The inflammatory response is amplified by the formation of blood clots, creating a vicious cycle between the two processes. Prevention of DVT is crucial, and understanding the mechanisms behind venous thrombosis, as well as identifying proteins that regulate inflammatory factors, is of significant value for delaying or reducing its formation and improving prognostic predictions. To identify potential biomarkers, an integrated bioinformatics approach was employed, utilizing self-sequencing and three machine learning algorithms (Lasso, RF, and SVM). Five feature genes—RGS1, HSF2, ARL4A, AAED1, and WDR81—were identified through this analysis. The model demonstrated a high diagnostic performance, with an AUC value of 0.85 in the training set and individual AUC values greater than 0.75 for each of the five genes, indicating their strong potential as biomarkers. Experimental validation was conducted to further substantiate these findings.

RGS1, a GTPase activator, serves as a negative regulator of chemokine receptor signaling in B and T cells. It was found to be enriched in the NF-κB signaling pathway, suggesting that RGS1 may contribute to DVT through this pathway. NF - κ B signaling plays a crucial role in the vascular system and various cell types involved in thrombotic inflammatory responses. It disrupts the balance between normal coagulation and fibrinolysis by regulating the interactions between endothelial cells, platelets, and inflammatory cells, thereby triggering the formation of blood clots [24, 25]. The role of HSF2 in neurological protection, inflammation, and other aspects has gradually been revealed [26]. Studies have found that HSF2 can reduce the secretion of IL-1 β and TNF - α, and promote TGF - β by enhancing Smad2/3 phosphorylation, thereby exerting anti-inflammatory effects [26]. The IPA enrichment analysis results of this study showed that the HSF2 gene is involved in regulating TNF and other related pathways, which is consistent with previous research findings. WDR81 is associated with vesicle transport and autophagy, and can inhibit the activity of PI3K, increasing exosome secretion by preventing the formation of class III PI3K complexes [27]. Extracellular vesicles have the function of clotting agents. Studies have found that the formation of extracellular vesicles leads to the externalization of phosphatidylserine (PS). In the presence of calcium ions, externalized PS can promote the activation of prothrombin, thereby promoting the formation of thrombin [28, 29, 30]. This study found that, HSF2, ARL4A, AAED1,and WDR81 are associated with oxidative phosphorylation pathways, which are closely related to cellular energy metabolism [31]. Abnormal oxidative phosphorylation system can cause mitochondrial dysfunction [32]. Bing Bo Yu’s study found that mitochondrial dysfunction in endothelial progenitor cells weakens their migration ability, leading to a decrease in angiogenesis [33]. However, their specific roles in deep vein thrombosis still need further experimental exploration to clarify.

DVT is primarily a disease of the elderly, with rare occurrences before late adolescence [34]. Advanced age is considered an independent risk factor for preoperative DVT, often associated with a prothrombotic state and impaired vascular function [35, 36, 37, 38]. In this study, DVT was predominantly found in elderly patients with low expression levels of ARL4A.

ARL4A, a member of the ADP-ribosylation factor (ARF) family, includes six ARF proteins (ARF1-6) and more than 20 ARL proteins. These proteins are regulated by developmental processes and show decreased expression levels as the body matures. In adults, ARL4A expression is restricted to specific tissues but has been implicated in various cellular events [39]. Gene expression profiling in zebrafish has demonstrated that ARL4A is expressed in hematopoietic tissues during development [40]. Furthermore, research by Lin et al. confirms that phosphorylation of ARL4A promotes its binding with the HYPK chaperone. This interaction stabilizes ARL4A’s targeting to the plasma membrane, thus promoting cell migration. Interestingly, targeting HYPK in the context of ARL4A overexpression enhances cell migration, while knocking down HYPK inhibits this process [39]. Guo et al. demonstrated that the expression of ARL4aa, a member of the ARL4A subfamily, in hematopoietic tissues during embryonic development helps maintain Golgi complex integrity, which in turn facilitates the maturation of the Notch receptor in hematopoietic stem cells [41]. Enhancing Notch receptor expression in venous endothelial cells has been shown to facilitate angiogenesis and the maturation of the vascular system [42].The results of this study showed that ARL4A is enriched in the oxidative phosphorylation pathway, and the ARL4A gene is highly expressed in the non thrombotic group. We speculate that ARL4A may activate the oxidative phosphorylation pathway, promote oxidative phosphorylation transport, enhance glucose metabolism to promote energy generation, and thereby enhance the migration ability of EPCs, thereby playing a role in repairing endothelial damage. In addition, we also observed a negative correlation between ARL4A expression and platelet distribution width (PDW), which is a marker of platelet activation and can quantify changes in platelet size [43]. Sonay O ğ uz et al. closely monitored the PDW values of patients with deep vein thrombosis and found that elevated PDW levels were associated with poor prognosis [21]. In summary, ARL4A may play an important role in the formation and development of deep vein thrombosis. In the future, it is necessary to further investigate the function and mechanism of action of ARL4A to explore its potential applications in clinical diagnosis and treatment.

VTE has a high recurrence rate, with approximately 30% of patients experiencing relapse within the next 10 years [44]. Inflammation plays a key role in converting the coagulation system into a procoagulant state, thereby increasing the risk of thrombosis [45]. Traumatic fractures, which are closely linked to inflammation, further elevate the risk of VTE [46]. In this study, an IPA was performed on ARL4A characteristic genes, revealing their involvement in pathways such as inflammatory diseases, inflammatory responses, body damage and abnormalities, cell movement, blood system development and function, and humoral immunity. According to Virchow’s classical theory, venous wall injury, slow blood flow, and a hypercoagulable state are the three primary factors contributing to DVT [47]. These factors are closely related to inflammation, where inflammatory substances trigger a hypercoagulant state and cause endothelial damage, which further exacerbates clot formation [48]. While slow venous blood flow is considered the most significant factor, blood stasis alone is insufficient to initiate clot formation [49]. However, the combination of slow venous blood flow with vascular injury or a hypercoagulable state significantly increases the risk of clot formation [50]. For DVT patients, anticoagulant therapy addresses the hypercoagulable state of the blood but does not alter slow blood flow or vascular damage. This limitation contributes to the recurrence of DVT and highlights the inadequacy of anticoagulants alone in preventing venous thrombosis.

Immune deficiency factors have also been identified as contributors to DVT [51]. Several studies have shown that immunological disorders, compared to healthy populations, are significant risk factors for VTE [52]. Our research shows that the ssGSEA algorithm identified differences in the functional scores of aDCs and MHC class I between the fracture complicated thrombus group and the fracture group. Mature dendritic cells can produce IL-23 and activate Th17 cells, thereby promoting the production of IL-17 [53]. Studies have shown that an increase in IL-17 in serum is positively correlated with the formation of DVT [54]. This study found a positive correlation between RGS1 and HSF2 and aDCs. In addition, qRT PCR results showed that the expression of RGS1 and HSF2 genes was upregulated in patients with deep vein thrombosis. Therefore, we speculate that RGS1 and HSF2 may play an important role as key regulatory factors in dendritic cell-mediated immune response and participate in the pathological mechanism of deep vein thrombosis by regulating the production of IL-17. ARL4A and RGS1 are positively correlated with MHC class I, indicating that both fracture related thrombosis and fractures may affect disease progression through MHC class I. Our research indicates that the expression level of MHC class I molecules is lower in the fracture thrombus group, which may lead to a decrease in CD8 + T cells and immune disorders, and even further exacerbate the pathogenesis of deep vein thrombosis. In summary, ARL4A and other genes provide new strategic directions for the prevention and treatment of DVT, which have important clinical significance and potential application value.

In this study, biomarkers such as ARL4A were associated with fractures and thrombosis. In clinical applications, postoperative detection of these biomarkers can identify risks early and reduce the risk of DVT through early use of anticoagulant therapy. At the same time, after DVT occurs, by detecting the expression of these markers, timely follow-up of the disease situation, optimization of treatment plans, and reduction of drug abuse and costs can be achieved. In addition, these findings lay the foundation for new strategies and target research for DVT treatment. However, given the complexity of the pathology of deep vein thrombosis in the lower limbs, it is crucial to explore the specific roles of these biomarkers in DVT, and continuous in vitro and in vivo experimental research is needed. In this regard, we plan to further explore the biological characteristics of biomarkers through more in-depth functional validation in the future, such as constructing overexpression cell models to investigate their biological effects in cases of abnormal expression. In addition, gene editing technology will be used to overexpress or knock out genes in specific animals to more directly observe their functions and effects in the living environment. In addition, this study lacks external validation to confirm the reliability and stability of its conclusions. We plan to explore possible external validation opportunities through multi institutional collaboration and other means, with the aim of strengthening our conclusions through broader external validation. At the same time, by comparing the actual number of immune cells in the whole blood cell count, the accuracy of predicting the content of immune cells can be verified and evaluated, enhancing the credibility of the prediction results.Finally, in this study, the extensive use of R language software packages for data analysis may face difficulties such as version compatibility, parameter optimization, and computational resource limitations. In the future, we will need to continuously optimize parameters, adjust models, and improve research efficiency and accuracy.

Conclusions

This study identifies RGS1, HSF2, ARL4A, AAED1, and WDR81 as central genes associated with DVT. Notably, ARL4A is linked to PDW and may play a role in autophagosome signaling pathways related to DVT onset. Therefore, ARL4A may serve as a key biomarker in the development of HF-induced DVT.