OSCR

Plasma proteomic signatures of early retinal neurodegeneration in diabetes: A multi-cohort study.

A correction to this paper has been published: the notice, 42461828, from Europe PMC.

Code ↔ Paper

1 match between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 1 match
  1. [1] § Methods › Statistical analysis ↔ R code for Pro-DRN.R, lines 693–707 · score 0.52 · age sex, Aspelund, Dagliati, Hippisley, ISDR, Tarasewicz

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,144 lines · 44 KB · no license · 1 match

  2. #———————————————————————————————————————————————————————————————————————————————————————————————————————————####
  3. # 线性回归（RNFL—蛋白）------------------------------------------------------------------------------------------------
  4. ------------------------------------------------------------------------------------------------------------
  5. rm(list = ls())
  6. gc()
  8. library(haven)
  9. library(readr)
  10. library(tidyr)
  11. library(dplyr)
  12. library(readxl)
  13. library(openxlsx)
  14. library(survival)
  15. library(fdrtool)
  16. library(progress)
  17. library(openxlsx)
  19. data_GDES <- read_dta("XXX.dta")
  20. data_pro <- read.xlsx("XXX.xlsx")
  21. # 存储要拟合的蛋白变量名
  22. proteins <- setdiff(names(data_pro), "id")
  23. proname <- read_tsv("XXX.tsv")
  24. proname <- proname %>%
  25. select(-coding) %>% # 删除 coding 列
  26. separate(meaning, into = c("abbreviation", "name"), sep = ";") %>% # 分离 meaning 列
  27. mutate(pro = tolower(abbreviation)) # 增加 pro 列，将 abbreviation 转为小写
  29. panel <- read_excel("XXX.xlsx")
  30. colnames(panel) <- panel[1,]
  31. panel <- panel[-1,-c(5,6)]
  32. colnames(panel) <- c("uniprot","name","pro","panel")
  33. panel <- panel %>%
  34. mutate(
  35. pro = tolower(pro), # 改小写
  36. pro = gsub("-", "_", pro) # "-"改"_"
  37. )
  38. panel <- panel[,c("pro","panel")]
  41. variables <- c("y0_prnfl_aver")
  43. # 年龄、性别、dm_duration、HbA1c、SBP、吸烟 ####
  44. # 创建一个列表，用于存储每个变量的回归结果
  45. lm_results_list <- list()
  46. pb <- progress_bar$new(
  47. format = " Progress [:bar] :percent Elapsed: :elapsed ETA: :eta",
  48. total = length(variables) * length(proteins),
  49. clear = FALSE,
  50. width = 60
  51. )
  53. # 循环进行亚场的线性回归
  54. for (var in variables) {
  55. # 创建一个空的数据框来存储回归结果
  56. lm_result <- data.frame(variable = character(),
  57. coef = numeric(),
  58. tvalue = numeric(),
  59. std = numeric(),
  60. pvalue = numeric(),
  61. stringsAsFactors = FALSE)
  63. # 循环拟合每个模型并存储结果
  64. for (var_pro in proteins) {
  65. formula <- paste0(var,"~", var_pro, "+ age + as.factor(sex) + dmdura_year + hba1c + as.factor(smoking) + sbp")
  66. lm_model <- lm(formula, data = data_GDES)
  68. coef <- as.numeric(summary(lm_model)$coefficients[, 1][2])
  69. std <- as.numeric(summary(lm_model)$coefficients[, 2][2])
  70. tvalue <- as.numeric(summary(lm_model)$coefficients[, 3][2])
  71. pvalue <- as.numeric(summary(lm_model)$coefficients[, 4][2])
  73. lm_result <- rbind(lm_result, c(var_pro, coef, tvalue, std, pvalue))
  74. # ✅ 更新进度条
  75. pb$tick()
  76. }
  78. colnames(lm_result) <- c("pro", "coef", "tvalue", "std", "pvalue")
  79. lm_result$coef <- as.numeric(lm_result$coef)
  80. lm_result$std <- as.numeric(lm_result$std)
  81. lm_result$tvalue <- as.numeric(lm_result$tvalue)
  82. lm_result$pvalue <- as.numeric(lm_result$pvalue)
  83. lm_result$low_limit <- as.numeric(lm_result$coef) - 1.96 * as.numeric(lm_result$std)
  84. lm_result$up_limit <- as.numeric(lm_result$coef) + 1.96 * as.numeric(lm_result$std)
  85. lm_result <- lm_result[c("pro", "coef", "low_limit", "up_limit", "std", "tvalue", "pvalue")]
  87. # 整理结果数据框
  88. lm_result$pvalue_sig <- ifelse(lm_result$pvalue <= 0.05, "significant", "-")
  90. # BH校正
  91. lm_result$BH <- p.adjust(lm_result$pvalue, "BH")
  92. lm_result$BH_sig <- ifelse(lm_result$BH <= 0.05, "significant", "-")
  93. fdr <- fdrtool(lm_result$pvalue, statistic="pvalue", plot=FALSE) # 关闭自动绘图
  94. lm_result$qval <- fdr$qval
  95. lm_result$lfdr <- fdr$lfdr
  97. # Bonferroni 校正与显著性标注
  98. lm_result$bonf <- p.adjust(lm_result$pvalue, method = "bonferroni")
  99. lm_result$bonf_sig <- ifelse(lm_result$bonf <= 0.05, "significant", "-")
  101. lm_result$BY <- p.adjust(lm_result$pvalue, method = "BY")
  102. lm_result$BY_sig <- ifelse(lm_result$BY <= 0.05, "significant", "-")
  104. lm_result$holm <- p.adjust(lm_result$pvalue, method = "holm")
  105. lm_result$holm_sig <- ifelse(lm_result$holm <= 0.05, "significant", "-")
  107. lm_result$hochberg <- p.adjust(lm_result$pvalue, method = "hochberg")
  108. lm_result$hochberg_sig <- ifelse(lm_result$hochberg <= 0.05, "significant", "-")
  111. # 注释蛋白
  112. lm_result <- left_join(lm_result, proname, by = "pro")
  113. lm_result <- left_join(lm_result, panel, by = "pro")
  116. #将每一个lm_result保存到相应的lm_results_list中
  117. lm_results_list[[var]] <- lm_result
  118. }
  120. # 输出结果
  121. lm_results_list
  123. # 写入到Excel
  124. library(openxlsx)
  125. wb <- createWorkbook()
  127. for (var in variables) {
  128. # 创建工作表
  129. addWorksheet(wb, var)
  130. # 写入数据到工作表
  131. writeData(wb, var, lm_results_list[[var]])
  132. }
  135. # Save the Excel file
  136. saveWorkbook(wb, "XXX.xlsx", overwrite = TRUE)
  140. #———————————————————————————————————————————————————————————————————————————————————————————————————————————####
  141. # 线性回归（蛋白-prnfl下降速率）------------------------------------------------------------------------------------------------
  142. ------------------------------------------------------------------------------------------------------------
  143. rm(list = ls())
  144. gc()
  146. library(haven)
  147. library(readr)
  148. library(tidyr)
  149. library(dplyr)
  150. library(readxl)
  151. library(openxlsx)
  152. library(survival)
  153. library(fdrtool)
  154. library(progress)
  155. library(openxlsx)
  157. data_GDES <- read_dta("XXX.dta")
  159. # 存储要拟合的蛋白变量名
  161. proname <- read_tsv("XXX.tsv")
  162. proname <- proname %>%
  163. select(-coding) %>% # 删除 coding 列
  164. separate(meaning, into = c("abbreviation", "name"), sep = ";") %>% # 分离 meaning 列
  165. mutate(pro = tolower(abbreviation)) # 增加 pro 列，将 abbreviation 转为小写
  167. panel <- read_excel("E:/UKB数据/olink数据/3072蛋白列表.xlsx")
  168. colnames(panel) <- panel[1,]
  169. panel <- panel[-1,-c(5,6)]
  170. colnames(panel) <- c("uniprot","name","pro","panel")
  171. panel <- panel %>%
  172. mutate(
  173. pro = tolower(pro), # 改小写
  174. pro = gsub("-", "_", pro) # "-"改"_"
  175. )
  176. panel <- panel[,c("pro","panel")]
  177. data_pro <- read_excel("XXX.xlsx")
  178. proteins <- data_pro$pro
  180. variables <- c("decline_rate")
  182. # 年龄、性别、dm_duration、HbA1c、SBP、吸烟 ####
  183. # 创建一个列表，用于存储每个变量的回归结果
  184. lm_results_list <- list()
  185. pb <- progress_bar$new(
  186. format = " Progress [:bar] :percent Elapsed: :elapsed ETA: :eta",
  187. total = length(variables) * length(proteins),
  188. clear = FALSE,
  189. width = 60
  190. )
  192. # 循环进行亚场的线性回归
  193. for (var in variables) {
  194. # 创建一个空的数据框来存储回归结果
  195. lm_result <- data.frame(variable = character(),
  196. coef = numeric(),
  197. tvalue = numeric(),
  198. std = numeric(),
  199. pvalue = numeric(),
  200. stringsAsFactors = FALSE)
  202. # 循环拟合每个模型并存储结果
  203. for (var_pro in proteins) {
  204. formula <- paste0(var,"~", var_pro, "+ age + as.factor(sex) + dmdura_year + hba1c + as.factor(smoking) + sbp")
  205. lm_model <- lm(formula, data = data_GDES)
  207. coef <- as.numeric(summary(lm_model)$coefficients[, 1][2])
  208. std <- as.numeric(summary(lm_model)$coefficients[, 2][2])
  209. tvalue <- as.numeric(summary(lm_model)$coefficients[, 3][2])
  210. pvalue <- as.numeric(summary(lm_model)$coefficients[, 4][2])
  212. lm_result <- rbind(lm_result, c(var_pro, coef, tvalue, std, pvalue))
  213. # ✅ 更新进度条
  214. pb$tick()
  215. }
  217. colnames(lm_result) <- c("pro", "coef", "tvalue", "std", "pvalue")
  218. lm_result$coef <- as.numeric(lm_result$coef)
  219. lm_result$std <- as.numeric(lm_result$std)
  220. lm_result$tvalue <- as.numeric(lm_result$tvalue)
  221. lm_result$pvalue <- as.numeric(lm_result$pvalue)
  222. lm_result$low_limit <- as.numeric(lm_result$coef) - 1.96 * as.numeric(lm_result$std)
  223. lm_result$up_limit <- as.numeric(lm_result$coef) + 1.96 * as.numeric(lm_result$std)
  224. lm_result <- lm_result[c("pro", "coef", "low_limit", "up_limit", "std", "tvalue", "pvalue")]
  226. # 整理结果数据框
  227. lm_result$pvalue_sig <- ifelse(lm_result$pvalue <= 0.05, "significant", "-")
  229. # BH校正
  230. lm_result$BH <- p.adjust(lm_result$pvalue, "BH")
  231. lm_result$BH_sig <- ifelse(lm_result$BH <= 0.05, "significant", "-")
  232. fdr <- fdrtool(lm_result$pvalue, statistic="pvalue", plot=FALSE) # 关闭自动绘图
  233. lm_result$qval <- fdr$qval
  234. lm_result$lfdr <- fdr$lfdr
  236. # 注释蛋白
  237. lm_result <- left_join(lm_result, proname, by = "pro")
  238. lm_result <- left_join(lm_result, panel, by = "pro")
  241. #将每一个lm_result保存到相应的lm_results_list中
  242. lm_results_list[[var]] <- lm_result
  243. }
  245. # 输出结果
  246. lm_results_list
  248. # 写入到Excel
  249. library(openxlsx)
  250. wb <- createWorkbook()
  252. for (var in variables) {
  253. # 创建工作表
  254. addWorksheet(wb, var)
  255. # 写入数据到工作表
  256. writeData(wb, var, lm_results_list[[var]])
  257. }
  260. # Save the Excel file
  261. saveWorkbook(wb, "XXXX.xlsx", overwrite = TRUE)
  265. # 通路分析 #####
  266. rm(list = ls())
  267. gc()
  269. ####加载包
  270. library(zlibbioc)
  271. library(GO.db)
  272. library(org.Hs.eg.db)
  273. library(clusterProfiler)
  274. library(tidyverse)
  275. library(ggplot2)
  276. library(forcats)
  277. library(readxl)
  278. library(writexl)
  280. ####加载数据
  281. data <- read_excel("XXXX.xlsx")
  283. ##基因ID转换（genesymbol → ENTREZID)
  284. IDs <- bitr(data$abbreviation,
  285. fromType = 'SYMBOL',
  286. toType = c('ENTREZID'),
  287. OrgDb = 'org.Hs.eg.db') #不同的物种对应不同的数据库，详见https://bioconductor.org/packages/release/BiocViews.html#___OrgDb
  289. ##基于ENTREZID进行GO富集分析
  290. enrichGO_result <- enrichGO(IDs$ENTREZID,
  291. OrgDb="org.Hs.eg.db",
  292. qvalueCutoff=0.05,
  293. pvalueCutoff=0.05,
  294. ont="all")
  296. enrichGO_result <- setReadable(enrichGO_result, OrgDb="org.Hs.eg.db", keyType = 'ENTREZID')
  298. enrichGO_result_data <- data.frame(enrichGO_result)
  300. writexl::write_xlsx(enrichGO_result_data,
  301. path = "XXXX.xlsx")
  303. ##基于ENTREZID进行KEGG富集分析
  304. enrichKEGG_result <- enrichKEGG(gene=IDs$ENTREZID,
  305. organism ="hsa",#物种
  306. keyType ="kegg",
  307. pvalueCutoff = 0.05,
  308. qvalueCutoff = 1
  309. )
  311. enrichKEGG_result <- setReadable(enrichKEGG_result,
  312. OrgDb="org.Hs.eg.db",
  313. keyType ="ENTREZID"
  314. )
  316. enrichKEGG_result_data <- data.frame(enrichKEGG_result)
  319. writexl::write_xlsx(enrichKEGG_result_data,
  320. path = "XXXX.xlsx")
  324. #———————————————————————————————————————————————————————————————————————————————————————————————————————————####
  325. # 4. 预测分析 #####
  326. rm(list = ls())
  327. gc()
  329. library(haven)
  330. library(readr)
  331. library(tidyr)
  332. library(dplyr)
  333. library(readxl)
  334. library(openxlsx)
  335. library(survival)
  336. library(fdrtool)
  337. library(progress)
  338. library(openxlsx)
  340. #整理数据
  341. data_GDES <- read_dta("XXXX.dta")
  342. data_train <- read_dta("XXXX.dta")
  343. data_test <- read_dta("XXXX.dta")
  345. # 计算 Q1 阈值
  346. q1_cut <- quantile(data_GDES$decline_rate, 0.25, na.rm = TRUE)
  348. # 根据 Q1 定义 drn
  349. data_GDES <- data_GDES %>%
  350. mutate(
  351. drn = ifelse(decline_rate <= q1_cut, 1, 0)
  352. )
  354. # 蛋白组模型(xgboost) ####
  355. outcome_list <- c("drn")
  356. results_pro_all <- data.frame()
  357. original_seed_results <- list()
  358. seed_list <- createWorkbook()
  360. # 用全量人群 id 作为对齐基准
  361. stopifnot("id" %in% names(data_GDES))
  362. data_score <- data.frame(id = data_GDES$id)
  363. set.seed(2025)
  365. params <- list(
  366. booster = "gbtree",
  367. objective = "binary:logistic",
  368. eval_metric = "auc",
  369. eta = 0.05,
  370. max_depth = 7,
  371. min_child_weight = 1,
  372. gamma = 2,
  373. subsample = 0.9,
  374. colsample_bytree = 1,
  375. colsample_bylevel = 0.5,
  376. lambda = 1,
  377. alpha = 0.5,
  378. max_delta_step = 0,
  379. tree_method = "exact",
  380. scale_pos_weight = 1,
  381. seed = 2025
  382. )
  385. for (outcome in outcome_list) {
  387. # --- 清理NA（保持最小必要清理） ---
  388. stopifnot(all(c("id", outcome) %in% colnames(data_train)))
  389. stopifnot(all(c("id", outcome) %in% colnames(data_test)))
  390. data_train <- data_train[!is.na(data_train[[outcome]]), ]
  391. data_test <- data_test[ !is.na(data_test[[outcome]]), ]
  393. # --- 特征对齐：仅保留在 train/test 都存在的蛋白列 ---
  394. proteins <- intersect(proteins, intersect(colnames(data_train), colnames(data_test)))
  395. stopifnot(length(proteins) > 0)
  397. # --- 标签与特征矩阵 ---
  398. train_y <- as.numeric(data_train[[outcome]])
  399. test_y <- as.numeric(data_test[[outcome]])
  400. train_x <- as.matrix(data_train[, proteins, drop = FALSE])
  401. test_x <- as.matrix(data_test[, proteins, drop = FALSE])
  403. dtrain <- xgb.DMatrix(data = train_x, label = train_y)
  404. dtest <- xgb.DMatrix(data = test_x, label = test_y)
  406. # --- 用训练集做交叉验证，寻找最佳轮次 ---
  407. val.cv <- xgb.cv(
  408. params = params,
  409. data = dtrain,
  410. nfold = 10,
  411. nrounds = 1000,
  412. early_stopping_rounds = 100,
  413. metrics = "auc",
  414. maximize = TRUE,
  415. verbose = 0
  416. )
  417. best_ntree <- val.cv$best_iteration
  419. # --- 用最佳轮次在训练集上拟合最终模型 ---
  420. model <- xgboost(
  421. data = dtrain,
  422. params = params,
  423. nrounds = best_ntree,
  424. verbose = 0
  425. )
  427. # --- 训练/测试集预测与ROC ---
  428. pred_train <- predict(model, dtrain)
  429. pred_test <- predict(model, dtest)
  431. roc_pro_all_train <- pROC::roc(as.factor(train_y) ~ pred_train, ci = TRUE, quiet = TRUE)
  432. roc_pro_all <- pROC::roc(as.factor(test_y) ~ pred_test, ci = TRUE, quiet = TRUE)
  434. # --- cutoff & 混淆矩阵（基于测试集） ---
  435. cutoff <- cutoff::roc(pred_test, as.factor(test_y))$cutoff
  436. prediction <- ifelse(pred_test > cutoff, 1, 0)
  437. cm <- caret::confusionMatrix(
  438. data = factor(prediction, levels = c(0,1)),
  439. reference = factor(test_y, levels = c(0,1)),
  440. positive = "1", mode = "everything"
  441. )
  442. TN <- cm$table[1]; FP <- cm$table[2]; FN <- cm$table[3]; TP <- cm$table[4]
  443. Accuracy <- cm$overall[["Accuracy"]]
  444. Sensitivity <- cm$byClass[["Sensitivity"]]
  445. Specificity <- cm$byClass[["Specificity"]]
  446. PPV <- cm$byClass[["Precision"]]
  447. NPV <- cm$byClass[["Neg Pred Value"]]
  448. Recall <- cm$byClass[["Recall"]]
  449. F1 <- cm$byClass[["F1"]]
  450. DR <- TP/(FN+TP)
  451. FPR <- FP/(TN+FP)
  452. LR <- ifelse(FPR == 0, Inf, DR/FPR)
  454. # --- 结果表 ---
  455. test_results <- data.frame(
  456. outcome = outcome,
  457. seed_pro_all = NA_real_, # 占位
  458. train_roc_pro_all = as.numeric(roc_pro_all_train$auc),
  459. train_low_pro_all = roc_pro_all_train$ci[1],
  460. train_up_pro_all = roc_pro_all_train$ci[3],
  461. test_roc_pro_all = as.numeric(roc_pro_all$auc),
  462. test_low_pro_all = roc_pro_all$ci[1],
  463. test_up_pro_all = roc_pro_all$ci[3],
  464. cutoff = cutoff,
  465. Accuracy = Accuracy, Sensitivity = Sensitivity, Specificity = Specificity,
  466. PPV = PPV, NPV = NPV, Recall = Recall, F1 = F1,
  467. DR = DR, FPR = FPR, LR = LR,
  468. stringsAsFactors = FALSE
  469. )
  471. original_seed_results[[outcome]] <- test_results
  472. addWorksheet(seed_list, sheetName = outcome)
  473. writeData(seed_list, outcome, original_seed_results[[outcome]])
  475. # --- 保存逐个体 score 到全量框架（按 id 对齐；score 列名改为 score_drn） ---
  476. score_col <- "score_drn" # ← 你要求的命名
  477. dset_col <- "dataset_drn" # 数据集标记列，便于区分
  479. data_score_temp <- data_GDES %>%
  480. dplyr::select(id) %>%
  481. dplyr::left_join(data.frame(id = data_test$id, pred = pred_test), by = "id") %>%
  482. dplyr::left_join(data.frame(id = data_train$id, pred_train = pred_train), by = "id") %>%
  483. dplyr::mutate(dataset = dplyr::case_when(!is.na(pred) ~ "test",
  484. !is.na(pred_train) ~ "train",
  485. TRUE ~ NA_character_),
  486. pred_final = ifelse(!is.na(pred), pred, pred_train)) %>%
  487. dplyr::transmute(id, !!score_col := pred_final, !!dset_col := dataset)
  489. data_score <- dplyr::left_join(data_score, data_score_temp, by = "id")
  490. }
  492. saveWorkbook(seed_list, "XXXX.xlsx", overwrite = TRUE)
  493. write.csv(data_score, "XXXXX.csv", row.names = FALSE)
  495. # 保存环境
  496. save.image("XXXXX.RData")
  499. # ——4.2.2 变量重要性 ------------------------------------------------------------------------------------------------
  500. rm(list = ls())
  501. load("XXXXX.RData")
  503. library(openxlsx)
  504. library(DT)
  505. variable_importance <- createWorkbook()
  508. # ———4.2.2.1 xgboost包自带VIP ####
  509. importance <- xgb.importance(feature_names = colnames(train), model = xgb_model)
  510. #1 Gain：使用某个特征进行拆分时，获得的平均训练损失减少量
  511. #2 Cover：首先得到某个特征被用于在所有树中拆分数据的次数，然后要利用经过这些拆分点的训练数据数量赋予权重
  512. #3 frequency：某个特征被用于在所有树中拆分数据的次数
  513. addWorksheet(variable_importance, sheetName = "importance")
  514. writeData(variable_importance, "importance", as.data.frame(importance))
  515. # 图示
  516. xgb.plot.importance(importance, measure = "Gain")
  517. pdf(file="XXXXX.pdf",onefile = FALSE,
  518. width = 8,
  519. height =18)
  520. xgb.plot.importance(importance, measure = "Gain")
  521. dev.off()
  522. # ggplot可视化排名前10变量
  523. library(ggplot2)
  524. xgb.ggplot.importance(importance_matrix = importance, top_n = 10)
  527. # ———4.2.2.2 SHAP ####
  528. library(SHAPforxgboost)
  529. library(dplyr)
  530. # SHAP data (individual-unit)
  531. shap_data <- shap.prep(xgb_model, X_train = train_x)
  532. addWorksheet(variable_importance, sheetName = "SHAP data")
  533. writeData(variable_importance, "SHAP data", as.data.frame(shap_data))
  534. # SHAP summary
  535. shap_summary <- shap_data %>%
  536. group_by(variable) %>%
  537. summarise(AverageMeanValue = mean(mean_value, na.rm = TRUE))
  538. addWorksheet(variable_importance, sheetName = "SHAP summary")
  539. writeData(variable_importance, "SHAP summary", as.data.frame(shap_summary))
  543. # *4.4.3 总人群预测 —— 汇总预测分数 ####
  544. data_GDES <- read_dta("XXXXX.dta")
  545. data_score_temp <- read_csv("XXXXX.csv")
  547. data_GDES <- data_GDES %>%
  548. left_join(
  549. data_score_temp %>%
  550. select(id, score_drn_XGB, dataset_drn),
  551. by = "id"
  552. )
  555. data_train <- data_GDES %>% filter(dataset_drn == "train")
  556. data_test <- data_GDES %>% filter(dataset_drn == "test")
  558. write_dta(data_train, "XXXXX.dta")
  559. write_dta(data_test, "XXXXX.dta")
  560. write_dta(data_GDES, "XXXXX.dta")
  562. # 单独predictor比较 ####
  564. rm(list = ls())
  565. gc()
  567. library(haven)
  568. library(readr)
  569. library(tidyr)
  570. library(dplyr)
  571. library(readxl)
  572. library(openxlsx)
  573. library(survival)
  574. library(fdrtool)
  575. library(progress)
  576. library(pROC)
  578. # 整理数据
  579. data_train <- read_dta("XXXXX.dta")
  580. data_test <- read_dta("XXXXX.dta")
  581. data_GDES <- read_dta("XXXXX.dta")
  583. outcomes <- c("drn")
  584. auc_results <- data.frame(stringsAsFactors = FALSE)
  586. for (outcome in outcomes) {
  588. # -------- 基准：蛋白（XGB 分数） --------
  589. model_pro <- glm(as.formula(paste(outcome, "~ score_drn_XGB")),
  590. data = data_test, family = binomial)
  591. data_test[[paste0(outcome, "_riskscore_pro")]] <-
  592. predict(model_pro, newdata = data_test, type = "response")
  594. roc_result_pro <- pROC::roc(
  595. response = data_test[[outcome]],
  596. predictor = data_test[[paste0(outcome, "_riskscore_pro")]],
  597. ci = TRUE, ci.method = "delong", na.rm = TRUE, quiet = TRUE
  598. )
  599. auc_pro <- as.numeric(roc_result_pro$auc)
  600. ci_lower_pro <- roc_result_pro$ci[1]
  601. ci_upper_pro <- roc_result_pro$ci[3]
  603. fit_and_roc <- function(rhs) {
  604. m <- glm(as.formula(paste(outcome, "~", rhs)), data = data_test, family = binomial)
  605. sc <- predict(m, newdata = data_test, type = "response")
  606. r <- pROC::roc(response = data_test[[outcome]], predictor = sc,
  607. ci = TRUE, ci.method = "delong", na.rm = TRUE, quiet = TRUE)
  608. list(model = m, roc = r, auc = as.numeric(r$auc), ci = r$ci)
  609. }
  611. res_age <- fit_and_roc("age")
  612. res_sex <- fit_and_roc("as.factor(sex)")
  613. res_income <- fit_and_roc("as.factor(income)")
  614. res_smoking <- fit_and_roc("as.factor(smoking)")
  615. res_drinking <- fit_and_roc("as.factor(drinking)")
  616. res_education <- fit_and_roc("as.factor(education)")
  617. res_bmi <- fit_and_roc("as.factor(bmi_c)")
  618. res_statin <- fit_and_roc("as.factor(med_statin)")
  619. res_hbp <- fit_and_roc("as.factor(med_hbp)")
  622. P_age <- pROC::roc.test(roc_result_pro, res_age$roc, method = "delong")$p.value
  623. P_sex <- pROC::roc.test(roc_result_pro, res_sex$roc, method = "delong")$p.value
  624. P_income <- pROC::roc.test(roc_result_pro, res_income$roc, method = "delong")$p.value
  625. P_smoking <- pROC::roc.test(roc_result_pro, res_smoking$roc, method = "delong")$p.value
  626. P_drinking <- pROC::roc.test(roc_result_pro, res_drinking$roc, method = "delong")$p.value
  627. P_education <- pROC::roc.test(roc_result_pro, res_education$roc, method = "delong")$p.value
  628. P_bmi <- pROC::roc.test(roc_result_pro, res_bmi$roc, method = "delong")$p.value
  629. P_statin <- pROC::roc.test(roc_result_pro, res_statin$roc, method = "delong")$p.value
  630. P_hbp <- pROC::roc.test(roc_result_pro, res_hbp$roc, method = "delong")$p.value
  632. lab <- function(pv, auc_other) ifelse(pv >= 0.05, "Comparable",
  633. ifelse(auc_pro > auc_other, "Superior", "-"))
  635. auc_age <- res_age$auc; ci_lower_age <- res_age$ci[1]; ci_upper_age <- res_age$ci[3]
  636. auc_sex <- res_sex$auc; ci_lower_sex <- res_sex$ci[1]; ci_upper_sex <- res_sex$ci[3]
  637. auc_income <- res_income$auc; ci_lower_income <- res_income$ci[1]; ci_upper_income <- res_income$ci[3]
  638. auc_smoking <- res_smoking$auc; ci_lower_smoking <- res_smoking$ci[1]; ci_upper_smoking <- res_smoking$ci[3]
  639. auc_drinking <- res_drinking$auc; ci_lower_drinking <- res_drinking$ci[1]; ci_upper_drinking <- res_drinking$ci[3]
  640. auc_education <- res_education$auc; ci_lower_education <- res_education$ci[1]; ci_upper_education <- res_education$ci[3]
  641. auc_bmi <- res_bmi$auc; ci_lower_bmi <- res_bmi$ci[1]; ci_upper_bmi <- res_bmi$ci[3]
  642. auc_med_statin <- res_statin$auc; ci_lower_med_statin <- res_statin$ci[1]; ci_upper_med_statin <- res_statin$ci[3]
  643. auc_med_hbp <- res_hbp$auc; ci_lower_med_hbp <- res_hbp$ci[1]; ci_upper_med_hbp <- res_hbp$ci[3]
  645. compare_age <- lab(P_age, auc_age)
  646. compare_sex <- lab(P_sex, auc_sex)
  647. compare_income <- lab(P_income, auc_income)
  648. compare_smoking <- lab(P_smoking, auc_smoking)
  649. compare_drinking <- lab(P_drinking, auc_drinking)
  650. compare_education <- lab(P_education, auc_education)
  651. compare_bmi <- lab(P_bmi, auc_bmi)
  652. compare_med_statin<- lab(P_statin, auc_med_statin)
  653. compare_med_hbp <- lab(P_hbp, auc_med_hbp)
  655. auc_results <- rbind(auc_results, c(
  656. outcome,
  657. auc_pro, ci_lower_pro, ci_upper_pro,
  659. auc_age, ci_lower_age, ci_upper_age, P_age, compare_age,
  660. auc_sex, ci_lower_sex, ci_upper_sex, P_sex, compare_sex,
  661. auc_income, ci_lower_income, ci_upper_income, P_income, compare_income,
  662. auc_smoking, ci_lower_smoking, ci_upper_smoking, P_smoking, compare_smoking,
  663. auc_drinking, ci_lower_drinking, ci_upper_drinking, P_drinking, compare_drinking,
  664. auc_education, ci_lower_education, ci_upper_education, P_education, compare_education,
  665. auc_bmi, ci_lower_bmi, ci_upper_bmi, P_bmi, compare_bmi,
  666. auc_med_statin, ci_lower_med_statin, ci_upper_med_statin, P_statin, compare_med_statin,
  667. auc_med_hbp, ci_lower_med_hbp, ci_upper_med_hbp, P_hbp, compare_med_hbp
  668. ))
  670. colnames(auc_results) <- c(
  671. "outcome", "PRO", "pro_lower", "pro_upper",
  672. "AGE", "age_lower", "age_upper", "Pvalue_AGE", "PRO_vs._AGE",
  673. "SEX", "sex_lower", "sex_upper", "Pvalue_SEX", "PRO_vs._SEX",
  674. "INCOME", "income_lower", "income_upper", "Pvalue_INCOME", "PRO_vs._INCOME",
  675. "SMOKING", "smoking_lower", "smoking_upper", "Pvalue_SMOKING", "PRO_vs._SMOKING",
  676. "DRINKING", "drinking_lower", "drinking_upper", "Pvalue_DRINKING", "PRO_vs._DRINKING",
  677. "EDUCATION", "education_lower", "education_upper", "Pvalue_EDUCATION", "PRO_vs._EDUCATION",
  678. "BMI", "bmi_lower", "bmi_upper", "Pvalue_BMI", "PRO_vs._BMI",
  679. "ANTI_LIPID", "anti_lipid_lower", "anti_lipid_upper", "Pvalue_ANTI_LIPID", "PRO_vs._ANTI_LIPID",
  680. "ANTI_BP", "anti_bp_lower", "anti_bp_upper", "Pvalue_ANTI_BP", "PRO_vs._ANTI_BP"
  681. )
  682. }
  684. # 转数值
  685. cmp_cols <- grep("^PRO_vs\\._", colnames(auc_results), value = TRUE) # 比较标签列
  686. num_cols <- setdiff(colnames(auc_results), c("outcome", cmp_cols)) # 其余才转数值
  688. auc_results[cmp_cols] <- lapply(auc_results[cmp_cols], as.character) # 显式设为字符
  689. auc_results[num_cols] <- lapply(auc_results[num_cols], function(x) suppressWarnings(as.numeric(x)))
  691. write.xlsx(auc_results, "XXXXX.xlsx")
  693. # 模型和结局函数 ####
  694. create_models <- function() {
  695. models <- list("Age_Sex" = "age + as.factor(sex)",
  696. "Aspelund_model" = "age + as.factor(sex) + as.factor(sbp_c) + dmdura_year + hba1c",
  697. "Hippisley_model" = "age + as.factor(sex) + as.factor(bmi_c) + as.factor(sbp_c) + TC_divide_HDL + hba1c",
  698. "Dagliati_model" = "age + as.factor(sex) + dmdura_year + as.factor(bmi_c) + hba1c + as.factor(prev_hbp) + as.factor(smoking)",
  699. "ISDR_model" = "age + as.factor(sex) + as.factor(sbp_c) + dmdura_year + hba1c + chol",
  700. "JJR_model" = "age + as.factor(sbp_c) + dmdura_year + hba1c + Egfr",
  701. "Tarasewicz_model" = "age + as.factor(sbp_c) + as.factor(insulin) + hba1c + as.factor(bmi_c) + ldlc + crea + mau",
  702. "All_model" = "age + as.factor(sex) + as.factor(sbp_c) + dmdura_year + hba1c + as.factor(bmi_c) + TC_divide_HDL + chol + Egfr + as.factor(insulin) + ldlc + crea + mau + as.factor(prev_hbp) + as.factor(smoking)"
  704. )
  706. return(list(models = models))
  707. }
  709. models <- create_models()$models
  712. # 人群预测 ####
  713. library(pROC)
  714. library(openxlsx)
  715. library(dplyr)
  717. score_vars <- c(
  718. "score_drn_XGB","score_drn_CART","score_drn_KNN","score_drn_LGBM",
  719. "score_drn_LogReg","score_drn_NN","score_drn_RF","score_drn_SVM"
  720. )
  721. outcomes <- c("drn")
  722. auc_results_list <- list()
  724. wb <- createWorkbook()
  726. for (model_name in names(models)) {
  727. auc_results <- data.frame()
  729. # —— 把 outcome 循环改成 score 循环
  730. for (score_var in score_vars) {
  731. outcome <- outcomes[1]
  732. covariate <- models[[model_name]]
  734. # 传统模型
  735. formula1 <- as.formula(paste(outcome, "~", covariate, sep = ""))
  736. model1 <- glm(formula1, data = data_test, family = binomial)
  737. data_test[[paste(outcome, "_", model_name, "_1", sep = "")]] <- predict(model1, newdata = data_test)
  738. roc_result1 <- roc(
  739. as.formula(paste(outcome, "~", paste(outcome, "_", model_name, "_1", sep = ""))),
  740. data = data_test, print.thres = TRUE, print.auc = TRUE, ci = TRUE, plot = FALSE, smooth = FALSE
  741. )
  742. auc_value1 <- roc_result1$auc
  743. ci_lower1 <- roc_result1$ci[1]
  744. ci_upper1 <- roc_result1$ci[3]
  746. # 加入蛋白（这里用当前 score_var）
  747. formula2 <- as.formula(paste(outcome, "~", covariate, "+", score_var))
  748. model2 <- glm(formula2, data = data_test, family = binomial)
  749. data_test[[paste(outcome, "_", model_name, "_2", sep = "")]] <- predict(model2, newdata = data_test)
  750. roc_result2 <- roc(
  751. as.formula(paste(outcome, "~", paste(outcome, "_", model_name, "_2", sep = ""))),
  752. data = data_test, print.thres = TRUE, print.auc = TRUE, ci = TRUE, plot = FALSE, smooth = FALSE
  753. )
  754. auc_value2 <- roc_result2$auc
  755. ci_lower2 <- roc_result2$ci[1]
  756. ci_upper2 <- roc_result2$ci[3]
  757. test <- roc.test(roc_result1, roc_result2) # DeLong's test
  758. p <- test[["p.value"]]
  759. auc_diff <- auc_value2 - auc_value1
  760. auc_improve_proportion <- auc_diff / auc_value1
  761. sig <- ifelse(p >= 0.05, "",
  762. ifelse(p < 0.001, "***",
  763. ifelse(p >= 0.001 & p < 0.01, "**", "*")))
  765. # 单独蛋白（当前 score_var）
  766. roc_result3 <- roc(
  767. as.formula(paste(outcome, "~", score_var)),
  768. data = data_test, print.thres = TRUE, print.auc = TRUE, ci = TRUE, plot = FALSE, smooth = FALSE
  769. )
  770. auc_value3 <- roc_result3$auc
  771. ci_lower3 <- roc_result3$ci[1]
  772. ci_upper3 <- roc_result3$ci[3]
  773. test_met <- roc.test(roc_result1, roc_result3) # DeLong's test
  774. p_met <- test_met[["p.value"]]
  775. compare_met <- ifelse(p_met >= 0.05, "Comparable", ifelse(auc_value3 > auc_value1, "Superior", "-"))
  777. # 结果入表
  778. auc_results <- rbind(
  779. auc_results,
  780. c(score_var, outcome,
  781. auc_value1, ci_lower1, ci_upper1,
  782. auc_value3, ci_lower3, ci_upper3, p_met, compare_met,
  783. auc_value2, ci_lower2, ci_upper2, auc_diff, auc_improve_proportion, p, sig)
  784. )
  785. colnames(auc_results) <- c("score", "outcome",
  786. "CONV", "conv_lower", "conv_upper",
  787. "PRO", "pro_lower", "pro_upper", "Pvalue_PRO_vs_CONV", "PRO_vs._CONV",
  788. "COMBINED", "combined_lower", "combined_upper", "Diff", "Diff_proportion",
  789. "Pvalue_COMBINED_vs_CONV", "Significance")
  790. }
  792. # 显著标记放在最后一列（保留 score 列）
  793. auc_results <- auc_results %>%
  794. select("score", which(names(auc_results) == "outcome"):which(names(auc_results) == "Pvalue_COMBINED_vs_CONV"), "Significance")
  796. # as.numeric
  797. auc_results$CONV <- as.numeric(auc_results$CONV)
  798. auc_results$conv_lower <- as.numeric(auc_results$conv_lower)
  799. auc_results$conv_upper <- as.numeric(auc_results$conv_upper)
  800. auc_results$PRO <- as.numeric(auc_results$PRO)
  801. auc_results$pro_lower <- as.numeric(auc_results$pro_lower)
  802. auc_results$pro_upper <- as.numeric(auc_results$pro_upper)
  803. auc_results$Pvalue_PRO_vs_CONV <- as.numeric(auc_results$Pvalue_PRO_vs_CONV)
  804. auc_results$COMBINED <- as.numeric(auc_results$COMBINED)
  805. auc_results$combined_lower <- as.numeric(auc_results$combined_lower)
  806. auc_results$combined_upper <- as.numeric(auc_results$combined_upper)
  807. auc_results$Diff <- as.numeric(auc_results$Diff)
  808. auc_results$Diff_proportion <- as.numeric(auc_results$Diff_proportion)
  809. auc_results$Pvalue_COMBINED_vs_CONV <- as.numeric(auc_results$Pvalue_COMBINED_vs_CONV)
  811. # 汇总与写表
  812. auc_results_list[[model_name]] <- auc_results
  813. addWorksheet(wb, sheetName = model_name)
  814. writeData(wb, model_name, auc_results_list[[model_name]])
  815. }
  816. auc_results_list
  817. # ===== 按 score 汇总成“图2样子”的表（每个 score 一个 sheet）=====
  818. # 先把所有模型的结果合并，并带上 model 列
  819. all_res <- dplyr::bind_rows(lapply(names(auc_results_list), function(m) {
  820. df <- auc_results_list[[m]]
  821. df$model <- m
  822. df
  823. }))
  825. # 列顺序
  826. col_order <- c("model","score","outcome",
  827. "CONV","conv_lower","conv_upper",
  828. "PRO","pro_lower","pro_upper","Pvalue_PRO_vs_CONV","PRO_vs._CONV",
  829. "COMBINED","combined_lower","combined_upper",
  830. "Diff","Diff_proportion","Pvalue_COMBINED_vs_CONV","Significance")
  832. # 逐个 score 输出一个 sheet
  833. for (sv in score_vars) {
  834. tag <- sub("^score_drn_", "", sv) # e.g. "XGB"
  835. sheet_name <- substr(paste0("Summary_", tag), 1, 31) # Excel 名称≤31字符
  837. df_out <- all_res %>%
  838. dplyr::filter(score == sv) %>%
  839. # 模型顺序排
  840. dplyr::mutate(model = factor(model, levels = names(models))) %>%
  841. dplyr::arrange(model)
  843. # 只保留需要的列
  844. keep_cols <- intersect(col_order, names(df_out))
  845. df_out <- df_out[, keep_cols, drop = FALSE]
  847. # 写入新工作表
  848. if (sheet_name %in% openxlsx::sheets(wb)) {
  849. openxlsx::removeWorksheet(wb, sheet_name)
  850. }
  851. openxlsx::addWorksheet(wb, sheetName = sheet_name)
  852. openxlsx::writeData(wb, sheet_name, df_out)
  853. }
  854. # 保存 Excel
  855. saveWorkbook(wb, "XXXXX.xlsx", overwrite = TRUE)
  857. # 保存预测label
  858. data_label <- data_test[, 2845:ncol(data_test)]
  859. data_id <- data_test[, 1]
  860. data_predict_label <- cbind(data_id, data_label)
  861. write.xlsx(data_predict_label, "XXXXX.xlsx")
  866. # =========================== RF训练 ===========================
  867. rm(list = ls()); gc()
  869. library(haven)
  870. library(dplyr)
  871. library(caret)
  872. library(ranger)
  873. library(pROC)
  874. library(cutoff)
  875. library(openxlsx)
  876. library(readxl)
  877. library(readr)
  878. library(tidyr)
  880. data_train <- read_dta("XXXXX.dta")
  881. data_test <- read_dta("XXXXX.dta")
  882. data_GDES <- read_dta("XXXXX.dta")
  884. data_pro2 <- read_excel("XXXXX.xlsx")
  885. proteins <- data_pro2$pro
  887. # 传统模型集合
  888. create_models <- function() {
  889. models <- list(
  890. "Age_Sex" = "age + as.factor(sex)",
  891. "Aspelund_model"= "age + as.factor(sex) + as.factor(sbp_c) + dmdura_year + hba1c",
  892. "Hippisley_model"= "age + as.factor(sex) + as.factor(bmi_c) + as.factor(sbp_c) + TC_divide_HDL + hba1c",
  893. "Dagliati_model"= "age + as.factor(sex) + dmdura_year + as.factor(bmi_c) + hba1c + as.factor(prev_hbp) + as.factor(smoking)",
  894. "ISDR_model" = "age + as.factor(sex) + as.factor(sbp_c) + dmdura_year + hba1c + chol",
  895. "JJR_model" = "age + sbp + dmdura_year + hba1c + Egfr",
  896. "Tarasewicz_model" = "age + as.factor(sbp_c) + as.factor(insulin) + hba1c + as.factor(bmi_c) + ldlc + crea + mau",
  897. "All_model" = "age + as.factor(sex) + as.factor(sbp_c) + dmdura_year + hba1c + as.factor(bmi_c) + TC_divide_HDL + chol + Egfr + as.factor(insulin) + ldlc + crea + mau + as.factor(prev_hbp) + as.factor(smoking)"
  898. )
  899. models
  900. }
  901. models <- create_models()
  904. outcome <- "drn"
  905. data_score <- data.frame(id = data_GDES$id)
  906. seed_level_perf <- data.frame()
  907. wb_seed_summary <- createWorkbook()
  908. wb_model_compare <- createWorkbook()
  911. # 仅删除结局缺失（ranger 不接受 y 中有 NA）；不做其它预处理/填补
  912. stopifnot(all(c("id","drn") %in% colnames(data_train)))
  913. stopifnot(all(c("id","drn") %in% colnames(data_test)))
  915. data_train <- data_train %>% filter(!is.na(drn))
  916. data_test <- data_test %>% filter(!is.na(drn))
  918. # 与训练/测试集取交集，确保特征存在
  919. proteins <- intersect(proteins, intersect(colnames(data_train), colnames(data_test)))
  920. stopifnot(length(proteins) > 0)
  922. # （沿用你原来做法）仅在训练集上去近零方差，锁定特征清单
  923. nzv <- nearZeroVar(data_train[, proteins, drop = FALSE])
  924. proteins_use <- if (length(nzv) > 0) proteins[-nzv] else proteins
  925. stopifnot(length(proteins_use) > 0)
  927. # 因变量转为 0/1 因子（仅供 ranger 用；不改变 data_train/data_test 中原始 drn）
  928. train_y_fac <- factor(data_train$drn, levels = c(0, 1), labels = c("0", "1"))
  929. test_y_fac <- factor(data_test$drn, levels = c(0, 1), labels = c("0", "1"))
  930. # 训练集必须包含两个类别
  931. stopifnot(length(unique(train_y_fac)) == 2)
  933. # 特征矩阵（不做填补/缩放等）
  934. train_x <- data_train[, proteins_use, drop = FALSE]
  935. test_x <- data_test[, proteins_use, drop = FALSE]
  938. set.seed(2025)
  940. y_fac <- train_y_fac
  941. X <- train_x
  943. # ——1) 构造分层 k 折（类别平衡）——
  944. K <- 5
  945. folds <- caret::createFolds(y_fac, k = K, list = TRUE, returnTrain = FALSE)
  947. # ——2) 参数网格（可按需精简/扩展）——
  948. p <- ncol(X)
  949. grid <- expand.grid(
  950. num.trees = c(500), # 可加 800, 1000
  951. mtry = unique(pmax(1, round(c(sqrt(p), 0.05*p, 0.1*p)))),
  952. min.node.size = c(5, 10, 20),
  953. sample.fraction = c(0.6, 0.8) # <1.0 可缓解过拟合
  954. )
  956. cv_res <- grid
  957. cv_res$mean_auc <- NA_real_
  958. cv_res$sd_auc <- NA_real_
  960. # ——3) 逐参数组合跑 k 折，记录验证 AUC——
  961. for (gi in seq_len(nrow(grid))) {
  962. g <- grid[gi, ]
  963. aucs <- numeric(K)
  965. for (k in seq_len(K)) {
  966. idx_te <- folds[[k]]
  967. idx_tr <- setdiff(seq_len(nrow(X)), idx_te)
  969. rf_k <- ranger::ranger(
  970. formula = as.formula("y ~ ."),
  971. data = data.frame(y = y_fac[idx_tr], X[idx_tr, , drop = FALSE]),
  972. num.trees = g$num.trees,
  973. mtry = g$mtry,
  974. min.node.size = g$min.node.size,
  975. sample.fraction = g$sample.fraction,
  976. replace = TRUE,
  977. probability = TRUE,
  978. importance = "impurity",
  979. oob.error = TRUE
  980. )
  982. pred_val <- predict(rf_k, data = data.frame(X[idx_te, , drop = FALSE]))$predictions[, "1"]
  983. aucs[k] <- as.numeric(pROC::roc(response = y_fac[idx_te], predictor = pred_val, quiet = TRUE)$auc)
  984. }
  986. cv_res$mean_auc[gi] <- mean(aucs, na.rm = TRUE)
  987. cv_res$sd_auc[gi] <- stats::sd(aucs, na.rm = TRUE)
  988. }
  990. # ——4) 选择最佳参数（按平均 AUC 最大；若并列取 sd 最小）——
  991. best_idx <- order(-cv_res$mean_auc, cv_res$sd_auc)[1]
  992. best_par <- cv_res[best_idx, , drop = FALSE]
  993. print(best_par)
  995. # ——5) 用最佳参数在「全部训练集」上重训最终模型，并拿到 OOB 预测/AUC——
  996. rf_model <- ranger::ranger(
  997. formula = as.formula(paste(outcome, "~ .")),
  998. data = data.frame(drn = y_fac, X),
  999. num.trees = best_par$num.trees,
  1000. mtry = best_par$mtry,
  1001. min.node.size = best_par$min.node.size,
  1002. sample.fraction = best_par$sample.fraction,
  1003. replace = TRUE,
  1004. probability = TRUE,
  1005. importance = "impurity"
  1006. )
  1008. # ——7) 训练/测试集预测与ROC ——
  1009. pred_train_prob <- predict(rf_model, data = data.frame(train_x))$predictions[, "1"]
  1010. pred_test_prob <- predict(rf_model, data = data.frame(test_x))$predictions[, "1"]
  1012. roc_train <- pROC::roc(response = train_y_fac, predictor = pred_train_prob, ci = TRUE, quiet = TRUE)
  1013. roc_test <- pROC::roc(response = test_y_fac, predictor = pred_test_prob, ci = TRUE, quiet = TRUE)
  1015. # cutoff & 混淆矩阵与派生指标
  1016. if (length(unique(test_y_fac)) == 2) {
  1017. cutoff_val <- cutoff::roc(pred_test_prob, test_y_fac)$cutoff
  1018. pred_label <- ifelse(pred_test_prob > cutoff_val, 1, 0)
  1020. cm <- caret::confusionMatrix(
  1021. data = factor(pred_label, levels = c(0,1), labels = c("0","1")),
  1022. reference = test_y_fac,
  1023. positive = "1", mode = "everything"
  1024. )
  1025. TN <- cm$table[1]; FP <- cm$table[2]; FN <- cm$table[3]; TP <- cm$table[4]
  1026. Accuracy <- cm$overall[["Accuracy"]]
  1027. Sensitivity <- cm$byClass[["Sensitivity"]]
  1028. Specificity <- cm$byClass[["Specificity"]]
  1029. PPV <- cm$byClass[["Precision"]]
  1030. NPV <- cm$byClass[["Neg Pred Value"]]
  1031. Recall <- cm$byClass[["Recall"]]
  1032. F1 <- cm$byClass[["F1"]]
  1033. DR <- TP/(FN+TP) # detection rate
  1034. FPR <- FP/(TN+FP) # false positive rate
  1035. LR <- ifelse(FPR == 0, Inf, DR/FPR)
  1037. } else {
  1038. # 单一类别时占位（只保留AUC；其余置NA）
  1039. cutoff_val <- NA
  1040. Accuracy <- Sensitivity <- Specificity <- PPV <- NPV <- Recall <- F1 <- DR <- FPR <- LR <- NA
  1041. }
  1043. seed_level_perf <- rbind(seed_level_perf, data.frame(
  1044. method = "RandomForest",
  1045. train_auc = as.numeric(roc_train$auc),
  1046. train_ci_low = roc_train$ci[1],
  1047. train_ci_up = roc_train$ci[3],
  1048. test_auc = as.numeric(roc_test$auc),
  1049. test_ci_low = roc_test$ci[1],
  1050. test_ci_up = roc_test$ci[3],
  1051. cutoff = cutoff_val,
  1052. Accuracy, Sensitivity, Specificity, PPV, NPV, Recall, F1, DR, FPR, LR
  1053. ))
  1055. # =========================== 保存逐个体score到全量框架（附dataset标记）
  1056. score_col <- "score_drn_RF"
  1057. dset_col <- "dataset_drn_RF"
  1059. data_score_temp <- data_GDES %>%
  1060. select(id) %>%
  1061. left_join(data.frame(id = data_test$id, pred = pred_test_prob), by = "id") %>%
  1062. left_join(data.frame(id = data_train$id, pred_train = pred_train_prob), by = "id") %>%
  1063. mutate(dataset = case_when(!is.na(pred) ~ "test",
  1064. !is.na(pred_train) ~ "train",
  1065. TRUE ~ NA_character_),
  1066. pred_final = ifelse(!is.na(pred), pred, pred_train)) %>%
  1067. transmute(id, !!score_col := pred_final, !!dset_col := dataset)
  1069. data_score <- left_join(data_score, data_score_temp, by = "id")
  1071. # 供下面模型比较使用
  1072. data_test[[score_col]] <- pred_test_prob
  1074. # =========================== 测试集：传统模型 vs RF score vs 组合模型
  1075. auc_rows_all_models <- data.frame()
  1077. for (model_name in names(models)) {
  1078. covars <- models[[model_name]]
  1080. # 传统模型
  1081. f1 <- as.formula(paste(outcome, "~", covars))
  1082. m1 <- glm(f1, data = data_test, family = binomial)
  1083. lp1 <- predict(m1, newdata = data_test) # logit 线性预测器即可；AUC对单调变换不敏感
  1084. roc1 <- pROC::roc(response = data_test[[outcome]],
  1085. predictor = as.numeric(lp1),
  1086. ci = TRUE, quiet = TRUE)
  1088. # 单独 RF score
  1089. roc3 <- pROC::roc(response = data_test[[outcome]],
  1090. predictor = data_test[[score_col]],
  1091. ci = TRUE, quiet = TRUE)
  1093. # 组合模型（传统 + RF分数）
  1094. f2 <- as.formula(paste(outcome, "~", covars, "+", score_col))
  1095. m2 <- glm(f2, data = data_test, family = binomial)
  1096. lp2 <- predict(m2, newdata = data_test)
  1097. roc2 <- pROC::roc(response = data_test[[outcome]],
  1098. predictor = as.numeric(lp2),
  1099. ci = TRUE, quiet = TRUE)
  1101. # DeLong检验（同样显式指向 pROC）
  1102. p_met <- as.numeric(pROC::roc.test(roc1, roc3)$p.value) # RF vs CONV
  1103. p_comb <- as.numeric(pROC::roc.test(roc1, roc2)$p.value) # COMBINED vs CONV
  1105. auc_diff <- as.numeric(roc2$auc) - as.numeric(roc1$auc)
  1106. auc_impr <- auc_diff / as.numeric(roc1$auc)
  1107. sig <- ifelse(p_comb >= 0.05, "",
  1108. ifelse(p_comb < 0.001, "***",
  1109. ifelse(p_comb < 0.01, "**", "*")))
  1110. compare_met <- ifelse(p_met >= 0.05, "Comparable",
  1111. ifelse(as.numeric(roc3$auc) > as.numeric(roc1$auc), "Superior", "-"))
  1113. auc_rows_all_models <- rbind(auc_rows_all_models, data.frame(
  1114. method = "RandomForest",
  1115. outcome = outcome,
  1116. model = model_name,
  1117. CONV = as.numeric(roc1$auc), conv_lower = roc1$ci[1], conv_upper = roc1$ci[3],
  1118. RF = as.numeric(roc3$auc), rf_lower = roc3$ci[1], rf_upper = roc3$ci[3],
  1119. Pvalue_RF_vs_CONV = p_met, RF_vs_CONV = compare_met,
  1120. COMBINED = as.numeric(roc2$auc), combined_lower = roc2$ci[1], combined_upper = roc2$ci[3],
  1121. Diff = auc_diff, Diff_proportion = auc_impr,
  1122. Pvalue_COMBINED_vs_CONV = p_comb, Significance = sig,
  1123. stringsAsFactors = FALSE
  1124. ))
  1125. }
  1127. # 1) RF汇总
  1128. addWorksheet(wb_seed_summary, "RF_Summary")
  1129. writeData(wb_seed_summary, "RF_Summary", seed_level_perf)
  1130. saveWorkbook(wb_seed_summary, "XXXX.xlsx", overwrite = TRUE)
  1132. # 2) 按模型分sheet的AUC对比结果
  1133. for (model_name in names(models)) {
  1134. sub <- auc_rows_all_models %>% filter(model == model_name)
  1135. addWorksheet(wb_model_compare, sheetName = model_name)
  1136. writeData(wb_model_compare, model_name, sub)
  1137. }
  1138. saveWorkbook(wb_model_compare, "XXXX.xlsx", overwrite = TRUE)
  1140. # 3) 个体层score（包含score_drn_RF与dataset_drn_RF）
  1141. write.csv(data_score, "XXXX.csv", row.names = FALSE)
  1143. # # 4) 保存模型与对象
  1144. save.image("XXXX.RData")

