巨大数据集的热力图

Question

我有一个制表符分隔的文件，其中包含区域和这些区域中发现的相应生物实体（我已经检查了67个，因此您可以说每个区域都被检查了这67个实体的存在或缺失以及它们的频率）。

我将所有这些数据以表格格式呈现。

以下是示例数据：

Region  ATF3    BCL3    BCLAF1  BDP1    BRF1    BRF2    Brg1    CCNT2   CEBPB   CHD2    CTCF    CTCFL   E2F6    ELF1
chr1:109102470:109102970    0   0   1   0   0   0   0   1   0   0   4   1   4   1
chr1:110526886:110527386    0   0   0   0   0   0   0   1   1   0   4   1   0   1
chr1:115300671:115301171    0   0   1   0   0   0   0   0   1   1   4   1   1   1
chr1:115323308:115323808    0   0   0   0   0   0   0   1   0   0   2   1   1   0
chr1:11795641:11796141  1   0   0   0   0   0   0   1   2   0   0   0   1   0
chr1:118148103:118148603    0   0   0   0   0   0   0   1   0   0   0   0   0   1
chr1:150521397:150521897    0   0   0   0   0   0   0   2   2   0   6   2   4   0
chr1:150601609:150602109    0   0   0   0   0   0   0   0   3   2   0   0   1   0
chr1:150602098:150602598    0   0   0   0   0   0   0   0   1   1   0   0   0   0
chr1:151119140:151119640    0   0   0   0   0   0   0   1   0   0   0   0   1   0
chr1:151128604:151129104    0   0   0   0   0   0   0   0   0   0   3   0   0   0
chr1:153517729:153518229    0   0   0   0   0   0   0   0   0   0   0   0   0   0
chr1:153962738:153963238    0   0   0   0   0   0   0   1   1   0   0   0   0   1
chr1:154155682:154156182    0   0   0   0   0   0   0   1   0   0   0   0   1   1
chr1:154155725:154156225    0   0   0   0   0   0   0   1   0   0   0   0   1   1
chr1:154192154:154192654    0   0   0   0   0   0   0   0   0   0   0   0   0   0
chr1:154192824:154193324    1   0   0   0   0   0   0   1   0   1   0   0   1   1
chr1:154192943:154193443    1   0   0   0   0   0   0   1   0   2   0   0   1   1
chr1:154193273:154193773    1   0   0   0   0   0   0   1   0   2   0   0   2   1
chr1:154193313:154193813    0   0   0   0   0   0   0   1   0   2   0   0   2   1
chr1:155904188:155904688    0   0   0   0   0   0   0   1   0   0   0   0   1   1
chr1:155947966:155948466    0   0   0   0   0   0   0   1   0   0   3   0   0   1
chr1:155948336:155948836    0   0   0   0   0   0   0   1   0   0   5   1   0   1
chr1:156023516:156024016    0   0   0   0   0   0   0   1   0   1   4   1   1   1
chr1:156024016:156024516    0   1   1   0   0   0   0   0   0   2   0   0   1   1
chr1:156163229:156163729    0   0   0   0   0   0   0   0   0   0   2   0   0   1
chr1:160990902:160991402    0   0   0   0   0   0   0   0   0   1   0   0   1   2
chr1:160991133:160991633    0   0   0   0   0   0   0   0   0   1   0   0   1   2
chr1:161474704:161475204    0   0   0   0   0   0   0   0   0   0   0   0   0   0
chr1:161509530:161510030    0   0   1   1   1   0   0   0   1   0   1   0   0   1
chr1:161590964:161591464    0   0   0   1   1   0   0   0   0   0   0   0   0   0
chr1:169075446:169075946    0   0   0   0   0   0   0   2   0   0   4   0   3   0
chr1:17053279:17053779  0   0   0   1   0   0   0   0   0   1   0   0   0   0
chr1:1709909:1710409    0   0   0   0   0   0   0   2   0   1   0   0   3   1
chr1:1710297:1710797    0   0   0   0   0   0   0   0   0   1   6   0   1   1

