单细胞交响乐
搜索文档…
⌃K
Links

3.4 降维

刘小泽写于2020.6.29

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如果很小,那么分辨率更高。
  • 有的时候会出现同一集群被分为两半的情况,但群间的距离并不能说明什么,解决这个问题,只需要跑多次找出效果最好的就可以了
下面看看不同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)