rna_seq统计方法试读

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

2021-05-17  本文已影响0人  单细胞空间交响乐

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，因为其能最大程度的保留原始数据的特征同时大幅度的降低特征维数。

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

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

图片.png

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

图片.png

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

图片.png
(1)定义了高维空间中任意两点之间观测距离的高斯概率，满足对称性规则。Eq。
(2)引入了困惑的概念作为一个约束,确定最优σ为每个样本。
(3)声明了低维嵌入中点对之间距离的学生t分布。学生t分布是为了克服嵌入低维时的拥挤问题。
(4)给出了Kullback-Leibler散度损失函数，将高维概率投影到低维概率上，并给出了梯度下降优化中使用的梯度的解析形式。
图片.png

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

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

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

图片.png

  1. ρ是一个重要的参数,它代表了每个i数据点距离第一个最近邻。这保证了流形的局部连接性。换句话说，这为每个数据点提供了一个局部自适应指数核，因此距离度量在点与点之间变化。

  2. ρ参数是唯一在UMAP部分2和3之间的桥梁。否则,我也看不出什么模糊单形建筑,即从第二节的拓扑数据分析,与从第三节UMAP的算法实现,似乎在一天结束的时候模糊单形集导致最近邻图施工。

  3. UMAP对高维或低维概率都不应用标准化，这与tSNE非常不同，感觉很奇怪。然而,高或低维函数形式的概率可以看到他们已经扩展[0,1],事实证明,缺乏规范化的情况下,类似于情商分母。

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

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

图片.png

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

图片.png

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

图片.png

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

图片.png

为了更好地理解曲线族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()`
图片.png

我们可以看到曲线族对参数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()
图片.png
UMAP使用二元交叉熵(CE)作为成本函数，而不是像tSNE那样使用kl -散度。
图片.png

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

图片.png

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


图片.png

我们可以使用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()`
图片.png

最后，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()
图片.png

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

图片.png

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

图片.png

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

图片.png

根据kl -散度的定义：

图片.png

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

图片.png

这是一个看起来很奇怪的函数，让我们画出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()
图片.png

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

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

图片.png

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

图片.png

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

图片.png

这里，如果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()
图片.png

在这里，我们可以看到图的“右”部分看起来与上面的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基本上都包含两个步骤：

  1. 建立一个图在高维度的带宽和计算指数概率,σ,使用二进制搜索和固定数量的最近的邻居需要考虑。

  2. 通过梯度下降优化低维表示。第二步是算法的瓶颈，它是连续的，不能是多线程的。由于tSNE和UMAP都执行第二步，所以并不清楚为什么UMAP比tSNE更有效。

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

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

  1. 采用随机梯度下降法(SGD)代替常规梯度下降法(GD)，如tSNE或FItSNE。这提高了速度，因为对于SGD，您可以从一个随机样本子集计算梯度，而不是像常规GD那样使用所有样本子集。除了速度之外，这还减少了内存消耗，因为您不再需要为内存中的所有样本保持梯度，而只需要为一个子集保持梯度。

  2. 我们不仅跳过了高维和低维概率的标准化。这也省去了第二阶段(优化低维嵌入)的昂贵求和。

  3. 由于标准的tSNE使用基于树的算法进行最近邻搜索，由于基于树的算法随着维数的增加呈指数增长，所以生成超过2-3个嵌入维数的速度太慢。这个问题在UMAP中通过去掉高维和低维概率的标准化得到了解决。

  4. 增加原始数据集的维数，我们引入稀疏性，即我们得到越来越多的碎片化流形，即有时存在稠密区域，有时存在孤立点(局部破碎流形)。UMAP解决这个问题通过引入本地连接ρ参数融合在一起(某种程度上)通过引入自适应稀疏区域指数内核,考虑了本地数据连接。这就是为什么UMAP(理论上)可以处理任意数量的维度，并且在将其插入主维度缩减过程之前不需要预降维步骤(自动编码器、PCA)的原因。

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

UMAP在数据上的运用

image

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

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

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

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

image

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

image

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

2.UMAP具有更好的可重复性
image

对三组数据（Samusik、Wong、Han_400k）分别进行数据随机降低至100-200,000之间不同的数量级，形成小数据集。纵轴为小数据集与原始数据集的相关性，代表降维方法在不同数据量上的可重复性。UMAP表现最好，数据集越大，优势越明显。

简单比较一下tSNE和UMAP

下图是UMAP和t-SNE对一套784维Fashion MNIST高维数据集降维到3维的效果的比较。


图片.png

虽然这两种算法都表现出强大的局部聚类并将相似的类别分组在一起，但UMAP还将这些相似类别的分组彼此分开。另外，UMAP降维用了4分钟，而多核t-SNE用了27分钟。

2、UMAP参数

UMAP的两个最常用的参数：n_neighbors 和 min_dist，它们可有效地用于控制最终结果中局部结构和全局结构之间的平衡。


图片.png

最重要的参数是 n_neighbors ，近似最近邻居数。它有效地控制了UMAP局部结构与全局结构的平衡，数据较小时，UMAP会更加关注局部结构，数据较大时，UMAP会趋向于代表大图结构，丢掉一些细节。

第二个参数是 min_dist，点之间的最小距离。此参数控制UMAP聚集在一起的紧密程度，数据较小时，会更紧密。较大的值会更松散，而将重点放在保留广泛的拓扑结构上。

图片.png
上图可以通过https://pair-code.github.io/understanding-umap/自己调整一下参数看一下。

3、进一步比较UMAP与t-SNE

t-SNE和UMAP大部分的表现非常相似，但以下示例明显例外：宽而稀疏的cluster中有密集的cluster(如下图所示)。UMAP无法分离两个嵌套的群集，尤其是在维数较高时。

图片.png

UMAP在初始图形构造中局部距离的使用可以解释该算法无法处理情况的原因。由于高维点之间的距离趋于非常相似(维数的诅咒)，所以可能会因此将其混合在一起。

算法很难，所以懂的人才显得牛

天行健，君子以自强不息

上一篇下一篇

猜你喜欢

热点阅读