文章:PHATE:Visualizing structure and transitions in high-dimensional biological data

1 篇文章 0 订阅
订阅专栏
1 篇文章 0 订阅
订阅专栏

本文发表于《nature biotechnology》2019年12月。文章整体质量颇高,介绍了一种研究细胞发育的算法——PHATE。虽然作者提出的新算法是 扩散映射多维标度分析 的结合改良版,但作者这种将 恰当的方法 应用于重要的实际问题上的能力是十分值得学习的。要想实现这一目标需要一定的数学功底。

文章链接:https://www.nature.com/articles/s41587-019-0336-3
软件链接: https://github.com/KrishnaswamyLab/PHATE
软件使用介绍: Tutoriall , Introduction, Demonstration

一、PHATE 目的

通过对复杂的生物数据降维,帮助研究人员理解数据结构。

二、PHATE 优点

  1. 无需假设数据服从何种分布,即无需先验分布。假设的先验分布对符合其假设的数据很有用,但也可能产生误导性结果,并且通常不适合用于提出假设或数据探索。
  2. 降低对噪声的敏感性,提高局部结构的准确性。注意这里说的是噪声不是异常值。
  3. 高效利用计算资源,可以处理大规模数据集。 t-SNE 等方法的计算包只是对算法的简单实现,没有对算法的计算进行优化,计算资源利用效率低,无法处理大规模数据集。
  4. 兼顾局部结构和全局结构。t-SNE 通常会为了保证数据的局部结构而破坏全局结构。PCA 牺牲局部结构展示群体结构,并且 PCA 可视化效果的好坏很依赖与角度的选择,并非 PC1&PC2 图总能展现出我们最想要的结构。下图展示了不同将降维算法的降维结果。可以看出 t-SNE 可以较好的展示数据的局部结构,各组分之间重叠的较少,但整体结构会被破坏,如图 b 中的紫色、前黄色、草绿色被割离出了整体;PCA 展示了全局结构,但小类之间存在明显的重合,并且其他角度下的结构难以展示;PHATE 可以有效地对数据去噪并保存数据局部和全局结构。

在这里插入图片描述

  1. 算法专门为可视化而设计。常见的降维方法在设计时并没有将可视化作为其参考指标之一,所以在可视化上效果不佳。
  2. 算法具有鲁棒性,对用户的参数配置不敏感。

三、PHATE 降维步骤

在这里插入图片描述

  1. 测定到的原始数据转化为 “样本N × 特征M” 的矩阵形式,如图 a。
  2. 以特征为维度,特征值为坐标,构建样本空间,计算样本之间的欧式距离。此距离代表了样本之间的差异,得到一个 N × N N×N N×N 的样本差异矩阵 Q Q Q,其中每个元素代表了两个样本间的距离, Q i j = ∣ ∣ c e l l i − c e l l j ∣ ∣ 2 Q_{ij} = \sqrt{||cell_i-cell_j||^2} Qij=cellicellj2 ,如图 b。我们可以注意到,目前为止样本间差异的度量与常用方法并无二致。
  3. 使用作者设计的 α − d e c a y α-decay αdecay 核函数将样本的欧式距离矩阵转化为核函数矩阵 K K K。进而转化为变换概率矩阵 P P P。作者称 K K K 为亲和力矩阵, P P P 为马尔可夫亲和力矩阵(Markov-normalized affinity matrix),如图 c。
    在这里插入图片描述
  4. 通过冯诺依曼熵(Von Neumann entropy)确定最适状态转换次数 t t t,计算经过 t t t 次转换后的变换概率矩阵 P t P^t Pt,如图 d。此时的矩阵不仅表现了局部信息、全局信息,同时还降低了噪声对数据的影响。
  5. 将变换概率矩阵 P t P^t Pt 中的元素取 l o g log log,用于计算样本之间的扩散距离(文中称为信息距离,informational distances),得到扩散(信息)距离矩阵 D D D,如图 e。
  6. 使用多维标度法 (Multidimensional Scaling, MDS),根据扩散距离矩阵 D D D,将样本嵌入到低维空间中,如图 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天收集一次样品)。

  1. 通过 PHATE 算法对 scRNA-seq 数据降维及可视化,并根据固定维度对嵌入到低维度的群体进行分类(如图 a.ii)得到 10 个分枝(如图 a.iii,不同分枝用不同颜色表示)
  2. 根据已有文献,作者将 ES -> EB 过程中的细胞分为 27 个类别,如图 a.i、b、c(图中是 30 个类,其中有 3 个类是本次分析新建立的,详见图 b 中红字表明的类名)。在图 b 的谱系树中,点是细胞类别,箭头表示类别的分化路线。谱系树详细 ES 细胞如何通过连续的中间状态来产生胚层衍生物。同时根据文献与 scRNA-seq 数据找出与此过程密切相关的 70 个基因的表达量数据。
  3. 将 70 个基因的表达量绘制在降维后的分布图中。可以直观的看出不同时间、不同类型的细胞受到不同基因的调控。
  4. 计算 70 个基因在 30 个类中的 EMD 值,聚类后绘制出热图 d。EMD 值用于比较两个分布之间的差异,分布差异越大,EMD值越大。 此处 EMD 值体现的是基因在类内表达量的频率分布直方图与在30个类中的总表达量的频率分布直方图之间的差异,EMD 值越大说明基因越有可能对此类细胞有调控作用。热图 d 横坐标是根据文献划分出的 30 个类及其所属的细胞大类,纵坐标为 70 个基因,颜色表示 EMD 值的大小,如基因 SOX10 在 NC 细胞中 EMD 较大,说明 SOX10 在 NC 中表达量分布情况与总分布相差较大,SOX10 可能对 NC 细胞有特殊的调控作用。
  5. 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维空间的函数,设 x = ( x 1 , x 2 , x 3 , x 4 ) x = (x_1, x_2, x_3, x_4) x=(x1,x2,x3,x4),有 f ( x ) = ( x 1 x 1 , x 1 x 2 , x 1 x 3 , x 1 x 4 , x 2 x 1 , x 2 x 2 , x 2 x 3 , x 2 x 4 , x 3 x 1 , x 3 x 2 , x 3 x 3 , x 3 x 4 , x 4 x 1 , x 4 x 2 , x 4 x 3 , x 4 x 4 ) f(x) = (x_1x_1, x_1x_2, x_1x_3, x_1x_4, x_2x_1, x_2x_2, x_2x_3, x_2x_4, x_3x_1, x_3x_2, x_3x_3, x_3x_4, x_4x_1, x_4x_2, x_4x_3, x_4x_4) f(x)=(x1x1,x1x2,x1x3,x1x4,x2x1,x2x2,x2x3,x2x4,x3x1,x3x2,x3x3,x3x4,x4x1,x4x2,x4x3,x4x4),使用内积度量样本间的距离,请计算样本 a = ( 1 , 2 , 3 , 4 ) , b = ( 5 , 6 , 7 , 8 ) a = (1, 2, 3, 4), b = (5, 6, 7, 8) a=(1,2,3,4),b=(5,6,7,8) 在 16 维空间中的距离 < f ( a ) , f ( b ) > <f(a), f(b)> <f(a),f(b)>

