Diabetic retinopathy (DR) is a major cause of blindness in adults^[1]. The global diabetic patients’ number was estimated to be 536.6 million in 2021, and is expected to rise to be 783.2 million in 2045^[2]. Approximately one-third of patients develop DR^[3]. Therefore, developing a reliable risk assessment method for the development of DR to facilitate early intervention and prevention of the disease is necessary.

Previous development of risk assessment models for DR relied on classical statistical methods to integrate multiple risk factors associated with the condition^[4,5]. However, these risk-prediction methods for DR exhibit limited efficacy due to their suboptimal predictive performance. In recent years, machine learning (ML) methods have demonstrated powerful self-learning abilities and improved disease prediction accuracy^[6]. Previously reported ML prediction models primarily include demographic factors^[7,8].

Genetic factors are widely recognized as important contributors to the progression of DR^[9-11]. However, only a few studies have employed ML methods that incorporate identified genetic risk factors. In the case of age-related macular degeneration, the inclusion of genetic factors enhances the accuracy of predictive models^[12,13].

Thus, the present study aimed to investigate the effect of an ML prediction model for DR while including the genetic risk factors identified in a DR pedigree using whole-exome sequencing (WES).

The DR pedigree, including probands and the type 2 diabetes mellitus (T2DM) cohort, was selected from the Desheng Diabetic Eye Study, the details of which have been previously described^[14]. The study protocol was approved by the Ethics Committees of Beijing Chao-Yang Hospital and Beijing Tongren Hospital and adhered to the tenets of the Declaration of Helsinki. Written informed consent was obtained from all participants before enrolment in the study. Diabetes was defined as a self-reported history of physician-diagnosed T2DM treated with insulin, oral hypoglycemic agents, or diet only; or by a fasting plasma glucose (FPG) concentration of 7.0 mmol/L (126 mg/dL) or more in at least two previous examinations; or a random plasma glucose concentration of ≥ 11.1 mmol/L (200 mg/dL). An overview of the study design is shown in Figure 1.

Figure 1.  Overview of the study design. The probands in the diabetic retinopathy (DR) pedigree and the type 2 diabetes mellitus (T2DM) cohort were derived from the Desheng Diabetic Eye Study. We collected clinical data through questionnaires and obtained venous blood samples from all participants. Genomic DNA was extracted from leukocytes. We identified seven single nucleotide polymorphisms (SNPs) and mutations associated with DR in the pedigree using whole-exome sequencing (WES). The dataset’s participants were genotyped for these SNPs employing Sequenom MassARRAY technology (Bioyong Technologies, Beijing, China). For the T2DM cohort, 70% of the data was allocated for training and 30% for validation. We assessed various machine learning algorithms, including the linear support vector machine (SVM), random forest, logistic regression, K-nearest neighbor, decision tree, multi-layer perceptron, linear discriminant analysis, multinomial naive Bayes, and AdaBoost. Model efficacy was evaluated by accuracy, sensitivity, specificity, the area under the receiver operating characteristic curve (AUC), and the kappa statistic.

All the participants underwent a standardized evaluation consisting of a questionnaire, ocular and anthropometric examinations, and laboratory tests. Laboratory tests, including fasting plasma glucose (FPG), glycosylated hemoglobin A1c (HbA1c), creatinine, uric acid, lipid profile (levels of total cholesterol, triglycerides, and high-density and low-density lipoprotein cholesterol), and urine albuminuria (UMALB) were conducted for all 420 participants. Venous blood samples were obtained from all participants, stored at −80 °C until DNA extraction.

The DR family was selected based on the following criteria: (1) more than two generations with T2DM members, (2) three or more DR members, (3) inclusion of at least one diabetic-without retinopathy (DWR) patient who had T2DM for more than four years, (4) healthy members older than 45 years, and (5) the fundus photograph of all included family members should be graded. Sixty Chinese families with DR were recruited for the study. One family was selected for further analysis (Table 1). The family pedigree is shown in Figure 2. Patients with one non-T2DM, one DWR, and four DR of the second-generation DR pedigree were selected for whole exome sequencing (WES).

