差异分析 | 当ComBat-seq校正后的Count数据仍非整数或含负值，如何精准进行下游分析？

1. 问题来了：ComBat-seq校正后，数据怎么还是“怪怪的”？

相信很多刚开始用ComBat-seq处理RNA-seq数据的小伙伴，都抱着一个美好的期待：把原始计数矩阵扔进去，它就能吐出来一份“干干净净”、全是整数的数据，然后直接喂给DESeq2或者edgeR做差异分析，一气呵成。我之前也是这么想的，直到在实际项目里踩了几个不大不小的坑。

有一次，我手头有一个整合了三个不同测序平台的数据集，批次效应非常明显。我兴冲冲地跑完了ComBat-seq，满心欢喜地检查结果，却傻眼了：调整后的数据矩阵里，居然零零散散地出现了一些负值，而且很多数值带着长长的小数点，根本不是整数。我当时的第一反应是：“是不是我代码写错了？参数调错了？” 反复检查了好几遍，流程完全按照官方文档来的。这下麻烦了，因为DESeq2和edgeR这些差异分析的金标准工具，明确要求输入的是原始计数数据，也就是非负整数。出现负值或者小数，它们要么直接报错罢工，要么给你一个警告，然后分析结果的可靠性就大打折扣了。

这种情况其实并不少见，尤其是在数据本身比较“脏”或者批次效应非常复杂的时候。ComBat-seq虽然是为计数数据设计的，但其核心的负二项回归模型在拟合和调整过程中，为了达到最佳的校正效果，可能会进行一些连续的数值运算。理论上，模型会尽力将结果向整数靠拢，但实际计算中，四舍五入、模型拟合的残差、极端值的处理等因素，都可能导致输出不再是完美的整数，甚至产生理论上不应出现的负值。这并不意味着ComBat-seq失败了，而是告诉我们，“校正完成”不等于“数据就绪”，中间还有一个重要的数据校验与后处理环节。

所以，如果你也遇到了同样的问题，先别慌，这很可能不是你的操作失误。接下来，我们就得像个侦探一样，一步步排查问题，并找到让数据“重回正轨”的方法，确保下游分析稳稳当当。

2. 诊断第一步：你的数据到底出了什么问题？

在急着动手修复之前，我们得先搞清楚数据“怪”在哪里，程度有多深。盲目处理可能会引入新的偏差。这里我分享几个我常用的诊断步骤，你可以像做体检一样，给校正后的数据来一套全面检查。

### 2.1 检查异常值的分布与比例

首先，最直观的就是看看负值和小数到底有多少。你可以写一段简单的R代码来快速统计：

# 假设 adjusted_counts 是 ComBat-seq 校正后的矩阵
# 检查负值
num_negative <- sum(adjusted_counts < 0, na.rm = TRUE)
total_values <- length(adjusted_counts)
prop_negative <- num_negative / total_values

cat(sprintf("负值数量：%d\n", num_negative))
cat(sprintf("负值比例：%.4f%%\n", prop_negative * 100))

# 检查非整数（这里定义绝对值大于1e-10的小数部分为非整数）
is_non_integer <- abs(adjusted_counts - round(adjusted_counts)) > 1e-10
num_non_int <- sum(is_non_integer, na.rm = TRUE)
prop_non_int <- num_non_int / total_values

cat(sprintf("非整数值数量：%d\n", num_non_int))
cat(sprintf("非整数值比例：%.4f%%\n", prop_non_int * 100))

跑完这段代码，你心里就有数了。如果负值比例低于0.1%，非整数比例也不高，那问题可能比较轻微，或许是一些极端基因或样本导致的。如果比例很高，比如超过1%，那就需要警惕了，可能意味着批次效应非常剧烈，或者数据中存在其他未考虑的混杂因素，ComBat-seq的模型假设在这种情况下受到了挑战。

### 2.2 追溯问题样本与基因

光知道有多少问题不够，我们还得知道是“谁”出了问题。是集中在某几个样本里，还是遍布所有样本？是某些低表达基因特有，还是高表达基因也未能幸免？

