10X单细胞（10X空间转录组）降维分析之UMAP

追风少年ii

于 2024-04-26 10:05:07 发布

阅读量791

点赞数 13

文章标签：数据分析 python

本文链接：https://blog.csdn.net/weixin_53637133/article/details/138210252

版权

hello，大家好，关于降维算法，在之前分享的文章10X单细胞（10X空间转录组）降维分析之tSNE（算法基础知识）详细介绍了tSNE的降维原理，但是大家都知道，在我们分析10X单细胞或者10X空间转录组数据的时候，有两种降维方法，一种是tSNE，这个已经分享过，今天我们就来分享另外一个，UMAP。

什么是UMAP？

UMAP ，全称uniform manifold approximation and projection，统一流形逼近与投影，是基于黎曼几何和代数拓扑的理论框架结构构建的。在处理大数据集时，UMAP优势明显，运行速度更快，内存占用小。Etienne Becht等人2019年在Nature Biotechnology上发表一篇文章将其应用在生物学数据上并阐述了UMAP在处理单细胞数据方面的应用和优势。

UMAP应该说是目前最好的降维算法了，目前大部分10X单细胞的降维图谱都选择了UMAP，因为其能最大程度的保留原始数据的特征同时大幅度的降低特征维数。

这是《生命科学的数理统计和机器学习》的相关探讨，我们来深入探讨这一个降维算法，UMAP，它主导了当今的单细胞转录组学。

如果你不知道tSNE是什么，它是如何工作的，也没有读过2008年的革命性的van der Maaten & Hinton原稿，可以参考我的那文章10X单细胞（10X空间转录组）降维分析之tSNE（算法基础知识）。尽管tSNE对一般的单细胞基因组学和数据科学产生了巨大的影响，但人们普遍认为它有一些缺点，这些缺点很快将得到解决。（tSNE的缺点在上次分享的文章中也做过详细的介绍）。

究竟是什么让我们想放弃使用tSNE单细胞转录组学?这里我用简短的评论总结了几点:

在scRNAseq中，tSNE并不能很好地扩展快速增加的样本量。试图用FItSNE来加速它会导致大量的内存消耗，使得在计算机集群之外无法进行分析。（集群之痛），这个大家用过Seurat分析过的同学应该都知道。
tSNE没有保留全局数据结构，这意味着只有在集群距离内才有意义，而集群之间的相似性并不能保证，因此使用tSNE作为集群不是一个很好的主意。

tSNE实际上只能嵌入到2维或3维中，即仅用于可视化目的。

请注意，这对于更现代的FItSNE算法来说仍然是一个问题。
tSNE执行从高维到低维的非参数映射，这意味着它不利用驱动观察到的集群的特性(也称为PCA维度)。
tSNE不能直接处理高维数据，在将其插入tSNE之前，通常使用自动编码器或PCA执行预降维
tSNE的计算占用了大量的内存，在使用大型perplexity超参数时，这一点变得尤为明显，因为k近邻的初始步长(就像Barnes-Hut过程)变得低效，而且对时间减少也很重要。更现代的FItSNE算法也不能解决这个问题。

tSNE是一种比较简单的机器学习算法，可以用以下四个公式表示:

(1)定义了高维空间中任意两点之间观测距离的高斯概率，满足对称性规则。Eq。

(2)引入了困惑的概念作为一个约束,确定最优σ为每个样本。

(3)声明了低维嵌入中点对之间距离的学生t分布。学生t分布是为了克服嵌入低维时的拥挤问题。

(4)给出了Kullback-Leibler散度损失函数，将高维概率投影到低维概率上，并给出了梯度下降优化中使用的梯度的解析形式。

看看上面的图，我想说的是 t分布应该提供全局距离信息，因为它们将高维空间中相距较远的点推到低维空间中更远的点。

然而，这种良好的意愿被成本函数(KL-divergence)的选择所扼杀，我们将在后面看到其原因。

当了解UMAP时，我的第一印象是这是一种全新的、有趣的降维技术，它基于坚实的数学原理，因此与纯机器学习半经验算法tSNE非常不同。我的生物学同事告诉我，最初的UMAP论文“太数学化了”，看着论文的第二部分，我很高兴地看到，严格而准确的数学终于进入生命和数据科学。然而，在阅读UMAP文档和观看Leland McInnes在2018年SciPy大会上的演讲时，我感到困惑，觉得UMAP是另一种邻居图技术，它与tSNE如此相似，以至于我很难理解UMAP与tSNE究竟有什么不同。
从UMAP论文中，虽然Leland McInnes试图在附录c中总结UMAP和tSNE之间的差异，但是它们之间的差异并不是很明显。在这里，我将首先总结我所注意到的UMAP和tSNE之间的不同之处，然后尝试解释为什么这些不同是重要的，并找出它们的影响有多大。