Figure 2.  Pedigree of the family showing diabetes status of each member. Circles represent female participants and squares male participants. A slash through the symbol indicates that the family member is deceased. Black symbols indicate patients with diabetic retinopathy (DR). Dashed symbol represented individuals with diabetes mellitus (DM) but without DR.

The ML dataset was collected from the Desheng Diabetic Eye Study cohort. Based on the duration of diabetes and the grading of fundus photographs, patients from the cohort were assigned to the DWR group if they had more than 10 years of T2DM with no signs of DR (microaneurysms, hemorrhages, and exudates) or if they had more than 15 years of T2DM with fewer than five microaneurysms. Patients with five or more microaneurysms in at least one eye were assigned to the DR group. Patients who did not meet the DWR or DR criteria were excluded from the study. The duration of diabetes was defined as the interval between the first diagnosis of diabetes and the time of enrollment in the present study. One hundred ninety-eight patients with DWR and 222 patients with DR were randomly recruited.

Genomic DNA was isolated from the leukocytes using a TIANamp Blood DNA Kit (Tiangen Biotech, Beijing, China). DNA concentration was measured using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, United States). All DNA samples passing initial quality checks at a concentration of ≥ 50 ng/μL were collected for further sequencing.

Genomic DNA from six selected individuals from the DR pedigree was subjected to WES using a commercial service provided by the Beijing Genomics Institute, Shenzhen, China.

Exonic DNA libraries were enriched by hybridizing genomic DNA samples with the NimbleGen 2.1 M-probe sequence capture array^[15]. The libraries were sequenced on an Illumina HiSeq 2000 platform, generating 90-bp reads for each sample. Sequence reads were aligned to the human genome using the SOAP Aligner (v2.21) software. SOAPsnp software generated genotypes for each position in the target region. High confidence single nucleotide polymorphism (SNP) sites (quality ≥ 20 and depth ≥ 8X) were selected for further analysis.

SNP sites, coding DNA sequences, and indels were identified using SAMtools and BCF tools. Variants with a minor allele frequency (MAF) lower than 1% or absent in the East Asian population from gnomAD and EXAC databases were retained.

The Polyphen2 uses sequence- and structure-based predictive algorithms, and predicts if a change is damaging (> 0.85), possibly damaging (0.85–0.15) or benign (< 0.15). The SIFT tool mainly assesses a change being damaging (≥ 0.95) or tolerated (< 0.95), through analyzing the conserved sequence in different species. Variants identified by WES in the Chinese family with DR were validated by polymerase chain reaction (PCR) and Sanger sequencing using specific primers. Mutations detected in the DR family were subjected to co-segregation analysis to determine their causality.

Study participants in the dataset were genotyped for SNP sites using the Sequenom Mass ARRAY (Bioyong Technologies, Beijing, China). Quality control measures included MAF < 0.05, genotyping success < 80%, and failed Hardy-Weinberg equilibrium (HWE) test in control samples (P < 0.001, Pearson goodness-of-fit test).

Feature selection was performed using least absolute shrinkage and selection operator (LASSO) regression algorithms. The model training features included specific SNP loci and ten clinical variables: age at examination, biological sex (male or female), FPG, HbA1c, UMALB, insulin usage, oral glucose tolerance test (OGTT) results, therapeutic approach, monthly income, and educational level. Age was recorded on the day of the examination, and sex was recorded as assigned at birth. Insulin use and OGTT results were categorized into two groups (usage/no use and tested/not tested, respectively). The therapeutic approaches were classified into three groups: dietary control, exercise, and pharmacotherapy. Monthly income was divided into six brackets ranging from no income to more than 4,000 Chinese Yuan (CY). Educational level was divided into junior school or lower and high school or higher.

