This study had two independent datasets. Dataset 1, a baseline dataset without follow-up data, comprised 690 patients with IgAN. These patients received the kidney biopsy in our center but returned to local for follow-up. Dataset 2, a follow-up dataset (PKU-IgAN cohort), included 1,808 patients with IgAN who were registered and with long-term follow-up in the Peking University First Hospital IgAN database from 1997 to 2020 (29). All patients with IgAN were diagnosed based on the histologic and immunofluorescence study of the renal biopsy, and those with <8 glomeruli per biopsy section were excluded (29). After excluding 243 patients without blood lipid data, 28 patients presented at younger than 16 years of age, and 14 patients presented acute kidney failure, 1,523 patients in dataset 2 were finally enrolled in this study, consisting of 1,168 patients with Oxford MEST-C scores and 355 patients lacking Oxford MEST-C scores.

Finally, a total of 690 patients in dataset 1 and 1,168 patients with Oxford MEST-C scores in dataset 2 were enrolled in our study as the modeling group, and 355 patients without Oxford MEST-C scores in dataset 2 were enrolled in this study as the external validation group ( Figure 1 ).

All clinical characteristics were collected at the time of the renal biopsy. The estimated glomerular filtration rate (eGFR) was calculated using the Chronic Kidney Disease Epidemiology Collaboration (CKD-EPI) formula (30). Renal biopsies were categorized according to established criteria for the Oxford MEST-C scoring system (24, 26, 31). Mean arterial pressure (MAP, mm Hg) was defined as diastolic pressure plus a third of the pulse pressure. The end-stage kidney disease (ESKD) was defined as eGFR <15 ml/min/1.73 m², dialysis, or kidney transplantation. Our study was approved by the Ethics Committee of Peking University First Hospital (IRB number 2020Y197). Written informed consent was provided by all participants.

The pathology T-score prediction (T _pre) model, constructed by the stacking algorithm, was used to predict whether IgAN patients would have T lesions (yes or no). The stacking algorithm is an integrated machine learning algorithm that can summarize several models and predict new observations. It utilizes the prediction of a collection of models as input for training a second-level model. This second-level model aims to find the best combination of the prediction of first-level models. Stacking can shield the capabilities of a range of well-performing models so that a better output prediction model can be achieved (32). In our study, we combined three machine learning algorithms, namely, support vector machine (SVM), random forest (RF), and logistic regression as first-level models, and then logistic regression as the second-level model to output the final probability of the binary T-score (with or without tubular atrophy/interstitial fibrosis, T _pre).

The input variables used in this model were chosen by AUCRF (33), a method using the random forest to find the optimal set for prediction. Variables entered into the AUCRF included age, sex, body mass index, systolic arterial pressure, diastolic arterial pressure, mean arterial pressure, hypertension, eGFR, proteinuria, microhematuria, history of gross hematuria, serum IgA, serum uric acid, serum triglycerides, total cholesterol, high-density lipoprotein, and low-density lipoprotein.

Several studies have demonstrated the value of tubular atrophy/interstitial fibrosis (T) in predicting ESKD in patients with IgAN (16, 19, 24, 34, 35). To evaluate whether the predicted T-score could help predict ESKD and how effective it was, we constructed a 5-year ESKD prediction model based on the XGBoost algorithm. To illustrate the significance of tubular atrophy/interstitial fibrosis in predicting ESKD, we first constructed a 5-year ESKD prediction model with only clinical variables as input variables (base model). Then, the 5-year ESKD prediction model using clinical variables and the real pathological T lesions score (T _bio, T0 was assigned 0, T1 and T2 were assigned 1) was also developed (base model plus T _bio) to evaluate the additive value of atrophy/interstitial fibrosis (T) in predicting ESKD. Finally, to evaluate whether the value of T _pre in predicting ESKD of patients with IgAN was consistent with real pathological T lesions (T _bio) when the base model plus T _bio was trained in the training set, the T _bio of the testing set was replaced by the corresponding T _pre predicted by the pathology T-score prediction model and then the testing set was used to evaluate the model performance (the base model plus T _pre). For the base model plus T _pre, the purpose of training the model using real pathological T-score (T _bio) was for the model to learn the true value of T for predicting ESKD.