UMAP在高维中使用指数概率分布，但不一定是像tSNE那样的欧氏距离，而是任何距离都可以代入。另外，概率没有归一化:

ρ是一个重要的参数,它代表了每个i数据点距离第一个最近邻。这保证了流形的局部连接性。换句话说，这为每个数据点提供了一个局部自适应指数核，因此距离度量在点与点之间变化。
ρ参数是唯一在UMAP部分2和3之间的桥梁。否则,我也看不出什么模糊单形建筑,即从第二节的拓扑数据分析,与从第三节UMAP的算法实现,似乎在一天结束的时候模糊单形集导致最近邻图施工。
UMAP对高维或低维概率都不应用标准化，这与tSNE非常不同，感觉很奇怪。然而,高或低维函数形式的概率可以看到他们已经扩展[0,1],事实证明,缺乏规范化的情况下,类似于情商分母。

(1),可以显著降低计算时间高维图像由于求和或集成是一个代价高昂的计算过程。想想马尔可夫链蒙特卡罗(MCMC)它基本上是试图近似地计算在贝叶斯规则的分母上的积分（UMAP使用最近邻居的数量而不是perplexity）

(2)定义perplexity, UMAP则定义了没有log2函数的最近邻居k的个数，即:

UMAP使用稍微不同的高维概率对称

symmterization是必要的因为UMAP融合在一起的点与本地不同的指标(通过参数ρ),它可能发生图A和B节点之间的重量不等于B之间的权重和节点。为什么UMAP使用这种对称而不是tSNE使用的对称还不清楚。我将在下一篇文章(从头开始编写UMAP)中展示我对不同的对称化规则的实验，这并没有使我相信这是如此重要的一步，因为它对最终的低维嵌入式产生了很小的影响。
UMAP使用曲线族1 / (1+a*y^(2b))在低维中建模距离概率，不是完全的学生t分布，但非常非常相似，请注意再次没有应用标准化:

其中，对于默认UMAP超参数a≈1.93,b≈0.79(实际上，对于min_dist = 0.001)。在实践中，UMAP从非线性最小二乘拟合到带有min_dist超参数的分段函数中找到a和b:

为了更好地理解曲线族1 / (1+a*y^(2b))的行为，让我们画出不同a和b的曲线:

