# 复现 Nature Communications 2026 (s41467-026-74406-0) Fig.1a（R/ggplot2 版）
# 同源 E 蛋白结合抗体纵向衰减轨迹：原发（P，虚线）vs 继发（S，实线）+ 95% CI 带
# 数据为依据原图轨迹模拟（明确标注）
try(invisible(Sys.setlocale("LC_ALL", "C.utf8")), silent = TRUE)
library(ggplot2)
library(ragg)
library(systemfonts)
register_font("NotoSansSC",
              plain = "/home/hatch/workspace/arbovirus-site/day1-figure/NotoSansSC-Regular.ttf")

set.seed(20261020)
isotypes <- c("IgG", "IgA", "IgM", "IgG1", "IgG2", "IgG3", "IgG4")
pos  <- c("30" = 0, "90" = 1, "180" = 2, "540" = 4)
kpos <- c(0, 1, 2, 4)
p_mean <- list(
  IgG  = c(5.0, 5.5, 4.9, 5.1), IgA = c(4.5, 2.5, 1.7, 2.2),
  IgM  = c(5.0, 0.0, -0.1, -0.2), IgG1 = c(5.5, 6.3, 5.6, 5.6),
  IgG2 = c(3.2, 5.0, 3.5, 4.0),  IgG3 = c(4.6, 3.0, 0.8, 0.7),
  IgG4 = c(-1.0, -0.7, -0.9, 1.0))
s_mean <- list(
  IgG  = c(6.4, 6.1, 5.3, 5.2), IgA = c(2.9, 2.3, 1.6, 1.5),
  IgM  = c(1.2, -0.5, -0.6, -0.5), IgG1 = c(7.1, 6.9, 5.9, 5.7),
  IgG2 = c(7.2, 6.5, 4.3, 4.0),  IgG3 = c(3.3, 1.8, 0.4, 0.2),
  IgG4 = c(0.3, 0.4, 0.9, 1.6))

xp <- seq(0, 4, length.out = 220)
interp1 <- function(v) approx(kpos, v, xout = xp)$y
ci_band <- function(mu) {
  n <- 400
  w <- 0.22 + 0.10 * abs(mu) / 8.0
  lo <- matrix(NA, n, length(xp)); hi <- matrix(NA, n, length(xp))
  for (i in seq_len(n)) {
    jit <- rnorm(length(xp), 0, w) + rnorm(1, 0, 0.08) * sin(xp / 2)
    lo[i, ] <- mu + jit - 0.18; hi[i, ] <- mu + jit + 0.18
  }
  data.frame(x = xp,
             lo = apply(lo, 2, quantile, 0.025),
             hi = apply(hi, 2, quantile, 0.975))
}
df <- do.call(rbind, lapply(isotypes, function(iso) {
  pm <- interp1(p_mean[[iso]]); sm <- interp1(s_mean[[iso]])
  cb_p <- ci_band(pm); cb_s <- ci_band(sm)
  data.frame(
    x = c(xp, xp, xp, xp),
    y = c(pm, sm, rep(NA_real_, 2 * length(xp))),
    ymin = c(cb_p$lo, cb_s$lo, cb_p$lo, cb_s$lo),
    ymax = c(cb_p$hi, cb_s$hi, cb_p$hi, cb_s$hi),
    group = factor(rep(c("P", "S"), each = length(xp)), levels = c("P", "S")),
    panel = factor(iso, levels = isotypes))
}))
line_df <- df[!is.na(df$y), ]
band_df <- df
band_df$y <- NULL
band_df <- band_df[!duplicated(band_df[, c("x", "group", "panel")]), ]

COL_P <- "#C9A227"; COL_S <- "#5A5A5A"

p <- ggplot() +
  geom_rect(data = data.frame(panel = factor(isotypes, levels = isotypes)),
            aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0),
            fill = "#E8E8E8", alpha = 0.6) +
  geom_ribbon(data = band_df, aes(x = x, ymin = ymin, ymax = ymax, fill = group),
              alpha = 0.30) +
  geom_line(data = line_df, aes(x = x, y = y, color = group, linetype = group),
            linewidth = 0.9) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.4) +
  facet_wrap(~ panel, nrow = 1, scales = "free_x") +
  scale_x_continuous(breaks = c(0, 1, 4), labels = c("30", "90", "540"),
                     limits = c(0, 4)) +
  scale_y_continuous(limits = c(-2, 9)) +
  scale_color_manual(values = c(P = COL_P, S = COL_S)) +
  scale_fill_manual(values = c(P = COL_P, S = COL_S)) +
  scale_linetype_manual(values = c(P = "dashed", S = "solid")) +
  labs(x = "Days post symptoms onset", y = "Log2(MFI) - BG",
       color = "Infection group", fill = "Infection group", linetype = "Infection group") +
  theme_bw(base_family = "NotoSansSC") +
  theme(panel.grid = element_blank(), legend.position = "right",
        strip.background = element_blank(), strip.text = element_text(size = 12),
        plot.margin = margin(20, 10, 10, 10))

agg_png("/home/hatch/workspace/arbovirus-site/day-2026-10-20-figure/repro_r.png",
        width = 2000, height = 560, res = 120)
print(p + labs(title = "a"))
invisible(dev.off())
cat("saved repro_r.png\n")