#!/usr/bin/env Rscript
# 复现 Telles-de-Deus et al., Nature Microbiology 2026 (doi:10.1038/s41564-026-02302-w)
# Fig. 2a: RT-PCR Ct 值 vs YFV 病毒载量 (log10 RPM) 散点图 + 线性拟合与 Spearman 相关。
# 数据：ct_rpm.csv（论文补充材料 Supplementary Table 8 的真实数据）。
# 运行：~/miniforge3/bin/Rscript plot_ct_rpm.R  ->  输出 repro_r.png（300 dpi）

# 沙箱默认 locale=C 会把 UTF-8 按字节渲染，必须先切到 C.utf8
try(invisible(Sys.setlocale("LC_ALL", "C.utf8")), silent = TRUE)

suppressPackageStartupMessages({
  library(ggplot2)
  library(ragg)
})

df <- read.csv("ct_rpm.csv", stringsAsFactors = FALSE)
df$group <- factor(df$group, levels = c("I", "II", "III"))

COLS <- c(I = "#DE8A7B", II = "#8FBCD8", III = "#2F6B9A")
GLAB <- c(I = "I (MO)", II = "II (NHP)", III = "III (NHP)")

# Spearman 相关（与原文一致的统计量）
sp <- suppressWarnings(cor.test(df$ct, df$log10_rpm, method = "spearman"))
r_s <- unname(sp$estimate); p_s <- sp$p.value

# P 值排版为 "1.3 × 10⁻¹⁴" 形式
sup_digits <- c("-" = "\u207B", "0" = "\u2070", "1" = "\u00B9", "2" = "\u00B2",
                "3" = "\u00B3", "4" = "\u2074", "5" = "\u2075", "6" = "\u2076",
                "7" = "\u2077", "8" = "\u2078", "9" = "\u2079")
p_fmt <- sprintf("%.1e", p_s)
parts <- strsplit(p_fmt, "e")[[1]]
p_str <- paste0(parts[1], " \u00D7 10",
                paste0(sup_digits[strsplit(parts[2], "")[[1]]], collapse = ""))
annot <- paste0("Spearman's r = ", sprintf("%.2f", r_s), "\nP value = ", p_str)

p <- ggplot(df, aes(x = ct, y = log10_rpm)) +
  geom_smooth(aes(group = 1), method = "lm", se = TRUE,
              color = "black", fill = "#D9D9D9", linewidth = 0.8) +
  geom_point(aes(fill = group), shape = 21, size = 2.6,
             stroke = 0.7, color = "#333333", alpha = 0.9) +
  scale_fill_manual(values = COLS, labels = GLAB, name = "Group:") +
  scale_x_continuous(breaks = c(10, 20, 30)) +
  scale_y_continuous(breaks = c(4, 5, 6, 7)) +
  coord_cartesian(xlim = c(4.5, 33.5), ylim = c(3.2, 7.3)) +
  labs(x = "Ct value", y = "Viral load (log10 RPM)", tag = "a") +
  annotate("text", x = 5.6, y = 3.55, label = annot, hjust = 0, vjust = 0, size = 4) +
  guides(fill = guide_legend(title.position = "top", nrow = 1)) +
  theme_classic(base_size = 13) +
  theme(
    legend.position = c(0.42, 0.93),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.title = element_text(size = 11),
    legend.text = element_text(size = 11),
    plot.tag = element_text(size = 20, face = "bold"),
    plot.tag.position = c(0.0, 1.0),
    axis.line = element_line(linewidth = 0.8)
  )

agg_png("repro_r.png", width = 6.6, height = 5.6, units = "in", res = 300)
print(p)
invisible(dev.off())
cat(sprintf("n=%d, Spearman r=%.3f, P=%.2e -> repro_r.png OK\n", nrow(df), r_s, p_s))