OSCR

XGBoost-SHAP Interpretable Modeling Identifies and Validates an Eight-Gene Biomarker for Hepatic Encephalopathy Risk Prediction in Cirrhosis.

Code ↔ Paper

18 matches 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 18 matches
  1. [1] § 4. Materials and Methods › 4.6. Construction and Validation of XGBoost–SHAP Models ↔ r.07_xgboost.R, lines 29–136 · score 0.89 · colsample_bytree, XGBoost regression, cross validation, squarederror, subsampling, depth
  2. [2] § 4. Materials and Methods › 4.2. Differential Expression Analysis ↔ r.01_DEG.R, lines 15–100 · score 0.86 · contrasts.fit, eBayes, lmFit, linear model, limma, volcano
  3. [3] § 2. Results › 2.1. Workflow of the Study ↔ r.02_Mfuzz.R, lines 15–56 · score 0.79 · chronic liver failure, eCLD, decompensated cirrhosis, ACLF, acute, CC
  4. [4] § 2. Results › 2.1. Workflow of the Study ↔ r.11_PCA.R, lines 17–63 · score 0.78 · eCLD, chronic liver failure, decompensated cirrhosis, ACLF, CC, healthy
  5. [5] § 4. Materials and Methods › 4.5. Machine Learning-Based Marker Gene Selection ↔ r.06_RF.R, lines 17–75 · score 0.74 · randomForest, IncMSE, IncNodePurity, regression, modeling, gene
  6. [6] § 4. Materials and Methods › 4.8. Gene Set Enrichment Analysis (GSEA) ↔ r.10_GSEA.R, lines 68–119 · score 0.74 · gradient permutation, minSize, maxSize, Spearman, NES, thresholds
  7. [7] § 2. Results › 2.5. Multi-Algorithm Cross-Validated Screening for HE-Specific Marker Genes ↔ r.06_RF.R, lines 17–75 · score 0.70 · IncMSE, IncNodePurity, random forest, genes
  8. [8] § 4. Materials and Methods › 4.3. Functional Enrichment Analysis ↔ r.03_GOKEGG.R, lines 39–92 · score 0.69 · enrichGO, enrichKEGG, bitr, BH, enrichment, gene
  9. [9] § 4. Materials and Methods › 4.6. Construction and Validation of XGBoost–SHAP Models ↔ r.12_ROC.R, lines 150–213 · score 0.64 · fusion weights, Risk scores, w1, w2, AUC, models
  10. [10] § 4. Materials and Methods › 4.8. Gene Set Enrichment Analysis (GSEA) ↔ R/fgsea.R, lines 316–376 · score 0.63 · minSize, maxSize, fgsea, BH, NES, permutation
  11. [11] § 4. Materials and Methods › 4.7. Structural Equation Modeling and Mediation Effect Analysis ↔ r.09_SEM.R, lines 305–351 · score 0.60 · TUBA1C, mediation, mediated, bootstrap, indirect, SEM
  12. [12] § 2. Results › 2.4. Identification of Prognosis-Associated DEGs ↔ r.01_DEG.R, lines 102–148 · score 0.57 · good prognosis, poor prognosis, limma, DEGs, cirrhotic, GSE15654
  13. [13] § 4. Materials and Methods › 4.7. Structural Equation Modeling and Mediation Effect Analysis ↔ r.09_SEM.R, lines 70–93 · score 0.56 · ComBat, batch correction, SEM, Modeling
  14. [14] § 4. Materials and Methods › 4.6. Construction and Validation of XGBoost–SHAP Models ↔ r.08_model_validation.R, lines 141–200 · score 0.56 · KM survival, model validation, probability, median, risk
  15. [15] § 4. Materials and Methods › 4.1. Data Acquisition and Cohort Construction ↔ r.11_PCA.R, lines 17–63 · score 0.56 · eCLD, ACLF, DC, CC, healthy, clustering
  16. [16] § 2. Results › 2.5. Multi-Algorithm Cross-Validated Screening for HE-Specific Marker Genes ↔ r.04_LASSO.R, lines 122–164 · score 0.55 · partial likelihood deviance, cross validation, optimal, LASSO
  17. [17] § 2. Results › 2.6. Dual-Model Construction Based on XGBoost–SHAP and Multi-Level Clinical Validation ↔ packaging.R, lines 1–47 · score 0.54 · SHapley, exPlanations, Additive, XGBoost, model
  18. [18] § 2. Results › 2.7. Structural Equation Modeling-Based Analyses of Characteristic Gene Regulatory Networks ↔ r.09_SEM.R, lines 305–351 · score 0.51 · TUBA1C, mediated, indirect, SEM, mediation, LRRC32

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 · 484 lines · 17 KB · MIT · 3 matches

  1. rm(list = ls()); gc()
  2. ORIGINAL_DIR <- ""
  3. output <- file.path(ORIGINAL_DIR, "09_SEM")
  5. if (!dir.exists(output)) {
  6. dir.create(output, recursive = TRUE)
  7. }
  8. setwd(ORIGINAL_DIR)
  10. library(tidyverse)
  11. library(mediation)
  12. library(limma)
  13. library(sva)
  14. library(corrplot)
  15. library(pheatmap)
  16. library(lavaan)
  17. library(semPlot)
  18. library(dplyr)
  19. library(ggplot2)
  20. library(scales)
  21. library(grid)
  22. library(ggrepel)
  24. # Helper function for normality testing
  25. test_normality <- function(exp_data, by_group = FALSE, group_info = NULL) {
  26. if (is.data.frame(exp_data)) exp_data <- as.matrix(exp_data)
  28. if (!by_group) {
  29. pvals <- apply(exp_data, 1, function(x) {
  30. x <- as.numeric(x)
  31. if (all(is.na(x)) || sum(!is.na(x)) < 3) return(NA_real_)
  32. res <- tryCatch(shapiro.test(x), error = function(e) NULL)
  33. if (is.null(res)) NA_real_ else as.numeric(res$p.value)
  34. })
  36. df <- data.frame(gene = rownames(exp_data), p.value = as.numeric(pvals), stringsAsFactors = FALSE)
  37. df$adj.p <- p.adjust(df$p.value, method = "BH")
  38. return(df[order(df$adj.p), ])
  39. } else {
  40. if (is.null(group_info)) stop("group_info required for by_group = TRUE")
  41. groups <- levels(factor(group_info))
  42. res_list <- list()
  44. for (g in groups) {
  45. cols <- which(group_info == g)
  46. if (length(cols) < 3) {
  47. tmp <- data.frame(gene = rownames(exp_data), group = g, p.value = NA_real_, stringsAsFactors = FALSE)
  48. } else {
  49. pvals <- apply(exp_data[, cols, drop = FALSE], 1, function(x) {
  50. x <- as.numeric(x)
  51. if (all(is.na(x)) || sum(!is.na(x)) < 3) return(NA_real_)
  52. res <- tryCatch(shapiro.test(x), error = function(e) NULL)
  53. if (is.null(res)) NA_real_ else as.numeric(res$p.value)
  54. })
  55. tmp <- data.frame(gene = rownames(exp_data), group = g, p.value = as.numeric(pvals), stringsAsFactors = FALSE)
  56. tmp$adj.p <- p.adjust(tmp$p.value, method = "BH")
  57. }
  58. res_list[[g]] <- tmp
  59. }
  61. df_all <- do.call(rbind, lapply(res_list, function(x) {
  62. if (!"adj.p" %in% colnames(x)) x$adj.p <- NA_real_
  63. x
  64. }))
  65. rownames(df_all) <- NULL
  66. return(df_all[order(df_all$group, df_all$adj.p), ])
  67. }
  68. }
  70. # Helper function for batch correction
  71. perform_batch_correction <- function(exp_data, batch_info, group_info = NULL) {
  72. sample_order <- colnames(exp_data)
  73. batch_factor <- batch_info[match(sample_order, names(batch_info))]
  74. unique_batches <- unique(batch_factor)
  76. cat("Unique batches:", unique_batches, "\n")
  77. cat("Number of unique batches:", length(unique_batches), "\n")
  79. if (length(unique_batches) > 1 && !any(is.na(unique_batches))) {
  80. if (!is.null(group_info)) {
  81. group_factor <- group_info[match(sample_order, names(group_info))]
  82. mod <- model.matrix(~ group_factor)
  83. exp_corrected <- ComBat(dat = exp_data, batch = batch_factor, mod = mod)
  84. } else {
  85. exp_corrected <- ComBat(dat = exp_data, batch = batch_factor)
  86. }
  87. } else {
  88. warning("Batch correction not performed - insufficient batches")
  89. exp_corrected <- exp_data
  90. }
  92. return(exp_corrected)
  93. }
  95. # Helper function for PCA plot
  96. create_pca_plot <- function(exp_data, group_info, project_info = NULL, title = "PCA Plot") {
  97. pca <- prcomp(t(exp_data), scale. = TRUE)
  98. pca_data <- data.frame(
  99. PC1 = pca$x[, 1],
  100. PC2 = pca$x[, 2],
  101. Sample = colnames(exp_data),
  102. Group = group_info
  103. )
  105. if (!is.null(project_info)) {
  106. pca_data$Project <- project_info
  107. }
  109. var_explained <- pca$sdev^2 / sum(pca$sdev^2)
  110. pc1_var <- round(var_explained[1] * 100, 2)
  111. pc2_var <- round(var_explained[2] * 100, 2)
  113. p <- ggplot(pca_data, aes(x = PC1, y = PC2, color = factor(Group))) +
  114. geom_point(size = 2, alpha = 0.7) +
  115. labs(
  116. x = paste0("PC1 (", pc1_var, "%)"),
  117. y = paste0("PC2 (", pc2_var, "%)"),
  118. title = title,
  119. color = "Group"
  120. ) +
  121. theme_minimal()
  123. return(list(plot = p, pca = pca, data = pca_data))
  124. }
  126. # Helper function for correlation heatmap
  127. create_correlation_heatmap <- function(exp_data, output_dir, filename_prefix) {
  128. cor_matrix <- cor(t(exp_data), use = "pairwise.complete.obs")
  129. p_matrix <- cor.mtest(cor_matrix)$p
  131. genes <- colnames(cor_matrix)
  133. cor_df <- as.data.frame(as.table(cor_matrix), stringsAsFactors = FALSE) %>%
  134. rename(gene_y = Var1, gene_x = Var2, r = Freq)
  136. p_df <- as.data.frame(as.table(p_matrix), stringsAsFactors = FALSE) %>%
  137. rename(gene_y = Var1, gene_x = Var2, p = Freq)
  139. plot_df <- cor_df %>%
  140. left_join(p_df, by = c("gene_y", "gene_x")) %>%
  141. mutate(
  142. sig = case_when(
  143. p < 0.001 ~ "***",
  144. p < 0.01 ~ "**",
  145. p < 0.05 ~ "*",
  146. TRUE ~ ""
  147. ),
  148. label = ifelse(gene_x == gene_y, sprintf("%.2f", r), sprintf("%.2f%s", r, sig))
  149. )
  151. plot_df$gene_x <- factor(plot_df$gene_x, levels = genes)
  152. plot_df$gene_y <- factor(plot_df$gene_y, levels = rev(genes))
  154. p <- ggplot(plot_df, aes(x = gene_x, y = gene_y, fill = r)) +
  155. geom_tile(color = "white", linewidth = 0.4) +
  156. geom_text(aes(label = label), size = 2.8, color = "black") +
  157. scale_fill_gradient2(
  158. low = "#8491B4", mid = "white", high = "#F39B7F",
  159. midpoint = 0, limits = c(-1, 1), name = "Correlation"
  160. ) +
  161. scale_x_discrete(position = "top") +
  162. coord_fixed() +
  163. theme_minimal(base_size = 11) +
  164. theme(
  165. panel.grid = element_blank(),
  166. axis.title = element_blank(),
  167. axis.text.x.top = element_text(face = "italic", color = "black",
  168. angle = 45, hjust = 0, vjust = 0,
  169. margin = margin(b = 0)),
  170. axis.text.y.left = element_text(face = "italic", color = "black",
  171. margin = margin(r = 0))
  172. )
  174. ggsave(file.path(output_dir, paste0(filename_prefix, "_corrplot_full.pdf")),
  175. p, width = 4.5, height = 4.5, dpi = 300)
  176. ggsave(file.path(output_dir, paste0(filename_prefix, "_corrplot_full.png")),
  177. p, width = 4.5, height = 4.5, dpi = 300)
  179. return(list(plot = p, cor_matrix = cor_matrix, p_matrix = p_matrix))
  180. }
  182. #### Load and prepare data ####
  183. # Load expression data
  184. exp01 <- read.csv(file.path("00_rawdata", "00.rawdata_GSE139602_exp.csv"), row.names = 1)
  185. group01 <- read.csv(file.path("00_rawdata", "00.rawdata_GSE139602_group.csv"), row.names = 1)
  186. group01 <- group01[colnames(exp01), , drop = FALSE]
  188. # Convert group labels
  189. group01$group <- ifelse(group01$characteristics_ch1 == "disease state: Healthy", 0,
  190. ifelse(group01$characteristics_ch1 == "disease state: eCLD", 1,
  191. ifelse(group01$characteristics_ch1 == "disease state: Compensated Cirrhosis", 2,
  192. ifelse(group01$characteristics_ch1 == "disease state: Decompesated Cirrhosis", 3,
  193. ifelse(group01$characteristics_ch1 == "disease state: Acute-on-chronic liver failure", 4, NA)))))
  195. exp02 <- read.csv(file.path("00_rawdata", "00.rawdata_GSE15654_exp.csv"), row.names = 1)
  196. group02 <- read.csv(file.path("00_rawdata", "00.rawdata_GSE15654_group.csv"), row.names = 1)
  198. # Load models and gene lists
  199. model1 <- readRDS(file.path("07_xgboost", "07_xgboost_GSE139602_final_xgboost_model.Rdata"))
  200. model2 <- readRDS(file.path("07_xgboost", "07_xgboost_GSE15654_final_xgboost_model_cox.Rdata"))
  202. model1_gene <- read.csv(file.path("07_xgboost", "07_xgboost_GSE139602_shap_mat.csv"), row.names = 1)
  203. model2_gene <- read.csv(file.path("07_xgboost", "07_xgboost_GSE15654_shap_mat.csv"), row.names = 1)
  205. all_features <- c(colnames(model1_gene), colnames(model2_gene))
  207. # Find common genes
  208. common_genes <- intersect(rownames(exp01), rownames(exp02))
  209. exp01 <- exp01[common_genes, ]
  210. exp02 <- exp02[common_genes, ]
  211. exp_merged <- cbind(exp01, exp02)
  213. # Prepare group information
  214. group01$project <- "GSE139602"
  215. group02$project <- "GSE15654"
  216. group_merged <- data.frame(
  217. row.names = c(rownames(group01), rownames(group02)),
  218. group = c(group01$group, group02$group),
  219. project = c(group01$project, group02$project)
  220. )
  221. group_merged <- group_merged[colnames(exp_merged), ]
  223. #### Batch correction ####
  224. # Pre-correction visualization
  225. png(file.path(output, "09_SEM_pre_merged.png"), width = 5, height = 4, res = 300, units = "in")
  226. boxplot(exp_merged, xaxt = "n", col = "lightblue",
  227. main = "Gene Expression Distribution (Pre-batch Correction)",
  228. ylab = "Expression Value")
  229. dev.off()
  231. pdf(file.path(output, "09_SEM_pre_merged.pdf"), width = 5, height = 4)
  232. boxplot(exp_merged, xaxt = "n", col = "lightblue",
  233. main = "Gene Expression Distribution (Pre-batch Correction)",
  234. ylab = "Expression Value")
  235. dev.off()
  237. # Perform batch correction
  238. group_merged$sample <- rownames(group_merged)
  239. sample_order <- colnames(exp_merged)
  240. batch_info <- group_merged$project[match(sample_order, group_merged$sample)]
  241. group_info <- group_merged$group[match(sample_order, group_merged$sample)]
  242. names(batch_info) <- sample_order
  243. names(group_info) <- sample_order
  245. exp_corrected <- perform_batch_correction(exp_merged, batch_info, group_info)
  247. # Post-correction visualization
  248. png(file.path(output, "09_SEM_post_merged.png"), width = 5, height = 4, res = 300, units = "in")
  249. boxplot(exp_corrected, xaxt = "n", col = "lightblue",
  250. main = "Gene Expression Distribution (Post-batch Correction)",
  251. ylab = "Expression Value")
  252. dev.off()
  254. pdf(file.path(output, "09_SEM_post_merged.pdf"), width = 5, height = 4)
  255. boxplot(exp_corrected, xaxt = "n", col = "lightblue",
  256. main = "Gene Expression Distribution (Post-batch Correction)",
  257. ylab = "Expression Value")
  258. dev.off()
  260. # Save corrected data
  261. write.csv(group_merged, file.path(output, "09_SEM_group_merged.csv"))
  262. write.csv(exp_merged, file.path(output, "09_SEM_exp_merged.csv"))
  263. write.csv(exp_corrected, file.path(output, "09_SEM_exp_corrected.csv"))
  265. # PCA plots
  266. pca_pre <- create_pca_plot(exp_merged, group_merged$group, group_merged$project,
  267. "PCA of Pre-corrected Expression Data")
  268. ggsave(file.path(output, "09_SEM_pca_pre.png"), pca_pre$plot, width = 6, height = 5)
  269. ggsave(file.path(output, "09_SEM_pca_pre.pdf"), pca_pre$plot, width = 6, height = 5)
  271. pca_post <- create_pca_plot(exp_corrected, group_merged$group, group_merged$project,
  272. "PCA of Post-corrected Expression Data")
  273. ggsave(file.path(output, "09_SEM_pca_post.png"), pca_post$plot, width = 6, height = 5)
  274. ggsave(file.path(output, "09_SEM_pca_post.pdf"), pca_post$plot, width = 6, height = 5)
  276. #### Normality testing ####
  277. exp_data <- exp_corrected[all_features, ]
  278. norm_res <- test_normality(exp_data, by_group = FALSE)
  279. write.csv(norm_res, file.path(output, "09_SEM_normality_shapiro_overall.csv"))
  281. #### Correlation analysis ####
  282. cor_results <- create_correlation_heatmap(exp_data, output, "09_SEM")
  284. #### Mediation analysis ####
  285. # Prepare data for mediation
  286. X <- as.data.frame(t(exp_corrected[all_features, ]))
  287. X <- scale(X)
  288. X <- as.data.frame(X)
  290. # Predict progression and survival
  291. progression <- predict(model1, scale(t(exp_corrected[colnames(model1_gene), ])))
  292. survival <- predict(model2, scale(t(exp_corrected[colnames(model2_gene), ])))
  294. X$progression <- progression
  295. X$survival <- survival
  297. # Define mediator and outcome models
  298. mediator_vars <- c("NPC2", "TLN1", "TUBA1C", "LRRC32", "PRB2")
  299. outcome_vars <- c("SOX9", "RNASE4", "SERPINA3")
  301. # Linear regression models for mediation
  302. model_a <- lm(progression ~ NPC2 + TLN1 + TUBA1C + LRRC32 + PRB2, data = X)
  303. model_bc <- lm(survival ~ NPC2 + TLN1 + TUBA1C + LRRC32 + PRB2 + SOX9 + RNASE4 + SERPINA3 + progression, data = X)
  305. # Perform mediation analysis
  306. mediation_results <- list()
  307. for (treat in mediator_vars) {
  308. med_result <- mediate(model.m = model_a,
  309. model.y = model_bc,
  310. treat = treat,
  311. mediator = "progression",
  312. boot = TRUE,
  313. sims = 1000)
  314. mediation_results[[treat]] <- summary(med_result)
  315. }
  317. #### Structural Equation Modeling (SEM) ####
  318. # Define SEM model
  319. sem_model <- '
  320. # Direct effects on progression
  321. progression ~ a_NPC2*NPC2 + a_TLN1*TLN1 + a_TUBA1C*TUBA1C + a_LRRC32*LRRC32 + a_PRB2*PRB2
  323. # Effects on outcome mediators
  324. SOX9 ~ b_SOX9*progression
  325. RNASE4 ~ b_RNASE4*progression
  326. SERPINA3 ~ b_SERPINA3*progression
  328. # Effects on survival
  329. survival ~ c_SOX9*SOX9 + c_RNASE4*RNASE4 + c_SERPINA3*SERPINA3 + c_progression*progression
  331. # Indirect effects
  332. indirect_NPC2_SOX9 := a_NPC2 * b_SOX9 * c_SOX9
  333. indirect_NPC2_RNASE4 := a_NPC2 * b_RNASE4 * c_RNASE4
  334. indirect_NPC2_SERPINA3 := a_NPC2 * b_SERPINA3 * c_SERPINA3
  336. indirect_TLN1_SOX9 := a_TLN1 * b_SOX9 * c_SOX9
  337. indirect_TLN1_RNASE4 := a_TLN1 * b_RNASE4 * c_RNASE4
  338. indirect_TLN1_SERPINA3 := a_TLN1 * b_SERPINA3 * c_SERPINA3
  340. indirect_TUBA1C_SOX9 := a_TUBA1C * b_SOX9 * c_SOX9
  341. indirect_TUBA1C_RNASE4 := a_TUBA1C * b_RNASE4 * c_RNASE4
  342. indirect_TUBA1C_SERPINA3 := a_TUBA1C * b_SERPINA3 * c_SERPINA3
  344. indirect_LRRC32_SOX9 := a_LRRC32 * b_SOX9 * c_SOX9
  345. indirect_LRRC32_RNASE4 := a_LRRC32 * b_RNASE4 * c_RNASE4
  346. indirect_LRRC32_SERPINA3 := a_LRRC32 * b_SERPINA3 * c_SERPINA3
  348. indirect_PRB2_SOX9 := a_PRB2 * b_SOX9 * c_SOX9
  349. indirect_PRB2_RNASE4 := a_PRB2 * b_RNASE4 * c_RNASE4
  350. indirect_PRB2_SERPINA3 := a_PRB2 * b_SERPINA3 * c_SERPINA3
  351. '
  353. # Fit SEM model
  354. set.seed(123)
  355. fit <- sem(sem_model, data = X, se = "bootstrap", bootstrap = 2000, missing = "FIML")
  357. # Extract results
  358. fit_summary <- summary(fit, standardized = TRUE, rsquare = TRUE)
  359. fit_measures <- fitmeasures(fit, c("cfi", "tli", "rmsea", "srmr"))
  360. parameter_estimates <- parameterEstimates(fit, standardized = TRUE)
  362. # Save results
  363. write.csv(parameter_estimates, file.path(output, "09_SEM_sem_results.csv"))
  364. write.csv(fit_measures, file.path(output, "09_SEM_fit_measures.csv"))
  366. # Extract indirect effects
  367. indirect_effects <- parameter_estimates[grepl("^indirect_", parameter_estimates$label), ]
  368. indirect_effects <- indirect_effects %>%
  369. mutate(p_adj = p.adjust(pvalue, method = "fdr"))
  371. write.csv(indirect_effects, file.path(output, "09_SEM_indirect_effects.csv"))
  372. write.csv(indirect_effects[indirect_effects$pvalue < 0.05, ],
  373. file.path(output, "09_SEM_significant_indirect.csv"))
  375. #### Create SEM path diagram ####
  376. # Prepare node positions
  377. nodes <- bind_rows(
  378. tibble(node = c("NPC2", "TLN1", "TUBA1C", "LRRC32", "PRB2"),
  379. x = c(-4, -2, 0, 2, 4), y = 4),
  380. tibble(node = "progression", x = 0, y = 3),
  381. tibble(node = c("SOX9", "RNASE4", "SERPINA3"),
  382. x = c(-2, 0, 2), y = 2),
  383. tibble(node = "survival", x = 0, y = 1)
  384. )
  386. # Prepare edges
  387. pe <- parameter_estimates %>%
  388. filter(op == "~") %>%
  389. transmute(
  390. from = rhs,
  391. to = lhs,
  392. est = std.all,
  393. pval = pvalue
  394. )
  396. edges <- pe %>%
  397. inner_join(nodes %>% rename(from = node, x = x, y = y), by = "from") %>%
  398. inner_join(nodes %>% rename(to = node, xend = x, yend = y), by = "to") %>%
  399. mutate(
  400. dx = xend - x,
  401. curve_type = case_when(
  402. dx > 0 ~ "right",
  403. dx < 0 ~ "left",
  404. TRUE ~ "vertical"
  405. ),
  406. xm = (x + xend) / 2 + ifelse(dx == 0, 0.18, 0),
  407. ym = (y + yend) / 2 + ifelse(dx == 0, 0, 0.12 * sign(dx))
  408. )
  410. # Prepare node labels with italics for genes
  411. gene_nodes <- c("NPC2", "TLN1", "TUBA1C", "LRRC32", "PRB2", "SOX9", "RNASE4", "SERPINA3")
  412. nodes <- nodes %>%
  413. mutate(
  414. label_expr = ifelse(
  415. node %in% gene_nodes,
  416. paste0("italic('", node, "')"),
  417. paste0("'", node, "'")
  418. )
  419. )
  421. # Create SEM path diagram
  422. p_sem <- ggplot() +
  423. geom_curve(
  424. data = edges %>% filter(curve_type == "right"),
  425. aes(x = x, y = y, xend = xend, yend = yend, color = est, linewidth = abs(est)),
  426. curvature = 0.25, alpha = 0.9, lineend = "round",
  427. arrow = arrow(length = unit(0.18, "cm"), type = "closed")
  428. ) +
  429. geom_curve(
  430. data = edges %>% filter(curve_type == "left"),
  431. aes(x = x, y = y, xend = xend, yend = yend, color = est, linewidth = abs(est)),
  432. curvature = -0.25, alpha = 0.9, lineend = "round",
  433. arrow = arrow(length = unit(0.18, "cm"), type = "closed")
  434. ) +
  435. geom_curve(
  436. data = edges %>% filter(curve_type == "vertical"),
  437. aes(x = x, y = y, xend = xend, yend = yend, color = est, linewidth = abs(est)),
  438. curvature = 0.20, alpha = 0.9, lineend = "round",
  439. arrow = arrow(length = unit(0.18, "cm"), type = "closed")
  440. ) +
  441. geom_point(
  442. data = nodes,
  443. aes(x = x, y = y),
  444. shape = 21, size = 14, stroke = 0.9,
  445. fill = alpha("#F2F2F2", 0.45), color = "#4D4D4D"
  446. ) +
  447. geom_label_repel(
  448. data = edges,
  449. aes(x = xm, y = ym, label = sprintf("%.2f", est)),
  450. size = 3.0, color = "black", fill = "#F7F7F7",
  451. alpha = 0.95, label.size = 0.15, label.r = unit(0.08, "lines"),
  452. box.padding = 0.12, point.padding = 0.05,
  453. min.segment.length = 0, segment.color = NA,
  454. segment.size = 0.25, force = 1.2,
  455. max.overlaps = Inf, show.legend = FALSE
  456. ) +
  457. geom_text(
  458. data = nodes,
  459. aes(x = x, y = y, label = label_expr),
  460. parse = TRUE, size = 4.3, fontface = "bold", color = "#1F1F1F"
  461. ) +
  462. scale_color_gradient2(
  463. low = "#2C7BB6", mid = "#BDBDBD", high = "#D7191C", midpoint = 0,
  464. name = "Standardized coefficient"
  465. ) +
  466. scale_linewidth(range = c(0.7, 2.8), name = "|Standardized coefficient|") +
  467. coord_cartesian(xlim = c(-5, 5), ylim = c(0.6, 4.4), clip = "off") +
  468. theme_void(base_size = 13) +
  469. theme(
  470. legend.position = "right",
  471. plot.margin = margin(10, 20, 10, 20)
  472. )
  474. ggsave(file.path(output, "09_SEM_semplot.png"), p_sem, width = 7.5, height = 6, dpi = 300)
  475. ggsave(file.path(output, "09_SEM_semplot.pdf"), p_sem, width = 7.5, height = 6, dpi = 300)
  477. message("SEM analysis completed successfully!")
  479. # Print summary
  480. cat("\n=== SEM Fit Measures ===\n")
  481. print(fit_measures)
  483. cat("\n=== Significant Indirect Effects ===\n")
  484. print(indirect_effects[indirect_effects$pvalue < 0.05, c("label", "est", "pvalue", "p_adj")])

