AI 辅助数据分析与可视化🔗

import pandas as pd
import numpy as np

# ============ 1. 读入数据 ============
df = pd.read_csv('enzyme_activity.csv', na_values=['#N/A', 'N/A', ''])

print("=== 原始数据概况 ===")
print(f"数据量: {df.shape[0]} 行, {df.shape[1]} 列")
print(f"\n数据类型:\n{df.dtypes}")
print(f"\n缺失值统计:\n{df.isnull().sum()}")
print(f"\n数值列描述性统计:\n{df.describe()}")

# ============ 2. 统一日期格式 ============
df['date'] = pd.to_datetime(df['date'], format='mixed')
print(f"\n日期范围: {df['date'].min()} ~ {df['date'].max()}")

# ============ 3. 处理负值（仪器错误） ============
n_negative = (df['absorbance_450'] < 0).sum()
print(f"\n吸光度负值数量: {n_negative}")
df.loc[df['absorbance_450'] < 0, 'absorbance_450'] = np.nan

# ============ 4. 处理缺失值 ============
n_missing_before = df['absorbance_450'].isnull().sum()
# 对每个 (sample_id, concentration) 组，用同组均值填充
df['absorbance_450'] = df.groupby(['sample_id', 'concentration'])[
    'absorbance_450'
].transform(lambda x: x.fillna(x.mean()))
n_missing_after = df['absorbance_450'].isnull().sum()
print(f"\n缺失值处理: {n_missing_before} → {n_missing_after}")

# 仍然为空的行（整组都缺失），直接删除
df = df.dropna(subset=['absorbance_450'])

# ============ 5. 验证 temperature 值 ============
invalid_temp = df[~df['temperature'].isin([25, 37])]
if len(invalid_temp) > 0:
    print(f"\n异常温度值: {invalid_temp['temperature'].unique()}")
    df = df[df['temperature'].isin([25, 37])]

# ============ 6. 保存清洗后数据 ============
df.to_csv('enzyme_activity_cleaned.csv', index=False)
print(f"\n=== 清洗完成 ===")
print(f"最终数据量: {df.shape[0]} 行")

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def detect_outliers(series, name="data"):
    """用三种方法检测异常值"""
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))

    # --- 方法 1: Z-score ---
    z_scores = np.abs((series - series.mean()) / series.std())
    outliers_z = z_scores > 3

    # --- 方法 2: IQR ---
    Q1, Q3 = series.quantile(0.25), series.quantile(0.75)
    IQR = Q3 - Q1
    lower, upper = Q1 - 1.5 * IQR, Q3 + 1.5 * IQR
    outliers_iqr = (series < lower) | (series > upper)

    # --- 可视化 ---
    axes[0].boxplot(series.dropna())
    axes[0].set_title(f'{name} - Boxplot')

    axes[1].hist(series.dropna(), bins=30, edgecolor='black')
    axes[1].axvline(lower, color='r', linestyle='--', label=f'IQR bounds')
    axes[1].axvline(upper, color='r', linestyle='--')
    axes[1].set_title(f'{name} - Histogram + IQR')
    axes[1].legend()

    axes[2].scatter(range(len(series)), series, c=outliers_z, 
                    cmap='coolwarm', alpha=0.5, s=10)
    axes[2].set_title(f'{name} - Z-score Outliers')

    plt.tight_layout()
    plt.savefig(f'outlier_detection_{name}.png', dpi=150, bbox_inches='tight')
    plt.show()

    print(f"Z-score 方法检测到 {outliers_z.sum()} 个异常值")
    print(f"IQR 方法检测到 {outliers_iqr.sum()} 个异常值")

    return outliers_z, outliers_iqr

# 使用
outliers_z, outliers_iqr = detect_outliers(df['absorbance_450'], 'Absorbance')

import pandas as pd
import numpy as np
from scipy import stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import matplotlib.pyplot as plt
import seaborn as sns

# ============ 1. 描述性统计 ============
desc = df.groupby('group')['expression'].agg(['mean', 'std', 'sem', 'count'])
print("=== 描述性统计 ===")
print(desc.round(3))
print()

