20
$\begingroup$

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?

$\endgroup$
1
  • $\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, 2017 at 2:03

4 Answers 4

17
$\begingroup$

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.

$\endgroup$
2
  • 10
    $\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, 2017 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, 2017 at 11:47
4
$\begingroup$

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 ]]
  then 
  echo '[ERROR]: Requested number of reads exceeds total read count in' $1 '-- exiting' && exit 1
fi

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

It have used both sambamba and samtools. Either one works great to downsample to a desired read depth after calculating the downsample factor.

Sambamba can use multiple threads by using -t INT, where INT is the number of threads. It can also set the compression level for the BAM file using -l.

BAM="Sample.bam"
READS="10000000"
fraction=$(samtools idxstats $BAM | cut -f3 | awk -v ct=$myReads 'BEGIN {total=0} {total += $1} END {print ct/total}')

# With samtools
samtools view -b -s  ${fraction} ${BAM} > output.bam

# With sambamba
sambamba view -f bam -l 5 -s  ${fraction} ${BAM} > output.bam
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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