单细胞交响乐
  • 前言:我与《单细胞交响乐》的缘分
  • 1 准备篇:背景知识
    • 1.1 数据结构
    • 1.2 总览 | 从实验到分析
  • 2 积累篇:文献阅读
    • 2.1.1 综述 | 2019-单细胞转录组分析最佳思路
    • 2.1.2 综述 | 2018-单细胞捕获平台
    • 2.1.3 综述 | 2017-scRNA中的细胞聚类分群
    • 2.1.4 综述 | scRNA已经开发出超过1000款工具了,你用过几种?
    • 2.1.5 综述 | 2021-单细胞测序的微流控技术应用
    • 2.2.1 研究 | 2018-单细胞转录组探索癌症免疫治疗获得性抗性机理
    • 2.2.2 研究 | 2018-人类结直肠癌单细胞多组学分析
    • 2.2.3 研究 | 2020-单细胞分析揭示葡萄膜黑色素瘤新的进化复杂性
    • 2.2.4 研究 | 2020-COVID-19病人支气管免疫细胞单细胞测序分析
    • 2.2.5 研究 | 2020-原汁原味读--单细胞肿瘤免疫图谱
    • 2.2.6 研究 | 2021-多发性骨髓瘤发展过程中肿瘤和免疫细胞的共同进化
    • 2.2.7 研究 | 2021-多个组织的成纤维细胞图谱
    • 2.2.8 研究 | 2021-多组学分析肺结核队列的记忆T细胞状态
    • 2.2.9 研究 | 2021-CancerSCEM: 人类癌症单细胞表达图谱数据库
    • 2.2.10 研究| 2021-单细胞转录组分析COVID-19重症患者肺泡巨噬细胞亚型
    • 2.2.11 研究 |2021-单细胞转录组揭示肺腺癌特有的肿瘤微环境
    • 2.2.12 研究 | 2021-单细胞转录组揭示乳头状甲状腺癌起始与发展
    • 2.2.13 研究 | 2021-解析食管鳞癌化疗病人的单细胞转录组
    • 2.2.14 研究 | 2021-单细胞水平看骨髓瘤的细胞状态和基因调控
    • 2.3.1 算法|2020-BatchBench比较scRNA批次矫正方法
    • 2.3.2 算法 | 2021-scPhere——用地球仪来展示降维结果
    • 2.3.3 算法 | 2021-单细胞差异分析方法评测
    • 2.3.4 算法 | 2021-细胞分群新方法——CNA(co-varying neighborhood analysis)
    • 2.3.5 工具 | 2018-iSEE:单细胞数据可视化辅助网页工具
    • 2.3.6 工具 | 2021-MACA: 一款自动注释细胞类型的工具
    • 2.3.7 工具 | 2021-一个很有想法的工具——Ikarus,想要在单细胞水平直接鉴定肿瘤细胞
  • 3 流程篇:分析框架
    • 3.1 质控
    • 3.2 归一化
    • 3.3 挑选表达量高变化基因
    • 3.4 降维
    • 3.5 聚类
    • 3.6 Marker/标记基因检测
    • 3.7 细胞类型注释
    • 3.8 批次效应处理
    • 3.9 多样本间差异分析
    • 3.10 检测Doublet
    • 3.11 细胞周期推断
    • 3.12 细胞轨迹推断
    • 3.13 与蛋白丰度信息结合
    • 3.14 处理大型数据
    • 3.15 不同R包数据的相互转换
  • 4 实战篇:活学活用
    • 4.1 实战一 | Smart-seq2 | 小鼠骨髓
    • 4.2 实战二 | STRT-Seq | 小鼠大脑
    • 4.3 实战三 | 10X | 未过滤的PBMC
    • 4.4 实战四 | 10X | 过滤后的PBMC
    • 4.5 实战五 | CEL-seq2 | 人胰腺细胞
    • 4.6 实战六 | CEL-seq | 人胰腺细胞
    • 4.7 实战七 | SMARTer | 人胰腺细胞
    • 4.8 实战八 | Smart-seq2 | 人胰腺细胞
    • 4.9 实战九 | 不同技术数据整合 | 人胰腺细胞
    • 4.10 实战十 | CEL-seq | 小鼠造血干细胞
    • 4.11 实战十一 | Smart-seq2 | 小鼠造血干细胞
    • 4.12 实战十二 | 10X | 小鼠嵌合体胚胎
    • 4.13 实战十三 | 10X | 小鼠乳腺上皮细胞
    • 4.14 | 实战十四 | 10X | HCA计划的38万骨髓细胞
  • 5 补充篇:开拓思路
    • 5.1 10X Genomics概述
      • 5.1.1 10X Genomics 问题集锦
    • 5.2 CellRanger篇
      • 5.2.1 CellRanger实战(一)数据下载
      • 5.2.2 CellRanger实战(二) 使用前注意事项
      • 5.2.3 CellRanger实战(三) 使用初探
      • 5.2.4 CellRanger实战(四)流程概览
      • 5.2.5 CellRanger实战(五) 理解count输出的结果
    • 5.3 Seurat的使用
      • 5.3.1 Seurat V3 | 实战之2700 PBMCs分析
      • 5.3.2 Seurat V3 | 如何改造Seurat包的DoHeatmap函数?
      • 5.3.3 scRNA的3大R包对比
      • 5.3.4 Seurat两种数据比较:integrated vs RNA assay
      • 5.3.5 seurat 的几种findmaker比较
    • 5.4 Monocle的使用
      • 5.4.1 Monocle V3实战
    • 5.5 多个数据集的整合
      • 5.5.1 使用Seurat的merge功能进行整合
      • 5.5.2 如何使用sctransform去除批次效应
