Python绘制顶刊同款散点拟合回归组合图(附数据和代码)

数据准备

# ==========================================
# 2. 数据处理：读取数据 + 计算残差（用于小图）
# ==========================================
def prepare_data(file_path):
    """把原始数据读进来，并算出拟合线需要的残差"""
    try:
        df = pd.read_csv(file_path, index_col=0)
    except FileNotFoundError:
        # 如果你还没准备好 csv，这里会生成假数据让你先能跑通
        data = {
            'length': np.random.uniform(10 ** 5, 10 ** 7, 500),
            'number_of_cds': np.random.uniform(500, 10000, 500),
            'type': np.random.choice(CONFIG["hue_order"], 500)
        }
        df = pd.DataFrame(data)

    # 线性回归公式：y = ax + b (在 log 空间里跑)
    def linear_func(x, a, b):
        return a * x + b

    # 算对数，因为生物数据通常服从幂律分布
    log_x = np.log(df["length"])
    log_y = np.log(df["number_of_cds"])

    # 拟合：得到斜率 a 和 截距 b
    popt, _ = sciopt.curve_fit(linear_func, log_x, log_y)

    # 计算残差：实际值和预测值差了多少
    df["log_cds"] = log_y
    df["log_length"] = log_x
    df["cds_fit"] = np.exp(linear_func(log_x, *popt))
    df["residuals_log"] = log_y - np.log(df["cds_fit"])

    return df, popt

# ==========================================
# 3. 绘图：主图（散点+边缘密度+回归线）
# ==========================================
def plot_main_joint(df, popt):
    g = sns.jointplot(
        data=df, x="length", y="number_of_cds", hue="type",
        kind="scatter", alpha=0.7, height=6,
        palette=CONFIG["colors"], hue_order=CONFIG["hue_order"],
        joint_kws=dict(s=25, edgecolor='w', linewidth=0.5),
        marginal_kws=dict(bw_adjust=0.2, fill=True)
    )

    # 消除顶部边缘图 Y 轴的科学计数法（如 1e-7）
    g.ax_marg_x.ticklabel_format(axis='y', style='plain')

    g.ax_joint.set(xscale="log", yscale="log", **CONFIG["plot_limits"])
    g.ax_joint.grid(True, which="major", ls="--", alpha=0.5, zorder=0)

    # 画回归线：在 log 空间生成点，再转回原始坐标
    x_range = np.exp(np.linspace(df["log_length"].min(), df["log_length"].max(), 100))
    y_fit = np.exp(popt[0] * np.log(x_range) + popt[1])
    g.ax_joint.plot(x_range, y_fit, color="#34495e", ls="--", lw=1.5, zorder=5)

    # 修改图例：加入样本量 N=xxx
    counts = df["type"].value_counts()
    handles, labels = g.ax_joint.get_legend_handles_labels()
    new_labels = [f"MHQ {l} (N={counts.get(l, 0)})" for l in labels if l in counts]
    g.ax_joint.legend(handles=handles, labels=new_labels, loc="upper left", fontsize=9)

    return g

# ==========================================
# 4. 绘图：右下角箱线图（显示残差分布 + 显著性标注）
def add_inset_boxplot(joint_grid, df):
    ax_ins = inset_axes(
        joint_grid.ax_joint,
        width="40%", height="22%",
        loc="lower right",
        bbox_to_anchor=(-0.03, 0.05, 1, 1),
        bbox_transform=joint_grid.ax_joint.transAxes
    )

    box_colors = [CONFIG["colors"][2], CONFIG["colors"][1], CONFIG["colors"][0]]
    box_args = dict(
        data=df, x="residuals_log", y="type", orient="h",
        order=["WGS", "MAG", "SAG"], palette=box_colors,
        linewidth=1, fliersize=1.5
    )
    sns.boxplot(**box_args, ax=ax_ins)

    pairs = [("WGS", "MAG"), ("WGS", "SAG"), ("MAG", "SAG")]
    annot = Annotator(ax_ins, pairs, **box_args)
    annot.configure(test="Mann-Whitney", text_format="star", loc="inside", verbose=0)
    annot.apply_and_annotate()

    ax_ins.set(xlabel=None, ylabel=None, yticks=[])
    ax_ins.set_title("Residuals (log-scale)", fontsize=9, fontweight='bold', pad=12)
    for s in ax_ins.spines.values():
        s.set_visible(False)

    rect = mpl.patches.FancyBboxPatch(
        (0, 0), 1, 1, transform=ax_ins.transAxes,
        facecolor="white", edgecolor="#bdc3c7", alpha=0.8,
        boxstyle="round,pad=0.1,rounding_size=0.05", clip_on=False, zorder=0
    )
    ax_ins.add_patch(rect)