ScRNA-seq data from five proliferative diabetic retinopathy (PDR) patients (GSE165784) (12) and three control samples (15) were processed alongside two bulk RNA PDR patient cohorts from the GEO database (cohort 1: GSE160306 (16), 76 samples; cohort 2: GSE102485 (17), 25 samples). Oxidative stress-related genes (ORGs) were identified from the Gene Ontology and PathCards databases.

The single-cell data of five PDR patients (12) and three healthy controls (15) were analyzed via Seurat (18). We filtered cells based on mitochondrial content (<10%), cell count (>300), and gene number (1000-5000). The t-distributed stochastic neighbor embedding (t-SNE) (19) and ‘RunHarmony’ functions (20) facilitated visualization and batch effect correction. Cell subtypes were annotated according to cell markers from the original study (12, 15). In the differential expression analysis between microglia and mesenchymal cells in PDR versus control samples, the mitochondrial and ribosomal genes were removed. We used the Wilcoxon signed-rank test to identify significant genes (adjusted P value <0.05, absolute log₂FC >1).

NMF analysis, specifically consensus NMF (cNMF), was applied to microglia in PDR samples, standardizing negative values to zero. After more than 100 iterations, we explored the components (k) ranging from 2 to 10 signatures, determining the optimal component number via a diagnostic plot from the provided tutorial (https://github.com/dylkot/cNMF) (21). A two-step gene ranking algorithm was used to identify nonoverlapping gene modules, which were further analyzed for expression patterns using Pearson correlations and hierarchical clustering, revealing three distinct meta-programs.

Twelve machine learning algorithms, including (1) Least Absolute Shrinkage and Selection Operator (LASSO), (2) Ridge, (3) Elastic network (Enet), (4) Stepglm, (5) Support Vector Machines (SVM), (6) GlmBoost, (7) Linear Discriminant Analysis (LDA), (8) Partial Least Squares Regression for Generalized Linear Models (plsRglm), (9) Random Forest (RSF), (10) Generalized Boosted Regression Models (GBMs), (11) XGBoost, (12) Naive Bayes, were utilized to develop a predictive ORG model for PDR. We constructed 109 model combinations, trained initial models with the GSE160306 cohort and validated them with the GSE102485 cohort. Model performance was assessed using the AUC.

We utilized AlphaFold 2, a tool that has achieved remarkably accurate levels comparable to those obtained through human observation via advanced techniques such as cryoelectron microscopy, for the prediction of protein structures (22). For our specific study objectives, we used AlphaFold 2 to predict the structures of select proteins relevant to our research. We focused on the proteins ALKBH1, PSIP1, and ATP13A2, which play significant roles in the context of PDR. The sequences of these proteins were meticulously retrieved from the NCBI database (23).

To investigate the binding affinities and interaction patterns of the drug candidates with their targets, we utilized AutoDock Vina 1.2.2, a software designed for in silico protein–ligand docking (24). The molecular structure of resveratrol was obtained from the PubChem Compound database (https://pubchem.ncbi.nlm.nih.gov/) (25). AlphaFold 2 was used to generate the 3D coordinates for ALKBH1. Before docking analysis, all proteins and ligand files were prepared by converting them into PDBQT format. This preparation involved the removal of water molecules and the addition of polar hydrogen atoms to ensure accurate docking simulations. The docking grid box was strategically positioned to encompass the target protein’s domain, allowing for unhindered molecular movement within the simulation. The dimensions of the grid box were set to 30 Å × 30 Å × 30 Å, with a grid point spacing of 0.05 nm to capture detailed interaction data. The molecular docking studies were conducted using AutoDock Vina 1.2.2 (http://autodock.scripps.edu/).

In our study, we investigated protein interactions involving three pivotal ORGs. The use of the STRING database (https://string-db.org/) (26), a comprehensive resource, enabled us to compile and amalgamate data on protein−protein interactions (PPIs). Our focus was directed toward interactions with confidence scores surpassing 0.7, a threshold set to ensure the biological relevance and significance of these interactions.

To deepen our analysis and improve its visualization, we transferred the relevant data into Cytoscape (version 3.8.2) (27). Within the Cytoscape environment, we leveraged the capabilities of the cytoHubba plugin. This allowed us to pinpoint and rank the top 10 nodes in the PPI network. The ranking process utilized seven distinct algorithms, each contributing a unique perspective to the analysis. These algorithms included the following: Radiality, which measures the centrality of a node; Maximum Neighborhood Component (MNC), which assesses the largest connected component around a node; Maximum Clique Centrality (MCC), which focuses on the largest clique a node belongs to; Edge Percolated Component (EPC), which evaluates the connectivity and clustering; DMNC, which is the Maximum Neighborhood Component Centrality, a derivative of the MNC; Degree, which counts the number of edges linked to a node; and Closeness, which measures the average distance to other nodes. To synthesize and present our findings, we utilized an UpSet diagram.

To identify the hub genes associated with PDR, we used the Comparative Toxicogenomics Database (CTD, http://ctdbase.org/) (28). Utilizing the CTD, we conducted an in-depth analysis to unravel the connections between potential key genes and a spectrum of relevant conditions. This included not only PDR but also a broader scope of related health issues, such as other eye diseases, retinal disorders, vascular diseases, complications arising from diabetes, and diabetes mellitus itself.

All the statistical analyses of single cells were performed with R (version 4.3.1). A P value less than 0.05 was considered to indicate statistical significance if not explicitly stated.

In this study, we conducted an in-depth analysis of single-cell RNA sequencing data from five patients with proliferative diabetic retinopathy (PDR) (GSE165784) (12) and three healthy controls (15), implementing t-distributed stochastic neighbor embedding (t-SNE) for visualization post-quality control and data normalization. This approach effectively distinguished between cellular clusters of the PDR and control cohorts.

Figure 1A displays the range and individual RNA counts per cell, reflecting successful quality control measures for our sample analysis. We highlighted the 2000 genes with the highest variability across samples in Figure 1B . To further dissect this complexity, we applied linear dimensionality reduction to compute principal components (PCs), as illustrated in Figure 1C . The determination of significant PCs was aided by the integration of both ElbowPlot ( Figure 1D ) and JackStrawPlot ( Figure 1E ), setting the stage for more nuanced analyses. The distribution of cells across the PDR and control groups is presented in Figure 1F , with the study encompassing 5 PDR and 3 control samples, as depicted in Figure 1G . A total of 26 clusters were identified across the samples ( Figure 1H ). Through marker analysis from the original study (12), we classified cells into five primary types: microglia, lymphocytes, myeloid cells, endothelial cells, and mesenchymal cells ( Figure 1I ). This categorization was further refined, resulting in the identification of 14 distinct subcell types ( Figure 1J ).

Our investigation of differential gene expression and intercellular communication within the retinal microenvironment of PDR patients highlighted important findings. We observed pronounced ligand-receptor interactions among various cell types, with notable interactions between microglia and mesenchymal cells, as well as between microglia and endothelial cells ( Figures 2A, B ). These interactions shed light on the intricate signaling pathways that could be instrumental in the development and progression of PDR, suggesting potential therapeutic targets.

In a detailed analysis of gene expression between microglia and mesenchymal cells in PDR versus control samples, 40 genes were upregulated, and 111 genes were downregulated in both cell types ( Figures 2C, D ; Supplementary Tables 1 – 4 ). The upregulated genes, including FN1, ATP5F1E, B2M, MALAT1, and ATP5MG, and downregulated genes, such as ATP5E, ALDOA, ATP5L, ATP5I, and C14orf2, indicate a complex regulatory landscape. Furthermore, we revealed nuanced gene expression patterns: GLUL, DAB2, SELENOP, and ALDH1A1 were downregulated in mesenchymal cells but upregulated in microglia ( Figure 2E ; Supplementary Tables 5 , 6 ), while C12orf75, ITM2C, and CCND1 showed the opposite pattern ( Figure 2F ; Supplementary Tables 7 , 8 ).

In our detailed investigation of specific microglia within PDR samples, we employed the sophisticated technique of NMF to determine the unique transcriptional landscape of these cells. This advanced approach allowed us to systematically catalog various gene modules, which are fundamentally crucial in defining the distinct states of cells. Through this meticulous process, we were able to identify and analyze patterns of gene coexpression within individual PDR samples.

Our comparative analysis across multiple PDR samples was instrumental in revealing recurring gene modules. This aspect of our study was particularly significant because it effectively minimized the impact of technical variations, thereby enhancing the reliability and accuracy of our findings. By focusing on these gene modules, we gained valuable insights into the transcriptional intricacies inherent in PDR microglia.

One of the most noteworthy outcomes of our analysis was the identification of three distinct meta-programs. These meta-programs were discerned and clustered based on their correlation coefficients, providing a clear representation of the transcriptional synergy within microglia. The top-scoring genes of these meta-programs were characterized, as depicted in Figure 3A . Notably, Meta-Program 1 emerged as particularly prominent, exhibiting the highest level of correlation among the three. These results suggest that the genes within Meta-Program 1 are potentially central to the transcriptional identity and function of microglia in the context of PDR.

Elucidation of these meta-programs is important to our understanding of PDR. This study provides a novel perspective on the transcriptional dynamics of microglia, a critical component of disease pathology. This insight not only enhances our understanding of the molecular mechanisms underlying PDR but also opens up new directions for targeted therapeutic strategies aimed at modulating these specific transcriptional programs.

To develop a predictive model for PDR based on ORGs, we examined a subset of 15 genes that notably intersected within Meta-Program 1, as illustrated in Figure 3B . We constructed ORG models by twelve diverse machine learning algorithms, including (1) Least Absolute Shrinkage and Selection Operator (LASSO), (2) Ridge, (3) Elastic network (Enet), (4) Stepglm, (5) Support Vector Machines (SVM), (6) GlmBoost, (7) Linear Discriminant Analysis (LDA), (8) Partial Least Squares Regression for Generalized Linear Models (plsRglm), (9) Random Forest (RSF), (10) Generalized Boosted Regression Models (GBMs), (11) XGBoost, (12) Naive Bayes ( Figure 4A ). Among the 109 models constructed, the cream of the crop emerged in the form of models based on a sophisticated Stepglm [backward]+RF approach. These standout models prominently featured three key ORGs, ALKBH1, PSIP1, and ATP13A2, as delineated in Figure 4B . The importance of this model was unmistakably demonstrated in the training cohort (GSE160306), in which an outstanding area under the curve (AUC) of 0.989 was achieved. This exceptional level of predictive accuracy underlines the model’s formidable potential as a tool for diagnosing PDR.

To further validate the model’s applicability in a clinical setting, we undertook a validation study using an external cohort (GSE102485). The results were encouraging, as the model retained a significant level of diagnostic accuracy, as evidenced by an AUC of 0.833. This performance in an external cohort not only reinforces the model’s robustness but also underscores its potential utility as an early detection and ongoing monitoring tool for PDR.

In our study, we harnessed the ability of AlphaFold 2 technology to predict the complex structures of three pivotal ORGs, ALKBH1, PSIP1, and ATP13A2, as shown in Figure 4C . The predictive confidence for ALKBH1 and ATP13A2 was notably high, whereas PSIP1 demonstrated lower confidence levels and was subsequently not included in further analysis.

Furthermore, to assess the binding affinity of potential therapeutic agents for these targets, we conducted a molecular docking analysis. Specifically, we explored the interaction between ALKBH1 and the candidate drug resveratrol utilizing AutoDock Vina v.1.1.2 for this purpose. The analysis provided insights into the binding modes and calculated the binding energies for the interactions ( Figures 5A–C ). The derived binding energy for the ALKBH1-resveratrol complex was -6.471 kcal/mol, suggesting a highly stable interaction. This strong and stable binding affinity further underscores the potential therapeutic relevance of targeting ALKBH1 with resveratrol in the context of oxidative stress-related conditions.

Next, we explored the protein-protein interactions (PPIs) of these ORGs. For this purpose, we utilized the STRING database, which is renowned for its extensive protein interaction data. Our focus was on interactions with confidence scores exceeding 0.7, ensuring that only biologically significant and reliable interactions were considered. This selective approach was instrumental in sifting through vast data to identify meaningful connections that could be crucial in the context of PDR.

The PPI network enriched with these curated data was then intricately analyzed using Cytoscape. This platform enabled us to visualize and dissect the complex web of interactions. Using the cytoHubba plugin within Cytoscape, we systematically identified the top 10 nodes in the network utilizing a suite of seven sophisticated algorithms. These included Radiality, Maximum Neighborhood Component (MNC), Maximum Clique Centrality (MCC), EPC (Edge Percolated Component), DMNC (Maximum Neighborhood Component Centrality), Degree, and Closeness, each offering a unique lens to view and understand the network’s structure. Figures 6A, B depict these findings, revealing a comprehensive map of the interactions.

Furthermore, to emphasize the interconnected nature of these interactions, we constructed an UpSet diagram ( Figure 6C ). This visualization succinctly highlighted the convergence of hub genes across different algorithms, revealing key proteins such as H2AC8, H2BC12, H2AC13, and H2AC16 that were consistently central across all algorithms, as detailed in Supplementary Table 9 . This representation was instrumental in highlighting the core genes within the network, thereby elucidating their potential collective role in the pathophysiology of PDR.

To complement our protein interaction analysis, we utilized CTD as an instrumental resource. The CTD facilitated the expansion of our study to investigate the connections between our identified hub ORGs and a range of conditions associated with PDR, such as diabetic retinopathy, various eye and retinal diseases, vascular complications, and diabetes mellitus itself. Figures 7A–F display these connections, emphasizing the marked correlation between genes ALKBH1, PSIP1, ATP13A2, and the aforementioned conditions, validated by substantial reasoning scores within the CTD. Additionally, we incorporated a “negative control” gene, PXDNL, an unrelated ORG, to bolster the conclusiveness of our analysis.