由 GitBook 提供支持
在本页
  • 1 前言
  • 数据准备
  • 为了方便下面的批次矫正,先做几件事
  • 2 检查批次效应
  • 合并数据
  • 进行PCA
  • 再看聚类分群结果
  • 看t-SNE结果
  • 3 矫正批次效应之线性回归
  • 4 MNN矫正
  • 4.1 先了解一下这个方法
  • 4.2 还是应用到PBMC数据
  • 4.3 矫正后的检查
  • 5 更为细致的检查
  • 5.1 同一批次内clusters的比较
  • 5.2 利用marker基因检查
  • 6 矫正后数据的应用
  • TODO:对13.7 Using the corrected values继续补充

这有帮助吗?

  1. 3 流程篇:分析框架

3.8 批次效应处理

刘小泽写于2020.7.7

上一页3.7 细胞类型注释下一页3.9 多样本间差异分析

最后更新于4年前

这有帮助吗?

1 前言

大型的scRNA分析一般都需要整合多个批次的数据集,但每个批次的质量可能存在不可调和的差异,例如仪器的调整、试剂质量的差异。结果就导致了不同批次细胞中基因表达量的系统差异,称之为“批次效应”。之前在 这一部分中也介绍了一些批次效应的事情,它最大的隐患就是:如果它这个差异成为整个数据的主力军,那么真正的生物学差异就难以检测,对结果的解读变得更难。

因此现在出了一些多数据整合时批次效应处理工具,早期的方法 (; )是基于线性模型,它假设细胞群体的组成要么是已知的,要么是在批次间也是同类型的。但很明显这个假设有局限,我们不能保证这个假设成立,于是定制的算法产生(; ; ),它不再需要对细胞组成有一个先验知识,保证了对未知的探索过程。

数据准备

将使用两个不同批次的10X PBMC数据进行整合,它们都是从 这个数据包获得的,分别进行了前期的数据处理。各自前期处理的一个好处,就是可以根据各自的特点先过滤掉一些低质量数据(比如各自根据QC结果对低质量细胞进行过滤)

常规的前期处理步骤还是:质控=》归一化=》找HVGs=》降维=》聚类 只是注意下面的lapply的使用,多个数据的批量处理用列表很高效

# 数据下载(数据我已做好,链接在:https://share.weiyun.com/iE0VjVD0)
library(TENxPBMCData)
all.sce <- list(
    pbmc3k=TENxPBMCData('pbmc3k'),
    pbmc4k=TENxPBMCData('pbmc4k'),
    pbmc8k=TENxPBMCData('pbmc8k')
)

# 批量质控
library(scater)
stats <- high.mito <- list()
for (n in names(all.sce)) {
    current <- all.sce[[n]] #接下来的每步还是和之前一样
    is.mito <- grep("MT", rowData(current)$Symbol_TENx)
    stats[[n]] <- perCellQCMetrics(current, subsets=list(Mito=is.mito))
    high.mito[[n]] <- isOutlier(stats[[n]]$subsets_Mito_percent, type="higher")
    all.sce[[n]] <- current[,!high.mito[[n]]]
}

# 批量归一化
all.sce <- lapply(all.sce, logNormCounts)

