提高Python脚本效率的方法

Question

我有一个基因列表、它们的坐标和它们的表达（现在只查看前 500 个最高表达的基因）和对应于 DNA 读数的 12 个文件。我有一个 python 脚本，用于搜索与每个基因的坐标重叠的读数并将值存储在字典中。然后我使用这个字典创建一个 Pandas 数据框并将其保存为 csv。 （我将使用这些来创建散点图。）

RNA 文件看起来像这样（标题是基因名称、染色体、开始、停止、基因覆盖/富集）：

MSTRG.38        NC_008066.1     9204    9987    48395.347656
MSTRG.36        NC_008066.1     7582    8265    47979.933594
MSTRG.33        NC_008066.1     5899    7437    43807.781250
MSTRG.49        NC_008066.1     14732   15872   26669.763672
MSTRG.38        NC_008066.1     8363    9203    19514.273438
MSTRG.34        NC_008066.1     7439    7510    16855.662109

DNA 文件看起来像这样（标题是染色体、开始、停止、基因名称、覆盖范围、链）：

JQ673480.1      697     778     SRX6359746.5505370/2    8       +
JQ673480.1      744     824     SRX6359746.5505370/1    8       -
JQ673480.1      1712    1791    SRX6359746.2565519/2    27      +
JQ673480.1      3445    3525    SRX6359746.7028440/2    23      -
JQ673480.1      4815    4873    SRX6359746.6742605/2    37      +
JQ673480.1      5055    5092    SRX6359746.5420114/2    40      -
JQ673480.1      5108    5187    SRX6359746.2349349/2    24      -
JQ673480.1      7139    7219    SRX6359746.3831446/2    22      +

RNA 文件有 >9,000 行，DNA 文件有 >12,000,000 行。

我最初有一个 for 循环，可以一次性生成包含所有 12 个文件的所有值的字典，但它运行得非常慢。由于我可以访问具有多核的计算系统，因此我决定运行一个脚本，该脚本一次只计算一个 DNA 文件的覆盖率，如下所示：

#import modules
import csv
import pandas as pd
import matplotlib.pyplot as plt

#set sample name
sample='CON-2'
#set fraction number
f=6

#dictionary to store values
d={}

#load file name into variable
fileRNA="top500_R8_7-{}-RNA.gtf".format(sample)
print(fileRNA)
#read tsv file
tsvRNA = open(fileRNA)
readRNA = csv.reader(tsvRNA, delimiter="\t")
expGenes=[]
#convert tsv file into Python list
for row in readRNA:
    gene=row[0],row[1],row[2],row[3],row[4]
    expGenes.append(row)
#print(expGenes)

#establish file name for DNA reads
fileDNA="D2_7-{}-{}.bed".format(sample,f)
print(fileDNA)
tsvDNA = open(fileDNA)
readDNA = csv.reader(tsvDNA, delimiter="\t")
#put file into Python list
MCNgenes=[]
for row in readDNA:
    read=row[0],row[1],row[2]
    MCNgenes.append(read)

#find read counts
for r in expGenes:
    #include FPKM in the dictionary
    d[r[0]]=[r[4]]
    regionCount=0
    #set start and stop points based on transcript file
    chr=r[1]
    start=int(r[2])
    stop=int(r[3])
    #print("start:",start,"stop:",stop)
    for row in MCNgenes:
        if start < int(row[1]) < stop:
            regionCount+=1
    d[r[0]].append(regionCount)
n+=1
df=pd.DataFrame.from_dict(d)
#convert to heatmap
df.to_csv("7-CON-2-6_forHeatmap.csv")

不过，这个脚本也运行得很慢。我可以进行任何更改以使其更有效地运行吗？

Answer 1

如果我理解正确并且您尝试在不同文件中的基因坐标之间进行匹配，我相信最好的选择是使用 KDTree 分区算法之类的东西。

您可以使用 KDtree 对 DNA 和 RNA 数据进行分区。我假设您使用“开始”和“停止”作为“坐标”：

import pandas as pd
import numpy as np
from sklearn.neighbors import KDTree

dna = pd.DataFrame() # this is your dataframe with DNA data
rna = pd.DataFrame() # Same for RNA

# Let's assume you are using 'start' and 'stop' columns as coordinates
dna_coord = dna.loc[:, ['start', 'stop']]
rna_coord = rna.loc[:, ['start', 'stop']]

dna_kd = KDTree(dna_coord)
rna_kd = KDTree(rna_coord)

# Now you can go through your data and match with DNA:
my_data = pd.DataFrame()

for start, stop in zip(my_data.start, my_data.stop):
    coord = np.array(start, stop)
    dist, idx = dna_kd.query(coord, k=1)

    # Assuming you need an exact match
    if np.islose(dist, 0):
        # Now that you have the index of the matchin row in DNA data 
        # you can extract information using the index and do whatever
        # you want with it

        dna_gene_data = dna.loc[idx, :]

您可以调整搜索参数以获得所需的结果，但这会比每次搜索要快得多。

Answer 2

一般来说，Python 非常容易使用，但代价是效率低下！科学库（例如 pandas 和 numpy）在这里提供帮助，只需支付 Python 开销最少的有限次数，以将工作映射到方便的空间，然后用更高效的语言进行“繁重的工作”（这可能会很痛苦） /不方便使用）。

一般建议

• 尽可能尝试将数据放入数据帧并将其保留在其中（不要将数据转换为一些中间 Python 对象，例如 list 或 dict）
• 尝试使用数据框的方法或部分数据框来完成工作（例如.apply() 和.map()-like 方法）
• 当您必须在原生 Python 中进行迭代时，请在数据框的较短一侧进行迭代（即，可能只有 10 列，但有 10,000 行；遍历列）

这里有更多关于这个主题的信息：

How to iterate over rows in a DataFrame in Pandas?
Answer: DON'T*!

---

一旦您有了一个程序，您就可以通过收集运行时信息来对其进行基准测试。为此有很多库，但也有一个名为 cProfile 的内置库，它可能对您有用。

文档：https://docs.python.org/3/library/profile.html

python3 -m cProfile -o profile.out myscript.py