It”s often a good idea (especially for large datasets) to use a subset of the reads for some preliminary analyses. A smaller dataset can be useful for figuring out the alignment parameters or checking for alignment problems.
The easiest way to get a subset of a fastq file is to take the reads in the beginning of the file. Of course, this is not a random subset of the reads, but it can still be useful and it”s very easy to do.
Just open a linux terminal and use the “head” command to get the first N lines of a file. Make sure to use a number that is divisible by four, otherwise you”ll end up with an incomplete record for one of the reads.
head -n 400 SRR515927.fastq > SRR515927_subset.fastq
Now, if you check the length of the new file (remember the word count [wc] command?), you can see that it contains 400 mobile casinos lines (i.e. 100 reads).
$ wc -l SRR515927_subset.fastq 400 SRR515927_subset.fastq
Remember, that this is often a bad representation of your dataset!
There are quite a few possible methods for getting a random subsample from fastq files. I will show you a few:
Subsampling using the HTSeq Python package (source SeqAnswers):
After installing HTSeq, you have to create a python script file containing the following lines:
import sys, random, itertools import HTSeq fraction = float( sys.argv ) in1 = iter( HTSeq.FastqReader( sys.argv ) ) in2 = iter( HTSeq.FastqReader( sys.argv ) ) out1 = open( sys.argv, "w" ) out2 = open( sys.argv, "w" ) for read1, read2 in itertools.izip( in1, in2 ): if random.random() < fraction: read1.write_to_fastq_file( out1 ) read2.write_to_fastq_file( out2 ) out1.close() out2.close()
To get approximately 10% of the reads, you can use the following command:
python htseq_subsample.py 0.1 SRR022913_1.fastq SRR022913_2.fastq SRR022913_1_htseq.fastq SRR022913_2_htseq.fastq
Let”s check the length of the original and new files!
$wc -l *_htseq.fastq 2479600 SRR022913_1_htseq.fastq 2479600 SRR022913_2_htseq.fastq
$wc -l reads/SRR022913_1.fastq 24818636 reads/SRR022913_1.fastq
You can also use the R BioConductor ShortRead package (source):
library(ShortRead) f1 <- FastqSampler("SRR022913_1.fastq", n=1000) f2 <- FastqSampler("SRR022913_2.fastq", n=1000) set.seed(123L); p1 <- yield(f1) set.seed(123L); p2 <- yield(f2) writeFastq(p1,"f1.fastq") writeFastq(p2,"f2.fastq")