# 批量找HVGs
library(scran)
all.dec <- lapply(all.sce, modelGeneVar)
all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1)
> length(all.hvgs[[1]])
[1] 1001
> length(all.hvgs[[2]])
[1] 1352

# 批量降维(使用了一个更高级的mapply:Apply a Function to Multiple List or Vector Arguments)
library(BiocSingular)
set.seed(10000)
all.sce <- mapply(FUN=runPCA, x=all.sce, subset_row=all.hvgs, 
    MoreArgs=list(ncomponents=25, BSPARAM=RandomParam()), 
    SIMPLIFY=FALSE)
set.seed(100000)
all.sce <- lapply(all.sce, runTSNE, dimred="PCA")

set.seed(1000000)
all.sce <- lapply(all.sce, runUMAP, dimred="PCA")

# 批量聚类
for (n in names(all.sce)) {
    g <- buildSNNGraph(all.sce[[n]], k=10, use.dimred='PCA')
    clust <- igraph::cluster_walktrap(g)$membership
    colLabels(all.sce[[n]])  <- factor(clust)
}

# 把数据提取出来
pbmc3k <- all.sce$pbmc3k
dec3k <- all.dec$pbmc3k
pbmc4k <- all.sce$pbmc4k
dec4k <- all.dec$pbmc4k

为了方便下面的批次矫正,先做几件事

1 选出两个数据集的基因交集,并各自取子集

这里两个数据都是Ensembl ID,因此这个过程很简单

universe <- intersect(rownames(pbmc3k), rownames(pbmc4k))
> length(rownames(pbmc3k));length(rownames(pbmc4k));length(universe)
[1] 32738
[1] 33694
[1] 31232

然后取子集

pbmc3k <- pbmc3k[universe,]
pbmc4k <- pbmc4k[universe,]

dec3k <- dec3k[universe,]
dec4k <- dec4k[universe,]

2 矫正批次间测序深度的影响

利用multiBatchNorm()函数,将两个批次放在一起,重新计算log-count,可以提高下游批次矫正的质量

还记得之前的size factor这个名词吗? 它是为了去除一个批次数据中的各个细胞之间的文库差异 而这里,是为了去除两个批次之间的差异

library(batchelor)
rescaled <- multiBatchNorm(pbmc3k, pbmc4k)
> class(rescaled)
[1] "list"
pbmc3k <- rescaled[[1]]
pbmc4k <- rescaled[[2]]

3 将之前单独挑出来的模型构建结果再次整合

当然还也可以把每个批次的HVGs取交集或并集,选多少基因这个没有成文规定。

这里只是因为scran提供了这么一个函数,并可以提供更多一些的基因,毕竟现在基因多一点没坏处,省的后面捉襟见肘

library(scran)
combined.dec <- combineVar(dec3k, dec4k)
chosen.hvgs <- combined.dec$bio > 0
sum(chosen.hvgs)
## [1] 13431

2 检查批次效应

在进行批次矫正之前,肯定要先看看有没有批次效应以及有多严重吧

最简单的方法是把两个数据合并,然后PCA+t-SNE

合并数据

如果简单的cbind(pbmc3k, pbmc4k),会报错

> cbind(pbmc3k, pbmc4k)
Error in FUN(X[[i]], ...) : 
  column(s) 'Symbol_TENx' in ‘mcols’ are duplicated and the data do not match
# 看似问题出现在mcols中的Symbol_TENx,看一眼

看样子的确存在不匹配的问题,第一行就暴露了,同一个Ensembl ID的Symbol就不同。这个问题先不深究。先看看两个数据的Ensembl 是否一致吧

> identical(mcols(pbmc3k)[,1],mcols(pbmc4k)[,1])
[1] TRUE
# 一致,那么就把它们的变成一样就好了,symbol的问题我们不考虑
rowData(pbmc3k) <- rowData(pbmc4k)
pbmc3k$batch <- "3k"
pbmc4k$batch <- "4k"
uncorrected <- cbind(pbmc3k, pbmc4k)

进行PCA

使用的HVGs也是合并两个数据集后得到的chosen.hvgs

library(scater)
set.seed(0010101010)
uncorrected <- runPCA(uncorrected, subset_row=chosen.hvgs,
    BSPARAM=BiocSingular::RandomParam())

再看聚类分群结果

如果说,这两个批次是真的重复,差异不大的话,那么分群的cluster应该在两个批次中包含差不多数量的细胞。但是下面发现基本clusters主要是仅仅包含了一个批次

