find coverage of specific exomes from a bam file

Question

[lait]

Hi

I know that

bedtools genomecov -ibam xx.bam  -bg

would report genome coverage in BEDGRAPH format

What if I need the coverage not for the whole genome, but just for specific regions (specific exomes) ? how can I achieve this ? my input files are a bam file and a .bed file containing the start and end positions of exomes.

[ATpoint]

Have a look at bedtools coverage.

Answer 1

[Pierre Lindenbaum]

I wrote : http://lindenb.github.io/jvarkit/BamStats04.html

$ java -jar dist/bamstats04.jar -B src/test/resources/toy.bed.gz src/test/resources/toy.bam 2> /dev/null | column -t 

#chrom  start  end  length  sample  mincov  maxcov  meancov  mediancov  nocoveragebp  percentcovered
ref     10     13   3       S1      3       3       3.0      3.0        0             100
ref2    1      2    1       S1      2       2       2.0      2.0        0             100
ref2    13     14   1       S1      6       6       6.0      6.0        0             100
ref2    16     17   1       S1      6       6       6.0      6.0        0             100

[lait]

thanks for sharing!

[lait]

Hi Pierre, why do I always get the error: Unknown contig ? which refers to the first line in my bed file ? I tried to play around with the bed file , use different coordinates, but I always get this error.

My bed file is something similar to:

ch1 115256528   115256530
ch2 .......

[Pierre Lindenbaum]

and i suppose your bam file use '1' and '2'... ?

[lait]

when replacing chx with x , I get the message "ignoring" for all the entries in the bed file

[lait]

just figured out the error, the bed file contains locations of length = 1 (a single nucleotide). it should be of length 2 or more , because of the condition present in your java file.

[Pierre Lindenbaum]

hum.. this ignoring is displayed when start > end. https://github.com/lindenb/jvarkit/blob/master/src/main/java/com/github/lindenb/jvarkit/tools/bamstats04/BamStats04.java#L274

what is the output of

file your.bed

and

grep X -m1 your.bed |tr "\t" "#"

please

[lait]

ASCII text

and the second command returned nothing

So in my case, the start and end coordinates are equal, I dont know why the condition turned to be TRUE

[Pierre Lindenbaum]

ha, it's x in lowercase grep x -m1 your.bed |tr "\t" "#"

[Pierre Lindenbaum]

otherwise send a bug report with a minimal bed and a minimal bam please https://github.com/lindenb/jvarkit/issues/

Answer 2

[trausch]

Alfred can be used for that:

alfred count_dna -i exome.bed.gz -o coverage.gz input.bam

[lait]

thanks! do you have an idea how to solve this error when installing Alfred:

sudo conda install -c bioconda alfred  
Fetching package
metadata .............


 PackageNotFoundError: Package not found: '' Package missing in current
 osx-64 channels: 
   - alfred

[trausch]

Osx builds are skipped for performance reasons in the Bioconda versions, linux-64 is the default platform. You can either build from source on osx (requires boost library) or use this minimal docker container of Alfred (if you are familiar with docker).