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

Have a look at bedtools coverage.

ADD REPLYlink written 10 weeks ago by ATpoint12k
3
gravatar for Pierre Lindenbaum
10 weeks ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum116k 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 10 weeks ago by Pierre Lindenbaum116k

thanks for sharing!

ADD REPLYlink written 10 weeks 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 10 weeks ago by lait120

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

ADD REPLYlink written 10 weeks ago by Pierre Lindenbaum116k

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

ADD REPLYlink modified 10 weeks ago • written 10 weeks 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 10 weeks 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 10 weeks ago by Pierre Lindenbaum116k

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 10 weeks ago • written 10 weeks ago by lait120

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

ADD REPLYlink written 10 weeks ago by Pierre Lindenbaum116k

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

ADD REPLYlink modified 10 weeks ago • written 10 weeks ago by Pierre Lindenbaum116k
2
gravatar for trausch
10 weeks 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 10 weeks ago • written 10 weeks 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 10 weeks ago • written 10 weeks 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 9 weeks 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: 1596 users visited in the last hour