Sampling random reads from BAM files with allowing same reads

Question

[SSK]

Hi.

I'm trying to randomly sample reads with allowing same reads.

I tried to random sampling by using "samtools view -s", but it seems like samtools samples without duplicate.

How can I do this ?

Thanks!

[WouterDeCoster]

I don't think that will give you a valid bam then if you reuse read names.

[jkbonfield]

The aim of samtools view -s is to simulate sequencing at a reduced depth.

It sub-samples based on a hash of the read name. Thus all reads in a pair are kept or none of the reads in that pair. This is better than a naive sampling which ends up breaking most read pairs.

I don't understand what you're asking for beyond how it already works.

Answer 1

[benformatics]

Please use the search function:

Question: Sampling with replacement from a bam file

However, here is a solution in R:

library(GenomicAlignments)
library(rtracklayer)

## your bam file
bam.file <- 'your_bam_file.bam'

## is it paired end or not
paired <- TRUE

## read in the bam file
if(paired){
  bam.input <- readGAlignmentPairs(bam.file)
}else{
  bam.input <- readGAlignments(bam.file)
}

## select your number of subsets
num.samples <- 100

## subsample WITH replacement
bam.subsample <- sample(bam.input,num.samples,replace=TRUE)

## export new bam file
export(bam.subsample,BamFile('bam_subsampled.bam'))

[benformatics]

You lose a lot of metadata this way so if you want anything specific you should look at using the ScanBamParam() parameter during the 'read in the bam file' step.

[SSK]

Thanks benformatics! I'ts OK if I lose a lot of meta data. I try using your script. Thanks!

[benformatics]

Fixed a typo in case you are using unpaired reads.