经计算得到:
f ( a ) f(a) f(a) = (1, 2, 3, 4, 2, 4, 6, 8, 3, 6, 9, 12, 4, 8, 12, 16)
f ( b ) f(b) f(b) = (25, 30, 35, 40, 30, 36, 42, 48, 35, 42, 49, 56, 40, 48, 56, 64)
< f ( a ) , f ( b ) > <f(a), f(b)> <f(a),f(b)> = 25 + 60 + 105 + 160 + 60 + 144 + 252 + 384 + 105 + 252 + 441 + 672 + 160 + 384 + 672 + 1024 = 4900

是否我们能够找到一个函数,使我们无需将样本先投射到高维空间再计算距离,可以直接在低维空间中就可以计算距离?设函数 k ( x , y ) = ( < x , y > ) 2 k(x, y) = (<x, y>)^2 k(x,y)=(<x,y>)2,有 k ( a , b ) = ( 5 + 12 + 21 + 32 ) 2 = 7 0 2 = 4900 k(a, b) = (5+12+21+32)^2 = 70^2 = 4900 k(a,b)=(5+12+21+32)2=702=4900 。发现 k ( a , b ) = < f ( a ) , f ( b ) > k(a, b) = <f(a), f(b)> k(a,b)=<f(a),f(b)>,使用 k ( x , y ) k(x, y) k(x,y) 可以简便计算高维空间中样本 x , y x, y x,y 之间的距离。核函数正是这样一类用于简化计算高维空间样本间距离的函数,所以核函数也称为 kernel trick。一般常用高斯核函数 K ( x , y ) = e x p ( − ( x − y ) 2 2 σ 2 ) K(x, y)=exp(-\dfrac{(x-y)^2}{2σ^2}) K(x,y)=exp(2σ2(xy)2) 计算样本在高维空间中的距离。本文使用的 α-decay 核函数是以高斯核函数为基础修改得到的。

2. 扩散映射与转换概率矩阵

度量两个样本之间差异大小的方法有很多,最基础的方法是计算两个样本间的欧式距离,用欧式距离的大小反映样本间的差异。此外还有用 pearson 相关系数等其他方法得到的值来表示样本间差异。扩散映射 (Diffusion Maps) 是于 2006 年提出的一种基于流形学习的非线性降维方法,其主要思想来自于动力系统。作为一种新的流形学习框架,扩散映射通过在扩散过程中尽可能地保持扩散距离来进行降维,即保持样本点的局部结构不变,通过局部关系定义全局关系,使样本点在低维空间中仍保持这种稳定的全局关系。运用的是图论的思想,假设所有样本之间均相互连接,样本可以随机转化为集合中的任意其他样本,转化的过程称之为扩散(diffusing)。扩散映射的思想即是从般含在核矩阵 K K K 中的局部几何信息中构建数据的全局几何信息。在扩散映射中,使用扩散距离来描述样本间差异。

