GEO (http://www.ncbi.nlm.nih.gov/geo), a public genomics data repository, is created and maintained by the National Center for Biotechnology Information (NCBI) and includes high-throughput gene expression data, chips, and microarrays. The gene expression profiles of MDS, namely, GSE4619, GSE19429, and GSE58831, were sourced and downloaded from GEO. The aforementioned validation set dataset comprises 397 MDS patients and 35 normal healthy controls. GSE2779, utilized as a validation set, includes 159 MDS patients and 17 normal healthy controls.

Initially, the GEOquery package was used to convert the probe matrix to the gene matrix, incorporating the probe annotation file in the process. In instances where multiple probes corresponded to the same gene, the expression value of the gene was determined by calculating the average value of the probes. To mitigate the batch effects arising from different platforms, the sva package was employed, given that the datasets originated from various sources.

Analysis of DEGs between MDS patients and normal controls was conducted using the limma package. The selection criteria for DEGs included a p-value <0.05 and an absolute fold change (FC) > 1.

The enrichment of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was analyzed using the R package ‘clusterProfiler’. Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA) were utilized to investigate differences in biological functions among various expression groups of hub genes, thereby revealing their underlying mechanisms.

LASSO regression, a machine learning algorithm, is commonly used for fitting generalized linear models. This algorithm is distinguished by its ability to perform variable selection and complexity regularization simultaneously (Cheung-Lee and Link, 2019). The parameter λ adjusts the complexity level, whereby higher values impose a greater penalty on linear models with many variables. This approach yields a smaller number of selected genes, producing a more concise and representative set of key genes. The glmnet package in R was utilized in this study for LASSO regression analysis of candidate hub genes. The optimal value of λ was determined via 10-fold cross-validation, by selecting the value that produced the smallest criterion.

The random forest model, a machine learning technique, employs multiple independent decision trees for predicting classification or regression (Rigatti, 2017). In this study, the R package ‘randomforest’ was utilized to construct the random forest model. The optimal number of variables was determined by calculating the average error rate of candidate hub genes. Subsequently, the error rate was evaluated for tree numbers ranging from 1 to 500, selecting the number with the lowest error rate. Following the determination of the parameters, the random forest tree model was constructed. The feature importance score for each candidate hub gene was then determined, selecting genes with an importance value greater than 1.

SVM-RFE is a widely utilized supervised machine learning protocol for classification and regression (Frost and Amos, 2017). The “e1071” package in R was used, applying the support vector machine recursive feature elimination (SVM-RFE) algorithm based on nonlinear support vector machines to identify hub genes.

To assess the accuracy of the hub genes identified by machine learning, ROC curves were generated for MDS patients and normal controls in the training group. A larger area under the curve (AUC) indicates higher accuracy of the gene as a hub gene in MDS. This method’s effectiveness was further validated in the validation group.

The CIBERSORT algorithm was employed to calculate the differential abundance of 22 types of immune infiltrating cells. Heatmaps and violin plots displaying the correlation of immune cells were prepared using the R packages “corrplot” and “ggplot2”. This study utilized the CIBERSORT algorithm to determine the relative proportion of different immune cells in MDS and normal control groups.

Analysis of the correlation between hub genes and immune cells was conducted using the Spearman correlation coefficient.

The human MDS cell lines, SKM-1, were acquired from Yuchi Biotechnology Ltd. (Shanghai, China) and cultured in Dulbecco’s modified Eagle medium supplemented with 10% fetal bovine serum. These cells were maintained at 37 °C in a cell culture incubator with 5% CO₂.

Total RNA was extracted using TRIzol reagent and processed for cDNA synthesize, adhering strictly to the manufacturer’s protocol (Cat: 11119ES60, YEASEN Biotechnology Ltd., Shanghai, China). The resulting cDNA was then subjected to RT-PCR using QuantStudio 6 Pro (Thermo Fisher Scientific, Waltham, MA, United States) in conjunction with SYBR Green PCR Master Mix. Fold changes in mean values were determined using the double delta CT method. Each experiment yielded three independent datasets. The primer sequences were as follows: LDLRAD4: 5′- GTT​GCA​CTT​AGG​CTG​GGT​CT -3′ (F); 5′- AGG​TGA​GGG​GCA​GAG​AGA​AA-3′ (R). GADPH:5′-GGAGCGAGATCCCTCCAAAAT-3′(F); 5′-GGC​TGT​TGT​CAT​ACT​TCT​CAT​GG-3′ (R).