现在我该如何将这个数据用热力图表示，颜色从浅红到深红（根据频率），没有数据的地方用白色表示？

还有其他更好的方式来呈现这种类型的数据吗？

- Angelo

---

除了你的数据，你能展示一下你尝试过的内容吗？ - Inbar Rose

嗯，这是一个实用编程的问答网站。我觉得在信息图表方面可能没有专家。 - Denis

你所说的“huge”是指什么？是指很多行还是很多列？如果其中一个数字接近你要创建图表的像素数，那将不是一个好的方法。 - Thorsten Kranz

@ThorstenKranz 67列和1100行 - Angelo

@Angelo，你应该将图像保存为大图像（使用plt.savefig("output.png")或类似方法），这样你的1100行才会有意义。 - Thorsten Kranz

Answer 1

使用 Matplotlib

import pylab as plt
import numpy as np

data = np.loadtxt("14318737.txt", skiprows=1, converters={0:lambda x: 0})
plot_data = np.ma.masked_equal(data[:,1:], 0)

plt.imshow(plot_data, cmap=plt.cm.get_cmap("Reds"), interpolation="nearest")
plt.colorbar()

plt.show()

我忽略了第一行和第一列（如果您需要它们作为标签，我们需要更改这个）。对于剩余的数据，所有零值都被屏蔽（因此在图中显示为白色），然后将这些数据绘制为彩色编码的图形。 imshow有许多其他参数可控制结果，例如起点（下/上）、方面（自动/相等/某些比率）。

你提到了地区 - 你是指地理区域吗？那么你可能想看看Basemap Toolkit for Matplotlib来创建彩色编码地图。 编辑 新要求，新示例。

import pylab as plt
import numpy as np

fn = "14318737.txt"
with open(fn, "r") as f:
    labels = f.readline().rstrip("\n").split()[1:]
data = np.loadtxt(fn, skiprows=1, converters={0:lambda x: 0})
plot_data = np.ma.masked_equal(data[:,1:], 0)

plt.subplots_adjust(left=0.1, bottom=0.15, right=0.99, top=0.95)
plt.imshow(plot_data, cmap=plt.cm.get_cmap("Reds"), interpolation="nearest", aspect = "auto")
plt.xticks(range(len(labels)), labels, rotation=90, va="top", ha="center")
plt.colorbar()

plt.show()

现在我首先从第一行读取标签。我添加了关键字参数aspect到imshow调用中。我为每个因素创建标签。

此外，我使用subplots_adjust调整图的位置。您可以根据需要调整这些参数。

现在的结果是： 如果您想要y轴的其他刻度，请使用plt.yticks，就像我例子中的xticks一样。

- Thorsten Kranz

---

这个我怀疑过，所以你不需要它。 - Thorsten Kranz

Thorsten Kranz：结果也可以进行聚类吗？ - Angelo

你能解释一下你所说的“cluster as well the results”是什么意思吗？你是指在这个二维图像中找到聚类吗？这是可能的，但我担心你可能指的是其他的东西。 - Thorsten Kranz

是的，我是指“在这个二维图像中找到聚类”，并且从零开始标尺而不是1。你能否请修改脚本并将这些内容整合进去，因为我没有成功调整它。谢谢。 - Angelo

问题2很简单 - 只需在imshow中使用一个额外的关键字参数，即vmin=0。问题1稍微复杂一些。你如何定义一个聚类？所有连接的值都高于某个阈值吗？那个阈值是多少？是0吗？我会发另一个答案，附上一些代码示例。 - Thorsten Kranz

Answer 2

由于我的另一个答案引起了评论，OP有关于搜索2D聚类的另一个问题。这里是一些解答。

从我的库eegpy中获取，我使用一个名为find_clusters的方法。它在2D数组上执行遍历，找到所有高于/低于给定阈值的聚类。

