Extract fastq reads from reads list

Question

[felipelira3]

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[1], "r")
fastq_file = open(sys.argv[2], "r")

recover = set()

for line in list_of_reads:
        if line not in recover:
                read_name = line.split()[0]
                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()

[dschika]

First guess: line in your first for loop ends with '\n', while title in the second for loop does not: add:

line = line.strip()

in the first loop before the if ... .

[felipelira3]

Wow a beginner mistake. Thank you @dschika

[WouterDeCoster]

I think (but this should be tested) that it would be quicker to just loop over the list_of_reads and create a list, and create a set out of that later. Something like

recover = set([line.split()[0] for line in list_of_reads])

[Brian Bushnell]

If performance is important, there are existing tools that will do the job:

A: Extracting fastq files, based on their fasta counterparts

filterbyname will alternatively accept a list of read names as the argument to "names", rather than a fasta file.

Answer 1

[felipelira3]

improved script

#usage: script [list_of_reads] [fastq_file] [-output_file_name]

import sys
from Bio import SeqIO

list_of_reads= open(sys.argv[1], "r")
recover = set([line.strip().split(' ')[0] for line in list_of_reads])

output = open(sys.argv[3], "w")         # uncomment to write the output file
for record in SeqIO.parse(sys.argv[2], "fastq"):
        if record.id in recover:
                output.write(record.format("fastq"))
output.close()

[WouterDeCoster]

Looks good, just another piece of advice: it's better (and more convenient) to use:

with open(sys.argv[1], 'r') as list_of reads:
    recover = set([line.strip().split(' ')[0] for line in list_of_reads])

The with statement takes care of opening and closing the file (which you didn't do for sys.argv[1])