第3章 数据基础：高通量测序数据格式与模拟批次效应数据构建
做批次效应可视化之前，你得先搞清楚数据长什么样。说白了，高通量测序出来的一堆 FASTQ 文件，经过一系列处理，最终会落到两种最常见的表格格式里——OTU 表和丰度矩阵。我刚开始接触这批数据时，差点被「OTU」「ASV」「feature table」这些名词绕晕。后来才明白，它们本质上就是同一个东西的不同叫法。
3.1 OTU 表：最原始的计数派
OTU（Operational Taxonomic Unit）表，行是物种分类单元（通常按 97% 相似度聚类），列是样本，交叉格子里是序列条数。你看到的就是一个整数矩阵，每个数字代表「这个 OTU 在这个样本里被检测到了多少次」。举个例子：
| OTU ID | Sample_A | Sample_B | Sample_C |
|---|---|---|---|
| OTU_1 | 1250 | 872 | 2034 |
| OTU_2 | 346 | 128 | 576 |
| OTU_3 | 89 | 245 | 12 |
嗯，这看起来挺清爽对吧？但实际项目里，OTU 表动辄上千行、几百列，数字分布极不均匀——有的 OTU 只在少数样本里出现，有的则是「核心物种」遍布所有样本。我曾在处理一批婴儿肠道样本时，发现一个 OTU 从 0 直接蹦到 50000，差点以为是程序 bug，后来一查是某个极端样本的稀释深度没做好。所以拿到 OTU 表，第一件事就是统计列和（library size），看看有没有量级差异巨大的样本。
df.sum() 扫一遍列总和，如果最大值/最小值 >10 倍，就得考虑 rarefaction 或 CSS 归一化，否则后续 PCoA/UMAP 会被深度最大的样本带偏。
3.2 丰度矩阵：归一化后的战场
OTU 表是绝对的计数，但不同样本的测序深度不一样，直接比较不合理。所以我们会把计数转成相对丰度（每个 OTU 占该样本总读数的百分比），或者用 TMM、CLR 等方法归一化。得到的矩阵就是丰度矩阵，数值通常在 0~1 之间（或经过 log 变换后可为负值）。
我个人的偏好是：做 PCoA 前尽量用 CLR 变换（Centered Log-Ratio），因为它是基于成分数据的，特别适合微生物组这种「和恒为 1」的数据结构。而 UMAP 对归一化方式没那么敏感，但建议至少做一下总读数归一化，否则极端值会让 UMAP 变成「一个样本一个孤岛」。
3.3 模拟批次效应数据：为什么需要自己造？
很多时候你拿不到真实的批次标签——实验记录丢了、样本来自不同项目、或者你只是想练练手。怎么办？自己造。模拟批次效应数据的好处是：你知道「真批次」是什么，可以清楚评估你的降维方法能否分离出批次。
怎么模拟？说白了就是在原本无批次的丰度矩阵上，给不同「批次组」的样本施加一个偏移量。这个偏移量可以加在 OTU 层级（所有物种被放大/缩小），也可以只加在部分物种上（模拟技术偏好）。下面我给出一个 Python 示例——用 numpy 和 pandas 构建模拟数据：
import numpy as np
import pandas as pd
# 设置随机种子，保证可复现
np.random.seed(42)
# 模拟 100 个 OTU × 30 个样本的原始计数
n_otu = 100
n_sample = 30
base_abund = np.random.negative_binomial(n=5, p=0.3, size=(n_otu, n_sample))
# 给前 10 个样本打上「批次1」标签，后 20 个样本打上「批次2」
batch_labels = np.array(['batch1']*10 + ['batch2']*20)
# 批次效应：给批次2的样本整体乘以一个系数 1.5
effect = np.ones((n_otu, n_sample))
effect[:, 10:] = 1.5
# 添加随机噪声（log-normal 噪声更真实）
noise = np.random.lognormal(mean=0, sigma=0.2, size=(n_otu, n_sample))
# 最终的模拟计数
sim_counts = (base_abund * effect * noise).astype(int)
# 构建 OTU 表
otu_ids = [f'OTU_{i+1}' for i in range(n_otu)]
sample_ids = [f'Sample_{j+1}' for j in range(n_sample)]
otu_table = pd.DataFrame(sim_counts, index=otu_ids, columns=sample_ids)
# 保存为 tsv
otu_table.to_csv('simulated_otu_table.tsv', sep='\t')
print(otu_table.iloc[:5, :5])上面这段代码的核心逻辑：先生成无批次基准数据，再给第二批次乘以 1.5 倍，加上对数正态噪声。这种「乘法偏移」在真实测序中很常见——同一个批次扩增效率不同，导致所有物种都按比例放大或缩小。你仔细看 effect 矩阵的构造方式，就知道批次效应只作用在列上，行（OTU）之间没有差别，这叫全局批次效应。
effect 改成稀疏矩阵，只随机选一些行做偏移。我在实际项目里见过更复杂的：不同批次里某些属的丰度直接翻倍或减半，而且只在特定 pH 条件下……
为什么这里用负二项分布模拟原始计数？因为微生物组数据过离散（方差 > 均值），负二项比泊松更合适。嗯，这个细节可能有点超纲，但你记住一句话：模拟数据必须模仿真实数据的统计特征，否则你的方法论验证是无效的。
3.4 整个流程：从原始测序到模拟批次数据
下面这张图把我们从 raw data 到可操作数据矩阵的步骤串起来了。核心就两步：先得到干净的丰度矩阵，再人为加入批次效应。
图里有一条虚线分支——你可以从 OTU 表直接进入模拟批次效应模块，也可以先归一化再加批次效应。我个人建议先归一化再加批次，这样能模拟出「不同批次使用了不同实验操作导致的技术偏差」这一真实场景；如果你直接对原始计数加偏移，测序深度的差异会和批次效应混在一起，实操时很难分离。
3.5 小练习：给你的模拟数据添个「真批次」标签
上面代码里我已经用 batch_labels 存好了批次信息。但如果你用的是真实数据，怎么知道哪些样本属于同一批次？这往往需要去查实验记录——样本编号、处理日期、提取试剂盒批次。我建议你在构建模拟数据时，顺便存一个批次注释文件，像这样：
sample_id,batch
Sample_1,batch1
Sample_2,batch1
...
Sample_11,batch2
Sample_12,batch2后面用 PCoA 或 UMAP 着色时，这个文件就是你的金标准。我曾见过一个团队，辛辛苦苦降维完发现聚类结果和批次标签高度吻合，但因为没有提前记录批次，愣是浪费了两周去猜「是不是样品类型导致的差异」——你猜怎么着？最后发现就是同一个批次的试剂污染。所以记录源信息，从模拟数据开始就要养成习惯。