0

我正在阅读格式的大型基因文件

>GeneID
ACTCTCTCTATATATATATAT\n
GCTCTGCTTCTAGAGAGAGTG\n
TCTATTTGTTTATATATCTTT\n
>GeneID
GCTCTGCTTCTAGAAATTCCC\n
ACTCTGTATATATTTTCAAAA\n
GCTCTGCTTCTAGAGAGAGTG\n

每个基因都以 > 开头，然后是唯一的 ID。之后是该基因的核苷酸线。不幸的是，这个文件的生成使得每行序列之间都有换行符。

我需要在每个序列中读取一个连续的字符串。所以，我一直在使用下一种方法（如下图）。

for line in filer:
    if line.startswith(">"):

        # Find Sequences
        seq_seg = next(filer)
        seq = ""

        # Concatenate lines until find next gene
        while not (seq_seg.startswith(">")):
            seq += seq_seg.strip()  # Get rid of '\n'
            seq_seg = next(filer)

我发现我的脚本只提取了文件中的一半基因，因为当在导致 while 循环失败的条件下调用 next 时，文件指针指向下一个基因 ID，然后当 for 的下一次迭代时循环执行，它移动到下一个文件。

有没有办法将文件指针倒回上一行，所以我的 for 循环将其作为新基因捕获？

我见过类似的问题，但没有一个解决我通过文件读取的具体方式

  for line in file:
        #do stuff
4

4 回答 4

3

我会使用生成器而不是跳过行（有人告诉我这可以大大简化）：

def parse_file(file):
    id = ''
    gene = ''

    for line in file:
        if line.startswith('>'):
            if gene:
                yield id, gene

            id = line[1:]
            gene = ''
        else:
            gene += line.strip()
    else:
        yield id, gene # Final gene

现在，您只需几行代码就可以安全地迭代整个事情：

with open('file.txt', 'r') as handle:
    for gene_id, nucleotides in parse_file(handle):
        print gene_id, nucleotides

还有pyfasta。

或更通用的功能itertools：

def grouper(line):
    return line.startswith('>') and line[1:]

def itersplit(it, pred):
    groups = (list(group) for key, group in itertools.groupby(it, pred))
    yield from zip(groups, groups)

def parse(file):
    for key, group in itersplit(file, grouper):
        yield key[0], ''.join(group)
于 2013-07-16T21:04:55.863 回答
1

这是另一种使用reand的方法mmap：

import mmap, re

with open(your_file) as fin:
    mm = mmap.mmap(fin.fileno(), 0, access=mmap.ACCESS_READ)
    for match in re.finditer('>([^\n]+)([^>]*)', mm, flags=re.DOTALL):
        print match.group(1), match.group(2).replace('\n', '')

#GeneID1 ACTCTCTCTATATATATATATGCTCTGCTTCTAGAGAGAGTGTCTATTTGTTTATATATCTTT
#GeneID2 GCTCTGCTTCTAGAAATTCCCACTCTGTATATATTTTCAAAAGCTCTGCTTCTAGAGAGAGTG

这样做可以将整个文件视为一个字符串，但将利用操作系统按需提供文件的一部分来完成正则表达式。当它使用时finditer，我们也没有在内存中建立结果集。

于 2013-07-16T21:47:21.890 回答
1

有没有办法将文件指针倒回上一行，所以我的 for 循环将其作为新基因捕获？

在 Python 3 中，没有。您不能将文件迭代与对文件指针的显式操作混合使用。

在 Python 2 中，也许。但这只是偶然起作用，这就是它在 3.0 中被禁止的原因，并且不能保证在每种情况下都能正常工作。所以，你不应该这样做。

一个更好的方法是询问如何回退迭代器。答案是itertools。您可以使用tee. 您可以将迭代器重新绑定到chain([pushed_back_value], iterator). 等等。

但是，正如其他人所指出的，有一种更好的方法可以做到这一点。您实际上不需要向前窥视或向后寻找，您只需要对事物进行分组。你也可以这样做itertools，但在这种情况下，这很简单，你不妨明确地这样做，正如 Blender 所示。

于 2013-07-16T21:10:29.690 回答
0

有很多更简单的方法可以读取 fasta 文件，例如：

entries = []
for line in filer:
    if line.startswith('>'):
        entries.append((line.rstrip()[1:], []))
    else:
        entries[-1][1].append(line.rstrip())

这将为您提供一个元组列表。第一个元素是序列 ID，第二个元素是序列列表。

在此之后加入序列很容易：

entries = [(x, "".join(y)) for x,y in entries]
于 2013-07-16T21:02:04.420 回答