# 找出包含负值的基因（行）和样本（列）
negative_gene_indices <- which(rowSums(adjusted_counts < 0) > 0)
negative_sample_indices <- which(colSums(adjusted_counts < 0) > 0)

cat("包含负值的基因行索引：\n")
print(head(negative_gene_indices))
cat("包含负值的样本列索引：\n")
print(negative_sample_indices)

# 可以进一步查看这些基因的表达量分布（原始数据）
if(length(negative_gene_indices) > 0){
  problematic_genes_raw <- count_matrix[negative_gene_indices, ]
  summary_row_means <- rowMeans(problematic_genes_raw)
  cat("问题基因（校正前）的平均表达量摘要：\n")
  print(summary(summary_row_means))
}

这个分析非常关键。如果问题只出现在极少数超低表达（比如在所有样本中原始计数都为0或1）的基因上，那么这些基因很可能本身噪音就极大，信息量有限，在后续分析中直接过滤掉它们可能是最简洁高效的做法。如果发现是某个特定批次的样本校正后普遍出现异常，那可能需要回头审视批次信息的定义是否正确，或者考虑这个批次的数据质量是否本身就有问题。

### 2.3 可视化审视数据分布

数字是冰冷的，图形是直观的。在修复前后，我都强烈建议画几个图对比一下。

library(ggplot2)
library(reshape2)
# 将校正后的数据“融化”，便于绘图
df_melted <- melt(adjusted_counts)
colnames(df_melted) <- c("Gene", "Sample", "Adjusted_Count")

# 绘制校正后值的分布直方图（聚焦于小范围）
p1 <- ggplot(subset(df_melted, Adjusted_Count < 10 & Adjusted_Count > -5), aes(x=Adjusted_Count)) +
  geom_histogram(bins=50, fill="steelblue", alpha=0.7) +
  geom_vline(xintercept = 0, color="red", linetype="dashed") +
  labs(title="校正后计数分布（-5到10区间）", x="校正后计数", y="频数") +
  theme_minimal()

# 绘制样本间相关性热图（校正后）
cor_matrix <- cor(adjusted_counts, method="spearman")
library(pheatmap)
pheatmap(cor_matrix,
         main="样本间相关性热图（ComBat-seq校正后）",
         clustering_distance_rows="euclidean",
         clustering_distance_cols="euclidean")

直方图能让你一眼看出负值和小数聚集在哪个区间。热图则可以帮助你判断批次效应是否被有效移除——理想情况下，同一批次内的样本高相关性“块状结构”应该消失，样本聚类应更多由生物学分组驱动，而不是批次。如果校正后热图依然显示清晰的批次区块，那说明校正可能不充分，需要调整模型参数或考虑其他方法。

3. 实战修复：让数据“规整”起来的几种策略

诊断清楚后，就到了最关键的一步：怎么修？这里没有一刀切的“最佳”方案，需要根据问题的严重程度和下游分析工具的要求来权衡选择。我把自己用过、试过的几种策略整理了一下，你可以按顺序尝试。

### 3.1 策略一：温和的数值修整（四舍五入与截断）

这是最直接、最快速的方法，适用于异常值比例很低（比如<0.5%）的情况。核心思想是，把这些微小的偏差给“抹平”。

• 四舍五入取整：对于非整数值，直接使用 round() 函数。这是最常用的方法，因为它简单且符合计数数据的直觉。

  adjusted_counts_rounded <- round(adjusted_counts)
  

  但要注意，对于绝对值小于0.5的负值，round() 后会变成0，这是可以接受的（视为表达量极低或检测不到）。但对于比如-1.4这样的值，round()后会变成-1，这仍然是一个负的“计数”，不符合逻辑。所以单纯四舍五入可能不够。

