为什么使用PHATE进行降维

博文名称：Why you should be using PHATE for dimensionality reduction
博文链接：https://towardsdatascience.com/why-you-should-be-using-phate-for-dimensionality-reduction-f202ef385eb7
发表时间：Jul 20, 2021

---

作为数据科学家，我们经常要处理超过3个感兴趣的特征或维度的高维数据。在监督机器学习中，我们可能会使用这些数据进行训练和分类，并且可能会减少维度以加快训练速度。在无监督学习中，我们使用这种类型的数据进行数据可视化和聚类。例如，在单细胞RNA测序(scRNA-seq)中，我们会检测超过一百万个细胞数据，每个细胞拥有数万个基因的测量值。大量数据给我们提供了了解细胞身份、状态和其他属性的窗口。更重要的是，这些属性使它和数据集中的无数其他细胞相关。尽管如此，这会创建一个包含1,000,000个细胞和10,000个基因的庞大矩阵，每个基因代表一个维度。我们如何解释这些数据？作为3-D世界中的人类，我们无法看到超过三个物理维度之外的任何东西，需要一种方法来捕捉这些数据集的本质，同时又不失去任何价值。

单细胞数据分析当中，最关键的步骤可以说是细胞分群聚类。我们如何将这个数据集压缩到2维或3维并且保留基本信息？你的第一直觉可能是使用PCA或tSNE将数据集投影到二维嵌入空间中。但是，这些方法中的每一种都有关键的取舍（权衡），可能会导致在下游分析期间对相关数据集得出错误的结论。使用PCA，会以捕获全局特征为代价，丢失数据点（或我们案例中的细胞）之间的局部关系。

此外，PCA是一种线性方法，它不能准确地捕获复杂的数据集，例如单细胞RNA测序数据，其中我们可以拥有由不同基因表达组合定义的无数细胞类型，同时还会经历各种过程，例如细胞分裂、分化（将干细胞特化为更成熟的细胞类型，如神经元）、新陈代谢等。tSNE是一种非线性方法，它虽然在保持局部关系方面做得更好，但代价是破坏了全局图景。例如，在干细胞分化中，干细胞不会轻按一下开关（或本例中的一组基因）自动变成神经元，而是经历在转录谱上逐渐改变引起的连续变化。我们经常将这种连续体描述为轨迹（trajectory）。tSNE的问题在于它打破了这些轨迹，导致细胞群脱节，几乎没有关于一种细胞类型如何与另一种相关的信息。

既然上面已经说明了这种有利弊的下游问题，我们该如何解决呢？PHATE（Potential of Heat diffusion for Affinity-based Transition Embedding）一种新的降维方式，就登场了，它是降维世界的新手。与tSNE一样，它是一种非线性、无监督的技术。然而，与tSNE不同，tSNE通常以牺牲全局结构为代价保留高维数据的局部结构，PHATE兼具了PCA和tSNE的优点，保留了数据点之间的局部和全局关系，准确地反映所讨论的高维数据集。

我们将继续以scRNA-seq数据为例。有了这个数据集，我们基本上有一个由m个细胞（行）和n个基因（列）组成的m×n矩阵，它们代表了信使RNA分子的离散计数。

我们首先计算这些细胞之间的欧几里得距离（Euclidean distance）的方阵，它只是基于笛卡尔坐标的两点之间线段的长度。在scRNA-seq的背景下，这些笛卡尔坐标是基因表达测量值，因此直观地，我们预计相距很近的细胞在基因表达方面非常相似，因此细胞类型相似。相距较远的细胞具有非常不同的基因表达模式，因此反映了截然不同的细胞类型（例如，神经元和红细胞）。然而，这个距离指标并不总是适合这种解释。维度诅咒（Curse of Dimensionality），其中太多维度可能导致数据集中的数据点与所有其他数据点的距离相等，从而难以得出数据趋势的结论，推导出有意义的聚类和局部邻域以及确定其他类型的模式。

