# -*- coding: utf-8 -*-
# R 复现脚本：Nature Microbiology 2026, doi:10.1038/s41564-026-02465-6, Fig. 2b
# Serum anti-ZIKV E IgG (ELISA, log10) 与 NT50 (log10)，CON / WT / FLM 分组柱状图+散点。
#
# 数据说明（模拟数据！）：论文未随刊发布该图 Source Data（PMC13652353 无附件），
# zika_vaccine_fig2b.csv 由 Python 脚本按原图散点分布手工估读+固定种子生成（模拟数据），
# 仅用于教学复现作图代码，不代表论文原始数据；显著性 p 值引用自论文原文。
#
# 运行：~/miniforge3/bin/Rscript plot_zika_vaccine_fig2b.R
# Tested: R 4.5.3 + ggplot2 + ragg — sandbox run OK.
try(invisible(Sys.setlocale("LC_ALL", "C.utf8")), silent=TRUE)

library(ggplot2)
library(ragg)

HERE <- dirname(normalizePath("plot_zika_vaccine_fig2b.R"))
FONT <- normalizePath(file.path(HERE, "..", "day1-figure", "NotoSansSC-Regular.ttf"))
systemfonts::register_font("Noto Sans SC", FONT)

df <- read.csv(file.path(HERE, "zika_vaccine_fig2b.csv"), stringsAsFactors = FALSE)
df$group <- factor(df$group, levels = c("CON", "WT", "FLM"))
df$panel <- factor(df$panel, levels = c("igg", "nt50"),
                   labels = c("Anti-ZIKV E IgG (log10)", "NT50 (log10)"))
set.seed(20261019)
df$xj <- as.integer(df$group) + runif(nrow(df), -0.22, 0.22)

stats <- aggregate(value_log10 ~ panel + group, data = df,
                   FUN = function(x) c(mean = mean(x), sem = sd(x) / sqrt(length(x))))
stats <- data.frame(panel = stats$panel, group = stats$group,
                    mean = stats$value_log10[, "mean"], sem = stats$value_log10[, "sem"])
stats$fill <- ifelse(stats$group == "CON", "white", "black")

lod <- data.frame(panel = levels(df$panel), y = c(1.45, 1.00))
brackets <- data.frame(
  panel = rep(levels(df$panel), each = 2),
  x1 = c(1, 1, 1, 1), x2 = c(2, 3, 2, 3),
  ytop = c(5.10, 5.60, 3.95, 4.45), ytxt = c(4.82, 5.32, 3.68, 4.18),
  label = c("<0.0001", "<0.0001", "0.0010", "<0.0001"))

ymax <- data.frame(panel = levels(df$panel), y = c(5.9, 5.0))
br_df <- merge(brackets, ymax, by = "panel")
br_long <- do.call(rbind, lapply(seq_len(nrow(br_df)), function(i) {
  r <- br_df[i, ]
  data.frame(panel = r$panel, bid = i,
             x = c(r$x1, r$x1, r$x2, r$x2),
             y = c(r$ytxt + 0.18, r$ytop, r$ytop, r$ytxt + 0.18))
}))

p <- ggplot() +
  geom_blank(data = ymax, aes(x = 1, y = y)) +
  geom_col(data = stats, aes(x = as.integer(group), y = mean, fill = fill),
           width = 0.72, color = "black", linewidth = 0.5) +
  geom_errorbar(data = stats, aes(x = as.integer(group), ymin = mean - sem, ymax = mean + sem),
                width = 0.14, linewidth = 0.5) +
  geom_point(data = df, aes(x = xj, y = value_log10,
             fill = ifelse(point_type == "pooled", "red", "white")),
             shape = 21, size = 2.6, stroke = 0.7, color = "black") +
  scale_fill_identity() +
  geom_hline(data = lod, aes(yintercept = y), linetype = "dotted", linewidth = 0.5) +
  geom_line(data = br_long, aes(x = x, y = y, group = bid), linewidth = 0.5) +
  geom_text(data = brackets, aes(x = (x1 + x2) / 2, y = ytop + 0.06, label = label), size = 3.2) +
  facet_wrap(~panel, scales = "free_y", strip.position = "left") +
  scale_x_continuous(breaks = 1:3, labels = c("CON", "WT", "FLM")) +
  scale_y_continuous(breaks = 0:5, expand = expansion(mult = c(0, 0.02))) +
  labs(x = NULL, y = NULL, title = "b") +
  theme_classic(base_family = "Noto Sans SC", base_size = 12) +
  theme(plot.title = element_text(face = "bold", size = 20, hjust = 0, margin = margin(b = -8)),
        strip.background = element_blank(),
        strip.text = element_text(size = 11),
        axis.text.x = element_text(angle = 45, hjust = 1))

ggsave(file.path(HERE, "reproduce_r.png"), p, width = 8.4, height = 4.4, dpi = 300,
       device = agg_png)
cat("saved:", file.path(HERE, "reproduce_r.png"), "\n")