> For the complete documentation index, see [llms.txt](https://jieandze1314.osca.top/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://jieandze1314.osca.top/04/04-9.md).

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

## 1 前言

2016年很多研究都对人类胰腺细胞的scRNA很感兴趣，也因此发表了很多文章：([Muraro et al. 2016](https://pubmed.ncbi.nlm.nih.gov/27693023/); [Grun et al. 2016](https://pubmed.ncbi.nlm.nih.gov/27345837/); [Lawlor et al. 2017](https://pubmed.ncbi.nlm.nih.gov/27864352/); [Segerstolpe et al. 2016](https://pubmed.ncbi.nlm.nih.gov/27667667/))。这些不同作者不同技术手段得到的数据，也给数据整合带来了不小的挑战。相比于之前的PBMC数据整合，这里更为复杂，因为包含的建库、测序方法多样，供体的类型、数量更是不一致。

## 2 简单一点的试验

首先拿技术相似的两套数据来做：分别是 [Muraro et al. 2016](https://pubmed.ncbi.nlm.nih.gov/27693023/); [Grun et al. 2016](https://pubmed.ncbi.nlm.nih.gov/27345837/)，采用了CEL-seq 和 CEL-seq2

> [sce.grun分享数据下载](<https://share.weiyun.com/iO5BCEc9 >)\
> [sce.muraro分享数据链接](https://share.weiyun.com/gvnArtAf)

```r
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 取两个数据的交集子集**

**首先获得交集基因**

```r
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
```

**对数据集取子集**

```r
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.

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

```r
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
```

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

```r
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 数据整合后**找表达量高变化基因

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

```r
## 原来的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.

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

### **2.4 矫正批次效应**

之前在：[3.8 批次效应处理](https://jieandze1314.osca.top/03/03-8) 中介绍过：

> bulk mRNA转录组中常用的矫正批次效应方法就是线性回归，对每个基因表达量拟合一个线性模型。例如limma的`removeBatchEffect()` ([Ritchie et al. 2015](https://pubmed.ncbi.nlm.nih.gov/25605792/)) 、sva的`comBat()` ( [Leek et al. 2012](https://pubmed.ncbi.nlm.nih.gov/22257669/))。如果要使用这类方法，就需要假设：批次间的细胞组成相同。另外的一个假设是：批次效应的累积的，对于任何给定的基因，在不同亚群中经过任何因素诱导的表达变化倍数是相同的。（其实，**从这两个假设就看出来，这个方法不适合我们的单细胞数据**，但还是要继续了解下去）

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

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

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

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

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-21-032715.png)

**再来使用fastMNN()**

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

可能之前听过：`mnnCorrect()`这个方法，它是[Haghverdi et al. (2018) ](https://pubmed.ncbi.nlm.nih.gov/29608177/)提出来的，之前也介绍过：[单细胞转录组数据校正批次效应实战](https://www.jianshu.com/p/b7f6a5efed85)

它和`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.

```r
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
```

再做个图

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

效果有明显改进

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-21-033244.png)

## 3 更具挑战性的操作

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

但这里，会再加上两个数据，分别来自[Lawlor et al. 2017](https://pubmed.ncbi.nlm.nih.gov/27864352/) 和 [Segerstolpe et al. 2016](https://pubmed.ncbi.nlm.nih.gov/27667667/) 的数据进行整合，这样四个数据就会包括不同的技术、不同的UMI、不同的表达量、更加不同的供体

> [sce.lawlor数据分享链接](<https://share.weiyun.com/mic3m90k >)\
> [sce.seger数据分享链接](https://share.weiyun.com/7w2vBGdC)

```r
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 取四个数据的子集**

**首先获得交集基因**

```r
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 #共有基因
```

**再分别取子集**

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

### **3.2 数据整合后的归一化**

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

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

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

```r
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)
```

**再取子集**

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

**再组合四组的结果**

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

### **3.4 矫正批次效应**

`fastMNN`也包含了PCA的操作

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

### **3.5 聚类**

```r
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()`

```r
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
)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-21-045902.png)

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

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

### **3.6 对批次效应的检查**

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

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

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-21-084814.png)

```r
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的供体各个细胞最为分散，因此它的供体批次效应是最强的

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-21-085745.png)

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

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

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

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

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

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

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-21-092211.png)

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

```r
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**

```r
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没有了

```r
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
```

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

```r
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
)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-21-092759.png)

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

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

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

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

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-21-093233.png)

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

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

最后检查一下

```r
table(proposed, clusters)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-21-093723.png)

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