6

我对 Python 很陌生，并且我编写了一个（可能非常难看的）脚本，该脚本应该从 fastq 文件中随机选择一个序列子集。fastq 文件将信息存储在每块四行的块中。每个块中的第一行以字符“@”开头。我用作输入文件的 fastq 文件为 36 GB，包含大约 14,000,000 行。

我试图重写一个已经存在的使用太多内存的脚本，并且我设法减少了很多内存使用量。但是脚本需要永远运行，我不明白为什么。

parser = argparse.ArgumentParser()
parser.add_argument("infile", type = str, help = "The name of the fastq input file.", default = sys.stdin)
parser.add_argument("outputfile", type = str, help = "Name of the output file.")
parser.add_argument("-n", help="Number of sequences to sample", default=1)
args = parser.parse_args()


def sample():
    linesamples = []
    infile = open(args.infile, 'r')
    outputfile = open(args.outputfile, 'w')
    # count the number of fastq "chunks" in the input file:
    seqs = subprocess.check_output(["grep", "-c", "@", str(args.infile)])
    # randomly select n fastq "chunks":
    seqsamples = random.sample(xrange(0,int(seqs)), int(args.n))
    # make a list of the lines that are to be fetched from the fastq file:
    for i in seqsamples:
        linesamples.append(int(4*i+0))
        linesamples.append(int(4*i+1))
        linesamples.append(int(4*i+2))
        linesamples.append(int(4*i+3))
    # fetch lines from input file and write them to output file.
    for i, line in enumerate(infile):
        if i in linesamples:
            outputfile.write(line)

grep 步骤几乎不需要任何时间，但 500 多分钟后，脚本仍未开始写入输出文件。所以我想这是 grep 和最后一个 for 循环之间的步骤之一，需要这么长时间。但我不明白究竟是哪一步，以及我能做些什么来加快它。

4

4 回答 4

2

根据 的大小linesamples，if i in linesamples将需要很长时间，因为您正在搜索列表中的每个迭代infile。您可以将其转换为 aset以缩短查找时间。此外，enumerate效率不是很高——我已经用line_num我们在每次迭代中递增的构造替换了它。

def sample():
    linesamples = set()
    infile = open(args.infile, 'r')
    outputfile = open(args.outputfile, 'w')
    # count the number of fastq "chunks" in the input file:
    seqs = subprocess.check_output(["grep", "-c", "@", str(args.infile)])
    # randomly select n fastq "chunks":
    seqsamples = random.sample(xrange(0,int(seqs)), int(args.n))
    for i in seqsamples:
        linesamples.add(int(4*i+0))
        linesamples.add(int(4*i+1))
        linesamples.add(int(4*i+2))
        linesamples.add(int(4*i+3))
    # make a list of the lines that are to be fetched from the fastq file:
    # fetch lines from input file and write them to output file.
    line_num = 0
    for line in infile:
        if line_num in linesamples:
            outputfile.write(line)
        line_num += 1
    outputfile.close()
于 2014-12-04T22:17:00.033 回答
1

你说 grep 运行得很快，所以在这种情况下，而不是仅仅使用 grep 来计算 @ 的出现次数，让 grep 输出它看到的每个 @ 字符的字节偏移量（使用-bgrep 的选项）。然后，用于random.sample选择您想要的任何块。一旦你选择了你想要的字节偏移量，使用infile.seek去每个字节偏移量并从那里打印出 4 行。

于 2014-12-04T22:21:21.650 回答
0

尝试并行化您的代码。我的意思是这个。您有 14,000,000 行输入。

  1. 处理你的 grep 并首先过滤你的行并将其写入 filtersInput.txt
  2. 将您的过滤输入拆分为 10.000-100.000 行文件，例如过滤输入001.txt、过滤输入002.txt
  3. 在这个拆分文件上工作我们的代码。将您的输出写入不同的文件，例如 output001.txt、output002.txt
  4. 合并您的结果作为最后一步。

由于您的代码根本不起作用。您也可以在这些过滤后的输入上运行您的代码。您的代码将检查过滤输入文件的存在，并将了解他所处的步骤，并从该步骤恢复。

您还可以使用 shell 或 python 线程以这种方式（在步骤 1 之后）使用多个 python 进程。

于 2014-12-04T22:16:52.520 回答
0

您可以使用Reservoir Sampling算法。使用此算法，您只需读取一次数据（无需提前计算文件的行数），因此您可以通过脚本传输数据。维基百科页面中有示例 python 代码。

Heng Li 的seqtk中还有一个用于 fastq 采样的 C 实现。

于 2015-03-20T23:28:27.440 回答