扩散距离的计算方法:

  1. 选择计算距离的核函数 k ( x , y ) k(x, y) k(x,y)。从上面对核函数的介绍中我们知道,使用核函数可以更为准确的度量样本间的距离。
  2. 用核函数计算各样本之间的距离,得到核函数矩阵 K K K K i , j = k ( x i , x j ) K_{i, j}=k(x_i, x_j) Ki,j=k(xi,xj) K K K描述了样本之间的局部连接情况。
  3. 构建变换概率矩阵 P P P P i , j = K i , j ∑ i = 1 K i , j P_{i, j} = \dfrac{K_{i, j}}{\displaystyle\sum_{i=1}K_{i, j}} Pi,j=i=1Ki,jKi,j P i , j P_{i, j} Pi,j 表示样本 i i i 变换成 j j j 的概率,变换概率的大小与距离正相关。如果只是将距离转化为概率,没有什么实际意义,只是换了种说法,下面的步骤是扩散映射的关键。
  4. 将变换概率矩阵 P P P 视为 Markov 状态转移矩阵,设矩阵 P n P^n Pn 表示对 P P P 进行 n n n 次叠乘,其中元素 p i , j n p^n_{i, j} pi,jn 表示样本 i i i 经过 n n n 次转移最终变换为 j j j 的概率。设群体总共有 3 个样本,通过步骤 1-3 得到矩阵 P = [ p 11 p 12 p 13 p 21 p 22 p 23 p 31 p 32 p 33 ] P=\begin{bmatrix}p_{11}&p_{12}&p_{13}\\p_{21}&p_{22}&p_{23}\\p_{31}&p_{32}&p_{33} \end{bmatrix} P=p11p21p31p12p22p32p13p23p33
    对矩阵进行 2 次叠乘,有 P 2 = [ p 11 p 12 p 13 p 21 p 22 p 23 p 31 p 32 p 33 ] [ p 11 p 12 p 13 p 21 p 22 p 23 p 31 p 32 p 33 ] = [ p 11 p 11 + p 12 p 21 + p 13 p 31 p 11 p 12 + p 12 p 22 + p 13 p 32 p 11 p 13 + p 12 p 23 + p 13 p 33 p 21 p 11 + p 22 p 21 + p 23 p 31 p 21 p 12 + p 22 p 22 + p 23 p 32 p 21 p 13 + p 22 p 23 + p 23 p 33 p 31 p 11 + p 32 p 21 + p 33 p 31 p 31 p 12 + p 32 p 22 + p 33 p 32 p 31 p 13 + p 32 p 23 + p 33 p 33 ] P^2=\begin{bmatrix}p_{11}&p_{12}&p_{13}\\p_{21}&p_{22}&p_{23}\\p_{31}&p_{32}&p_{33} \end{bmatrix}\begin{bmatrix}p_{11}&p_{12}&p_{13}\\p_{21}&p_{22}&p_{23}\\p_{31}&p_{32}&p_{33} \end{bmatrix}=\begin{bmatrix}p_{11}p_{11}+p_{12}p_{21}+p_{13}p_{31}&p_{11}p_{12}+p_{12}p_{22}+p_{13}p_{32}&p_{11}p_{13}+p_{12}p_{23}+p_{13}p_{33}\\p_{21}p_{11}+p_{22}p_{21}+p_{23}p_{31}&p_{21}p_{12}+p_{22}p_{22}+p_{23}p_{32}&p_{21}p_{13}+p_{22}p_{23}+p_{23}p_{33}\\p_{31}p_{11}+p_{32}p_{21}+p_{33}p_{31}&p_{31}p_{12}+p_{32}p_{22}+p_{33}p_{32}&p_{31}p_{13}+p_{32}p_{23}+p_{33}p_{33} \end{bmatrix} P2=p11p21p31p12p22p32p13p23p33p11p21p31p12p22p32p13p23p33=p11p11+p12p21+p13p31p21p11+p22p21+p23p31p31p11+p32p21+p33p31p11p12+p12p22+p13p32p21p12+p22p22+p23p32p31p12+p32p22+p33p32p11p13+p12p23+p13p33p21p13+p22p23+p23p33p31p13+p32p23+p33p33,其中 P 1 , 2 2 = p 11 p 12 + p 12 p 22 + p 13 p 32 P^2_{1,2} = p_{11}p_{12}+p_{12}p_{22}+p_{13}p_{32} P1,22=p11p12+p12p22+p13p32,表示样本1->样本1->样本2的概率 + 样本1->样本2->样本2的概率 + 样本1->样本3->样本2的概率,即从样本 1 出发先变换成 1、2、3 然后在变换为样本 2 的概率,变换路径用矩阵表示为 [ 1 1 2 1 2 2 1 3 2 ] \begin{bmatrix}1&1&2\\1&2&2\\1&3&2\end{bmatrix} 111123222
    对矩阵进行 3 次叠乘,有 P 3 = [ p 11 p 12 p 13 p 21 p 22 p 23 p 31 p 32 p 33 ] [ p 11 p 12 p 13 p 21 p 22 p 23 p 31 p 32 p 33 ] [ p 11 p 12 p 13 p 21 p 22 p 23 p 31 p 32 p 33 ] P^3=\begin{bmatrix}p_{11}&p_{12}&p_{13}\\p_{21}&p_{22}&p_{23}\\p_{31}&p_{32}&p_{33} \end{bmatrix}\begin{bmatrix}p_{11}&p_{12}&p_{13}\\p_{21}&p_{22}&p_{23}\\p_{31}&p_{32}&p_{33} \end{bmatrix}\begin{bmatrix}p_{11}&p_{12}&p_{13}\\p_{21}&p_{22}&p_{23}\\p_{31}&p_{32}&p_{33} \end{bmatrix} P3=p11p21p31p12p22p32p13p23p33p11p21p31p12p22p32p13p23p33p11p21p31p12p22p32p13p23p33,其中
    P 1 , 2 3 = ( p 11 p 11 + p 12 p 21 + p 13 p 31 ) p 12 + ( p 11 p 12 + p 12 p 22 + p 13 p 32 ) p 22 + ( p 11 p 13 + p 12 p 23 + p 13 p 33 ) p 32 P^3_{1,2} = (p_{11}p_{11}+p_{12}p_{21}+p_{13}p_{31})p_{12}+(p_{11}p_{12}+p_{12}p_{22}+p_{13}p_{32})p_{22}+(p_{11}p_{13}+p_{12}p_{23}+p_{13}p_{33})p_{32} P1,23=(p11p11+p12p21+p13p31)p12+(p11p12+p12p22+p13p32)p22+(p11p13+p12p23+p13p33)p32
    = p 11 p 11 p 12 + p 12 p 21 p 12 + p 13 p 31 p 12 + p 11 p 12 p 22 + p 12 p 22 p 22 + p 13 p 32 p 22 + p 11 p 13 p 32 + p 12 p 23 p 32 + p 13 p 33 p 32 =p_{11}p_{11}p_{12}+p_{12}p_{21}p_{12}+p_{13}p_{31}p_{12}+p_{11}p_{12}p_{22}+p_{12}p_{22}p_{22}+p_{13}p_{32}p_{22}+p_{11}p_{13}p_{32}+p_{12}p_{23}p_{32}+p_{13}p_{33}p_{32} =p11p11p12+p12p21p12+p13p31p12+p11p12p22+p12p22p22+p13p32p22+p11p13p32+p12p23p32+p13p33p32,表示 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 1 1 2 1 1 2 2 1 1 3 2 1 2 1 2 1 2 2 2 1 2 3 2 1 3 1 2 1 3 2 2 1 3 3 2 ] \begin{bmatrix}1&1&1&2\\1&1&2&2\\1&1&3&2\\1&2&1&2\\1&2&2&2\\1&2&3&2\\1&3&1&2\\1&3&2&2\\1&3&3&2\end{bmatrix} 111111111111222333123123123222222222
    可以发现,随着叠乘的增加,1->2 的变换轨迹经过各样本的次数增加,变换的总概率受到群体中其他位点的影响越大,即受到样本 1、2 所处全局几何位置的影响越大,实现了在局部信息中引入了全局信息。同时,两个样本之间距离不再是由两个样本所决定,而是由全部样本共同参与决定,所以距离对噪声更具有鲁棒性。
  5. 下面计算扩散距离 D n ( i , j ) 2 = ∑ u ∈ X ∣ P i , u n − P j , u n ∣ 2 D_n(i, j)^2=\displaystyle\sum_{u\in{X}}|P^n_{i, u}-P^n_{j, u}|^2 Dn(i,j)2=uXPi,unPj,un2,其中 n n n 是变换次数。扩散距离使用的也是样本整体来度量两个样本间的距离。总结来说,扩散映射的思想是从局部信息中构建出全局信息。

