Truncate unevenly covered BAM file at certain coverage

Question

[predeus]

Hello all,

This is basic enough question, but I can't seem to find the answer online.

I have a BAM from RNA-seq alignment which is very unevenly covered. However the mean coverage is huge (~100,000x - it's a mitochondrion). I want to use it for polishing, so I would like to retain about 100x "sliding window" (or simply at any given nucleotide) coverage across the whole chromosome, and lose the rest.

How can I achieve this?

Thank you in advance, as always

UPD: Apparently, there's a discussion here that directly answers my question.

Answer 1

[Pierre Lindenbaum]

use the option -s of samtools view

[predeus]

Nice, my current version of samtools (1.7) doesn't even have that yet! Thank you very much.

[predeus]

All right, I've read about the option and still can't figure out how to use the "filter expressions" for coverage-based filtering. I might be missing something obvious...

[Pierre Lindenbaum]

my bad. I meant option -s not option -e. (fixed)

  -s FLOAT subsample reads (given INT.FRAC option value, 0.FRAC is the
           fraction of templates/read pairs to keep; INT part sets seed)

[predeus]

I see. I knew about -s but this is not what I'm after - I don't want to proportionally downsample the reads, I only want to remove reads in the areas where coverage is > 100, for example.

[Pierre Lindenbaum]

ok, so you can try my tool: http://lindenb.github.io/jvarkit/Biostar154220.html

see also Capping coverage in bam file

[predeus]

Brilliant, I didn't find that thread somehow (it's tricky to select the right keywords). Thank you!

Brian's explanation about normalization via k-mers also makes a lot of sense - it has long puzzled me what Trinity was doing for example.

Answer 2

[QLFblaireau]

You could use samtools depth

for example

samtools depth -a my.bam |awk '$3 >=100' > file.coverage

Then make a bedfile out of the file.coverage and use it in combination with samtools view to only extract position where there is the expected coverage.

Or easier, you could use mosdepth that will directly output a bed file that you can filter.

mosdepth

[predeus]

Thank you! However, I don't want to find regions with some type of coverage - I want to remove the excessive coverage, retaining reads up to a certain maximum, e.g. 100.

This approach also can work (together with samtools view -s), but it becomes dreadfully convoluted - you need to

• split bam into windows;
• downsample each file according to mean coverage in each window;
• put them back together.

Seems like there should be a tool that's doing it...