# ============ 2. 正态性检验 ============
print("=== Shapiro-Wilk 正态性检验 ===")
groups = df['group'].unique()
normal_flag = True
for g in groups:
    data = df[df['group'] == g]['expression']
    stat, p = stats.shapiro(data)
    result = "正态" if p > 0.05 else "非正态"
    print(f"  {g}: W={stat:.4f}, p={p:.4f} → {result}")
    if p <= 0.05:
        normal_flag = False
print()

# ============ 3. 方差齐性检验 ============
group_data = [df[df['group'] == g]['expression'].values for g in groups]
stat, p = stats.levene(*group_data)
equal_var = p > 0.05
print(f"=== Levene 方差齐性检验 ===")
print(f"  F={stat:.4f}, p={p:.4f} → {'方差齐性' if equal_var else '方差不齐'}")
print()

# ============ 4. 组间比较 ============
if normal_flag and equal_var:
    # 满足参数检验条件 → One-way ANOVA
    stat, p = stats.f_oneway(*group_data)
    print(f"=== One-way ANOVA ===")
    print(f"  F={stat:.4f}, p={p:.4f}")
    test_name = "ANOVA"
else:
    # 不满足 → Kruskal-Wallis
    stat, p = stats.kruskal(*group_data)
    print(f"=== Kruskal-Wallis 检验 ===")
    print(f"  H={stat:.4f}, p={p:.4f}")
    test_name = "Kruskal-Wallis"

# ============ 5. 事后多重比较 ============
if p < 0.05:
    print(f"\n=== Tukey HSD 事后比较 ===")
    tukey = pairwise_tukeyhsd(df['expression'], df['group'], alpha=0.05)
    print(tukey)

import matplotlib.pyplot as plt
import matplotlib as mpl

# ============ 期刊级图表全局设置 ============
def set_publication_style():
    """设置适合期刊发表的 matplotlib 全局参数"""
    mpl.rcParams.update({
        # --- 字体 ---
        'font.family': 'Arial',       # 大多数期刊要求 Arial 或 Helvetica
        'font.size': 8,                # 基础字号 8pt
        'axes.titlesize': 9,           # 子图标题
        'axes.labelsize': 8,           # 坐标轴标签
        'xtick.labelsize': 7,          # 刻度标签
        'ytick.labelsize': 7,
        'legend.fontsize': 7,          # 图例

        # --- 线条 ---
        'axes.linewidth': 0.8,         # 坐标轴线宽
        'xtick.major.width': 0.8,      # 刻度线宽
        'ytick.major.width': 0.8,
        'xtick.major.size': 3,         # 刻度线长
        'ytick.major.size': 3,
        'lines.linewidth': 1.0,        # 数据线宽
        'lines.markersize': 4,         # 标记大小

        # --- 图像 ---
        'figure.dpi': 150,             # 屏幕显示 DPI
        'savefig.dpi': 300,            # 保存 DPI（照片用 300，线图用 600）
        'savefig.bbox': 'tight',       # 自动裁剪白边
        'savefig.pad_inches': 0.05,    # 留极小的边距

        # --- 其他 ---
        'axes.spines.top': False,      # 隐藏上边框
        'axes.spines.right': False,    # 隐藏右边框（更简洁）
        'legend.frameon': False,       # 图例无边框
        'pdf.fonttype': 42,            # PDF 中嵌入字体（期刊要求）
        'ps.fonttype': 42,
    })

# 在脚本开头调用
set_publication_style()

# Nature 单栏宽度
FIG_WIDTH_SINGLE = 89 / 25.4   # 89mm → inches
# Nature 双栏宽度
FIG_WIDTH_DOUBLE = 183 / 25.4  # 183mm → inches

import matplotlib.pyplot as plt
import numpy as np

set_publication_style()  # 使用上面定义的全局设置

