183_RNA-seq批次效应校正¶

一句话概述¶

RNA-seq批次效应是由实验操作、测序时间、试剂批号等非生物学因素引起的系统性技术偏差，通过ComBat/ComBat-seq、limma::removeBatchEffect、Harmony等方法校正后才能进行跨批次样本的合理比较。

核心知识点表格¶

知识点	说明
批次效应来源	建库时间、操作员、测序lane、试剂批号、RNA提取方法
检测方法	PCA着色、RLE图、层次聚类
ComBat	sva包，经验贝叶斯方法，最经典的批次校正
ComBat-seq	专为RNA-seq counts设计的ComBat变体
removeBatchEffect	limma包，在差异分析前校正可视化用
SVA	隐变量分析，无需已知批次信息
RUVSeq	利用对照基因去除不需要的变异
注意事项	批次与实验设计完全混淆时无法校正

步骤详解¶

第一步：检测批次效应¶

白话解释：在校正之前，首先需要确认批次效应是否存在。最直观的方法是做PCA，看样本是否按批次分离而非按实验分组分离。

技术细节：检测方法包括：(1)PCA/t-SNE按批次着色；(2)相对对数表达（RLE）图检查分布一致性；(3)层次聚类看样本是否先按批次聚类。

library(DESeq2)
library(ggplot2)
library(pheatmap)

# 读取数据
counts <- read.csv("count_matrix.csv", row.names = 1)
coldata <- read.csv("sample_info.csv", row.names = 1)

# 确保有batch信息
head(coldata)
#          condition batch
# sample1  Control   batch1
# sample2  Control   batch1
# sample3  Treatment batch1
# sample4  Control   batch2
# sample5  Treatment batch2

# DESeq2标准化
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData = coldata,
                               design = ~ condition)
vsd <- vst(dds, blind = TRUE)  # 方差稳定化变换

# PCA检查批次效应
pca_data <- plotPCA(vsd, intgroup = c("condition", "batch"), returnData = TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"))

ggplot(pca_data, aes(x = PC1, y = PC2, color = batch, shape = condition)) +
  geom_point(size = 4) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  ggtitle("PCA: Before Batch Correction") +
  theme_minimal()

# 样本间相关性热图
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)

annotation_col <- data.frame(
  Batch = coldata$batch,
  Condition = coldata$condition,
  row.names = rownames(coldata)
)

pheatmap(sampleDistMatrix,
         annotation_col = annotation_col,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         main = "Sample Distance (Before Correction)")

第二步：使用ComBat-seq校正counts数据¶

白话解释：ComBat-seq是专门为RNA-seq原始count数据设计的批次校正方法。它保持数据的count性质（整数、非负），校正后的数据可以直接用于DESeq2或edgeR分析。

技术细节：ComBat-seq基于负二项分布建模，使用经验贝叶斯方法估计和去除批次效应。与原始ComBat不同，它直接在count层面操作，不需要先做log变换，输出仍为整数counts。

library(sva)

# ComBat-seq校正
adjusted_counts <- ComBat_seq(
  counts = as.matrix(counts),
  batch = coldata$batch,
  group = coldata$condition,  # 保护的生物学变量
  covar_mod = NULL            # 额外协变量
)

# 验证校正效果
dds_adj <- DESeqDataSetFromMatrix(countData = adjusted_counts,
                                    colData = coldata,
                                    design = ~ condition)
vsd_adj <- vst(dds_adj, blind = TRUE)

pca_adj <- plotPCA(vsd_adj, intgroup = c("condition", "batch"), returnData = TRUE)
percentVar_adj <- round(100 * attr(pca_adj, "percentVar"))

ggplot(pca_adj, aes(x = PC1, y = PC2, color = batch, shape = condition)) +
  geom_point(size = 4) +
  xlab(paste0("PC1: ", percentVar_adj[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar_adj[2], "% variance")) +
  ggtitle("PCA: After ComBat-seq Correction") +
  theme_minimal()

# 使用校正后的counts进行差异分析
dds_adj <- DESeq(dds_adj)
res <- results(dds_adj, contrast = c("condition", "Treatment", "Control"))

第三步：使用limma::removeBatchEffect¶

白话解释：limma的removeBatchEffect是最简单的批次校正方法，直接从表达矩阵中减去批次效应。但注意，它只用于可视化和探索性分析，差异分析时应在模型中加入batch作为协变量。

技术细节：removeBatchEffect使用线性模型拟合批次效应并移除，同时保留design中指定的生物学效应。它作用于连续值矩阵（log-counts或VST变换后的数据）。

library(limma)

# 对VST变换后的数据校正（用于可视化）
mat <- assay(vsd)
mat_corrected <- removeBatchEffect(
  mat,
  batch = coldata$batch,
  design = model.matrix(~ condition, data = coldata)
)

# 可视化校正后的PCA
pca_result <- prcomp(t(mat_corrected))
pca_df <- data.frame(
  PC1 = pca_result$x[, 1],
  PC2 = pca_result$x[, 2],
  Batch = coldata$batch,
  Condition = coldata$condition
)
var_pct <- round(summary(pca_result)$importance[2, 1:2] * 100, 1)

ggplot(pca_df, aes(x = PC1, y = PC2, color = Condition, shape = Batch)) +
  geom_point(size = 4) +
  labs(x = paste0("PC1 (", var_pct[1], "%)"),
       y = paste0("PC2 (", var_pct[2], "%)"),
       title = "After removeBatchEffect") +
  theme_minimal()

# 差异分析中正确的做法：在模型中加入batch
design <- model.matrix(~ batch + condition, data = coldata)
dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ batch + condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "Treatment", "Control"))

第四步：使用SVA发现隐藏的批次效应¶