R code for Pro-DRN.R at commit 459bdfe, no license · at the source

Overview

Authors: Huangdong Li1,2, Ziyu Zhu1,2, Shaopeng Yang1,2, Weijing Cheng1,2, Shaoying Tan3,4,5, Zhuoyao Xin6, Lei Zhang7,8,9, Zhuoting Zhu10, Shida Chen1,2, Wenyong Huang1,2, Wei Wang1,2,11
ORCID iDs: Ziyu Zhu, Shaopeng Yang, Shaoying Tan, Lei Zhang, Shida Chen, Wei Wang
  1. State Key Laboratory of Ophthalmology, Zhongshan Ophthalmic Center, Sun Yat-Sen University, Guangdong Provincial Key Laboratory of Ophthalmology and Visual Science, Guangdong Provincial Clinical Research Center for Ocular Diseases, Guangzhou, China
  2. Guangdong Basic Research Center of Excellence for Major Blinding Eye Diseases Prevention and Treatment, Guangzhou, China
  3. School of Optometry, The Hong Kong Polytechnic University, Kowloon, Hong Kong SAR, China
  4. Research Centre for SHARP Vision (RCSV), The Hong Kong Polytechnic University, Hong Kong, Hong Kong SAR, China
  5. Centre for Eye and Vision Research (CEVR), Hong Kong, Hong Kong SAR, China
  6. Department of Biomedical Engineering, Columbia University, New York, New York, United States of America
  7. Clinical Medical Research Center, Children’s Hospital of Nanjing Medical University, Nanjing, Jiangsu Province, China
  8. Artificial Intelligence and Modelling in Epidemiology Program, Melbourne Sexual Health Centre, Alfred Health, Melbourne, Australia
  9. Central Clinical School, Faculty of Medicine, Nursing and Health Sciences, Monash University, Melbourne, Australia
  10. Centre for Eye Research Australia, Royal Victorian Eye and Ear Hospital, Melbourne, Australia
  11. Hainan Eye Hospital and Key Laboratory of Ophthalmology, Zhongshan Ophthalmic Center, Sun Yat-sen University, Haikou, China
