超越forlooping：一个大的、格式良好的数据fi的高性能解析

Question

我正在寻找使用python优化大数据解析问题的性能。如果有人感兴趣：下面的数据是六种灵长类动物全基因组DNA序列比对的片段。在

目前，我知道如何处理这类问题的最好方法是打开我的每个~250（大小为20-50MB）的文件，逐行循环并提取我想要的数据。格式（如示例所示）是相当规则的，尽管在每10-10万行线段上都有重要的变化。循环工作很好，但很慢。在

我最近一直在使用numpy来处理大量（>10gb）的数值数据集，我对我能够在数组上执行不同计算的速度印象深刻。我想知道是否有一些处理格式化文本的高性能解决方案可以避免循环的繁琐？在

我的文件包含多个具有模式的段

<MULTI-LINE HEADER>  # number of header lines mirrors number of data columns
<DATA BEGIN FLAG>  # the word 'DATA'
<DATA COLUMNS>  # variable number of columns
<DATA END FLAG>  # the pattern '//'
<EMPTY LINE>

示例：

^{pr2}$

我将使用段，其中在头中同时存在homo_sapiens和{}，字段6（我在上面的注释中称为质量标志）等于每个物种的“1”。由于macaca_mulatta没有出现在第二个示例中，我将完全忽略此段。在

我只关心segment_start和segment_end的坐标，因此在存在homo_sapiens的段中，我将记录这些字段并将它们用作dict()的键。segment_start还告诉我homo_sapiens的第一个位置坐标，对于当前段中的每一行数据，它严格增加1。在

我想比较homo_sapiens和{}的字母（DNA碱基）。出现homo_sapiens和{}的标题行（即第一示例中的1和5）对应于表示它们各自序列的数据列。在

重要的是，这些列并不总是相同的，因此我需要检查页眉以获得每个片段的正确索引，并检查两个物种在当前片段中是否相等。

看一下示例1中的两行数据，对我来说相关的信息是

# homo_sapiens_coordinate homo_sapiens_base macaca_mulatta_base
11388669 G G
11388670 C T

对于每个包含homo_sapiens和macaca_mulatta信息的段，我将从头和两个不匹配的每个位置记录{}的开始和结束。最后，有些职位有“差距”或质量较低的数据，即

aaa--A

我只记录来自homo_sapiens和macaca_mulatta都有有效基的位置（必须在集合ACGT）中，因此我考虑的最后一个变量是每个段的有效基的计数器。在

给定文件的最终数据结构是一个字典，如下所示：

{(segment_start=i, segment_end=j, valid_bases=N): list(mismatch positions), 
    (segment_start=k, segment_end=l, valid_bases=M): list(mismatch positions), ...}

这里写了一个循环，用于执行这个函数：

