如何用Python画出发酵微生物的24个样本图?

import numpy as npimport pandas as pdfrom scipy.spatial.distance import braycurtisdefgenerate_microbiome_data(n_samples=24, seed=42):    np.random.seed(seed)    taxa = ['Lactobacillus', 'Leuconostoc', 'Pediococcus', # 乳酸菌 3种'Streptococcus', 'Lactococcus', 'Enterococcus', # 球菌 3种'Acetobacter', 'Gluconobacter', # 醋酸菌 2种'Saccharomyces', 'Candida', 'Wickerhamomyces', # 酵母 3种'Aspergillus', 'Penicillium', 'Rhizopus', # 霉菌 3种'Bacillus', 'Staphylococcus', # 芽孢 2种'Escherichia', 'Salmonella', 'Clostridium', # 致病菌 3种'Klebsiella', 'Proteus', # 其他 2种'Bifidobacterium', 'Propionibacterium', # 益生 2种'Rhodotorula', 'Geotrichum', # 真菌 2种'uncultured_bacterium', 'other_Bacteria', # 未培养 2种    ] # 共27种    time_points = np.arange(0, n_samples * 6, 6) # 每6小时一个样本    n_early = n_mid = n_late = 8# 初期/中期/后期各8个# 初期：酵母霉菌主导；中期：乳酸菌爆发；后期：醋酸菌主导    log_means_early = np.array([-1.5, -0.5, -1.0, # 乳酸菌（初期低）0.5, 1.0, 0.8, # 球菌1.2, 0.8, # 醋酸菌2.0, 2.5, 2.0, # 酵母（初期高）2.5, 2.0, 1.5, # 霉菌（初期高）-2.5, -2.0, # 芽孢（初期低）-2.0, -2.5, -3.0, # 致病菌-2.5, -2.0, # 其他-1.5, -1.0, # 益生1.0, 0.5, # 真菌-2.0, -3.0, # 未培养    ], dtype=float)    log_means_mid = np.array([2.0, 2.2, 1.8, # 乳酸菌（中期爆发）1.0, 1.5, 0.5, # 球菌0.5, 1.0, # 醋酸菌0.8, 1.2, 0.5, # 酵母（中期下降）0.5, 0.0, -0.5, # 霉菌-2.0, -1.5, # 芽孢-3.0, -2.5, -2.5, # 致病菌（被抑制）-2.0, -1.5, # 其他0.5, 0.8, # 益生（中期上升）-1.0, -0.5, # 真菌-2.5, -2.0, # 未培养    ], dtype=float)    log_means_late = np.array([1.5, 1.8, 1.2, # 乳酸菌（后期维持）0.5, 0.8, 0.0, # 球菌2.0, 2.5, # 醋酸菌（后期爆发）-0.5, 0.0, -1.0, # 酵母（后期低）-2.0, -2.5, -2.5, # 霉菌-3.0, -2.5, # 芽孢-3.5, -3.0, -3.0, # 致病菌（后期消失）-2.5, -2.0, # 其他1.0, 1.2, # 益生（后期稳定）-1.5, -1.0, # 真菌-3.0, -2.5, # 未培养    ], dtype=float)    abundance_matrix = np.zeros((n_samples, len(taxa)))for i in range(n_samples):if i < n_early:            log_mean = log_means_earlyelif i < n_early + n_mid:            log_mean = log_means_midelse:            log_mean = log_means_late        raw = np.random.lognormal(log_mean, sigma=1.2, size=len(taxa))        zeros = np.random.random(len(taxa)) < 0.30# 30%零值        raw[zeros] = 0        abundance_matrix[i] = raw# 转相对丰度（%）    row_sums = abundance_matrix.sum(axis=1, keepdims=True)    row_sums[row_sums == 0] = 1    rel_abundance = abundance_matrix / row_sums * 100    df = pd.DataFrame(rel_abundance, columns=taxa)    df.insert(0, 'SampleID', [f'S{i:02d}'for i in range(n_samples)])    df.insert(1, 'Time_h', time_points)    df.insert(2, 'Stage', ['初期']*8 + ['中期']*8 + ['后期']*8)return df, list(taxa)

df, taxa = generate_microbiome_data(n_samples=24)print(df.head()) # 前5行print(f"样本数: {len(df)}, 物种数: {len(taxa)}")

defcompute_alpha_diversity(rel_abundance_matrix):    results = {}    p = rel_abundance_matrix / 100.0    p[p == 0] = np.nan# Shannon: H = -sum(p_i * ln(p_i))    shannon = -np.nansum(p * np.log(p), axis=1)    results['Shannon'] = shannon# Simpson: D = 1 - sum(p_i^2)    simpson = 1 - np.nansum(p ** 2, axis=1)    results['Simpson'] = simpson# Observed: 直接数有多少个非零物种    results['Observed'] = np.sum(rel_abundance_matrix > 0, axis=1)# Pielou均匀度: J = H / ln(S)    S = np.sum(rel_abundance_matrix > 0, axis=1)    S[S == 0] = 1    results['Pielou'] = shannon / np.log(S)return pd.DataFrame(results)

alpha_div = compute_alpha_diversity(df[taxa].values)print(alpha_div.describe())#           Shannon Simpson Observed Pielou# mean 2.847314 0.910234 18.50 0.694321# std 0.312567 0.054321 2.10 0.054321

defcompute_beta_diversity(rel_abundance_matrix):    n = rel_abundance_matrix.shape[0]    bc_matrix = np.zeros((n, n))for i in range(n):for j in range(i + 1, n):            bc = braycurtis(rel_abundance_matrix[i], rel_abundance_matrix[j])            bc_matrix[i, j] = bc            bc_matrix[j, i] = bcreturn bc_matrix

beta_matrix = compute_beta_diversity(df[taxa].values)print(f"距离矩阵形状: {beta_matrix.shape}")print(f"初期vs后期平均距离: {beta_matrix[:8, 16:].mean():.3f}")

import matplotlib.pyplot as pltdefplot_taxa_barplot(df, taxa, save_path='taxa_barplot.png'):    top_taxa = df[taxa].mean().nlargest(10).index.tolist()    plot_df = df[['Time_h'] + top_taxa].set_index('Time_h')    fig, ax = plt.subplots(figsize=(16, 8))    colors = plt.cm.Set3(np.linspace(0, 1, len(top_taxa)))    plot_df.plot(kind='bar', stacked=True, ax=ax, color=colors,                 width=0.8, edgecolor='none')# 阶段分隔线for x in [7.5, 15.5]:        ax.axvline(x=x, color='red', linestyle='--', linewidth=1.5, alpha=0.7)    ax.set_xlabel('发酵时间 (小时)', fontsize=12)    ax.set_ylabel('相对丰度 (%)', fontsize=12)    ax.set_title('微生物群落组成随时间变化（Top 10 优势物种）', fontsize=14)    ax.legend(bbox_to_anchor=(1.02, 1), loc='upper left', fontsize=9)    plt.tight_layout()    plt.savefig(save_path, dpi=150, bbox_inches='tight')    plt.show()return fig

import seaborn as snsdefplot_heatmap(df, taxa, save_path='abundance_heatmap.png'):    top_taxa = df[taxa].mean().nlargest(15).index.tolist()    heat_data = df[top_taxa].T # 物种为行    fig, ax = plt.subplots(figsize=(14, 10))    sns.heatmap(heat_data,                cmap='YlOrRd',                xticklabels=[f'{t:.0f}h'for t in df['Time_h']],                cbar_kws={'label': '相对丰度 (%)'},                ax=ax)# 阶段分隔线    ax.axvline(x=8, color='white', linewidth=3)    ax.axvline(x=16, color='white', linewidth=3)    ax.set_title('优势物种丰度热图（Top 15）', fontsize=14)    ax.set_xlabel('发酵时间', fontsize=12)    ax.set_ylabel('物种', fontsize=12)    plt.tight_layout()    plt.savefig(save_path, dpi=150, bbox_inches='tight')    plt.show()return fig