Finding High Coverage Regions In A Bam File?

Question

[2184687-1231-83-]

I would like to have a list of regions that are in the top X% coverage in a bam file, centered at the peak of coverage for a window of size W. For example, something like:

~/find_top_X_perc -perc 10 -window 250 file.bam

That will produce a list of regions ready to be used on samtools tview:

1:1000000-1000250
2:2000250-2000300
[...]

Any suggestions?

Answer 1

[Jeremy Leipzig]

i would go with R & Bioconductor:

library("Rsamtools")
library("rtracklayer")
rbga<-readBamGappedAlignments("myBamfile.bam")
bamCov<-coverage(rbga)
#quantile returns a list of the ith % for each chromosome, you could weight the mean properly if you are ambitious
cutoff<-mean(quantile(bamCov,.9))
hotslices<-slice(bamCov,lower=cutoff)
hotwindows<-resize(hotslices,width=250,fix="center")
dummy<-file("mybed.bed")
export.bed(hotwindows,dummy)

[Istvan Albert]

This is neat - declarative rather than imperative programming - I need to learn the NGS packages for R

[Jeremy Leipzig]

you could actually combine lines 3-9 into one line but it would be unreadable

Answer 2

[Ketil]

As it happens, I recently wrote a small program to output coverage in user-defined buckets from a BAM file. Currently you need to get it via darcs

darcs get http://malde.org/~ketil/biohaskell/bamcov

And compile it yourself (e.g. per instructions here: http://biohaskell.org/Installation)

The program outputs five columns per bucket: chromosome, position, and paired, orphan, and link coverage (the latter being paired reads mapping to different chromosomes).

[2184687-1231-83-]

I tried to install bamcov, but failed. As a reference, here are the installation steps: [?] sudo apt-get install darcs sudo apt-get install cabal-install darcs get http://malde.org/~ketil/biohaskell/bamcov cabal update cabal install [...] cabal: Error: some packages failed to install: bamcov-0 depends on seqloc-0.2.1 which failed to install. samtools-0.1.2 depends on seqloc-0.2.1 which failed to install. seqloc-0.2.1 failed during the configure step. The exception was: ExitFailure 1 [?]

[Ketil]

Basically, that shoulda worked. You'd need to check out why seqloc failed, it should be in the output somewhere. (The seqloc library is used by the samtools bindings, so this is an unavoidable requirement, I'm afraid.)

Answer 3

[Zhidkov]

You can start with turning your bam file to pileup, perform summary statistics on coverage column (to establish what is high coverage in current bam), then extract coordinates position with high coverage, then again run the samtools this time extract reads which correspond to high coverage regions.

Ilia