将fastq信息复制到新fastq文件中
我正在尝试写一段代码,这段代码会打开一个fasta文件,从另一个fastq文件中提取读取名称(标题)、序列(seq)和质量分数(qual),但只有当这些信息在fasta文件中找到时,才会提取,并把这些fastq的信息写入一个新的fastq文件。不过,我在写这最后一步的时候遇到了困难(我在代码中用粗体标出了我遇到问题的地方)。有没有人知道怎么写这一部分,或者我可以在哪里找到关于如何在python中输入这些内容的信息?
到目前为止,我有:
from sys import argv
from Bio.SeqIO.QualityIO import FastqGeneralIterator
script, merged_seqs, raw_seqs = argv
merged_from_raw = "merged_only.fastq"
merged_names = set()
for line in open(merged_seqs):
if line[0] == ">":
read_name = line.split()[0][1:]
merged_names.add(read_name)
raw_fastq = raw_seqs
temp_handle = open(merged_from_raw, "w")
for title, seq, qual in FastqGeneralIterator(open(raw_fastq)) :
if title in merged_names:
**handle.write() #this is where I don't know how to write what I need in python**
1 个回答
1
除非你有特别的理由自己去实现文件解析,否则使用SeqIO解析器来处理你的输入和输出文件会更好。你可以试试下面这样的代码(注意:我之前没有用过Bio库,也没有测试过这段代码):
from sys import argv
from Bio import SeqIO
output_filename = 'merged_only.fastq'
merged_seqs, raw_seqs = argv[1:2]
# Get fasta iterator, and read source fastq file into a dict-like object
merged_names = SeqIO.parse(merged_seqs, 'fasta')
source_seqs = SeqIO.index(raw_seqs, 'fastq')
filtered_seqs = (source_seqs[record.id] for record in merged_names if record.id in source_seqs)
SeqIO.write(filtered_seqs, output_filename, 'fastq')