3. 高斯核函数

高斯核函数: K ( x , y ) = e x p ( − ( x − y ) 2 / ε ) K(x, y)=exp(-(x-y)^2/ε) K(x,y)=exp((xy)2/ε)

以高斯核函数的值作为距离的度量时,参数 ε ε ε 的取值反映了概率矩阵 P P P 在全局信息和局部信息之间的权衡。
∣ ∣ x − y ∣ ∣ 2 = 10 , ∣ ∣ x − z ∣ ∣ 2 = 5 ||x-y||_2 = 10, ||x-z||_2 = 5 xy2=10,xz2=5
ε = 5 ε=5 ε=5 时, K ( x , y ) = e − 2 , K ( x , z ) = e − 1 , K ( x , z ) / K ( x , y ) = e K(x, y)=e^{-2}, K(x, z)=e^{-1}, K(x, z)/K(x, y)=e K(x,y)=e2,K(x,z)=e1,K(x,z)/K(x,y)=e
ε = 1 ε=1 ε=1 时, K ( x , y ) = e − 10 , K ( x , z ) = e − 5 , K ( x , z ) / K ( x , y ) = e 5 K(x, y)=e^{-10}, K(x, z)=e^{-5}, K(x, z)/K(x, y)=e^5 K(x,y)=e10,K(x,z)=e5,K(x,z)/K(x,y)=e5
ε = 10 ε=10 ε=10 时, K ( x , y ) = e − 1 , K ( x , z ) = e − 1 / 2 , K ( x , z ) / K ( x , y ) = e 1 / 2 K(x, y)=e^{-1}, K(x, z)=e^{-1/2}, K(x, z)/K(x, y)=e^{1/2} K(x,y)=e1,K(x,z)=e1/2,K(x,z)/K(x,y)=e1/2
可以发现, ε ε ε 的减小 ( 5 → 1 ) (5\rarr1) (51) 使 K ( x , z ) K(x, z) K(x,z) K ( x , y ) K(x, y) K(x,y) 之间的差异显著扩大, K ( x , z ) / K ( x , y ) = P ( x , z ) / P ( x , y ) K(x, z)/K(x, y) = P(x, z)/P(x, y) K(x,z)/K(x,y)=P(x,z)/P(x,y) e → e 5 e\rarr e^5 ee5,样本x向邻近样本z转换的概率显著增加,而样本 x x x 向较远样本 y y y 转换显著降低,使样本变换的发生更趋向于邻近样本。而 ε ε ε 的增加 ( 1 → 10 ) (1\rarr10) (110) 使 K ( x , z ) K(x, z) K(x,z) K ( x , y ) K(x, y) K(x,y) 之间的差异显著减少, K ( x , z ) / K ( x , y ) = P ( x , z ) / P ( x , y ) K(x, z)/K(x, y) = P(x, z)/P(x, y) K(x,z)/K(x,y)=P(x,z)/P(x,y) e 5 → e 1 / 2 e^5\rarr e^{1/2} e5e1/2,样本 x x x 向邻近样本 z z z 转换的概率与向较远样本 y y y 转换的概率差异逐渐减少。所以当 ε ε ε 较小时,样本间的变换主要发生最邻近样本间,处于稀疏采样区的样本将会随着转移的迭代而逐渐被孤立;当 ε ε ε 较大时,样本与邻近样本发生变换的概率与较远样本相差不大,样本的变换不主要集中于最邻近样本间,随着 ε ε ε 的增加,样本变换的主要发生范围逐渐扩大。可以试想,当 ε → ∞ , K ( x , z ) / K ( x , y ) = e x p ( − 1 ε ∗ [ ( x − z ) 2 − ( x − y ) 2 ) ] → e x p ( 5 / ε ) → 1 ε\rarr \infin, K(x, z)/K(x, y)=exp(-\dfrac{1}{ε} * [(x-z)^2-(x-y)^2)]\rarr exp(5/ε)\rarr1 ε,K(x,z)/K(x,y)=exp(ε1[(xz)2(xy)2)]exp(5/ε)1,样本 x x x 向其他样本转移的概率均相等,不再受到样本间距离的影响,样本间的转移概率不再体现局部信息,仅体现全局信息;当 ε → 0 , K ( x , z ) / K ( x , y ) = e x p ( 5 / ε ) → e x p ( ∞ ) → ∞ ε\rarr0, K(x, z)/K(x, y)=exp(5/ε)\rarr exp( \infin) \rarr \infin ε0,K(x,z)/K(x,y)=exp(5/ε)exp(),样本 x x x 有且仅向距离最近的样本点转移,样本间的转移概率仅体现局部信息,不再体现全局信息。所以 ε ε ε 的取值对概率矩阵 P P P 至关重要,它反映了 P P P 在全局信息和局部信息之间的权衡。

在本文中,作者将 ε ε ε 变成一个自适应参数, ε ε ε 值会随着样本点的变化而变化, ε = ε= ε= 样本到最近邻近 k k k 个样本的距离均值(k-NN)。在这种算法下,处于密集区的样本的 ε ε ε 值较小,其转移概率更多的反映局部信息;处于稀疏区的样本的 ε ε ε 值较大,其转移概率更多的反映全局信息。这种算法相比于 ε ε ε 值固定,可以在 n n n 次迭代后的变换概率矩阵中展现更多的数据信息。对于空间中样本密集区,局部信息十分丰富,如果 ε ε ε 较大,空间中样本距离较近的点无法准确区分, ε ε ε 应当较小,如 ∣ ∣ x − y ∣ ∣ 2 = 10 , ∣ ∣ x − z ∣ ∣ 2 = 5 ||x-y||_2 = 10, ||x-z||_2 = 5 xy2=10,xz2=5,若 ε = 50 ε=50 ε=50,则 K ( x , z ) / K ( x , y ) = e 1 / 10 K(x, z)/K(x, y)=e^{1/10} K(x,z)/K(x,y)=e1/10,相对于 x x x 而言, y , z y, z y,z之间的差异难以体现;对于空间中样本稀疏区,局部信息基本没有,如果 ε ε ε 较小,概率矩阵 P P P 中稀疏区样本所在的行、列中将有大量元素为 0,矩阵对空间结构的记录效率大幅降低, ε ε ε 应当较大。总体来看自适应参数 ε ε ε 使距离远的样本被拉近,距离近的样本被拉远。

