Nature 期刊图复现 | Python 实现多情景生态抗性六边形频数回归图

import osimport numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom scipy import stats# ==============================================================================# 主流期刊标准学术配色方案字典 (20种)# ==============================================================================COLOR_PALETTES = {    1:  {"cmap": "viridis",       "line": "#FF2E17", "ci": "#CCCCCC"},  # 鲜红拟合线 + 经典浅灰置信区间阴影，最符合您最初的视觉样式}# 设定当前使用的方案索引（选择1）SELECTED_PALETTE_INDEX = 1def plot_from_excel_only(excel_path='data.xlsx'):    output_dir = '图表'    output_image_name = 'functional_resistance_vs_redundancy.png'    output_path = os.path.join(output_dir, output_image_name)    if not os.path.exists(excel_path):        raise FileNotFoundError(f"未找到 '{excel_path}' 文件，请确保该文件已存在。")    os.makedirs(output_dir, exist_ok=True)    print("正在从 Excel 读取数据并实时计算统计指标...")    try:        sheets = {            'a': pd.read_excel(excel_path, sheet_name='Panel_a'),            'b': pd.read_excel(excel_path, sheet_name='Panel_b'),            'c': pd.read_excel(excel_path, sheet_name='Panel_c')        }    except Exception as e:        print(f"读取 Excel 失败，请确保工作表名称精确为: Panel_a, Panel_b, Panel_c。")        print(f"具体错误信息: {e}")        return    # 应用选定的期刊配色方案    palette = COLOR_PALETTES[SELECTED_PALETTE_INDEX]    cmap_chosen = palette["cmap"]    line_color = palette["line"]    ci_color = palette["ci"]    plt.rcParams['font.sans-serif'] = ['DejaVu Sans', 'Arial', 'Microsoft YaHei']    plt.rcParams['axes.unicode_minus'] = False    fig, axes = plt.subplots(1, 3, figsize=(16, 5.5), sharey=True)    fig.patch.set_facecolor('white')    fig.subplots_adjust(wspace=0.06)    for i, label in enumerate(['a', 'b', 'c']):        ax = axes[i]        df = sheets[label]        if 'Relative_redundancy' not in df.columns or 'Functional_resistance' not in df.columns:            raise KeyError(f"Sheet 'Panel_{label}' 中必须包含 'Relative_redundancy' 和 'Functional_resistance' 两列。")        x_val = df['Relative_redundancy'].dropna().values        y_val = df['Functional_resistance'].dropna().values        n_samples = len(x_val)        if n_samples < 2:            print(f"警告：Sheet 'Panel_{label}' 中的有效数据量过少，无法进行回归拟合。")            continue        ax.set_facecolor('white')        ax.grid(True, color='#E5E5E5', linestyle='-', linewidth=1.0)        ax.set_axisbelow(True)        # 使用选定的色系绘制六边形频数图        hb = ax.hexbin(x_val, y_val, gridsize=28, cmap=cmap_chosen, mincnt=1, edgecolors='none')        # 基于真实数据计算回归指标        slope, intercept, r_value, p_value, std_err = stats.linregress(x_val, y_val)        x_range = np.linspace(x_val.min(), x_val.max(), 100)        y_pred = slope * x_range + intercept        # 计算 95% 置信区间        y_actual_pred = slope * x_val + intercept        residual_std = np.sqrt(np.sum((y_val - y_actual_pred) ** 2) / (n_samples - 2))        x_mean = np.mean(x_val)        ss_x = np.sum((x_val - x_mean) ** 2)        se_y = residual_std * np.sqrt(1 / n_samples + (x_range - x_mean) ** 2 / ss_x)        ci = 1.96 * se_y        # 使用选定方案的置信区间颜色与回归线颜色进行绘制        ax.fill_between(x_range, y_pred - ci, y_pred + ci, color=ci_color, alpha=0.5, zorder=3)        ax.plot(x_range, y_pred, color=line_color, linewidth=2.5, zorder=4)        ax.text(0.05, 0.94, f"Coefficient estimate: {slope:.2f}", transform=ax.transAxes, fontsize=12, va='top')        p_text = "p < 0.001" if p_value < 0.001 else f"p = {p_value:.3f}"        ax.text(0.05, 0.86, p_text, transform=ax.transAxes, fontsize=12, va='top')        ax.text(0.92, 0.94, label, transform=ax.transAxes, fontsize=20, va='top', fontweight='bold')        ax.set_xlim(x_val.min() - 0.5, x_val.max() + 0.5)        ax.tick_params(axis='both', which='major', labelsize=12, colors='#333333')        for spine in ax.spines.values():            spine.set_color('#CCCCCC')            spine.set_linewidth(1)    all_y = np.concatenate([sheets['a']['Functional_resistance'].dropna(),                            sheets['b']['Functional_resistance'].dropna(),                            sheets['c']['Functional_resistance'].dropna()])    if len(all_y) > 0:        axes[0].set_ylim(all_y.min() - 0.5, all_y.max() + 0.5)    fig.text(0.5, 0.02, 'Relative redundancy', ha='center', va='center', fontsize=15)    fig.text(0.015, 0.5, 'Functional resistance', ha='center', va='center', rotation='vertical', fontsize=15)    plt.tight_layout()    fig.subplots_adjust(bottom=0.12, left=0.06, right=0.98, top=0.95)    plt.savefig(output_path, dpi=300)    plt.close()    print(f"绘图成功，最终图表已导出至: {output_path}")if __name__ == '__main__':    plot_from_excel_only('data.xlsx')