Question: Extracting randomly subset of fastq reads from a huge file
0
gravatar for KVC_bioinfo
15 months ago by
KVC_bioinfo350
Boston
KVC_bioinfo350 wrote:

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 extract fastq • 1.9k views
ADD COMMENTlink modified 14 months ago by Bioinformatics_NewComer320 • written 15 months ago by KVC_bioinfo350
1

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.

ADD REPLYlink modified 15 months ago • written 15 months ago by Devon Ryan86k

long time: more than 2 hours

ADD REPLYlink written 15 months ago by KVC_bioinfo350

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

ADD REPLYlink modified 14 months ago • written 14 months ago by KVC_bioinfo350

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).

ADD REPLYlink modified 14 months ago • written 14 months ago by genomax58k

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

ADD REPLYlink written 14 months ago by KVC_bioinfo350

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).

ADD REPLYlink modified 14 months ago • written 14 months ago by Devon Ryan86k

Got it. thank you very much.

ADD REPLYlink written 14 months ago by KVC_bioinfo350

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

ADD REPLYlink written 14 months ago by KVC_bioinfo350
3
gravatar for genomax
15 months ago by
genomax58k
United States
genomax58k wrote:

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.
ADD COMMENTlink written 15 months ago by genomax58k
1

Ah, I got ninja'd :)

ADD REPLYlink written 15 months ago by Brian Bushnell16k

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

ADD REPLYlink written 14 months ago by KVC_bioinfo350
1

"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.

ADD REPLYlink written 14 months ago by Brian Bushnell16k

Thank you for the explanation.

ADD REPLYlink written 14 months ago by KVC_bioinfo350
1
gravatar for Brian Bushnell
15 months ago by
Walnut Creek, USA
Brian Bushnell16k wrote:

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).

ADD COMMENTlink written 15 months ago by Brian Bushnell16k
1
gravatar for Bioinformatics_NewComer
14 months ago by
Genomic Island
Bioinformatics_NewComer320 wrote:

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

ADD COMMENTlink written 14 months ago by Bioinformatics_NewComer320

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

ADD REPLYlink written 14 months ago by h.mon21k

It allows:

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

I'm using: vsearch v2.0.3

ADD REPLYlink modified 14 months ago • written 14 months ago by Bioinformatics_NewComer320
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1614 users visited in the last hour