白话解释：有时批次效应的来源不清楚，或者存在未记录的技术混淆因素。SVA（替代变量分析）可以自动发现这些隐藏的变异来源。

library(sva)

# 准备模型矩阵
mod <- model.matrix(~ condition, data = coldata)
mod0 <- model.matrix(~ 1, data = coldata)

# 估计替代变量数量
n_sv <- num.sv(as.matrix(counts), mod, method = "be")
cat("估计的替代变量数:", n_sv, "\n")

# 计算替代变量
svobj <- svaseq(as.matrix(counts), mod, mod0, n.sv = n_sv)

# 将SV加入DESeq2模型
coldata$SV1 <- svobj$sv[, 1]
coldata$SV2 <- svobj$sv[, 2]

dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ SV1 + SV2 + condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "Treatment", "Control"))

第五步：Python实现批次校正¶

import scanpy as sc
import numpy as np
import pandas as pd
from combat.pycombat import pycombat
import matplotlib.pyplot as plt

# 读取数据
counts = pd.read_csv("count_matrix.csv", index_col=0)
metadata = pd.read_csv("sample_info.csv", index_col=0)

# 使用pycombat校正
corrected = pycombat(
    counts,
    batch=metadata['batch'].values,
    mod=pd.get_dummies(metadata['condition'])
)

# 可视化
from sklearn.decomposition import PCA

pca = PCA(n_components=2)
log_counts = np.log2(counts.T + 1)
pca_before = pca.fit_transform(log_counts)

log_corrected = np.log2(corrected.T + 1)
pca_after = pca.fit_transform(log_corrected)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
for batch in metadata['batch'].unique():
    mask = metadata['batch'] == batch
    axes[0].scatter(pca_before[mask, 0], pca_before[mask, 1], label=batch, alpha=0.7)
    axes[1].scatter(pca_after[mask, 0], pca_after[mask, 1], label=batch, alpha=0.7)
axes[0].set_title('Before Correction')
axes[1].set_title('After Correction')
for ax in axes:
    ax.legend()
    ax.set_xlabel('PC1')
    ax.set_ylabel('PC2')
plt.tight_layout()
plt.savefig('batch_correction.png', dpi=300)

实战命令速查¶

# ComBat-seq一行校正
adjusted <- sva::ComBat_seq(as.matrix(counts), batch = coldata$batch, group = coldata$condition)

# limma校正可视化矩阵
corrected <- limma::removeBatchEffect(assay(vsd), batch = coldata$batch)

# DESeq2模型中加入批次
dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ batch + condition)

# SVA自动发现批次
svobj <- sva::svaseq(as.matrix(counts), mod, mod0)

# RUVSeq
library(RUVSeq)
set <- newSeqExpressionSet(as.matrix(counts), phenoData = coldata)
set <- RUVg(set, control_genes, k = 1)

面试常问点¶

Q1: 批次效应和生物学差异如何区分？ A: 关键在于实验设计。如果所有Treatment样本在batch1、所有Control在batch2，则批次效应和生物学差异完全混淆，无法区分。好的设计应该让每个批次都包含所有实验组。可以通过PCA查看样本是按batch还是condition分离来初步判断。

Q2: ComBat和ComBat-seq的区别是什么？ A: ComBat作用于连续值（如log-CPM），输出也是连续值，适合芯片数据。ComBat-seq专为RNA-seq设计，直接作用于count数据，基于负二项分布建模，输出整数counts，可直接输入DESeq2/edgeR。

Q3: removeBatchEffect可以用于差异分析吗？ A: 不推荐。removeBatchEffect直接修改了表达矩阵，这会改变数据的方差结构，导致差异分析的统计检验失效（p值不可信）。正确做法是在差异分析模型中加入batch作为协变量（如~ batch + condition）。

Q4: 批次和实验设计完全混淆时怎么办？ A: 无法用统计方法解决。任何校正方法都可能在去除批次效应的同时去除真实的生物学信号。唯一的解决方案是改进实验设计（平衡设计）或增加独立实验重复。

Q5: SVA和已知批次校正哪个更好？ A: 如果批次信息明确可靠，直接使用ComBat-seq或在模型中加入batch更简单有效。SVA适用于批次信息未知或存在多个未记录的技术混淆时。两者也可以结合使用（先校正已知batch，再用SVA发现残余的隐藏因素）。

Q6: 批次校正会不会去除真实的生物学信号？ A: 有可能。特别是当batch和condition存在一定程度的混淆时，校正可能部分移除生物学信号。ComBat和ComBat-seq通过group参数保护已知的生物学变量来缓解这个问题。

易错点¶

batch和condition完全混淆仍强行校正：结果不可信，应从实验设计层面解决
对原始counts使用ComBat：应使用ComBat-seq，ComBat需要log变换后的数据
removeBatchEffect后做差异分析：应在模型中加入batch协变量
过度校正：加入太多协变量可能移除真实信号
忽略batch平衡性：校正前应检查每个batch中各组的样本数是否平衡

补充知识¶

批次校正方法对比¶

方法	输入	原理	适用场景
ComBat-seq	Counts	经验贝叶斯+负二项	已知batch的RNA-seq
ComBat	连续值	经验贝叶斯	芯片数据
removeBatchEffect	连续值	线性回归	仅可视化
SVA	Counts/连续	替代变量分析	未知batch来源
RUVSeq	Counts	利用对照基因	有对照基因时
Harmony	降维空间	迭代聚类	单细胞数据

实验设计建议¶

每个batch包含所有实验组（平衡设计）
记录所有可能的技术变量（日期、操作员、试剂批号）
加入技术重复检测batch效应大小
避免将对照和处理分在不同批次