为了解决这个问题，我们将我们的距离转换为亲和力指标，这量化了我们数据中观察值之间的局部相似性（PHATE中“A”中“A”的来源，表示“Affinity-based”）。这些亲和力与距离成反比，因此两个观测值在欧几里德空间中距离越远，它们的亲和力就越小；同样，距离越近，亲和力越大。亲和力通常是通过使用核函数来转换欧几里得距离来计算的。用非常简单的话说，这是你的概率质量或密度函数减去对其进行归一化以确保概率介于0和1之间的因子/系数。它们经常用于其他机器学习问题，例如支持向量机。里面的内核是主流的高斯内核（Gaussian kernel）：

image.png

其中x和y是高维空间X中的坐标,ε是带宽度量（bandwidth），用于测量该内核捕获的邻域的“扩展”或半径。PHATE论文的作者使用了一个稍微高级一点的核函数，它在量化相似性方面做得更好，为简洁起见，我将省略其详细信息。虽然这是保留数据集局部结构的一个方便的技巧，但与tSNE的情况一样，仅嵌入这些将破坏全局结构。因此，除了保留局部结构之外，PHATE的另一个目标是维护跨数据的全局关系。为了实现这一点，PHATE使用亲和力通过马尔可夫随机游走（Markov random-walk）在数据中“扩散”。当我们说扩散时，我们指的是从高浓度区域到低浓度区域的净扩散。在亲和力的背景下，这可以被认为是从高亲和力（即我们高维数据集中的一组细胞）到低亲和力（即更分散的细胞）。更直观地，我们可以认为这是房间内的热量从温暖的来源（例如壁炉）传播到较不温暖的区域（例如，您在沙发上），这可以在数学上建模为热量方程，其解是热核。这就是PHATE（Heat diffusion，热扩散）中的“H”的来源。

至于随机游走，这表示通过我们的高维空间（即从一个细胞到另一个细胞的转换）的连续随机步骤的轨迹，其中每个可能的步骤或转换都有一个确定的沿着所述路线的概率。你可以将其视为房间内的热量最初向房间一个角落的方向随机传播，然后切换方向。这被更广泛地称为随机过程。在我们的单细胞数据中，我们有细胞i去细胞j的概率取决于我们访问的最后一个细胞。现在抛开所有这些术语，让我们看看它是如何运作的！

首先我们通过对先前计算过的亲和力进行归一化来计算随机游走的初始概率。这会产生以下结果：

image.png

在这里，

image.png

这为我们提供了在单个时间步长内从细胞x移动到细胞y的N×N转移概率矩阵。

PHATE论文中使用的这个矩阵的一个更简单的术语是扩散算子。从数学上讲，为了获得有效的扩散，我们将扩散算子提高到最佳步数t以学习我们数据的全局结构。这给了我们在t时间步长内从细胞x转换到细胞y的概率。使用更大的t，我们可以在高维空间中覆盖更多距离，并了解更多关于全局结构的信息，而不会被我们基于单时间步移亲和性概率的局部性所困扰。你可以将其视为寻找远足小径以构建周围环境的地图：每走几步，您就放置一个标记以记录您之前的位置，通常位于明显的区域（例如，一棵大树、河岸等），而不是在每一步都设置标记。这使你可以构建该区域的一般地图，而不会被每根树枝和树枝所困扰。为简洁起见，我将省略如何计算最佳步骤t的技术细节，这在PHATE论文中使用冯诺依曼熵进行了描述，但这都将输入我们的“转换”算法的“T”。

好的，我们有了我们的亲和力，一方面可以捕捉附近细胞之间的局部关系，另一方面可以用强力扩散算子来学习我们的全局空间。怎么办？嵌入这个强大的扩散算子？还没那么快。一个限制是，通过这个算子（或作者定义的扩散距离，diffusion distances）得到的基于概率的细胞之间的距离对远点之间的距离不是很敏感，当我们考虑我们的高维空间边界点时，可能会遇到稳定性问题（本文更详细地介绍了这些缺点）。然而，这可以通过我们首字母缩略词“P”的第一个字母来解决，“P”代表潜力（Potential）。作者受信息论和随机动力学的启发，定义了这个称为潜在距离的巧妙度量，我们从幂扩散算子测量对数转换概率之间的距离。这增加了我们得到的距离的敏感性，并使PHATE实现数据可视化保留局部和全局结构特征。这在数学上定义如下：

image.png

