python - 如何创建索引来解析大文本文件

Question

我有两个 FASTQ 格式的文件 A 和 B，它们基本上是几亿行文本，以 @ 开头，每行 4 行，如下所示:

@120412_SN549_0058_BD0UMKACXX:5:1101:1156:2031#0/1
GCCAATGGCATGGTTTCATGGATGTTAGCAGAAGACATGAGACTTCTGGGACAGGAGCAAAACACTTCATGATGGCAAAAGATCGGAAGAGCACACGTCTGAACTCN
+120412_SN549_0058_BD0UMKACXX:5:1101:1156:2031#0/1
bbbeee_[_ccdccegeeghhiiehghifhfhhhiiihhfhghigbeffeefddd]aegggdffhfhhihbghhdfffgdb^beeabcccabbcb`ccacacbbccB

我需要比较

5:1101:1156:2031#0/

将文件 A 和 B 之间的部分分开，并将与新文件匹配的 4 行组写入文件 B 中。我在 python 中得到了一段代码，它可以做到这一点，但只适用于小文件，因为它为文件 A 中的每个 @ 行解析文件 B 的整个 @ 行，并且两个文件都包含数亿行。

有人建议我应该为文件B创建一个索引；我用谷歌搜索但没有成功，如果有人能指出如何做到这一点或让我知道一个教程以便我可以学习，我将非常感激。谢谢。

==编辑== 理论上，每组 4 行在每个文件中只应存在一次。如果在每次匹配后中断解析，它是否会提高足够的速度，或者我是否需要完全不同的算法？

Answer 1

索引只是您正在使用的信息的缩写版本。在这种情况下，您将需要“键” - @ 行上的第一个冒号 (':') 和靠近末尾的最后一个斜杠 ('/') 之间的文本 - 以及某种值。

由于本例中的“值”是 4 行 block 的全部内容，并且由于我们的索引将为每个 block 存储一个单独的条目，因此如果我们使用指数中的实际值。

相反，我们使用 4 行 block 开头的文件位置。这样，您可以移动到该文件位置，打印 4 行，然后停止。总成本是存储整数文件位置所需的 4 或 8 或任意字节，而不是实际基因组数据的任意字节。

这里有一些代码可以完成这项工作，但也进行了大量的验证和检查。您可能想扔掉不使用的东西。

import sys

def build_index(path):
    index = {}
    for key, pos, data in parse_fastq(path):
        if key not in index:
            # Don't overwrite duplicates- use first occurrence.
            index[key] = pos

    return index

def error(s):
    sys.stderr.write(s + "\n")

def extract_key(s):
    # This much is fairly constant:
    assert(s.startswith('@'))
    (machine_name, rest) = s.split(':', 1)
    # Per wikipedia, this changes in different variants of FASTQ format:
    (key, rest) = rest.split('/', 1)
    return key

def parse_fastq(path):
    """
    Parse the 4-line FASTQ groups in path.
    Validate the contents, somewhat.
    """
    f = open(path)
    i = 0
    # Note: iterating a file is incompatible with fh.tell(). Fake it.
    pos = offset = 0
    for line in f:
        offset += len(line)
        lx = i % 4
        i += 1
        if lx == 0:     # @machine: key
            key = extract_key(line)
            len1 = len2 = 0
            data = [ line ]
        elif lx == 1:
            data.append(line)
            len1 = len(line)
        elif lx == 2:   # +machine: key or something
            assert(line.startswith('+'))
            data.append(line)
        else:           # lx == 3 : quality data
            data.append(line)
            len2 = len(line)

            if len2 != len1:
                error("Data length mismatch at line "
                        + str(i-2)
                        + " (len: " + str(len1) + ") and line "
                        + str(i)
                        + " (len: " + str(len2) + ")\n")
            #print "Yielding @%i: %s" % (pos, key)
            yield key, pos, data
            pos = offset

    if i % 4 != 0:
        error("EOF encountered in mid-record at line " + str(i));

def match_records(path, index):
    results = []
    for key, pos, d in parse_fastq(path):
        if key in index:
            # found a match!
            results.append(key)

    return results

def write_matches(inpath, matches, outpath):
    rf = open(inpath)
    wf = open(outpath, 'w')

    for m in matches:
        rf.seek(m)
        wf.write(rf.readline())
        wf.write(rf.readline())
        wf.write(rf.readline())
        wf.write(rf.readline())

    rf.close()
    wf.close()

#import pdb; pdb.set_trace()
index = build_index('afile.fastq')
matches = match_records('bfile.fastq', index)
posns = [ index[k] for k in matches ]
write_matches('afile.fastq', posns, 'outfile.fastq')

请注意，此代码返回到第一个文件以获取数据 block 。如果文件之间的数据相同，则当发生匹配时，您将能够从第二个文件复制该 block 。

另请注意，根据您尝试提取的内容，您可能希望更改输出 block 的顺序，并且您可能希望确保键是唯一的，或者可能确保键不是唯一的，但按匹配顺序重复。这取决于你 - 我不确定你要如何处理这些数据。