使用kBET检测批次效应 by 单细胞天地

生信技能树的老师指出来了：本来就不应该从数据分析角度解决的批次效应，其实并不完全是这样的，也有些时候实验设计没办法克服困难，数据分析人员就只能说是来收拾烂摊子！

以往的推文中我们没有对批次效应有个量化的概念，往往是画个PCA或者单细胞中做个UMAP、tSNE肉眼看看，上周组会我注意到师兄讲文献提到了这个（使用kBET检测批次效应）方法但没有解释，刚好我就自己学习学习，这周周更补全我们这个板块的空白

github仓库：https://github.com/theislab/kBET

本文主要参考仓库示例Rmd文档和原始文献 A test metric for assessing single-cell RNA-seq batch correction [https://www.nature.com/articles/s41592-018-0254-1]

使用这篇文章的数据集作为示例：

Batch effects and the effective design of single-cell gene expression studies [https://www.nature.com/articles/srep39921]

单细胞RNA测序（scRNA-seq）可用于以高分辨率表征基因表达水平的变化。然而，scRNA-seq中的实验噪声源尚不清楚。我们研究了与使用单细胞Fluidigm C1平台进行样品处理相关的技术变化。为此，我们处理了来自三个人诱导多能干细胞（iPSC）系的三个C1重复。我们向所有样品添加了唯一分子标识符（UMI），以解决扩增偏差。我们发现基因表达数据变异的主要来源是由基因型驱动的，但我们也观察到技术重复之间的巨大差异。我们观察到，使用UMIs将读段转换为分子受到生物学和技术差异的影响，这表明UMI计数不是基因表达水平的无偏估计指标。根据我们的结果，我们提出了一个有效的scRNA-seq研究框架。

https://github.com/jdblischak/singleCellSeq/tree/master/data

kBET原理：

A test metric for assessing single-cell RNA-seq batch correction [https://www.nature.com/articles/s41592-018-0254-1]

用于评估单细胞RNA-seq批次校正的测试指标

(Abstract)
使用kBET来评估常用的批次回归和归一化方法，并量化它们在保留生物变异性的同时消除批次效应的程度
...
(Main)
在这项研究中，我们将kBET应用于使用基于微孔板和基于液滴的方法（每批100-3,000个细胞）分析来自研究的四个小鼠单细胞数据集，并评估了11种归一化和7种批次效应回归方法的性能和准确性（图1e）。基于对数（计数 + 1）、对数（每百万计数（CPM） + 1）或 scran 池的批次校正，以及 ComBat 或 limma 回归，在保留所有数据集的生物结构的同时降低了批次效应（表 1）。最后，我们探索了kBET评估独立研究整合的潜力，并确定kBET还允许人们研究复杂人体组织数据中的个体间变异性

从图a,b中我们可以看到因为技术偏差导致的批次效应对实验设计的影响

我们之前无论在单细胞还是bulk中都提到过这个问题

单细胞参考：

在harmony、不harmony，这是个问题这篇中我们着重讨论了harmony以及单细胞何时需要处理批次效应

在多分组单细胞测序数据第一层次未整合和整合分析对B细胞细分的分群有何影响？这篇中我们也顺带提了一下CCA方法

bulk参考：

在奇怪的转录组差异表达矩阵之实验分组这篇中，我们强调了并不是所有的批次效应都可以被矫正

比如图b右边的confounded实验设计，批次效应和contrl vs condition混在了一起就不可以使用我们过往介绍的方法去除，这会导致生物学差异也被去除

作者这里提出kBET方法来通过总体拒绝率来评估批次效应的大小

kBET方法，简单直观来说，就是在降维聚类图中选取固定大小的随机邻域，基于卡方分布看这个随机领域是否混合良好（如上图b中左边和图c所示），因为随机邻域如果具有与完整数据集相同的批次标签分布则能说明混合良好，获得每个邻域的二元测试结果，然后对其进行平均以计算总体拒绝率。

对二元测试结果（binary test results）的解释：
"二元测试结果"是指kBET方法中固定大小的随机邻域的基于χ2的测试结果。
基于χ2的检验将相邻样本中批次标签的分布与整个数据集进行比较。测试结果是二元结果，表明相邻样本是否混合良好（表示低批量效应）或混合不好（表示高批量效应）。
获得每个邻域的二元测试结果，然后对其进行平均以计算总体拒绝率。低拒绝率表明相邻样本的批次标签分布与完整数据集相似，表明重复混合良好。

此外，原文中介绍了许多normalization、batch regression和其他批次效应评估的方法，并基于kBET和其他批次效应评估对不同normalization、batch regression处理效果进行了比较，很值得一读

kBET

准备数据集

library(devtools)
# install_github('theislab/kBET')
# install_local("kBET-master.zip",repos = F, dependencies = T)

#########准备数据集############

# 数据来源：
umi <- read.table("data/molecules.txt")
# 表达矩阵 行为umi基因 列为barcode细胞
# > umi[1:4,1:4]
# NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03 NA19098.r1.A04
# ENSG00000237683              0              0              0              1
# ENSG00000187634              0              0              0              0
# ENSG00000188976              3              6              1              3
# ENSG00000187961              0              0              0              0
anno <- read.table("data/annotation.txt",sep = "\t",header = TRUE)
# 注释信息
# > head(anno)
# individual replicate well      batch      sample_id
# 1    NA19098        r1  A01 NA19098.r1 NA19098.r1.A01
# 2    NA19098        r1  A02 NA19098.r1 NA19098.r1.A02
# 3    NA19098        r1  A03 NA19098.r1 NA19098.r1.A03
# 4    NA19098        r1  A04 NA19098.r1 NA19098.r1.A04
# 5    NA19098        r1  A05 NA19098.r1 NA19098.r1.A05
# 6    NA19098        r1  A06 NA19098.r1 NA19098.r1.A06

# 通过质控的细胞barcode id
quality_single_cells <- scan("data/quality-single-cells.txt", what = "character")
# > head(quality_single_cells)
# [1] "NA19098.r1.A01" "NA19098.r1.A02" "NA19098.r1.A04" "NA19098.r1.A05" "NA19098.r1.A06" "NA19098.r1.A07"

# 过滤掉低质量细胞
umi.qc <- umi[, colnames(umi) %in% quality_single_cells]
anno.qc <- anno[anno$sample_id %in% quality_single_cells, ]

# 在原始数据集中作者三个样本每个都有3个重复,但是由于NA19098的r2重复没有通过质控所以被排除
table(anno.qc$individual, anno.qc$replicate)
# r1 r2 r3
# NA19098 85  0 57
# NA19101 80 70 51
# NA19239 74 68 79

exp.design <- as.matrix(table(anno.qc$individual, anno.qc$replicate))
knitr::kable(exp.design, booktabs=TRUE, caption='Experimental design', row.names =TRUE)
# Table: Experimental design
# |        | r1| r2| r3|
# |:-------|--:|--:|--:|
# |NA19098 | 85|  0| 57|
# |NA19101 | 80| 70| 51|
# |NA19239 | 74| 68| 79|

# 去除读数为0的基因
umi.qc <- umi.qc[rowSums(umi.qc) > 0, ]
dim(umi.qc)
# 18464   564
# 可以看到剩下了564个细胞 18464个基因

# Finally, we create a vector with the IDs of *ERCC spike-in* genes. We keep them in the dataset, but it is prudent to remove them before using kBET or some normalisation approaches.
spikes <- grep('ERCC', rownames(umi.qc))

测试与比较

library(kBET)
patients <- unique(anno.qc$individual)
kBET.umi.counts <- list()
for (patient in patients){
  kBET.umi.counts[[patient]] <- kBET(df = umi.qc[-spikes,anno.qc$individual==patient],   # 非棘突蛋白
                                     batch = anno.qc$replicate[anno.qc$individual==patient],  # 当前样本各细胞对应的批次 r1-3
                                     plot = FALSE)
}

kBET 函数为这个包的核心函数，“kBET runs a chi square test to evaluate the probability of a batch effect.”，它返回kBET结果的summary信息：

summary - a rejection rate for the data, an expected rejection rate for random labeling and the significance for the observed result
results - detailed list for each tested cells; p-values for expected and observed label distribution
average.pval - significance level over the averaged batch label distribution in all neighbourhoods
stats - extended test summary for every sample
params - list of input parameters and adapted parameters, respectively
outsider - only shown if adapt=TRUE. List of samples without mutual nearest neighbour:

index - index of each outsider sample)
categories - tabularised labels of outsiders
p.val - Significance level of outsider batch label distribution vs expected frequencies. If the significance level is lower than alpha, expected frequencies will be adapted

# kBET返回每个患者样本的测试摘要summary和以下结果图：
library(ggplot2)
library(gridExtra)

kBET.plots <- list()
for (patient in patients){
  plot.data <- data.frame(class = rep(c('observed', 'expected'), each=100), 
                          data = c(kBET.umi.counts[[patient]]$stats$kBET.observed,
                                   kBET.umi.counts[[patient]]$stats$kBET.expected))
  kBET.plots[[patient]] <- ggplot(plot.data, aes(class, data)) + 
    geom_boxplot() + 
    theme_bw() + 
    labs(x='Test', y='Rejection rate',title=patient) + 
    scale_y_continuous(limits=c(0,1))
}

n <- length(kBET.plots)
do.call("grid.arrange", c(kBET.plots, ncol=n))

该算法运行一个null model，其中批次标签是随机排列的。使用null model，我们估计了混合良好的数据集的预期拒绝率。观察到的拒绝率使用样品的实际批次标签，并描述了批次效应引起的偏差。默认情况下，kBET只测试样本的一个子集的良好混合性，并重复过程“n_repeat”次以创建显示的统计数据。我们使用统计数据来计算拒绝率的显著性，并将其添加到kBET摘要中。在这里，我们进一步显示一名患者的摘要信息：

kBET中的零模型（null model）包括随机排列批次标签以估计预期的拒绝率，而观察到的拒绝率使用实际的批次标签来测量批次效应引起的偏差
我个人理解就是非参数方法，可以参考基因集分析的前世今生（附进行通路富集分析的9个tips）这篇推文中的gene sampling 和 phenotype permutation

summary.kBET <- kBET.umi.counts[[patients[1]]]$summary
summary.kBET$kBET.expected <- signif(summary.kBET$kBET.expected, digits = 2)
summary.kBET$kBET.observed <- signif(summary.kBET$kBET.observed, digits = 2)
colnames(summary.kBET) <- c('kBET (null model)', 'kBET (observed)', 'kBET p_value')
knitr::kable(summary.kBET, 
             booktabs=TRUE, caption=paste0('Summary for ', patient[[1]]), row.names=TRUE)
# Table: Summary for NA19239
# |      | kBET (null model)| kBET (observed)| kBET p_value|
#   |:-----|-----------------:|---------------:|------------:|
#   |mean  |                 0|            0.55|            0|
#   |2.5%  |                 0|            0.38|            0|
#   |50%   |                 0|            0.56|            0|
#   |97.5% |                 0|            0.72|            0|

从summary结果和箱型图我们可以看到，kBET方法认为这三个样本都存在批次效应的影响

这里其实需要注意一下这个单细胞实验分组是没有对照的，所以不存在一个根据condition分组差异分析，生物学差异和批次效应混在一起的情况，这里就是单纯地看批次效应

library(ggplot2)
pca.umis <- prcomp(log10(1+t(umi.qc[-spikes,])))
pca.df <- data.frame(PC1=pca.umis$x[,1], 
                     PC2=pca.umis$x[,2], 
                     replicate= anno.qc$replicate, 
                     patient=anno.qc$individual)
ggplot(pca.df, aes(x=PC1, y=PC2, shape=replicate, colour=patient)) + geom_point() + theme_bw() +
  ggtitle(expression(PCA * ' '* of* ' '* log[10]-normalised * ' '*raw* ' '* counts.))

从老方法PCA图中我们也可以看到确实每个样本都存在批次效应，样本重复没有很好地混在一起

这个示例其实简化了问题，因为单细胞大部分时候也是有实验分组的，这时就不能仅仅考虑批次效应，还需要考虑对生物学差异的影响

后面的部分作者使用kBET来评估常用的批次回归和归一化方法，并量化它们在保留生物变异性的同时消除批次效应的程度

得出normalization方法使得表达谱同分布来解决批次效应会对下游分析产生不利影响，而batch regression方法可以更加有效地处理批次效应的结论

小结

我们简单地介绍了kBET方法并演示了基础代码

需要注意的是，前提条件要知道批次标签才能计算kBET，这个我们可以简单地根据样本重复来打标签，也可以根据自己掌握的技术细节来进行更有意义的探索

原文中介绍了许多normalization、batch regression和其他批次效应评估的方法，并基于kBET和其他批次效应评估对不同normalization、batch regression处理效果进行了比较，很值得一读

另外这个方法的提出目的是用于scRNAseq检测批次效应，我个人想了想可能是因为bulk我们是以每个个体为单位，落到降维聚类图上就没单细胞以细胞为单位那么多，所以不用这个方法，但如果思路打开面对其它高维度数据，kBET背后的统计检验方法不失为一种简单但合理的思路

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

ixxmu / mp_duty

使用kBET检测批次效应 #4035