本文发表于《nature biotechnology》2019年12月
链接:https://www.nature.com/articles/s41587-019-0336-3
一、PHATE 目的
通过对复杂的生物数据降维,帮助研究人员理解数据结构。
二、PHATE 优点
- 无需假设数据服从何种分布,即无需先验分布。假设的先验分布对符合其假设的数据很有用,但也可能产生误导性结果,并且通常不适合用于提出假设或数据探索。
- 降低对噪声的敏感性,提高局部结构的准确性。注意这里说的是噪声不是异常值。
- 高效利用计算资源,可以处理大规模数据集。 t-SNE 等方法的计算包只是对算法的简单实现,没有对算法的计算进行优化,计算资源利用效率低,无法处理大规模数据集。
- 兼顾局部结构和全局结构。t-SNE 通常会为了保证数据的局部结构而破坏全局结构。PCA 牺牲局部结构展示群体结构,并且 PCA 可视化效果的好坏很依赖与角度的选择,并非 PC1&PC2 图总能展现出我们最想要的结构。下图展示了不同将降维算法的降维结果。可以看出 t-SNE 可以较好的展示数据的局部结构,各组分之间重叠的较少,但整体结构会被破坏,如图 b 中的紫色、前黄色、草绿色被割离出了整体;PCA 展示了全局结构,但小类之间存在明显的重合,并且其他角度下的结构难以展示;PHATE 可以有效地对数据去噪并保存数据局部和全局结构。
- 算法专门为可视化而设计。常见的降维方法在设计时并没有将可视化作为其参考指标之一,所以在可视化上效果不佳。
- 算法具有鲁棒性,对用户的参数配置不敏感。
三、PHATE 降维步骤
- 测定到的原始数据转化为 “样本N × 特征M” 的矩阵形式,如图 a。
- 以特征为维度,特征值为坐标,构建样本空间,计算样本之间的欧式距离。此距离代表了样本之间的差异,得到一个 的样本差异矩阵 ,其中每个元素代表了两个样本间的距离, ,如图 b。我们可以注意到,目前为止样本间差异的度量与常用方法并无二致。
- 使用作者设计的
核函数将样本的欧式距离矩阵转化为核函数矩阵
。进而转化为变换概率矩阵
。作者称
为亲和力矩阵,
为马尔可夫亲和力矩阵(Markov-normalized affinity matrix),如图 c。
- 通过冯诺依曼熵(Von Neumann entropy)确定最适状态转换次数 ,计算经过 次转换后的变换概率矩阵 ,如图 d。此时的矩阵不仅表现了局部信息、全局信息,同时还降低了噪声对数据的影响。
- 将变换概率矩阵 中的元素取 ,用于计算样本之间的扩散距离(文中称为信息距离,informational distances),得到扩散(信息)距离矩阵 ,如图 e。
- 使用多维标度法 (Multidimensional Scaling, MDS),根据扩散距离矩阵 ,将样本嵌入到低维空间中,如图 f。
关于上述内容的详细介绍见 补充说明 。
四、PHATE 效果
上图显示了在 5 个具有已知轨迹(a,d,e)和簇(b,c)结构的单细胞数据集上 PHATE 方法与其他 7 个方法的降维比较。PHATE为数据提供了干净且相对去噪的可视化,突出了局部和全局结构。这些分支中的许多都与已有文献验证的细胞类型或簇相一致。如图 a、b 体现了 PHATE 兼顾了局部(与t-SNE比较)、全局(与PCA比较)、去噪(lsomap比较)的特性。可以发现 t-SNE 倾向于将轨迹粉碎为簇,从而产生错误的印象,即数据包含自然簇。
五、PHATE 有效性
我们在研究降维聚类问题时,必定要考虑一个问题:聚成一类的样本是否因为我们想要性质的相似而被聚集? 如果聚类不是因为各阶段调控发育的关键基因的表达量,那么用聚类后的结果研究发育是不合理的。在聚类时,如果我们对特征不赋予权重,意味着所有特征被等价对待 。所以当我们处理细胞发育的转录组数据时,我们必须要证明聚类及变迁是因为与发育有关的基因的表达量发生改变。为此我们需要对基因表达量进行聚类,观察不同类别间,与发育有关的基因的表达量是否有显著差异。结果证明,类别间的差异确实与发育相关的基因的表达量有关。这一点从数据中未包含时间序列信息,但降维后样本的排布有明显的时间线中也可以证明,因为时间代表的就是发育。表达量的聚类热图如下所示 (来自附图7) 。
下图展示了对特征不加权和多种加权后的降维可视化结果。其中对 stem cell markers 的加权感觉效果最好,展示了细胞随时间的分化。
下图详细的展示了各种加权下与加权相关的基因的表达了的变化,基因的表达量多→少用颜色表示黄→紫表示。如 i 表示对细胞周期相关基因进行加权,其中 ki67、ps6 分别展示了 ki67、ps6 基因的表达量变化规律,可以发现两个基因的表达规律有明显的区别,阐述了相关基因在时空上的配合。
六、PHATE 应用
1. 研究物种进化
上图展示在文章附录中,展示了以 SNP 作为特征进行降维与可视化,可以发现相比 PCA,使用 PHATE 可以更有效划分群体。PHATE 善于解决既需要局部结构有需要全局结构的实际问题,如文中列举的发育分化、遗传进化等包含时间线的问题。
2. 研究细胞分化
文章详细阐述了使用 PHATE 分析人类 ES 细胞分化数据。
作者分析人胚胎干细胞(ES)分化为胚状体(EB)过程中的单细胞 RNA 测序 (scRNA-seq) 数据。测量了约 31,000 个细胞,在 27 天的分化时间过程中平均分布 (每隔3天收集一次样品)。
- 通过 PHATE 算法对 scRNA-seq 数据降维及可视化,并根据固定维度对嵌入到低维度的群体进行分类(如图 a.ii)得到 10 个分枝(如图 a.iii,不同分枝用不同颜色表示)
- 根据已有文献,作者将 ES -> EB 过程中的细胞分为 27 个类别,如图 a.i、b、c(图中是 30 个类,其中有 3 个类是本次分析新建立的,详见图 b 中红字表明的类名)。在图 b 的谱系树中,点是细胞类别,箭头表示类别的分化路线。谱系树详细 ES 细胞如何通过连续的中间状态来产生胚层衍生物。同时根据文献与 scRNA-seq 数据找出与此过程密切相关的 70 个基因的表达量数据。
- 将 70 个基因的表达量绘制在降维后的分布图中。可以直观的看出不同时间、不同类型的细胞受到不同基因的调控。
- 计算 70 个基因在 30 个类中的 EMD 值,聚类后绘制出热图 d。EMD 值用于比较两个分布之间的差异,分布差异越大,EMD值越大。 此处 EMD 值体现的是基因在类内表达量的频率分布直方图与在30个类中的总表达量的频率分布直方图之间的差异,EMD 值越大说明基因越有可能对此类细胞有调控作用。热图 d 横坐标是根据文献划分出的 30 个类及其所属的细胞大类,纵坐标为 70 个基因,颜色表示 EMD 值的大小,如基因 SOX10 在 NC 细胞中 EMD 较大,说明 SOX10 在 NC 中表达量分布情况与总分布相差较大,SOX10 可能对 NC 细胞有特殊的调控作用。
- EMD 展现的是两个分布之间的差异,但分布的差异是否能有与调控等价?调控的充分不必要条件是基因的表达量与细胞的分化具有很高的相关性,而通过 PHATE 和 EMD 热图找到是与 batch 有关的基因,并不能证明相关性。为此作者详细分析了在 batch iii、vii 里 EMD 值高而在其他 batch 中低的基因,根据热图 e 中 EMD 值及抗体可用性选取 ITGA4、F3、CD82 作为研究对象,如果根据选取基因的表达量对细胞排序后其他基因的表达量,可以与根据时间排序后其他基因的表达量相吻合,则能证明选取的基因与 batch 内细胞的发育有很高的相关性且可能是位于上游的分化启动基因。如图 f,其中相关系数为 Spearman 相关系数。(Pearson 相关系数只能计算 ITGA4 与其他基因之间的相关性,或 ITGA4 与时间之间的相关性,无法结合起来考虑,但Spearman可以)
补充说明
1. 核函数
核函数:一类用于简化计算高维空间样本间距离的函数。
首先我们要了解为什么要将数据从原有的维度空间投射到高维空间中。将低维数据投射到高维空间中,可以使数据分布的更为离散。如支持向量机 SVM 通过将线性不可分数据投射到高维空间中从而实现线性可分,如下图所示:
将数据投射到高维空间中可能会使样本间的距离被更好的度量,但高维空间中各样本距离计算是复杂的。如定义我们定义f是一个将样本从4维空间映射到16维空间的函数,设 ,有 ,使用内积度量样本间的距离,请计算样本 在 16 维空间中的距离 。
经计算得到:
= (1, 2, 3, 4, 2, 4, 6, 8, 3, 6, 9, 12, 4, 8, 12, 16)
= (25, 30, 35, 40, 30, 36, 42, 48, 35, 42, 49, 56, 40, 48, 56, 64)
= 25 + 60 + 105 + 160 + 60 + 144 + 252 + 384 + 105 + 252 + 441 + 672 + 160 + 384 + 672 + 1024 = 4900
是否我们能够找到一个函数,使我们无需将样本先投射到高维空间再计算距离,可以直接在低维空间中就可以计算距离?设函数 ,有 。发现 ,使用 可以简便计算高维空间中样本 之间的距离。核函数正是这样一类用于简化计算高维空间样本间距离的函数,所以核函数也称为 kernel trick。一般常用高斯核函数 计算样本在高维空间中的距离。本文使用的 α-decay 核函数是以高斯核函数为基础修改得到的。
2. 扩散映射与转换概率矩阵
度量两个样本之间差异大小的方法有很多,最基础的方法是计算两个样本间的欧式距离,用欧式距离的大小反映样本间的差异。此外还有用 pearson 相关系数等其他方法得到的值来表示样本间差异。扩散映射 (Diffusion Maps) 是于 2006 年提出的一种基于流形学习的非线性降维方法,其主要思想来自于动力系统。作为一种新的流形学习框架,扩散映射通过在扩散过程中尽可能地保持扩散距离来进行降维,即保持样本点的局部结构不变,通过局部关系定义全局关系,使样本点在低维空间中仍保持这种稳定的全局关系。运用的是图论的思想,假设所有样本之间均相互连接,样本可以随机转化为集合中的任意其他样本,转化的过程称之为扩散(diffusing)。扩散映射的思想即是从般含在核矩阵 中的局部几何信息中构建数据的全局几何信息。在扩散映射中,使用扩散距离来描述样本间差异。
扩散距离的计算方法:
- 选择计算距离的核函数 。从上面对核函数的介绍中我们知道,使用核函数可以更为准确的度量样本间的距离。
- 用核函数计算各样本之间的距离,得到核函数矩阵 , 。 描述了样本之间的局部连接情况。
- 构建变换概率矩阵 , 。 表示样本 变换成 的概率,变换概率的大小与距离正相关。如果只是将距离转化为概率,没有什么实际意义,只是换了种说法,下面的步骤是扩散映射的关键。
- 将变换概率矩阵
视为 Markov 状态转移矩阵,设矩阵
表示对
进行
次叠乘,其中元素
表示样本
经过
次转移最终变换为
的概率。设群体总共有 3 个样本,通过步骤 1-3 得到矩阵
。
对矩阵进行 2 次叠乘,有 ,其中 ,表示样本1->样本1->样本2的概率 + 样本1->样本2->样本2的概率 + 样本1->样本3->样本2的概率,即从样本 1 出发先变换成 1、2、3 然后在变换为样本 2 的概率,变换路径用矩阵表示为 。
对矩阵进行 3 次叠乘,有 ,其中
,表示 1→1→1→2 的概率 + 1→2→1→2 的概率 + 1→3→1→2 的概率 + 1→1→2→2 的概率 + 1→2→2→2 的概率 + 1→3→2→2 的概率 + 1→1→3→2 的概率 + 1→2→3→2 的概率 + 1→3→3→2 的概率,即从样本 1 出发先变换成 1、2、3,然后在此基础上变换成 1、2、3,最后变换为样本 2 的概率,变换路径用矩阵表示为 。
可以发现,随着叠乘的增加,1->2 的变换轨迹经过各样本的次数增加,变换的总概率受到群体中其他位点的影响越大,即受到样本 1、2 所处全局几何位置的影响越大,实现了在局部信息中引入了全局信息。同时,两个样本之间距离不再是由两个样本所决定,而是由全部样本共同参与决定,所以距离对噪声更具有鲁棒性。 - 下面计算扩散距离 ,其中 是变换次数。扩散距离使用的也是样本整体来度量两个样本间的距离。总结来说,扩散映射的思想是从局部信息中构建出全局信息。
3. 高斯核函数
高斯核函数:
以高斯核函数的值作为距离的度量时,参数
的取值反映了概率矩阵
在全局信息和局部信息之间的权衡。
设
,
当
时,
,
当
时,
,
当
时,
。
可以发现,
的减小
使
与
之间的差异显著扩大,
由
,样本x向邻近样本z转换的概率显著增加,而样本
向较远样本
转换显著降低,使样本变换的发生更趋向于邻近样本。而
的增加
使
与
之间的差异显著减少,
由
,样本
向邻近样本
转换的概率与向较远样本
转换的概率差异逐渐减少。所以当
较小时,样本间的变换主要发生最邻近样本间,处于稀疏采样区的样本将会随着转移的迭代而逐渐被孤立;当
较大时,样本与邻近样本发生变换的概率与较远样本相差不大,样本的变换不主要集中于最邻近样本间,随着
的增加,样本变换的主要发生范围逐渐扩大。可以试想,当
,样本
向其他样本转移的概率均相等,不再受到样本间距离的影响,样本间的转移概率不再体现局部信息,仅体现全局信息;当
,样本
有且仅向距离最近的样本点转移,样本间的转移概率仅体现局部信息,不再体现全局信息。所以
的取值对概率矩阵
至关重要,它反映了
在全局信息和局部信息之间的权衡。
在本文中,作者将 变成一个自适应参数, 值会随着样本点的变化而变化, 样本到最近邻近 个样本的距离均值(k-NN)。在这种算法下,处于密集区的样本的 值较小,其转移概率更多的反映局部信息;处于稀疏区的样本的 值较大,其转移概率更多的反映全局信息。这种算法相比于 值固定,可以在 次迭代后的变换概率矩阵中展现更多的数据信息。对于空间中样本密集区,局部信息十分丰富,如果 较大,空间中样本距离较近的点无法准确区分, 应当较小,如 ,若 ,则 ,相对于 而言, 之间的差异难以体现;对于空间中样本稀疏区,局部信息基本没有,如果 较小,概率矩阵 中稀疏区样本所在的行、列中将有大量元素为 0,矩阵对空间结构的记录效率大幅降低, 应当较大。总体来看自适应参数 使距离远的样本被拉近,距离近的样本被拉远。
对于高斯核函数,还需要考虑函数的拖尾问题。因为核函数图像是曲线的,随着距离的增加,函数值减少的速度不断降低,到达一定距离后,附近横坐标对应的函数值基本不再改变,造成概率矩阵 中信息的冗余。我们可以通过函数图像发现这一问题。设有核函数 ,其函数图像如下:
当 时核函数 等价于函数 ,随着带宽 (bandwidth) 的增加,拖尾现象逐渐加重,即达到一定距离后,后续距离的增加对函数值的影响不大。所以对于处在空间稀疏区的样本,我们既不希望 过小导致概率矩阵 中有大量元素为 0,也不希望 过大导致概率矩阵P中有大量元素值差异微小,这两种情况都会导致矩阵中包含大量的冗余信息。自适应参数的引入避免处在空间稀疏区样本的 过小,但可能会导致 过大,为了解决这一问题,作者引入了参数 、双向平均距离代替单向距离、对概率取对数三种方法。下面公式是作者对核函数的改进,涉及参数 α、双向平均距离代替单向距离两种方法,而概率取对数方法在计算扩散距离时使用。
观察 的函数图像可以发现随着 的增加,拖尾现象得到缓解,但 附近图像逐渐平坦,即 的值也反映出了概率矩阵 在全局信息和局部信息之间的权衡。使用双向距离既可以有效缓解拖尾现象,又可以包含更多的局部信息。
设样本 处在稀疏区,样本 处在密集区,样本 到 的距离相等且略小于到 。因为 的值较大,对距离变化不敏感, 与 的差异很小;但使用双向距离后, 的值较小,对距离变化敏感, 与 有一定差异。两者总和的双向距离可以有效缓解拖尾现象。另一方面, 与 的距离虽然相等,但 所处的空间环境是不同的,单向距离只能反映出发点的局部信息,无法反映到达点的局部信息。使用双向距离则可以反映出发点、到达点的局部信息,如 不等于 。作者实践发现,对于本文中介绍的所有数据集, 效果较好。PS:作者通过定性和定量的实验证明 PHATE 对于参数 具有鲁棒性。
我们知道随着距离 的增加,函数 值不断减小,并且减小的速度不断降低,使得距离较远的样本函数值一致,造成信息的冗余。如果我们能找到一个函数 ,其函数值随着自变量的减小不断减小,但减小的速度不断提高,将函数 套在函数 外,形成复合函数 ,使其随着 增加函数值减小的过程中,减小速度增加、不变或降低的不那么快。作者使用的 为对数函数,将发生 次转换后的转换概率矩阵 中的元素转化为 ,并计算扩散距离,进一步解决拖尾现象。为了与原始的扩散距离 相区别,作者将以 为基础计算出的扩散距离称之为势距离(Potential distance)。
总的来说,上述方法、改进等都是优化“距离”矩阵,使矩阵包含更少的冗余信息,记录更多有用的信息。
4. MDS 多维标度分析
MDS,multidimensional scaling,多维标度分析的算法具体使用有空再补充。
5. DEMaP 度量标准
作者认为先前量化可视化降维效果的度量标准衡量的仅是降维保持成对的距离或本地邻域的能力,而没有考虑降维方法保持全局结构的能力。从 t-SNE 等方法中我们可以看到,仅考虑局部结构是不足够的,全局结构对我们分析数据也十分重要。作者还提出:原始数据是存在噪音的,如何能够衡量出降维算法的去噪能力?如何定义噪声?
因此,作者提出了 DEMaP 度量标准,同时衡量降维方法保持局部、全局结构的能力以及去噪的能力。DEMaP 通过将无噪声数据上的测地距离与从噪声数据中提取的嵌入的欧几里德距离进行比较,从而封装了这些属性中的每一个。
DEMaP 的具体算法有空再补充。
PS
数据分布在三维的一个瑞士卷上,那么数据的本身是在一个二维矩形上,但通过折叠旋转后嵌入进了三维空间中。我们期望的对瑞士卷的降维结果是2维结构并且能够合理反映原空间中各点之间的相对距离。PCA是线性降维,什么是线性降维?线性降维是指使用线性变换的方式进行的降维方法。什么是线性变换?当对坐标系的变换符合一下两条性质时,我们就称它为线性变换:1. 直线在变换后仍保持为直线,不能有所弯曲(Lines remain lines)或者说原空间中的网格线在变换后仍然等距且平行;2. 原点必须固定。所有线性降维,都满足线性变换的两个性质,即新的坐标轴是通过原坐标轴平移、伸缩、旋转得到,而不能对原坐标轴扭曲。所以PCA的降维后的结果只能展示原空间中图形的各个截面,而无法将瑞士卷展开拉伸。
文章里的continual progression翻译过来是连续的过程,具体可以指胚胎干细胞随时间分化过程。
启发 & 应用
- 扩散距离是否可以替代我们实验室中常用的一些度量距离的方法,如欧式距离、pearson相关系数等,用于自学习等。
- 扩散距离矩阵是否可以替代kinship从而更好的度量个体之间的相关性。
- 利用 PHATE 是否可以判断植株的长势分化,是否有按照我们期望的方向发展。
- 类比细胞的发育,PHATE 可以尝试用于解释物种的遗传进化路线。
- 通过 batch 研究基因在时空上的配合。