exon coverage by coverageBed

[frankyue50]

Hello everyone,

So I just want to compute the average coverage per base for our exome-seq. First I collected all the exons in a bed file, say exon.bed and it looks like this:

chr1 11873 12227
chr1 12594 12721

I also have an exam seq bam file, and what I did is:

samtools view -b bamfile.bam | coverageBed -abam stdin -b exon.bed -d > report, which looks like this:
chr1 161480623 161480746 1 0
chr1 161480623 161480746 2 0
chr1 161480623 161480746 3 0
chr1 161480623 161480746 4 0
chr1 161480623 161480746 5 0
chr1 161480623 161480746 6 0

But the thing is the last column is always zero. What did I do wrong? Thank you guys for your help!

[GenoMax]

Your chromosome names seem to match between the bed file and bam file so that can't be the problem. That was not the case.

Did you create/edit the bed file on a PC/Mac and then move it to unix server? I wonder if the file format differences is causing the problem you are seeing.

[frankyue50]

this is from samtools view reads.bam:
name 163 1 10001 29 39S25M37S = 10026 126 ...
Looks like a normal bam file to me

Your chromosome names seem to match between the bed file and bam file so that can't be the problem.

Did you create/edit the bed file on a PC/Mac and then move it to unix server? I wonder if the file format differences is causing the problem you are seeing.

[GenoMax]

this is from samtools view reads.bam:
name 163 1 10001 29 39S25M37S = 10026 126 ...
Looks like a normal bam file to me

Your reference file seems to have chromosome names as "1" as opposed to chr1. If all the chromosomes do not have "chr" in front then you could edit your bed file accordingly.

[frankyue50]

I guess that's the reason... I print out the bed file from bam and things worked out.

Thanks, GenoMax!

Your reference file seems to have chromosome names as "1" as opposed to chr1. If all the chromosomes do not have "chr" in front then you could edit your bed file accordingly.