This is my code to extract the reads from a raw fastq file. I have one list with the headers of the reads and I want to retreive them from the original file. My output is an empty file. Maybe I did it not so well. Any help, improving this script will be appreciate. Less is more.
import sys from Bio import SeqIO from Bio.SeqIO.QualityIO import FastqGeneralIterator list_of_reads = open(sys.argv, "r") fastq_file = open(sys.argv, "r") recover = set() for line in list_of_reads: if line not in recover: read_name = line.split() recover.add(read_name) # add value to set output = open("recovered.fastq", "w") for title, seq, qual in FastqGeneralIterator(fastq_file): if title in recover: output.write("@%s\n%s\n+\n%s\n" % (title, seq, qual)) output.close()