对于高斯核函数,还需要考虑函数的拖尾问题。因为核函数图像是曲线的,随着距离的增加,函数值减少的速度不断降低,到达一定距离后,附近横坐标对应的函数值基本不再改变,造成概率矩阵 P P P 中信息的冗余。我们可以通过函数图像发现这一问题。设有核函数 K α , σ ( x ) = e x p ( − ( ∣ x ∣ / σ ) α ) K_{α, σ}(x) = exp(-(|x|/σ)^α) Kα,σ(x)=exp((x/σ)α),其函数图像如下:

加粗样式

α = 2 α=2 α=2 时核函数 K K K 等价于函数 e x p ( − ( x − y ) 2 / ε ) exp(-(x-y)^2/ε) exp((xy)2/ε),随着带宽 (bandwidth) σ σ σ 的增加,拖尾现象逐渐加重,即达到一定距离后,后续距离的增加对函数值的影响不大。所以对于处在空间稀疏区的样本,我们既不希望 ε ε ε 过小导致概率矩阵 P P P 中有大量元素为 0,也不希望 ε ε ε 过大导致概率矩阵P中有大量元素值差异微小,这两种情况都会导致矩阵中包含大量的冗余信息。自适应参数的引入避免处在空间稀疏区样本的 ε ε ε 过小,但可能会导致 ε ε ε 过大,为了解决这一问题,作者引入了参数 α α α、双向平均距离代替单向距离、对概率取对数三种方法。下面公式是作者对核函数的改进,涉及参数 α、双向平均距离代替单向距离两种方法,而概率取对数方法在计算扩散距离时使用。

观察 K α , σ ( x ) = e x p ( − ( ∣ x ∣ / σ ) α ) K_{α, σ}(x) = exp(-(|x|/σ)^α) Kα,σ(x)=exp((x/σ)α) 的函数图像可以发现随着 α α α 的增加,拖尾现象得到缓解,但 x = 0 x=0 x=0 附近图像逐渐平坦,即 α α α 的值也反映出了概率矩阵 P P P 在全局信息和局部信息之间的权衡。使用双向距离既可以有效缓解拖尾现象,又可以包含更多的局部信息。

设样本 w , x w, x w,x处在稀疏区,样本 y , z y, z y,z 处在密集区,样本 w w w x , y x, y x,y 的距离相等且略小于到 z z z。因为 ε k ( w ) ε_k(w) εk(w) 的值较大,对距离变化不敏感, K ( w , y ) K(w, y) K(w,y) K ( w , z ) K(w, z) K(w,z) 的差异很小;但使用双向距离后, ε k ( y ) , ε k ( z ) ε_k(y), ε_k(z) εk(y),εk(z) 的值较小,对距离变化敏感, K ( y , w ) K(y, w) K(y,w) K ( z , w ) K(z, w) K(z,w) 有一定差异。两者总和的双向距离可以有效缓解拖尾现象。另一方面, w w w x , y x, y x,y 的距离虽然相等,但 x , y x, y x,y 所处的空间环境是不同的,单向距离只能反映出发点的局部信息,无法反映到达点的局部信息。使用双向距离则可以反映出发点、到达点的局部信息,如 K k , α ( w , x ) K_{k, α}(w, x) Kk,α(w,x) 不等于 K k , α ( w , y ) K_{k, α}(w, y) Kk,α(w,y)。作者实践发现,对于本文中介绍的所有数据集, k = 5 , α = 10 k=5, α=10 k=5,α=10 效果较好。PS:作者通过定性和定量的实验证明 PHATE 对于参数 k , α k, α k,α 具有鲁棒性。

我们知道随着距离 x x x 的增加,函数 K α , σ ( x ) K_{α, σ}(x) Kα,σ(x) 值不断减小,并且减小的速度不断降低,使得距离较远的样本函数值一致,造成信息的冗余。如果我们能找到一个函数 g g g,其函数值随着自变量的减小不断减小,但减小的速度不断提高,将函数 g g g 套在函数 K α , σ K_{α, σ} Kα,σ 外,形成复合函数 g ( K α , σ ( x ) ) g(K_{α, σ}(x)) g(Kα,σ(x)),使其随着 x x x 增加函数值减小的过程中,减小速度增加、不变或降低的不那么快。作者使用的 g g g 为对数函数,将发生 t t t 次转换后的转换概率矩阵 P t P^t Pt 中的元素转化为 l o g P i , j t logP^t_{i, j} logPi,jt,并计算扩散距离,进一步解决拖尾现象。为了与原始的扩散距离 D n ( i , j ) 2 = ∑ u ∈ X ∣ P i , u n − P j , u n ∣ 2 D_n(i, j)^2=\displaystyle\sum_{u\in{X}}|P^n_{i, u}-P^n_{j, u}|^2 Dn(i,j)2=uXPi,unPj,un2相区别,作者将以 l o g P i , j t logP^t_{i,j} logPi,jt 为基础计算出的扩散距离称之为势距离(Potential distance)。

总的来说,上述方法、改进等都是优化“距离”矩阵,使矩阵包含更少的冗余信息,记录更多有用的信息。

4. MDS 多维标度分析

参考: https://zhuanlan.zhihu.com/p/50715681

MDS,multidimensional scaling,多维标度分析,是一种经典的数据降维方法,同时也是数据可视化的一种手段。方法最初提出时所要解决的问题为:当我们仅能获得物体之间的距离 (相似性) 矩阵时,如何由此来重构它们的欧几里德坐标

