单细胞交响乐
  • 前言:我与《单细胞交响乐》的缘分
  • 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 简单一点的试验
  • 2.1 取两个数据的交集子集
  • 2.2 数据整合后的归一化
  • 2.3 数据整合后找表达量高变化基因
  • 2.4 矫正批次效应
  • 3 更具挑战性的操作
  • 3.1 取四个数据的子集
  • 3.2 数据整合后的归一化
  • 3.3 数据整合后找表达量高变化基因
  • 3.4 矫正批次效应
  • 3.5 聚类
  • 作图
  • 3.6 对批次效应的检查
  • 3.7 进行一次更严格的批次矫正

这有帮助吗?

  1. 4 实战篇:活学活用

4.9 实战九 | 不同技术数据整合 | 人胰腺细胞

刘小泽写于2020.7.21

上一页4.8 实战八 | Smart-seq2 | 人胰腺细胞下一页4.10 实战十 | CEL-seq | 小鼠造血干细胞

最后更新于4年前

这有帮助吗?

1 前言

2016年很多研究都对人类胰腺细胞的scRNA很感兴趣,也因此发表了很多文章:(; ; ; )。这些不同作者不同技术手段得到的数据,也给数据整合带来了不小的挑战。相比于之前的PBMC数据整合,这里更为复杂,因为包含的建库、测序方法多样,供体的类型、数量更是不一致。

2 简单一点的试验

首先拿技术相似的两套数据来做:分别是 ; ,采用了CEL-seq 和 CEL-seq2

load('final.sce.grun.Rdata')
load('final.sce.muraro.RData')
final.sce.grun
# class: SingleCellExperiment 
# dim: 17548 1063 
# metadata(0):
#   assays(2): counts logcounts
# rownames(17548): ENSG00000268895 ENSG00000121410 ...
# ENSG00000074755 ENSG00000036549
# rowData names(2): symbol chr
# colnames(1063): D2ex_1 D2ex_2 ... D17TGFB_94
# D17TGFB_95
# colData names(3): donor sample sizeFactor
# reducedDimNames(0):
#   altExpNames(1): ERCC

final.sce.muraro
# class: SingleCellExperiment 
# dim: 16940 2299 
# metadata(0):
#   assays(2): counts logcounts
# rownames(16940): ENSG00000268895 ENSG00000121410 ...
# ENSG00000159840 ENSG00000074755
# rowData names(2): symbol chr
# colnames(2299): D28-1_1 D28-1_2 ... D30-8_93
# D30-8_94
# colData names(4): label donor plate sizeFactor
# reducedDimNames(0):
#   altExpNames(1): ERCC

2.1 取两个数据的交集子集

首先获得交集基因

universe <- intersect(rownames(final.sce.grun), rownames(final.sce.muraro))

> nrow(final.sce.grun);nrow(final.sce.muraro);length(universe)
[1] 17548
[1] 16940
[1] 15974

对数据集取子集

sce.grun2 <- final.sce.grun[universe,]
sce.muraro2 <- final.sce.muraro[universe,]

既然是经过处理后的数据,那么就略过了之前介绍的质控步骤

2.2 数据整合后的归一化

首先测序深度导致的文库大小差异是批次效应的一个重要来源,因此可以先对不同的批次进行文库矫正。会以文库最小的批次为基准,对其他批次进行文库归一化。最后返回一个列表

使用一个归一化函数:multiBatchNorm ,它应用的就是最简单的library size normalization归一化方法: Perform scaling normalization within each batch to provide comparable results to the lowest-coverage batch.

既然是要处理文库差异,那就先看看各自原本的文库大小

summary(colSums(logcounts(sce.muraro2)))
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 4433    8041    8680    8558    9193   10594 
summary(colSums(logcounts(sce.grun2)))
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 2347    4141    4565    4538    4969    5992

然后进行处理,可以看看前后的变化,就明白了这个函数做了什么事情

library(batchelor)
normed.pancreas <- multiBatchNorm(sce.grun2, sce.muraro2)
sce.grun3 <- normed.pancreas[[1]]
sce.muraro3 <- normed.pancreas[[2]]

