The question is whether there is ever a use case to have a simple table.
This is analogous to base R's table and data.frame. For example, if you
call xtabs(), you get a table, then you have to call as.data.frame to get
back into the data.frame context. This is sort of clean, and we could
create an extension of table that for efficiency stores the associated
GRanges along with the counts in the .Data slot. Then as(x, "GRanges") on
that would generate the GRanges with the count column. That would be
complicated though.

Another issue is that table() cannot be used in the general way, due to
restrictions on dispatch with "...".

So I think I'm in favor of a new "count" generic. That naming is consistent
with countOverlaps, countSubjects, countQueries, etc.

Michael



On Wed, Dec 5, 2012 at 12:55 PM, Hervé Pagès <hpages@fhcrc.org> wrote:

> Hi,
>
> As Tim pointed out, there is the question of what the output should
> be.
>
> Proposal 1: Output is a "table" object.
>
> tableGenomicRanges1 <- function(x)
> {
>     x_levels <- sort(unique(x))
>     y <- IRanges:::matchIntegerQuads(**as.factor(seqnames(x)),
>                                      as.factor(strand(x)),
>                                      start(x), width(x),
>                                      as.factor(seqnames(x_levels)),
>                                      as.factor(strand(x_levels)),
>                                      start(x_levels), width(x_levels))
>     ans_names <- paste(seqnames(x_levels),
>                        "[", start(x_levels), ",", end(x_levels), "]",
>                        strand(x_levels), sep="")
>     ans <- table(factor(y, levels=seq_along(x_levels)))
>     names(ans) <- ans_names
>     ans
> }
>
> > tableGenomicRanges1(gr)
>    chr1[6,10]+     chr1[1,5]-    chr1[1,10]-     chr1[2,6]-     chr1[3,7]-
>            303            333            345            337            345
> chr1[101,110]- chr1[102,111]- chr1[103,112]-    chr1[5,10]*    chr2[2,10]+
>            327            308            343            308            338
>    chr2[3,10]+     chr2[4,8]-     chr2[5,9]-    chr2[6,10]- chr2[104,113]-
>            324            344            350            355            342
> chr2[105,114]- chr2[106,115]- chr2[107,116]-    chr2[4,10]*    chr3[7,10]+
>            354            338            318            336            351
>    chr3[8,10]+    chr3[7,11]-    chr3[8,12]-    chr3[9,10]-    chr3[9,13]-
>            355            360            357            324            325
>   chr3[10,10]-   chr3[10,14]- chr3[108,117]- chr3[109,118]- chr3[110,119]-
>            346            307            277            320            330
>
> but this output is not very convenient...
>
> Proposal 2: Output is a GenomicRanges object containing the
> sorted "levels". The counts are stored as a metadata col.
>
> tableGenomicRanges2 <- function(x)
> {
>     ans <- sort(unique(x))
>     names(ans) <- mcols(ans) <- NULL
>     y <- IRanges:::matchIntegerQuads(**as.factor(seqnames(x)),
>                                      as.factor(strand(x)),
>                                      start(x), width(x),
>                                      as.factor(seqnames(x_levels)),
>                                      as.factor(strand(x_levels)),
>                                      start(x_levels), width(x_levels))
>     mcols(ans)$count <- tabulate(y, nbins=length(x_levels))
>     ans
> }
>
> > tableGenomicRanges2(gr)
> GRanges with 30 ranges and 1 metadata column:
>        seqnames     ranges strand   |     count
>           <Rle>  <IRanges>  <Rle>   | <integer>
>    [1]     chr1 [  6,  10]      +   |       303
>    [2]     chr1 [  1,   5]      -   |       333
>    [3]     chr1 [  1,  10]      -   |       345
>    [4]     chr1 [  2,   6]      -   |       337
>    [5]     chr1 [  3,   7]      -   |       345
>    [6]     chr1 [101, 110]      -   |       327
>    [7]     chr1 [102, 111]      -   |       308
>    [8]     chr1 [103, 112]      -   |       343
>    [9]     chr1 [  5,  10]      *   |       308
>
>    ...      ...        ...    ... ...       ...
>   [22]     chr3 [  7,  11]      -   |       360
>   [23]     chr3 [  8,  12]      -   |       357
>   [24]     chr3 [  9,  10]      -   |       324
>   [25]     chr3 [  9,  13]      -   |       325
>   [26]     chr3 [ 10,  10]      -   |       346
>   [27]     chr3 [ 10,  14]      -   |       307
>   [28]     chr3 [108, 117]      -   |       277
>   [29]     chr3 [109, 118]      -   |       320
>   [30]     chr3 [110, 119]      -   |       330
>   ---
>   seqlengths:
>    chr1 chr2 chr3
>    1000 2000 1500
>
> Since the "table" method for GenomicRanges objects should really return
> a "table" object it should do tableGenomicRanges1().
> So what do we do for tableGenomicRanges2? Should we provide this
> functionality thru a new generic function e.g. count()?
> Note that tableGenomicRanges2() not only produces more convenient
> output, it's also 4x faster or more than tableGenomicRanges1() on
> big objects.
>
> Thanks,
> H.
>
>
>
> On 12/05/2012 11:18 AM, Michael Lawrence wrote:
>
>> I've often thought that this would be a really useful feature, that is,
>> counting the combination of seqname, start, end and strand. This would
>> follow the convention of the duplicated() method. All someone (Herve)
>> needs
>> to do is make an IRanges:::tableIntegerQuads. In the past, I have hacked
>> this together by tabulating a string representation. Given the way R
>> internalizes strings, that's not an ideal solution.
>>
>> Michael
>>
>> On Tue, Dec 4, 2012 at 4:25 PM, Vedran Franke <vfranke@bioinfo.hr> wrote:
>>
>>  Dear All,
>>>
>>> Is there an equivalent of the table function for a GenomicRanges object?
>>>
>>> Best,
>>>
>>> Vedran
>>>
>>>
>>> --
>>> Vedran Franke
>>>
>>> Bioinformatics Group,
>>> Department of Molecular Biology,
>>> Faculty of Science, Zagreb
>>>
>>> ______________________________**_________________
>>> Bioconductor mailing list
>>> Bioconductor@r-project.org
>>> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https://stat.ethz.ch/mailman/listinfo/bioconductor>
>>> Search the archives:
>>> http://news.gmane.org/gmane.**science.biology.informatics.**conductor<http://news.gmane.org/gmane.science.biology.informatics.conductor>
>>>
>>>
>>         [[alternative HTML version deleted]]
>>
>>
>> ______________________________**_________________
>> Bioconductor mailing list
>> Bioconductor@r-project.org
>> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https://stat.ethz.ch/mailman/listinfo/bioconductor>
>> Search the archives: http://news.gmane.org/gmane.**
>> science.biology.informatics.**conductor<http://news.gmane.org/gmane.science.biology.informatics.conductor>
>>
>>
> --
> Hervé Pagès
>
> Program in Computational Biology
> Division of Public Health Sciences
> Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N, M1-B514
> P.O. Box 19024
> Seattle, WA 98109-1024
>
> E-mail: hpages@fhcrc.org
> Phone:  (206) 667-5791
> Fax:    (206) 667-1319
>

	[[alternative HTML version deleted]]