# --- 数据 ---
groups = ['Control', 'Treatment A', 'Treatment B']
means = [1.0, 1.45, 1.82]
sems = [0.05, 0.06, 0.08]
# 原始数据点（用于散点叠加）
raw_data = [
    np.random.normal(1.0, 0.15, 10),
    np.random.normal(1.45, 0.20, 10),
    np.random.normal(1.82, 0.25, 10)
]

# --- Nature 配色 ---
colors = ['#BFBFBF', '#4DBBD5', '#E64B35']

# --- 绑定 ---
fig, ax = plt.subplots(figsize=(89/25.4, 60/25.4))  # 单栏宽，适中高度

# 柱状图
bars = ax.bar(range(3), means, yerr=sems, capsize=3,
              color=colors, edgecolor='black', linewidth=0.5,
              width=0.6, error_kw={'linewidth': 0.8})

# 散点叠加
for i, data in enumerate(raw_data):
    jitter = np.random.uniform(-0.15, 0.15, len(data))
    ax.scatter(i + jitter, data, color='black', s=8, alpha=0.6,
               zorder=5, linewidths=0)

# 显著性标注
def add_significance(ax, x1, x2, y, p_text, h=0.03):
    """添加显著性标注横线和星号"""
    ax.plot([x1, x1, x2, x2], [y, y+h, y+h, y], lw=0.8, c='black')
    ax.text((x1+x2)/2, y+h, p_text, ha='center', va='bottom', fontsize=7)

add_significance(ax, 0, 1, 1.7, '*')
add_significance(ax, 0, 2, 1.9, '***')

# --- 格式调整 ---
ax.set_xticks(range(3))
ax.set_xticklabels(groups)
ax.set_ylabel('Relative Expression')
ax.set_ylim(0, 2.2)

# 保存
fig.savefig('barplot.pdf')
fig.savefig('barplot.tiff', dpi=600)
plt.show()

import seaborn as sns
import pandas as pd
import numpy as np

# 使用 seaborn clustermap
g = sns.clustermap(
    df_expression,
    cmap='RdBu_r',
    center=0,
    vmin=-3, vmax=3,
    figsize=(90/25.4, 120/25.4),
    dendrogram_ratio=(0.15, 0.15),
    cbar_kws={'label': r'$\log_2$(Fold Change)', 'shrink': 0.6},
    linewidths=0.5,
    linecolor='white',
    xticklabels=True,
    yticklabels=True,
    col_colors=col_color_map,  # 列注释颜色
    method='ward',
    metric='euclidean'
)

# 调整字体大小
g.ax_heatmap.tick_params(axis='both', labelsize=6)
g.ax_heatmap.set_xlabel('')
g.ax_heatmap.set_ylabel('')

plt.savefig('heatmap.pdf', bbox_inches='tight')

import matplotlib.pyplot as plt
import numpy as np
from adjustText import adjust_text

set_publication_style()

fig, ax = plt.subplots(figsize=(89/25.4, 75/25.4))

# 分类标注
conditions = [
    (df['log2FC'] > 1) & (df['padj'] < 0.05),   # 显著上调
    (df['log2FC'] < -1) & (df['padj'] < 0.05),   # 显著下调
]
colors = ['#E64B35', '#4DBBD5', '#BFBFBF']  # 红、蓝、灰
labels = [f'Up ({conditions[0].sum()})', 
          f'Down ({conditions[1].sum()})', 
          'NS']

# 计算 -log10(padj)
df['neg_log10p'] = -np.log10(df['padj'].clip(lower=1e-300))

# 不显著的灰色点先画（在底层）
ns = ~conditions[0] & ~conditions[1]
ax.scatter(df.loc[ns, 'log2FC'], df.loc[ns, 'neg_log10p'],
           c='#BFBFBF', s=3, alpha=0.5, label=f'NS ({ns.sum()})')

# 显著上调
ax.scatter(df.loc[conditions[0], 'log2FC'], 
           df.loc[conditions[0], 'neg_log10p'],
           c='#E64B35', s=5, alpha=0.7, label=labels[0])

# 显著下调
ax.scatter(df.loc[conditions[1], 'log2FC'], 
           df.loc[conditions[1], 'neg_log10p'],
           c='#4DBBD5', s=5, alpha=0.7, label=labels[1])