以下是我的代码：

import pylab as plt
import numpy as np
from Queue import Queue


def find_clusters(ar,thres,cmp_type="greater"):
    """For a given 2d-array (test statistic), find all clusters which
are above/below a certain threshold.
"""
    if not cmp_type in ["lower","greater","abs_greater"]:
        raise ValueError("cmp_type must be in [\"lower\",\"greater\",\"abs_greater\"]")
    clusters = []
    if cmp_type=="lower":
        ar_in = (ar<thres).astype(np.bool)
    elif cmp_type=="greater":
        ar_in = (ar>thres).astype(np.bool)
    else: #cmp_type=="abs_greater":
        ar_in = (abs(ar)>thres).astype(np.bool)

    already_visited = np.zeros(ar_in.shape,np.bool)
    for i_s in range(ar_in.shape[0]): #i_s wie i_sample
        for i_f in range(ar_in.shape[1]):
            if not already_visited[i_s,i_f]:
                if ar_in[i_s,i_f]:
                    #print "Anzahl cluster:", len(clusters)
                    mask = np.zeros(ar_in.shape,np.bool)
                    check_queue = Queue()
                    check_queue.put((i_s,i_f))
                    while not check_queue.empty():
                        pos_x,pos_y = check_queue.get()
                        if not already_visited[pos_x,pos_y]:
                            #print pos_x,pos_y
                            already_visited[pos_x,pos_y] = True
                            if ar_in[pos_x,pos_y]:
                                mask[pos_x,pos_y] = True
                                for coords in [(pos_x-1,pos_y),(pos_x+1,pos_y),(pos_x,pos_y-1),(pos_x,pos_y+1)]: #Direct Neighbors
                                    if 0<=coords[0]<ar_in.shape[0] and 0<=coords[1]<ar_in.shape[1]:
                                        check_queue.put(coords)
                    clusters.append(mask)
    return clusters

fn = "14318737.txt"
with open(fn, "r") as f:
    labels = f.readline().rstrip("\n").split()[1:]
data = np.loadtxt(fn, skiprows=1, converters={0:lambda x: 0})

clusters = find_clusters(data, 0, "greater")

plot_data = np.ma.masked_equal(data[:,1:], 0)

plt.subplots_adjust(left=0.1, bottom=0.15, right=0.99, top=0.95)
plt.imshow(plot_data, cmap=plt.cm.get_cmap("Reds"), interpolation="nearest", aspect = "auto", 
           vmin=0, extent=[0.5,plot_data.shape[1]+0.5, plot_data.shape[0] - 0.5, -0.5])
plt.colorbar()

for cl in clusters:
    plt.contour(cl.astype(np.int),[0.5], colors="k", lw=2)
plt.xticks(np.arange(1, len(labels)+2), labels, rotation=90, va="top", ha="center")


plt.show()

这句话的意思是：“它提供了一个形式的图像：”，其中保留了HTML标记。

clusters是一个布尔值2D数组（True / False）的列表。每个数组表示一个聚类，其中每个布尔值指示特定“点”是否属于此聚类。您可以在任何进一步的分析中使用它。

编辑

现在有一些关于聚类的有趣信息

import pylab as plt
import numpy as np
from Queue import Queue