• 非负截断与取整组合拳：我更喜欢用这个组合策略，分两步走：

  # 第一步：将所有负值设置为0（或一个很小的正数，如0.001，但通常0更安全）
  adjusted_counts_nonneg <- adjusted_counts
  adjusted_counts_nonneg[adjusted_counts_nonneg < 0] <- 0
  
  # 第二步：将整个矩阵四舍五入为整数
  adjusted_counts_final <- round(adjusted_counts_nonneg)
  
  # 检查一下还有没有漏网之鱼
  sum(adjusted_counts_final < 0) # 应该为0
  any(adjusted_counts_final != round(adjusted_counts_final)) # 应该为FALSE
  

  这个方法粗暴但有效，特别适合那些由计算精度引起的、非常接近0的负值和小数。它的优点是最大限度地保留了ComBat-seq校正后的数据格局，同时满足了下游工具对整数输入的要求。风险在于，如果原始负值或小数值很大，这种截断可能会扭曲部分基因的表达信息，尤其是低表达基因。

### 3.2 策略二：基于分布的重新缩放（Scaling）

当数值偏离整数较远，或者你希望处理方式更“统计”一点时，可以考虑缩放。ComBat-seq校正可能会改变数据整体的尺度（scale）。我们可以通过缩放，让每个样本的库大小（library size）回归到一个合理范围，同时促使数值向整数靠拢。

一个常见的方法是，以校正后数据的中位数或75%分位数（更稳健）为基准进行缩放：

# 计算每个样本的缩放因子（以75%分位数为基准）
scale_factors <- apply(adjusted_counts, 2, function(x) quantile(x[x > 0], probs = 0.75))
# 或者使用中位数 median(x[x>0])
# 进行缩放，使所有样本的75%分位数对齐
adjusted_counts_scaled <- t(t(adjusted_counts) / scale_factors) * median(scale_factors)

# 缩放后，再进行四舍五入和截断
adjusted_counts_scaled_nonneg <- adjusted_counts_scaled
adjusted_counts_scaled_nonneg[adjusted_counts_scaled_nonneg < 0] <- 0
adjusted_counts_scaled_final <- round(adjusted_counts_scaled_nonneg)

这种方法在去除批次效应的同时，还能标准化样本间的测序深度差异，相当于把 ComBat-seq 和 TMM（edgeR用的方法）或 median-of-ratios（DESeq2用的方法）的归一化思想结合了一下。需要注意，缩放可能会改变离散度，对于依赖精确离散度估计的 DESeq2 来说，最好在其内部进行归一化，所以使用此策略后，在运行 DESeq2 时可能需要设置 normalizationFactor 参数。

### 3.3 策略三：转换数据格式，适配广义线性模型

如果我们不想动 ComBat-seq 产出的数据本身，还有一个思路是“曲线救国”：转换下游分析工具的数据输入格式。DESeq2 和 edgeR 虽然最擅长处理原始计数，但它们的内核是广义线性模型（GLM）。对于非负的连续数据（比如我们的校正后数据，在将负值设为0后），理论上是可以尝试的，但需要格外小心。

• 对于DESeq2：它最核心的 DESeqDataSetFromMatrix 函数严格要求整数。一个变通的方法是，使用 DESeqDataSetFromTximport 如果你有转录本丰度数据，但这不适用这里。更实际的 hack 方法是，将处理后的连续数据乘以一个很大的数（比如1e6），再取整，模拟“伪计数”。但这会极大扭曲离散度估计，非常不推荐，除非万不得已且数据异常值极少。

  # 警告：此方法风险极高，仅作知识性展示
  pseudo_counts <- round(adjusted_counts_final * 1e6)
  # 然后使用 pseudo_counts 创建 DESeqDataSet
  

• 对于limma-voom：这是更好的选择！limma 配合 voom 转换，本身就是为处理归一化后的连续表达量数据（log-CPM）而设计的。我们的处理流程可以调整为：

  1. 对原始计数数据使用 ComBat-seq 校正。
  2. 对校正后的数据（进行非负截断和取整后）使用 edgeR 的 calcNormFactors 计算归一化因子。
  3. 使用 voom 函数将归一化的计数转换为log2-CPM值及其权重，这部分数据允许是连续的。
  4. 将 voom 转换后的数据（E 列表元素）输入 limma 进行线性模型拟合和差异分析。

  library(limma)
  library(edgeR)
  # 假设 adjusted_counts_final 是经过策略一处理后的非负整数矩阵
  dge <- DGEList(counts=adjusted_counts_final)
  dge <- calcNormFactors(dge) # TMM归一化
  design <- model.matrix(~ your_group_variable) # 你的实验设计矩阵
  v <- voom(dge, design, plot=TRUE) # voom转换，并绘制均值-方差趋势图检查
  fit <- lmFit(v, design) # 线性拟合
  fit <- eBayes(fit) # 经验贝叶斯收缩
  topTable(fit, coef=2) # 查看差异结果
  

  这条路走通了，就完美避开了整数约束问题。voom 的权重计算能很好地处理连续数据中的异方差性。

