# 4.4 实战四 | 10X | 过滤后的PBMC

## 1 前言

这次使用的数据是来自[Zheng et al. 2017](https://pubmed.ncbi.nlm.nih.gov/28091601/) 的三个PBMC数据，而且这个数据是经过前期过滤的

### **准备数据**

```r
library(TENxPBMCData)
all.sce <- list(
    pbmc3k=TENxPBMCData('pbmc3k'),
    pbmc4k=TENxPBMCData('pbmc4k'),
    pbmc8k=TENxPBMCData('pbmc8k')
)

all.sce
# $pbmc3k
# class: SingleCellExperiment 
# dim: 32738 2700 
# metadata(0):
#   assays(1): counts
# rownames(32738): ENSG00000243485 ENSG00000237613 ...
# ENSG00000215616 ENSG00000215611
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames: NULL
# colData names(11): Sample Barcode ... Individual
# Date_published
# reducedDimNames(0):
#   altExpNames(0):
#   
#   $pbmc4k
# class: SingleCellExperiment 
# dim: 33694 4340 
# metadata(0):
#   assays(1): counts
# rownames(33694): ENSG00000243485 ENSG00000237613 ...
# ENSG00000277475 ENSG00000268674
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames: NULL
# colData names(11): Sample Barcode ... Individual
# Date_published
# reducedDimNames(0):
#   altExpNames(0):
#   
#   $pbmc8k
# class: SingleCellExperiment 
# dim: 33694 8381 
# metadata(0):
#   assays(1): counts
# rownames(33694): ENSG00000243485 ENSG00000237613 ...
# ENSG00000277475 ENSG00000268674
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames: NULL
# colData names(11): Sample Barcode ... Individual
# Date_published
# reducedDimNames(0):
#   altExpNames(0):
```

**多个数据集的批量处理，重点就是列表list和for循环的熟练使用，还有相关的apply家族函数**。而且每个结果数据也要对应放在一个新列表中，比如下面质控使用的`stats <- high.mito <- list()` ，就是新建了两个空列表，然后把结果放进去

## 2 批量质控

### 数据备份

把unfiltered数据主要用在质控的探索上

```r
unfiltered <- all.sce
```

还是先通过线粒体含量计算质控结果，然后根据这个结果进行过滤，一个for循环搞定

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

### **看一下根据线粒体过滤的结果**

```r
> lapply(high.mito, summary)
$pbmc3k
   Mode   FALSE    TRUE 
logical    2609      91 

$pbmc4k
   Mode   FALSE    TRUE 
logical    4182     158 

$pbmc8k
   Mode   FALSE    TRUE 
logical    8157     224
```

### **批量作图（也是把作图结果放进list，方便后期批量导出）**

```r
qcplots <- list()

for (n in names(all.sce)) {
  current <- unfiltered[[n]]
  colData(current) <- cbind(colData(current), stats[[n]])
  current$discard <- high.mito[[n]]
  qcplots[[n]] <- plotColData(current, x="sum", y="subsets_Mito_percent",
                              colour_by="discard") + scale_x_log10()
}
# do.call也是list处理中的常用函数
do.call(gridExtra::grid.arrange, c(qcplots, ncol=3))
```

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

## 3 批量归一化

这里使用的是最简单的方法`logNormCounts()`，就是用某个细胞中每个基因或spike-in转录本的表达量除以这个细胞计算的size factor，最后还进行一个log转换，得到一个新的矩阵：`logcounts` 【不过这个名字并不是很贴切，只是因为拼写简单，真实应该是：log-transformed normalized expression values。而不是单纯字面意思取个log】

```r
all.sce <- lapply(all.sce, logNormCounts)
```

**看到这里，可能会想，为什么没计算size factor就直接进行了logNormCounts？**

前面提到的操作，一般都是：

```r
# 常规方法
lib.sf.zeisel <- librarySizeFactors(sce.zeisel)
lib.sf.zeisel <- logNormCounts(lib.sf.zeisel)
# 去卷积方法
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
sce.pbmc <- logNormCounts(sce.pbmc)
```

**看一下`logNormCounts`的帮助文档就能明白了，逻辑很清楚：**

* 函数默认的参数是：`size_factors=NULL`，如果没有计算size factor更新给函数，那么函数会执行`normalizeCounts` 的操作
* 再来看`normalizeCounts`会有什么操作：如果没有提供size factor，它会根据数据类型去自己计算size factor
  * 对于count矩阵和SummarizedExperiment数据类型，会通过`librarySizeFactors`计算
  * 对于SingleCellExperiment这种数据类型，它会首先在数据中寻找size factor是否存在，如果找不到，也是会使用`librarySizeFactors`

也就是说，**如果我们不提前计算，就会自动帮我们用最简单的`librarySizeFactors`计算，并添加到我们的数据中**。正是因为我们只需要最简单的方法，所以才可以不提供。如果要使用去卷积方法，还是要自己先计算好

**最后看下结果**

```r
> lapply(all.sce, function(x) summary(sizeFactors(x)))
$pbmc3k
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2338  0.7478  0.9262  1.0000  1.1571  6.6042 

$pbmc4k
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3155  0.7109  0.8903  1.0000  1.1272 11.0267 

$pbmc8k
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2963  0.7043  0.8772  1.0000  1.1177  6.7942
```

## 4 批量找表达量高变化基因

同样也是使用了最简单的计算方法

```r
library(scran)
all.dec <- lapply(all.sce, modelGeneVar)
all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1)
```

作图

```r
par(mfrow=c(1,3))
for (n in names(all.dec)) {
    curdec <- all.dec[[n]]
    plot(curdec$mean, curdec$total, pch=16, cex=0.5, main=n,
        xlab="Mean of log-expression", ylab="Variance of log-expression")
    curfit <- metadata(curdec)
  # 可视化一条线（下图的蓝线），这条线指所有的基因都会存在的一种偏差
    curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
}
```

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

* 每个点表示一个基因
* 图中蓝线指的是：技术因素导致的偏差
* 纵坐标表示总偏差：它等于技术偏差+生物因素偏差

因此，要衡量一个基因的生物因素偏差大小，就看对应的纵坐标减去对应的蓝线的值

## 5 批量降维

这里将每个PBMC数据单独进行降维，而并没有把它们混合起来再分析

关于降维方法，这里选择是与PCA近似的**SVD算法（singular value decomposition，奇异值分解）**，scater或scran都可以直接通过函数计算SVD，利用参数`BSPARAM=`传递一个`BiocSingularParam`对象到`runPCA`中

* [SVD与PCA的解释](https://www.cnblogs.com/bjwu/p/9280492.html)
* SVD是一种矩阵分解方法，相当于因式分解，目的纯粹就是将一个矩阵拆分成多个矩阵相乘的形式
* **对于稀疏矩阵来说，SVD算法更适用**，这样对于大数据来说节省了很大空间

```r
library(BiocSingular)
set.seed(10000)
# 这里的runPCA是一个降维的函数名字，它可以用本身的PCA算法，还可以用SVD算法
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")
```

这里使用SVD其实还是为了帮助更好地进行PCA，支持4种方式：

* ExactParam: exact SVD with runExactSVD.
* IrlbaParam: approximate SVD with irlba via runIrlbaSVD.
* RandomParam: approximate SVD with rsvd via runRandomSVD.
* FastAutoParam: fast approximate SVD, chosen based on the matrix representation.

## 6 批量聚类

使用基于图形的聚类，最基础的想法是：我们首先构建一个图，其中每个节点都是一个细胞，它与高维空间中最邻近的细胞相连。连线基于细胞之间的相似性计算权重，权重越高，表示细胞间关系更密切。如果一群细胞之间的权重高于另一群细胞，那么这一群细胞就被当做一个群体 “communiity”。

```r
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)
}
# 看看各自分了多少群
lapply(all.sce, function(x) table(colLabels(x)))
## $pbmc3k
## 
##   1   2   3   4   5   6   7   8   9  10 
## 487 154 603 514  31 150 179 333 147  11 
## 
## $pbmc4k
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13 
##  497  185  569  786  373  232   44 1023   77  218   88   54   36 
## 
## $pbmc8k
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
## 1004  759 1073 1543  367  150  201 2067   59  154  244   67   76  285   20   15 
##   17   18 
##   64    9
```

作图

```r
# 还是先新建一个空列表，方便保存数据
all.tsne <- list()

for (n in names(all.sce)) {
    all.tsne[[n]] <- plotTSNE(all.sce[[n]], colour_by="label") + ggtitle(n)
}
do.call(gridExtra::grid.arrange, c(all.tsne, list(ncol=3)))
```

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

## 7 数据整合

> 前面只是批量进行了各个数据集的分析，现在要把它们整合起来再分析一下

### **找共有基因**

```r
# 先看看各自多少基因
> lapply(all.sce, function(x) length(rownames(x)))
$pbmc3k
[1] 32738

$pbmc4k
[1] 33694

$pbmc8k
[1] 33694

# 找共同基因
universe <- Reduce(intersect, lapply(all.sce, rownames))
> length(universe)
[1] 31232
```

### **对每个数据批量取子集**

```r
# 这个操作可以记下来，lapply批量取子集
all.sce2 <- lapply(all.sce, "[", i=universe,)

> lapply(all.sce2, dim)
$pbmc3k
[1] 31232  2609

$pbmc4k
[1] 31232  4182

$pbmc8k
[1] 31232  8157

# 同样的，对找高变异基因的结果取子集
all.dec2 <- lapply(all.dec, "[", i=universe,)
```

### **把三个数据当做一个数据的三个批次，重新进行归一化**

```r
library(batchelor)
normed.sce <- do.call(multiBatchNorm, all.sce2)
```

### **根据重新归一化的结果，再次找HVGs**

这次是把3个批次放在一起再找的

```r
combined.dec <- do.call(combineVar, all.dec2)
combined.hvg <- getTopHVGs(combined.dec, n=5000)
```

### **对一个大数据进行降维**

> 结合：[单细胞交响乐10-数据集整合后的批次矫正](https://www.jianshu.com/p/4b10cff43920) 中的【第4部分 MNN矫正】

`fastMNN`先进行PCA降维，以加速下面的聚类环节

```r
set.seed(1000101)
merged.pbmc <- do.call(fastMNN, c(normed.sce, 
    list(subset.row=combined.hvg, BSPARAM=RandomParam())))

merged.pbmc
# class: SingleCellExperiment 
# dim: 5000 14948 
# metadata(2): merge.info pca.info
# assays(1): reconstructed
# rownames(5000): ENSG00000090382 ENSG00000163220 ...
# ENSG00000122068 ENSG00000011132
# rowData names(1): rotation
# colnames: NULL
# colData names(2): batch label
# reducedDimNames(1): corrected
# altExpNames(0):

# 这时就出现了批次信息
> table(merged.pbmc$batch)
pbmc3k pbmc4k pbmc8k 
  2609   4182   8157
```

检查一下结果，使用`lost.var` ，值越大表示丢失的真实生物异质性越多

* It contains a matrix of the **variance lost in each batch** (column) at each merge step (row).
* Large proportions of lost variance **(>10%)** suggest that correction is **removing genuine biological heterogeneity.**&#x20;

```r
metadata(merged.pbmc)$merge.info$lost.var
##         pbmc3k    pbmc4k   pbmc8k
## [1,] 7.003e-03 3.126e-03 0.000000
## [2,] 7.137e-05 5.125e-05 0.003003
```

### **对一个大数据进行聚类**

```r
g <- buildSNNGraph(merged.pbmc, use.dimred="corrected")
colLabels(merged.pbmc) <- factor(igraph::cluster_louvain(g)$membership)
table(colLabels(merged.pbmc), merged.pbmc$batch)
##     
##      pbmc3k pbmc4k pbmc8k
##   1     113    387    825
##   2     507    395    806
##   3     175    344    581
##   4     295    539   1018
##   5     346    638   1210
##   6      11      3      9
##   7      17     27    111
##   8      33    113    185
##   9     423    754   1546
##   10      4     36     67
##   11    197    124    221
##   12    150    180    293
##   13    327    588   1125
##   14     11     54    160
```

### **可视化**

```r
set.seed(10101010)
merged.pbmc <- runTSNE(merged.pbmc, dimred="corrected")
# 查看3个数据混合后的聚类结果，以及有没有批次效应
gridExtra::grid.arrange(
    plotTSNE(merged.pbmc, colour_by="label", text_by="label", text_colour="red"),
    plotTSNE(merged.pbmc, colour_by="batch")
)
```

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


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://jieandze1314.osca.top/04/04-4.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.