coverageBed counts not the same as input bed

[cmccabe]

When I use a sorted bed file with unique 56,548 targets in it and process it in bedtools 2.25.0 with the below command:

coverageBed -d -g /home/cmccabe/Desktop/NGS/bedtools2-25.0/genomes/human.hg19.Ion.genome.txt -a /home/cmccabe/Desktop/NGS/bed/bedtools/bedtools_unix_sorted_xgen_targets.bed -b /home/cmccabe/Desktop/NGS/IonXpress_004.R_2015_10_29_13_11_15_user_Proton-32-Lurie_Inh_Disease_151029_newheader.bam > /home/cmccabe/Desktop/NGS/IonXpress_004.R_2015_10_29_13_11_15_user_Proton-32-Lurie_Inh_Disease_151029_newheader_basecounts.bed

I get 56,543 unique targets returned. The sort order is lexical in the bed and human in the bam so I use the -g option to specific the ordering of the bam. The 56,543 was confirmed in perl and awk. Any ideas why the two #'s are not the same? Thank you .

[GenoMax]

Perhaps those 5 targets are not covered by reads in your BAM? Did you look at the BAM in IGV?

[cmccabe]

I will try that, but shouldn't bedtools output that region with all zero's.

There are other regions in the coverageBed output where the regions is zero.

For example each one of my targets is 120 bases and the coverage per base is calculated. So if there was no coverage in those 120 bases, there is no output in the file? Thank you .

[GenoMax]

I see.

Just for the heck of it can you do a "sort/uniq" on the names to see if those 5 are for some reason repeated in the file?

Edit: Just thought that the intervals may be repeated (instead of just names) as another possibility.

[cmccabe]

If I do

wc -l

on the original bed, I get 56,548, but if I do:

awk '$4!=d{c++;d=$4}END{print c}' file

or perl -nae '$seen{$F[3]}++;
END{
print "There are ", scalar keys %seen, " unique fourth fields\n";
}' originalbed [/CODE]

I get 56,543

However, in excel there are 56,548. Since the file was sent to me I think I am just going to re-do it without using excel (which is very problematic). Thank you .

[GenoMax]

5 blank lines at the end of that file?

[cmccabe]

I have put the bed file on box.net as it is too large to attach:

https://app.box.com/s/8fuk40t7fn77jikymullug3ueeenk2j3

When I do a wc-l test.bed I get 56,548

However, when I use awk or perl I get 56,543 and that is also what bedtools outputs.

cmccabe@DTV-A5211QLM:~/Desktop/NGS/bed/bedtools$ wc -l test.bed
56548 test.bed
cmccabe@DTV-A5211QLM:~/Desktop/NGS/bed/bedtools$ awk '$4!=d{c++;d=$4}END{print c}' test.bed
56543
cmccabe@DTV-A5211QLM:~/Desktop/NGS/bed/bedtools$ perl -nae '$seen{$F[3]}++;
>     END{
>         print "There are ", scalar keys %seen, " unique fourth fields\n";
>     }' test.bed
There are 56543 unique fourth fields

I can not seem to figure out why there are 5 missing. Thank you .

edit:
you were right @Genomax they were non-unique and excel didn't mark them:

4748d4747
< chr11	47270217	47270425	chr11:47270217-47270425	unknown-1062|gc=64.9
4970d4968
< chr11	5248271	5248449	chr11:5248271-5248449	HBB-283|gc=55.1
24883d24880
< chr19	13010118	13010237	chr19:13010118-13010237	SYCE2-864|gc=47.9
33027d33023
< chr22	38153605	38154160	chr22:38153605-38154160	TRIOBP-610|gc=68.6
54957d54952
< chrX	33357316	33359011	chrX:33357316-33359011	DMD-581|gc=33.7

but fgrep did

$ fgrep 'chr11:47270217-47270425' test.bed
chr11	47270217	47270425	chr11:47270217-47270425	ACP2-1062|gc=64.9
chr11	47270217	47270425	chr11:47270217-47270425	unknown-1062|gc=64.9

Thank you

[GenoMax]

Ah. It is good to have a logical explanation.