find coverage of specific exomes from a bam file
2
0
Entering edit mode
3.8 years ago
lait ▴ 160

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.

coverage bedtools genomecov bam bed • 1.9k views
ADD COMMENT
1
Entering edit mode

Have a look at bedtools coverage.

ADD REPLY
3
Entering edit mode
3.8 years ago

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
ADD COMMENT
0
Entering edit mode

thanks for sharing!

ADD REPLY
0
Entering edit mode

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 .......
ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
2
Entering edit mode
3.8 years ago
trausch ★ 1.8k

Alfred can be used for that:

alfred count_dna -i exome.bed.gz -o coverage.gz input.bam
ADD COMMENT
0
Entering edit mode

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
ADD REPLY
0
Entering edit mode

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).

ADD REPLY

Login before adding your answer.

Traffic: 1842 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6