[BioC] reducing hits from countGenomicOverlaps()

[Robert Castelo]

hi,

On 10/27/11 12:07 AM, Martin Morgan wrote:
> On 10/25/2011 05:35 PM, Robert Castelo wrote:
>> dear list,
>>
>> the following three lines allow one to count overlaps of aligned
>> short-reads with annotations:
>>
>> aln <- readGappedAlignments("somebamfile.bam")
>> txdb <- makeTranscriptFromUCSC(genome="hg19", tablename="ensGene")
>> ensGenes <- exonsBy(txdb, by="gene")
>> ov <- countGenomicOverlaps(aln, ensGenes)
>>
>> then i want to get read-counts per gene and the first thing that comes
>> to my head is doing:
>>
>> counts <- sapply(ov, function(x) sum(values(x)[["hits"]]))
>>
>> which goes through every gene and adds up the "hits" of its exons.
>> however, this latter step of "just adding" takes longer than the actual
>> calculation of the hits with countGenomicOverlaps() and i guess that
>> there are more efficient ways to approach this, probably something
>> around "reducing the hits value column". i've been looking at rdapply()
>> and reduce() and googled too, but couldn't find anything, so i look
>> forward to your suggestions.
>
> Hi Robert -- A strategy here is along the lines (untested!) of
>
> hits <- values(unlist(ov)))[["hits"]]
> genes <- rep(names(ov), elementLengths(ov))
> counts <- sapply(split(hits, genes), sum)

beautiful and fast, i've just tested it with the ~ 50,000 Ensembl genes 
and the execution time is nearly negligible, less than half a second in 
my laptop.

> but you'll want to make sure that this is a conceptually sensible way of
> counting hits per gene.

i'm aiming at counting exonic-reads and adding them up at gene level. 
the function countGenomicOverlaps() allows me to tune the part of the 
logic related to counting exonic reads, but is there any function that 
would allow me to tune the summing of exonic reads?

> The countGenomicOverlaps function has been heavily updated to
> summarizeOverlaps in 'devel'; you might want to check that out, or at
> least the vignette for summarizeOverlaps at
>
> http://bioconductor.org/packages/devel/bioc/html/GenomicRanges.html

that sounds great, i'll upgrade next week when it becomes release.

thanks again!
robert.