python - 提高循环数学脚本的效率

Question

我编写了一个脚本来计算分数，以及遗传谱列表的频率。

这里的遗传图谱由 SNP 的组合组成。每个 SNP 有两个等位基因。3 个 SNP 的输入文件如下所示：

    AA   CC   TT
    AT   CC   TT
    TT   CC   TT
    AA   CG   TT
    AT   CG   TT
    TT   CG   TT
    AA   GG   TT
    AT   GG   TT
    TT   GG   TT
    AA   CC   TA
    AT   CC   TA
    TT   CC   TA
    AA   CG   TA
    AT   CG   TA
    TT   CG   TA
    AA   GG   TA
    AT   GG   TA
    TT   GG   TA
    AA   CC   AA
    AT   CC   AA
    TT   CC   AA
    AA   CG   AA
    AT   CG   AA
    TT   CG   AA
    AA   GG   AA
    AT   GG   AA
    TT   GG   AA

然后我有以下代码，它可以接收上面的输入文件和一个包含权重和频率的表，如下所示：

(SNP        RiskAll  RefAll         OR            log(OR)    RiskAllFreq) # example header, not in file
SNP1             A       T       1.25    0.223143551314     0.97273 
SNP2             C       G       1.07    0.0676586484738    0.3     
SNP3             T       A       1.08    0.0769610411361    0.1136  

然后，它根据遗传谱中每个 SNP 的每个风险等位基因的对数优势比的总和以及基于乘以等位基因频率的频率计算分数，假设 Hardy Weinberg 平衡。

import sys

snp={}
riskall={}
weights={}
freqs={}    # effect allele, *MAY NOT BE MINOR ALLELE

pop = int(int(sys.argv[4]) + 4) # for additional columns due to additional populations. the example table given only has one population (column 6)

# read in OR table
pos = 0
with open(sys.argv[1], 'r') as f:
    for line in f:
        snp[pos]=(line.split()[0])
        riskall[line.split()[0]]=line.split()[1]
        weights[line.split()[0]]=line.split()[4]
        freqs[line.split()[0]]=line.split()[pop]

        pos+=1



### compute scores for each combination
with open(sys.argv[2], 'r') as f:
    for line in f:
        score=0
        freq=1
        for j in range(len(line.split())):
            rsid=snp[j]
            riskallele=riskall[rsid]
            frequency=freqs[rsid]
            wei=weights[rsid]
            allele1=line.split()[j][0]
            allele2=line.split()[j][1]
            if allele2 != riskallele:      # homozygous for ref
                score+=0
                freq*=(1-float(frequency))*(1-float(frequency))
            elif allele1 != riskallele and allele2 == riskallele:  # heterozygous, be sure that A2 is risk allele!
                score+=float(wei)
                freq*=2*(1-float(frequency))*(float(frequency))
            elif allele1 == riskallele: # and allele2 == riskall[snp[j]]:      # homozygous for risk, be sure to limit risk to second allele!
                score+=2*float(wei)
                freq*=float(frequency)*float(frequency)

            if freq < float(sys.argv[3]):   # threshold to stop loop in interest of efficiency 
                break

        print(','.join(line.split()) + "\t" + str(score) + "\t" + str(freq))

到目前为止，我已经设置了一个变量，我可以在其中指定一个阈值以在频率变得极低时打破循环，例如大约 1e-10。我希望将其扩大到至少包含 20 个 SNP。可以做哪些改进来加快脚本速度？

编辑：添加了表格文件的示例。

编辑 2：我现在尝试以 1e-4 的频率阈值运行脚本。到目前为止已经大约六天了，它仍在运行，而且太慢了，所以我正在寻找更多建议！

编辑3：澄清输入文件，标题实际上不在输入文件中，它们只是一个指示。

编辑 4：尝试使用 Pandas，但速度很慢，并且不确定是否可以为此目的进行矢量化。Dask 在我的 Unix 服务器上安装有问题。我现在已尽可能将所有数据结构更改为字典。我还可以做些什么？

标签： pythonperformancemathmemory