Counting number of reads mapped to an allele

[jyu429]

Hi,

I have a sam file and I want to count how many queries mapped to a specific reference in the reference file. Is there a simple way to do this? Thanks!

[blancha]

One way would be with bedtools multicov.

You'll need to convert your SAM file to a BAM file, sort the BAM file and index it.
If you put your reference in the form of a BED file, you can then easily count the number of alignments for each interval defined in the BED file with bedtools multicov.

[jyu429]

So if I want to count the number of alignments to the reference with rname "IGHV01-02" how would I define this in the BED file?

[blancha]

Sorry, I misunderstood your question.
I thought reference meant a chromosomal position.
BED files work with chromosomal positions.
You could still answer your question with bedtools multicov though, if you defined the start as 0 and the end as a position larger than the largest position in IGHV01-02.
E. g.
IGHV01-02 0 10000000

However, if your question is that simple, you could just use grep.
grep -c IGHV01-02 yoursamfile

[jyu429]

Hi,

Yes, I did try grep but when I do grep -c "IGHV7-81*01" myfile.sam the count is 0 and when I grep for "IGHV7-81*" the count is 119, but I'm very sure there is at least one IGHV7-81*01.

[blancha]

I don't know.

If your SAM file has a header, you could check the reference names in the header.
You could also check the reference names in the reference genome, as well as test the grep syntax on the reference genome.

You could also convert the SAM file to a BAM file, sort it and index it.
You could then load the BAM file and the reference genome in IGV to visualize the alignments.