r.09_SEM.R at commit b1fc48c, under MIT · at the source

Overview

Authors: Yuanfeng Lan1,2,3,4, Tian Zhao1,2,3, Ying Xu4, Haihong Ye1,2,3
ORCID iDs: Yuanfeng Lan, Haihong Ye
  1. Department of Medical Genetics and Developmental Biology, School of Basic Medical Sciences, Capital Medical University, Beijing 100069, China; (Y.L.); (T.Z.)
  2. Laboratory for Clinical Medicine, Capital Medical University, Beijing 100069, China
  3. Beijing Key Laboratory of Cell and Gene Therapy in Otology, Beijing 100069, China
  4. Department of Human Cell Biology and Genetics, SUSTech Homeostatic Medicine Institute, School of Medicine, Southern University of Science and Technology, Shenzhen 518055, China
Institutions: Capital Medical University (China); Southern University of Science and Technology (China)
Journal: International journal of molecular sciences, volume 27, issue 15, article 6925
Dates: received 1 July 2026; accepted 30 July 2026; published online 1 August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3390/ijms27156925 · PMID 42589578 · PMCID PMC13467542 · OpenAlex W7172247900
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), other condition (population)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing
Keywords: cirrhosis, hepatic encephalopathy, machine learning, prediction model, prognostic analysis
MeSH: Biomarkers*, Hepatic Encephalopathy*, Liver Cirrhosis*, Boosting Machine Learning Algorithms, Gene Expression Profiling, Humans, Prognosis, Transcriptome (* major topic)
Topic: Liver Disease and Transplantation (Hepatology, Medicine), according to OpenAlex
Funding: Beijing Municipal Natural Science Foundation (7262004)
Citations: not cited yet (Europe PMC); 39 references in the paper

