OSCA单细胞数据分析笔记8—Dimensionality reduction

生信技能树jimmy

发布于 2021-07-02 19:04:54

1.1K0

发布于 2021-07-02 19:04:54

对应原版教程第9章:http://bioconductor.org/books/release/OSCA/overview.html 在scRNA-seq中，根据成千上万个基因表达信息(维度)定义细胞间的距离是令人头痛的，最大程度不丢失有效信息的前提下，进行降维处理对于后续的cluster分群非常有必要；尤其对于我们只能观察到低维信息，提供可视化的方法。

笔记要点

1、关于降维的背景知识
2、PCA降维的简单理解与应用
3、选择最佳PCs数量的思路
4、降维可视化

1、关于降维的背景知识

（1）在单细胞表达矩阵中，细胞的维度定义就是：有多少个基因表达数据，就有多少维度，以表示每个细胞间的相对位置。根据维度信息，可计算细胞间的距离，用于分群；
（2）但是对于具有成千上万个维度信息的单细胞间距离计算是十分低效的；
（3）而且考虑到在细胞的生命活动中，多个基因的表达量是高度相关的，即可以用少数特征值来代表多数基因的表达量；
（4）基于上述因素，单细胞数据降维就是使用几十个维度的特征信息，来衡量细胞间的距离，大大减少计算量；并且可一定程度上去除技术误差，以及对细胞间相对位置的二维可视化提供便利。

2、PCA降维的简单理解与应用

（1）简单理解

PCA降维是针对多维复杂数据常用的线性降维手段
可以简单理解为是基于原始数据中心点的相同维度坐标系的重构，新的坐标系的坐标轴就称之为主成分(PC, principal componets)。一般PC1代表最能衡量细胞差异性（最大方差捕获）的新变量，PC2、PC3....次之。不过一般来说至少前10个新变量(主成分)就可以解释/代表原始细胞上千维度表达矩阵的90%的信息，实现的信息的高度压缩与有效利用。
参考教程视频（statquest）：【中字】主成分分析法（PCA）| 分步步骤解析看完你就懂了！(https://www.bilibili.com/video/BV1C7411A7bj?from=search&seid=1673872467165998571)

（2）PCA降维与scRNA

对scRNA进行PCA降维的前提假设是多数基因的表达是高度相关的，可以用少数特征维度“概括”多数基因相对冗余的高维数据。
scRNA降维，产生的排在前面的若干个主成分往往代表有生物意义的主成分指标。而排在后面的，仅捕捉到微小波动性的主成分往往代表着技术误差引起的转录水平扰动等。
综上，对于scRNA的降维结果，选取Top主成分，在尽可能不损失细胞异质性信息的前提下，大大降低的高维数据的复杂度，而且减少了技术误差的干扰。

（2）scate包的`runPCA()`函数

示例数据（获取方式见原教程）

sce.zeisel
# class: SingleCellExperiment 
# dim: 19839 2816 

#获得高变基因
library(scran)
dec.zeisel <- modelGeneVarWithSpikes(sce.zeisel, "ERCC")
top.zeisel <- getTopHVGs(dec.zeisel, n=2000)

根据高变基因，进行降维，以避免技术误差导致的“异质性”被降维时所捕获；

library(scater)
sce.zeisel <- runPCA(sce.zeisel, subset_row=top.zeisel)

降维结果默认计算前50个主成分，保存在sce对象的ReduceDims的slot里

reducedDimNames(sce.zeisel)
#[1] "PCA"
dim(reducedDim(sce.zeisel, "PCA"))
#[1] 2816   50
reducedDim(sce.zeisel, "PCA")[1:10,1:6]

10个细胞的前6个主成分的指标

观察每个主成分的细胞异质性(方差解释)的捕获比例

percent.var <- attr(reducedDim(sce.zeisel), "percentVar")
# [1] 24.5181077  7.1739169  4.8484962  2.7507716  2.3263866  1.4646539  1.0064506
# [8]  0.9452909  0.8281589  0.7215028  0.5226110  0.4959496  0.4708611  0.4485736
# [15]  0.4059012  0.3924625  0.3512433  0.3220025  0.3009666  0.2864609  0.2743744
# [22]  0.2523282  0.2360735  0.2284612  0.2105227  0.2088941  0.1931187  0.1874229
# [29]  0.1809175  0.1704797  0.1636240  0.1594421  0.1556688  0.1519302  0.1479599
# [36]  0.1436477  0.1391318  0.1368290  0.1341435  0.1314207  0.1300718  0.1237573
# [43]  0.1208721  0.1195724  0.1170765  0.1139158  0.1131542  0.1120158  0.1114791
# [50]  0.1105582

sum(percent.var)
#[1] 55.35963
sum(percent.var[1:10])
#[1] 46.58374

plot(percent.var, log="y", xlab="PC", ylab="Variance explained (%)")

结合上述统计与下图所示，基本前10个主成分的方差解释率远高于剩余的40个主成分的总和。

主成分的方差解释率

这就引出了下面一个问题：选择多少个PC用于接下来的PC合适？这类似上一节的问题(选择多少个hvg合适？)
一般情况下，选择的范围在10~50之间。事实上，由于后面的主成分的方差解释率很低，所以在不需要考虑计算量的情况下，PC选择的多少（10~50）不会特别影响后面的分群结果。
但还是有多种思路去提供一个最佳PC数量选择的参考。

3、选择最佳PCs数量的思路

3.1 scree plot for elbow point 碎石图找拐点

假设：PC1的方差解释率应该远大于PC2，PC2的方差解释率应远大于PC3...以此类推，直到发现连续两个PC的方差解释率十分接近的前一个PC，认为就是最后一个能够捕获生物异质性的主成分。表现在图中就是一个“拐点”
可通过PCAtools包的findElbowPoint()函数自动鉴别

chosen.elbow <- PCAtools::findElbowPoint(percent.var)
chosen.elbow
# 7
plot(percent.var, xlab="PC", ylab="Variance explained (%)")
abline(v=chosen.elbow, col="red")

3.2 基于技术误差的阈值

假设：Top PCs尽可能多的捕捉到是生物水平异质性，而后面PCs则多代表技术噪音。
基于此，我们将hvg的方差分解为生物水平误差a与技术噪音b两部分(a+b=1)。然后，当Top n个主成分的方差解释率之和达到a时，即为我们希望保留的PCs
scran包的denoisePCA()函数，需要提供计算hvg时的方差分解结果以及hvg

sce.pbmc  #获取见原教程
#class: SingleCellExperiment 
#dim: 33694 3985
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop=0.1)
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop=0.1)

