#!/usr/bin/env Rscript
# 复现 Nature Microbiology 2026 (doi:10.1038/s41564-026-02303-9) Fig. 6a:
# 慢性期 CHIKV 感染小鼠经 SRI-42718 治疗 7 天后踝关节组织 CHIKV 基因组
# 拷贝数（RT-qPCR，log10 转换）。
# 数据: chikv_genomes_fig6a.csv —— 论文补充材料 Source Data (MOESM10, sheet
# 'Fig 6A-Viral RNA Copies') 的 16 个真实样本原始值，非模拟。
try(invisible(Sys.setlocale("LC_ALL", "C.utf8")), silent = TRUE)
args <- commandArgs(trailingOnly = FALSE)
file_arg <- grep("^--file=", args, value = TRUE)
HERE <- if (length(file_arg)) dirname(normalizePath(sub("^--file=", "", file_arg[[1]]))) else getwd()

library(ggplot2)
library(ragg)
library(systemfonts)

register_font("NotoSansSC",
              plain = file.path(HERE, "NotoSansSC-Regular.ttf"))

df <- read.csv(file.path(HERE, "chikv_genomes_fig6a.csv"), stringsAsFactors = FALSE)
df$log10_copies <- log10(df$copies_per_ug_rna)
df$group <- factor(df$group, levels = c("Vehicle", "SRI-42718"))
df$x <- as.numeric(df$group) - 1  # 0, 1

grp_mean <- aggregate(log10_copies ~ group + x, df, mean)
p_val <- t.test(log10_copies ~ group, data = df, var.equal = TRUE)$p.value
cat(sprintf("Vehicle mean=%.4f, SRI-42718 mean=%.4f, t-test p=%.3e\n",
            grp_mean$log10_copies[1], grp_mean$log10_copies[2], p_val))

BLACK <- "#000000"; WHITE <- "#FFFFFF"

set.seed(42)
df$xj <- df$x + runif(nrow(df), -0.16, 0.16)

p <- ggplot(df, aes(x = xj, y = log10_copies)) +
  geom_point(aes(shape = group), size = 4.2, colour = BLACK,
             fill = WHITE, stroke = 1.1) +
  scale_shape_manual(values = c("Vehicle" = 21, "SRI-42718" = 22),
                     name = NULL,
                     labels = c("Vehicle", "SRI-42718")) +
  guides(shape = guide_legend(override.aes = list(fill = WHITE, size = 4.5),
                              nrow = 1)) +
  geom_segment(data = grp_mean,
               aes(x = x - 0.28, xend = x + 0.28,
                   y = log10_copies, yend = log10_copies),
               colour = BLACK, linewidth = 1.0) +
  # 显著性括号
  annotate("segment", x = 0.12, xend = 0.88, y = 6.92, yend = 6.92,
           colour = BLACK, linewidth = 0.8) +
  annotate("segment", x = 0.12, xend = 0.12, y = 6.86, yend = 6.92,
           colour = BLACK, linewidth = 0.8) +
  annotate("segment", x = 0.88, xend = 0.88, y = 6.86, yend = 6.92,
           colour = BLACK, linewidth = 0.8) +
  annotate("text", x = 0.5, y = 7.04, label = "****",
           family = "NotoSansSC", size = 5.5, colour = BLACK) +
  scale_x_continuous(limits = c(-0.55, 1.55), breaks = NULL) +
  scale_y_continuous(limits = c(4, 7.5), breaks = c(4, 5, 6, 7)) +
  labs(title = "CHIKV genomes",
       y = "CHIKV genomes per \u00b5g RNA\n(log\u2081\u2080)") +
  theme_minimal(base_family = "NotoSansSC", base_size = 13) +
  theme(
    panel.grid = element_blank(),
    axis.line.x = element_line(colour = BLACK, linewidth = 0.7),
    axis.line.y = element_line(colour = BLACK, linewidth = 0.7),
    axis.ticks.y = element_line(colour = BLACK),
    axis.text.y = element_text(colour = BLACK, size = 12),
    axis.title.y = element_text(colour = BLACK, size = 13),
    axis.title.x = element_blank(),
    plot.title = element_text(colour = BLACK, size = 15, hjust = 0.5,
                              margin = margin(b = 14)),
    plot.margin = margin(t = 20, r = 20, b = 10, l = 10),
    legend.position = c(0.5, 0.94),
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.text = element_text(colour = BLACK, size = 12.5),
    legend.key.size = unit(1.4, "lines")
  )

agg_png(file.path(HERE, "repro_r.png"), width = 4.0, height = 5.2, units = "in",
        res = 300, background = "white")
print(p)
# 面板标号 a
grid::grid.text("a", x = 0.02, y = 0.965, gp = grid::gpar(
  fontfamily = "NotoSansSC", fontsize = 22, fontface = "bold", col = BLACK))
invisible(dev.off())
cat("saved repro_r.png\n")