genomeCoverageBed out > mapping out

Question

[Siva]

Hello guys,

An other doubt today ;)

Got a bwa sam file results from a mapping of single reads, I have in my sam file 33 results, means 33 of my reads into 33 reads have been mapped, no overlap between reads and just uniquely mapped reads.

So I expected to have 33 lines in my output file from "genomeCoverageBed" but it's not I got 42 lines ! O_o

$ samtools view -b mapped_bowtie_ref-33CDS_local.sorted.bam | genomeCoverageBed -ibam stdin -bga -split -g chr_size_HA412.tab | awk '$4 > 0' > cov_map_bowtl.txt

I've done a -bga -split http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html so should have only 33, and moreover no overlap so I have the value "1" in the 4th column that's good, but why 42 lines instead of 33?

Thanks,
Siva

Answer 1

[James Ashmore]

You have specified -bga flag and that will report regions between your reads which have zero coverage. That might be why you have more BED records than BAM alignments. Also, there is no need to pipe samtools output into bedtools, just specify your bam file name after the -ibam flag.