I would use python (no dependencies):
1. read you read names into list1 and change list to set (it's hashable, so checking for present of element is much faster than in list)
2. parse you fastq file and check if each name is present in set of read names of interest

#!/usr/bin/python
import sys

#get fname from parameter
idfile=sys.argv[1]

#load ids
ids = set( x.strip() for x in open(idfile) )

#read stdid 
handle=sys.stdin  

while ids:
  #parse fastq
  idline=handle.readline()
  seq   =handle.readline()
  spacer=handle.readline()
  quals =handle.readline()
  #check
  id=idline[:-1]
  if id in ids:
    #print fastq
    sys.stdout.write( '%s%s%s%s' % ( idline, seq, spacer, quals ) )
    #update ids
    ids.remove( id )

My python implementation is over 3x faster than c++ by Pierre. It's due to implementation of set (hash-able) instead of usual list.

#c++
time cat PL429_q10_filtered_1.fastq | ./parse PL429_q10_filtered_1.fastq.lids > PL429_q10_filtered_1.fastq.lids.fqcpp

real    4m21.137s
user    4m15.950s
sys 0m10.510s

#python (3.4x faster than c++)
time cat PL429_q10_filtered_1.fastq | python fastq_id_parser.py PL429_q10_filtered_1.fastq.lids > PL429_q10_filtered_1.fastq.fqSet

real    1m11.155s
user    1m5.270s
sys 0m10.920s

#python list instead of set (56x slower than python using set!)
time cat PL429_q10_filtered_1.fastq | python fastq_id_parser.noBioPython.noSet.py PL429_q10_filtered_1.fastq.lids > PL429_q10_filtered_1.fastq.lids.fqNoSet

real    70m54.534s
user    70m18.650s
sys 0m16.740s

#parsing file
time cat PL429_q10_filtered_1.fastq | wc -l
115959528

real    0m30.373s
user    0m1.990s
sys 0m8.250s

My test set was 10k: first and last 5k headers of the fq. Fq was 4.9GB, 115,959,528 lines, seq between 31-75bases. I used 15k SAS drive - read speed is ~200MB/s, access time ~5ms.

if you have only 10,000 reads you want to pull out, you can just iterate over the file once. For each record you check if the current read header is in the list (or hash) of 10,000. Something like this:

use Bio::SeqIO;
use Bio::Seq::Quality;

# read in 10K headers into a hash with value of 1 or whatever.
# you'll have to add your own code to do this.
my %headers = ();

$seqio  = Bio::SeqIO->new('-format'=>'fastq' , '-file'=>'some.fasq');

my $out_fastq = Bio::SeqIO->new( -format => 'fastq', '-file'=> 'subset.fastq');

while((my $rec = $seqio->next_seq())) {
    # if the current record doesn't exist in the 10,000, skip it.
    if(! $headers{$rec->display_id} ) { next; }
    # otherwise, write it to the output file.
    $out_fastq->write_seq($rec);
}

Here is my C++ solution:

#include <iostream>
#include <set>
#include <string>
#include <fstream>
#include <cstdlib>
using namespace std;


int main(int argc,char **argv)
    {
    string line;
    set<string> names;
    /* reads your names */
    ifstream in(argv[1],ios::in);
    if(!in.is_open()) return EXIT_FAILURE;
    while(getline(in,line,'\n'))
        {
        names.insert(line);
        }
    in.close();
    /* reads the fastq */
    string name;
    string seq;
    string name2;
    string qual;
    set<string>::iterator r;
    /* assuming 4 lines / read */
    while(!names.empty())
        {
        if(!getline(cin,name,'\n')) break;
        if(!getline(cin,seq,'\n')) break;
        if(!getline(cin,name2,'\n')) break;
        if(!getline(cin,qual,'\n')) break;
        r=names.find(name);
        if(r==names.end()) continue;
        names.erase(r);//names are unique we don't need this one anymore
        cout << name << "\n" << seq << "\n" << name2 << "\n" << qual << endl;
        }
    return 0;
    }

Compile:

g++ -Wall -O3 -o parse parser.cpp

Example:

$ cat names.txt
@EAS54_6_R1_2_1_540_792
@EAS54_6_R1_2_1_443_348

$ gunzip -c my.fastq.gz 
@EAS54_6_R1_2_1_413_324
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;88
@EAS54_6_R1_2_1_540_792
TTGGCAGGCCAAGGCCGATGGATCA
+
;;;;;;;;;;;7;;;;;-;;;3;83
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+EAS54_6_R1_2_1_443_348
;;;;;;;;;;;9;7;;.7;393333

$ gunzip -c  my.fastq.gz | ./parse names.txt
@EAS54_6_R1_2_1_540_792
TTGGCAGGCCAAGGCCGATGGATCA
+
;;;;;;;;;;;7;;;;;-;;;3;83
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+EAS54_6_R1_2_1_443_348
;;;;;;;;;;;9;7;;.7;393333

If all you need to do is extract the sequences, you could also just use grep, which will likely be about as fast as writing your own parser.

You can tell it to print out the line after the read name matched (-A 1) and can specify the list of patterns to search for in a file.

This doesn't really answer your question, but anyway: WTF would anyone store 70M reads in a text file and then try to access it as if it was a database?

I would convert the humongous FastQ to BAM. The conversion is lossless, saves space, and BAM/SAM are in fact easier to parse than FastQ. I don't know of a tool that would index BAM to retrieve reads by name, but if people start asking for that functionality, maybe someone will write such a tool.

Hi, I have written a fast parser for FastQ files. It's in Ruby but with a C extension underneath, so it's considerably faster then pure Ruby/Python/Perl code https://github.com/fstrozzi/bioruby-faster .

Look at the Wiki for simple examples on how to use it https://github.com/fstrozzi/bioruby-faster/wiki . It's easy and you get back an array for each entry in the FastQ file, with the complete sequence header (ID + comment), the sequence and quality values.

You can then parse the file once and just filter the sequence ID of each record by searching it against a Set of IDs using for example Ruby methods like Set#include?

my_ids = Set.new ["X,"Y","Z"]
fastq = Bio::Faster.new("sequences.fastq")
fastq.each_record do |sequence_header, sequence, quality|
     seq_id = sequence_header.split("\s").first
     puts sequence_header, sequence, quality if my_ids.include? seq_id
end

Sorry for a separate answer. It shows the table in my previous answer, which is malformated during the stackexchange transition. I cannot edit my old answer either as it exceeds the 5000 character limit. I will remove my this answer when I can reformat my old answer properly.

=========================================================================================
 Program        Lang       User  Sys   Mem/MB    Comments
-----------------------------------------------------------------------------------------
 cat            C           0.4  10.1     0.4    Unix cat
 wc -l          C           4.3  10.4     0.4    Unix wc -l
 fqextract.c    C          17.3  10.8     2.0    Hardcoded maximum line length: 4096
 Pierre-r2.cc   C++        19.5  12.2     3.8    Use tr1::unordered_set
 seqtk          C          19.9   8.9     4.1    Multi-line fastq; subregion; gzip input
 Pierre-r1.cc   C++        22.8  12.6     3.9    Replace cin with reading /dev/stdin
 wc             C          32.0  10.6     0.4    Unix wc (no command-line options)
 fqextract.py   Python     54.0   7.9     8.0    Same algorithm as fqextract.c
 fqextract.java Java       55.8   8.9   562.7    Same algorithm as fqextract.c
 fqextract.pl   Perl       71.5  10.5     5.5    Same algorithm as fqextract.c
 Leszek.py      Python     78.0  13.3     8.0    Infinite loop given improper input
 fqextract.lua  LuaJIT    109.6   8.6     6.1    Same algorithm as fqextract.c
 readfq.py      Python    122.9  11.3     8.1    Multi-line fastq
 bp-fqiter      Python    124.5   9.9    15.2    Multi-line fastq; biopython
 drio.rb        Ruby      130.1   6.9    10.0    Same algorithm as fqextract.c
 readfq.lua     Lua       189.5  11.5     6.0    Multi-line fastq
 readfq.pl      Perl      203.1  69.6     5.7    Multi-line fastq
 Pierre.cc      C++       323.6  13.3     3.9    Same algorithm as Leszek.py
 biopython.py   Python    928.8  34.2    15.6    Multi-line fastq
 bioperl.pl     Perl    >2400.0         >15.0    Multi-line fastq; killed
 grep -f -A     C       >2400.0        >414.0    Killed; possible false positives
=========================================================================================

Did anyone try do it out with:

• Go
• D
• Vala

You could also use a binary search. I have a library in Python that will only require you to write a small function to adapt to your file type...

Just a thought. I think Map/reduce in Hadoop might be super fast in doing this. Conceptually

Very simple solution using biopython:

def batch_iterator(iterator, batch_size:)
    entry = True 
    while entry :
        batch = []
        while len(batch) < batch_size :
            try :
                entry = iterator.next()
            except StopIteration :
                entry = None
            if entry is None :
                break
            batch.append(entry)
        if batch :
            yield batch

from Bio import SeqIO
record_iter = SeqIO.parse(open("filename.fastq"), "fastq")
for i, batch in enumerate(batch_iterator(record_iter, no.of records)) :
    filename = "group_%i.fastq" % (i+1)
    handle = open(filename, "w")
    count = SeqIO.write(batch, handle, "fastq")
    handle.close()
    print "Wrote %i records to %s" % (count, filename)

source: BiopythonWiki

> Other than this, there's also a BioPerl module (Bio::Index::fastq) which I didn't try because I read somewhere that you can run into the same problem when dealing with huge fastq files.

The idea is that you should never try to load that huge file in RAM. (you could do it if the file would be only half of installed amount of RAM). Just enumerate each record and see if it meets your criteria (the name you want). When you have found one record don't keep it in RAM. Write it to disk in a new file. An app built like this (C: Efficiently process (view, analize, clip ends, convert, demultiplex, dereplicate) would require way under 10MB or RAM (with GUI included in those 10MB).

If you still need this, I can build it for you.

Another really fast Python implementation which can deal also with very huge input lists of reads ids (e.g. 200 million read ids) which do not fit at once in the memory. It is using also frozenset for speed/memory reasons.

https://code.google.com/p/fusioncatcher/source/browse/extract_short_reads.py

What about AWK?

Sorry for necroposting, but I like puzzles, and had some free time on a vacation.

I've put together a short script that does the task at least for 4-line for even multiline FASTQ, using the simple idea of setting record separator to @ (i.e. RS = "@").

I want to find ~10,000 reads (by looking for their name) and pull out their sequence and quality

Also, I'm gonna assume you just want a regular FASTQ-formated output. The script is as follows:

#!/usr/bin/awk -f

## read IDs to extract
FNR == NR {
    ids[$0]++
    next
}

## parse FASTQ
FNR != NR {
    ## set record separator
    RS = "@"
    OFS = "\n"

    if ($1 in ids) {
        record = ""; sep = ""
        for (i = 1; i <= NF; i++) {record = record sep $i; sep = OFS}
        print "@"record
    }
}

You save it e.g. as fq.filter.awk and invoke as:

awk -f fq.filter.awk ids.txt input.fq > output.fq

It works with mawk and gawk, so the speed should be workable. Of course, bioawk is even faster, and more robust, as it works with multiline FASTQ (if you make it read the ids.txt)...