Institutions: State Key Laboratory of Ophthalmology; Zhongshan Ophthalmic Center, Sun Yat-sen University; The Royal Victorian Eye & Ear Hospital (Australia); Centre for Eye Research Australia (Australia); Hong Kong Polytechnic University (Hong Kong SAR China); Columbia University (United States); Melbourne Sexual Health Centre (Australia); Second Affiliated Hospital of Nanjing Medical University (China); Alfred Health (Australia); Monash University (Australia); Nanjing Medical University (China); Hainan Eye Hospital (China)
Journal: PLoS medicine, volume 23, issue 6, article e1004868
Dates: received 30 November 2025; accepted 20 April 2026; published online 2 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pmed.1004868 · PMID 42228639 · PMCID PMC13229346 · OpenAlex W7163180447
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), other condition (population), clinical / translational (subfield)
Methods: Statistics, Machine learning, Physiology & signal measures
MeSH: Diabetes Mellitus, Type 2*, Diabetic Retinopathy*, Proteomics*, Retinal Degeneration*, Biomarkers, Cohort Studies, Female, Humans, Male, Middle Aged, Prospective Studies, Tomography, Optical Coherence (* major topic)
Journal subjects: Medicine and Health Sciences, Endocrinology, Endocrine Disorders, Diabetes Mellitus, Medical Conditions, Metabolic Disorders, Epidemiology, Medical Risk Factors, Biology and Life Sciences, Biochemistry, Proteomics, Research and Analysis Methods, Database and Informatics Methods, Biological Databases, Proteomic Databases, Proteins, Plasma Proteins, Mathematical and Statistical Techniques, Statistical Methods, Forecasting, Physical Sciences, Mathematics, Statistics, Diagnostic Medicine, Diagnostic Radiology, Tomography, Imaging Techniques, Radiology and Imaging, Biomarkers
Topic: Retinal Diseases and Treatments (Ophthalmology, Medicine), according to OpenAlex
Funding: Guangdong Basic Research Center of Excellence (N/A); Hainan Province Clinical Medical Center (N/A); Guangzhou Municipal Science and Technology Program key projects (2025A04J7150); National Natural Science Foundation of China (82371086, 82301253, 82571271); Natural Science Foundation of Guangdong Province (2026A1515010675); Guangdong Basic and Applied Basic Research Foundation (2022A151511); Projects of Research Center for Sharp Vision at The Hong Kong Polytechnic University (P0057931); Health and Medical Research Fund-Research Fellowship Scheme (07210207); Lumitin Vision to Brightness Research Funding for the Young and middle-aged Ophthalmologists (BCF-KH-YK-20230803-03)
Citations: cited by 2 papers (Europe PMC); 85 references in the paper
Notices: A correction to this paper has been published (42461828, from Europe PMC)