Cells were harvested and incubated with RIPA lysis buffer on ice for 30 min to extract total protein. Following centrifugation at 12,000 rpm for 15 min at a temperature of 4 °C, the protein-rich supernatant was gathered. The samples were then applied to a gel for sodium dodecyl-sulfate polyacrylamide gel electrophoresis and subsequently transferred onto a 0.22-μm polyvinylidene difluoride membrane. To block the membrane, 5% skimmed milk was used at a temperature range of 24 °C–30 °C for 30 min, followed by incubation with primary antibody at 4 °C. After three washes with PBS-Tween 20 (PBST), the membrane was incubated with a secondary antibody (Cat: SSA004, 1:3000, Sino Biological Inc., China) for 1 h at 24 °C–30 °C. After additional washes with PBST, the protein bands were detected using a Tanon-5200 image analyzer (Thermo Fisher Scientific, Waltham, MA, United States). ImageJ software was employed to analyze the film strips, normalizing the protein intensities to the corresponding β-actin bands. The antibodies used and their dilution ratios were as follows: anti-LDLRAD4 (Cat:PA5-70568, 1.0 μg/mL, Thermo Fisher Scientific, Waltham, MA, United States), anti-β-Tubulin (Cat: E021040-01, 1:5000, EarthOx, United States).

Transfected cells were distributed into 96-well plates at a density of 1,000 cells per well, containing 100 μL of medium in each. Following incubation periods of 24, 48, and 72 h, 10 μL of CCK-8 solution was introduced to each well, and the plates were incubated in the dark at 37 °C for 2 h. Subsequently, the absorbance at 450 nm was measured to assess cell viability.

Trypsin-treated cells were harvested, rinsed with PBS, suspended in pre-cooled 70% ethanol, and fixed at 4 °C for 2 h. Subsequently, the supernatant was aspirated, and the cells were washed again with PBS. Propidium iodide staining solution was then added, and the cells were incubated in the dark at 37 °C for 30 min.

Cells were digested and resuspended in 100 μL 1× binding buffer. Stained with 5 μL Annexin V and 10 μL PI (Cat:40302ES20,YEASEN, Shanghai, China) for 15 min at room temperature (37 °C) in the dark, and finally 400 μL 1× binding buffer was added. The number of apoptotic cells was analyzed by using FACS Aria flow cytometer with CellQuest software and the data were analyzed with FlowJo software.

SKM-1 were plated into 6-well plates and allowed to proliferate to a 50%–60% confluence prior to transfection. Cells were divided into negative control (pcDNA3.1-NC) and transfection (pcDNA3.1-LDLRAD4) groups. Transfection was performed following the instructions for Lipo8000 protocol (Beyotime Biotechnology Ltd., Shanghai, China). Eight hours later, the serum-free Opti-MEM was substituted with complete medium containing 10% serum. Cells were further incubated for 48 h at 37 °C with 5% CO2) after transfection.

Statistical analysis of the data in this study was conducted using R (version 4.3.2). For continuous variables assumed to follow a normal distribution between the two groups, t-tests were applied. The exploration of the correlation between gene expression and immune cell components was conducted using the Spearman rank correlation test. The experimental data were plotted and analyzed using GraphPad Prism 8.0, employing the Student’s t-test for statistical significance. The threshold for statistical significance was established at a p-value <0.05.

The datasets GSE4619, GSE19429, and GSE58831 were selected from the GEO database to constitute the training group, including 397 MDS patients and 35 normal controls. Principal component analysis (PCA) was performed on these datasets to facilitate subsequent differential gene analysis, ensuring accuracy and comparability post de-batch processing, as demonstrated in the figures below (Figures 1A,B).

Analysis of the merged and normalized datasets revealed 108 significantly differentiated genes in MDS. Volcano plots visualized the expression changes of these genes (Figure 1C), featuring 14 upregulated genes such as IFITM1, KCNK5, TRBV27, PABPC4L, MAMDC2, and 94 significantly downregulated genes including GPR176, OR7A5, ARPP21, LDLRAD4, IGHV5-78 (Figure 1D).

Subsequent KEGG analysis compared the MDS group with the normal control group, highlighting three strongly relevant signaling pathways: the primary immunodeficiency-related pathway, the hematopoietic cell regulatory signaling pathway, and the B-cell receptor signaling pathway (Figure 1E). Additionally, GO analysis identified biological processes closely associated with immune responses in terms of BP, including lymphocyte differentiation, monocyte differentiation, and regulation of cell-cell adhesion (Figure 1F). Notably, MDS represents a group of clonal diseases originating from the malignant transformation of hematopoietic stem cells (HSCs). A distinguishing feature of MDS is its ineffective hematopoiesis (Bazinet and Bravo, 2022). Furthermore, research into the bone marrow microenvironment and immune abnormalities forms a significant area in the study of MDS pathogenesis (Lee et al., 2021). The results of this study align closely with the current state of research, reinforcing the validity and credibility of our findings.