library(scran)
snn.gr <- buildSNNGraph(uncorrected, use.dimred="PCA")
clusters <- igraph::cluster_walktrap(snn.gr)$membership
tab <- table(Cluster=clusters, Batch=uncorrected$batch)
tab
##        Batch
## Cluster   3k   4k
##      1     1  781
##      2     0 1309
##      3     0  535
##      4    14   51
##      5     0  605
##      6   489    0
##      7     0  184
##      8  1272    0
##      9     0  414
##      10  151    0
##      11    0   50
##      12  155    0
##      13    0   65
##      14    0   61
##      15    0   88
##      16   30    0
##      17  339    0
##      18  145    0
##      19   11    3
##      20    2   36

看t-SNE结果

set.seed(1111001)
uncorrected <- runTSNE(uncorrected, dimred="PCA")
plotTSNE(uncorrected, colour_by="batch")

虽然说t-SNE图上的点大小、远近说明不了问题,但是两个批次分的这么开,也提供了两个批次的证据

当然,还可以继续找证据。利用上一次的细胞类型注释,看看两个批次中的细胞类型,是不是差异很大。

3 矫正批次效应之线性回归

下面将使用batchelor 包中的rescaleBatches()函数进行处理

它也是对每个基因的log表达量进行了线性回归,并提高了一些运行性能。另外与removeBatchEffect()不同的是,rescaleBatches()保持了数据的稀疏性,而removeBatchEffect()会破坏稀疏性

library(batchelor)
rescaled <- rescaleBatches(pbmc3k, pbmc4k)
# 依然进行PCA、聚类
set.seed(1010101010) 
rescaled <- runPCA(rescaled, subset_row=chosen.hvgs, exprs_values="corrected")

snn.gr <- buildSNNGraph(rescaled, use.dimred="PCA")
clusters.resc <- igraph::cluster_walktrap(snn.gr)$membership
tab.resc <- table(Cluster=clusters.resc, Batch=rescaled$batch)
tab.resc
##        Batch
## Cluster    1    2
##      1   278  525
##      2    16   23
##      3   337  606
##      4    43  748
##      5   604  529
##      6    22   71
##      7   188   48
##      8    25   49
##      9   263    0
##      10  123  135
##      11   16   85
##      12   11   57
##      13  116    6
##      14  455 1035
##      15    6   31
##      16   89  187
##      17    3   36
##      18    3    8
##      19   11    3

效果还是有的,现在大部分clusters都包含了两个批次的细胞,但是依然有存在批次差异的cluster(看cluster9,在batch2中细胞数为0),表明处理有效果但不彻底 。影响效果的原因是:我们的数据违背了这个方法的假设

再做个图看看

rescaled <- runTSNE(rescaled, dimred="PCA")
rescaled$batch <- factor(rescaled$batch)
plotTSNE(rescaled, colour_by="batch")

注意

除了这个函数,还可以尝试更常规的regressBatches() ,只不过它的表现可能会更差,因为会丢失数据稀疏性

4 MNN矫正

4.1 先了解一下这个方法

试想batchA中有一个细胞a,然后想根据挑选的特征基因(如HVGs)去在batchB中鉴定与a细胞空间相近的细胞们。同样的,对batchB中的细胞b 重复这个过程,也鉴定它在batchA中的近邻。于是这个MNN就是:Mutual nearest neighbors ,是对一个批次的细胞找到它们在另一个批次中对应的最相邻细胞集合。

在进行批次矫正前,MNN找到的一对对细胞都是生物学状态相似的,因此如果再有不同,那么就姑且认为是外部的批次效应导致的,也就方便了去除。

与线性回归方法相比,MNN方法不会假设细胞群组成相同或者事先已知。MNN会自己学习细胞群的结构并进行估计。

4.2 还是应用到PBMC数据

batchelor 这个包除了线性模型,还提供了MNN 的函数fastMNN(),但这个函数与单纯的MNN算法还有不同。fastMNN()先进行PCA降维,以加速下面的检测细胞近邻。

set.seed(1000101001)
# 这里的d指的是前多少个主成分,k指近邻的数量(k增大会导致数据集合并更激进;一般不要设置太大,如果看到相同的细胞类型没有在批次之间充分合并,可以适当增加k)
mnn.out <- fastMNN(pbmc3k, pbmc4k, d=50, k=20, subset.row=chosen.hvgs,
    BSPARAM=BiocSingular::RandomParam(deferred=TRUE))
