我正在尝试使用 python 在两个单独的文件中找到感兴趣的四行 block ,然后按受控顺序打印出其中一些行。下面是两个输入文件和所需输出文件的示例。请注意,Input.fasta 中的 DNA 序列与 Input.fastq 中的 DNA 序列不同,因为 .fasta 文件已被读取更正。
输入.fasta
>read1
AAAGGCTGT
>read2
AGTCTTTAT
>read3
CGTGCCGCT
输入.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 操作”
请您参考如下方法:
我建议你使用 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