Abstract

Cirrhosis, accounting for 2.4% of global mortality in 2019, represents a leading cause of death in chronic liver disease. Hepatic encephalopathy (HE), a decompensated complication of cirrhosis, is associated with a median survival of only 0.92 years post-diagnosis. Current screening methods relying on neuropsychological tests (e.g., Psychometric Hepatic Encephalopathy Score, PHES) have limitations such as time-consuming procedures and subjective interpretation, potentially delaying diagnosis. To address this, we integrated four cirrhotic transcriptomic cohorts (GSE41919, GSE57193, GSE139602, and GSE15654) and employed an integrated algorithm (LASSO [Least Absolute Shrinkage and Selection Operator]–RFE [Recursive Feature Elimination]–random forest) to identify HE-specific biomarker genes. Ultimately, we developed an HE risk-prediction system centered on eight HE-specific marker genes, namely, PRB2, TUBA1C, NPC2, LRRC32, TLN1, SOX9, SERPINA3 and RNASE4. Based on these genes, an XGBoost (eXtreme Gradient Boosting)-based HE risk stratification model was constructed, and SHAP (SHapley Additive exPlanations) analysis was further introduced to address the “black-box” limitation of conventional machine learning models and to improve the interpretability. The finalized eight-gene system enables accurate, efficient, and interpretable HE risk assessment in patients with cirrhosis. Functional characterization through gene set enrichment analysis and structural equation modeling further revealed that these marker genes converge on four interconnected biological processes, namely, metabolic homeostasis, synaptic and neural transmission, immune inflammatory signaling, and hepatic detoxification, which collectively reflect the gut–liver–brain axis disruption central to HE pathogenesis. This dual-model system, incorporating both cirrhosis progression and survival prognosis, provides a reliable and clinically applicable tool for early HE risk warning and stratification, reducing the limitations of traditional neuropsychological screening and offering a translational foundation for timely intervention and prognostic optimization in high-risk cirrhotic patients.

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

