Differential gene expression profiling and machine learning-based discovery of key genetic markers in VTE and CKD

Gene expression data associated with CKD and VTE were obtained from the Gene Expression Omnibus (GEO) database. For CKD, the GSE66494 dataset, based on the GPL6480 platform, included placenta samples from 48 CKD patients and 5 healthy controls, and the GSE37171 dataset, based on GPL570, contained 75 CKD samples and 40 healthy controls, which was used to validate diagnostic efficiency and the robustness of simulated gene expression results. For VTE, the GSE19151 dataset, based on the GPL571 platform, comprised 70 VTE patient and 63 healthy control blood samples, and the GSE48000 dataset, based on GPL10558, including 15 VTE and 11 healthy control samples, served as an independent validation set. Additionally, single-cell RNA sequencing data for CKD (GSE198621) incorporating samples from 3 CKD cases and 3 normal controls were utilized to validate the expression levels of key genes at the single-cell level. All datasets were retrieved in accordance with the GEO repository protocols. Prior to analysis, expression data were processed for quality control and normalized, and batch effect correction was performed using the ComBat method to reduce heterogeneity arising from different platforms and experimental conditions.

Differentially expressed genes (DEGs) were identified using the "limma" R package for both CKD and VTE datasets. For the GSE19151 dataset, genes with |log2 fold change (FC)| ≥ 0.5 and p < 0.05 were considered differentially expressed, while for the GSE66494 dataset, the threshold was set at |log2FC| ≥ 1 and p < 0.05. The different thresholds were selected considering the sample sizes and variability in each dataset, consistent with previously published studies. To minimize potential bias due to platform differences, all datasets were subjected to batch effect correction prior to DEG analysis. Functional enrichment analysis was conducted using the "clusterProfiler" package, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment. Pathways with p < 0.05 were considered significantly enriched. For validation at the single-cell level, key DEGs were mapped to the CKD single-cell dataset to confirm their expression patterns across different cell types.

The crosstalk genes for the two diseases were obtained by hybridizing the DEGs of CKD with those of VTE. To reduce the high false discovery rate, these genes were further screened using the Wilcoxon test (p < 0.05) between CKD and control samples, and a similar analysis was carried out between VTE and control samples.

To identify key gene features for the diagnostic model, we utilized three machine learning algorithms: LASSO (Least Absolute Shrinkage and Selection Operator), Random Forest, and SVM-RFE (Support Vector Machine Recursive Feature Elimination). For LASSO, we employed Cox proportional hazards regression using the "glmnet" R package. Cross-validation was applied to optimize the regularization parameter (λ) via 10-fold cross-validation, where the λ value that minimized the mean cross-validation error was selected. The coxnet option was used to perform survival analysis, with 100 iterations to ensure model stability.

For Random Forest, we used the randomForest R package to rank the importance of CKD-related marker genes. The model was configured with 500 trees, which is the default setting, and the number of variables at each split was set to the square root of the total number of features. To minimize overfitting, a 5-fold cross-validation was performed, and genes with an importance score greater than 0.8, as determined by the mean decrease in Gini index, were considered for further analysis.

SVM-RFE was implemented using the e1071 R package with a linear kernel. The cost parameter (C) was set to 1 to balance the trade-off between maximizing the margin and minimizing classification errors. The algorithm was run for 10 iterations to identify the most relevant features, with genes ranked based on the average performance across all iterations. A 10-fold cross-validation approach was used to assess the optimal number of features to retain, focusing on accuracy and error rates.

In the final step, a hybrid approach was adopted, combining the results from LASSO, Random Forest, and SVM-RFE. This consensus approach helped to identify a set of robust and reliable gene features for further validation. The diagnostic performance of these marker genes was evaluated using Receiver Operating Characteristic (ROC) curve analysis, and the Area Under the Curve (AUC) was calculated to assess their ability to distinguish CKD samples. To validate the diagnostic potential, the expression profiles of the selected genes were analyzed in external datasets, including GSE37171, GSE19151, and GSE48000.

The diagnostic nomogram was created using the "rms" R package. Patient scores were calculated based on the expression of individual core genes, with the total risk score defined as the sum of risk scores for all individual genes. The dCA curve, calibration curve, and ROC curve were utilized to evaluate the diagnostic value of the nomogram for CKD.