plt.figure(figsize=(20, 15))
y = np.linspace(0, 10, 1000)
my_prob = lambda y, a, b: np.power(1 + a*y**(2*b), -1)
plt.plot(y, my_prob(y, a = 1, b = 1))``plt.plot(y, my_prob(y, a = 1.93, b = 0.79))
plt.plot(y, my_prob(y, a = 1.93, b = 5))
plt.gca().legend(('a = 1, b = 1', 'a = 1.93, b = 0.79', 'a = 1.93, b = 5'), fontsize = 20)
plt.title('Low-Dimensional Probability of Pairwise Euclidean Distances', fontsize = 20)
plt.xlabel('Y', fontsize = 20);
plt.ylabel('Q(Y)', fontsize = 20)
plt.show()`

我们可以看到曲线族对参数b非常敏感，在大的参数b处，在小的参数y处，曲线族形成了一种高峰。这意味着在UMAP超参数min_dist之下，所有的数据点都是同样紧密相连的。由于Q(Y)函数的行为几乎像一个Heaviside阶跃函数，这意味着UMAP为所有在低维空间中相互靠近的点分配了几乎相同的低维坐标。min_dist正是导致在UMAP维数降低图中经常观察到的超密集集群的原因。

为了演示如何准确地找到a和b参数，让我们展示一个简单的分段函数(其中高峰部分是通过min_dist参数定义的)，并使用函数族1 / (1+ay^(2b))通过优化来拟合它。curve_fit来自Scipy Python库。作为拟合的结果，我们得到了函数1 / (1+ay^(2b))的初值a和初值b。

from scipy import optimize
import matplotlib.pyplot as plt
import numpy as np
MIN_DIST = 1
x = np.linspace(0, 10, 300)
def f(x, min_dist):    
        y = []    
        for i in range(len(x)):        
                if(x[i] <= min_dist):            
                          y.append(1)        
                 else:            
                          y.append(np.exp(- x[i] + min_dist))    
        return y
dist_low_dim = lambda x, a, b: 1 / (1 + a*x**(2*b))    
p , _ = optimize.curve_fit(dist_low_dim, x, f(x, MIN_DIST))
print(p)
plt.figure(figsize=(20,15))
plt.plot(x, f(x, MIN_DIST), 'o')
plt.plot(x, dist_low_dim(x, p[0], p[1]), c = 'red')
plt.title('Non-Linear Least Square Fit of Piecewise Function', fontsize = 20)
plt.gca().legend(('Original', 'Fit'), fontsize = 20)
plt.xlabel('X', fontsize = 20)
plt.ylabel('Y', fontsize = 20)
plt.show()

UMAP使用二元交叉熵(CE)作为成本函数，而不是像tSNE那样使用kl -散度。

图片.png

由于我们需要知道交叉熵的梯度，以便以后实现梯度下降，让我们快速计算它。忽略只包含p(X)的常数项，我们可以将交叉熵重新写一下，并将其微分如下:

UMAP使用图拉普拉斯变换分配初始的低维坐标，与tSNE使用的随机正常初始化形成对比。

然而，这应该会对最终的低维表示产生较小的影响，至少在tSNE中是这样的。

但是，这应该会减少UMAP每次运行时的更改，因为它不再是一个随机初始化。

Linderman和Steinerberger提出了一个有趣的假设，即在tSNE的初始阶段最小化kl散度，并在早期进行放大，这一假设的动机是通过图拉普拉斯变换来进行初始化。

图拉普拉斯、谱聚类、拉普拉斯Eignemaps、扩散图、谱嵌入等，实际上是指将矩阵分解和邻接图方法结合起来解决降维问题的同一种有趣的方法。在这种方法中，我们首先构造一个图(或knn图)，然后通过构造拉普拉斯矩阵用矩阵代数(邻接矩阵和度矩阵)将其形式化，最后分解拉普拉斯矩阵，即求解特征值分解问题。

我们可以使用scikit-learn Python库，并使用spectralembedded函数在演示数据集(即与癌症相关的成纤维细胞(CAFs) scRNAseq数据)上轻松地显示初始的低维坐标:

from sklearn.manifold import SpectralEmbedding
model = SpectralEmbedding(n_components = 2, n_neighbors = 50)
se = model.fit_transform(np.log(X_train + 1))
plt.figure(figsize=(20,15))
plt.scatter(se[:, 0], se[:, 1], c = y_train.astype(int), cmap = 'tab10', s = 50)
plt.title('Laplacian Eigenmap', fontsize = 20)
plt.xlabel('LAP1', fontsize = 20)
plt.ylabel('LAP2', fontsize = 20)
plt.show()`

最后，UMAP使用随机梯度下降(SGD)代替常规梯度下降(GD)，如tSNE / FItSNE，这既加快了计算速度，又减少了内存消耗。

现在让我们简要地讨论一下为什么他们说tSNE只保留数据的局部结构。可以从不同的角度来理解tSNE的局部性。首先,我们有σ参数Eq。(1)本地数据点集这样互相“感觉”。因为成对欧几里得距离衰减指数的概率,在小的σ值,它基本上是零遥远的点(大型X)和快速增长仅为最近的邻居(小X)。相比之下,在大的σ,遥远而近点的概率成为限制可比和σ→∞,概率就等于1为所有任何一对点之间的距离,即成为等距点。

plt.figure(figsize=(20, 15))
x = np.linspace(0, 10, 1000)
sigma_list = [0.1, 1, 5, 10]
for sigma in sigma_list:    
        my_prob = lambda x: np.exp(-x**2 / (2*sigma**2))    

plt.plot(x, my_prob(x))
plt.gca().legend(('$\sigma$ = 0.1','$\sigma$ = 1', '$\sigma$ = 5', '$\sigma$ = 10'), fontsize = 20)
plt.title('High-Dimensional Probability of Pairwise Euclidean Distances', fontsize = 20)
plt.xlabel('X', fontsize = 20); 
plt.ylabel('P(X)', fontsize = 20)
plt.show()

有趣的是,如果我们扩大成对欧几里得距离的概率高维度成泰勒级数在σ→∞,我们会在第二近似幂律:

关于两两欧几里得距离的幂律类似于多维定标法(MDS)的成本函数，MDS是通过保存每对点之间的距离来保存全局距离，而不管它们是相距很远还是很近。一个可以解释这个大的σtSNE远程数据点之间的相互作用,所以是不完全正确的说tSNE只能处理当地的距离。然而，我们通常会受到perplexity有限值的限制，Laurens van der Maaten建议perplexity的取值范围在5到50之间，尽管在局部信息和全局信息之间可能会有一个很好的折衷，那就是使用平方根≈N^(1/2)来选择perplexity，其中N为样本量。相反的极限,σ→0,我们最终的极端“局部性”高维概率的行为类似于狄拉克δ函数的行为。

