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分布

大家在概率与统计课程中都接触过分布的概念，从正态总体中抽取容量为的随机样本，若该正态总体的均值为，方差为。随机样本均值为，方差为，随机变量可表示为：

此时我们称服从自由度为的分布，即

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

正态分布与t分布

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

kNN图

kNN图(k-Nearest Neighbour Graph)实际上是在经典的kNN(k-Nearest Neighbor)算法上增加了一步构图过程。假设空间中有个节点，对节点，通过某种距离度量方式（欧式距离、编辑距离）找出距离它最近的个邻居，然后分别将与这个邻居连接起来，形成条有向边。对空间中所有顶点均按此方式进行，最后就得到了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把这种距离关系转换为一种条件概率来表示相似性，在这里，邻居关系通过概率体现。假设在高维空间中有两个样本点x_i和x_j,x_j以p_i|j的概率作为x_i的邻居，将样本之间的欧氏距离转化成概率值，借助于正态分布，此概率的计算公式为

图片.png
其中，σ_i表示以x_i为中心的正态分布的标准差，上式除以分母的值是为了将所有值归一化成概率。由于不关系一个点与他自身的相似度，因此p_j|i = 0 。投影到低维空间后仍然要保持这个概率关系。假设x_i和x_j投影之后对应的点为y_i和y_j,在低维空间中对应的近邻概率记为q_j|i,计算公式与上面的相同，但标准差统一设为1/√2 ，即：
图片.png
上面定义的是点x_i与它一个邻居点的概率关系，如果考虑其他所有点，这些点构成一个离散型概率分布P_i，是所有其他样本点成为x_i的邻居的概率。在低维空间中对应的概率分布为Q_i,投影的目标是将两个概率分布尽可能接近，因此，需要衡量两个概率分布的之间的相似度或差异。在机器学习中一般用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散度，可以将目标函数写成所有y_i的函数。该问题可以用梯度下降法求解，目标函数对y_i的梯度为：
图片.png

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

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

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

拥挤问题(The Crowding Problem)

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

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

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

代码如下所示：

# -*- 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()

结果如下图所示：

高维空间距离分布

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

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

神奇的长尾——t分布

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

异常点与t分布

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

神奇的长尾

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

在SNE中和是不对等的，因此概率值不对称。可以用两个样本点的联合概率替代它们之间的条件概率解决此问题。在高维空间中两个样本点的联合概率定义为：

图片.png

显然这个定义是对称的，即=。同样地，在低维空间两个点的联合概率为：

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

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

图片.png

其中，n为样本点总数。这样能确保对所有的有：

图片.png

因为：

图片.png
因此每个样本点都对目标函数有显著的贡献。目标函数对的梯度为：
图片.png

这中方法称为堆成SNE。

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

图片.png
目标函数同样采用KL散度。采用梯度下降法求解，目标函数对的梯度为：
图片.png

同样地，求解梯度之后可以用梯度下降法进行迭代从而得到的最优解。

熟悉了这么多之后，很多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

生活很好，有你更好