We utilized Python 3.9 and the SKLearn ML package for model development. The feature contributions were assessed using the Shapley Additive explanation (SHAP) approach. A linear support vector machine (SVM) served as the primary algorithm, partitioning the data into 70% training and 30% validation sets. We also examined additional machine learning algorithms, including random forest, logistic regression, K-nearest neighbor, decision tree, multilayer perceptron, linear discriminant analysis, multinomial naive Bayes, and AdaBoost. The model efficacy was measured using accuracy, sensitivity, specificity, area under the receiver operating characteristic curve (AUC), and kappa statistics.

Statistical analyses were performed using R (version 4.2.2) and R Studio (version 2023.06.1). The package used was EpiCalc (version 2.8.1.1, 2009). Chi-square test or Fisher’s exact test was used to compare categorical data between the two groups. Parametric variables were analyzed using t-tests. Nonparametric data were assessed using the Kruskal-Wallis rank-sum test. Binomial logistic regression was used for the multivariate analysis. The odds ratio (OR) and 95% confidence intervals (CI) were calculated. The statistical significance level was set at P < 0.05.

Within the DR pedigree comprising ten members, we identified seven SNP sites and mutations (rs116911833 in TRIM7, 1997T>C in LRBA, 1643T>C in PRMT10, rs117858678 in C9orf152, rs201922794 in CLDN25, rs146694895 in SH3GLB2, and rs201407189 in FANCC). These genetic variations were present in DR-afflicted individuals, but absent in DWR or non-T2DM subjects. The characteristics of these SNP sites and mutations, such as their location in the coding DNA sequence, SNP sites, allele frequencies, and predictive impact scores from SIFT and PolyPhen2, are detailed in Table 2.

The percentages of calorie sources (carbohydrate, fat, and protein) were analyzed using multiple logistic regression. After adjusting for multiple factors, including sex, ethnicity, use of insulin, glycosylated hemoglobin, hypertension, and exercise, dietary vitamin E [OR (95% CI): 0.97 (0.95, 1.00), P = 0.036] and selenium [OR (95% CI): 0.98 (0.96, 1.00), P = 0.017] intakes were inversely associated with the risk of DR (Table 2).

The demographic and clinical features selected for ML training in the T2DM cohort of 420 patients are shown in Table 3. These features, including age, FPG, HbA1c, UMALB, monthly income, insulin use, therapeutic interventions, and educational level, demonstrated statistically significant differences between the DR and non-DR groups (P < 0.05, T-test, Fisher’s exact test, and Kruskal-Wallis rank-sum test, respectively). Conversely, sex showed no significant difference (P = 0.65, chi-squared test). SNP rs146694895 exhibited significant variation (P = 0.02, Fisher’s exact test), whereas rs201407189 did not (P = 0.82, Fisher’s exact test). All features have all been selected through LASSO regularization, ensuring no collinearity issues are present among them; the lambda figures are shown in Figure 3.

Figure 3.  The lambda figures (coefficients and MSE) of Lasso features selection procedures.

During the evaluation of different ML algorithms, the SVM algorithm outperformed the others, achieving the highest accuracy, sensitivity, specificity, AUC, and kappa values for the test set (80.2%, 83.3%, 76.7%, 80.0%, and 60.1%, respectively), as presented in Table 4.

The rs146694895 and rs201407189 were observed enhanced the performance of ML models for DR prediction (Table 5). Other identified SNP sites and mutations did not contribute significantly to predictive accuracy and were thus excluded from the feature set.

The final model for predicting DR incorporated both demographic data and selected SNP sites (rs146694895 and rs201407189). The contribution of each parameter as quantified by the SHAP values is shown in Figure 4.

Figure 4.  SHAP analysis illustrating the importance of each feature in the proposed model. A: Variable importance plot displaying the impact of various predictors on the model’s output. B: Bar chart showing the mean decrease in the magnitude of each variable’s importance. HbA1c, glycosylated hemoglobin; UMALB, urine microalbumin; OGTT, oral glucose tolerance test; FPG, fasting plasma glucose; therapy, therapeutic approaches.

