第4章 数据探索与可视化：用箱线图和密度图揪出批次差异

数据加载进来之后，别急着跑模型。我个人习惯先拿眼睛“看一看”数据长什么样。这一步叫数据探索，说白了就是让数据自己说话。批次效应检测的第一步，往往就是靠视觉发现端倪。嗯，今天咱们就用 matplotlib 和 seaborn 这两个工具，画出箱线图和密度图，直观地看看不同批次之间到底有没有差异。

4.1 准备数据：模拟一个有批次效应的表达矩阵

为了演示方便，我们先造一批数据。我假设你有3个批次（Batch1、Batch2、Batch3），每个批次里各有10个样本。每个样本测了5个基因的表达值。我会让Batch1和Batch2比较像，Batch3则人为加上一个偏移——这就是我们想检测的批次效应。

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

# 设置随机种子，保证结果可复现
np.random.seed(42)

# 定义批次和样本数
batches = ['Batch1', 'Batch2', 'Batch3']
samples_per_batch = 10
genes = ['Gene_A', 'Gene_B', 'Gene_C', 'Gene_D', 'Gene_E']

# 生成表达矩阵
data = []
batch_labels = []

for i, batch in enumerate(batches):
    # 每个批次的基础表达水平，Batch3 整体上调 2 个单位
    base_mean = 5.0 if i < 2 else 7.0
    for _ in range(samples_per_batch):
        # 每个样本的基因表达值，加入随机噪声
        expression = np.random.normal(loc=base_mean, scale=1.0, size=len(genes))
        data.append(expression)
        batch_labels.append(batch)

# 构建 DataFrame
df = pd.DataFrame(data, columns=genes)
df['Batch'] = batch_labels

print(df.head())

4.2 箱线图：一眼看出批次间的中心偏移

箱线图能告诉我们每个批次里基因表达值的中位数、四分位数范围，还有异常值。我最喜欢用它来找“整体偏移”。

如果你想把所有基因合在一起看，直接画一个图就够了。但如果你想看看每个基因在不同批次的表现，那就得用“分面”或者分组。

4.2.1 所有基因合在一起画

plt.figure(figsize=(8, 6))
sns.boxplot(x='Batch', y='value', data=df.melt(id_vars='Batch', var_name='Gene', value_name='value'))
plt.title('所有基因表达值的批次间箱线图')
plt.show()

你会看到：Batch1和Batch2的箱子高度和中位数差不多，但Batch3整个箱子都往上提了一截。这就是批次效应的典型表现。我当年做RNA-seq项目时，第一次看到这种图，心里就有底了——后续的校正方向基本确定。

4.2.2 分基因观察

不过，合在一起画有时会掩盖一些基因特异性的问题。比如某个基因在Batch1里表达特别高，但平均值被拖平了。这时候，我建议你分基因画箱线图。

plt.figure(figsize=(12, 6))
sns.boxplot(x='Gene', y='value', hue='Batch', data=df.melt(id_vars='Batch', var_name='Gene', value_name='value'))
plt.title('按基因拆分，观察各批次表达差异')
plt.show()

这里 hue='Batch' 会把每个基因的箱子按批次分开。如果某个基因上三个批次的箱子高度差异特别大，那就值得重点排查。

4.3 密度图：看清数据分布的“形状”

箱线图只能看中位数和范围，但分布的形状它给不出来。比如说，两个批次的中位数一样，但一个分布很宽、另一个很窄，这其实也是批次效应的表现。密度图正好补上这个缺口。

4.3.1 核密度估计图（KDE）

plt.figure(figsize=(8, 6))
sns.kdeplot(data=df, x='Gene_A', hue='Batch', fill=True, common_norm=False, alpha=0.5)
plt.title('Gene_A 在不同批次的密度分布')
plt.show()

这段代码会画出Gene_A在每个批次里的概率密度曲线。如果三条曲线基本重合，说明这个基因不受批次影响。反之，如果Batch3的曲线明显整体右移，那就坐实了批次效应的存在。

我记得有一次，客户死活不承认自己的数据有批次问题。我什么话都没说，给他画了一张密度图——Batch1和Batch2的曲线几乎重叠，Batch3的曲线孤零零地飘在右边。他看完就沉默了。数据往往比语言更有说服力。

4.3.2 结合直方图的联合展示

有时候密度图太光滑，你可能会错过一些数据的细节。我个人喜欢把直方图和密度图叠在一起用。

plt.figure(figsize=(8, 6))
sns.histplot(data=df, x='Gene_B', hue='Batch', kde=True, stat='density', alpha=0.4)
plt.title('Gene_B 直方图 + 核密度估计')
plt.show()

这样做的好处是：直方图告诉你原始数据的“柱子”分布，KDE曲线则给你一个平滑的趋势。两者互相验证，不容易被假象误导。

4.4 多基因表达谱的热图+聚类

光看单个基因还不够。有时候我们需要看一个批次内所有基因的整体表达模式。热图加上层次聚类，就是一个很好的组合。

# 构造一个按批次排序的数据矩阵
heatmap_data = df.set_index('Batch').T  # 转置成基因为行、样本为列
sns.clustermap(heatmap_data, col_colors=df['Batch'].map({'Batch1': 'blue', 'Batch2': 'green', 'Batch3': 'red'}),
               cmap='vlag', linewidths=0.5, figsize=(10, 8))
plt.suptitle('表达数据热图 + 层次聚类', y=1.02)
plt.show()

如果你看到样本的聚类树把Batch1和Batch2聚在一起，Batch3单独聚成一簇，那就说明：批次效应比你想象的要大。有一次我做的项目里，热图上的聚类结果甚至跟生物学分组完全反了——那是我第一次深刻意识到批次效应校正的必要性。

4.5 本章知识结构图

下面的SVG图帮你理清本章的数据探索流程。我们从原始表达矩阵出发，分别进行箱线图分析和密度图分析，再整合成热图+聚类，最终判断是否存在批次效应。

4.6 一些避坑提醒

• 样本量不均衡：如果一个批次有50个样本，另一个批次只有3个，箱线图里的箱子大小会差很多，视觉上容易产生误导。我曾经吃过这个亏，最后用“小提琴图”替代了箱线图。
• 不要只盯着一个基因看：批次效应往往是一个全局现象。如果你只看一个基因，可能会误以为是生物学变异。建议至少看5-10个基因的表达分布，或者直接用所有基因的平均值作图。
• 密度图的带宽参数：sns.kdeplot 里的 bw_method 参数会影响曲线的平滑程度。默认值在大多数情况下够用，但如果你的样本量很小（比如每个批次只有10个样本），可以适当调大带宽避免过度拟合。
• 注意图的可读性：当基因数量超过10个时，分面箱线图会变得很拥挤。我在一个750个基因的数据集上试过，结果X轴的标签完全搅在一起。这时候我一般会先筛选出差异最大的几个基因，或者用热图替代。