在python循环中从交替文件中打印行答案

【问题标题】：in python loop print lines from alternating files在python循环中从交替文件中打印行
【发布时间】：2018-03-01 14:44:36
【问题描述】：

我正在尝试使用 python 在两个单独的文件中查找感兴趣的四行块，然后以受控顺序打印出其中一些行。下面是两个输入文件和所需输出文件的示例。请注意，Input.fasta 中的 DNA 序列与 Input.fastq 中的 DNA 序列不同，因为 .fasta 文件已被读取更正。

Input.fasta

>read1
AAAGGCTGT
>read2
AGTCTTTAT
>read3
CGTGCCGCT

Input.fastq

@read1
AAATGCTGT
+
'(''%$'))
@read2
AGTCTCTAT
+
&---+2010
@read3
AGTGTCGCT
+
0-23;:677

DesiredOutput.fastq

@read1
AAAGGCTGT
+
'(''%$'))
@read2
AGTCTTTAT
+
&---+2010
@read3
CGTGCCGCT
+
0-23;:677

基本上我需要序列行“AAAGGCTGT”，来自“input.fasta”的“AGTCTTTAT”和“CGTGCCGCT”以及来自“input.fastq”的所有其他行。这允许将质量信息恢复为读取更正的 .fasta 文件。

这是我最近失败的尝试：

fastq = open(Input.fastq, "r")
fasta = open(Input.fasta, "r")

ReadIDs = []
IDs = []

with fastq as fq:
    for line in fq:
        if "read" in line:  
            ReadIDs.append(line)
            print(line.strip())
            for ID in ReadIDs:
                IDs.append(ID[1:6])
            with fasta as fa:
                for line in fa:
                    if any(string in line for string in IDs):
                        print(next(fa).strip())
            next(fq)
            print(next(fq).strip())
            print(next(fq).strip())

我认为我在尝试将“with”调用嵌套到同一个循环中的两个不同文件时遇到了麻烦。这会正确打印 read1 所需的行，但不会继续遍历剩余的行并抛出错误“ValueError: I/O operation on closed file”

【问题讨论】：

为什么只有读取 1 和 3 的序列 AAAGGCTGT 和 CGTGCCGCT 被复制到新文件中？
很抱歉给您带来了困惑。需要针对其各自读取的每个序列。我会更新以反映这一点。
with 语句是上下文管理器。当with 完成时，它会关闭文件。将with fasta as fa: 替换为with open(Input.fasta, "r") as fa:。好吧，也许这会解决你的问题。你想要的输出没有任何意义。

标签： python bioinformatics biopython fasta fastq

【解决方案1】：

我建议你使用Biopython，因为它为这些文件格式提供了很好的解析器，它不仅可以处理标准情况，还可以处理多行fasta。

这里是一个用对应的fasta序列行替换fastq序列行的实现：

from Bio import SeqIO

fasta_dict = {record.id: record.seq for record in
              SeqIO.parse('Input.fasta', 'fasta')}

def yield_records():
    for record in SeqIO.parse('Input.fastq', 'fastq'):
        record.seq = fasta_dict[record.id]
        yield record

SeqIO.write(yield_records(), 'DesiredOutput.fastq', 'fastq')

如果您不想使用标头而只依赖顺序，则解决方案更简单，内存效率更高（只需确保记录的顺序和数量相同），无需定义字典首先，一起迭代记录：

fasta_records = SeqIO.parse('Input.fasta', 'fasta')
fastq_records = SeqIO.parse('Input.fastq', 'fastq')

def yield_records():
    for fasta_record, fastq_record in zip(fasta_records, fastq_records):
        fastq_record.seq = fasta_record.seq
        yield fastq_record

【讨论】：

我看到了这在原则上和提供的示例序列上是如何工作的。但在实际大型 fasta 和 fastq 文件的实践中，我遇到了以下错误... raise ValueError("You must empty the letter annotations first!")。另外需要注意的是，此解决方案需要与示例中提供的相同的读取名称。
@Paul 我添加了一个使用顺序而不是读取名称的解决方案。关于你的错误，我认为这意味着 fastq 和 fasta 序列行的长度不同！这说明了为什么 Biopython 很有用，我猜你在某个地方有一个错误，如果序列长度不同，就不可能匹配质量分数
克里斯，感谢您指出这一点。你是对的，我的分析中有一个与这个脚本无关的皱纹。读取校正的 fasta 序列可能包含添加或删除的碱基，导致 fasta 和 fastq 文件之间的碱基数量不相等。感谢您的帮助和回答！

【解决方案2】：

## Open the files (and close them after the 'with' block ends)
with open("Input.fastq", "r") as fq, open("Input.fasta", "r") as fa:

    ## Read in the Input.fastq file and save its content to a list
    fastq = fq.readlines()

    ## Do the same for the Input.fasta file
    fasta = fa.readlines()


## For every line in the Input.fastq file
for i in range(len(fastq)):
    print(fastq[i]))
    print(fasta[2 * i])
    print(fasta[(2 * i) + 1])

【讨论】：

@Harvey 也可以，但我已将数据存储到 lists 以防 OP 稍后需要使用数据。

【解决方案3】：

我更喜欢 @Chris_Rands 的 Biopython solution 更适合小文件，但这里有一个解决方案，它只使用 Python 附带的电池并且内存高效。它假定 fasta 和 fastq 文件以相同的顺序包含相同数量的读取。

with open('Input.fasta') as fasta, open('Input.fastq') as fastq, open('DesiredOutput.fastq', 'w') as fo:
    for i, line in enumerate(fastq):
        if i % 4 == 1:
            for j in range(2):
                line = fasta.readline()
        print(line, end='', file=fo)

【讨论】：