extracting reads from fastq file based on read_id

Question

[sumithrasank75]

I have a file with a set of read ids called tmp. I want to extract reads with these read ids from another fastq file called reads.fastq. When I use grep -f tmp reads.fastq I get only the first line (header) of the read. How can I get the output in full fastq format.

Answer 1

[GenoMax]

grep -A 3 seq_ID your file should get you the full fastq record. Iterate over the ID's you have.

filterbyname.sh from BBMap would also work with your list of ID's in a file.

[poisonAlien]

Today I learnt about -A. Thanks.

Answer 2

[Charles Plessy]

If your tmp and reads.fastq have their reads in the same order, you can try our syncpairs tools.

Answer 3

[brent_wilson]

For bash, I would use a nice one-liner:

for x in $(cat tmp); do reads.fastq | grep $x -A 3; done

Of course you can redirect to wherever you want:

for x in $(cat tmp); do reads.fastq | grep $x -A 3; done > tmp2

I hope that this helps!

Brent Wilson, PhD | Project Scientist | Cofactor Genomics

4044 Clayton Ave. | St. Louis, MO 63110 | tel. [phone]

Catch the latest from Cofactor on our blog.

[Devon Ryan]

For m desired reads out of a fastq file with n total reads, this takes O(mn) time, as opposed to the -f flag given to grep (or just using BBMap), which is likely O(n).