使用PHATE进行单细胞高维数据的可视化
简介
PHATE[1] (Potential of Heat-diffusion for Affinity-based Transition Embedding) 是Krishnaswamy实验室开发的一种用于可视化具有自然进程或轨迹的高维单细胞数据的工具。PHATE 使用一种新颖的概念框架来学习和可视化生物系统中固有的流形。其中，平滑过渡标志着细胞从一种状态到另一种状态的进展。PHATE生成一个低维嵌入，使用数据点之间的几何距离信息（即：样品表达谱之间的距离）来捕获数据集中的局部和全局非线性结构。
目前该文章发表在Nature Biotechnology顶级期刊上：Visualizing Structure and Transitions in High-Dimensional Biological Data. 2019\. Nature Biotechnology.[2]
以下是它的原理图：
数据简介
在本教程中，我们将演示如何使用 PHATE 来分析胚胎体 (EB) 在27天不同时间点分化的 31,000 个单细胞的数据。运行本教程大约需要 15 分钟（不包括 t-SNE 比较），或 25 分钟（包括比较）。
人胚胎体不同时间点的分化进程
首先，将低代 H1 hESCs 细胞培养在基质胶涂层的培养皿中，并在 DMEM/F12-N2B27 培养基中补充 FGF2成分。对于 EB 的形成，细胞用 Dispase 处理，解离成小团块并铺在非粘附板中，培养基中补充有 20% FBS。在为期 27 天的分化过程中，每隔 3 天收集一次样品，还包括未分化的 hESC 样本。通过 qPCR 验证这些 EB 培养物中关键胚层标记基因的表达。对于单细胞文库的构建，将EB 培养物分离，FACS 分选以去除双联体细胞和死细胞，并使用10x Genomics仪器生成 cDNA 文库，然后对其进行测序。
安装PHATE及其依赖包
使用pip直接进行安装
pip install --user phate
#安装MAGIC和scprep依赖包
pip install --user magic-impute
pip install --user scprep
下载源码进行安装
git clone --recursive git://github.com/KrishnaswamyLab/PHATE.git
cd PHATE/Python
python setup.py install --user
下载 10x 单细胞示例数据
EB示例数据集scRNAseq.zip可以从（https://data.mendeley.com/datasets/v6n743h5ng/） Mendelay Datasets 中下载使用。scRNA-seq文件夹中包含有5个子目录，每个子目录中包含有CellRanger处理生成的三个文件：barcodes.tsv，genes.tsv，和matrix.mtx。
# 导入所需的python包
import os
import scprep
download_path = os.path.expanduser("~")
print(download_path)
>>> /home/user
# 下载数据
if not os.path.isdir(os.path.join(download_path, "scRNAseq", "T0_1A")):
# need to download the data
scprep.io.download.download_and_extract_zip(
"https://md-datasets-public-files-prod.s3.eu-west-1.amazonaws.com/"
"5739738f-d4dd-49f7-b8d1-5841abdbeb1e",
download_path)
这是目录的结构：
使用scprep包导入数据并进行数据质控预处理
我们使用scprep包来进行数据的导入， 通过scprep.io.load_10X函数将会自动加载 10x scRNAseq 数据集到Pandas的 DataFrame 类型中。当然，我们也可以使用scprep.io.load_tsv()函数直接读取基因表达矩阵的TSV文件。
注意：默认情况下，scprep.io.load_10X函数使用 Pandas 的 SparseDataFrame 格式（参见 Pandas 文档）[3]加载 scRNA-seq 数据以最大化内存效率。但是，这将比加载dense矩阵的速度要慢一些。如果要加载dense矩阵，可以将sparse=False参数传递给load_10X函数。同时，我们还设置gene_labels = 'both'参数，这样就可以同时保留基因symbol和基因 ID 的信息。
# 导入所需的python包
import pandas as pd
import numpy as np
import phate
import scprep
# matplotlib settings for Jupyter notebooks only
%matplotlib inline
# 使用load_10X函数分别导入每个样本文件
sparse=True
T1 = scprep.io.load_10X(os.path.join(download_path, "scRNAseq", "T0_1A"), sparse=sparse, gene_labels='both')
T2 = scprep.io.load_10X(os.path.join(download_path, "scRNAseq", "T2_3B"), sparse=sparse, gene_labels='both')
T3 = scprep.io.load_10X(os.path.join(download_path, "scRNAseq", "T4_5C"), sparse=sparse, gene_labels='both')
T4 = scprep.io.load_10X(os.path.join(download_path, "scRNAseq", "T6_7D"), sparse=sparse, gene_labels='both')
T5 = scprep.io.load_10X(os.path.join(download_path, "scRNAseq", "T8_9E"), sparse=sparse, gene_labels='both')
# 查看数据信息
T1.head()
文库数据质量控制
首先，我们会过滤掉具有非常大或非常小的library size的细胞。对于这个数据集，文库的大小与样本有一定的相关性，因此我们会对每个样本单独进行过滤。这里，我们使用filter_library_size()函数过滤掉每个样本顶部和底部 20% 的细胞。
# 查看文库大小分布
scprep.plot.plot_library_size(T1, percentile=20)
plt.show()
# 分别对每个文库进行过滤
filtered_batches = []
for batch in [T1, T2, T3, T4, T5]:
batch = scprep.filter.filter_library_size(batch, percentile=20, keep_cells='above')
batch = scprep.filter.filter_library_size(batch, percentile=75, keep_cells='below')
filtered_batches.append(batch)
del T1, T2, T3, T4, T5 # removes objects from memory
合并所有样本的数据
接下来，我们将使用combine_batches()函数合并所有样本的数据，并构建一个表示每个样本时间点的向量。
EBT_counts, sample_labels = scprep.utils.combine_batches(
filtered_batches,
["Day 00-03", "Day 06-09", "Day 12-15", "Day 18-21", "Day 24-27"],
append_to_cell_names=True
)
del filtered_batches # removes objects from memory
EBT_counts.head()
去除稀有基因
接下来，我们会去除那些在 10 个或更少细胞中表达的基因。
EBT_counts = scprep.filter.filter_rare_genes(EBT_counts, min_cells=10)
数据归一化
为了校正文库大小的差异，我们使用library_size_normalize()函数将每个细胞除以其库大小，然后按文库大小的中值进行重新缩放。
EBT_counts = scprep.normalize.library_size_normalize(EBT_counts)
去除线粒体含量较高的死细胞
通常，死细胞中线粒体 RNA 的表达水平可能会高于正常活细胞。因此，我们通过消除平均线粒体 RNA 表达水平最高的细胞来去除可疑的死细胞。
首先让我们看看线粒体基因的分布。
mito_genes = scprep.select.get_gene_set(EBT_counts, starts_with="MT-") # Get all mitochondrial genes. There are 14, FYI.
scprep.plot.plot_gene_set_expression(EBT_counts, genes=mito_genes, percentile=90)
plt.show()
这里，我们看到在前 90% 以上，线粒体 RNA 的表达急剧增加。因此，我们将在后续的分析中删除这些细胞。
数据转换
这里，我们将使用scprep.transform.sqrt()函数进行平方根转换。
EBT_counts = scprep.transform.sqrt(EBT_counts)
使用PHATE进行低维嵌入降维可视化
首先，我们实例化一个 PHATE 估计器对象，用于将 PHATE 低维嵌入拟合到给定数据集中。接下来，使用fit和fit_transform函数生成低维嵌入。有关更多信息，请查看PHATE 阅读文档[4]。
这里，我们将使用默认参数，但可以调整以下参数（阅读我们在phate.readthedocs.io[5]上的文档以了解更多信息）：
knn：最近邻居的数量（默认值：5）。如果得到的 PHATE 低维嵌入看起来非常不连贯，我们可以增加此值（例如，增加到 20）。如果您的数据集非常大（例如 >100k 细胞），您还应该考虑增加该值。decay：Alpha 衰减值（默认值：15）。减少decay值可以增加knn图的连通性，增加decay值会降低图的连通性。通常情况下，很少需要调整该值。t：Number of times to power the operator（默认值：'auto'）。这相当于对数据进行平滑调整的量。默认情况下会自动选择，但如果您的低维嵌入缺乏结构，可以考虑增加它，或者如果结构看起来太紧凑，则减少它。gamma：Informational distance constant（默认值：1）。gamma=1gives the PHATE log potential, but other informational distances can be interesting. If most of the points seem concentrated in one section of the plot, you can trygamma=0。
这是应用 PHATE 的最简单方法：
# 实例化phate估计器对象
phate_operator = phate.PHATE(n_jobs=-2)
# 使用fit_transform函数进行低维嵌入拟合
Y_phate = phate_operator.fit_transform(EBT_counts)
Calculating PHATE...
Running PHATE on 16821 cells and 17845 genes.
Calculating graph and diffusion operator...
Calculating PCA...
Calculated PCA in 32.16 seconds.
Calculating KNN search...
Calculated KNN search in 14.15 seconds.
Calculating affinities...
Calculated affinities in 0.26 seconds.
Calculated graph and diffusion operator in 52.76 seconds.
Calculating landmark operator...
Calculating SVD...
Calculated SVD in 2.15 seconds.
Calculating KMeans...
Calculated KMeans in 36.05 seconds.
Calculated landmark operator in 40.53 seconds.
Calculating optimal t...
Automatically selected t = 21
Calculated optimal t in 2.98 seconds.
Calculating diffusion potential...
Calculated diffusion potential in 2.12 seconds.
Calculating metric MDS...
Calculated metric MDS in 13.67 seconds.
Calculated PHATE in 112.09 seconds.
接下来，我们可以使用scprep.plot.scatter2d()函数对 PHATE 低维嵌入后的结果进行可视化展示。
scprep.plot.scatter2d(Y_phate, c=sample_labels, figsize=(12,8), cmap="Spectral",
ticks=False, label_prefix="PHATE")
plt.show()
由于我们希望寻找出EB在不同时间点分化过程中的微妙结构，并且我们期望一些分化轨迹是稀疏的。因此，我们可以将knn值从默认的 5 减少，并将t值从自动检测的21减少（在上面的输出结果中）。对于单细胞 RNA-seq数据，如果您想寻找数据集中微妙结构的变化，可以尝试将knn值降低至 3 或 4，如果您有数十万个细胞，则可以尝试将该值增加到30 或 40。此处，我们还将alpha值减少到 15。
phate_operator.set_params(knn=4, decay=15, t=12)
# We could also create a new operator:
# phate_operator = phate.PHATE(knn=4, decay=15, t=12, n_jobs=-2)
Y_phate = phate_operator.fit_transform(EBT_counts)
Calculating PHATE...
Running PHATE on 16821 cells and 17845 genes.
Calculating graph and diffusion operator...
Calculating PCA...
Calculated PCA in 30.84 seconds.
Calculating KNN search...
Calculated KNN search in 17.06 seconds.
Calculating affinities...
Calculated affinities in 8.93 seconds.
Calculated graph and diffusion operator in 62.44 seconds.
Calculating landmark operator...
Calculating SVD...
Calculated SVD in 10.38 seconds.
Calculating KMeans...
Calculated KMeans in 35.15 seconds.
Calculated landmark operator in 47.79 seconds.
Calculating diffusion potential...
Calculated diffusion potential in 1.45 seconds.
Calculating metric MDS...
Calculated metric MDS in 7.24 seconds.
Calculated PHATE in 118.93 seconds.scprep.plot.scatter2d(Y_phate, c=sample_labels, figsize=(12,8), cmap="Spectral",
ticks=False, label_prefix="PHATE")
plt.show()
当然，我们还可以使用scprep.plot.scatter3d()函数进行三维可视化展示。
phate_operator.set_params(n_components=3)
Y_phate_3d = phate_operator.transform()
scprep.plot.scatter3d(Y_phate_3d, c=sample_labels, figsize=(8,6), cmap="Spectral",
ticks=False, label_prefix="PHATE")
plt.show()
我们甚至还可以生成一个3D旋转的动态gif图。
scprep.plot.rotate_scatter3d(Y_phate_3d, c=sample_labels,
figsize=(8,6), cmap="Spectral",
ticks=False, label_prefix="PHATE")
# to save as a gif:
# >>> scprep.plot.rotate_scatter3d(Y_phate_3d, c=sample_labels,
# ... figsize=(8,6), cmap="Spectral",
# ... ticks=False, label_prefix="PHATE", filename="phate.gif")
# to save as an mp4:
# >>> scprep.plot.rotate_scatter3d(Y_phate_3d, c=sample_labels,
# ... figsize=(8,6), cmap="Spectral",
# ... ticks=False, label_prefix="PHATE", filename="phate.mp4")
与其他降维可视化工具的比较
接下来，我们将对PCA和t-SNE降维可视化的结果进行比较。
import sklearn.decomposition # PCA
import sklearn.manifold # t-SNE
import time
start = time.time()
# PCA降维可视化
pca_operator = sklearn.decomposition.PCA(n_components=2)
Y_pca = pca_operator.fit_transform(np.array(EBT_counts))
end = time.time()
print("Embedded PCA in {:.2f} seconds.".format(end-start))
start = time.time()
pca_operator = sklearn.decomposition.PCA(n_components=100)
# tSNE降维可视化
tsne_operator = sklearn.manifold.TSNE(n_components=2)
Y_tsne = tsne_operator.fit_transform(pca_operator.fit_transform(np.array(EBT_counts)))
end = time.time()
print("Embedded t-SNE in {:.2f} seconds.".format(end-start))# plot everything
import matplotlib.pyplot as plt
fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(16, 5))
#plotting PCA
scprep.plot.scatter2d(Y_pca, label_prefix="PC", title="PCA",
c=sample_labels, ticks=False, cmap='Spectral', ax=ax1)
#plotting tSNE
scprep.plot.scatter2d(Y_tsne, label_prefix="t-SNE", title="t-SNE", legend=False,
c=sample_labels, ticks=False, cmap='Spectral', ax=ax2)
#plotting PHATE
scprep.plot.scatter2d(Y_phate, label_prefix="PHATE", title="PHATE", legend=False,
c=sample_labels, ticks=False, cmap='Spectral', ax=ax3)
plt.tight_layout()
plt.show()
可以看到，PHATE可以很好的展示出EB在不同分化时间点的连续变化状态，同时也能保留它们全局的整体结构。
参考资料
PHATE: https://github.com/KrishnaswamyLab/PHATE
[2]
Visualizing Structure and Transitions in High-Dimensional Biological Data. 2019. Nature Biotechnology.: https://www.nature.com/articles/s41587-019-0336-3
[3]
Pandas 文档: https://pandas.pydata.org/pandas-docs/stable/sparse.html
[4]
PHATE 阅读文档: https://phate.readthedocs.io/
[5]
phate.readthedocs.io: https://phate.readthedocs.io/
[6]
参考: https://nbviewer.org/github/KrishnaswamyLab/PHATE/blob/master/Python/tutorial/EmbryoidBody.ipynb