mnn.out
## class: SingleCellExperiment 
## dim: 13431 6791 
## metadata(2): merge.info pca.info
## assays(1): reconstructed
## rownames(13431): ENSG00000239945 ENSG00000228463 ... ENSG00000198695
##   ENSG00000198727
## rowData names(1): rotation
## colnames: NULL
## colData names(1): batch
## reducedDimNames(1): corrected
## altExpNames(0):

其中mnn.out的每列表示某个批次中的一个细胞(具体存储在batch),行表示chosen.hvgs中的基因

head(mnn.out$batch) 
## [1] 1 1 1 1 1 1

看到降维相关模块reducedDimNames中多了一个corrected,它是根据50个主成分得到的所有6791在低维空间的坐标

dim(reducedDim(mnn.out, "corrected"))
## [1] 6791   50

> reducedDim(mnn.out, "corrected")[1:3,1:3]
            [,1]       [,2]         [,3]
[1,] -0.12783777  0.0409469 -0.000116194
[2,] -0.03364614 -0.1466529  0.161690348
[3,] -0.09895631  0.1219669  0.010321992

看到表达量模块assays() 中多了一个reconstructed矩阵,是每个细胞中每个基因矫正后的表达量

> dim(assay(mnn.out, "reconstructed"))
[1] 13431  6791

> assay(mnn.out, "reconstructed")[1:4,1:4]
<4 x 4> matrix of class LowRankMatrix and type "double":
                         [,1]          [,2]
ENSG00000239945 -2.522191e-06 -1.851424e-06
ENSG00000228463 -6.626821e-04 -6.724341e-04
ENSG00000237094 -8.077231e-05 -8.038006e-05
ENSG00000229905  3.838135e-06  6.179994e-06
                         [,3]          [,4]
ENSG00000239945 -1.198984e-05 -3.192269e-06
ENSG00000228463 -4.820230e-04 -6.330560e-05
ENSG00000237094 -9.630608e-05 -6.855333e-05
ENSG00000229905  5.432122e-06 -2.118592e-05

4.3 矫正后的检查

这次函数已经做过PCA了,那就继续分群

library(scran)
snn.gr <- buildSNNGraph(mnn.out, use.dimred="corrected")
clusters.mnn <- igraph::cluster_walktrap(snn.gr)$membership
tab.mnn <- table(Cluster=clusters.mnn, Batch=mnn.out$batch)
tab.mnn
##        Batch
## Cluster    1    2
##      1   337  606
##      2   289  542
##      3   152  181
##      4    12    4
##      5   517  467
##      6    17   19
##      7   313  661
##      8   162  118
##      9    11   56
##      10  547 1083
##      11   17   59
##      12   16   58
##      13  144   93
##      14   67  191
##      15    4   36
##      16    4    8

可以看到这里的效果不错,没有哪个cluster在某个batch中细胞数为0

再作图看看

library(scater)
set.seed(0010101010)
mnn.out <- runTSNE(mnn.out, dimred="corrected")

mnn.out$batch <- factor(mnn.out$batch)
plotTSNE(mnn.out, colour_by="batch")

注意

fastMNN()函数其实也自带了检测数据

metadata(mnn.out)$merge.info$lost.var
##             [,1]        [,2]
## [1,] 0.006617087 0.003315395

它的含义是:the proportion of variance within each batch that is lost during MNN correction。如果这个值很大(比如>10%),说明有些”矫枉过正“,把一些内在的生物学差异也给矫正了。

5 更为细致的检查

5.1 同一批次内clusters的比较

两个批次数据混合后进行矫正,还是应该保留每个批次的特性。一个很不想看到的情况是:看上去细胞混合效果很好,以为矫正的效果不错,但没想到矫正过了,细胞中生物学异质性也被抹除了。

我们想在整合+矫正批次后,看到矫正后的一个cluster中既包含了原来的batch1信息,也包含batch2。当然表现的差异越大越好,体现在下图中就是:after1(即矫正后的cluster1)对应的PBMC3k 的clusters与PBMC 4k的clusters最好别重复,并且差别越大越好

library(pheatmap)
# 画出二者混合后的cluster与各自之前的cluster对比
# batch 1
tab <- table(paste("after", clusters.mnn[rescaled$batch==1]),
    paste("before", colLabels(pbmc3k)))
heat3k <- pheatmap(log10(tab+10), cluster_row=FALSE, cluster_col=FALSE,
    main="PBMC 3K comparison", silent=TRUE)

# batch 2
tab <- table(paste("after", clusters.mnn[rescaled$batch==2]),
    paste("before", colLabels(pbmc4k)))
