5.6 处理混杂因素
5.6.1 简介
在上一章中，我们根据文库大小进行了标准化，有效地将其作为混杂因素去除。现在，我们将从数据中去除其他定义不太明确的混杂因素。技术混杂因素（又称批次效应）可能来自试剂、分离方法、进行实验的实验室/实验者，甚至进行实验的日期或时间的差异。考虑技术混杂因素，特别是批次效应，是一个大课题，且涉及实验设计原理。这里我们讨论当实验设计合理时处理混杂因素的方法。
从根本上来说，考虑技术混杂因素涉及识别并去除表达数据中与感兴趣的生物信号不相关（即混淆）的变异源。有多种方法，其中一些使用ERCC或管家基因，一些使用内源基因。
使用ERCC作为控制基因很有优势，因为在实验中向每个细胞添加了相同数量的ERCC。原则上，我们观察到的这些基因的所有变异都是由于技术噪声造成的；而内源基因则受到技术噪声和生物变异的共同影响。可以通过对ERCC进行模型拟合并将其从内源基因中“减去”来消除技术噪声。基于此前提，有几种方法可用，每种方法都使用不同的模型和不同的拟合程序。或者识别超出技术噪声的显著变化的基因。但使用ERCC进行标准化存在一些问题（特别是来自细菌序列的ERCC），出于各种原因，它们的变异实际上可能高于内源基因。
由于使用ERCC的诸多问题，使用内源基因通常可以获得更好的结果。当大量内源基因在细胞间不会发生系统性变化，并且我们预期技术影响会影响大多数基因时，此类方法就能表现良好。
ScRNA-seq数据整合有两种场景。在第一种场景中，预期细胞组成是相同的，这时为Bulk RNA-seq开发的方法（例如ComBat）就能表现出良好的性能。对于同一实验的生物学重复，这通常是可行的；对于Tung数据中的批次，这也是可行的。在第二种情况下，数据之间的重叠是部分的——例如，如果数据代表健康和患病组织，其细胞类型组成就存在很大差异。在这种情况下，基于相互最近邻（MNN）的方法往往表现更好。
在这里，我们使用两种方法消除批次效应——基于经验贝叶斯框架的ComBat和来自batchelor包的基于MNN的方法fastMNN。
5.6.2 加载并标准化Tung数据
> library(scRNA.seq.funcs)
> library(scater)
> library(scran)
> library(sva)
> library(batchelor)
> library(kBET)
> set.seed(1234567)
读取预处理的数据集并使用scran包中的logNormCounts对其进行标准化。
在umi.qc对象中，除了之前存在的counts和logcounts_raw之外，还会出现一个名为logcounts的assay：
> umi <- readRDS("data/tung/umi.rds")
> umi.qc <- umi[! rowData(umi)$discard, ! colData(umi)$discard]
> qclust <- quickCluster(umi.qc, min.size = 30)
> umi.qc <- computeSumFactors(umi.qc, clusters = qclust)
Warning message:
In .guessMinMean(x, min.mean = min.mean, BPPARAM = BPPARAM) :
assuming UMI data when setting 'min.mean'
> umi.qc <- logNormCounts(umi.qc)
5.6.3 Combat
如果有一个平衡设计的实验，可以使用ComBat来消除批次效应，同时使用mod参数来保留生物学影响。然而，Tung数据包含多个实验重复，不是平衡设计，因此使用mod1来保留生物变异将导致错误。
> assay(umi.qc, "combat") <- ComBat(logcounts(umi.qc),batch = umi.qc$replicate)
Found 40 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
Found3batches
Adjusting for0covariate(s) or covariate level(s)
Standardizing Data across genes
Fitting L/S model and finding priors
Finding parametric adjustments
Adjusting the Data
5.6.4 mnnCorrect（batchelor）
fastMNN函数可以实现基于MNN的标准化。
> mnn_out <- fastMNN(umi.qc,batch = umi.qc$replicate)
> assay(umi.qc, "mnn") <- assay(mnn_out,'reconstructed')
5.6.5 消除批次效应方法的评估与比较
在考虑消除混杂因素的不同方法时，一个关键问题是如何定量地确定哪一种方法最有效。具有挑战的主要原因是，通常很难知道什么对应于技术混杂因素及什么是感兴趣的生物变异。这里，我们考虑三种不同的指标，这些指标根据我们对实验设计的了解都是合理的。根据想要解决的生物学问题，选择一个评估在特定情况下可能最受关注的混杂因素的指标非常重要。
5.6.5.1 评估1
通过查看PCA图来评估标准化的有效性，其中颜色对应技术重复、形状对应不同的生物样本（个体）。生物样品和批次的分离表明技术变异已被消除。我们始终使用log2-cpm归一化数据来匹配PCA的假设。
> for(n in assayNames(umi.qc)) {
tmp <- runPCA(umi.qc, exprs_values = n, ncomponents = 20)
print(
plotPCA(
tmp,
colour_by = "batch",
size_by = "detected",
shape_by = "individual"
) +
ggtitle(n)
)
}
5.6.5.2 评估2
还可以使用细胞的相对对数表达（RLE）检查批次校正的有效性，以确认技术噪声已从数据集中消除。注意：RLE仅评估每个细胞中高于和低于平均值的基因数量是否相等——即系统偏差。RLE可能无法检测到批次之间的随机技术噪声。
> res <- list()
> for(n in assayNames(umi.qc)) {
res[[n]] <- suppressWarnings(calc_cell_RLE(assay(umi.qc, n)))
> }
> par(mar=c(6,4,1,1))
> boxplot(res, las=2)
5.6.5.3 评估3
检查批次效应校正效果的另一种方法是考虑数据局部子样本中不同批次的点的混合。如果没有批次效应，那么任何局部区域的每个批次的细胞比例应该等于每个批次中细胞的全局比例。
kBET（Büttner等，2018）围绕随机细胞采用kNN网络，并根据二项分布测试每批细胞的数量。这些测试的拒绝率表明数据中仍然存在批次效应的严重性（高拒绝率=强批次效应）。kBET假设每个批次组成相同，因此只有在使用完美的平衡设计时它才能应用于整个数据集。但如果将kBET分别应用于每个分组，它也可以应用于重复数据。对于Tung数据，我们对每个个体独立应用kBET，以检查残留的批次效应。然而，这种方法无法识别与生物条件相混淆的残留批次效应。此外，kBET无法确定生物信号是否被保留。
> compare_kBET_results <- function(sce){
sce <- umi.qc
indiv <- unique(as.character(sce$individual))
norms <- assayNames(sce) # Get all normalizations
results <- list()
for (i in indiv){
for (j in norms){
tmp <- kBET(
df = t(assay(sce[,sce$individual== i], j)),
batch = sce$batch[sce$individual==i],
heuristic = TRUE,
verbose = FALSE,
addTest = FALSE,
plot = FALSE)
results[[i]][[j]] <- tmp$summary$kBET.observed[1]
}
}
return(do.call(rbind.data.frame, results))
}
> eff_debatching <- compare_kBET_results(umi.qc)
> eff_debatching
counts logcounts_raw logcounts combat mnn
NA19098 1.0000000 0.7340000 0.7324000 0.7924000 0.5024000
NA19101 0.4400000 0.4572000 0.8012000 0.7924000 0.7520000
NA19239 0.7140741 0.6740741 0.8811111 0.8703704 0.8737037
最后，让我们直观地看一下kBET计算的输出：
往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（1）
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（2）
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（3）
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（4）
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（5）