GSEA was used to analyze the biological function of the key genes. Based on the gene sets of 28 immune-related cells, the ssGSEA algorithm from the R package "GSVA" was employed to assess the immune activity for each sample. Differences in immune infiltration between the CKD and VTE groups were analyzed. Additionally, the correlation between immune infiltration levels and the expression of key genes in the two diseases was investigated.

For single-cell characterization studies, the scRNA-seq dataset GSE198621 was analyzed using the standard protocol of "Seurat". Cells with fewer than 100 genes, more than 5000 total genes, and mitochondrial gene content exceeding 20% were further filtered out. The R package "Harmony" was employed to mitigate batch effects between samples. The "FindVariableFeatures" function was used to identify the top 2000 variably expressed genes. Cells were annotated using "SingleR", and the expression of key genes in different cell types was subsequently validated.

Peripheral blood samples were collected from 15 patients with clinically diagnosed CKD (stage 3–4) and 15 age- and sex-matched healthy controls at Huizhou First People's Hospital, following approval by the institutional Ethics Committee and acquisition of written informed consent from all participants.

Peripheral blood mononuclear cells (PBMCs) were isolated from whole blood using Ficoll-Paque density gradient centrifugation according to standard protocols. Total RNA was extracted from PBMCs using the FastPure Cell/Tissue Total RNA Isolation Kit (Vazyme, China) following the manufacturer's instructions. RNA concentration and purity were assessed by NanoDrop spectrophotometry. Reverse transcription was performed using the ReverTra Ace qPCR RT Master Mix with gDNA Remover (Toyobo, Japan) to synthesize complementary DNA (cDNA).

qRT-PCR was conducted using SYBR Premix Ex Taq II (Takara, Japan) on a real-time PCR system. The thermocycling conditions were as follows: initial denaturation at 95 °C for 10 minutes, followed by 45 cycles of 95 °C for 5 seconds and 60 °C for 30 seconds. The housekeeping gene GAPDH was used as the internal control. Primer sequences are listed in. 1

All statistical analyses were performed using R (version 4.3.2) and GraphPad Prism 9.0. Differential expression was determined with the "limma" package, and p values were adjusted using the Benjamini–Hochberg method. Comparisons between two groups were analyzed using Student's t-test or the Wilcoxon rank-sum test as appropriate. Correlations were assessed by Spearman's method. Diagnostic performance was evaluated by ROC curve analysis, and qRT-PCR results were calculated using the 2^–ΔΔCt method. A two-tailed p < 0.05 was considered statistically significant. In figures, significance is indicated as follows: p < 0.05 (*), p < 0.01 (**), p < 0.001 (***), and p < 0.0001 (****).

Differential analysis revealed that between VTE patients and the control group, there were 637 differentially expressed genes, of which 224 were downregulated and 413 were upregulated (Figures 1A, C). Analysis of differences between CKD patients and the control group showed 671 differentially expressed genes, including 572 downregulated genes and 99 upregulated genes (Figures 1B, D).

The biological functions of differentially expressed genes in different diseases were explored through GO and KEGG enrichment analyses. Differentially expressed genes in VTE patients were mainly enriched in cytoplasmic translation, cell activation involved in the immune response, and ribosomal subunits (Figure 2A). KEGG enrichment analysis showed that VTE differentially expressed genes were mainly enriched in oxidative phosphorylation, reactive oxygen species associated with chemical carcinogenesis, and signaling pathways of the 2019 coronavirus disease (Figure 2C). Differentially expressed genes in CKD patients were primarily enriched in regulation of muscle contraction, cation-transporting ATPase complex, and endopeptidase inhibitor activity (Figure 2B). KEGG enrichment analysis showed that CKD differentially expressed genes were primarily enriched in signaling pathways such as pancreatic secretion, protein digestion and absorption, and vascular smooth muscle contraction (Figure 2D).

Venn diagram results identified 23 overlapping genes differentially expressed in both CKD and VTE diseases, including CCNL2, HNRNPA0, PI4KA, FOS, HBD, TSC22D3, DUSP1, and ZNF692 (Figure 3A). Results from the Wilcoxon test indicated that these 23 genes were significantly differentially expressed in both CKD and VTE diseases (Figures 3B, C).