def human_macaque_divergence(chromosome):

    """
    A function for finding the positions of human-macaque divergent sites within segments of species alignment tracts
    :param chromosome: chromosome (integer:
    :return div_dict: a dictionary with tuple(segment_start, segment_end, valid_bases_in_segment) for keys and list(divergent_sites) for values
    """

    ch = str(chromosome)
    div_dict = {}

    with gz.open('{al}Compara.6_primates_EPO.chr{c}_1.emf.gz'.format(al=pd.align, c=ch), 'rb') as f:

        # key to the header fields:
        # header_flag chromosome segment_start segment_end quality_flag chromosome_info
        # SEQ homo_sapiens 1 14163 24841 1 (chr_length=249250621)

        # flags, containers, counters and indices:
        species   = []
        starts    = []
        ends      = []
        mismatch  = []

        valid        = 0
        pos          = -1
        hom          = None
        mac          = None

        species_data = False  # a flag signalling that the lines we are viewing are alignment columns

        for line in f:

            if 'SEQ' in line:  # 'SEQ' signifies a segment info field

                assert species_data is False
                line = line.split()

                if line[2] == ch and line[5] == '1':  # make sure that the alignment is to the desired chromosome in humans quality_flag is '1'

                    species += [line[1]]  # collect each species in the header
                    starts  += [int(line[3])]  # collect starts and ends
                    ends    += [int(line[4])]

            if 'DATA' in line and {'homo_sapiens', 'macaca_mulatta'}.issubset(species):

                species_data = True

                # get the indices to scan in data columns:
                hom       = species.index('homo_sapiens') 
                mac       = species.index('macaca_mulatta')
                pos       = starts[hom]  # first homo_sapiens positional coordinate

                continue

            if species_data and '//' not in line:

                assert pos > 0

                # record the relevant bases:
                human   = line[hom]
                macaque = line[mac]

                if {human, macaque}.issubset(bases):
                    valid += 1

                if human != macaque and {human, macaque}.issubset(bases):
                    mismatch += [pos]

                pos += 1

            elif species_data and '//' in line:  # '//' signifies segment boundary

                # store segment results if a boundary has been reached and data has been collected for the last segment:
                div_dict[(starts[hom], ends[hom], valid)] = mismatch

                # reset flags, containers, counters and indices
                species   = []
                starts    = []
                ends      = []
                mismatch  = []

                valid        = 0
                pos          = -1
                hom          = None
                mac          = None
                species_data = False

            elif not species_data and '//' in line:

                # reset flags, containers, counters and indices
                species   = []
                starts    = []
                ends      = []

                pos       = -1
                hom       = None
                mac       = None

    return div_dict

这段代码运行得很好（也许需要一些调整），但我真正的问题是，是否有一种更快速的方法可以在不运行for循环和检查每一行的情况下提取这些数据？例如，使用f.read()加载整个文件只需不到一秒钟的时间，尽管这会创建一个相当复杂的字符串。（原则上，我假设我可以使用正则表达式来解析至少一些数据，比如头信息，但是我不确定如果不使用批量方法来处理每个段中的每个数据列，这是否一定会提高性能）。在

有人对我如何绕过数十亿行的循环并以更大的方式解析这种文本文件有什么建议吗？在

请让我知道，如果有什么不清楚的评论，高兴编辑或直接回复改善帖子！在

Answer 1

当您有工作代码并且需要提高性能时，请使用探查器，一次测量一个优化的效果。（即使不使用profiler，也一定要使用后者。）您当前的代码看起来很合理，也就是说，从性能上看，我没有看到任何“愚蠢”的地方。在

尽管如此，对所有字符串匹配使用预编译正则表达式可能是值得的。通过使用re.MULTILINE，您可以将整个文件作为字符串读入并拉出部分行。例如：

s = open('file.txt').read()
p = re.compile(r'^SEQ\s+(\w+)\s+(\d+)\s+(\d+)\s+(\d+)', re.MULTILINE)
p.findall(s)

产生：

^{pr2}$

然后，您将需要对这些数据进行后处理，以处理代码中的特定条件，但总体结果可能更快。在

Answer 2

您的代码看起来不错，但是还有一些特殊的地方需要改进，比如使用map，等等

有关良好性能的提示，请参见Python指南：

https://wiki.python.org/moin/PythonSpeed/PerformanceTips

我已经用上面的技巧让代码运行得几乎和C代码一样快。基本上，尽量避免for循环（使用map），尝试使用find内置函数，等等。通过使用Python的内置函数（大部分是用C编写的），使Python尽可能地为您工作

获得可接受的性能后，可以使用以下方法并行运行：

https://docs.python.org/dev/library/multiprocessing.html#module-multiprocessing

编辑：

我还刚刚意识到您正在打开一个压缩的gzip文件。我怀疑花了很多时间去减压。您可以尝试使用多线程来加快速度：

https://code.google.com/p/threadzip/

Answer 3

是的，您可以使用一些正则表达式一次性提取数据；这可能是工作/性能的最佳比率。在

如果您需要更多的性能，您可以使用mx.TextTools来构建一个有限状态机；我很有信心这将大大加快速度，但是编写规则和学习曲线所需的工作量可能会使您望而却步。在

您还可以将数据分成块并并行处理，这可能会有所帮助。在