### 3.4 策略四：回退与替代方案（重新审视预处理）

如果以上方法都感觉不理想，或者数据异常情况太严重，我们可能需要退一步，思考是不是前置步骤或工具选择有问题。

• 检查输入数据：ComBat-seq 要求输入是原始计数。你是否不小心输入了经过CPM、TPM或log转换的数据？这会导致模型拟合出错。确保你的 count_matrix 是纯粹的整数矩阵。
• 调整ComBat-seq参数：ComBat_seq 函数有一些参数可以微调。例如，shrink 参数控制是否收缩批次效应估计值（默认为TRUE），对于小样本数据集，设为 FALSE 可能更稳定。shrink.disp 参数控制是否收缩离散度。可以尝试不同的参数组合，看看输出是否改善。

  adjusted_v2 <- ComBat_seq(count_matrix, batch=batch, group=NULL, shrink = FALSE)
  

• 考虑其他批次校正工具：如果 ComBat-seq 始终产出不理想的数据，可以试试其他专门为计数数据设计或对整数要求不那么严格的方法：
  • RUVseq：利用控制基因（如管家基因或通过经验确定的基因）来估计和移除不需要的变异。它对输出数据的整数性没有保证，但后续通常接 edgeR 或 DESeq2 的GLM模式，兼容性较好。
  • sva 包的 svaseq：基于奇异值分解（SVD）估计替代变量（SVs），然后将SVs作为协变量加入 DESeq2 或 edgeR 的模型中。这是一种“校正于模型内”的策略，完全避免了直接修改计数数据，是我在棘手情况下非常偏爱的方法。

  library(sva)
  # 使用svaseq估计替代变量
  mod <- model.matrix(~ your_condition_of_interest)
  mod0 <- model.matrix(~1, data=your_pheno_data)
  svseq <- svaseq(count_matrix, mod, mod0, n.sv=num_sv) # num_sv需要估计
  # 然后将 svseq$sv 作为协变量加入DESeq2或edgeR的设计矩阵中
  

4. 下游分析适配：与DESeq2、edgeR无缝对接

数据修复好了，最终目的是要能顺利、可靠地跑差异分析。这里分别说说怎么和我们最常用的 DESeq2 和 edgeR 对接。

### 4.1 与DESeq2的适配要点

DESeq2 对输入数据要求最严格，也最敏感。我们的目标是创建一个能被 DESeqDataSetFromMatrix 接受的“干净”的整数矩阵。

1. 数据准备：采用 3.1节 的“非负截断+四舍五入”策略，得到 adjusted_counts_final。
2. 创建DESeqDataSet：这是最关键的一步，要确保 colData 中的批次信息和你用来做 ComBat-seq 校正的批次信息一致。注意：此时批次效应已经被“消除”在数据中了，所以设计公式（design）中不应该再包含批次变量，否则等于校正了两次。设计公式应只包含你感兴趣的生物学条件。

   library(DESeq2)
   # 假设 pheno_data 是样本信息数据框，包含‘condition’列
   dds <- DESeqDataSetFromMatrix(countData = adjusted_counts_final,
                                 colData = pheno_data,
                                 design = ~ condition) # 注意：没有batch！
   

3. 运行DESeq：之后按标准流程进行即可。因为数据已经过校正，DESeq2 内部的估计离散度和对数折叠变化收缩可能会受到轻微影响。对于特别严谨的分析，可以考虑在运行 DESeq 时关闭自动的 Cook's 距离过滤（cooksCutoff = FALSE），因为校正可能改变了异常值的分布。

   dds <- DESeq(dds, cooksCutoff = FALSE)
   resultsNames(dds)
   res <- results(dds, contrast=c("condition", "treated", "control"))
   