# 阈值线
ax.axhline(-np.log10(0.05), ls='--', lw=0.5, c='grey')
ax.axvline(-1, ls='--', lw=0.5, c='grey')
ax.axvline(1, ls='--', lw=0.5, c='grey')

# 标注 top 基因
sig = df[conditions[0] | conditions[1]].nlargest(10, 'neg_log10p')
texts = []
for _, row in sig.iterrows():
    texts.append(ax.text(row['log2FC'], row['neg_log10p'], row['gene'],
                         fontsize=5, style='italic'))
adjust_text(texts, arrowprops=dict(arrowstyle='-', color='grey', lw=0.5))

ax.set_xlabel(r'$\log_2$(Fold Change)')
ax.set_ylabel(r'$-\log_{10}$(adjusted p-value)')
ax.legend(fontsize=6, loc='upper right', markerscale=1.5)

fig.savefig('volcano_plot.pdf')
plt.show()

import seaborn as sns
import matplotlib.pyplot as plt

set_publication_style()
fig, ax = plt.subplots(figsize=(89/25.4, 65/25.4))

# 小提琴图
parts = ax.violinplot([data_ctrl, data_trtA, data_trtB], 
                       positions=[0, 1, 2], showextrema=False)
for i, pc in enumerate(parts['bodies']):
    pc.set_facecolor(colors[i])
    pc.set_alpha(0.3)

# 箱线图叠加
bp = ax.boxplot([data_ctrl, data_trtA, data_trtB],
                positions=[0, 1, 2], widths=0.15,
                patch_artist=True, showfliers=False,
                medianprops=dict(color='black', linewidth=1),
                boxprops=dict(facecolor='white', linewidth=0.8))

# 散点叠加
for i, data in enumerate([data_ctrl, data_trtA, data_trtB]):
    jitter = np.random.uniform(-0.05, 0.05, len(data))
    ax.scatter(i + jitter, data, color=colors[i], s=10, alpha=0.6,
               edgecolors='none', zorder=5)

ax.set_xticks([0, 1, 2])
ax.set_xticklabels(['Control', 'Treatment A', 'Treatment B'])
ax.set_ylabel('Expression Level')

fig.savefig('violin_box.pdf')

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

set_publication_style()

fig = plt.figure(figsize=(183/25.4, 120/25.4))

# 创建不等宽网格
gs = gridspec.GridSpec(2, 3, 
                       height_ratios=[1, 1],
                       width_ratios=[1, 1, 1],
                       hspace=0.35, wspace=0.35)

ax_a = fig.add_subplot(gs[0, 0])     # A: 第一行第一列
ax_b = fig.add_subplot(gs[0, 1])     # B: 第一行第二列
ax_c = fig.add_subplot(gs[0, 2])     # C: 第一行第三列
ax_d = fig.add_subplot(gs[1, 0:2])   # D: 第二行跨前两列
ax_e = fig.add_subplot(gs[1, 2])     # E: 第二行第三列

# 面板标签
for ax, label in zip([ax_a, ax_b, ax_c, ax_d, ax_e], 
                      ['a', 'b', 'c', 'd', 'e']):
    ax.text(-0.15, 1.10, label, transform=ax.transAxes,
            fontsize=10, fontweight='bold', va='top')

# [在每个 ax 上绑定对应内容...]

fig.savefig('Figure1.pdf')
fig.savefig('Figure1.tiff', dpi=600)

#!/usr/bin/env python3
"""
蛋白质表达量分析与绘图脚本
分析三组 Western Blot 定量数据，生成 Nature 风格柱状图
"""

import pandas as pd
import numpy as np
from scipy import stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import matplotlib.pyplot as plt
import matplotlib as mpl
from itertools import combinations