Repositories

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

seandavi/GEOquery

License: MIT
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: cb12423bf2691af82f742fc10f891835e8bab3ce, 17 August 2026
Languages: R (53), Quarto (6)
Size: 154 files, 59 scripts
Software Heritage: not archived
Found in: the text, “4.1. Data Acquisition and Cohort Construction”
Holds: README, license file, CITATION.cff, environment (DESCRIPTION, Dockerfile), tests, continuous integration, documentation, 6 notebooks
Tools: tidyverse (6 files), data.table (3 files), limma (2 files), SingleCellExperiment (2 files), Seurat (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
62 files
At the source: github.com/seandavi/GEOquery

ModelOriented/shapviz

License: GPL-2.0
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: a3ead2491c66b7345eaf23e212eb1edda03b4697, 31 August 2026
Languages: R (28)
Size: 101 files, 28 scripts
Software Heritage: not archived
Found in: the text, “4.6. Construction and Validation of XGBoost–SHAP”
Holds: README, license file, environment (DESCRIPTION), tests, continuous integration, documentation, 4 notebooks
Not found: CITATION.cff
Tools: ggplot2 (12 files), XGBoost (9 files), patchwork (7 files), LightGBM (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
30 files
At the source: github.com/ModelOriented/shapviz

alserglab/fgsea

License: MIT
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 570f5903545835cb340e699fcdfdf444b1ba07ec, 22 September 2026
Languages: R (33), C++ (10), C/C++ (8)
Size: 107 files, 51 scripts
Software Heritage: archived
Found in: the text, “4.8. Gene Set Enrichment Analysis (GSEA)”
Holds: README, license file, environment (DESCRIPTION), tests, continuous integration, documentation, 2 notebooks
Not found: CITATION.cff
Tools: data.table (12 files), ggplot2 (4 files), limma (4 files), cowplot (3 files), Seurat (2 files)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
53 files
At the source: github.com/alserglab/fgsea

yuanfeng-lan/HE-risk-Prediction

License: MIT
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: b1fc48c5cbe5977c65ab0cb4a0b9cccdeb51e2d4, 24 July 2026
Languages: R (13)
Size: 17 files, 13 scripts
Software Heritage: not archived
Found in: the text, “4.9. Statistical Analysis and Data Visualization”
Holds: README, environment (renv.lock)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: ggplot2 (11 files), tidyverse (10 files), survival (4 files), caret (3 files), cowplot (3 files), XGBoost (3 files), data.table (2 files), limma (2 files), patchwork (2 files), clusterProfiler (1 file), ggpubr (1 file), glmnet (1 file), lavaan (1 file), pheatmap (1 file), pROC (1 file), randomForest (1 file), reshape2 (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
14 files
At the source: github.com/yuanfeng-lan/HE-risk-Prediction

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:

  • 4 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 151 scripts, each with its path and the digest of its content;
  • 18 matches 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

Datasets cited

Other data links

Data Availability Statement

All data utilized in this study were obtained from the Gene Expression Omnibus (GEO) database, a publicly accessible repository maintained by the National Center for Biotechnology Information (NCBI). All datasets are freely available for download without any access restrictions. The specific GEO accession numbers and corresponding datasets are cited within the manuscript.

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, 4 authors, 5 keywords, 8 MeSH terms, 1 funder, 38 references.

Cite

This paper

Lan, Y., Zhao, T., Xu, Y., & Ye, H. (2026). XGBoost-SHAP Interpretable Modeling Identifies and Validates an Eight-Gene Biomarker for Hepatic Encephalopathy Risk Prediction in Cirrhosis. International journal of molecular sciences, 27(15), 6925. https://doi.org/10.3390/ijms27156925

BibTeX

@article{lan2026xgboost,
author = {Lan, Yuanfeng and Zhao, Tian and Xu, Ying and Ye, Haihong},
title = {{XGBoost-SHAP Interpretable Modeling Identifies and Validates an Eight-Gene Biomarker for Hepatic Encephalopathy Risk Prediction in Cirrhosis}},
journal = {International journal of molecular sciences},
year = {2026},
month = aug,
volume = {27},
number = {15},
pages = {6925},
publisher = {Multidisciplinary Digital Publishing Institute (MDPI)},
issn = {1422-0067},
doi = {10.3390/ijms27156925},
url = {https://doi.org/10.3390/ijms27156925},
pmid = {42589578},
pmcid = {PMC13467542}
}

RIS

TY - JOUR
AU - Lan, Yuanfeng
AU - Zhao, Tian
AU - Xu, Ying
AU - Ye, Haihong
TI - XGBoost-SHAP Interpretable Modeling Identifies and Validates an Eight-Gene Biomarker for Hepatic Encephalopathy Risk Prediction in Cirrhosis
T2 - International journal of molecular sciences
J2 - Int J Mol Sci
PY - 2026
DA - 2026/08/01
VL - 27
IS - 15
SP - 6925
SN - 1422-0067
PB - Multidisciplinary Digital Publishing Institute (MDPI)
DO - 10.3390/ijms27156925
UR - https://doi.org/10.3390/ijms27156925
LA - en
ER -

CSL-JSON

{
"id": "10.3390/ijms27156925",
"type": "article-journal",
"title": "XGBoost-SHAP Interpretable Modeling Identifies and Validates an Eight-Gene Biomarker for Hepatic Encephalopathy Risk Prediction in Cirrhosis",
"container-title": "International journal of molecular sciences",
"author": [
{
"family": "Lan",
"given": "Yuanfeng"
},
{
"family": "Zhao",
"given": "Tian"
},
{
"family": "Xu",
"given": "Ying"
},
{
"family": "Ye",
"given": "Haihong"
}
],
"container-title-short": "Int J Mol Sci",
"volume": "27",
"issue": "15",
"page": "6925",
"DOI": "10.3390/ijms27156925",
"PMID": "42589578",
"PMCID": "PMC13467542",
"ISSN": "1422-0067",
"publisher": "Multidisciplinary Digital Publishing Institute (MDPI)",
"URL": "https://doi.org/10.3390/ijms27156925",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
1
]
]
}
}

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.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: randomForest, glmnet, survival, 12 other tools
[2] doi:10.3390/ijms27104466 [code]
Uncovering the Key Circuit FOSL2/FOS/EGR3/EGR1, Contributing to the Hyperexcitability of Excitatory Neurons in the Epileptic Temporal Cortex and Hippocampus.
Journal: International journal of molecular sciences
In common: glmnet, pROC, caret, 11 other tools, genetics / omics, 2 references
[3] doi:10.1016/j.xcrm.2026.102682 [code]
TET CpG sequence-context-specific DNA demethylation shapes progression of IDH-mutant gliomas.
Journal: Cell reports. Medicine
In common: glmnet, pROC, survival, 11 other tools, genetics / omics, other condition, 1 reference
[4] doi:10.7717/peerj.21426 [code]
Integrated transcriptomic identification and validation reveal key autophagy-associated biomarkers in sleep deprivation.
Journal: PeerJ
In common: randomForest, glmnet, pROC, 11 other tools, genetics / omics
[5] 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: glmnet, pROC, survival, 11 other tools, other condition
[6] 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: randomForest, pROC, survival, 10 other tools, genetics / omics, other condition
[7] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: pROC, SingleCellExperiment, limma, 10 other tools, genetics / omics, other condition
[8] doi:10.3390/ijms27093997 [code]
Coordinated Multicellular Immune Programs and Drug Targets Revealed by Single-Cell Analysis in Driver-Mutated NSCLC.
Journal: International journal of molecular sciences
In common: glmnet, survival, clusterProfiler, 9 other tools, genetics / omics, other condition
[9] 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: glmnet, pROC, survival, 8 other tools, genetics / omics, other condition
[10] doi:10.1038/s41593-026-02367-0 [code]
A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: SingleCellExperiment, limma, clusterProfiler, 9 other tools, 1 reference

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.