XGBoost is a kind of ensemble of the decision tree, whose advantages include higher-order interactions and complex non-linear relationships between the model features and the outcome (21). It has been shown to achieve impressive performance in predicting renal failure risk and provide explanations for variables by ranking their importance (16, 34). We also applied other machine learning algorithms to our data set for evaluating whether the predicted T could be used in ESKD prediction models based on different algorithms, including RF, penalized regression, artificial neural network (ANN), and SVM.

Characteristics selected by the Cox proportional hazards model were collected at the time of the renal biopsy at enrollment [age, sex, systolic arterial pressure, diastolic arterial pressure, proteinuria, eGFR, serum IgA, serum uric acid, serum triglycerides, total cholesterol, low-density lipoprotein, and history of previous use of renin–angiotensin system (RAS) inhibitors and immunosuppressants as well as pathological T lesions], whereas the binary outcome (ESKD within 5 years after diagnostic kidney biopsy, yes or no) represented the output data. For these variables, we imputed missing values to the means for continuous characteristics and the mode for categorical characteristics. Because of missing information on serum triglycerides, total cholesterol, and low-density lipoprotein in some cases, 243 patients without blood lipid data were excluded to avoid inaccuracy due to missing value filling ( Figure 1 ).

To confirm that the T _pre can be used in the ESKD prediction model at multiple levels, we also constructed a lifetime ESKD prediction model based on XGBoost. The process and approach were the same as building the 5-year ESKD prediction model. The primary outcome was time-to-event ESKD. The survival time for the kidney without ESKD event was calculated from the kidney biopsy to the last follow-up.

The XGBoost was allowed to generate boosting trees at most 110 times, and the maximum depth of each tree was constrained to 5. To avoid overfitting, we further set the L2 regularization term on weights as 1 and stop training if the performance did not improve by more than 15 rounds. At last, the optimal prediction model parameters and architectures were selected by the five-fold cross-validation.

The patients of dataset 2 without Oxford MEST-C scores combined with the corresponding T _pre were used as an additional external validation set to evaluate the performance of the ESKD prediction model using T _pre.

The sociodemographic and clinical variables were calculated and expressed as the mean ± standard deviation for variables with approximately symmetrical distributions and as median (interquartile range 25th–75th percentile) for variables with skewed distribution. All categorical variables are expressed as frequencies and percentages. Univariate analyses based on the Cox proportional hazards model (36) were conducted to evaluate the association between the baseline clinical characteristics and ESKD event. Clinical characteristics associated with ESKD event in univariate analysis (P < 0.05) or if they were clinically relevant were used as input features of the 5-year ESKD prediction model.

For predicting 5-year ESKD status (yes or no) and T-score (0 or 1), the performance of the models was assessed by calculating the accuracy, sensitivity, specificity, and area under the receiver operating characteristic (ROC) curve (AUC). For predicting lifetime ESKD risk, we quantify the performance of the model by concordance statistic (C-statistic), which is a general concept of the area under the curve (AUC) for time-to-event survival data (37). The C-statistic compares the rank of predicting probability and the rank of the survival time in the real world. The calibration ability of the models was assessed by the Hosmer–Lemeshow test and calibration scatter plot, in which P-value >0.05 indicated no very significant difference between the predicted probability predicted by the model and the true outcome frequencies during a certain time period. SPSS version 26.0 software and R 3.6.3 were used for the statistical analysis. All P-values were two-tailed, and P <0.05 was considered statistically significant.

