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

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: 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 (how the code is found and 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 12 November 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: 2 checks, the latest on 12 November 2026: the link answers
  • 12 November 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.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: pROC, glmnet, caret, 11 other tools, genetics / omics, 2 references
[2] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: randomForest, survival, glmnet, 12 other tools
[3] doi:10.1038/s41698-025-00909-6 [code]
Diagnostic and prognostic potential of cell-free RNAs in cerebrospinal fluid and plasma for brain tumors.
Journal: NPJ precision oncology
In common: randomForest, pROC, survival, 10 other tools, genetics / omics, other condition, 1 reference
[4] doi:10.1038/s42003-025-09115-3 [code]
Characterizing and decoding ultraconserved regions uncovers their regulatory significance in human brain development and disorders.
Journal: Communications biology
In common: pROC, survival, caret, 12 other tools
[5] 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: pROC, survival, glmnet, 11 other tools, genetics / omics, other condition, 1 reference
[6] doi:10.1038/s41467-025-66109-9 [code]
Mid-gestational cell-type-specific transcriptomic signatures in the prefrontal and superior temporal cortex in Down syndrome.
Journal: Nature communications
In common: randomForest, pROC, caret, 11 other tools, genetics / omics, other condition
[7] doi:10.7717/peerj.21426 [code]
Integrated transcriptomic identification and validation reveal key autophagy-associated biomarkers in sleep deprivation.
Journal: PeerJ
In common: randomForest, pROC, glmnet, 11 other tools, genetics / omics
[8] 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, glmnet, 11 other tools, other condition
[9] doi:10.1126/sciadv.adu0726 [code]
CHCHD2 mutant mice link mitochondrial deficits to PD pathophysiology.
Journal: Science advances
In common: randomForest, pROC, caret, 11 other tools
[10] doi:10.1016/j.isci.2022.105276 [code]
Integrative genomic analysis facilitates precision strategies for glioblastoma treatment.
Journal: iScience
In common: pROC, survival, SingleCellExperiment, 10 other tools, other condition, 2 references

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.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.3390/ijms27156925.svg)](https://openscicode.org/badge/doi/10.3390/ijms27156925)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.3390/ijms27156925.svg
   :target: https://openscicode.org/badge/doi/10.3390/ijms27156925
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.3390/ijms27156925"><img src="https://openscicode.org/badge/doi/10.3390/ijms27156925.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

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