1 前言
基本上scRNA数据分析步骤都会包含基于基因表达量来对细胞进行比较。例如：细胞聚类就是将相似表达谱的细胞聚集在一起，而如何判断基因表达谱是否相似，就是计算基因间的欧氏距离。
scRNA分析就是在于数据打交道。本来看不见摸不着的基因，现在体现在数据上，虽然还是摸不着，但能看得到。每个基因现在都作为数据的一个维度存在。假如我们有一个scRNA数据集，其中只有两个基因，其中两个坐标轴表示两个基因的表达量，图中的点就表示细胞，每个细胞在图上都由x、y轴确定。也就是说，基因的表达决定了细胞的分布
如果有两个基因，那么细胞的位置就由两个维度决定；三个基因，就由 三维空间决定；而我们有2万多个基因，那么就存在2万多个维度（高维空间）。
高维空间我们无法理解，最常见的当属二维和三维，如果再加上时间线，就是四维空间。
降维，顾名思义，就是把各种离散的维度“去粗存精”。因为各个基因并不是毫无关联，同一生物学过程的各个基因是相互关联的，因此它们各自的维度就可以进行整合，结果可以用一个维度表示相互关联的多个维度，例如基因共表达网络中的eigengene
降维的作用非常重大，下游分析中的聚类只利用少量的维度进行计算，会比使用全部的维度计算更快，结果也更容易理解（比如我们目前最多只能做出三维的图片）。
数据准备之sce.zeisel
数据依然给大家准备好了，sce.zeisel链接：https://share.weiyun.com/mNwJS8U9 密码：g8379h
library(scRNAseq)
# 将每个细胞的所有基因表达量加起来，得到每个细胞的文库大小
library(scater)
sce.zeisel <- aggregateAcrossFeatures(sce.zeisel,
id=sub("_loc[0-9]+$", "", rownames(sce.zeisel)))
# 基因注释
library(org.Mm.eg.db)
rowData(sce.zeisel)$Ensembl <- mapIds(org.Mm.eg.db,
keys=rownames(sce.zeisel), keytype="SYMBOL", column="ENSEMBL")
# 质控（先perCellQCMetrics，后quickPerCellQC）
stats <- perCellQCMetrics(sce.zeisel, subsets=list(
Mt=rowData(sce.zeisel)$featureType=="mito"))
qc <- quickPerCellQC(stats, percent_subsets=c("altexps_ERCC_percent",
"subsets_Mt_percent"))
sce.zeisel <- sce.zeisel[,!qc$discard]
# 归一化（去卷积方法）
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.zeisel)
sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clusters)
sce.zeisel <- logNormCounts(sce.zeisel)
# 利用spike-in找HVGs（前10%）
dec.zeisel <- modelGeneVarWithSpikes(sce.zeisel, "ERCC")
top.hvgs <- getTopHVGs(dec.zeisel, prop=0.1)
sce.zeisel
## class: SingleCellExperiment
## dim: 19839 2816
## metadata(0):
## assays(2): counts logcounts
## rownames(19839): 0610005C13Rik 0610007N19Rik ... Zzef1 Zzz3
## rowData names(2): featureType Ensembl
## colnames(2816): 1772071015_C02 1772071017_G12 ... 1772063068_D01
## 1772066098_A12
## colData names(11): tissue group # ... level2class sizeFactor
## reducedDimNames(0):
## altExpNames(2): ERCC repeat
数据准备之PBMC
具体的操作在上一章也提到过
# 下载
library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
raw.path <- bfcrpath(bfc, file.path("http://cf.10xgenomics.com/samples",
"cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"))
untar(raw.path, exdir=file.path(tempdir(), "pbmc4k"))
# 读取
suppressMessages(library(DropletUtils))
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)
> dim(sce.pbmc)
[1] 33694 737280
# 基因注释之ID整合
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(
rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)
# 基因注释之添加位置信息
library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID,
column="SEQNAME", keytype="GENEID")
# 空液滴检测
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]
# 质控
stats <- perCellQCMetrics(sce.pbmc, subsets=list(Mito=which(location=="MT")))
high.mito <- isOutlier(stats$subsets_Mito_percent, type="higher")
sce.pbmc <- sce.pbmc[,!high.mito]
# 归一化（去卷积）
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
sce.pbmc <- logNormCounts(sce.pbmc)
# 利用数据分布来表示技术噪音，并挑选HVGs
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop=0.1)
2 主成分分析( PCA)概览
它的计算过程之前写过：StatQuest-PCA学习 以及 StatQuest--在R中拆解PCA，其中提到：
PCA的作用就是：对超过4维的数据降维到一个2D/3D平面图中，并且这个图中"相似相聚"
为什么可以通过PCA可以看批次效应？因为PCA图中的每个点都是一个样本，这个点中包含了大量的基因表达量信息；如果说本来属于生物学重复的几个样本在PCA图上离得很远，那么就意味着它们包含的基因表达量差异很大，这是不符合实际的，因此可能存在批次效应
PCA全称Principal components analysis，一般最常见的是二维的图，x轴是PC1，y轴是PC2。PC1捕获了细胞间最大的差异（暂且称它为V1），PC2与PC1正交，捕获了除V1外的最大差异。这样看来，前几个PCs就能捕获整个数据集主要的差异因素
我们降维也是想得到这几个能包含主要生物差异的PCs，方便下游分析【其实可以看到，一步步走来，都是上一步在为下一步省事省力做贡献，逐渐降低分析难度】。这几个包含主要生物差异的PCs其实也就是前几个PCs，PC越往后越集中在技术误差。
基本每个包都有自己PCA的函数，例如scater可以用runPCA ，seurat的RunPCA，monocle V3的reduce_dimension ，下面以scater包为例：
首先还是选择前多少的HVGs，然后给到runPCA默认计算前50个PCs，然后结果保存在SingleCellExperiment对象的reducedDims模块中
library(scran)
top.zeisel <- getTopHVGs(dec.zeisel, n=2000)
library(scater)
set.seed(100) #因为PCA是随机的
sce.zeisel <- runPCA(sce.zeisel, subset_row=top.zeisel)
reducedDimNames(sce.zeisel)
## [1] "PCA"
dim(reducedDim(sce.zeisel, "PCA")) # 50个PCs，2000多个细胞
## [1] 2816 50
如果是更大的数据，可以用近似的SVD算法（奇异值分解），scater或scran都可以直接通过函数计算SVD，利用参数BSPARAM=传递一个BiocSingularParam对象到runPCA中
SVD是一种矩阵分解方法，相当于因式分解，目的纯粹就是将一个矩阵拆分成多个矩阵相乘的形式
PCA从名字上就很直观，找到矩阵的主成分，也就意味这从一出生这就是个降维的方法。
对于稀疏矩阵来说，SVD更适用，这样对于大数据来说节省了很大空间
library(BiocSingular)
set.seed(1000)
sce.zeisel <- runPCA(sce.zeisel, subset_row=top.zeisel,
BSPARAM=RandomParam(), name="IRLBA")
reducedDimNames(sce.zeisel)
## [1] "PCA" "IRLBA"
> dim(reducedDim(sce.zeisel, "IRLBA"))
[1] 2816 50
3 选择合适的PCs数
看到这个问题，好像又回到了之前选择top HVGs的时候了，那么主成分数应该选多少合适呢？
如果选的PC多，可以避免丢弃一些生物信号，但同时又增加了噪音的风险
很多有经验的人会设置PC数在一个“合理”区间，例如10-50之间
下面我们来看看，如果经验不足，应该怎么帮助判断？
3.1 利用elbow point
elbow point就是在它之前变化幅度很大，之后变化幅度很小，属于一个转折点。
# 首先提取各个PCs对整体差异的贡献比例
percent.var <- attr(reducedDim(sce.zeisel), "percentVar")
可以看到前几个主成分的贡献最高，越往后越小
# 辅助选择
chosen.elbow <- PCAtools::findElbowPoint(percent.var)
chosen.elbow
## [1] 7
# 再画个图
plot(percent.var, xlab="PC", ylab="Variance explained (%)")
abline(v=chosen.elbow, col="red")
上面基于的假设是：每个PCs都能捕获一些生物差异，而且前面的PC比后面的PC包含的差异信息更多，更有价值。就像“二八准则”，仅有20%的变因操纵着80%的局面，PCA中也一样。于是当我们从头到尾观察每个PC的差异贡献时，会有一个明显的“落差”，“落差”之前的那些PC就是“二八准则”中20%的变因。
这个方法得到的PC数一般是比较少而精的，但我们能保证elbow近邻的那几个PCs没有感兴趣的生物差异信号吗？
3.2 利用技术噪音
这个方法是：设置一个差异贡献阈值，将每个PC的贡献率加起来，直到达到这个阈值，保留其中的所有PC。例如设置80%，也就是保留能够解释80%差异的PCs。
看到这里可能会奇怪，这个阈值的设置，不是也很主观吗？为了不显得这么主观，我们通过计算数据有多少比例的差异与生物因素相关，从而得到这个阈值。
更准确的描述是：Threshold is defined as the ratio of the sum of the biological components to the sum of total variances
整个计算是利用denoisePCA()函数完成，下面利用10X PBMC来展示
# 需要考虑技术噪音，通过参数technical指定
library(scran)
set.seed(111001001)
denoised.pbmc <- denoisePCA(sce.pbmc, technical=dec.pbmc, subset.row=top.pbmc)
ncol(reducedDim(denoised.pbmc))
## [1] 8
一般来说，denoisePCA()这种方法得到的PC数比elbow计算得到的要多，但是它在去除技术噪音的时候，并没有去除一些不感兴趣的生物噪音（例如transcriptional bursting）。它默认得到的PC数也是在5-50之间，取决于技术噪音与生物偏差的多少（例如：当技术误差占主导地位，就会得到很少的PC；当技术误差很小时，PC数又会非常多）
set.seed(001001001)
denoised.zeisel <- denoisePCA(sce.zeisel, technical=dec.zeisel,
subset.row=top.zeisel)
ncol(reducedDim(denoised.zeisel))
## [1] 50 #这个数据就达到了PC数的天花板
另外，modelGeneVarByPoisson() 和 modelGeneVarWithSpikes() 是考虑了技术误差，它们与denoisePCA的配合效果会更好。而基于log-counts估计的模型modelGeneVar() ，与denoisePCA的结果中，PC数就会减少
dec.pbmc2 <- modelGeneVar(sce.pbmc)
denoised.pbmc2 <- denoisePCA(sce.pbmc, technical=dec.pbmc2, subset.row=top.pbmc)
ncol(reducedDim(denoised.pbmc2))
## [1] 5
# 之前可是得到了8个PC
不论如何，想要更多的PCs（denoisePCA方法）还是保留更“精”的生物偏差（elbow方法），最终取决于自己的想法
3.3 利用细胞分群的推断
简言之就是根据亚群的数量来估计PC的数量，比如scRNA数据有10种不同的细胞类型，同时还要考虑一些噪音的存在，因此可以估计PC约等于10。不过这个想法有些理想，因为细胞分群的情况我们一开始不可能知晓，因此可以先进行预聚类
pcs <- reducedDim(sce.zeisel)
choices <- getClusteredPCs(pcs)
metadata(choices)$chosen
## [1] 17
plot(choices$n.pcs, choices$n.clusters,
xlab="Number of PCs", ylab="Number of clusters")
abline(a=1, b=1, col="red")
abline(v=metadata(choices)$chosen, col="grey80", lty=2)
图中红线是cluster数量的最大理论值，灰线就是getClusteredPCs()建议的PC数，可以看到如果使用推荐的17个PCs，基本能得到18个cluster
不过这个方法基于的假设是：整个数据可以分出许多亚群，而它们也是体现了不同的生物学意义。如果数据不能够分出这么多亚群，这种方法就不会得到太多的PCs，会导致信息丢失（例如研究连续生物学过程——分化时，可能就不会产生许多亚群）
3.4 自定义PC数量
如果之前获得了PC的数量，但还想再精简（比如之前得到了50个PC，现在可以根据自己需求来定义PC数量，比如取前20个）
可以直接在PCA结果上修改
reducedDim(sce.zeisel, "PCA") <- reducedDim(sce.zeisel, "PCA")[,1:20]
ncol(reducedDim(sce.zeisel, "PCA"))
## [1] 20
还可以新增一个PCA结果
# 赋值一个新的变量名
reducedDim(sce.zeisel, "PCA_20") <- reducedDim(sce.zeisel, "PCA")[,1:20]
reducedDimNames(sce.zeisel)
## [1] "PCA" "IRLBA" "PCA_20"
其实runPCA() 可以直接定义PC数，这个过程只是为了方便后期探索更多不同的组合，多设置几组PCs看看差异
4 了解非负矩阵分解
可以把它当成降维的高阶玩法
英文名是：Non-negative matrix factorization（NMF），非负矩阵分解 是一个用来替代主成分分析的方法。它利用两个低维的矩阵（W和H，一个代表基因，一个代表细胞）来估计整个表达矩阵，并且其中所有的数据都要求是非负。它的想法和PCA类似，都是“缩小”数据，NMF是利用小的矩阵来归纳大矩阵的主要特性，最终降噪、压缩数据。相比于PCA的坐标，NMF坐标更容易解释，因为更大的值表示相应因子中基因的表达量更高。
scRNA的数据使用NMF也是很常见的 (Shao and Höfer 2017; Kotliar et al. 2019) ，因此数据正好符合要求（虽然稀疏矩阵存在大量的0，但没有负数，并且即使采用log转换，还是会使用log(x+1) 的方法）
比如有一个LIGER算法，就首先使用综合非负矩阵分解（iNMF）来学习低维空间，获得一组特异因子，之后根据因子的maximum factor loading构建邻域图，从而达到降维目的
使用起来也不麻烦，scater包中有一个runNMF的函数（也是基于NNLM 包）。它的计算结果存储在reducedDims()的NMF模块中，之后就可以直接应用于分群。 相比于PCA的优点是：可以通过找到基因相关的小矩阵（W）每个因子中表达量高的基因，然后通过这些基因推测相应的细胞，再看细胞相关的小矩阵（H）中相应的高表达细胞，从而将细胞与细胞类型对应起来
# 运行很简单
set.seed(101001)
nmf.zeisel <- runNMF(sce.zeisel, ncomponents=10, subset_row=top.zeisel)
# 提取结果
nmf.out <- reducedDim(nmf.zeisel, "NMF")
nmf.basis <- attr(nmf.out, "basis")
colnames(nmf.out) <- colnames(nmf.basis) <- 1:10
# 画每个细胞与因子的关系
per.cell <- pheatmap::pheatmap(nmf.out, silent=TRUE,
main="By cell", show_rownames=FALSE,
color=rev(viridis::magma(100)), cluster_cols=FALSE)
# 画每个基因与因子的关系
per.gene <- pheatmap::pheatmap(nmf.basis, silent=TRUE,
main="By gene", cluster_cols=FALSE, show_rownames=FALSE,
color=rev(viridis::magma(100)))
# 组合
gridExtra::grid.arrange(per.cell[[4]], per.gene[[4]], ncol=2)
横坐标是定义好的10个因子，首先检查每个因子对应的高表达基因，看看其中有没有认识的marker基因。基于背景知识，Mog是少突细胞的marker基因，Gad1是中间神经元的marker基因。如果我们看到某个因子中Mog高表达（右图），那么这个因子就可能表示少突细胞。再看左图的细胞与因子的关系，找到对应的因子中高表达的细胞，那么这些细胞可能就是少突细胞。
因此，NMF在细胞类型识别中应该非常广。不过常规的流程还是基于PCA（入门选手必备），NMF可以用在更复杂的下游分析中（中高阶玩家选配）
5 降维后的可视化
一般的可视化算法都会基于10-50个PCs
5.1 PCA结果可视化
plotReducedDim(sce.zeisel, dimred="PCA", colour_by="level1class")
因为PCA是线性降维模型，每个PC捕获的是高维空间中线性的差异，不能够在前两个PC中压缩全部的维度。上图就能看到，前两个PC对亚群的显示效果不太理想，很多都没有分开。如果说第一个PC表示各个亚群之间最大的差异，PC2就是第二大差异，但是剩下的差异呢？没有显现！
虽然也可以在一张图画多个PC，但还是不好解释；另外即使用了这个方法，也不能将一些亚群分开
# 在一张图画多个PC
plotReducedDim(sce.zeisel, dimred="PCA", ncomponents=4,
colour_by="level1class")
PCA的优点就是：可预测、数据结构不引入人为因素、对输入数据的微小差异比较灵敏
5.2 t-SNE
t-SNE是由SNE(Stochastic Neighbor Embedding, SNE; Hinton and Roweis, 2002)发展而来。于2008年提出，全称是： t-stochastic neighbor embedding，可以说是现在scRNA分析的标准可视化方法。与PCA不同，它不局限于线性变换，也不需要精确地表示各个细胞群之间的距离。能在高维空间中直接捕获细胞非线性关系，它的分群更灵活，可以在更复杂细胞群体中区分亚群
（因此tSNE的点不表示数据的远近，相距较近的点也不意味着细胞相似）
set.seed(00101001101)
# runTSNE()的结果也是存储在reducedDims()中
sce.zeisel <- runTSNE(sce.zeisel, dimred="PCA")
plotReducedDim(sce.zeisel, dimred="TSNE", colour_by="level1class")
不过，tSNE的最大缺点是：计算量很大。不过可以在runtTSNE()中使用参数dimred="PCA" 来缓解一下；
缺点二：一般需要运行多次才能得到想要的结果，并且如果想要重复之前的结果，需要设置随机种子，因为tsne每次对细胞映射的坐标都不同（可能第一次运行这一团细胞在左下角，下一次又跑到右上角去了）
缺点三： 参数perplexity的使用。perplexity的含义是：the number of close neighbors each point has。 这个参数一般会多调整几次
关于tsne这个流行的算法，有必要了解一下：
tsne的作者Laurens强调，可以通过t-SNE的可视化图提出一些假设，但是不要用t-SNE来得出一些结论，想要验证你的想法，最好用一些其他的办法。
t-SNE中集群之间的距离并不表示相似度 ，同一个数据上运行t-SNE算法多次，很有可能得到多个不同“形态”的集群。但话说回来，真正有差异的群体之间，不管怎么变换形态，它们还是有差别
关于perplexity的使用：(默认值是30) 如果忽视了perplexity带来的影响，有的时候遇到t-SNE可视化效果不好时，对于问题无从下手。perplexity表示了近邻的数量，例如设perplexity为2，那么就很有可能得到很多两个一对的小集群，也就意味着：perplexity如果很小，那么分辨率更高。
有的时候会出现同一集群被分为两半的情况，但群间的距离并不能说明什么，解决这个问题，只需要跑多次找出效果最好的就可以了
引用自： https://bindog.github.io/blog/2018/07/31/t-sne-tips/
很好的tsne可视化：https://distill.pub/2016/misread-tsne/c
下面看看不同perplexity参数的影响：
set.seed(100)
sce.zeisel <- runTSNE(sce.zeisel, dimred="PCA", perplexity=5)
out5 <- plotReducedDim(sce.zeisel, dimred="TSNE",
colour_by="level1class") + ggtitle("perplexity = 5")
set.seed(100)
sce.zeisel <- runTSNE(sce.zeisel, dimred="PCA", perplexity=20)
out20 <- plotReducedDim(sce.zeisel, dimred="TSNE",
colour_by="level1class") + ggtitle("perplexity = 20")
set.seed(100)
sce.zeisel <- runTSNE(sce.zeisel, dimred="PCA", perplexity=80)
out80 <- plotReducedDim(sce.zeisel, dimred="TSNE",
colour_by="level1class") + ggtitle("perplexity = 80")
multiplot(out5, out20, out80, cols=3)
看到，t-SNE的结果看上去更像是一幅“地图”，但正是因为这样更容易引起误导，让我们根据clusters的大小和位置产生一些联想。事实上，上面的点不能说明什么问题，t -SNE的算法会让原本稠密的点变膨胀，让原本稀疏的点变紧凑，不可以用cluster的大小来衡量亚群的异质性。它并没有义务帮我们保留cluster的相对位置，因此我们不能根据点的位置远近来确定cluster之间的相似程度。
5.3 UMAP
2018年提出McInnes et al., (2018)，全称是：Uniform manifold approximation and projection，与t-SNE一样，也是一种非线性方法，都是在高维空间中寻找保持相邻关系的低维表示方法。不过它们的基础理论不同，因此图形展示也不同
set.seed(1100101001)
sce.zeisel <- runUMAP(sce.zeisel, dimred="PCA")
plotReducedDim(sce.zeisel, dimred="UMAP", colour_by="level1class")
可以看到，UMAP相比t-SNE，各个cluster更加紧凑，而且cluster之间的距离空间也更大。UMAP比t-SNE保留了更多的空间结构。运算速度也更快，因此对于大型数据是个福音（尽管如此，还是推荐使用前几个PC进行UMAP）。另外它也需要指定随机种子。
劣势一：
UMAP也有自己的复杂参数设定，例如近邻的数量n_neighbors 和 点之间最短距离min_dist 都会结果展示都比较大的影响。如果这些值太小，随机噪音就会显现出来；如果设置太大，就会丢弃一些原本不错的数据结构。所以建议也是：多试几组参数的设置
劣势二：
UMAP旨在保留更多的全局结构，但这必然会降低每个cluster中的分辨率。
5.4 最后，作者对于可视化的建议
Do not let the tail (of visualization) wag the dog (of quantitative analysis)