The clinical characteristics of 690 patients with IgAN in dataset 1 are shown in Table 1 . The mean age of these patients was 32.38 ± 11.32 years at the time of renal biopsy. The male-to-female ratio was 1.2:1. The mean arterial pressure was 94.44 ± 14.02 mm Hg. The median value of eGFR was 84.66 (range, 63.32–107.50) ml/min per 1.73 m², and daily proteinuria was 1.38 (range, 0.66–2.89) g/day.

For the 1,168 follow-up patients with Oxford MEST-C scores in dataset 2, the mean age was 35.10 ± 11.73 years at the time of renal biopsy. The male-to-female ratio was 1:1. The mean arterial pressure was 93.59 ± 11.42 mm Hg. The eGFR was 85.91 (range, 60.94–107.23) ml/min per 1.73 m², and daily proteinuria was 1.27 (range, 0.66–2.45) g/day ( Table 1 ). For the variables used to train the pathology T-score prediction (T _pre) model, there were no statistically significant differences in clinical parameters between dataset 1 and dataset 2 except for age (32.38 ± 11.32 vs. 35.10 ± 11.73, P = 1.00 × 10⁻⁶), serum IgA level (3.13 ± 1.21 vs. 3.29 ± 1.20, P = 0.01), and serum uric acid level (347.10 ± 114.95 vs. 367.63 ± 101.86, P = 1.52 × 10⁻⁴). Among these, 158 patients (13.53%) had reached the event of ESKD during the median 67.5-month follow-up. The unadjusted hazard ratios (HRs) between the different variables and ESKD are reported in Table 2 . The risk of ESKD significantly increased for every 10.0 mm Hg increase in the MAP [HR: 1.34, 95% confidence interval (CI): 1.18–1.53, P = 1.10 × 10⁻⁵] and increased for every 1.0 g/day in the daily proteinuria (HR: 1.10, 95% CI: 1.05–1.15, P = 1.60 × 10⁻⁵). For each ml/min per 1.73 m² decrease in eGFR, the risk of ESKD increased by 4% (HR: 0.96, 95% CI: 0.96–0.97, P = 1.24 × 10⁻²⁷). For each mg/dl increase in uric acid, the risk of ESKD increased by 38% (HR: 1.38, 95% CI: 1.29–1.49, P = 1.47 × 10⁻¹⁹). Moreover, there was the strongest association between the risk of ESKD and the presence of tubulointerstitial lesions (HR: 3.34, 95% CI: 2.73–4.07, P = 1.72 × 10⁻³²).

Feature reductions were conducted using the AUCRF algorithm, which was used to select the optimal random forest model with the least number of predictive variables to predict the presence or absence of T lesions. Clinical variables with a probability of selection higher than 0.7 were selected in repeated cross-validation of the optimal random forest model (optimal AUC = 0.82). Finally, the features selected by AUCRF for the T prediction model included age, systolic arterial pressure, diastolic arterial pressure, proteinuria, eGFR, serum IgA, and uric acid ( Figure 2 ). The 690 IgAN patients with Oxford MEST-C scores in dataset 1 as the training set were taken to develop a pathology T-score prediction model. The 1,168 IgAN patients with Oxford MEST-C scores in dataset 2 as the testing set were used only for reporting the performance of the model and were not used for development or fine-tuning.

If a predictive model has an AUC of higher than 0.75, it will be considered to have a good discriminating ability. The pathology T prediction model achieved a discrimination of 0.82 (95% CI: 0.80–0.85) [area under the receiver operating characteristic (ROC) curve (AUC)] in the testing set ( Figure 3A ). The ROC curve had 0.74 sensitivity and 0.77 specificity, which indicated that it had better clinical utility.