The performance of the SVM model is depicted through a confusion matrix and receiver operating characteristic curve in Figures 5. To facilitate the practical application of this model, we developed a web-based application for DR prediction that is accessible on DR Prediction Model Web Server. This tool allows users to input 12 features into designated fields and receive DR predictions (Figure 6).

Figure 5.  (A) Confusion matrix of the optimal model. “0” indicates diabetes without retinopathy (DWR), and “1” refers to diabetic retinopathy (DR). (B) The receiver operating characteristic curve of the proposed model. The area under the receiver operating characteristic (ROC) curve was 80.0%.

Figure 6.  Public internet calculator for diabetic retinopathy (DR) risk assessment using 12 features. The DR prediction model can be accessed via the public web server at DR Prediction Model, which utilizes 12 distinct features for risk calculation. Users can input 12 features to predict DR risk. The interface provides text boxes for the following variables: SH3GLB2 (SH3GLB2 rs146694895: GG = 0 and GA = 1), FANCC (FANCC rs201407189: GG = 0 and AG = 1), age (continuous variable), sex (biological sex: male = 1, female = 2), FPG (fasting plasma glucose [mmol/L], continuous variable), HbA1c (glycosylated hemoglobin [%], continuous variable), UMALB (microalbuminuria [mg/L], continuous variable), monthly income (no income = 1, basic living allowance = 2, < 2,000 Chinese Yuan [CY] = 3, 2,000–3,000 CY = 4, 3,000–4,000 CY = 5, and > 4,000 CY = 6), OGTT (oral glucose tolerance test: tested = 1 and not tested = 2), therapy (therapeutic approach: dietary control = 1, exercise = 2, and drug therapy = 3), insulin usage (with usage = 1 and without usage = 2), education (educational level: junior school or lower = 0, high school or above = 1). Upon submission, the model computes the results, and the resulting probability indicates the DR status of the subject.

In the present study, we identified seven novel SNP sites in seven genes associated with DR, using a DR pedigree. Subsequently, an ML model for predicting DR, including both demographic and genetic factors, was established and tested. The inclusion of rs146694895 in SH3GLB2 and rs201407189 in FANCC improved model performance.

The identified 7 SNP sites and mutations have not been reported previously. Multiple candidate genes have been proposed to be associated with the development of DR^[9-11]. However, previous studies involving linkage analyses, candidate gene studies, and underpowered genome-wide association studies have found limited or inconsistent genetic signals that cannot be replicated across multiple cohorts^[16-18]. The heritability of DR attributed solely to common genetic variants is estimated to be 7%, based on studies involving distantly related individuals^[10]. However, heritability estimates from specific family studies ranged from 18% to 52%^[19,20]. Thus, identifying DR-associated genes in the DR pedigree is more efficient. Therefore, we first attempted to identify DR-associated genes in a DR pedigree using WES.

To date, only a few studies have reported on the utilization of machine learning analysis for assessing DR risks^[7,8]. Moreover, these studies incorporated only the common risk factors associated with DR into their models. This study aimed to develop a robust DR risk prediction model by integrating genetic factors from the DR pedigree with demographic data.

SH3GLB2 rs146694895 and FANCC rs201407189 were identified as significant contributors to the model efficiency. These two genes were first reported to be associated with DR. Both SH3GLB2 and FANCC are expressed in the retina^[21,22]. Previous studies have shown that SH3GLB2 and SH3GLB1 play crucial roles in mitochondrial autophagy^[21]. Disruption of these genes can lead to weaker degradation of damaged mitochondrial membranes^[23]. In diabetes, problems with mitochondria can lead to the overproduction of harmful substances such as superoxide, which can trigger unhealthy cellular reactions, contributing to oxidative stress, a known factor in DR development^[24]. Furthermore, SH3GLB1 has been implicated in lipid metabolism in the adipose tissue and insulin resistance^[25]. Because SH3GLB2 can interact with SH3GLB1 and share significant sequences and similar structures, they may have a similar function. However, direct evidence linking SH3GLB2 with DR is lacking.