MDS 常用于通过样本间的距离矩阵,将样本从高维空间重构到低维空间中,实现数据的降维。 用机器学习角度的话说就是:MDS 的目标函数为 “低维空间下的距离矩阵与高维空间下的距离矩阵的差异”,最小化目标函数得到样本在低维空间中的坐标。 如果将高维数据重构到 2、3 维的空间中,则也可以视为一种数据可视化的手段。本研究是将细胞的高维转录组数据重构到 2、3 维空间中,实现单细胞数据的可视化,直观反映胚发育过程中各类细胞的关系。

可以发现,MDS 利用的仅是数据间距离关系。在 MDS 的视角中,样本的空间坐标仅是样本间距离的一种表示方式。坐标只要是能准确反映出样本间距离关系即可,维数不再重要。所以,若已知高维数据点间的距离,且数据点的绝对位置对我们的研究目的而言意义不大时,可以使用 MDS 方法对数据进行降维。

MDS 利用的是样本间距离关系,那其他降维方法利用的是什么呢?其实无论是 pca、tsne 还是 MDS,所有的降维方法本质上都是利用样本间的距离关系,但不同方法利用距离关系的方式有一定区别。其中 pca 只是在原空间中重新找一组正交坐标系来描述样本,不改变数据点的绝对位置。但如果只取 pca 的前 N 维的坐标,那么其含义则变成了样本从高维空间向这 N 维空间的投影,其样本间的距离自然也发生了改变。而 tsne、MDS 等方法都是利用低维空间中的距离矩阵来拟合高维空间中的距离矩阵,或多或少都会改变一些数据的相对位置。具体改变的程度与各算法的映射方法相关。总的来说,pca 是最 “原汁原味” 的降维方法。所以 如果 pca 降维后的空间维数低,不需要再度精简,没有可视化需求时,pca 无疑是一个非常优秀的方法。也可以说 “能用 pca 解决的降维问题,尽量都用 pca 来解决,别整那些花里胡哨的”。如处理实验室数据库中 368 个玉米的高光谱数据,降维的目的是去冗,并且 pca 降维后空间维数低。但 如果固定了维数(如可视化需要将数据降到 2、3 维),并且此维度下 pca 所能解释的方差较小,则不再适合使用 pca 降维方法,如本论文中的单细胞转录组数据。

相比 pca 直接利用坐标进行降维,MDS 这类利用距离矩阵 “拟合” 的降维方法还有一个优势:可以使用不同的距离度量方式得到具有不同侧重点的距离矩阵。如本研究使用扩散映射 (diffusion map) 度量样本间的距离,使其即包含局部信息又包含全局信息;如 tsne (t-distributed stochastic neighbor embedding,t 分布随机邻近嵌入) 使用 t 分布度量样本间的距离,使其更关注局部信息,较为忽略全局信息。这一特性也使得 距离矩阵类降维方法可以处理许多非线性(流形)结构的数据。所以 距离矩阵相比 pca 有更加广阔的使用范围

5. DEMaP 度量标准

DEMaP 为样本间在高维空间中的测地距离与在低维空间中的欧式距离的相关性。作者通过 DEMaP 度量标准,同时衡量降维方法保持局部、全局结构的能力以及去噪的能力。这里我对 DEMaP 度量标准的有效性保持怀疑,如果测地距离是黄金标准,那为什么 MDS 的距离矩阵不用测地距离计算得到的而用扩散映射方法计算得到的?

启发 & 应用

  1. 扩散距离是否可以替代我们实验室中常用的一些度量距离的方法,如欧式距离、pearson相关系数等,用于自学习等。
  2. 扩散距离矩阵是否可以替代kinship从而更好的度量个体之间的相关性。
  3. 利用 PHATE 是否可以判断植株的长势分化,是否有按照我们期望的方向发展。
  4. 类比细胞的发育,PHATE 可以尝试用于解释物种的遗传进化路线。
  5. 通过 batch 研究基因在时空上的配合。