The unadjusted Cox regression analysis suggested that sex, systolic arterial pressure, diastolic arterial pressure, proteinuria, eGFR, uric acid, triglycerides, and tubular atrophy/interstitial fibrosis (T) were risk factors for developing ESKD ( Table 2 ). A study supported elevated serum IgA as a causal factor in IgA nephropathy through Mendelian randomization (38). Some studies have suggested the association between the poor prognosis of renal disease and dyslipidemia. Higher triglycerides and cholesterol levels have been proven to be independent risk factors for the progression of kidney disease (39). Hence, clinical variables (age, sex, systolic arterial pressure, diastolic arterial pressure, proteinuria, eGFR, serum IgA, uric acid, triglycerides, total cholesterol, low-density lipoprotein, history of previous use of RAS inhibitors and immunosuppressants) and the pathology T lesions (T _bio, T0 was assigned 0, T1 and T2 were assigned 1) were used as the input variables of the 5-year ESKD prediction model.

To make the predictive model achieve a good performance, the 1,168 follow-up IgAN patients with Oxford MEST-C scores in dataset 2 were randomly divided into training and testing sets at a ratio of 8:2. The training set included 936 patients and the testing set included 232 patients. The training set was used to perform five-fold cross-validation to select the optimal prediction model. The testing set was used to assess the performance.

The performance value of the 5-year ESKD prediction model using only the above clinical variables as input variables (base model) was 0.86 (95% CI: 0.75–0.97) in the test set ( Figure 3B ). To test whether the T _bio could improve the predictive performance of the 5-year ESKD prediction model, we added T _bio to the base model. An increase in AUC [from 0.86 (95% CI: 0.75–0.97) to 0.92 (95% CI: 0.85–0.98); P = 0.03] showed a better discriminating ability, which indicated that the T was important for judging the prognosis of patients with IgAN ( Figure 3B ). To test whether T _pre had a similar effect on judging the prognosis of IgAN patients, after training the 5-year ESKD prediction model with the training set, we replaced the T _bio in the testing set with the corresponding T _pre to see the discrimination effect. The AUC was 0.90 (95% CI: 0.82–0.99) in the testing set ( Figure 3B ). The performance of the base model plus T _pre did not differ from that of the base model plus T _bio [AUC for the base model plus T _pre 0.90 (95% CI: 0.82–0.99) vs. AUC for the base model plus T _bio 0.92 (95% CI: 0.85–0.98), P = 0.52, Table 3 ], which showed that the value of the T _pre in predicting the ESKD of patients was comparable to that of T _bio. The calibration of the three prediction models is shown in Figures 4A–C . The P-values for the Hosmer–Lemeshow test of the base model, the base model plus T _bio, and the base model plus T _pre were 0.42, 0.79, and 0.92, respectively, which indicated that these models had a good calibration. These results suggested the importance of T in predicting ESKD, and T _pre can be used to assist clinicians in assessing the prognosis of patients without pathology reports.

Table 4 shows the performance of the 5-year ESKD prediction model based on different machine learning algorithms in the testing set using T _pre. All models have good prediction performance, which indicated that T _pre could be used in ESKD predictive models built on different algorithms.

For the lifetime ESKD prediction model based on XGBoost using only clinical variables (base model), the C-statistic was 0.82 (95% CI: 0.80–0.84) in the testing set. The discriminating ability of the base model plus T _pre was also comparable to the base model plus T _bio [C-statistic: 0.85 (95% CI: 0.83–0.86) vs. 0.85 (95% CI: 0.83–0.86), P = 0.11] in the testing set.

The 355 patients without MEST-C scores in dataset 2 were included as the external validation population for evaluating the performance of the 5-year ESKD prediction model. Because patients did not have MEST-C scores, the T _pre predicted by the pathology T-score prediction model was used in the 5-year ESKD prediction model. The AUC of the 5-year ESKD prediction model using T _pre based on XGBoost was 0.93 (95% CI: 0.87–0.99). We listed the AUC of the applied other machine learning algorithms in Table 5 .

In the lifetime ESKD prediction model using T _pre, the C-statistic was 0.92 (95% CI: 0.90–0.94). We have shown here that both models have a good performance in the external validation set, indicating the reliability of T _pre for assisting in evaluating the prognosis of IgAN.