summary(colSums(logcounts(sce.muraro3)))
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 3072    4735    4994    4954    5204    5864 
summary(colSums(logcounts(sce.grun3)))
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 2347    4141    4565    4538    4969    5992

看到混合以后的sce.muraro3相对于之前独立的sce.muraro2的变化了吧

2.3 数据整合后找表达量高变化基因

首先对表达量变化模型取子集

## 原来的sce.grun模型
library(scran)
block1 <- paste0(final.sce.grun$sample, "_", final.sce.grun$donor)
dec.grun <- modelGeneVarWithSpikes(final.sce.grun, spikes="ERCC", block=block1)
# 取子集
dec.grun2 <- dec.grun[universe,]

## 原来的sce.muraro模型
block2 <- paste0(final.sce.muraro$plate, "_", final.sce.muraro$donor)
dec.muraro <- modelGeneVarWithSpikes(final.sce.muraro, "ERCC", block=block2)
# 取子集
dec.muraro2 <- dec.muraro[universe,]

之后组合两组的结果

使用combineVar ,它的作用是:

Combine the results of multiple variance decompositions, usually generated for the same genes across separate batches of cells.

library(scran)
combined.pan <- combineVar(dec.grun2, dec.muraro2)
# 把更有可能代表生物差异的基因选出来,用于下游的PCA和聚类
chosen.genes <- rownames(combined.pan)[combined.pan$bio > 0]

2.4 矫正批次效应

先来看看基于线性回归的rescaleBatches()

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

library(scater)
rescaled.pancreas <- rescaleBatches(sce.grun2, sce.muraro2)

set.seed(100101)
rescaled.pancreas <- runPCA(rescaled.pancreas, subset_row=chosen.genes,
    exprs_values="corrected")

rescaled.pancreas <- runTSNE(rescaled.pancreas, dimred="PCA")
plotTSNE(rescaled.pancreas, colour_by="batch")

不过结果并不尽如人意,两个批次还是分得很开,表明处理有效果但不彻底 。影响效果的原因是:我们的数据违背了这个方法的假设

再来使用fastMNN()

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

它和fastMNN()原理类似,但速度会慢很多,总之它们的不同可以概括为:

For scRNA-seq data, fastMNN() tends to be both faster and better at achieving a satisfactory merge. mnnCorrect() is mainly provided here for posterity’s sake, though it is more robust than fastMNN() to certain violations of the orthogonality assumptions.

set.seed(1011011)
mnn.pancreas <- fastMNN(sce.grun2, sce.muraro2, subset.row=chosen.genes)

snn.gr <- buildSNNGraph(mnn.pancreas, use.dimred="corrected")
clusters <- igraph::cluster_walktrap(snn.gr)$membership
tab <- table(Cluster=clusters, Batch=mnn.pancreas$batch)
tab
##        Batch
## Cluster   1   2
##      1  239 281
##      2  312 257
##      3  200 837
##      4   56 193
##      5   37   1
##      6   24 108
##      7  109 391
##      8   63  80
##      9   18 115
##      10   0  17
##      11   5  19

再做个图

mnn.pancreas <- runTSNE(mnn.pancreas, dimred="corrected")
plotTSNE(mnn.pancreas, colour_by="batch")

效果有明显改进

3 更具挑战性的操作

前面是将一个CEL-seq和一个CEL-seq2数据整合,总的说来还不是很复杂

rm(list = ls())
load('final.sce.grun.RData')
load('final.sce.muraro.RData')
load('final.sce.lawlor.RData')
load('final.sce.seger.RData')
sce.grun=final.sce.grun
sce.muraro=final.sce.muraro

sce.grun
# class: SingleCellExperiment 
# dim: 17548 1063 
# metadata(0):
#   assays(2): counts logcounts
# rownames(17548): ENSG00000268895 ENSG00000121410 ...
# ENSG00000074755 ENSG00000036549
# rowData names(2): symbol chr
# colnames(1063): D2ex_1 D2ex_2 ... D17TGFB_94
# D17TGFB_95
# colData names(3): donor sample sizeFactor
# reducedDimNames(0):
#   altExpNames(1): ERCC