另一种理解tSNE“局部性”的方法是考虑KL-divergence函数。假设X是高维空间中点之间的距离Y是低维空间中点之间的距离

根据kl -散度的定义：

方程(9)的第一项对于X的大小都是趋近于0的，对于X的大小也是趋近于0的，因为指数趋近于1，而log(1)=0。对于大X，这一项仍然趋近于0因为指数前因子趋近于0的速度快于对数趋近于负无穷。因此，为了直观地理解kl散度，只考虑第二项就足够了:

这是一个看起来很奇怪的函数，让我们画出KL(X, Y)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import LinearLocator, FormatStrFormatter
fig = plt.figure(figsize=(20, 15))
ax = fig.gca(projection = '3d')# Set rotation angle to 30 degrees
ax.view_init(azim=30)
X = np.arange(0, 3, 0.1)
Y = np.arange(0, 3, 0.1)
X, Y = np.meshgrid(X, Y)
Z = np.exp(-X**2)*np.log(1 + Y**2)``# Plot the surface.
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,linewidth=0, antialiased=False)
ax.set_xlabel('X', fontsize = 20)
ax.set_ylabel('Y', fontsize = 20)
ax.set_zlabel('KL (X, Y)', fontsize = 20)
ax.set_zlim(0,2.2)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()

这个函数的形状非常不对称。如果点在高维度X之间的距离很小,指数因子变成1和对数项行为日志(1 + Y ^ 2)这意味着如果Y是在低维空间的距离大,将会有一个大的惩罚,因此tSNE试图减少Y在小X为了减少罚款。相反，对于高维空间中的长距离X, Y基本上可以是0到∞之间的任何值，因为指数项趋于0，并且总是胜过对数项。因此，在高维空间中相距遥远的点，在低维空间中可能会相互靠近。因此，换句话说，tSNE并不能保证高维空间中相距较远的点在低维空间中会保持较远的距离。然而，它确实保证了在高维空间中相邻的点在低维空间中保持相邻。所以tSNE不是很擅长远距离投射至低维,所以它只保留本地数据结构提供了σ不去∞。

与tSNE不同，UMAP使用交叉熵(CE)作为成本函数，而不是KL-divergence

这导致了地方-全球结构保护平衡的巨大变化。在X的小值处，我们得到了与tSNE相同的极限，因为第二项由于前因子和对数函数比多项式函数慢的事实而消失:

因此，为了使惩罚规则最小化，Y坐标必须非常小，即Y→0。这与tSNE的行为完全一样。但是，在大X的相反极限，即X→∞时，第一项消失，第二项的前因子为1，得到:

这里，如果Y很小，我们会得到一个很大的惩罚，因为Y在对数的分母上，因此，我们鼓励Y很大，这样，对数下的比率就变成了1，我们得到零惩罚。因此，我们在X→∞处得到Y→∞，所以从高维空间到低维空间的整体距离保持不变，这正是我们想要的。为了说明这一点，让我们绘制UMAP CE成本函数:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import LinearLocator, FormatStrFormatter
fig = plt.figure(figsize=(20, 15))
ax = fig.gca(projection = '3d') # Set rotation angle to 30 degrees
ax.view_init(azim=30)
X = np.arange(0, 3, 0.001)
Y = np.arange(0, 3, 0.001)
X, Y = np.meshgrid(X, Y)
Z = np.exp(-X**2)*np.log(1 + Y**2) + (1 - np.exp(-X**2))*np.log((1 + Y**2) / (Y**2+0.01))
# Plot the surface.
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,linewidth=0, antialiased=False)
ax.set_xlabel('X', fontsize = 20)
ax.set_ylabel('Y', fontsize = 20)
ax.set_zlabel('CE (X, Y)', fontsize = 20)
ax.set_zlim(0,4.3)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
# Add a color bar which maps values to colors.fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()

在这里，我们可以看到图的“右”部分看起来与上面的kl散度曲面非常相似。这意味着在X低的时候，为了减少损失，我们仍然想要Y低。然而，当X很大时，Y的距离也要很大，因为如果它很小，CE (X, Y)的损失将是巨大的。记住，之前，对于KL (X, Y)曲面，在X很大的情况下，我们在高Y值和低Y值之间没有差别，这就是为什么CE (X, Y)代价函数能够保持全局距离和局部距离。

我们知道UMAP是速度比tSNE担忧)时大量的数据点,b)嵌入维数大于2或3,c)大量环境维度的数据集。在这里,让我们试着了解UMAP要优于tSNE来自于数学和算法实现。