Mutations in FANCA and FANCC have been identified in more than 70% of patients with Fanconi anemia patients worldwide^[26]. Approximately half of patients with Fanconi anemia present with abnormalities in glucose metabolism^[22], and FANCC knockout mice display a higher propensity towards developing diabetes mellitus^[27]. Both patients with Fanconi anemia and knockout mouse models show increased levels of reactive oxygen species^[28,29]. This surge contributed to oxidative stress and triggered endothelial cell death. Endothelial dysfunction or cell death can cause vascular barrier disruption, vasoconstriction, diastolic dysfunction, vascular smooth muscle cell proliferation, migration, inflammatory responses, and thrombosis. Each outcome is closely associated with DR^[30].

In this study, the LASSO shrinkage and selection operator regression algorithm was used for feature selection. It has been proven to automatically select significant features that are relevant to the target variable, while setting the coefficients of irrelevant or redundant features to zero^[31]. Oh et al.^[7] embarked on a pioneering exploration of the potential of machine learning in the development of a risk prediction model for DR. This study demonstrated that the LASSO algorithm outperformed conventional indicators, delivering superior AUC values and accuracy rates. Using a similar sample size, our study achieved an improvement (5%) in accuracy. This may be attributed to the incorporation of both genetic factors and different ML algorithms.

Compared with prior studies, our research employed nine machine learning algorithms, facilitating flexibility in exploring diverse perspectives and techniques in model construction. Commonly used machine learning algorithms include the k-nearest neighbor, support vector machine, decision tree, random forest, linear discriminate analysis, AdaBoost, XGBoost, and CATBoost^[32]. Our study adopted all these algorithms and found that the SVM demonstrated the highest effectiveness, with superior AUC and accuracy scores. Hemoglobin A1c, duration of diabetes, cholesterol, and blood pressure only accounted for 10%–11% of the DR risk^[33,34]. Li et al.^[8] conducted an exhaustive analysis of the electronic medical records of a large cohort of 32,452 T2DM inpatients. Logistic regression, random forest, SVM, and XGBoost were adopted. The XGBoost algorithm exhibits the highest AUC for the validation set. This model relied on hospital data, which may introduce bias and limit its generalizability.

We also found that insulin usage, HbA1c, UMALB, age, FPG, OGTT, sex, therapeutic approach, monthly income, and educational level also influenced DR conditions. Previous studies have demonstrated that HbA1c, duration of diabetes, cholesterol, and blood pressure account for 10% to 11% of DR risk^[33,34]. The heritability of DR attributed solely to common genetic variants is estimated to be 7%, based on studies involving distantly related individuals^[10]. According to previous studies, the contribution of demographic data to DR risk is higher than that of heritability, which is consistent with our findings. We noticed that the contribution of SNP factors was less than that reported previously. The reason for this might be that DR is a multifactorial disease influenced by many genes, but we only included two SNP sites. Therefore, additional DR-associated genes should be included in prediction models.

The present study has some limitations. First, although our study marks an important step forward in the field, we acknowledge that the limited sample size used for model construction poses constraints on the overall accuracy of the model. Second, although we found that genetic factors identified in the DR pedigree played a crucial role in the performance of the ML model, the collection of adequate pedigree can pose a challenge due to the relatively advanced age at the onset of T2DM. Furthermore, although whole exome sequencing is a powerful tool, it is not as comprehensive as whole-genome sequencing for detecting a wide range of intronic sites in genes. This could potentially overlook some factors influencing gene function. Additionally, the absence of pathway validation for the relevant SNP sites constrains our capacity to substantiate the impact of these SNP sites on DR. Finally, the restricted sample size for model construction limits the model accuracy.

In conclusion, this study represents a crucial advancement in DR risk prediction modeling through the integration of multiple machine learning algorithms and comprehensive evaluation metrics. SNPs of SH3GLB2 and FANCC screened from the DR pedigree were instrumental in enhancing the performance of the DR prediction model, thus proving the effectiveness and accuracy of this screening method. Our study is the first to determine the importance of genetic factors in the DR risk prediction model. Further enrollment of other DR-associated genetic sites in the prediction model is required.