How to call "gaps" from spliced read alignment?

[sulicon]

Hi all,

I have aligned my long 454 reads to the reference genome by using GMAP, a spliced read aligner. Now I want to extract the best set of intronic regions (gaps) in the alignments. Is there any tool could generate such information from a SAM file?

I feel this task is somewhat similar with indel calling, if you considered intron regions as huge deletions. But I'm not sure whether I can use indel calling tools to fit my need.

A simple (and naive) way I can image is counting the coverage of the gaps. The chromosome regions with many gaps mapped should be reliable intronic regions. For this purpose, is there a convenient way to calculate the coverage of gaps from a SAM file?

Thanks a lot,
Shuli

[frozenlyse]

Are you familiar with R/Bioconductor at all? It should be pretty straightforward to make a table of how many times each intron is observed using the GenomicRanges package.

[sulicon]

Thanks. It seems this package doesn't output the gap regions directly:

The genomic intervals between the "start" and "end" of the query that are "covered" by the alignment. Saying that the full [start,end] interval is covered is the same as saying that the alignment has no gap (no N in the CIGAR). It is then considered a simple alignment. Note that a simple alignment can have ismatches or deletions (in the reference). In other words, a deletion, encoded with a D, is NOT considered a gap.

So I still need calculate the gap regions from the CIGAR strings...

And I have just found a solution this afternoon. BEDTools (http://code.google.com/p/bedtools/) could calculate the coverage for each nucleotide with/without intronic regions counted. Then the difference is the the coverage/depth in the intronic regions.

[malcook]

tabulating gaps in bam files using bioconductor's GenomicRanges IS POSSIBLE

Thanks. It seems this package doesn't output the gap regions directly:
...

I think you may misunderstand the manual. In any case, the key word in your comment is perhaps "directly"

I've just blogged on using GenomicRanges for this exact purpose here. It works quite nicely.

And I have just found a solution this afternoon. BEDTools (http://code.google.com/p/bedtools/) could calculate the coverage for each nucleotide with/without intronic regions counted. Then the difference is the the coverage/depth in the intronic regions.

Hmmm. I wondered about this approach too. At the very least it requires you to know beforehand the location of all your introns. Let us know how you fare if you take this route....