大数据集(>19000 行)上的 Python All 与 All 比率比较不断崩溃

Question

我一直在做一些关于基因表达和药物相互作用的研究，最近被要求对整个人类基因组(> 19000 行)进行分析。本质上，我的程序所做的就是获取所有表达数据的比率，并将其与一些药物信息(文件的最后 9 行)相关联。

发生的情况是，每次我运行代码时，我的计算机都会崩溃。有没有任何明显的方法可以使其工作而不必在全新的机器上运行它？我正在远程工作，因此目前不太可能。我正在使用相对较新的 macbook pro 在终端上运行 python。我的代码如下，如果有任何问题需要清理，请告诉我。感谢您的帮助!

编辑 1:一位同事提到使用 Amazon Web Services 来运行此程序。这听起来是一个可行的选择吗？从我目前所看到的来看，仅运行这个 Python 脚本的设置似乎有点复杂。

编辑 2:以下是一些输入文件的示例:行继续向下表示 1970 个基因，列继续表示 57 个不同的细胞系。

              WT     Mut    Mut     Mut
              WT     Mut    Mut     Mut
Cell lines  BR:MCF7 BR:MDA_MB_231   BR:HS578T
            BR:MCF7 BR:MDA_MB_231   BR:HS578T
5-HT3C2     -0.27   0.99    0.7 -0.42
A1BG-AS1    1.36    -0.15   0.87    1.7
A1CF        -0.14   -0.18   0.15    -0.1
A2LD1       0.62    -0.59   -0.29   2.45
A2M        -0.38    -0.4    -0.24   -0.39
A2ML1       -0.11   -0.13   -0.13   -0.12
A2MP1       0.31    0.65    1.2 0.03
A4GALT      1.99    0.41    -0.75   0.19
A4GNT       0.28    0.08    1.08    0.74
AAA1        -0.27   -0.25   -0.19   -0.16
AAAS        1.16    -1.46   -2.06   -1
AACS        0.73    0.11    -1.11   -2.08

代码

import csv
import numpy
from scipy.stats.stats import pearsonr

Reader = csv.reader(open('/Users/_57_genes.csv', 'rU'))

fout = "/Users/ratio_correlations_whole_genome_6232014.csv"

fileout = open(fout, 'w') #open file

#Create lists of lists including averages
Label = []
Resistance = []


for row in Reader:
    if len(row) == 0 or row[0] == '':
        continue
    else:
        Resistance.append(map(float, row[2:]))
        Label.append(row[0])



Ratios = []
Name = []

for index, i in enumerate(Resistance[0:-9]):
    gene = []
    Name1 = []
    for index2, j in enumerate(Resistance[0:-9]):
        r = []
        if j == i:
             continue
        for k in range(len(i)):
            if j[k] == 0:
                fraction = 0
            else:
                fraction = (i[k]/j[k])
            r.append(fraction)
        Name1.append('%s VS %s' %(Label[index], Label[index2]) )
        gene.append(r)
    Ratios.append(gene)
    Name.append(Name1)



for index, i in enumerate(Ratios):
    GeneName = Name[index]
    for index2, k in enumerate(i):
        Comparitor = GeneName[index2]
        #print k
        fileout.write('%s, ' % (Comparitor))
        for avg in (Resistance[-9:]):
            correlate = pearsonr(k, avg)
            fileout.write('%0.6f, %0.6f, ' % correlate)
        fileout.write('\n')

fileout.close()

Answer 1

(我在这里假设是 64 位架构；如此大的数据的 32 位代码会在内存耗尽之前耗尽地址空间......我已经使用 pympler 测量了对象的大小code> 可从 pypi 获取模块。)

您正在构建两个包含 19000*19888=360981000 个项目的二维数组。名称数组本身包含那么多字符串。 Python 中的每个字符串对象至少需要 40 个字节，再加上您需要为对象指针添加至少 8 个字节……这意味着这个数组至少需要 13GB。

更有趣的是另一个数组:另一个数组每个项目至少需要 112 个字节(至少一个 double 的列表需要 104 个字节，指针需要 8 个字节)，这又是 37GB。

您的计算机有 50GB RAM 吗？如果没有，您需要认真重新考虑您的代码。您真的需要同时存储所有值吗？也许您可以使用不同的算法，而不是从表中获取值，而是仅根据需要(当实际需要时)计算必要的值？

也许您可以使用 numpy 数组来代替 Python 列表，它在内存中的效率更高( double 组的每个项目 8 个字节，而 double 组的每个项目 32 个字节)。

假设您可以从他们的报价中购买最昂贵的机器，那么获得一台 AWS 机器似乎是一个快速的解决方案。请注意，确切的内存需求将会更大——它们取决于字符串长度、Python列表的实现细节(它们使用的指数缩放算法会增加额外的开销以使列表快速调整大小)等。

另外，还有一点与算法无关的评论。从统计角度来看，计算 360981000 个 p 值几乎没有意义。请记住，为了使此过程科学有效，您很可能需要对 p 值应用多重检验校正。对于如此多的 p 值，校正将是巨大的，并且可能会使您从该过程中获得的任何见解无效(例如，要获得 0.05 的调整 p 值，您需要观察低于 0.0000000002 的未调整 p 值)。