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.
2 answers
use the option -s of samtools view
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)
ok, so you can try my tool: http://lindenb.github.io/jvarkit/Biostar154220.html
see also Capping coverage in bam file
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.
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...
Log in to answer this question.