Python中基于基因表达矩阵的层次聚类

Question

我想知道如何在Python中进行层次聚类（这里是针对基因表达数据），并且能够同时显示基因表达值的矩阵和树状图。我的意思是像这里的例子：

http://www.mathworks.cn/access/helpdesk/help/toolbox/bioinfo/ug/a1060813239b1.html

在第6个要点之后展示的（图1），树状图被绘制在基因表达矩阵的左侧，行的顺序已经重新排列，以反映聚类的结果。

我该如何在Python中使用numpy/scipy或其他工具来实现这个呢？另外，使用欧几里得距离作为度量标准，处理大约11,000个基因的矩阵在计算上是否可行？

编辑：很多人建议使用聚类包，但我仍然不确定如何在Python中绘制我上面链接的那种图像。我该如何使用Matplotlib等工具将树状图叠加在热图矩阵旁边呢？

谢谢。

Answer 1

有几个人在使用scipy和matplotlib创建层次聚类和热图可视化的原型模块方面取得了一些不错的进展。

如何获取与scipy创建的树状图中的颜色簇对应的平面聚类

我一直在调整这段代码，想要做一个完整的层次聚类模块，这样我就可以把它集成到我的转录组分析软件包中。我对最终的产品很满意，它可以使用各种聚类指标和方法生成热图，并且有颜色渐变。这里展示了代码和一个示例输出：

http://altanalyze.blogspot.com/2012/06/hierarchical-clustering-heatmaps-in.html

Answer 2

你可以使用scipy的cluster.hierarchy模块来实现这个功能。其实这些命令非常相似。不过，你需要在pdist中用correlation来代替corr，而不是用cluster，scipy的聚类模块里对应的函数是fcluster。另外，关于树状图，scipy中的函数是dendrogram，而在Matlab中是clustergram。

你完全可以使用欧几里得距离（我认为这也是pdist的默认设置）。我觉得处理11,000个基因是可行的，因为这需要计算的距离是11000*(11000-1)/2 = 60494500（从11000中选2）的组合。这是个很大的数字，但我认为是可以做到的。

Answer 3

许多聚类方法，包括 scipy.cluster，首先会对所有成对的距离进行排序，像你这种情况大约有6000万条，算是比较小的数量。
你觉得下面这段代码运行需要多长时间呢？

import scipy.cluster.hierarchy as hier
import pylab as pl

def fcluster( pts, ncluster, method="average", criterion="maxclust" ):
    """ -> (pts, Y pdist, Z linkage, T fcluster, clusterlists)   
        ncluster = n1 + n2 + ... (including n1 singletons)
        av cluster size = len(pts) / ncluster
    """
    pts = np.asarray(pts)
    Y = scipy.spatial.distance.pdist( pts )  # ~ N^2 / 2
    Z = hier.linkage( Y, method )  # N-1                         
    T = hier.fcluster( Z, ncluster, criterion=criterion )
        # clusters = clusterlists(T)
    return (pts, Y, Z, T)

hier.dendrogram( Z )

关于如何重新排列矩阵并漂亮地绘图的问题，可以在这里找到，那个帖子是在三月份提问的，里面有部分答案。