Abstract

Background: Retinal neurodegeneration is an early and independent feature of diabetic retinal disease and has been proposed as a window into the systemic neural consequences of diabetes, yet accessible molecular biomarkers and individualized prediction tools remain scarce. We aimed to identify circulating plasma protein signatures of diabetic retinal neurodegeneration (DRN) and to translate them into a clinically usable risk prediction system.

Methods and findings: In this multi-cohort prospective observational study, we integrated high-throughput plasma proteomics with longitudinal optical coherence tomography (OCT) in two independent populations. The discovery cohort comprised 1,492 participants had baseline plasma proteomics and OCT, and 1,218 were followed with repeated OCT over 6 years in Guangzhou Diabetic Eye Study (GDES). DRN was quantified by the annualized OCT-derived retinal nerve fiber layer thinning rate. In multivariable analyses adjusted for age, sex, smoking, systolic blood pressure, HbA1c, and diabetes duration, we identified 71 plasma proteins associated with development and progression of DRN. These proteins mapped onto pathways governing inflammatory immune recruitment, extracellular matrix remodeling, and microvascular homeostasis, providing a plausible biological basis for DRN. We developed a proteomics-based DRN model (Pro-DRN) using eight machine learning (ML) algorithms, including XGBoost and LightGBM. In the independent test set, Pro-DRN achieved a C-index of 0.860, rising to 0.908 when integrated with clinical variables. Compared with six conventional models, Pro-DRN improved discrimination (ΔC-index 0.137 to 0.159; all P < 0.001), reclassification (IDI 0.212 to 0.245; NRI 0.226 to 0.452; all P < 0.05). In the Hippisley model, the C-index increased from 0.739 (95% CI [0.670, 0.808]) to 0.898 (95% CI [0.858, 0.937]), with IDI 0.245 (95% CI [0.177, 0.318]), NRI 0.452 (95% CI [0.222, 0.673]) (both P < 0.001), and higher net benefit. The proteins most consistently driving model performance included ACTA2, COL6A3, and HSPG2. For clinical translation, we deployed the locked model as an interactive, web-based risk-assessment tool to support early DRN screening and longitudinal monitoring. Cross-ethnic external validation in UK Biobank (n = 502; recruited 2006–2010) reproduced core protein signals and consistent effect directions, confirming robustness across populations. Principal methodological limitation lies in single time point proteomic assessment.