sce.muraro
# class: SingleCellExperiment 
# dim: 16940 2299 
# metadata(0):
#   assays(2): counts logcounts
# rownames(16940): ENSG00000268895 ENSG00000121410 ...
# ENSG00000159840 ENSG00000074755
# rowData names(2): symbol chr
# colnames(2299): D28-1_1 D28-1_2 ... D30-8_93
# D30-8_94
# colData names(4): label donor plate sizeFactor
# reducedDimNames(0):
#   altExpNames(1): ERCC

sce.lawlor
# class: SingleCellExperiment 
# dim: 26616 604 
# metadata(0):
#   assays(2): counts logcounts
# rownames(26616): ENSG00000229483 ENSG00000232849 ...
# ENSG00000251576 ENSG00000082898
# rowData names(2): SYMBOL SEQNAME
# colnames(604): 10th_C11_S96 10th_C13_S61 ...
# 9th-C96_S81 9th-C9_S13
# colData names(9): title age ... Sex sizeFactor
# reducedDimNames(0):
#   altExpNames(0):

sce.seger
# class: SingleCellExperiment 
# dim: 25454 2090 
# metadata(0):
#   assays(2): counts logcounts
# rownames(25454): ENSG00000118473 ENSG00000142920 ...
# ENSG00000278306 eGFP
# rowData names(2): symbol refseq
# colnames(2090): HP1502401_H13 HP1502401_J14 ...
# HP1526901T2D_N8 HP1526901T2D_A8
# colData names(4): CellType Donor Quality sizeFactor
# reducedDimNames(0):
#   altExpNames(1): ERCC

3.1 取四个数据的子集

首先获得交集基因

all.sce <- list(Grun=sce.grun, Muraro=sce.muraro, 
    Lawlor=sce.lawlor, Seger=sce.seger)
universe <- Reduce(intersect, lapply(all.sce, rownames))

nrow(sce.grun);nrow(sce.muraro);nrow(sce.lawlor);nrow(sce.seger);length(universe)
# [1] 17548
# [1] 16940
# [1] 26616
# [1] 25454
# [1] 15231 #共有基因

再分别取子集

all.sce <- lapply(all.sce, "[", i=universe,)

3.2 数据整合后的归一化

normed.pancreas <- do.call(multiBatchNorm, all.sce)

3.3 数据整合后找表达量高变化基因

首先获得各个数据集的表达量变化模型

library(scran)
# sce.grun
block1 <- paste0(sce.grun$sample, "_", sce.grun$donor)
dec.grun <- modelGeneVarWithSpikes(sce.grun, spikes="ERCC", block=block1)

# sce.muraro
block2 <- paste0(sce.muraro$sample, "_", sce.muraro$donor)
dec.muraro <- modelGeneVarWithSpikes(sce.muraro, spikes="ERCC", block=block2)

# sce.lawlor:没有ERCC信息,就用modelGeneVar()
dec.lawlor <- modelGeneVar(sce.lawlor, block=sce.lawlor$`islet unos id`)

# sce.seger
for.hvg <- sce.seger[,librarySizeFactors(altExp(sce.seger)) > 0
                     & sce.seger$Donor!="AZ"]
dec.seger <- modelGeneVarWithSpikes(for.hvg, "ERCC", block=for.hvg$Donor)

# 整合起来
all.dec <- list(Grun=dec.grun, Muraro=dec.muraro, 
                Lawlor=dec.lawlor, Seger=dec.seger)

再取子集

all.dec <- lapply(all.dec, "[", i=universe,)

再组合四组的结果

combined.pan <- do.call(combineVar, all.dec)
# 把更有可能代表生物差异的基因选出来,用于下游的PCA和聚类
chosen.genes <- rownames(combined.pan)[combined.pan$bio > 0]

3.4 矫正批次效应

fastMNN也包含了PCA的操作

set.seed(1011110)
mnn.pancreas <- fastMNN(normed.pancreas)

3.5 聚类

