python - python中大文件的高效文件缓冲&扫描方法

Question

我遇到的问题的描述有点复杂，我宁愿提供更完整的信息。对于不耐烦的人，这里是我可以总结的最简短的方式:

What is the fastest (least execution time) way to split a text file in to ALL (overlapping) substrings of size N (bound N, eg 36) while throwing out newline characters.

我正在编写一个模块来解析基于 FASTA ascii 的基因组格式的文件。这些文件包含所谓的“hg18”人类引用基因组，您可以从 UCSC genome browser 下载。 (去鼻涕虫!)如果你喜欢的话。

正如您所注意到的，基因组文件由 chr[1..22].fa 和 chr[XY].fa 以及一组未在该模块中使用的其他小文件组成。

已经存在多个用于解析 FASTA 文件的模块，例如 BioPython 的 SeqIO。 (抱歉，我会发布一个链接，但我还没有这样做的要点。)不幸的是，我能够找到的每个模块都没有执行我正在尝试执行的特定操作。

我的模块需要将基因组数据(例如，“CAGTACGTCAGACTATACGGAGCTA”可能是一行)拆分为每个重叠的 N 长度子字符串。让我举个例子，使用一个非常小的文件(实际的染色体文件长度在 355 到 2000 万个字符之间)并且 N=8

>>>import cStringIO>>>example_file = cStringIO.StringIO("""\>headerCAGTcagTFgcACF""")>>>for read in parse(example_file):...    print read...CAGTCAGTFAGTCAGTFGGTCAGTFGCTCAGTFGCACAGTFGCACAGTFGCACF

The function that I found had the absolute best performance from the methods I could think of is this:

def parse(file):
  size = 8 # of course in my code this is a function argument
  file.readline() # skip past the header
  buffer = ''
  for line in file:
    buffer += line.rstrip().upper()
    while len(buffer) >= size:
      yield buffer[:size]
      buffer = buffer[1:]

这可行，但不幸的是，以这种方式解析人类基因组仍需要大约 1.5 小时(见下面的注释)。也许这是我将看到的最好的这种方法(完整的代码重构可能是有序的，但我想避免它，因为这种方法在代码的其他领域有一些非常具体的优势)，但我我想我会把它交给社区。

谢谢!

• 请注意，这次包括很多额外的计算，例如计算反向链读取以及对大小约为 5G 的哈希进行哈希表查找。

回答后结论:事实证明，与程序的其余部分，因此我使用了该方法。谢谢大家!

Answer 1

你能映射文件并开始使用滑动窗口浏览它吗？我写了一个运行起来非常小的愚蠢的小程序:

USER       PID %CPU %MEM    VSZ   RSS TTY      STAT START   TIME COMMAND
sarnold  20919  0.0  0.0  33036  4960 pts/2    R+   22:23   0:00 /usr/bin/python ./sliding_window.py

处理一个 636229 字节的 fasta 文件(通过 http://biostar.stackexchange.com/questions/1759 找到)花费了 0.383 秒。

#!/usr/bin/python

import mmap
import os

  def parse(string, size):
    stride = 8
    start = string.find("\n")
    while start < size - stride:
        print string[start:start+stride]
        start += 1

fasta = open("small.fasta", 'r')
fasta_size = os.stat("small.fasta").st_size
fasta_map = mmap.mmap(fasta.fileno(), 0, mmap.MAP_PRIVATE, mmap.PROT_READ)
parse(fasta_map, fasta_size)