# ================================================================
# 全局设置
# ================================================================
mpl.rcParams.update({
    'font.family': 'Arial',
    'font.size': 8,
    'axes.labelsize': 8,
    'axes.titlesize': 9,
    'xtick.labelsize': 7,
    'ytick.labelsize': 7,
    'legend.fontsize': 7,
    'axes.linewidth': 0.8,
    'xtick.major.width': 0.8,
    'ytick.major.width': 0.8,
    'axes.spines.top': False,
    'axes.spines.right': False,
    'figure.dpi': 150,
    'savefig.dpi': 300,
    'savefig.bbox': 'tight',
    'pdf.fonttype': 42,
    'ps.fonttype': 42,
})

# Nature 风格配色
COLORS = {
    'Control': '#BFBFBF',
    'Drug_A': '#4DBBD5',
    'Drug_B': '#E64B35'
}

# ================================================================
# 1. 读入数据
# ================================================================
df = pd.read_csv('protein_expression.csv')
print("=" * 50)
print("1. 数据概况")
print("=" * 50)
print(df.head())
print(f"\n样本量: {df.shape[0]}")
print(f"各组样本量:\n{df['group'].value_counts()}")

# ================================================================
# 2. 描述性统计
# ================================================================
print("\n" + "=" * 50)
print("2. 描述性统计")
print("=" * 50)
desc = df.groupby('group')['expression'].agg(
    ['count', 'mean', 'std', 'sem', 'median']
).round(4)
print(desc)

# ================================================================
# 3. 正态性检验
# ================================================================
print("\n" + "=" * 50)
print("3. Shapiro-Wilk 正态性检验")
print("=" * 50)
groups = ['Control', 'Drug_A', 'Drug_B']
group_data = {}
all_normal = True
for g in groups:
    data = df[df['group'] == g]['expression'].values
    group_data[g] = data
    stat, p = stats.shapiro(data)
    normal = "✓ 正态" if p > 0.05 else "✗ 非正态"
    if p <= 0.05:
        all_normal = False
    print(f"  {g}: W={stat:.4f}, p={p:.4f} {normal}")

# ================================================================
# 4. 方差齐性检验
# ================================================================
print("\n" + "=" * 50)
print("4. Levene 方差齐性检验")
print("=" * 50)
stat, p = stats.levene(*group_data.values())
equal_var = p > 0.05
print(f"  F={stat:.4f}, p={p:.4f} {'✓ 方差齐性' if equal_var else '✗ 方差不齐'}")

# ================================================================
# 5. 组间比较
# ================================================================
print("\n" + "=" * 50)
print("5. 组间比较")
print("=" * 50)

if all_normal and equal_var:
    stat, p_overall = stats.f_oneway(*group_data.values())
    print(f"[One-way ANOVA]  F={stat:.4f}, p={p_overall:.6f}")
    test_method = "One-way ANOVA"
else:
    stat, p_overall = stats.kruskal(*group_data.values())
    print(f"[Kruskal-Wallis]  H={stat:.4f}, p={p_overall:.6f}")
    test_method = "Kruskal-Wallis"

# ================================================================
# 6. 事后多重比较
# ================================================================
pairwise_results = {}
if p_overall < 0.05:
    print("\n" + "=" * 50)
    print("6. 事后比较 (Tukey HSD)")
    print("=" * 50)
    tukey = pairwise_tukeyhsd(df['expression'], df['group'], alpha=0.05)
    print(tukey)

    # 提取 p 值
    for i, (g1, g2) in enumerate(combinations(groups, 2)):
        p_val = tukey.pvalues[i]
        pairwise_results[(g1, g2)] = p_val

# ================================================================
# 7. 效应量 (Cohen's d)
# ================================================================
print("\n" + "=" * 50)
print("7. 效应量 (Cohen's d)")
print("=" * 50)

def cohens_d(a, b):
    na, nb = len(a), len(b)
    pooled_std = np.sqrt(((na-1)*np.std(a,ddof=1)**2 + 
                           (nb-1)*np.std(b,ddof=1)**2) / (na+nb-2))
    return (np.mean(a) - np.mean(b)) / pooled_std

for g1, g2 in combinations(groups, 2):
    d = cohens_d(group_data[g1], group_data[g2])
    size = "小" if abs(d) < 0.5 else ("中" if abs(d) < 0.8 else "大")
    print(f"  {g1} vs {g2}: d={d:.3f} ({size}效应)")

