Question: find coverage of specific exomes from a bam file
0
gravatar for lait
9 days ago by
lait120
lait120 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 9 days ago by Pierre Lindenbaum114k • written 9 days ago by lait120
1

Have a look at bedtools coverage.

ADD REPLYlink written 9 days ago by ATpoint9.3k
3
gravatar for Pierre Lindenbaum
9 days ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum114k 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 9 days ago by Pierre Lindenbaum114k

thanks for sharing!

ADD REPLYlink written 8 days ago by lait120

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 8 days ago by lait120

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

ADD REPLYlink written 8 days ago by Pierre Lindenbaum114k

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

ADD REPLYlink modified 8 days ago • written 8 days ago by lait120

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 8 days ago by lait120

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 8 days ago by Pierre Lindenbaum114k

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 8 days ago • written 8 days ago by lait120

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

ADD REPLYlink written 8 days ago by Pierre Lindenbaum114k

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

ADD REPLYlink modified 8 days ago • written 8 days ago by Pierre Lindenbaum114k
2
gravatar for trausch
9 days ago by
trausch1.1k
Germany
trausch1.1k wrote:

Alfred can be used for that:

alfred count_dna -i exome.bed.gz -o coverage.gz input.bam
ADD COMMENTlink modified 9 days ago • written 9 days ago by trausch1.1k

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 8 days ago • written 8 days ago by lait120

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 5 days ago by trausch1.1k
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: 770 users visited in the last hour