def find_clusters(ar,thres,cmp_type="greater"):
    """For a given 2d-array (test statistic), find all clusters which
are above/below a certain threshold.
"""
    if not cmp_type in ["lower","greater","abs_greater"]:
        raise ValueError("cmp_type must be in [\"lower\",\"greater\",\"abs_greater\"]")
    clusters = []
    if cmp_type=="lower":
        ar_in = (ar<thres).astype(np.bool)
    elif cmp_type=="greater":
        ar_in = (ar>thres).astype(np.bool)
    else: #cmp_type=="abs_greater":
        ar_in = (abs(ar)>thres).astype(np.bool)

    already_visited = np.zeros(ar_in.shape,np.bool)
    for i_s in range(ar_in.shape[0]): #i_s wie i_sample
        for i_f in range(ar_in.shape[1]):
            if not already_visited[i_s,i_f]:
                if ar_in[i_s,i_f]:
                    #print "Anzahl cluster:", len(clusters)
                    mask = np.zeros(ar_in.shape,np.bool)
                    check_queue = Queue()
                    check_queue.put((i_s,i_f))
                    while not check_queue.empty():
                        pos_x,pos_y = check_queue.get()
                        if not already_visited[pos_x,pos_y]:
                            #print pos_x,pos_y
                            already_visited[pos_x,pos_y] = True
                            if ar_in[pos_x,pos_y]:
                                mask[pos_x,pos_y] = True
                                for coords in [(pos_x-1,pos_y),(pos_x+1,pos_y),(pos_x,pos_y-1),(pos_x,pos_y+1)]: #Direct Neighbors
                                    if 0<=coords[0]<ar_in.shape[0] and 0<=coords[1]<ar_in.shape[1]:
                                        check_queue.put(coords)
                    clusters.append(mask)
    return clusters

fn = "14318737.txt"
data = []
with open(fn, "r") as f:
    labels = f.readline().rstrip("\n").split()[1:]
    for line in f:
        data.append([int(v) for v in line.split()[1:]])
data = np.array(data) #np.loadtxt(fn, skiprows=1, usecols=range(1,15))#converters={0:lambda x: 0})

clusters = find_clusters(data, 0, "greater")
large_clusters = filter(lambda cl: cl.sum()>5, clusters) #Only take clusters with five or more items
large_clusters = sorted(large_clusters, key=lambda cl: -cl.sum())

plot_data = np.ma.masked_equal(data[:,:], 0)

plt.subplots_adjust(left=0.1, bottom=0.15, right=0.99, top=0.95)
plt.imshow(plot_data, cmap=plt.cm.get_cmap("Reds"), interpolation="nearest", aspect = "auto", 
           vmin=0, extent=[-0.5,plot_data.shape[1]-0.5, plot_data.shape[0] - 0.5, -0.5])
plt.colorbar()

for cl in large_clusters:
    plt.contour(cl.astype(np.int),[.5], colors="k", lw=2)
plt.xticks(np.arange(0, len(labels)+1), labels, rotation=90, va="top", ha="center")

print "Summary of all large clusters:\n"
print "#\tSize\tIn regions"
for i, cl in enumerate(large_clusters):
    print "%i\t%i\t" % (i, cl.sum()),
    regions_in_cluster = np.where(np.any(cl, axis=0))[0]
    min_region = labels[min(regions_in_cluster)]
    max_region = labels[max(regions_in_cluster)]
    print "%s to %s" % (min_region, max_region)

plt.xlim(-0.5,plot_data.shape[1]-0.5)

plt.show()

我过滤掉包含五个以上点的所有聚类，只绘制这些聚类。您也可以使用每个聚类内数据的总和。然后按大小降序排序这些大型聚类。

最后，我打印所有大型聚类的摘要，包括它们跨越的所有聚类的名称。

- Thorsten Kranz

---

感谢提供更新的代码。由于我的数据有2700个区域，是否可以基于行进行聚类呢？我得到了一个非常糟糕的图像。 - Angelo

"基于行进行聚类"是什么意思？我的算法同时在行和列中搜索。你是指只在行中搜索吗？ - Thorsten Kranz

好的，我明白了。我能否将这些聚类以某种变量的形式获取，并按照从最大聚类到最小聚类的顺序显示它们的名称（例如聚类1 --> E2F6和ELF1）？ - Angelo

是的，你可以用这些集群做任何你想做的事情。我建议过滤掉所有小的集群（统计波动）等等...我会给你一个例子。 - Thorsten Kranz