library(scran)
set.seed(111001001)
denoised.pbmc <- denoisePCA(sce.pbmc, technical=dec.pbmc, subset.row=top.pbmc)
ncol(reducedDim(denoised.pbmc))
# 9

因为考虑到需要将方差分解为biological与technical var两部分的正确划分。作者认为modelGeneVar()方法可能会低估了生物水平差异性，而modelGeneVarByPossion()/modelGeneVarWithSpikes()方法相对来说更能反应真实的分布情况。因此denoisePCA()与后两种方差分解方法搭配使用，效果最佳~

3.3 PCs=subclusters

假设：如果一堆细胞可分为2群，则最多只需要两个维度即可区分开；如果可分为3群，则最多只需要三个维度即可区分开... 以此类推，当已知存在n个细胞群时，通过n个维度指标信息可以确保区分开。
但是细胞含有多少个潜在的cluster是未知的，而分群操作需要指定使用的PCs，这是矛盾的。但可以通过逐个尝试Top n个主成分，进行分群，得到m个cluster。然后找到n=m(或者相似)的点，即我们希望保留的PCs。
scran包的getClusteredPCs()函数

pcs <- reducedDim(sce.zeisel)
choices <- getClusteredPCs(pcs)
val <- metadata(choices)$chosen

plot(choices$n.pcs, choices$n.clusters,
     xlab="Number of PCs", ylab="Number of clusters")
abline(a=1, b=1, col="red")  # n=m
abline(v=val, col="grey80", lty=2)

最后作者还介绍了基于随机矩阵理论(RMT)的两种方法，的确是看不太明白，暂时不做记录了。有兴趣的朋友可以参考原教程。

4、降维可视化

二维平面是对于我们人类可接受的表征细胞间距离的可视化方式。因此，尽管上一步PCA已经降至50个维度以内，但在可视化呈现方面，仍需采取一定手段。

4.1 基于PCA

采用Top 2 即前两个组成分作为坐标轴进行可视化。缺点就是抛弃了剩余其它所有PCs的信息，势必会丢弃一定部分的有效信息。

plotReducedDim(sce.zeisel, dimred="PCA", colour_by="level1class")
#利用已有的level1class注释来辅助佐证我们的可视化效果

如下图所示部分clutser的分离效果并不明显。

4.2 基于t-SNE

t-stochastic neighbor embedding（非线性降维）
用二维表示高维空间中细胞间的距离，保持邻居细胞间的距离尽可能不变，但是距离较远的细胞变换后二维距离相对来说不能精确表示。
此外这种方法较耗计算量，故一般在使用PCA降维的基础上进行t-SNE处理；并且为保证图形可重现性，最好设置下种子

sed.seed(1)
sce.zeisel <- runTSNE(sce.zeisel, dimred="PCA")
plotReducedDim(sce.zeisel, dimred="TSNE", colour_by="level1class")

如下图，根据cluster的定义为一群距离相近的细胞，所以t-SNE的结果对于1个cluster内的细胞距离有参考意义；而对cluster间的较远距离的衡量是没有意义的。
关于t-SNE的图形调整细节，详见：[https://distill.pub/2016/misread-tsne/]

4.3 基于UMAP

uniform manifold approximation and projection（非线性降维）
相对t-SNE来说，UMAP将cluster压的更'紧凑'；更能维持cluster分布的总体结构；并且速度更快。
同样需要设置随机种子，以保证可重现性。

set.seed(1)
sce.zeisel <- runUMAP(sce.zeisel, dimred="PCA")
plotReducedDim(sce.zeisel, dimred="UMAP", colour_by="level1class")

综合考虑t-SNE与UMAP，虽然实现算法不同，但二者都是旨在维持高位空间中，邻近细胞的距离不变，在二维空间中可视化。但对于相距较远的细胞的二维呈现是没有参考意义的，切记决不能用于衡量cluster的相似性。可用于在之后的分群结果中，观察同一cluster内的细胞是否为距离相近的细胞。至于选择哪一种方法，各有优劣，可以都尝试一下，看看哪种展示方式更适合叙述接下来的生物学故事。