The GEO database (https://www.ncbi.nlm.nih.gov/) was used to acquire mRNA expression profile of the GSE221521 and GSE189005 datasets with the GPL24676 and GPL23126 platforms, respectively. The study included 41 and 50 venous blood samples of DR and normal samples, respectively, in GSE221521 as the training set. GSE189005 consisted of 10 and 9 venous blood samples of DR and normal samples, respectively. The IEU OpenGWAS database (https://gwas.mrcieu.ac.uk/) was employed to obtain the GWAS ID and data for exposure factors and DR. DR was considered as the outcome (finn-b-DM_RETINOPATHY). The finn-b-DM_RETINOPATHY dataset comprised of 14,584 DR samples and 16,380,459 single nucleotide polymorphisms (SNPs).

In the GSE221521 dataset, DR-differentially expressed genes (DEGs) between DR and normal samples were identified conducted utilizing the DESeq2 package (v 1.34.0) (|Log2FC| > 0.5, p < 0.05) (10). The P value was corrected by Benjamini-Hochberg (BH) method. The volcano map and heat map were drawn employing the ggplot2 (v 3.4.1) and ComplexHeatmap (v 2.16.0) packages, respectively (11, 12).

The WGCNA was implemented to seek DR key module genes. The WGCNA package (v 1.71) was used to perform hierarchical clustering of all samples in the GSE221521 dataset, and outliers were removed (13). Then, the optimal soft threshold (β) was determined by realizing scale-free distribution and setting R² above 0.85. Based on the β value, the module was segmented by applying the standard of hybrid dynamic tree cutting algorithm (deepSplit = 2, mergeCutHeight = 0.3), with each module containing at least 100 genes. Modules with an absolute correlation greater than 0.6 with DR were further analyzed, and genes in these module were DR key module genes.

The ggVennDiagram package (v 1.2.2) was used to determine the overlap between DR-DEGs and DR key module genes, and candidate genes were identified if there was overlap (14). The clusterProfiler package (v 4.6.0) was used to perform functional annotation analysis, including the Gene Ontology (GO) functions, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and the disease ontology (DO) enrichment analysis (p < 0.05) (15). The STRING database (https://string-db.org/) was employed to predict interactions between proteins corresponding to candidate genes with a confidence score threshold of 0.4. The protein-protein interaction (PPI) network was subsequently displayed adopting Cytoscape (v 3.10.0) (16).

The candidate genes served as exposure factors, and DR was used as the outcome for MR analysis. In MR studies, the following three assumptions were made: (a) the presence of a significant correlation between instrumental variables (IVs) and exposure factors is imperative, (b) IVs should not be affected by confounding factors related to exposure factors or outcome, and (c) IVs can affect the outcome only through exposure factors.

The mv harmonize data function in the TwoSampleMR package (v 0.5.6) was used to unify effect alleles and effect sizes (17). Next, SNPs exhibiting significant correlation with candidate genes were selected as IVs (p<5×10^-8), and IVs for linkage disequilibrium (LD) were removed (clump=TRUE, R^{2 =} 0.001, kb=10000). The function extract instruments of TwoSampleMR package (v 0.5.8) was employed for this procedure (17). MR analysis of causality was carried out by five methods—MR Egger, Weighted median, Inverse variance weighted (IVW), Simple mode, Weighted mode), of which results of the IVW were the primary reference (p<0.05) (18–22). Scatter plots, forest plots, and funnel plots were prepared to visualize the results. An odds ratio (OR) greater than 1 indicated that the gene was a risk factor for DR, while the value was less than 1, the gene was considered a protective factor.

The reliability of MR analysis was assessed via sensitivity analysis, consisting of heterogeneity and horizontal pleiotropy tests, as well as Leave-One-Out (LOO) analysis. Initially, the Cochran’s Q test for heterogeneity was conducted, with p>0.05 indicating the absence of heterogeneity. Subsequently, a horizontal pleiotropy test was performed by employing the MR pleiotropy test function, with p>0.05 indicating that SNPs affected the outcome only through exposure factors. Finally, the LOO analysis was conducted using the MR leave one out function to determine whether a single SNP could significantly alter the overall effects. Genes exhibiting causal relationships with the outcome and passing the sensitivity analysis were identified as potential.

The pROC (v 1.18.0) (23) package was employed to plot ROC curves of potential biomarkers in the GSE221521 and GSE189005 datasets, and genes exhibiting area under the curve (AUC) values greater than 0.7 were considered reliable biomarkers. Subsequently, the biomarker expression levels in the GSE221521 and GSE189005 datasets were validated. Based on the expression of biomarkers, the rms package (v 6.5.0) was used to develop a diagnosis nomogram for DR patients (24). Calibration curves and ROC curves were constructed to evaluate the accuracy and reliability of the nomogram predictions. The calibration curve was plotted using the calibrate function and boot method.

Using c2.cp.kegg.v2023.1.Hs.symbols.gmt as the background gene set, GSEA of biomarkers was conducted with R clusterProfiler (v 4.6.0) (p<0.05) (15). Based on biomarker expression, all samples in the GSE221521 dataset were reorganized into high and low expression groups. The GSVA package (v 1.46.0) was used for GSVA of biomarkers in these two groups (25).

The single sample GSEA (ssGSEA) algorithm was used to assess the abundance of 28 immune cells in the samples of DR and control groups. The Wilcoxon test was used to compare the difference in immune cell infiltration between the two groups (p<0.05). The Spearman correlation coefficient between distinct immune cells, as well as between biomarkers and these cells was computed (|r|>0.3, p<0.05).

The miRNA targeting biomarkers were predicted using the miRNet database (https://www.mirnet.ca), while the transcription factors (TFs) regulating biomarkers were obtained from the networkanalyst database (http://www.networkanalyst.ca). Cytoscape (v 3.10.0) was then used to construct a regulatory network involving TF-mRNA-miRNA interactions (16). In order to explore potential drugs for the treatment of DR, targeted drugs for biomarkers were searched based on Drug Signatures Database (DSigDB) database (http://tanlab.ucdenver.edu/DSigDB), and results with p < 0.05 were selected.

Peripheral venous blood samples were collected from 20 patients in the Second affiliated Hospital of Kunming Medical University, including 10 patients with DR and 10 patients without DR. These samples were divided into two parts, each comprising of five DR samples and five control samples, and RT-qPCR was performed on each part to verify the screened biomarkers. This experiment was approved by the Institutional Review Board of the Second affiliated Hospital of Kunming Medical University (Review -PJ- Research -2024–134). TRIzol (Ambion, Austin, TX) was used to separate total RNA from 10 samples following the manufacturer’s instructions. The concentration of RNA was extracted using NanoPhotometer N50 (IMPLEN GmbH), and the purity of RNA was assessed by measuring the ratio of A260/A280. The SureScript-First-strand-cDRA-synthesis-kit (Servicebio, Wuhan, China) was used for reverse transcription of total RNA into cDNA as per the manufacturer’s instructions. The temperature was set to 25°C for 5 minutes, followed by 50°C for 15 minutes, then raised to 85°C for a duration of 5 seconds. Finally, maintain the temperature at 4°C during reverse transcription. Subsequently, qPCR analysis was performed using the 2xUniversal Blue SYBR Green qPCR Master Mix (Servicebio) according to the provided manual. The amplification conditions were 95°C for 1min, 95°C for 20s, 55°C for 20s, and 72°C for 30s. The primer sequences for PCR are listed in Supplementary Table 1 . Gene expression levels were normalized to GAPDH as an internal reference and calculated using the 2−ΔΔCq method (26). The characteristics of patients are listed in Supplementary Table 2.

R software (v 4.2.1) (https://www.R-project.org/.) was utilized to process and analyze data. Statistical analysis was performed using the wilcox.test method in R.

The differential expression analysis between DR and normal samples in GSE221521 yielded 2,283 DR-DEGs, comprising 1,304 upregulated DR-DEGs and 979 downregulated DR-DEGs. The volcano map illustrates both upregulated and downregulated DR-DEGs ( Figure 2A ). The heatmap displays the expression of the top 50 DR-DEGs, ranked based on p.adj values, in the two groups ( Figure 2B ).

Sample clustering analysis revealed two outliers (GSM6881234 and GSM6881340) in the GSE221521 dataset and were removed for subsequent analysis ( Figure 2C ). The β value was determined to be 13 by setting a threshold of scale-free R² above 0.85 to construct gene modules ( Figure 2D ). As a result, 10 modules were obtained by the hybrid dynamic tree cutting algorithm ( Figure 2E ). Subsequently, the correlation between the green (R=0.73), light green (R=0.66), and midnight blue (R=0.63) module genes with DR were greater than 0.6, and a total of 1,079 DR module genes were obtained from these modules ( Figure 2F ). Finally, 423 candidate genes were determined based on interaction between 2,283 DR-DEGs and 1,079 DR module genes ( Figure 2G ).

In order to understand the functions, diseases, and pathways of candidate genes, enrichment analysis was performed. The GO functions of candidate genes included semaphorin receptor complex, semaphorin-plexin signaling pathway involved in axon guidance, and regulation of cell shape ( Figure 3A ). The KEGG enrichment analysis showed that candidate genes were markedly enriched in endometrial cancer, basal cell carcinoma, fatty acid biosynthesis, and apoptosis ( Figure 3B ). The DO enrichment analysis showed that candidate genes were significantly related to meningioma, tuberous sclerosis, intracranial arterial disease, and cerebral arterial disease ( Figure 3C ).

The PPI network constructed to investigate the interactions among genes contained 324 nodes and 815 edges. AGO2 had a direct interaction with 18 genes (e.g., CLCN6, UBR4, SIDT2), whereas ACTG1 interacted with 14 genes (e.g., SRCAP, ITGA3, HCLS1) ( Figure 3D ).

The 423 candidate genes served as exposure factors, with 176 genes exhibiting SNPs, and DR was used as the outcome for the MR analysis. The IVW method was used to identify 13 candidate biomarkers (OSER1, HIPK2, DDRGK1, PCK2, IK, IRF5, COLGALT1, TRPM2, SLC38A10, TSNARE1, PAQR7, ZNF142, and ARID1A) that exhibited significant causal relationships with DR (p<0.05). Among them, eight candidate biomarkers (OSER1, HIPK2, DDRGK1, PCK2, TRPM2, SLC38A10, TSNARE1, and ZNF142) exhibited an OR greater than 1 and identified as risk factors for DR. Conversely, five candidate biomarkers (IK, IRF5, COLGALT1, PAQR7, and ARID1A) demonstrated an OR less than 1 and were considered protective factors for DR ( Table 1 ).

In the scatter plot, the effect of SNPs of OSER1, HIPK2, DDRGK1, PCK2,TRPM2, SLC38A10, TSNARE1, and ZNF142 on DR were positively correlated overall, while the effect of SNPs of IK, IRF5, COLGALT1, PAQR7, and ARID1A on DR were negatively correlated overall. These results confirmed the above conclusions ( Supplementary Figure 1 ). Consistent with the previous results, the forest map illustrated that the MR effect size of risk factors for DR exceeded 0 and the MR effect sizes of risk factor were less than 0, providing further evidence that IVs exhibit no or weak correlation with outcome ( Supplementary Figure 2 ). At last, the funnel plot illustrated that MR analysis of 13 candidate biomarkers and DR was consistent with Mendel’s second random law ( Supplementary Figure 3 ).

Sensitivity analysis results showed that all 13 candidate biomarkers passed the tests of horizontal pleiotropy and heterogeneity, and LOO analysis affirming the robustness and reliability of our MR analysis ( Supplementary Tables 3 , 4 , Supplementary Figure 4 ).

The AUC values in ROC curves of OSER1, HIPK2, and DDRGK1 were 0.868, 0.815, 0.806, respectively, in GSE221521, and 0.700, 0.833, and 0.722, respectively, in GSE189005 ( Figures 4A–F ). Given that all the AUC values exceeded 0.7, these three genes could serve as reliable biomarkers. The expression levels of OSER1, DDRGK1 and HIPK2 were consistent in GSE221521 and GSE189005 ( Figures 4G, H ). RT-qPCR results also demonstrated a significant high expression of HIPK2 in patients with DR, whereas OSER1 exhibited a significant lower expression level ( Figures 4I–K , Supplementary Figure S5 ). A nomogram was constructed based on the expression levels of OSER1, HIPK2, and DDRGK1 ( Figure 5A ). The AUC value for the nomogram was 0.942, and the calibration curve of the nomogram was almost a straight line, indicating a good predictive performance of the nomogram ( Figures 5B, C ).

GSEA analysis demonstrated significant associations between HIPK2 and DDRGK1 with ribosome, Parkinson’s disease, oxidative phosphorylation, and other pathways. OSER1 was enriched in spliceosome, neuroactive ligand receptor interaction, and olfactory transduction ( Figures 6A–C ). The top 20 pathways enriched by DDRGK1, HIPK2 and OSER1 in GSVA are shown in Figures 6D–F . The first three pathways (glycine serine and threonine metabolism, base excision repair, etc.) of DDRGK1 were enriched in the high expression group, while the remaining 17 pathways (such as inositol phosphate metabolism and glioma) were enriched in the low expression group. Similarly, the initial 11 pathways (lysine degradation, prostate cancer, etc.) of HIPK2 were enriched in the high expression group, whereas the remaining nine pathways (glycine serine and threonine metabolism, ribosome, etc.) were enriched in the low expression group. Among pathways related to OSER1, the first 14 (beta alanine metabolism, nitrogen metabolism, etc.) were enriched in the high expression group, and the remaining six (glyoxylate and dicarboxylate metabolism, primary immunodeficiency, etc.) were enriched in the low expression group.

The distributions of the 28 immune cells in the sample were visualized using a heatmap ( Figure 7A ). There were 10 immune cell types exhibiting notable disparities between the DR and control groups, such as central memory CD8⁺ T cells, activated CD8⁺ T cells, monocytes, effector memory CD8⁺ T cells, activated CD4⁺ T cells, and myeloid derived suppressor cells ( Figure 7B ). There was a significant positive association between activated CD8⁺ T cells and activated CD4⁺ T cells (R=0.60). Conversely, a significant negative correlation was identified between monocytes and activated CD4⁺ T cells (R=-0.52) ( Figure 7C ). DDRGK1 expression significantly and positively correlated with effector memory CD8⁺ T cells and negatively with plasmacytoid dendritic cells. Activated CD4⁺ T cells correlated significantly and positively with OSER1 expression. HIPK2 expression significantly and positively correlated with myeloid derived suppressor cells and negatively with activated CD8+ T cells ( Figures 7D–F ).

TF-mRNA-miRNA network was composed of 25 TFs, three biomarkers, and 191 miRNAs, with a total of 219 nodes and 233 edges ( Figure 8 ). The hsa-mir-92a-3p linked with all three biomarkers. GATA2 and RELA had strong associations with HIPK2 and OSER1, thereby highlighting their role in the pathogenesis of DR. The target drugs of HIPK2 mainly included GW5074 (Raf1 Kinase Inhibitor I) MRC, LY-317615 Kinome Scan, and GSK650394A MRC. Furthermore, OSER1 was primarily targeted by Thioguanosine PC3 UP, Elesclomol CTD 00004602, and Gedunin CTD 00003449 ( Table 2 ).