🤩 scRNA-seq | 吐血整理的单细胞入门教程（Normalization的影响因素）（十一）
🤩 scRNA-seq | 吐血整理的单细胞入门教程（Normalization的影响因素）（十一）
1写在前面
上期我们介绍了使用scater, scran以及scRNA.seq.funcs包进行Normalization的方法，这种Normalization主要是针对library大小差异。😘
而在实际分析中，scRNAseq的影响因素可能不仅仅是library大小问题，还包括由于试剂、分离方法、实验者不同而引起的batch effects。🤗
本期我们介绍一下如何去除这些因素导致的noise。
2用到的包
rm(list = ls())
library(scRNA.seq.funcs)
library(scater)
library(scran)
library(sva)
library(batchelor)
library(kBET)
library(RColorBrewer)
library(reshape2)
library(patchwork)
3示例数据
这里我们用一下之前介绍的counts文件和annotation文件，然后通过SingleCellExperiment创建SingleCellExperiment格式的文件，并且经过过滤，ID转换等。
这里我们用logNormCounts进行一下上次介绍的Normalization，以去除library大小的影响。😘
load("umi_umiqc.Rdata")
umi.qc <- umi[! rowData(umi)$discard, ! colData(umi)$discard]
qclust <- quickCluster(umi.qc, min.size = 30)
umi.qc <- computeSumFactors(umi.qc, clusters = qclust)
umi.qc <- logNormCounts(umi.qc)
umi.qc
4Combat
这里我们先介绍一下combat包去除批次效应。🧐
4.1 去除批次效应-replicate
我们的dataset里有多个replicate，我们用Combat函数来去除一下batch effets吧。
assay(umi.qc, "combat") <- ComBat(logcounts(umi.qc),batch = umi.qc$replicate😉)
4.2 去除批次效应-detected
我们再换个因素试试！~
assay(umi.qc, "combat_tf") <- ComBat(logcounts(umi.qc),batch = umi.qc$detected)
5mnnCorrect (batchelor)
MNN，(mutual nearest neighbor)，主要思想是找到不同批次中相同的细胞类型，然后计算同种细胞类型的gene表达的差异，然后去除批次效应。
5.1 原理
这里我们补充一下原理:👇
- 1️⃣ 对不同批次的基因表达谱信息按细胞进行
余弦标准化（cosine normalization）; - 2️⃣ 依次计算不同批次中每个细胞到另一批次中所有细胞的
欧式距离，其实际等同于表达数据标准化前的余弦距离; - 3️⃣ 这个时候我们就得到了一个不同批次间细胞互相配对的数据; 😘
- 4️⃣ 计算细胞间的
基因表达差值，即表达差异向量，也称为配对特异的批次效应校正向量（pair-specific batch convection vector）; - 5️⃣ 利用
高斯核函数，计算它们的加权平均数，即批次效应校正向量，最后将其应用到所有细胞中进行批次效应的校正。
5.2 具体操作
mnn_out <- fastMNN(umi.qc,batch = umi.qc$replicate)
assay(umi.qc, "mnn") <- assay(mnn_out,'reconstructed')
6效果评估-PCA
这里我们写个for循环对不同Normalization方法进行比较，以便我们在实际应用中选择合适的方法。
lis <- list()
for(i in assayNames(umi.qc)) {
tmp <- runPCA(umi.qc, exprs_values = i, ncomponents = 20)
lis[[i]] <- plotPCA(
tmp,
colour_by = "batch",
size_by = "detected",
shape_by = "individual"
) +
ggtitle(i) +
theme(plot.title = element_text(size = 26))
}
wrap_plots(lis)
7效果评估-RLE
不知道大家还记不记得relative log expression (RLE)。
res <- list()
for(i in assayNames(umi.qc)) {
res[[i]] <- suppressWarnings(calc_cell_RLE(assay(umi.qc, i)))
}
par(mar=c(6,4,1,1))
boxplot(res, las=2)
Note! RLE只评估每个细胞高于和低于平均水平的基因数量是否相等，而批次间的随机噪音可能无法识别。
8效果评估-kNN
这里我们再介绍一种方法，kBET包利用kNN networks进行评估。🤯
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
Note! 在应用到有replicates的dataset时，需要分别应用到每个生物学分组中，所以在这里，我们提取的是individual。
👀 可视化~
# Plot results
dod <- melt(as.matrix(eff_debatching), value.name = "kBET")
colnames(dod)[1:2] <- c("Normalisation", "Individual")
colorset <- c('gray', brewer.pal(n = 9, "RdYlBu"))
ggplot(dod, aes(Normalisation, Individual, fill=kBET)) +
geom_tile() +
scale_fill_gradient2(
na.value = "gray",
low = colorset[2],
mid=colorset[6],
high = colorset[10],
midpoint = 0.5, limit = c(0,1)) +
scale_x_discrete(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
theme(
axis.text.x = element_text(
angle = 45,
vjust = 1,
size = 12,
hjust = 1
)
) +
ggtitle("Effect of batch regression methods per individual")
最后祝大家早日不卷!~
腾讯云开发者