tSNE和UMAP基本上都包含两个步骤：

建立一个图在高维度的带宽和计算指数概率,σ,使用二进制搜索和固定数量的最近的邻居需要考虑。
通过梯度下降优化低维表示。第二步是算法的瓶颈，它是连续的，不能是多线程的。由于tSNE和UMAP都执行第二步，所以并不清楚为什么UMAP比tSNE更有效。

然而，我注意到UMAP的第一步比tSNE快得多。这有两个原因：

首先，我们去掉了最近邻数量定义中的log部分，即不像tSNE那样使用全熵:

由于在算法上，对数函数是通过泰勒级数展开来计算的，而且在线性项前面加上一个对数前因子实际上并没有增加多少，因为对数函数比线性函数慢，所以最好完全跳过这一步。
第二个原因是我们忽略了高维概率的归一化，即式(1)中用于tSNE的高维概率归一化。可以说，这一小步实际上对演出产生了戏剧性的影响。这是因为求和或积分是一个计算开销很大的步骤。

接下来，UMAP实际上在第二步中也变得更快了。这种改善也有几个原因:

采用随机梯度下降法(SGD)代替常规梯度下降法(GD)，如tSNE或FItSNE。这提高了速度，因为对于SGD，您可以从一个随机样本子集计算梯度，而不是像常规GD那样使用所有样本子集。除了速度之外，这还减少了内存消耗，因为您不再需要为内存中的所有样本保持梯度，而只需要为一个子集保持梯度。
我们不仅跳过了高维和低维概率的标准化。这也省去了第二阶段(优化低维嵌入)的昂贵求和。
由于标准的tSNE使用基于树的算法进行最近邻搜索，由于基于树的算法随着维数的增加呈指数增长，所以生成超过2-3个嵌入维数的速度太慢。这个问题在UMAP中通过去掉高维和低维概率的标准化得到了解决。
增加原始数据集的维数，我们引入稀疏性，即我们得到越来越多的碎片化流形，即有时存在稠密区域，有时存在孤立点(局部破碎流形)。UMAP解决这个问题通过引入本地连接ρ参数融合在一起(某种程度上)通过引入自适应稀疏区域指数内核,考虑了本地数据连接。这就是为什么UMAP(理论上)可以处理任意数量的维度，并且在将其插入主维度缩减过程之前不需要预降维步骤(自动编码器、PCA)的原因。

在这篇文章中，我们了解到尽管tSNE多年来一直服务于单细胞研究领域，但它有太多的缺点，如速度快、缺乏全球距离保存。UMAP总体上遵循了tSNE的哲学，但是引入了一些改进，例如另一个成本函数和缺少高维和低维概率的标准化。

UMAP在数据上的运用

UMAP在处理大数据时具有独特的优势

除了运行速度快，内存占用小等特点，UMAP在处理细胞学数据时还有一个大的优势，就是可以反映细胞群体之间分化的连续性和组织性。下面将通过文献中的数据【2】来为大家详细讲解。

1.UMAP体现细胞集群分化的连续性

对同一组数据分别进行tSNE和UMAP降维，该数据为多达30万个从8种不同组织富集得到的T细胞和NK细胞的样本，并使用Phenograph聚类把细胞分为6大类，每种颜色代表一种细胞。从图中可以看出，UMAP和tSNE都可以较好地把不同类别的细胞分开。但tSNE倾向于把相同细胞群划分为更多的群，如图显示，黑色圈中CD8 T细胞，在tSNE结果中，群数更多，距离更远。

同样这组数据用组织来源对UMAP和t-SNE图上细胞的进行颜色区分，可以观察到一个有意思的现象。与UMAP相比，t-SNE更加倾向于根据它们的来源来分离总体细胞。而 UMAP则会兼顾细胞群的类别和来源来排列，如图中在CD4 T细胞和CD8 T细胞群内，细胞的排列与来源也会有一定的规律性，都是大致从脐带血（CB）和外周血单核细胞（PBMC），到肝脏（Liver）和脾脏（Spleen），最后到一端的扁桃或另一端的皮肤（Skin）、肠道（Gut）和肺（Lung）。

通过驻留记忆T细胞标志物CD69/CD103、记忆T细胞标志物CD45 RO和naïve T细胞标志物CCR7表达群的分布，可以观察到UMAP可以展示出T细胞连续的分化阶段。而tSNE结果中，这些群之间也是连续的，但是却没有非常明显的沿轴结构。同样的现象也在造血细胞系统中被观察到。由此可见， UMAP在大数据集的处理时可以展现细胞集群的连续性。