Three machine learning methods were employed to identify hub genes closely related to MDS: LASSO regression analysis, SVM-RFE, and the random forest model. These methods were utilized to screen key candidates among differentially expressed genes. Initially, LASSO regression analysis was used to screen for differentially expressed genes, and the stability and reliability of the results were ensured through a 10-fold cross-validation method. Ultimately, 23 key genes were identified (Figures 2A,B). Subsequently, the SVM-RFE method was employed to identify 23 key genes, with an error rate of 0.0136 and an accuracy rate of 0.986 (Figures 2C,D). Ultimately, a random forest model was used to rank the importance of all genes, focusing on those with importance scores greater than 1 (Figures 2E,F). To further enhance the reliability of the results, an intersection analysis was performed on the key genes identified by each method, ultimately pinpointing three pivotal genes closely associated with MDS: LDLRAD4, KCNK5, and FAM43A (Figure 2G).

An in-depth analysis was conducted to validate the reliability of the three hub genes, focusing on their expression levels in the training group. It was found that the expression levels of LDLRAD4 and FAM43A were significantly downregulated in the MDS group compared to the normal control group, while KCNK5 showed a significant uptrend (Figure 2H). Subsequently, the area under the receiver operating characteristic curve (AUC-ROC) was calculated for these three genes. The AUC values were: 0.943 for LDLRAD4 (95% CI: 0.915-0.966), 0.854 for KCNK5 (95% CI: 0.812-0.892), and 0.845 for FAM43A (95% CI: 0.786-0.901) (Figure 2I). Given that the AUC values of these genes were all above 0.8, this indicates their high diagnostic efficiency in predicting MDS.

Based on the aforementioned analyses, LDLRAD4 was selected as the core hub gene for this study. To verify the stability and reliability of LDLRAD4, the GSE2779 dataset from the GEO database was employed as the validation set. The validation set included 17 normal controls and 159 MDS cases. Re-analysis of LDLRAD4’s expression levels and AUC-ROC values in the validation set revealed high consistency with the training set results. Specifically, LDLRAD4’s expression levels were significantly downregulated in MDS patients, with an AUC value of 0.984 (95% CI: 0.929-1.000) (Figure 2J). This result not only confirms LDLRAD4’s reliability as a hub gene but also suggests its potential applications in the diagnosis and prediction of MDS.

Having identified LDLRAD4 as the focal point of this study, its function and mechanism were explored in depth. LDLRAD4’s mRNA expression level in MDS was categorized based on the median value, dividing it into low and high expression groups. Subsequent differential analysis of these groups led to the identification of 26 differential genes, including MAMDC2, RPS4Y1, MME1, and EBF1 (Figure 3A). Co-expression analysis was conducted on the obtained differential bases, sorting them by correlation strength. Results indicated a positive correlation of LDLRAD4 with genes including MME, EBF1, VPREB1, KIAA0226L, and PAX5, and a negative correlation with MAMDC2 and RPS4 (Figure 3B).

To further elucidate LDLRAD4’s role in MDS, GSVA analysis of GO and KEGG pathways was performed for its high and low expression groups. In terms of GO, the high-expression group showed enrichment in: GOBP: Positive Regulation of Response to Extracellular Stimulus, GOBP: Negative Regulation of Histone H3 K9 Methylation, and GOMF: Cardiolipin Binding, among others. The low-expression group showed enrichment in: GOBP: Neuromuscular Junction Development, GOBP: Bone Trabecula Morphogenesis, and GOBP: Regulation of Membrane Repolarization During Cardiac Muscle Cell Action Potential, among others (Figure 3C). KEGG enrichment analysis revealed for the high-expression group: Vasopressin Regulated Water Reabsorption, Focal Adhesion, Arrhythmogenic Right Ventricular Cardiomyopathy (ARVC), GnRH Signaling Pathway, Other Glycan Degradation, and others; for the low-expression group: Primary Immunodeficiency, Cell Cycle, Non-Homologous End Joining, and others (Figure 3D).