在这里，

image.png

pₓᵗ指的是我们扩散算子的第x行，它是我们的转移概率矩阵的次方。

当我们说这个距离度量更敏感时，假设从细胞a到细胞b的转换概率为0.4，而对于细胞a到c的转换概率为0.5。在这些扩散距离下，它们对fold-changes不是很敏感。距离0.1在尝试将这种关系封装在较低维度的投影中时可能会遇到敏感性问题。但是，如果我们对这些概率进行对数变换并取它们的距离，我们会得到更大的距离0.223（log(0.5)-log(0.4)），这与这些概率分别为0.4和0.5相同（回想一下loga-logb=log(a/b)）。

好的，我们已经准备好PHATE的“E”部分，嵌入（embedding）。通常，使用这些扩散度量，我们通常会执行特征分解（即，通过相关矩阵的特征向量将其分解，在本例中为我们的幂扩散算子）以导出数据的扩散图，这是一种流行的方法用于研究scRNA-seq数据中的分化轨迹。然而，这种方法的问题在于它将轨迹打散成无数反映扩散分量的特征向量。因此，这种高内在维度使它们无法进行可视化。为了绕过这个限制，作者使用度量多维缩放（MDS度量）。这是一种嵌入方法，通过最小化所谓的“压力”函数，为距离矩阵作为输入量身定制：

image.png

虽然这个方程可能看起来很吓人，但它的本质是测量拟合优度，在这种情况下，嵌入坐标与我们寻求可视化的高维数据的拟合程度。应力越小，贴合度越高。因此，如果这些嵌入点的应力为零，则在此MDS嵌入中已成功捕获数据。由于噪声或少量MDS维数（即，m=2或3），可能会出现小的非零值。然而，只要保留轨迹和其他关键特征，这种程度的失真通常是可以容忍的。

到目前为止，我所讨论的所有内容都可以从PHATE论文作者的这张图中直观地总结出来：

Moonetal.,NatureBiotechnology2019，可视化高维生物数据中的结构和转换

我们现在准备学习如何在Python中实现PHATE并查看它的实际效果！为简单起见，我们将使用通用的MNIST数据集的较小版本，其中有8×8的手绘数字图像（你可以在此处找到论文作者将PHATE应用于真实scRNA-seq数据的示例））

github网址：https://github.com/KrishnaswamyLab/PHATE

使用pip或pip3安装 PHATE 库：

pip install phate

在你最喜欢的编辑器或Jupyter notebooks 中的notebook中创建一个新的Python文件并运行以下命令，保存为mnist_phate.py：

# Load necessary packages
import phate
import numpy as np
import pandas as pd
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
np.random.seed(26)

# Load our MNIST dataset
from sklearn import manifold, datasets

X, y = datasets.load_digits(return_X_y=True)
n_samples, n_features = X.shape

# This is a 1797 by 64 dataset
print(X.shape)
digits = load_digits()
# Convert to a pandas dataframe
digits_df = pd.DataFrame(X)

# Run PHATE: this will take a little under 10 seconds to run depending on your machine
phate_operator = phate.PHATE(n_components=2)
digits_phate = phate_operator.fit_transform(digits_df)

fig, ax = plt.subplots()
scatter = ax.scatter(digits_phate[:,0], digits_phate[:,1], c=digits.target, cmap = 'tab10')
ax.set_xlabel('PHATE_1')
ax.set_ylabel('PHATE_2')
ax.set_title("PHATE of MNIST")
legend1 = ax.legend(*scatter.legend_elements(),bbox_to_anchor=(1.01, 1),
                    loc="upper left", ncol = 1, title="Digits")
ax.add_artist(legend1)
fig.savefig("PHATE.png", bbox_inches="tight")

# Now run tSNE for comparison
digits_tsne = TSNE(n_components=2).fit_transform(digits_df)
fig, ax = plt.subplots()
scatter = ax.scatter(digits_tsne[:,0], digits_tsne[:,1], c=digits.target, cmap = 'tab10')
ax.set_xlabel('tSNE_1')
ax.set_ylabel('tSNE_2')
ax.set_title("tSNE of MNIST")
legend1 = ax.legend(*scatter.legend_elements(),bbox_to_anchor=(1.01, 1),
                    loc="upper left", ncol = 1, title="Digits")
