10X单细胞(10X空间转录组)降维分析之tSNE(算法基础知识)

hi,大家好,不知道大家对算法感不感兴趣,之前我分享了PCA的原理和一些分析技巧,大家可以参考单细胞PCA分析的降维原理,这篇主要讲PCA的降维原理,10X单细胞10X空间转录组降维分析之PCA轴的秘密,这篇主要讲PCA降维后的一些知识,10X单细胞(10X空间转录组)SeuratPCA分析之三---维度的选取,这篇主要讲我们在分析单细胞数据的时候,PCA维度的选取,而我们今天就是要在这个的基础上,进行下一步tSNE非线性降维,也就是数据的二维(或者三维)的降维可视化。

数据可视化是大数据领域非常倚重的一项技术,在数据可视化领域有很多值得尊敬的研究者,比如深度学习大牛Hinton(SNE)、Maaten(t-SNE)还有唐建大神,下面言归正传,我们从简单的基础知识开始。

说实在的,想要彻底搞清楚这些算法的原理并不轻松,需要长时间关注和积累。这里我把所需要知识和资料简要列出,供大家有针对性的了解。

降维

降维顾名思义就是把数据或特征的维数降低,一般分为线性降维和非线性降维,比较典型的如下:

  • 线性降维:PCA(Principal Components Analysis)、LDA(Linear Discriminant Analysis)、MDS(Classical Multidimensional Scaling),其中PCA我们分享的最多

  • 非线性降维:Isomap(Isometric Mapping)、LLE(Locally Linear Embedding)、LE(Laplacian Eigenmaps) 大家可能对线性降维中的一些方法比较熟悉了,但是对非线性降维并不了解,非线性降维中用到的方法大多属于流形学习方法。(但我们今天所分享的tSNE并不是流形学习的方法,UMAP是,PHATE也是

流形学习

流形学习(Manifold Learning)听名字就觉得非常深奥,涉及微分流行和黎曼几何等数学知识。当然,想要了解流形学习并不需要我们一行一行的去推导公式,通过简单的例子也能够有一个直观的认识。关于流行学习的科普文章首推pluskid写的《浅谈流行学习》,里面有很多通俗易懂的例子和解释。(这个对数学知识的要求更加苛刻,本人不是数学专业的人,但是我请教了一些数学专业的硕士,发现他们也讲不清楚)。

简单来说,地球表面就是一个典型的流形,在流形上计算距离与欧式空间有所区别。例如,计算南极与北极点之间的距离不是从地心穿一个洞计算直线距离,而是沿着地球表面寻找一条最短路径,这样的一条路径称为测地线。如下面所示的三幅图

流形与测地线

其中第一张图为原始数据点分布,红色虚线是欧式距离,蓝色实线是沿着流形的真实测地线距离。第二张图是在原始数据点的基础上基于欧式距离构造的kNN图(灰色线条,下面还会具体介绍kNN图),红色实线表示kNN图中两点之间的最短路径距离。第三张图是将流形展开后的效果,可以看到,kNN图中的最短路径距离(红色实线)要略长于真实测地线距离(蓝色实线)。

在实际应用中,真实测地距离较难获得,一般可以通过构造kNN图,在kNN图中寻找最短路径距离作为真实测地线距离的近似。这些知识大家应该都不陌生吧,单细胞数据分析常用的方法

t分布

大家在概率与统计课程中都接触过t分布的概念,从正态总体中抽取容量为N的随机样本,若该正态总体的均值为\mu,方差为{\sigma}^2。随机样本均值为\bar{x},方差为s^2=\frac{1}{N-1}\sum_{i=1}^{N}(x_i-\bar{x})^2,随机变量t可表示为:

此时我们称t服从自由度为n-1t分布,即t \sim t(n-1)

下图展示了不同自由度下的t分布形状与正态分布对比,其中自由度为1的t分布也称为柯西分布。自由度越大,t分布的形状越接近正态分布。

正态分布与t分布

从图中还可以看出,t分布比正态分布要“胖”一些,尤其在尾部两端较为平缓。t分布是一种典型的长尾分布。实际上,在稳定分布家族中,除了正态分布,其他均为长尾分布。长尾分布有什么好处呢?在处理小样本和一些异常点的时候作用就突显出来了。下文介绍t-sne算法时也会涉及到t分布的长尾特性。

kNN图

kNN图(k-Nearest Neighbour Graph)实际上是在经典的kNN(k-Nearest Neighbor)算法上增加了一步构图过程。假设空间中有n个节点,对节点v_i,通过某种距离度量方式(欧式距离、编辑距离)找出距离它最近的k个邻居{v_1,v_2,\cdots,v_k},然后分别将v_i与这k个邻居连接起来,形成k条有向边。对空间中所有顶点均按此方式进行,最后就得到了kNN图。
有关knn的知识可以参考我之前分享的文章机器学习的常用算法(生信必备)

当然,为方便起见,在许多场景中我们往往将kNN图中的有向边视为无向边处理。如下图是一个二维空间中以欧式距离为度量的kNN图。

kNN图样例

kNN图的一种用途上文已经提到过:在计算流形上的测地线距离时,可以构造基于欧式距离的kNN图得到一个近似。原因很简单,我们可以把一个流形在很小的局部邻域上近似看成欧式的,也就是局部线性的。这一点很好理解,比如我们所处的地球表面就是一个流形,在范围较小的日常生活中依然可以使用欧式几何。但是在航海、航空等范围较大的实际问题中,再使用欧式几何就不合适了,使用黎曼几何更加精确。

kNN图还可用于异常点检测。在大量高维数据点中,一般正常的数据点会聚集为一个个簇,而异常数据点与正常数据点簇的距离较远。通过构建kNN图,可以快速找出这样的异常点。

k-d树与随机投影树

刚才说到kNN图在寻找流形的过程中非常有用,那么如何来构建一个kNN图呢?常见的方法一般有三类:第一类是空间分割树(space-partitioning trees)算法,第二类是局部敏感哈希(locality sensitive hashing)算法,第三类是邻居搜索(neighbor exploring techniques)算法。其中k-d树和随机投影树均属于第一类算法。(这一部分我们可以简单看一下,如果对算法感兴趣可以深入了解一下)。

很多同学可能不太熟悉随机投影树(Random Projection Tree),但一般都听说过k-d树。k-d树是一种分割k维数据空间的数据结构,本质上是一棵二叉树。主要用于多维空间关键数据的搜索,如范围搜索、最近邻搜索等。那么如何使用k-d树搜索k近邻,进而构建kNN图呢?我们以二维空间为例进行说明,如下图所示:

k-d树示意图

上图是一个二维空间的k-d树,构建k-d树是一个递归的过程,根节点对应区域内所有点,将空间按某一维划分为左子树和右子树之后,重复根结点的分割过程即可得到下一级子节点,直到k-d树中所有叶子节点对应的点个数小于某个阈值。

有了k-d树之后,我们寻找k近邻就不用挨个计算某个点与其他所有点之间的距离了。例如寻找下图中红点的k近邻,只需要搜索当前子空间,同时不断回溯搜索父节点的其他子空间,即可找到k近邻点。

k-d树找k近邻

当然,搜索过程还有一些缩小搜索范围的方法,例如画圆判断是否与父节点的分割超平面相交等等,这里就不展开讨论了。

不过k-d树最大的问题在于其划分空间的方式比较死板,是严格按照坐标轴来的。对高维数据来说,就是将高维数据的每一维作为一个坐标轴。当数据维数较高时,k-d树的深度可想而知,维数灾难问题也不可避免。相比之下,随机投影树划分空间的方式就比较灵活,还是以二维空间为例,如下图所示:

随机投影树示意图

随机投影树的基本思路还是与k-d树类似的,不过划分空间的方式不是按坐标轴了,而是按随机产生的单位向量。有的同学说,这样就能保证随机投影树的深度不至于太深吗?随机产生的单位向量有那么靠谱吗?这里需要注意的是,我们所分析的数据处于一个流形上的,并非是杂乱无章的,因此从理论上讲,随机投影树的深度并不由数据的维数决定,而取决于数据所处的流形维数。(此处可参考Freund等人的论文《Learning the structure of manifolds using random projections》)

那么如何使用随机投影树寻找k近邻呢?当然可以采用和k-d树类似回溯搜索方法。但是当我们对k近邻的精确度要求不高时,可以采用一个更加简单巧妙的方式,充分利用随机投影树的特性。简单来说,我们可以并行的构建多个随机投影树,由于划分的单位向量都是随机产生的,因此每棵随机投影树对当前空间的划分都是不相同的,如下图所示

随机投影树找k近邻

例如我们想搜索红点的k近邻,只需要在不同的随机投影树中搜索其所处的子空间(或者仅回溯一层父结点),最后取并集即可。这样做虽然在构建随机投影树的过程中较为耗时耗空间,但是在搜索阶段无疑是非常高效的。

先了解SNE(stochastic neighbor embedding)

SNE即stochastic neighbor embedding,是Hinton老人家2002年提出来的一个算法,出发点很简单:在高维空间相似的数据点,映射到低维空间距离也是相似的。常规的做法是用欧式距离表示这种相似性,而SNE把这种距离关系转换为一种条件概率来表示相似性,在这里,邻居关系通过概率体现。假设在高维空间中有两个样本点xixj,xj以pi|j的概率作为xi的邻居,将样本之间的欧氏距离转化成概率值,借助于正态分布,此概率的计算公式为

图片.png

其中,σi表示以xi为中心的正态分布的标准差,上式除以分母的值是为了将所有值归一化成概率。由于不关系一个点与他自身的相似度,因此pj|i = 0 。投影到低维空间后仍然要保持这个概率关系。假设xixj投影之后对应的点为yiyj,在低维空间中对应的近邻概率记为qj|i,计算公式与上面的相同,但标准差统一设为1/√2 ,即:
图片.png

上面定义的是点xi与它一个邻居点的概率关系,如果考虑其他所有点,这些点构成一个离散型概率分布Pi,是所有其他样本点成为xi的邻居的概率。在低维空间中对应的概率分布为Qi,投影的目标是将两个概率分布尽可能接近,因此,需要衡量两个概率分布的之间的相似度或差异。在机器学习中一般用KL(Kullback-Leibler)散度衡量两个概率分布之间的距离(关于KL散度的大家可以参考文章KL散度,我们下一篇详细分享一下KL散度),在生成对抗网络,变分自动编码器中都有它的应用。假设x为离散型随机变量,p(x)和q(x)是它的两个概率分布,KL散度定义为:
图片.png

KL散度不具有对称性,即一般情况下KL(p||q)不等于 KL(q||p),因此不是距离。KL散度是非负的,如果两个概率分布完全相同,有极小值0。对于连续性随机变量,则将求和变成定积分。
由此得到投影的目标为最小化如下函数:
图片.png

这里对所有样本点的KL散度求和,l为样本数。把概率计算的公式代入KL散度,可以将目标函数写成所有yi的函数。该问题可以用梯度下降法求解,目标函数对yi的梯度为:
图片.png

计算出梯度之后可用梯度下降法迭代,得到的最优yi值即为xi投影后的结果。

接下来我们来学习t分布随机近邻嵌入

虽然SNE有较好的效果,但训练时难以优化,而且很容易导致拥挤的问题。

拥挤问题(The Crowding Problem)

所谓拥挤问题,顾名思义,看看SNE的可视化效果,不同类别的簇挤在一起,无法区分开来,这就是拥挤问题。有的同学说,是不是因为SNE更关注局部结构,而忽略了全局结构造成的?这的确有一定影响,但是别忘了使用对称SNE时同样存在拥挤问题。实际上,拥挤问题的出现与某个特定算法无关,而是由于高维空间距离分布和低维空间距离分布的差异造成的。

我们生活在一个低维的世界里,所以有些时候思维方式容易受到制约。比如在讨论流形学习问题的时候,总喜欢拿一个经典的“Swiss roll”作为例子,这只不过是把一个简单的二维流形嵌入到三维空间里而已。实际上真实世界的数据形态远比“Swiss roll”复杂,比如一个10维的流形嵌入到更高维度的空间中,现在我们的问题是把这个10维的流形找出来,并且映射到二维空间上可视化。在进行可视化时,问题就来了,在10维流形上可以存在11个点且两两之间距离相等。在二维空间中呢?我们最多只能使三个点两两之间距离相等,想将高维空间中的距离关系完整保留到低维空间是不可能的。

这里通过一个实验进一步说明,假设一个以数据点x_i为中心,半径为rm维球(二维空间就是圆,三维空间就是球),其体积是按r^m增长的,假设数据点是在m维球中均匀分布的,我们来看看其他数据点与x_i的距离随维度增大而产生的变化。

代码如下所示:

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import norm

npoints = 1000 # 抽取1000个m维球内均匀分布的点
plt.figure(figsize=(20, 4))
for i, m in enumerate((2, 3, 5, 8)):
    # 这里模拟m维球中的均匀分布用到了拒绝采样,即先生成m维立方中的均匀分布,再剔除m维球外部的点
    accepts = []
    while len(accepts) < 1000:
        points = np.random.rand(500, m)
        accepts.extend([d for d in norm(points, axis=1) if d <= 1.0]) # 拒绝采样
    accepts = accepts[:npoints]
    ax = plt.subplot(1, 4, i+1)
    ax.set_xlabel('distance') # x轴表示点到圆心的距离
    if i == 0:
        ax.set_ylabel('count') # y轴表示点的数量
    ax.hist(accepts, bins=np.linspace(0., 1., 50), color='green')
    ax.set_title('m={0}'.format(str(m)), loc='left')
plt.show()

结果如下图所示:

高维空间距离分布

从图中可以看到,随着维度的增大,大部分数据点都聚集在m维球的表面附近,与点x_i的距离分布极不均衡。如果直接将这种距离关系保留到低维,肯定会出现拥挤问题。如何解决呢?这个时候就需要请出t分布了。

t分布随机近邻嵌入是对SNE的改进。t-SNE采用对称的概率计算公式,另外在低维空间中计算样本点之间的概率时使用t分布代替了正态分布。

神奇的长尾——t分布

刚才预备知识部分说到,像t分布这样的长尾分布,在处理小样本和异常点时有着非常明显的优势,例如下面这个图:

异常点与t分布

从图中可以看到,在没有异常点时,t分布与高斯分布的拟合结果基本一致。而在第二张图中,出现了部分异常点,由于高斯分布的尾部较低,对异常点比较敏感,为了照顾这些异常点,高斯分布的拟合结果偏离了大多数样本所在位置,方差也较大。相比之下,t分布的尾部较高,对异常点不敏感,保证了其鲁棒性,因此其拟合结果更为合理,较好的捕获了数据的整体特征。 那么如何利用t分布的长尾性来改进SNE呢?我们来看下面这张图,注意这个图并不准确,主要是为了说明t分布是如何发挥作用的。

神奇的长尾

图中有高斯分布和t分布两条曲线,表示点之间的相似性与距离的关系,高斯分布对应高维空间,t分布对应低维空间。那么对于高维空间中相距较近的点,为了满足p_{ij}=q_{ij},低维空间中的距离需要稍小一点;而对于高维空间中相距较远的点,为了满足p_{ij}=q_{ij},低维空间中的距离需要更远。这恰好满足了我们的需求,即同一簇内的点(距离较近)聚合的更紧密,不同簇之间的点(距离较远)更加疏远。

在SNE中p<sub>i|j</sub>p<sub>j|i</sub>是不对等的,因此概率值不对称。可以用两个样本点的联合概率替代它们之间的条件概率解决此问题。在高维空间中两个样本点的联合概率定义为:

图片.png

显然这个定义是对称的,即=p<sub>ij</sub>=p<sub>ji</sub>。同样地,在低维空间两个点的联合概率为:

图片.png

目标函数采用KL散度,定义为:
图片.png

但这样定义联合概率会存在异常点问题。如果某一个样本x<sub>i</sub>是异常点,即离其他点很远,则所有的||x<sub>i</sub>-x<sub>j</sub>||2都很大,导致与x<sub>i</sub>有关的p<sub>ij</sub>很小,从而导致低维空间中的y<sub>i</sub>对目标函数的影响很小。解决方法是重新定义高维空间中的联合概率,具体为:

图片.png

其中,n为样本点总数。这样能确保对所有的x<sub>i</sub>有:

图片.png

因为:

图片.png

因此每个样本点都对目标函数有显著的贡献。目标函数对y<sub>i</sub>的梯度为:
图片.png

这中方法称为堆成SNE。

堆成SNE随便对SNE进行了改进,但还是存在拥挤的问题,各类样本降维后聚集在一起而缺乏区分度。解决方法是用t分布替代了高斯分布,计算低维空间的概率值。相比于正态分布,t分布更长尾,对异常点更健壮。如果在低维空间使用t分布,则在高维空间中距离近的点,在低维空间距离也要近;但在高维空间中距离远的点,在低维空间中距离要更远。因此,可以有效拉大各个类样本之间的距离。使用t分布后,低维空间的概率计算公式为:

图片.png

目标函数同样采用KL散度。采用梯度下降法求解,目标函数对y<sub>i</sub>的梯度为:
图片.png

同样地,求解梯度之后可以用梯度下降法进行迭代从而得到y<sub>i</sub>的最优解。

熟悉了这么多之后,很多Seurat分析单细胞数据的RunTSNE的参数我们就可以理解了:

library(Seurat)
help(RunTSNE)

参数的意义:

  • seed.use 这个参数最难就理解,我们知道计算机产生的随机数都是伪随机数,是利用算法产生的一系列数。因此,需要给函数一个随机值作为初始值,以此基准不断迭代得到一系列随机数。这个初始值就叫做随机种子。
  • tsne.method :Select the method to use to compute the tSNE. Available methods are
    tSNE和FIt-SNE(有关Flt-SNE也需要很多内容,我们有机会再分享一篇这个算法)。

分析起来是如此的简单,背后的原理确实十分的复杂,大家多多学习,不要仅仅局限于会跑别人的代码

图片.png

生活很好,有你更好

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 207,113评论 6 481
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 88,644评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 153,340评论 0 344
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 55,449评论 1 279
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 64,445评论 5 374
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,166评论 1 284
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,442评论 3 401
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,105评论 0 261
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,601评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,066评论 2 325
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,161评论 1 334
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,792评论 4 323
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,351评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,352评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,584评论 1 261
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,618评论 2 355
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,916评论 2 344

推荐阅读更多精彩内容