Report the number of alignments for every read

Question

[sorrymouse]

I have a standard bam file produced by bowtie using -m 50. I'd like to have a report of the number of possible alignments for every read in the bam file - for example if the read maps to 5 locations or 20. I have specified that one mismatch is allowed, so this should be a discrete number. Googling this problem leads me to so many posts about just counting the number of reads in a bam file that I cannot figure out how to effectively google this problem. Any suggestions would be welcome.

Answer 1

[ATpoint]

I would say it comes down to counting the number of read names. In this below example I made a dummy alignment where one read gets aligned to three different chromosomes:

  read1 0   chr2    1   255 48M *   0   0   TTAGCTACGTACGATCGATCATCGACTGATCATCGGATACTAGCTACG    TTAGCTACGTACGATCGATCATCGACTGATCATCGGATACTAGCTACG    XA:i:0  MD:Z:48 NM:i:0
  read1 0   chr3    1   255 48M *   0   0   TTAGCTACGTACGATCGATCATCGACTGATCATCGGATACTAGCTACG    TTAGCTACGTACGATCGATCATCGACTGATCATCGGATACTAGCTACG    XA:i:0  MD:Z:48 NM:i:0
  read1 0   chr1    1   255 48M *   0   0   TTAGCTACGTACGATCGATCATCGACTGATCATCGGATACTAGCTACG    TTAGCTACGTACGATCGATCATCGACTGATCATCGGATACTAGCTACG    XA:i:0  MD:Z:48 NM:i:0

So you would simply extract the first column that contains the read name and then count its occurrence. if you have a BAM file you can do:

samtools view foo.bam | cut -f1 | sort -k1,1 | uniq -c | sed 's/^ *//g' | tr " " "\t"

In this case that would be:

3   read1

As bowtie only does ungapped alignments I think this very simple approach would be fine, there is no need to handle split reads etc because these (in my understanding) all are not supported by bowtie and the read would go unmapped.

[sorrymouse]

Thank you, that worked!