降维及可视化
驽马十驾,功在不舍;锲而不舍,金石可镂。
07-28 1336
降维及可视化1.PCA2.T-SNE3.Unimap4.phate5.seaborn参考文献 1.PCA 2.T-SNE 3.Unimap 4.phate 5.seaborn 参考文献 [1] seaborn之kdeplot
10X单细胞降维分析之PHATE
最新发布
weixin_53637133的博客
04-23 900
10X单细胞降维分析之PHATE
PHATE:PHATE(基于亲和力的过渡嵌入的热扩散潜能)是用于可视化高维数据的工具
05-04
PHATE-可视化生物数据探索的过渡和结构 快速开始 如果您想开始使用PHATE,请查看以下教程。 介绍 PHATE(基于亲和力的轨迹嵌入的热扩散潜能)是一种可视化高维数据的工具。 PHATE使用新颖的概念框架来学习和可视化流形,以保留局部和全局距离。 要了解如何将PHATE应用于人类胚胎干细胞的面部图像和单细胞数据集等数据集,请查看我们在《自然生物技术》上的出版物。 PHATE已在 , MATLAB和R中实现。 目录 系统要求 Python 用pip安装 从源安装 快速开始 教程和参考 的MATLAB 安装 教程和参考 [R 从CRAN和PyPi安装 使用devtools安装并进行网状安装 从源安装 快速开始 教程和参考 帮助 系统要求 Windows(> = 7),Mac OS X(> = 10.8)或Linux Python> = 3.5或MATLAB (>
Nat. Biotechnol | PHATE:高维生物数据的可视化方法
DrugAI
12-20 1907
作者 | 涂心琪单位 | 湖南大学研究方向 | 高维数据可视化高维生物数据的可视化能帮助研究者以直观的方式了解数据。今天介绍2019年12月发表在Nature Biotechnology...
m-phate:用于张量嵌入的多层PHATE
05-18
邻苯二甲酸酯 多层PHATE(M-PHATE)是一种降维算法,用于可视化随时间变化的数据。 要了解有关M-PHATE的更多信息,您可以阅读arXiv上的预印本,在该预印本中,我们将其应用于训练过程中神经网络的演化。 上面我们展示了一个M-PHATE在300多个训练阶段应用于三层MLP的演示,该训练由一个时期(左),隐藏层(中)和最能激活每个隐藏单元的数字标签(右)进行了着色。 在下面,您会看到在3D嵌入的训练中应用了辍学的同一网络,该网络也由大多数活动单位着色。 目录 参数调整 人物复制 帮助 这个怎么运作 多层PHATE(M-PHATE)将新颖的多层内核构造与PHATE可视化结合在一起。 我们的内核捕获了不断变化的图结构的动态特性,当可视化时,它给出了有关系统演变的独特直觉; 在我们的预印本中,我们显示了在训练和再训练过程中将其应用于神经网络的过程。 我们将M-PHATE与其他降维技
phateR:R中实现PHATE降维的方法
05-22
磷酸盐 该R包提供了。 有关PHATE可视化方法的完整概述,请参见《。 有关我们的Python和Matlab实现,请参见 。 目录 快速开始 教程 问题 常问问题 帮助 安装 为了在R中使用PHATE,您还必须安装Python软件包。 如果未安装python或pip ,则需要安装它们。 我们建议Miniconda3一起安装Python和pip ,否则,您可以从https://pip.pypa.io/en/stable/installing/安装pip 。 从CRAN和PyPi安装 首先通过在终端上运行以下代码在Python中安装phate : pip install --user phate 然后通过在R中运行以下代码从CRAN安装phateR : install.packages( " phateR " ) 使用devtools安装并进行reticulate安装 可以
ppkg:Phate的软件包管理器,用Rust编写的二进制软件包管理器(计划为混合)
02-07
**ppkg:Phate的软件包管理器** `ppkg` 是一个专为Phate操作系统设计的软件包管理器,其核心是用编程语言Rust实现的。Rust以其内存安全性和高性能而闻名,这使得`ppkg`在处理软件包安装、升级和卸载时能提供高效且...
PackagesNotFoundError: The following packages are not available from current channels: - phate
07-15
PackagesNotFoundError 错误提示意味着当前的软件源中找不到所请求的软件包(如 `phate`)。 要解决这个问题,您可以尝试以下方法之一: 1. 检查拼写错误:确保您输入的软件包名称正确,没有拼写错误。 2. 更新...
PyPI 官网下载 | phate-0.4.5-py3-none-any.whl
02-07
在Python环境中,PhATE可以通过Python Package Index (PyPI) 官方平台获取,如标题所示的"PyPI 官网下载 | phate-0.4.5-py3-none-any.whl"文件。 PhATE的核心功能是生成热图表示,这种表示方法能够捕捉到细胞间的...
探索生物数据的新维度:PHATE
gitblog_00083的博客
04-19 292
探索生物数据的新维度:PHATE 项目地址:https://gitcode.com/KrishnaswamyLab/PHATE PHATE 是一个强大的可视化工具,全称是“潜在的高维异构嵌入”,旨在帮助研究人员在高复杂度的单细胞或多细胞数据中揭示出生物系统的结构和动态。作为一个基于JavaScript的开源库,PHATE 可以轻松地集成到任何Web应用或数据分析工作流中。 技术分析 PHATE 的...
《统计学习方法》笔记
Algernon98的博客
09-05 924
李航版《统计学习方法》学习笔记,随科研进度需要填充内容
Unsupervised Image Matching based on Manifold Alignment(笔记)
thystar的专栏
11-07 1395
基于流形对齐的无监督图像配准
示例代码-对称正定流形(Symmetric Positive Definite Manifold ,简称SPD流形)上距离计算常用的四种度量
Kai-Xuan
05-14 4566
示例代码:https://github.com/Kai-Xuan/MyNote/tree/master/ML/SPD-Metrics 对于任意两个SPD 矩阵 X,Y 放射不变黎曼度量 (Affine Invariant Riemannian Metric): dA2(X,Y)=∣∣log⁡(X−1/2Ylog⁡(X−1/2)∣∣F2d_A^2(X,Y)=|| \log(X^{-1/2}Y \log(X^{-1/2}) ||_F^2dA2​(X,Y)=∣∣log(X−1/2Ylog(X−1/2)∣∣F2​
流形学习(转)
aasys848082的博客
06-27 506
流形(manifold)的概念最早是在1854年由 Riemann 提出的(德文Mannigfaltigkeit),现代使用的流形定义则是由 Hermann Weyl 在1913年给出的。江泽涵先生对这个名词的翻译出自文天祥《正气歌》“天地有正气,杂然赋流形”,日本人则将之译为“多样体”,二者孰雅孰鄙,高下立判。 流形(Manifold),一般可以认为是局部具有欧氏空间性质的空间。而实...
微分流形与黎曼几何学习笔记(转自http://blog.sciencenet.cn/home.php?mod=space&uid=81613&do=blog&id=333317)
西西弗的石头
08-10 8096
由于种种原因要恶补一下微分流形和黎曼几何,吸取一下“前辈”们的经验,也希望大家能提供一些更好的经验!   1  自几何佳缘 在这方面我是很有感受的。我整理了一些心得笔记,打算以后给学生上课的时候,把这些内功心法传授给他们。  这里先随便讲两句。 如果楼主想聊聊的话,可
Task00 Pandas安装指令和初识
jjwan123的博客
12-13 117
Pandas 是python的一个数据分析包,该工具是为了解决数据分析任务而创建的! “七星联盟的小伙伴”,直接使用“pip install pandas”直接安装第三方库就可以了
多目标优化算法平台PlatEMO的基本使用方法
热门推荐
认知无线电
09-04 1万+
简介和使用方法PlatEMO是基于Matlab的一款拥有可视化界面的多目标优化算法平台。下载文件,下载链接:https://github.com/BIMK/PlatEMO文件下载完成后,...
python用tsne降维图像_python代码实现TSNE降维数据可视化教程
weixin_42385930的博客
02-10 3361
TSNE降维降维就是用2维或3维表示多维数据(彼此具有相关性的多个特征数据)的技术,利用降维算法,可以显式地表现数据。(t-SNE)t分布随机邻域嵌入 是一种用于探索高维数据的非线性降维算法。它将多维数据映射到适合于人类观察的两个或多个维度。python代码km.py#k_mean算法import pandas as pdimport csvimport pandas as pdimport nu...
杭电4 License Plate Recognition 降维映射
srh20的博客
07-29 256
https://acm.hdu.edu.cn/showproblem.php?pid=6993 License Plate Recognition Time Limit: 4000/2000 MS (Java/Others)Memory Limit: 262144/262144 K (Java/Others) Total Submission(s): 141Accepted Submission(s): 74 Problem Description Little Rabbit i...
写文章

热门文章

  • 利用 DIAMOND、MCScanX、TBtools 分析物种基因组间的共线性区段与基因复制事件 24022
  • 孟德尔随机化法(Mendelian Randomization,MR) 18775
  • MEME 使用简介 16325
  • 使用 Blastp 和 Hmmer 筛选出包含特定结构域的蛋白 16061
  • MUSCLE、IQtree 软件及使用简介 14318

分类专栏

  • 3D 基因组 1篇
  • 玉米基因组 7篇
  • 进化遗传学 11篇
  • 系统发育树 7篇
  • 基因组学 1篇
  • RNA甲基化 2篇
  • 转录组学 3篇
  • 降维 1篇
  • 可视化 1篇
  • 解释模型 2篇
  • 现代农业 1篇
  • 高光谱 1篇
  • Motif 2篇
  • github 1篇
  • LightGBM 1篇
  • 线性代数 2篇
  • 多组学分析 1篇
  • 自训练

最新评论

  • IQtree:使用 SNP 数据(vcf file)构建玉米群体的 无根 系统发育树

    2401_82438120: log for logging怎么解决啊?

  • IQtree:使用 SNP 数据构建 有根 系统发育树及踩坑

    2401_82438120: log for logging怎么解决啊?

  • MUSCLE、IQtree 软件及使用简介

    m0_75044134: 你好,我和你的情况一样,请问你的解决了吗

  • 利用 OrthoFinder、IQtree、Notung、iTOL 绘制基因树

    2301_78533177: 请问为什么用Notung基因树显示复制与丢失结果但在物种树上显示不出来呀

  • 利用 OrthoFinder、IQtree、Notung、iTOL 绘制基因树

    浓香鸭腿面: 当物种较多时,多生根点是无法避免的,一般选择距离外群最近的生根点

最新文章

  • DipC 构建基因组 3D 结构(学习笔记)
  • IQtree:使用 SNP 数据构建 有根 系统发育树及踩坑
  • IQtree:使用 SNP 数据(vcf file)构建玉米群体的 无根 系统发育树
2023年1篇
2022年11篇
2021年20篇
2020年13篇

目录

目录

评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43元 前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值

玻璃钢生产厂家玻璃钢雕塑的设计特别推荐南川玻璃钢雕塑公司玻璃钢卡通雕塑设计报价福建个性化玻璃钢雕塑定制四川特色商场美陈现价上海商场主题创意商业美陈报价萍乡景观玻璃钢雕塑市场亳州定做玻璃钢雕塑厂选哪家玻璃钢天壶流水景观雕塑江西省玻璃钢雕塑找哪家现代玻璃钢雕塑报价表新密价值观玻璃钢人物雕塑厂家景观玻璃钢雕塑单价山西个性化玻璃钢雕塑订做价格湛江玻璃钢造型雕塑邵阳人物玻璃钢雕塑价格义乌玻璃钢雕塑制造厂家最好的云南玻璃钢雕塑保定玻璃钢雕塑报价玻璃钢雕塑安装工艺玻璃钢雕塑制品图片泸州市玻璃钢雕塑定制浙江永康玻璃钢景观雕塑供应商云南景区玻璃钢雕塑制作惠州玻璃钢雕塑接单中原财富网创新玻璃钢雕塑厂广东大型商场美陈批发温州玻璃钢景观雕塑制造商抽象玻璃钢花盆儿赤峰玻璃钢卡通雕塑香港通过《维护国家安全条例》两大学生合买彩票中奖一人不认账让美丽中国“从细节出发”19岁小伙救下5人后溺亡 多方发声单亲妈妈陷入热恋 14岁儿子报警汪小菲曝离婚始末遭遇山火的松茸之乡雅江山火三名扑火人员牺牲系谣言何赛飞追着代拍打萧美琴窜访捷克 外交部回应卫健委通报少年有偿捐血浆16次猝死手机成瘾是影响睡眠质量重要因素高校汽车撞人致3死16伤 司机系学生315晚会后胖东来又人满为患了小米汽车超级工厂正式揭幕中国拥有亿元资产的家庭达13.3万户周杰伦一审败诉网易男孩8年未见母亲被告知被遗忘许家印被限制高消费饲养员用铁锨驱打大熊猫被辞退男子被猫抓伤后确诊“猫抓病”特朗普无法缴纳4.54亿美元罚金倪萍分享减重40斤方法联合利华开始重组张家界的山上“长”满了韩国人?张立群任西安交通大学校长杨倩无缘巴黎奥运“重生之我在北大当嫡校长”黑马情侣提车了专访95后高颜值猪保姆考生莫言也上北大硕士复试名单了网友洛杉矶偶遇贾玲专家建议不必谈骨泥色变沉迷短剧的人就像掉进了杀猪盘奥巴马现身唐宁街 黑色着装引猜测七年后宇文玥被薅头发捞上岸事业单位女子向同事水杯投不明物质凯特王妃现身!外出购物视频曝光河南驻马店通报西平中学跳楼事件王树国卸任西安交大校长 师生送别恒大被罚41.75亿到底怎么缴男子被流浪猫绊倒 投喂者赔24万房客欠租失踪 房东直发愁西双版纳热带植物园回应蜉蝣大爆发钱人豪晒法院裁定实锤抄袭外国人感慨凌晨的中国很安全胖东来员工每周单休无小长假白宫:哈马斯三号人物被杀测试车高速逃费 小米:已补缴老人退休金被冒领16年 金额超20万

玻璃钢生产厂家 XML地图 TXT地图 虚拟主机 SEO 网站制作 网站优化