Counting entries in bam-files: difficulties with samtools and bamToBed

[Azazel]

I have a bam-file 'stuff.bam' and want to know how many sequences are described in it. I am not getting the same results when using samtools and bamToBed:

[prompt]$ ~/samtools-0.1.12a/samtools view -c stuff.bam 
21054492

...ok, now I count the entries in a conversion of stuff in bed-format:

[prompt]$ bamToBed -i stuff.bam | wc -l
20567680

Why does this give different results, and how can I get the correct number?

--
I am using
samtools-0.1.12a
BEDTools-Version-2.9.0

[droog_22]

BAM files can contain mapped and unmapped reads

Dear Azazel,

Using samtools view you are counting the complete number of reads in your BAM file - be it mapped or unmapped reads.

$ samtools view data.bam | wc -l

Information about which read was mapped and how can be seen in the flag field of the BAM file. See more here http://samtools.sourceforge.net/samtools.shtml and here http://picard.sourceforge.net/explain-flags.html

You could use this to count the different flag tags in your file:

$ samtools view data.bam | awk '{print $2}' | sort | uniq -c
10976575 0
10968005 16
5483180 4

So 5.5M reads do not map the genome at all, and roughly 11M map to the forward respectively the reverse strand of the genome you used to map the reads.

Using BamToBed you just count the number of mapped reads. A which hasn't been mapped obviously does not have any chromosome, start, and end information and is thus not taken into account.

So the following command will give you the total 22M reads which map to the genome.

$ bamToBed -i data/bowtie_remapped/input_0.5-2h_1_o.bam | wc -l
21944580

Hope this helps,

Cheers, droog_22

[frozenlyse]

"samtools flagstat" is probably a better tool for the job, it will give a summary of all the reads in your bam file eg

$ samtools flagstat data_1.bam 
25038903 in total
0 QC failure
0 duplicates
16921771 mapped (67.58%)
0 paired in sequencing
0 read1
0 read2
0 properly paired (-nan%)
0 with itself and mate mapped
0 singletons (-nan%)
0 with mate mapped to a different chr
0 with mate mapped to a different chr (mapQ>=5)

[DerSeb]

What I like to di is to count the flags that are present:

samtools view *.bam | awk ' { per[$2] += 1 }END { for (i in per)print i, per[i] }

This command will show you which flags are present and how often they occur.

[Azazel]

OK! I understand. Just wanted to say thanks for the answers.