snn.gr <- buildSNNGraph(mnn.pancreas, use.dimred="corrected", k=20)
clusters <- igraph::cluster_walktrap(snn.gr)$membership
clusters <- factor(clusters)
tab <- table(Cluster=clusters, Batch=mnn.pancreas$batch)

tab
##        Batch
## Cluster Grun Lawlor Muraro Seger
##      1    77     33    677   211
##      2   311     26    254   383
##      3   104    244    390   180
##      4   225     16    246   138
##      5   125    203    158   108
##      6    56     17    196   109
##      7     0      0      0    43
##      8     0      0      0    42
##      9     0      2     17     4
##      10   69     12     82    16
##      11   24     19    108    55
##      12    0      0      0    50
##      13   27      0      1     0
##      14   22      6     34    49
##      15   18     18    117   157
##      16    0      0      0   208
##      17    0      0      0    26
##      18    0      0      0   108
##      19    0      0      0   186
##      20    5      8     19    17

看到一个批次中都包含了很多clusters,说明数据整合的效果还不错,批次效应没有很强;当然有些clusters只显示在了seger批次中(比如7、8、12、16、17、18、19),那究竟这些clusters到底是不是seger数据特有的细胞类型呢?这个还有待考证

作图

注意其中使用到了一个很有趣的函数I() ,简单的一个字母,它是base包里的函数

因为我们是根据mnn.pancreas进行作图的,但clusters这个向量是根据mnn.pancreas创建的,但又不直接存在于mnn.pancreas(不像batch一样存在于mnn.pancreas中)。

因此要从外部把它导入到作图函数中,就可以用这个I()

mnn.pancreas <- runTSNE(mnn.pancreas, dimred="corrected")
gridExtra::grid.arrange(
    plotTSNE(mnn.pancreas, colour_by="batch", text_by=I(clusters)),
    plotTSNE(mnn.pancreas, colour_by=I(clusters), text_by=I(clusters)),
    ncol=2
)

上面我们粗略根据四个数据集看了下批次效应,发现Seger这个数据还有点特殊,因为很多clusters只存在于Seger中

批次效应的来源除了表面上的4个数据集整合,还有一个重点考虑对象是:供体的种类

3.6 对批次效应的检查

看一下来自各个数据中供体的批次效应

首先检查一下各个数据的供体信息,看到Seger最多,因此它的风险也最大。但这个怀疑到底对不对,还要做个图看看

seger.donors <- donors
seger.donors[mnn.pancreas$batch!="Seger"] <- NA

grun.donors <- donors
grun.donors[mnn.pancreas$batch!="Grun"] <- NA

lawlor.donors <- donors
lawlor.donors[mnn.pancreas$batch!="Lawlor"] <- NA

muraro.donors <- donors
muraro.donors[mnn.pancreas$batch!="Muraro"] <- NA

gridExtra::grid.arrange(
  plotTSNE(mnn.pancreas, colour_by=I(muraro.donors))+ ggtitle('muraro.donors'),
  plotTSNE(mnn.pancreas, colour_by=I(lawlor.donors))+ ggtitle('lawlor.donors'),
  plotTSNE(mnn.pancreas, colour_by=I(grun.donors))+ ggtitle('grun.donors'),
  plotTSNE(mnn.pancreas, colour_by=I(seger.donors))+ ggtitle('seger.donors'),
  ncol=2
)

看到图中Seger的供体各个细胞最为分散,因此它的供体批次效应是最强的

虽说供体也是生物信息,但它对于后续的细胞类型注释没有直接的帮助,相反还会产生混淆的作用(比如前面看到很多clusters只存在于Seger中,说不定就是由于Seger的供体批次效应导致的)

因此,一个更为谨慎的操作是:除了去除数据集之间的批次效应以外,还要将每个数据内部的供体信息作为另一个批次效应,分别处理掉

3.7 进行一次更严格的批次矫正

将原来的4个分开的数据聚合在一起,使用noCorrect 进行简单的聚合

它的含义是:This function is effectively equivalent to cbinding the matrices together without any correction.

combined <- noCorrect(normed.pancreas)
assayNames(combined) <- "logcounts"
combined$donor <- donors

对这个数据进行批次矫正,但首先要把数据批次和供体批次分开

