[Bioc-sig-seq] GenomicFeatures package

[Martin Morgan]

Hi Steve --

Steve Lianoglou <mailinglist.honeypot at gmail.com> writes:

> Hi,
>
> On Oct 28, 2009, at 3:44 PM, Martin Morgan wrote:
> <snip>
>>>>
>>> It could be useful, but you'll need to get it into the shape
>>> expected by the
>>> GenomicFeatures functions. With specific regard to transcripts(),
>>> it will
>>> need a 'chrom' column with the chromosome names, 'name' column for
>>> the gene
>>> name, and 'txStart' and 'txEnd' for the start and end of the
>>> transcripts,
>>> using UCSC coordinate conventions.
>>>
>>> I'm now thinking that it would be more convenient for the user if
>>> there was
>>> a transcripts method on a RangesList object, which could provide the
>>> necessary information. This could be extracted from the RangedData,
>>> which
>>> rtracklayer can create from your GFF file, with one caveat: if the
>>> GFF file
>>> contains a mixture of features (genes, exons, etc) and relies on the
>>> hierarchical features of GFF, it will take more work to get things
>>> into the
>>> right shape.
>>>
>>> The question then is where would this new functionality belong? The
>>> future
>>> of the GenomicFeatures package was a bit uncertain, the last time I
>>> checked.
>>
>> The intention is that GenomicFeatures will mature to contain these
>> sorts of data structures and functions, so that it's easy to create or
>> retrieve transcript or exon information.
> </snip>

> I've actually been working on a package of my own that does a lot of
> this stuff, since I've been working lately with rna-seq data.
>
> It takes an input file similar to what you get from the ucsc table
> browser with gene, tx/cds -start/-end, exon, etc. columns and whips up
> an annotation package for the genome/annotation source you are using
> (info is stored in an SQLite database). You can then query the
> database for genes of interest by name, chromosome, or bounded (start/
> end) chromosome position for further downstream use.
>
> Usually results are returned to the user as lists of gene objects,
> which the user would have to whip into a data.frame format if that's
> what's desired.
>
> I'll post some code below to show you what I mean, and see if you all
> think it's helpful. I'm not sure that it's ready for public
> consumption, but I guess I just wanted to mention that it was my
> intention to develop this a bit further in the near future. If it's
> useful for other people, all the better.

Just wanted to acknowledge this email. I like how straight-forward
your interface is, and the use of appropriate classes. RTree and the
custom compile is a little problematic; I wonder whether it's really
necessary for data this size (i.e., not too large) or whether there
are other ways in which the query could be optimized (representing the
data in-memory and using IRanges::findOverlaps?).

What ever the case, it would be good to coordinate your activities
with what we do here in Seattle and with GenomicFeatures, so we'll try
and keep the channels of communication open.

Martin

> Unfortunately I need to put much of this on hold for the immediate
> future, since I have to prepare for my upcoming qualifying exam ...
>
> I realize it somehow does similar things to what the GenomicFeatures
> package might have done, but I needed this exact functionality for my
> immediate analysis needs, and thought I'd have a go at it. The SQLite
> database that's generated actually stores gene, transcript, and exon
> info in RTree indices[1], so ranged queries (like all genes/
> transcripts between some start/stop positions) are super fast. It
> requires a custom compile of SQLite, though
>
> Examples of how I use it are pasted below.
>
> -steve
>
> [1] RTree: http://www.sqlite.org/rtree.html
>
> ##############################
> # Using with refseq annotations for hg18
>
> R> library(GenomeAnnotations) # name of my library
> R> hg18.refseq <- GenomeDB("hg18", "refseq")
> R> d1r <- Gene("DICER1", genome=hg18.refseq)
>
> ## Getting transcript info
>
> R> transcripts(d1r)
> SimpleRangesList instance:
> $NM_030621
> IRanges instance:
>          start      end width
> [1]  94622320 94626752  4433
> [2]  94627120 94627200    81
> [3]  94627296 94627456   161
> [4]  94629976 94630248   273
> [5]  94631912 94632800   889
> [6]  94635872 94636024   153
> [7]  94639432 94640216   785
> [8]  94641160 94641336   177
> [9]  94641768 94641872   105
> ...       ...      ...   ...
> [20] 94660288 94660760   473
> [21] 94662672 94662840   169
> [22] 94665560 94665720   161
> [23] 94666144 94666280   137
> [24] 94667600 94667728   129
> [25] 94668608 94668768   161
> [26] 94669408 94669592   185
> [27] 94676752 94676848    97
> [28] 94677776 94677808    33
>
> <1 additional element>
>
> ## UTR Info
> R> getUtr3(d1r)
> $NM_030621
> IRanges instance:
>         start      end width
> [1] 94622320 94626587  4268
>
> $NM_177438
> IRanges instance:
>         start      end width
> [1] 94622320 94626587  4268
>
> R> getUtr5(d1r)
> $NM_030621
> IRanges instance:
>         start      end width
> [1] 94669548 94669592    45
> [2] 94676752 94676848    97
> [3] 94677776 94677808    33
>
> $NM_177438
> IRanges instance:
>         start      end width
> [1] 94669548 94669592    45
> [2] 94693320 94693512   193
>
>
> ## Get genes on chromosome 1
> R> genes.1 <- getGenesOnChromosome(hg18.refseq, 1)
> R> names(genes.1)[1:5]
> [1] "WASH5P"       "LOC100288778" "FAM138F"      "FAM138A"
> "FAM138C"
>
> ## Get genes on chr1 between pos 1e6 and 1e7
> R> genes.1.sub <- getGenesOnChromosome(hg18.refseq, 1, start=1e6,
> end=1e7)
> R> names(genes.1.sub)[1:5]
> [1] "C1orf159" "TTLL10"   "TTLL10"   "TNFRSF18" "TNFRSF18"
>
> #########################################
> # or, hg18 aceview annos
> R> hg18.ace <- GenomeDB("hg18", "aceview")
> R> d1a <- Gene("DICER1", genome=hg18.ace)
> R> transcripts(d1a)
> SimpleRangesList instance:
> $DICER1.cApr07
> IRanges instance:
>          start      end width
> [1]  94622320 94626752  4433
> [2]  94627120 94627200    81
> [3]  94627296 94627456   161
> [4]  94629976 94630248   273
> [5]  94631912 94632800   889
> [6]  94635872 94636024   153
> [7]  94639432 94640216   785
> [8]  94641160 94641336   177
> [9]  94641768 94641872   105
> ...       ...      ...   ...
> [19] 94653712 94653840   129
> [20] 94660288 94660760   473
> [21] 94662672 94662840   169
> [22] 94665560 94665720   161
> [23] 94666144 94666280   137
> [24] 94667600 94667728   129
> [25] 94668608 94668768   161
> [26] 94669408 94669592   185
> [27] 94693320 94693512   193
>
> <10 additional elements>
>
> R> getUtr3(d1a)
> $DICER1.cApr07
> IRanges instance:
>         start      end width
> [1] 94622320 94626587  4268
>
> $DICER1.bApr07
> IRanges instance:
>         start      end width
> [1] 94622320 94626587  4268
>
> .... etc ...

-- 
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793