I know how to downsample a BAM file to lower coverage. I know I can randomly select lines in SAM, but this procedure can't guarantee two reads in a pair are always sampled the same time. Is there a way to downsample BAM while keeping pairing information intact?

  • $\begingroup$ Could you please give a bit more explanation as to why it would be useful to downsample reads? There are a few different ways to do it, and the answers will depend on your goal. For example, you can do a digital normalisation to reduce the input read count for genome assembly, but that method would not be appropriate for a quantitative analysis of reads (e.g. for RNASeq). $\endgroup$ – gringer Jun 3 '17 at 2:03

samtools has a subsampling option:

-s FLOAT: Integer part is used to seed the random number generator [0]. Part after the decimal point sets the fraction of templates/pairs to subsample [no subsampling]

samtools view -bs 42.1 in.bam > subsampled.bam

will subsample 10 percent mapped reads with 42 as the seed for the random number generator.

| improve this answer | |
  • 9
    $\begingroup$ To avoid possible confusion in the future, it should be noted that this works by hashing the read name. If for some terrible reason the mates don't have the exact same name (e.g., /1 and /2 weren't stripped) then this will break. $\endgroup$ – Devon Ryan Jun 3 '17 at 11:13
  • 1
    $\begingroup$ And also due to the way this works, the number of subsampled entries is not guaranteed to be the exact required percentage: some transformation of the hash residing in the [0,1] interval is compared against the sampling proportion. Depending on the read names present in the file, the number of sam entries effectively sampled will vary. $\endgroup$ – bli Jun 3 '17 at 11:47

With this function, you can subsample a given BAM file to a user-specified number of reads using SamBamba. The function automatically calculates the scaling factor. Ideas in part taken from here

function SubSample {

## Calculate the sampling factor based on the intended number of reads:
FACTOR=$(samtools idxstats $1 | cut -f3 | awk -v COUNT=$2 'BEGIN {total=0} {total += $1} END {print COUNT/total}')

if [[ $FACTOR > 1 ]]
  echo '[ERROR]: Requested number of reads exceeds total read count in' $1 '-- exiting' && exit 1

sambamba view -s $FACTOR -f bam -l 5 $1


## Usage example, selecting 100.000 reads:
SubSample in.bam 100000 > subsampled.bam

Note that $2 is the number of reads. Say you choose 10000 reads but have paired-end data, you'll end up with 5000 fragments.

| improve this answer | |

A totally different tactic...it helps if you aren't super picky about how many reads you have at the end.

Illumina read names contain the lane, tile, and xy coordinates of the read in them. If you filtered by these, you would be guaranteed to get either both members of a pair, or neither. I would recommend trying to avoid the edges of the flowcell.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.