Downsampling a bam file for a specific number of reads

[kmac]

Is there an easy way to downsample a bam file for a specific number of reads? I know that picard and samtools can downsample to a proportion of reads (like 10% of your total reads) but I have three replicates for which I would like to extract an equal number of random reads.
Thanks in advance for any help!

[Brian Bushnell]

I have a tool, ReformatReads, that will extract an exact number of randomly distributed reads from a file. It's designed for fasta, fastq, and scarf, but it will process sam files (and bam files if samtools is installed) if the cigar strings are in 1.4 format (with = and X instead of M) and the reads are not paired. If the reads paired or in an old format, you would need to sample the fastq files first and remap them. It can convert sam/bam to fastq regardless of the cigar string version, but sam/bam have no guarantees about read order, so pairing information will be lost (so don't do that!)

reformat.sh in=reads.fq out=sampled.fq samplereads=1000

...would sample 1000 reads from a fastq file, or 1000 pairs if the file was interleaved.

reformat.sh in1=r1.fq in2=r2.fq out1=s1.fq out2=s2.fq samplereads=1000

...would sample 1000 pairs from paired read files.

reformat.sh in=mapped.bam out=sampled.bam samplereads=1000

...should sample 1000 reads from a bam file if all of the conditions are met (samtools installed, single-ended, sam 1.4 format reads).

You alternatively sample a specific number of bases worth of reads, if you have variable-length reads, with the "samplebases" flag.

[mgogol]

Here's a few ideas:

https://groups.google.com/forum/#!topic/bedtools-discuss/gf0KeAJN2Cw

Also look into bamtools random

[shenkers]

StreamSampler

I uploaded a command line utility I created to uniformly sample lines from a stream of input. Using samtools in a bash comand shell, you can use this command to uniformly downsample a bam file to a specific number of reads:

(
samtools view -H [bamfile];
samtools view -F 0x004 [bamfile] |
java -jar StreamSampler.jar [# of reads to sample] [total # reads]
) |
samtools -bS - > [sampled bam file]

It's important to keep in mind that this just does the downsampling, which as Brian mentions above, would result in a bam file with inconsistent flags if the data is paired. If this is important for your application, you would need another step to fix the flags on reads whose mates were discarded... Alternatively, you could avoid fixing flags by sampling the read names, and then filter out all reads that match one of the sampled names. If you want to try this out, you can download the jar from

https://github.com/shenkers/sampling/releases

[Chipper]

samtools sort -n bam namesorted
samtools view namesorted | head -n xxxxxx
or just calculate the percentage and use picard?