Given the GO analysis results, an in-depth exploration of the regulation of Smad protein complex assembly, a critical biological process, is warranted. Additionally, enrichment in the primary immunodeficiency-related pathway and B-cell receptor signaling pathway in KEGG analysis aligns with previous analyses of MDS samples. These pathways are particularly intriguing. A substantial proportion of genes correlating with LDLRAD4 have yielded significant research outcomes in immunity. This suggests the promise of focusing on immunity in this direction for future research endeavors.

To conduct a more comprehensive analysis of LDLRAD4’s function and mechanism, GSEA was again employed to analyze the high and low expression groups of LDLRAD4 in terms of GO and KEGG. Regarding GO, the high-expression group was enriched in GOBP DNA Templated and DNA Replication, GOBP Mitotic Sister Chromatid Segregation, and others; the low-expression group showed enrichment in GOBP Hemostasis, GOBP Regulation of Body Fluid Levels, and others (Figure 4A). In KEGG analysis, the high-expression group was enriched in Cell Cycle, DNA Replication, Oocyte Meiosis, and others; while the low-expression group showed enrichment in Drug Metabolism Cytochrome P450, Focal Adhesion, MAPK Signaling Pathway, and others (Figure 4B).

Many GO analysis results, including DNA replication, chromosome segregation, and mitosis, were found to be associated with cell proliferation processes. Results from KEGG analysis suggest that LDLRAD4 may influence cell cycle-related pathways, corroborated by the GO analysis. Thus, it is hypothesized that LDLRAD4 could regulate the cell cycle through its impact on DNA templating and replication processes. Additionally, pathways aligning with the GEVA analysis results were identified, including those related to primary immunodeficiency disease and cell cycle. This finding further underscores the significance of LDLRAD4 within these pathways.

Given the close association of MDS with immunity in terms of LDLRAD4’s function and pathway enrichment, the CIBERSORT algorithm was employed to investigate the infiltration abundance of 22 immune cell types in the training set samples. Comparisons between the MDS and normal groups revealed significant differences in the abundance of T cells CD4 memory resting and Macrophages M1 (P < 0.05) (Figure 5A). To elucidate the complex interactions among immune cells, correlations between various immune cell types were investigated. Findings indicated a strong positive correlation between M1 and M2 macrophages, and a significant negative correlation with resting dendritic cells. Additionally, the strongest positive correlation was found between T cells CD4 memory resting and activated NK cells, and the most significant negative correlation with resting mast cells (Figure 5B).

Furthermore, ssGSEA was utilized to analyze the degree of immune cell activation and function in the high and low expression groups of LDLRAD4. Results revealed significant elevations in B cells, pDCs, Th1 helper cells, TIL, and inflammation-promoting in the high-expression group of the MDS group. Conversely, APC_co_inhibition and T helper cells showed significantly higher expression in the low-expression group (Figure 5C).

The relationship between infiltrating immune cells and LDLRAD4 was evaluated, revealing positive correlations with activated NK cells, T cells CD4 memory resting, and memory B cells, and negative correlations with Plasma cells, CD8 T cells, Eosinophils, resting NK cells, and Tregs (Figure 5D). These findings complement the conclusions from GO, KEGG, and co-expression analyses of LDLRAD4, enhancing the understanding of its association with immunity. Additionally, they offer crucial insights into the mechanisms of immune dysfunction and disease progression in MDS patients.

The analysis indicates that LDLRAD4 is significantly downregulated in MDS patients. To further investigate the impact of LDLRAD4 on the biological functions of MDS cell lines, we selected SKM-1 as our in vitro cellular model. We established a control group (pcDNA3.1-NC) and an overexpression group (pcDNA3.1-LDLRAD4). Western blotting and qRT-PCR revealed that the LDLRAD4 levels in the overexpression group were significantly higher compared to the control group (Figures 6A,B). Initially, we employed the CCK8 assay to evaluate the impact of LDLRAD4 on the viability of SKM-1 cells. The results indicated that overexpression of LDLRAD4 in SKM-1 cells significantly inhibited their proliferative capacity compared to the control group (Figure 6C). Subsequently, using flow cytometry, we found that compared to the control group, the overexpression group had cells arrested in the S phase, with a reduction in G2/M phase cells, indicating inhibited cell proliferation (Figure 6D). Finally, flow cytometry was used to assess apoptosis, and the overexpression group showed a significantly higher rate of cell apoptosis compared to the control group (Figure 6E). These findings suggest that LDLRAD4 may be a gene that inhibits the proliferation of MDS cells and promotes apoptosis.