Question: find coverage of specific exomes from a bam file
0
gravatar for lait
4 months ago by
lait130
lait130 wrote:

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.

ADD COMMENTlink modified 4 months ago by Pierre Lindenbaum118k • written 4 months ago by lait130
1

Have a look at bedtools coverage.

ADD REPLYlink written 4 months ago by ATpoint14k
3
gravatar for Pierre Lindenbaum
4 months ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum118k wrote:

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 COMMENTlink written 4 months ago by Pierre Lindenbaum118k

thanks for sharing!

ADD REPLYlink written 4 months ago by lait130

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 REPLYlink written 4 months ago by lait130

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

ADD REPLYlink written 4 months ago by Pierre Lindenbaum118k

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

ADD REPLYlink modified 4 months ago • written 4 months ago by lait130

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 REPLYlink written 4 months ago by lait130

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 REPLYlink written 4 months ago by Pierre Lindenbaum118k

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 REPLYlink modified 4 months ago • written 4 months ago by lait130

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

ADD REPLYlink written 4 months ago by Pierre Lindenbaum118k

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

ADD REPLYlink modified 4 months ago • written 4 months ago by Pierre Lindenbaum118k
2
gravatar for trausch
4 months ago by
trausch1.2k
Germany
trausch1.2k wrote:

Alfred can be used for that:

alfred count_dna -i exome.bed.gz -o coverage.gz input.bam
ADD COMMENTlink modified 4 months ago • written 4 months ago by trausch1.2k

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 REPLYlink modified 4 months ago • written 4 months ago by lait130

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 REPLYlink written 4 months ago by trausch1.2k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 971 users visited in the last hour