#!/usr/bin/env python3
"""复现 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 相关。

数据来源：论文补充材料 Supplementary Table 8（ct_rpm.csv 由其清洗得到，
含 Ct、log10(RPM)、分组三列真实数据；原表 RPM 列即为 log10 尺度，
为 Smart9n 宏基因组测序的 YFV 病毒载量）。用真实数据计算得
Spearman r=-0.72, P=1.3e-14，与原文报道完全一致。

运行：python3 plot_ct_rpm.py   ->  输出 repro_python.png（300 dpi）
"""
import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy import stats

HERE = os.path.dirname(os.path.abspath(__file__))
df = pd.read_csv(os.path.join(HERE, "ct_rpm.csv"))

x = df["ct"].to_numpy(dtype=float)
y = df["log10_rpm"].to_numpy(dtype=float)

# --- Spearman 相关（与原文一致的统计量） ---
r_s, p_s = stats.spearmanr(x, y)

# --- OLS 线性拟合 + 95% 置信带 ---
slope, intercept, _, _, _ = stats.linregress(x, y)
n = len(x)
yhat = intercept + slope * x
s_err = np.sqrt(np.sum((y - yhat) ** 2) / (n - 2))
xbar = x.mean()
Sxx = np.sum((x - xbar) ** 2)
xg = np.linspace(5.0, 32.5, 300)  # 与原图一致：拟合线/置信带止于数据范围
yg = intercept + slope * xg
tcrit = stats.t.ppf(0.975, n - 2)
se_mean = s_err * np.sqrt(1.0 / n + (xg - xbar) ** 2 / Sxx)

# --- 绘图 ---
COLORS = {"I": "#DE8A7B", "II": "#8FBCD8", "III": "#2F6B9A"}  # MO 珊瑚 / NHP 浅蓝 / NHP 深蓝
LABELS = {"I": "I (MO)", "II": "II (NHP)", "III": "III (NHP)"}

fig, ax = plt.subplots(figsize=(6.6, 5.6))
ax.fill_between(xg, yg - tcrit * se_mean, yg + tcrit * se_mean,
                color="#D9D9D9", alpha=0.9, zorder=1, linewidth=0)
ax.plot(xg, yg, color="black", lw=1.6, zorder=2)
for g in ["III", "II", "I"]:  # 深色组先画，浅色/珊瑚在上层
    d = df[df["group"] == g]
    ax.scatter(d["ct"], d["log10_rpm"], s=30, c=COLORS[g], edgecolors="#333333",
               linewidths=0.7, alpha=0.9, zorder=3, label=LABELS[g])

ax.set_xlim(4.5, 33.5)
ax.set_ylim(3.2, 7.3)
ax.set_xticks([10, 20, 30])
ax.set_yticks([4, 5, 6, 7])
ax.set_xlabel("Ct value", fontsize=13)
ax.set_ylabel("Viral load (log10 RPM)", fontsize=13)

# 图例：顶部两行 —— "Group:" 独占一行，三组在下一行（与原图一致）
ax.text(0.30, 0.99, "Group:", transform=ax.transAxes, fontsize=11,
        va="top", ha="left")
leg = ax.legend(loc="upper left", bbox_to_anchor=(0.16, 0.94),
                ncol=3, frameon=False, fontsize=11,
                handletextpad=0.4, columnspacing=1.2,
                handlelength=1.2, borderpad=0.2)
for h in getattr(leg, "legend_handles", None) or leg.legendHandles:
    h.set_alpha(1.0)
    h.set_edgecolor("#333333")

# 左下角统计注释（数值由本脚本从真实数据算出）
mant, exp = f"{p_s:.1e}".split("e")
sup = str(int(exp)).translate(str.maketrans("-0123456789", "⁻⁰¹²³⁴⁵⁶⁷⁸⁹"))
ax.text(0.035, 0.06,
        f"Spearman's $r$ = {r_s:.2f}\n$P$ value = {mant} × 10{sup}",
        transform=ax.transAxes, fontsize=11, va="bottom", ha="left")

ax.text(-0.02, 1.04, "a", transform=ax.transAxes, fontsize=20,
        fontweight="bold", va="bottom", ha="right")

for spine in ["top", "right"]:
    ax.spines[spine].set_visible(False)
ax.tick_params(direction="out", length=4, labelsize=11)
ax.spines["left"].set_linewidth(1.0)
ax.spines["bottom"].set_linewidth(1.0)

fig.tight_layout()
fig.savefig(os.path.join(HERE, "repro_python.png"), dpi=300)
print(f"n={n}, Spearman r={r_s:.3f}, P={p_s:.2e} -> repro_python.png OK")