Upon performing 10-fold cross-validation, the optimal lambda for the LASSO algorithm was determined to be 0.002. The penalty parameter (l) corresponding to the minimal partial likelihood deviance was selected, yielding 5 feature genes: CCNL2, HNRNPA0, PI4KA, FOS, and HBD (Figures 4A). In contrast, for SVM-RFE, the model classifier with the highest accuracy and the lowest error had 2 features: HNRNPA0 and PI4KA (Figures 4B, C). The random forest algorithm selected 7 genes with significance greater than 0.8, including HNRNPA0, FOS, PI4KA, DLST, NGRN, CCNL2, and SEC16A (Figures 4D, E). Based on the aforementioned three machine learning algorithms, two key genes were selected: HNRNPA0 and PI4KA (Figure 4F). The AUC for the key genes HNRNPA0 and PI4KA were both 1.000 (Figures 4G, H).

The "rms" R package was used to construct a nomogram model for RM diagnosis based on the feature genes (HNRNPA0 and PI4KA) (Figure 5A). Calibration curves indicated a minimal discrepancy between the predicted CKD risk and the actual CKD risk, suggesting that the nomogram model has high accuracy (Figure 5B). Decision curve analysis (DCA) demonstrated that CKD patients could benefit from the nomogram model (Figure 5C). The validity of this model was also confirmed by gene and nomogram ROC curve analyses (Figures 5D, E).

GSEA results suggested that HNRNPA0 was primarily involved in pathways like Amoebiasis, Ascorbate and aldarate metabolism (Figure 6A). PI4KA was mainly enriched in pathways like Ascorbate and aldarate metabolism and Butanoate metabolism (Figure 6A). Results from ssGSEA for immune cell infiltration revealed that VTE patients exhibited higher infiltration levels of activated CD8 T cells, effector memory CD4 T cells, CD56bright natural killer cells, and eosinophils, while having lower infiltration of activated B cells, central memory CD4 T cells, and plasmacytoid dendritic cell-like cells (Figure 6B). For GSE19151, the immune cell infiltration analysis revealed that HNRNPA0 was significantly associated with 15 immune cell types, positively correlating with cell types such as Central memory CD4 T cell, Plasmacytoid dendritic cell, Effector memory CD8 T cell, and negatively correlating with 4 immune cell types including Type 17 T helper cell, Macrophage (Figure 6D). PI4KA showed positive correlations with 9 immune cell types and negative correlations with 8 cell types including Type 17 T helper cell, Effector memory CD4 T cell (Figure 6E). For GSE66494, ssGSEA analysis indicated a higher infiltration of Type 2 T helper cells and Immature B cells in CKD patients (Figure 6C). Within GSE66494, the correlation analysis between immune cells and key genes showed that HNRNPA0 negatively correlated with Central memory CD8 T cell, Neutrophil (Figure 6F). PI4KA exhibited a positive correlation with cell types including Type 17 T helper cell and negative correlations with cell types such as Effector memory CD4 T cell, Type 2 T helper cell (Figure 6G).

For the GSE37171 dataset, the expression levels of HNRNPA0 and PI4KA in CKD patients were lower than in normal samples (Figures 7A, B). The AUC for HNRNPA0 and PI4KA were 0.720 and 0.859, respectively (Figures 7E, F). Similarly, differential expression analysis of the key genes in the GSE48000 dataset produced consistent results (Figures 7C, D). Furthermore, the AUC for HNRNPA0 and PI4KA were 0.985 and 1.000, respectively (Figures 7G, H). ROC curve analysis for key genes in GSE19151 showed that the AUC for HNRNPA0 and PI4KA were 0.839 and 0.881, respectively (Figures 7I, J).

Following the data filtering and integration described in the methods, we obtained gene expression profiles for 25,449 cells from control samples and 29,463 cells from CKD samples (Figures 8A, B). Clustering all cells, we identified 11 cell clusters (Figure 8C). Annotation and visualization of the 11 cell clusters resulted in eight cell types, such as Epithelial cells and Macrophages (Figures 8D, E). Finally, the expression of key genes across different cell types was displayed (Figures 8F, G).

To validate the expression of key genes identified from the transcriptomic and single-cell analyses, we performed qRT-PCR on PBMCs collected from 15 CKD patients and 15 healthy controls. Consistent with the multi-omics and single-cell findings, both HNRNPA0 and PI4KA exhibited significantly lower mRNA expression levels in CKD patients compared to healthy individuals (Figures 9A, B).