Hello, what command could I use to extract every 3rd or every 4th pair of reads from two fastq files corresponding to pair-end reads (file read1 and file read2)? I want to use it for getting smaller raw data files, which would be representative of the larger files, since they would be reduced randomly.
Also, the same question regarding how to randomly reduce bam and bed files, such as accepted_hits.bam and junctions.bed generated by tophat?
Thank you,
Ephraim Trakhtenberg
•
written
by
trakhtenberg
2 answers
To uniformly sample a third of your paired reads (without replacement) from a paired fastq file(*), via bash and sample:
$ fq_count=`calc "$(wc -l < all_reads.fastq) / 8"`
$ samples=`calc "floor($fq_count / 3)"`
$ sample --lines-per-offset=8 --sample-size=${samples} all_reads.fastq > random_sample.fastq
(* - A paired fastq file would be one interleaved file that stores paired reads consecutively, one pair after its mate.)
Similarly, to sample a BED file:
$ sample --sample-size=12345 foo.bed > random_sample.bed
Or if you have sufficient memory:
$ shuf foo.bed | head -n12345 > random_sample.bed
The calc tool is not a builtin, but I am using bash notation for shell commands. Source is available here, though it is probably available in convenient package form with your favorite package manager.
A minor point, recent versions of shuf will do memory-efficient reservoir sampling with the -n argument, so you don't need to use head.
I think expr $fq_count / 3 will give you the floor by default, so you could replace all "calc" with "expr" in this case.
•
written
by
Matt Shirley
I believe a truly random approach is already implemented in samtools for bam file extraction.
The option -s enables you to subsample a bam file.
-s FLOAT Fraction of templates/pairs to subsample; the inte‐
ger part is treated as the seed for the random num‐
ber generator [-1]
samtools view -s 0.25 -b mymapping.bam > random_25%_of_mymapping.bam
I had 167115 reads in the input BAM, so 25% is 41779: using 0.25 gets 42134; using 23241.25 gets 41567 (which is closer to 41779 but still 212 off). I would appreciate a clarification on how could I produce output closer to the specified percent, and a brief explanation on why replacing 0 with 23241 made the output closer to 25%, and how could it increase randomness? Also, just to confirm, in both approaches if reads are paired, they are randomly reduced also in pairs (just to make sure the reads are not selected randomly regardless of their pairing). Thank you.
Log in to answer this question.
-
Tophat output accepted_hits.bam is missing
written by Biologist
-
【tophat】 RNA-seq mapping reference genome
written by bio_zhangxl •
-
need to separate read1 and read2 to run tophat2
written by Ann
-
Picard error Illegal Mate State in converting BAM to Fastq
written by trakhtenberg
-
Extracting from tophat outputs reads pairs and splice-junctions with a single best match
written by trakhtenberg
-
tophat output files containing the reads which mapped uniquely as a pair
written by trakhtenberg
-
tophat2 reported the same alignment twice for the same read (paired end data)
written by Ann
-
QC based on Samtools Flagstat
written by indranipoddar •
-
Questions About Pair-End Reads Which Are Not Match For Bwa
written by Tonyzeng •
-
Junctions.Bed File Produced From Tophat
written by Varun Gupta
The approach you mention is not technically random, but this task is a pretty common one. There has been plenty of good discussion of this site about sampling large files, for example, Selecting Random Pairs From Fastq? (and see the 'similar posts' section to the right of that page for more posts on the topic). Hope that helps.