python - 超越循环:高性能，大格式的数据文件解析

Question

我希望优化我使用python时遇到的大数据解析问题的性能。以防有人感兴趣：下面显示的数据是六种灵长类动物全基因组DNA序列比对的片段。
目前，我知道如何处理这类问题的最好方法是打开我的~250（大小20-50MB）文件，逐行循环并提取我想要的数据。格式（如示例所示）是相当规则的，尽管在每个10000到100000行的线段上都有重要的更改。循环工作很好，但速度很慢。
我最近一直在使用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>

例子：

# key to the header fields:
# header_flag chromosome segment_start segment_end quality_flag chromosome_data
SEQ homo_sapiens 1 11388669 11532963 1 (chr_length=249250621)
SEQ pan_troglodytes 1 11517444 11668750 1 (chr_length=229974691)
SEQ gorilla_gorilla 1 11607412 11751006 1 (chr_length=229966203)
SEQ pongo_pygmaeus 1 218866021 219020464 -1 (chr_length=229942017)
SEQ macaca_mulatta 1 14425463 14569832 1 (chr_length=228252215)
SEQ callithrix_jacchus 7 45949850 46115230 1 (chr_length=155834243)
DATA
GGGGGG
CCCCTC
......  # continue for 10-100 thousand lines
//

SEQ homo_sapiens 1 11345717 11361846 1 (chr_length=249250621)
SEQ pan_troglodytes 1 11474525 11490638 1 (chr_length=229974691)
SEQ gorilla_gorilla 1 11562256 11579393 1 (chr_length=229966203)
SEQ pongo_pygmaeus 1 219047970 219064053 -1 (chr_length=229942017)
DATA
CCCC
GGGG
....  # continue for 10-100 thousand lines
//

<ETC>

我将使用头中同时存在物种homo_sapiens和macaca_mulatta的段，并且字段6（我在上面的注释中称之为质量标志）对于每个物种都等于“1”。因为macaca_mulatta不会出现在第二个示例中，所以我将完全忽略此段。
我只关心segment_start和segment_end坐标，所以在存在homo_sapiens的段中，我将记录这些字段并将它们用作homo_sapiens的键。dict()还告诉我segment_start的第一个位置坐标，当前段中每行数据的第一个位置坐标严格增加1。
我想比较homo_sapiens和homo_sapiens的字母（DNA碱基）。出现macaca_mulatta和homo_sapiens的标题行（即第一个示例中的1和5）对应于表示其各自序列的数据列。
重要的是，这些列并不总是相同的，所以我需要检查头以获得每个段的正确索引，并检查两个物种是否都在当前段中。
看一下示例1中的两行数据，我的相关信息是

# homo_sapiens_coordinate homo_sapiens_base macaca_mulatta_base
11388669 G G
11388670 C T

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

aaa--A

我只从homo_sapiens和homo_sapiens都有有效基（必须在集合macaca_mulatta中）的位置进行记录，因此我考虑的最后一个变量是每段有效基的计数器。
给定文件的最终数据结构是一个字典，它如下所示：

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

下面是我编写的使用for循环执行此操作的函数：

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循环并检查每一行？例如，使用ACGT加载整个文件只需不到一秒钟的时间，尽管它创建了一个相当复杂的字符串。（原则上，我假设我可以使用正则表达式来解析至少一些数据，例如头信息，但我不确定如果没有某种批量方法来处理每个段中的每个数据列，这是否一定会提高性能）。
有人对我如何绕过数十亿行的循环并以更大的方式解析这种文本文件有什么建议吗？
请让我知道，如果有任何不清楚的意见，乐意编辑或直接回应，以改善职位！

Answer 1

是的，您可以使用一些正则表达式一次性提取数据；这可能是工作/性能的最佳比率。
如果您需要更多的性能，可以使用mx.TextTools来构建一个有限状态机；我很有信心这将大大加快速度，但编写规则和学习曲线所需的努力可能会使您气馁。
您还可以将数据分成块并并行处理，这可能会有所帮助。