Downsampling long-read BAM files

Question

[eebloom]

Im looking for the best way to downsample nanopore long-reads. Ideally I would want to downsample to a specific average read depth e.g. 5x

I have found a few methods e.g. using samtools but this seems to work for downsampling to a specific number or fraction of reads. Rasusa downsamples to a specific coverage /depth but requires fastq input and I would prefer to downsample my BAM files so I don't have to re-run my alignment on many samples again.

[GenoMax]

You can use reformat.sh from BBMap suite. But this may not satisfy the "specific avg depth".

reformat.sh -Xmx10g in=your.bam out=sampled.bam parameter_from_below

Sampling parameters:

reads=-1                Set to a positive number to only process this many INPUT reads (or pairs), then quit.
skipreads=-1            Skip (discard) this many INPUT reads before processing the rest.
samplerate=1            Randomly output only this fraction of reads; 1 means sampling is disabled.
sampleseed=-1           Set to a positive number to use that prng seed for sampling (allowing deterministic sampling).
samplereadstarget=0     (srt) Exact number of OUTPUT reads (or pairs) desired.
samplebasestarget=0     (sbt) Exact number of OUTPUT bases desired.

[eebloom]

Okay, thanks. Just double-checking this will work with long-reads?

[GenoMax]

Yes it should work with long read BAM's.

Answer 1

[eebloom]

random downsampling of bam files is now possible using rasusa v1.0.0 with the aln subcommand

[GenoMax]

Is this what you used ultimately? Did this tool satisfy the "specific average depth" requirement?

You could accept (green check mark) this answer to provide closure to this thread. If @Pierre's answer was also useful then it can also be accepted.

Answer 2

[Pierre Lindenbaum]

I'm not sure if this is what you want. I wrote Biostar154220

https://jvarkit.readthedocs.io/en/latest/Biostar154220/#biostar154220

(see Capping coverage in bam file Truncate unevenly covered BAM file at certain coverage )

[eebloom]

This is not what I needed for this particular use case, as capping the coverage would lose the information on regions of copy number amplification / the variation in coverage across the genome