ax.add_artist(legend1)
fig.savefig("tSNE.png", bbox_inches="tight")

# Now let's compare to PCA
pca = PCA(n_components=2)
digits_pca = pca.fit_transform(digits_df)
fig, ax = plt.subplots()
scatter = ax.scatter(digits_pca[:,0], digits_pca[:,1], c=digits.target, cmap = 'tab10')
ax.set_xlabel('PCA_1')
ax.set_ylabel('PCA_2')
ax.set_title("PCA of MNIST")
legend1 = ax.legend(*scatter.legend_elements(),bbox_to_anchor=(1.01, 1),
                    loc="upper left", ncol = 1, title="Digits")
ax.add_artist(legend1)
fig.savefig("PCA.png", bbox_inches="tight")

# Finally, let's show all three of these methods side-by-side
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(25,9))
fig.suptitle('PCA vs tSNE vs PHATE')
ax1.scatter(digits_pca[:,0], digits_pca[:,1], c=digits.target, cmap = 'tab10')
ax1.set_xlabel('PCA_1')
ax1.set_ylabel('PCA_2')
ax1.set_title("PCA of MNIST")
ax2.scatter(digits_tsne[:,0], digits_tsne[:,1], c=digits.target, cmap = 'tab10')
ax2.set_xlabel('tSNE_1')
ax2.set_ylabel('tSNE_2')
ax2.set_title("tSNE of MNIST")
ax3.scatter(digits_phate[:,0], digits_phate[:,1], c=digits.target, cmap = 'tab10')
ax3.set_xlabel('PHATE_1')
ax3.set_ylabel('PHATE_2')
ax3.set_title("PHATE of MNIST")
legend1 = ax3.legend(*scatter.legend_elements(),bbox_to_anchor=(1.01, 1),
                    loc="upper left", ncol = 1, title="Digits")
fig.savefig("MNIST_all_dim_red.png", bbox_inches="tight")

让我们运行这个程序看看结果

python3 mnist_phate.py

image.png

查看 PHATE 嵌入，你可以看到我们总体上对数字进行了清晰的分离，而cluster本身具有独特的形状和分布。让我们看看这与 PCA 和 tSNE 相比如何。

image.png

PCA（左）数据过于拥挤，难以得出具体结论。 另一方面，tSNE（中间）将它们分离成球状cluster这方面做得很好，保留了数字与另一个数字之间的清晰的全局结构关系，以及抑制cluster内数据的独特传播。 PHATE（右）协调了这些问题，你可以清楚地分辨出cluster并了解它们之间的关系（例如，3 和 9 有很多相似之处）。 对于数字来说，这不是什么大问题，但对于单细胞生物学而言，当我们研究连续过程时，例如将干细胞转化为神经元，增加有关细胞如何从一种细胞类型转变为另一种细胞类型的细节非常有用。

在本文中，我们了解了PHATE算法、它的工作原理以及如何在Python中实现它。 如你所见，它在保留数据的局部和全局方面做得很好。 有关PHATE算法背后的数学的更多信息，以及更多示例和与 tSNE 和 UMAP 等其他方法的比较，我鼓励你查看原始出版物以及作者编写的一些关于 PHATE 和单细胞数据分析的教程。

参考

[1] K. Moon, D. van Dijk, Z. Wang, S. Gigante, D. Burkhart, W. Chen, K. Yim, A. van den Elzen, M.J. Hirn, R.R. Coifman, N.B. Ivanova, G. Wolf, and S. Krishnaswamy, Visualizing structure and transitions in high-dimensional biological data (2019), Nature Biotechnology

---

案例实操
注意要加上这句代码，否则执行报错，我在前面的代码中补充了。
digits = load_digits()

image.png

---

问题：如何理解tSNE的局部结构和全局结构？

刚看作者的这篇技术博文，有个点很疑惑，tsne的特点中解释局部结构和全局结果，没理解清楚：
前面文章说明：

tSNE, on the other hand, is a non-linear method that does a better job preserving local relationships but at the cost of shattering the global picture.
tSNE是一种非线性方法，它在保持局部关系方面做得更好，但代价是破坏了全局图景。

