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")

---

最后祝大家早日不卷!~

---

需要示例数据的小伙伴，在公众号回复scRNAseq获取吧！

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰

本文由 mdnice 多平台发布