donors.per.batch <- split(combined$donor, combined$batch)

# 获得每个数据批次下的供体批次
donors.per.batch <- lapply(donors.per.batch, unique)
donors.per.batch
## $Grun
## [1] "D2"  "D3"  "D7"  "D10" "D17"
## 
## $Lawlor
## [1] "ACIW009"  "ACJV399"  "ACCG268"  "ACCR015A" "ACEK420A" "ACEL337"  "ACHY057" 
## [8] "ACIB065" 
## 
## $Muraro
## [1] "D28" "D29" "D31" "D30"
## 
## $Seger
##  [1] "HP1502401"    "HP1504101T2D" "AZ"           "HP1508501T2D" "HP1506401"   
##  [6] "HP1507101"    "HP1509101"    "HP1504901"    "HP1525301T2D" "HP1526901T2D"

依然是使用fastMNN

set.seed(1010100)
# batch信息使用全部的供体
# 增加一步指定weights:可以理解为哪些供体属于哪个数据集
multiout <- fastMNN(combined, batch=combined$donor, 
    subset.row=chosen.genes, weights=donors.per.batch)

# 将两大批次信息记录在新的矫正结果multiout中
multiout$dataset <- combined$batch
multiout$donor <- multiout$batch

检查一下聚类结果

从下面的结果中可以看到单独属于Seger的cluster没有了

library(scater)
g <- buildSNNGraph(multiout, use.dimred=1, k=20)
clusters <- igraph::cluster_walktrap(g)$membership
tab <- table(clusters, multiout$dataset)
tab
##         
## clusters Grun Lawlor Muraro Seger
##       1   246     20    278   187
##       2   200    239    835   862
##       3   171    254    473   294
##       4   315     27    260   387
##       5    57     17    193   108
##       6    24     18    107    55
##       7    26      0      0     0
##       8     5      9     19    17
##       9    19     19    118   176
##       10    0      1     16     4

作图结果也发现数据混合更理想了

multiout <- runTSNE(multiout, dimred="corrected")
gridExtra::grid.arrange(
    plotTSNE(multiout, colour_by="dataset", text_by=I(clusters)),
    plotTSNE(multiout, colour_by=I(seger.donors)),
    ncol=2
)

最后,和已发表的细胞类型做对比

由于这些数据都已发表,数据集中也包含了最后作者注释的细胞类型(sce.grun除外)

因此,可以将我们自己整合后又矫正的分群,与发表的细胞分群对比,来说明批次处理质量

# 获得已发表的细胞类型信息
proposed <- c(rep(NA, ncol(sce.grun)), 
    sce.muraro$label,
    sce.lawlor$`cell type`,
    sce.seger$CellType)

看到其中大小写参差不齐,可以全变成小写

proposed <- tolower(proposed)
# 并根据原文章修改一下细胞类型名称
proposed[proposed=="gamma/pp"] <- "gamma"
proposed[proposed=="pp"] <- "gamma"
proposed[proposed=="duct"] <- "ductal"
proposed[proposed=="psc"] <- "stellate"

最后检查一下

table(proposed, clusters)

看到,我们处理完批次效应后的分群结果可以比较好的匹配到真实的细胞类型,因此说明了这里使用的批次矫正的方法的力度刚刚好

之前在: 中介绍过:

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

可能之前听过:mnnCorrect()这个方法,它是提出来的,之前也介绍过:

但这里,会再加上两个数据,分别来自 和 的数据进行整合,这样四个数据就会包括不同的技术、不同的UMI、不同的表达量、更加不同的供体

Muraro et al. 2016
Grun et al. 2016
Lawlor et al. 2017
Segerstolpe et al. 2016
Muraro et al. 2016
Grun et al. 2016
sce.grun分享数据下载
sce.muraro分享数据链接
3.8 批次效应处理
Ritchie et al. 2015
Leek et al. 2012
Haghverdi et al. (2018)
单细胞转录组数据校正批次效应实战
Lawlor et al. 2017
Segerstolpe et al. 2016
sce.lawlor数据分享链接
sce.seger数据分享链接