### 4.2 与edgeR的适配要点

edgeR 相对灵活一些，它同样接受整数计数矩阵，但其强大的GLM框架也能在一定程度上处理经过合理处理的连续数据（尤其是结合 limma-voom 时）。

1. 标准流程（基于计数）：如果使用处理后的整数矩阵，流程和DESeq2类似，但模型构建方式不同。

   library(edgeR)
   y <- DGEList(counts=adjusted_counts_final, group=pheno_data$condition)
   y <- calcNormFactors(y) # 计算归一化因子（TMM）
   design <- model.matrix(~condition, data=pheno_data) # 设计矩阵同样不含batch
   y <- estimateDisp(y, design) # 估计离散度
   fit <- glmQLFit(y, design) # 拟合准似然广义线性模型
   qlf <- glmQLFTest(fit, coef=2) # 检验
   topTags(qlf)
   

2. limma-voom 流程（推荐用于处理后的数据）：正如 3.3节 所述，这是处理经过 ComBat-seq 校正（及后续处理）数据的非常稳健的流程。voom 转换将计数数据（或近似计数数据）转换为适用于线性模型的连续数据，并计算观测权重，完美地桥接了 edgeR 的归一化能力和 limma 的线性建模能力。当对 ComBat-seq 输出数据的整数性心存疑虑时，这是我首选的差异分析流程。

### 4.3 结果验证与敏感性分析

无论采用哪种修复策略和下游工具，做完差异分析后，验证是必不可少的。我通常会做以下几件事：

• 检查关键基因：看看你已知的、预期会差异表达的基因（比如阳性对照）是否在结果列表中，其变化方向是否符合预期。
• PCA/聚类图：对修复后的数据（或 voom 转换后的数据）做PCA，可视化检查批次效应是否真的被移除，样本是否按生物学条件聚类。
• 敏感性分析：尝试不同的修复策略（比如，将负值设为0 vs 设为0.001；使用四舍五入 vs 向下取整）。比较不同策略得到的差异基因列表的重合度（如韦恩图）。如果核心的、高显著性的差异基因集合在不同方法间高度一致，那说明你的分析结果是稳健的。如果差异很大，就需要深入挖掘原因，并谨慎解释结果。

5. 避坑指南与最佳实践建议

踩了这么多坑，我也总结出一些经验，希望能帮你少走弯路。

第一，预防优于治疗。 在运行 ComBat-seq 之前，务必确保输入数据是原始计数，并且已经过滤掉了在所有样本中表达量都极低（比如总和小于10）的基因。这些低表达基因噪音大，更容易在校正后产生异常值。可以使用 edgeR 的 filterByExpr 函数进行智能过滤。

第二，理解工具的输出。 ComBat-seq 的论文和文档虽然说了它输出整数，但那是在理想情况下。在实际的数值计算世界里，接近整数和等于整数是两回事。不要把它当成一个黑箱，要养成检查输出数据的习惯。

第三，策略选择有优先级。 我的个人经验是：如果异常值极少（<0.1%），直接使用“非负截断+四舍五入”，简单高效。如果异常值稍多，或者你想更保守一点，可以尝试“缩放+取整”。如果数据情况复杂，或者你对统计模型更放心，那么转向 limma-voom 管道是最稳健的选择，它避免了直接修改计数数据带来的理论尴尬。

第四，记录与透明。 在分析报告或代码注释中，务必清晰记录你使用了 ComBat-seq 进行批次校正，并详细说明后续处理非整数/负值所采取的具体步骤（例如：“将ComBat-seq输出中小于0的值设置为0，随后对所有值四舍五入取整，以符合DESeq2输入要求”）。这保证了分析的可重复性。

最后，保持怀疑精神。 批次效应校正是一个强大的技术，但也不是万能的。如果校正后的数据仍然“奇奇怪怪”，或者下游分析结果与生物学常识严重不符，不妨回过头来，重新审视实验设计、批次定义、甚至原始数据的质量。有时候，问题可能不在分析方法，而在数据本身。生信分析，一半是技术，一半是对数据的理解和洞察。