This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extracting randomly subset of fastq reads from a huge file

Hello,

I am trying to extract a random subset of paired end fastq reads from a huge 5Gb file. I tried using seqtk and zcat. However, it is taking long time.

The command I used:

zcat huge_file.fastq.gz | sort -R | head -1000 > reads.sample_random.fq

Could someone please help me? Thank you in advance.

random fastq extract

Ummm, won't sort -R completely screw everything up (and be slow and use a LOT of memory)? Define "a long time", you're mostly limited by IO and gzip.

long time: more than 2 hours

I found these commands. However, what is -s (seed)? And is it sampling randomly?

Thank you in advance.

Apply a seed to extract the same reads from two, paired end files:

seqtk sample -s 10 /ebs/ecoli/SRR001666_1.fastq.gz 1000 > SRR001666_1_1000.fastq

seqtk sample -s 10 /ebs/ecoli/SRR001666_2.fastq.gz 1000 > SRR001666_2_1000.fastq

In this particular case the same seed is needed to keep the pairing of reads (i.e. extract the same read pair from two files).

Thank you. But what does the number after s represent?

It's a random seed. It doesn't matter what number you use, you just have to use the same one (I usually use 1234 for things like this).

Got it. thank you very much.

Worked perfectly for me. Thank you, everyone, for responses.

4 answers

reformat.sh from BBMap suite.

reformat.sh in1=file_R1.fq.gz in2=file_R2.fq.gz out1=Sampled_R1.fq.gz out2=Sampled_R2.fq.gz parameters_below

Sampling parameters:

reads=-1                Set to a positive number to only process this many INPUT reads (or pairs), then quit.
skipreads=-1            Skip (discard) this many INPUT reads before processing the rest.
samplerate=1            Randomly output only this fraction of reads; 1 means sampling is disabled.
sampleseed=-1           Set to a positive number to use that prng seed for sampling (allowing deterministic sampling).
samplereadstarget=0     (srt) Exact number of OUTPUT reads (or pairs) desired.
samplebasestarget=0     (sbt) Exact number of OUTPUT bases desired.
                        Important: srt/sbt flags should not be used with stdin, samplerate, qtrim, minlength, or minavgquality.

So shall I use reads=100 (no of reads I want to extract?) parameter?

"reads=100" means it will sample from 100 reads. So if you set"reads=100 samplerate=0.5" you'll get approximately 50 reads, sampled from the first 100 reads in the file. Whereas if you set "samplereadstarget=100" you will get exactly 100 reads, sampled from the full file.

Note that if your reads are paired, these will be the number of PAIRS you get, so the number of reads would be twice that.

Thank you for the explanation.

I assume huge_file.fastq.gz is interleaved?

Using BBMap, you can do this:

reformat.sh in=huge_file.fastq.gz out=subsampled.fastq.gz samplerate=0.1 interleaved

With 2 files you can use in1= and in2= instead. There's also "samplereadstarget=X" if you want a specific number of reads (for paired data, it will give you that number of pairs).

See if you'd like to try this with VSEARCH.

I think VSEARCH subsampling doesn't work with fastq files, only with fasta.

It allows:

    Subsampling
  --fastx_subsample FILENAME  subsample sequences from given FASTA/FASTQ file

I'm using: vsearch v2.0.3

You could use sample:

$ sample -k 1000 -l 4 huge_file.fastq > 1k_sample.fq

This uses reservoir sampling on newline offsets to reduce memory overhead. If your system is thrashing for lack of system memory, this will help a great deal by reducing the memory overhead required for sampling with other tools (including GNU sort).

Specify -k N for the number N of samples, and -l L for number of lines per record-to-sample (L). In the case of fastq, that would be four lines per record, so -l 4.

Add -d S to specify a seed value S, if desired, or it is drawn from a Mersenne Twister PRNG.

Generating your own random seed and applying to sampling two fastq files would let you grab the same sample offsets from two files of paired reads. If you have one interleaved file, then you could specify -l 8 to sample a pair of records from every eight lines.

You can also add -r to sample with replacement; the default is to sample without replacement. Or add -p to print a sample in the order provided by the original input.

Does that work on compressed files?

The tradeoff for sampling on line endings (and having a much lower memory hit) is the need to seek through a regular file, not a file stream, so an uncompressed regular file is required. I haven't tested a named pipe, but I don't think that would work because a pass through a stream discards past data buffers.

Thank you. I have compressed files.

Log in to answer this question.