数据清洗与预处理：把脏数据变成分析就绪态
数据清洗这步，说实在的，占了我实战项目里差不多一半的时间。你想想看，拿到手的丰度表，十个样本里三四个有缺失，低丰度OTU一抓一大把，不做清洗直接跑差异性分析，结果根本没法看。我刚开始做微生物组项目的时候，就吃过这个亏——把原始数据直接塞进LEfSe，出来的结果花里胡哨，可老板一质疑就全垮了。从那以后，我养成了一个习惯：先花一整天清洗，再花一下午分析。
1. 缺失值处理——别让NA毁了你的统计
缺失值在丰度数据里太常见了。测序深度不够、PCR扩增失败、样本标签搞混……原因五花八门。我自己最头疼的是那种“间歇性缺失”：同一处理组的几个重复里，有的OTU有值，有的没有。这时候直接删掉整个OTU，信息损失太大；不删吧，统计检验又报错。
常见的处理策略就三种：
- 删除：样本缺失率超过50%，直接删样本；OTU在所有样本中缺失率超过80%，删OTU。这招简单粗暴，但我建议别一上来就用。先看看缺失是不是随机出现的。
- 填补：用均值、中位数、最小值（比如检测限的1/2）或kNN预测。我偏爱用
zCompositions包的cmultRepl，它能利用成分数据的比例关系做贝叶斯替换，效果比简单均值好得多。 - 模型估计：把缺失当作待估参数，用MCMC或极大似然来补。计算量大，但精度最高。小数据集慎用，容易过拟合。
下面是我常用的填补代码：
library(zCompositions)
# 假设你的丰度表是 otu_table，行是OTU，列是样本
# cmultRepl 会自动识别0值，也可以指定处理NA
otu_filled <- cmultRepl(otu_table, label = 0, method = "CZM")
# 如果是NA，可以先用0占位再填补：
otu_table[is.na(otu_table)] <- 0
otu_filled <- cmultRepl(otu_table, label = 0, method = "CZM")2. 低丰度过滤——把噪声踢出去
低丰度OTU（比如平均相对丰度 < 0.01%）多数是测序错误或罕见菌。留着它们，不仅增加多重比较的校正负担，还会让稀疏性加剧，很多多元统计方法（如CCA）都对零过多的数据不友好。我个人的经验是：先画一个丰度分布直方图，看看阈值大概在哪。
常见的过滤策略：
- 基于丰度值：删除所有样本中平均丰度低于某个阈值（如1e-4）的OTU。
- 基于出现频率：删除在少于X%样本中出现的OTU（比如只在5%以下样本中出现的去掉）。
- 两者结合：比如“平均丰度 > 1e-5 且出现频率 > 10%”。
给个代码片段：
# 假设 otu 是矩阵，行OTU，列样本
min_abundance <- 1e-4
min_prevalence <- 0.1 # 至少出现在10%的样本中
# 计算每个OTU的平均丰度
mean_abund <- rowMeans(otu)
# 计算每个OTU的出现频率（非零样本比例）
prevalence <- rowSums(otu > 0) / ncol(otu)
# 过滤
keep <- mean_abund > min_abundance & prevalence > min_prevalence
otu_filtered <- otu[keep, ]3. 标准化方法——让样本之间能比
标准化这一步，说白了就是把不同测序深度的样本拉到同一把尺子上。如果直接用原始counts做差异分析，测序深度高的样本天然会有更多高丰度OTU，这可不是生物学差异。
下面我挑三种（相对丰度、CSS、TMM）聊聊，这也是我项目里最常用的。
| 方法 | 核心思想 | 适用场景 | R包/函数 |
|---|---|---|---|
| 相对丰度 (RA) | 每个OTU count / 总count，得到比例 | 初步探索、可视化、非参数检验 | prop.table() / decostand(RA) |
| CSS (Cumulative Sum Scaling) | 用上分位数（如75%）归一化，抵抗离群样本的影响 | 有离群大总count的样本、准备做距离计算 | metagenomeSeq::cumNorm() |
| TMM (Trimmed Mean of M-values) | 基于加权对数的截尾均值，常用在差异丰度分析前 | edgeR差异分析、样本间分散度差异大 | edgeR::calcNormFactors() |
RA最简单，但有个大坑：总count的变动会将扰动传导给所有OTU，所以我只在探索性分析时用RA。CSS是Paulson在metagenomeSeq里提出的，它只用了前75%分位的count和来缩放，能有效避免一两个超高丰度OTU带来的偏差。我有个项目，样本总count从1万到10万不等，做RA后β多样性PCoA第一轴全是总count差异；换成CSS，立马洗出生物学信号。
TMM的话，我一般在做差异丰度分析时才请它出场。它通过计算样本间M值的截尾均值来得到缩放因子，比RA更稳健。不过注意，TMM假设大多数OTU是非差异的，如果你的数据里有一半以上的OTU都是差异丰度的，那TMM反而会失效。
代码示例：
# ---- 相对丰度 ----
ra <- prop.table(as.matrix(otu_filtered), 2) # 按列归一化
# ---- CSS 标准化 ----
library(metagenomeSeq)
# 创建 MRexperiment 对象
obj <- newMRexperiment(otu_filtered)
# 计算归一化因子（默认使用上分位数0.75）
obj <- cumNorm(obj, p = cumNormStat(obj))
# 提取标准化后的counts（不是相对丰度，是缩放后的counts）
css_counts <- cumNormMat(obj)
# ---- TMM 标准化 ----
library(edgeR)
y <- DGEList(counts = otu_filtered)
y <- calcNormFactors(y, method = "TMM")
# 得到的缩放因子存于 y$samples$norm.factors
# 若需要标准化后的logCPM，可以用 cpm(y, log = TRUE)- 缺失值处理 → 用
cmultRepl填补（成分数据）或kNN（非成分数据）。 - 低丰度过滤 → 平均丰度 > 1e-4 & 出现频率 > 10%。
- 标准化 → 如果后续要做差异分析（edgeR/DESeq2），用TMM；如果只是β多样性或聚类，用CSS；探索性可视化用RA。
最后，我想强调一点：标准化方法的选择不是一成不变的。我遇到过不同标准化方法导致截然不同的排序图——同一批数据，RA第一轴解释16%，TMM第一轴解释9%。你猜哪个更可靠？答案是：得结合你的生物学假设去验证，比如看分组是否在轴上有明显分离。标准化只是工具，别让它替你做决定。