Conclusion: In this multi-cohort study, we present a proteomics- and ML–based precision prediction system for DRN. Pro-DRN substantially enhanced early risk stratification beyond conventional clinical factors and may support targeted screening and timely neuroprotective interventions, advancing molecularly guided strategies for diabetic eye disease prevention.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 1 match between paragraphs and lines of code.

ZOC-skl/Pro-DRN

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 459bdfe6dacaa07319ae981066700c7caf894863, 24 February 2026
Languages: R (1)
Size: 2 files, 1 script
Software Heritage: not archived
Found in: “Data Availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: caret (1 file), clusterProfiler (1 file), ggplot2 (1 file), pROC (1 file), survival (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
2 files
At the source: github.com/ZOC-skl/Pro-DRN

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1 script, each with its path and the digest of its content;
  • 1 match between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

Data from the Guangzhou Diabetic Eye Study (GDES) are subject to Chinese data privacy regulations and are not publicly available. The data are managed by the Preventive Ophthalmology Data Management Unit, Zhongshan Ophthalmic Center, Sun Yat-sen University. Requests for access should be directed to the data manager via email at . Requests will be reviewed within 90 days according to Center policy, and, if approved, an inter-institutional data use agreement specifying use for non-commercial academic research will be required. UK Biobank data are available through application via the official platform (http://www.ukbiobank.ac.uk) under a material transfer agreement (application number: 105658). All software used in this study is publicly available. The code used in this study can be accessed at https://github.com/ZOC-skl/Pro-DRN.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 11 authors, 12 MeSH terms, 9 funders, 85 references, 1 integrity notice.

Cite

This paper

Li, H., Zhu, Z., Yang, S., Cheng, W., Tan, S., Xin, Z., Zhang, L., Zhu, Z., Chen, S., Huang, W., & Wang, W. (2026). Plasma proteomic signatures of early retinal neurodegeneration in diabetes: A multi-cohort study. PLoS medicine, 23(6), e1004868. https://doi.org/10.1371/journal.pmed.1004868

BibTeX

@article{li2026plasma,
author = {Li, Huangdong and Zhu, Ziyu and Yang, Shaopeng and Cheng, Weijing and Tan, Shaoying and Xin, Zhuoyao and Zhang, Lei and Zhu, Zhuoting and Chen, Shida and Huang, Wenyong and Wang, Wei},
title = {{Plasma proteomic signatures of early retinal neurodegeneration in diabetes: A multi-cohort study}},
journal = {PLoS medicine},
year = {2026},
month = jun,
volume = {23},
number = {6},
pages = {e1004868},
publisher = {PLOS},
issn = {1549-1277},
doi = {10.1371/journal.pmed.1004868},
url = {https://doi.org/10.1371/journal.pmed.1004868},
pmid = {42228639},
pmcid = {PMC13229346}
}

RIS

TY - JOUR
AU - Li, Huangdong
AU - Zhu, Ziyu
AU - Yang, Shaopeng
AU - Cheng, Weijing
AU - Tan, Shaoying
AU - Xin, Zhuoyao
AU - Zhang, Lei
AU - Zhu, Zhuoting
AU - Chen, Shida
AU - Huang, Wenyong
AU - Wang, Wei
TI - Plasma proteomic signatures of early retinal neurodegeneration in diabetes: A multi-cohort study
T2 - PLoS medicine
J2 - PLoS Med
PY - 2026
DA - 2026/06/02
VL - 23
IS - 6
SP - e1004868
SN - 1549-1277
PB - PLOS
DO - 10.1371/journal.pmed.1004868
UR - https://doi.org/10.1371/journal.pmed.1004868
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pmed.1004868",
"type": "article-journal",
"title": "Plasma proteomic signatures of early retinal neurodegeneration in diabetes: A multi-cohort study",
"container-title": "PLoS medicine",
"author": [
{
"family": "Li",
"given": "Huangdong"
},
{
"family": "Zhu",
"given": "Ziyu"
},
{
"family": "Yang",
"given": "Shaopeng"
},
{
"family": "Cheng",
"given": "Weijing"
},
{
"family": "Tan",
"given": "Shaoying"
},
{
"family": "Xin",
"given": "Zhuoyao"
},
{
"family": "Zhang",
"given": "Lei"
},
{
"family": "Zhu",
"given": "Zhuoting"
},
{
"family": "Chen",
"given": "Shida"
},
{
"family": "Huang",
"given": "Wenyong"
},
{
"family": "Wang",
"given": "Wei"
}
],
"container-title-short": "PLoS Med",
"volume": "23",
"issue": "6",
"page": "e1004868",
"DOI": "10.1371/journal.pmed.1004868",
"PMID": "42228639",
"PMCID": "PMC13229346",
"ISSN": "1549-1277",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pmed.1004868",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
2
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1182/bloodadvances.2025019294 [code]
Plasma proteome signatures are predictive of mortality in sickle cell disease.
Journal: Blood advances
In common: genetics / omics, 8 references
[2] doi:10.3390/ijms27156925 [code]
XGBoost-SHAP Interpretable Modeling Identifies and Validates an Eight-Gene Biomarker for Hepatic Encephalopathy Risk Prediction in Cirrhosis.
Journal: International journal of molecular sciences
In common: pROC, survival, caret, 3 other tools, genetics / omics, other condition
[3] doi:10.1073/pnas.2516601123 [code]
Unveiling the glymphatic system's role in brain aging: A comprehensive biomarker and modifiable intervention target.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: pROC, survival, ggplot2, 1 other tool, clinical / translational, 2 references
[4] doi:10.1093/neuonc/noag128 [code]
Spatially-resolved single-cell imaging of melanoma brain metastases identifies localized immune patterns predictive of immune checkpoint blockade response.
Journal: Neuro-oncology
In common: pROC, survival, caret, 2 other tools, clinical / translational, genetics / omics, other condition
[5] doi:10.21037/jtd-2026-0997 [code]
Machine learning models based on XGBoost algorithm to predict prognosis of lung cancer brain metastases.
Journal: Journal of thoracic disease
In common: pROC, survival, caret, 2 other tools, clinical / translational, other condition
[6] doi:10.1038/s41467-026-77170-3 [code]
DNA methylation profiling identifies long-range epigenetic silencing of clustered protocadherins as a key determinant of meningioma progression.
Journal: Nature communications
In common: pROC, survival, caret, 2 other tools, genetics / omics, other condition
[7] doi:10.1016/j.isci.2026.115657 [code]
Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.
Journal: iScience
In common: pROC, survival, clusterProfiler, 2 other tools, clinical / translational, other condition
[8] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: pROC, survival, clusterProfiler, 2 other tools, genetics / omics, other condition
[9] doi:10.1080/20002297.2026.2705667 [code]
Oral microbiota dysbiosis related to the cortical thinning and cognitive impairment in cerebral small vessel disease.
Journal: Journal of oral microbiology
In common: pROC, survival, caret, 2 other tools
[10] doi:10.7717/peerj.21426 [code]
Integrated transcriptomic identification and validation reveal key autophagy-associated biomarkers in sleep deprivation.
Journal: PeerJ
In common: pROC, caret, clusterProfiler, 2 other tools, genetics / omics

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.