# ================================================================
# 8. 绑定
# ================================================================
print("\n" + "=" * 50)
print("8. 生成图表")
print("=" * 50)

fig, ax = plt.subplots(figsize=(89/25.4, 65/25.4))

# 柱状图
means = [group_data[g].mean() for g in groups]
sems = [stats.sem(group_data[g]) for g in groups]
x_pos = np.arange(len(groups))

bars = ax.bar(x_pos, means, yerr=sems, capsize=3,
              color=[COLORS[g] for g in groups],
              edgecolor='black', linewidth=0.5,
              width=0.6, error_kw={'linewidth': 0.8})

# 散点叠加
for i, g in enumerate(groups):
    data = group_data[g]
    jitter = np.random.default_rng(42).uniform(-0.12, 0.12, len(data))
    ax.scatter(i + jitter, data, color='black', s=12, alpha=0.6,
               zorder=5, linewidths=0.3, edgecolors='white')

# 显著性标注
def add_sig(ax, x1, x2, y, p_val):
    if p_val >= 0.05:
        return  # 不显著不标注
    text = '***' if p_val < 0.001 else ('**' if p_val < 0.01 else '*')
    h = 0.02 * (ax.get_ylim()[1] - ax.get_ylim()[0])
    ax.plot([x1, x1, x2, x2], [y, y+h, y+h, y], lw=0.8, c='black')
    ax.text((x1+x2)/2, y+h, text, ha='center', va='bottom', fontsize=7)

# 动态计算标注位置
y_max = max(max(group_data[g]) for g in groups)
sig_y = y_max * 1.05
for (g1, g2), p_val in pairwise_results.items():
    i1, i2 = groups.index(g1), groups.index(g2)
    add_sig(ax, i1, i2, sig_y, p_val)
    sig_y += y_max * 0.12  # 逐级升高避免重叠

# 格式
ax.set_xticks(x_pos)
ax.set_xticklabels([g.replace('_', ' ') for g in groups])
ax.set_ylabel('Relative Protein Expression')
ax.set_ylim(0, sig_y + y_max * 0.1)

# 保存
fig.savefig('Figure_protein_expression.pdf')
fig.savefig('Figure_protein_expression.tiff', dpi=600)
print("  ✓ 已保存 Figure_protein_expression.pdf")
print("  ✓ 已保存 Figure_protein_expression.tiff (600 DPI)")
plt.show()

# ================================================================
# 9. 论文用统计描述
# ================================================================
print("\n" + "=" * 50)
print("9. 论文用统计描述（英文）")
print("=" * 50)

mean_ctrl = desc.loc['Control', 'mean']
sem_ctrl = desc.loc['Control', 'sem']
mean_a = desc.loc['Drug_A', 'mean']
sem_a = desc.loc['Drug_A', 'sem']
mean_b = desc.loc['Drug_B', 'mean']
sem_b = desc.loc['Drug_B', 'sem']

p_str = "p<0.001" if p_overall < 0.001 else f"p={p_overall:.3f}"

print(f"""
Protein expression levels were compared across three groups using 
{test_method}. Drug A treatment significantly increased target protein 
expression compared to the control group ({mean_a:.2f} ± {sem_a:.2f} vs 
{mean_ctrl:.2f} ± {sem_ctrl:.2f}, mean ± SEM; Tukey HSD, 
p={'<0.001' if pairwise_results.get(('Control','Drug_A'),1)<0.001 
   else f"={pairwise_results.get(('Control','Drug_A'),1):.3f}"}). 
Drug B treatment showed an even greater increase ({mean_b:.2f} ± {sem_b:.2f}; 
p={'<0.001' if pairwise_results.get(('Control','Drug_B'),1)<0.001 
   else f"={pairwise_results.get(('Control','Drug_B'),1):.3f}"} vs control). 
Overall {test_method}: {p_str}. All values are presented as mean ± SEM 
(n=6 per group).
""")