[BioC] reducing hits from countGenomicOverlaps()

[Valerie Obenchain]

On 10/27/2011 12:44 PM, Robert Castelo wrote:
> Hi Valerie,
>
> On 10/27/11 7:53 PM, Valerie Obenchain wrote:
>> Hi Robert,
>>
>>
>> On 10/26/2011 03:34 PM, Robert Castelo wrote:
>>> 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?
>>
>> How do you want to tune the summing? If you want to sum all exons by
>> gene the piece of code Martin provided should do it. Was there another
>> goal you had in mind?
>
> well, Martin was warning me about whether just adding up exon counts 
> per gene makes complete sense and i guess the problem might arise with 
> overlapping exons of different transcripts since then by just summing 
> i'd be adding some counts twice. so i was wondering whether something 
> existed to address this summing along the lines of what is available 
> about the decision logic options for counting reads overlapping 
> annotations with countGenomicOverlaps().

Both the counting and summarizing problems have been addressed in the 
summarizedOverlaps function Martin mentioned. The function has different 
"modes"  that determine how to assign a read that overlaps multiple 
features. This function does not double count so a read is either 
assigned according to the rules of the "mode" or it is discarded. In the 
end you have counts per feature. A feature can be represented by a 
single range (ie, single exon, transcript or gene) or by multiple ranges 
(ie, exons in a gene, transcripts in a gene, etc.).  I'd be interested 
in feedback once you've had a chance to use the function.

Valerie



>
> cheers,
> robert.