后面案例中说明：

image.png
总共有9个cluster，一个cluster代表一个数字，同一个颜色表示。我们可以看到9个颜色。
tSNE保留了数字与另一个数字之间的清晰的全局结构关系，以及抑制cluster内数据的独特传播。优点保留全局结构，缺点丢失局部结构。
我的感觉是这段话，前后矛盾，前面指出tSNE保留局部结构，缺点丢失全局结构。后面却说清晰的全局结构关系。

那么，什么是全局结构和局部结构？
找到一段英文

global structure of the data (connections between all the data points)
local structure of the data (connections between neighboring points)
local structure that is to group neighboring data points together
within clusters (local structure) and between clusters (global structure)

全局结构是指所有数据点的连接，特指cluster之间的连接和距离关系；
局部结构是指邻近数据点的连接，指某个cluster内部数据点之间的连接程度和距离远近；

中间的tSNE图，每个cluster内部都聚集得很好，成一个小球状，说明局部结构很好，但是9个cluster之间却间隔甚远，看不出cluster彼此的关系，即全局结构不够好。是我误读了作者的本意。

PCA 的目的是通过降维，捕获数据的全局结构。但是， PCA只能捕获特征中的线性结构。而t-SNE 算法，它侧重于在某些映射到低维数据时保留高维数据的局部距离。而PHATE兼顾了全局结构和局部结构，作者的意思是更佳。

2.PCA，UMAP，t-SNE，PHATE对比

• PCA（ Principal Component Analysis ）
  是传统的一种线性降维方法，做法是将数据点投射在一个新的维度，在此空间找出能够解释最大变异程度新数据点的代表性组合（即为主成分），藉此数学计算，故能将差异较大的数据点分开，找出群聚间的关系；然而 PCA 却不太适用于复杂的数据（例如单细胞数据，一个基因代表一个维度空间），线性计算方法不足以在人肉眼可见的 2D/3D 图片当中，将如此复杂且多维度的数据点做区隔，因此细胞点若以 PCA 方式呈现，会发现细胞会有重叠的现象。
  为了解决这个问题，后来才衍生出 t-SNE、UMAP 方法，对于细胞聚类展现较佳的视觉呈现。

• t-SNE（t-distributed stochastic neighbour embedding）
  是一种非线性降维方法，此方法是改良自 SNE （stochastic neighbour embedding）的计算，相对于SNE主要是改变了 (1)代价函数 cost function 的计算（考虑对称性） (2)分布规则（改为t-分布假设）。计算概念是取得每个数据点跟其他数据点的距离，平均两点之间的相互距离之后做 scaling，并假设数据点为 t-分布，最后尽可能将相近分数的数据点排列在一起。

• UMAP （Uniform Manifold Approximation and Projection）
  计算大量数据的时候 t-SNE 可能会耗费许多时间，所以后来提出 UMAP的方法，除了节省计算时间，效果甚至可能比 t-SNE 好。 UMAP 跟 t-SNE 一样是一种非线性降维方法，主要是引入 manifold 的概念，并且有下列几种假设 ：

1. The data is uniformly distributed on Riemannian manifold;
2. The Riemannian metric is locally constant (or can be approximated as such);
3. The manifold is locally connected.
   上面这些假设是为了让整个数据 manifold 形成一个类似 topology 的结构，相似的数据点在低维度空间当中会有相似的topology。
   t-SNE 或是 UMAP 貌似解决了 PCA 的问题，但是，这两种方法可以达到聚类良好的效果的同时仅保留了局部结构（local structure），忽略了全局结构（global structure）：
   只找出数据点远近的「关系」，但不考虑远近的「距离」
   以单细胞分析为例，不论t-SNE或是UMAP图，我们只能得知某群细胞是相似的聚类在一起，但无法从图中比较这群细胞跟其他细胞、或这群细胞里各点的远近关系。

• PHATE（Potential of Heat diffusion for Affinity-based Transition Embedding ）
  是2019 年底发表于Nature的文章，这个方法宣称可以同时保留局部结构（local structure）以及全局结构（global structure），并且对于隐含树状结构的数据会有更好的解释。