heat4k <- pheatmap(log10(tab+10), cluster_row=FALSE, cluster_col=FALSE,
    main="PBMC 4K comparison", silent=TRUE)

gridExtra::grid.arrange(heat3k[[4]], heat4k[[4]])

就以after5来说,分别对应3k混合之前的cluster3、4、7;以及4k混合前的cluster1、10、3,差别蛮大,说明混合后依然保持了各自的个性

除了看图,还能通过数据体现——计算Rand index

这个指标是矫正后各个batch原来差异的保留度,这个数值越大越好(但如果某个批次矫正方法不给力,这个值也是大的,所以还是要慎重对待)

library(fossil)
# batch1
ri3k <- rand.index(as.integer(clusters.mnn[rescaled$batch==1]),
    as.integer(colLabels(pbmc3k)))
ri3k
## [1] 0.9335226

# batch2
ri4k <- rand.index(as.integer(clusters.mnn[rescaled$batch==2]),
    as.integer(colLabels(pbmc4k)))
ri4k
## [1] 0.9575746

5.2 利用marker基因检查

目的就是看数据整合+矫正后,原始batch的内部结构是否还完整

首先分别从原来的两个批次数据中寻找marker基因,因为marker基因的存在就保证了这个批次数据集中”生物学的有趣性“存在

stats3 <- pairwiseWilcox(pbmc3k, direction="up")
markers3 <- getTopMarkers(stats3[[1]], stats3[[2]], n=10)

stats4 <- pairwiseWilcox(pbmc4k, direction="up")
markers4 <- getTopMarkers(stats4[[1]], stats4[[2]], n=10)

然后用双方的marker基因合集,作为矫正的HVGs,看看算法会不会把这部分”有趣性“去除

marker.set <- unique(unlist(c(unlist(markers3), unlist(markers4))))
length(marker.set) 
## [1] 314

set.seed(1000110)
mnn.out2 <- fastMNN(pbmc3k, pbmc4k, subset.row=marker.set,
    BSPARAM=BiocSingular::RandomParam(deferred=TRUE))

最后作图看看

mnn.out2 <- runTSNE(mnn.out2, dimred="corrected")
gridExtra::grid.arrange(
    plotTSNE(mnn.out2[,mnn.out2$batch==1], colour_by=I(colLabels(pbmc3k))),
    plotTSNE(mnn.out2[,mnn.out2$batch==2], colour_by=I(colLabels(pbmc4k))),
    ncol=2
)

看样子根据这些marker基因进行的批次矫正,并没有去掉marker基因背后的生物学差异,因为两个批次数据还是能分出来群的,并且分群结果看似还不错

6 矫正后数据的应用

数据整合并矫正后,对合并后数据的cluster分析即可,不用再对每个batch的各个cluster单独分析。另一个好处是:batch整合后作为一整个数据集,细胞数量增加,相当于增加了群体结构的分辨率,方便下游的细胞层面分析(比如后面的轨迹推断)。

可能你还想利用合并后的数据进行基因层面的分析,例如差异分析marker基因鉴定。但一般不推荐,因为算法对多个批次进行合并时,不会保留每个基因表达的差异。如果要探索基因表达相关,最好还是要在未矫正数据基础上进行,并且还要把批次信息封锁掉(还记得之前findMarkers中的block参数吗?)

# 例如
m.out <- findMarkers(uncorrected, clusters.mnn, block=uncorrected$batch,
    direction="up", lfc=1, row.data=rowData(uncorrected)[,3,drop=FALSE])

TODO:对13.7 Using the corrected values继续补充

bulk mRNA转录组中常用的矫正批次效应方法就是线性回归,对每个基因表达量拟合一个线性模型。例如limma的removeBatchEffect() () 、sva的comBat() ( )。如果要使用这类方法,就需要假设:批次间的细胞组成相同。另外的一个假设是:批次效应的累积的,对于任何给定的基因,在不同亚群中经过任何因素诱导的表达变化倍数是相同的。(其实,从这两个假设就看出来,这个方法不适合我们的单细胞数据,但还是要继续了解下去)

Mutual nearest neighbors are pairs of cells from different batches that belong in each other s set of nearest neighbors.

3.2 归一化
Ritchie et al. 2015
Leek et al. 2012
Haghverdi et al. 2018
Butler et al. 2018
Lin et al. 2019
TENxPBMCData
Ritchie et al. 2015
Leek et al. 2012
Haghverdi et al. (2018)