Thanks so much Valerie.
This sounds like a good approach. I did not even know there is a way to
read GFFs directly from Granges. I was using rtracklayer before.

I will follow this advice and see how it goes.

I really appreciate your help.

Gowthaman

On Tue, May 29, 2012 at 9:03 AM, Valerie Obenchain <vobencha@fhcrc.org>wrote:

> Hi Gowthaman,
>
> Here is one approach,
>
> (1) create the GRanges or GRangesList from gff
>
> Using a gff file in GenomicFeatures as an example,
>
> >  gffFile <- system.file("extdata","a.gff3"**,package="GenomicFeatures")
> >      txdb <- makeTranscriptDbFromGFF(file=**gffFile,
>   +               format="gff3",
>   +               dataSource="partial gtf file for Tomatoes for testing",
>   +               species="Solanum lycopersicum",
>   +               circ_seqs=DEFAULT_CIRC_SEQS,
>   +               miRBaseBuild=NULL)
>
> If you are interested in genes either exons by genes or transcripts by
> genes might be useful,
>
>  exbygene <- exonsBy(txdb, "gene")
>  txbygene <- transcriptsBy(txdb, "gene")
>
> The gene names are associated with each outer list element of the
> GRangesList,
>
> > exbygene
>  GRangesList of length 488:
>  $gene:Solyc00g005000.2
>  GRanges with 2 ranges and 2 elementMetadata cols:
>          seqnames         ranges strand |   exon_id
> exon_name
> <Rle> <IRanges> <Rle> | <integer> <character>
>    [1] SL2.40ch00 [16437, 17275]      + |         1
> exon:Solyc00g005000.2.1.1
>    [2] SL2.40ch00 [17336, 18189]      + |         2
> exon:Solyc00g005000.2.1.2
>
>  $gene:Solyc00g005020.1
>  GRanges with 3 ranges and 2 elementMetadata cols:
>          seqnames         ranges strand | exon_id                 exon_name
>    [1] SL2.40ch00 [68062, 68211]      + |       3 exon:Solyc00g005020.1.1.1
>    [2] SL2.40ch00 [68344, 68568]      + |       4 exon:Solyc00g005020.1.1.2
>    [3] SL2.40ch00 [68654, 68764]      + |       5 exon:Solyc00g005020.1.1.3
>
>
> (2) findOverlaps()
>
> For the sake of example, I'll take a subset of one of the GRangesLists and
> overlap it with the full list,
>
>  subset <- exbygene[c(2,4,6)]
>  olap <- findOverlaps(subset, exbygene)
>
> > olap
> Hits of length 3
> queryLength: 3
> subjectLength: 488
>  queryHits subjectHits
> <integer> <integer>
>  1         1           2
>  2         2           4
>  3         3           6
>
> The information in the 'Hits' object can be used to map back to the query
> and subject used in findOverlaps().
> See ?queryHits and ?subjectHIts.
>
> Next identify which subjects were hit by the queries and extract the gene
> names,
>
>  names(exbygene[subjectHits(**olap)])
>
>
> (3) reconsiling gene sets
>
> Take a look at the 'type' argument to findOverlaps(). To find exact
> matches you'll want to find ranges with type='equal', to find any overlaps
> use type='any' (default).
>
> You may also want to look at summarizeOverlaps(). It is an extension of
> findOverlaps() that counts a read a maximum of once (i.e., no double
> counting). If a read hits multiple features different rule sets are applied
> to assign the read to a single feature.
>
> Valerie
>
>
>
>
> On 05/24/2012 10:05 AM, gowtham wrote:
>
>> Hi Everyone,
>> This is my first attempt in R/bioconductor outside of plot()/hist()
>> function. Thanks to Martin Morgan, I was impressed with what
>> bioconductor can do for me ( after attending bioc workshop last week).
>>
>> When I try to findOverlap on two GFF files I get information on which
>> feature is overlapping with another feature. But, instead of using
>> "feature
>> names" it just indicates 1,2, 3...etc. Is there a way, i can make it
>> output
>> gene names. (my GRanges objects are read from gff files using Rtracklayer&
>> later have names (geneids) assigned to them using names() function).
>>
>>  findOverlaps(au, am)
>>>
>> Hits of length 1740
>> queryLength: 790
>> subjectLength: 1054
>>      queryHits subjectHits
>>       <integer>    <integer>
>>  1            1           1
>>  2            1         528
>>  3            2          13
>>  4            2         540
>>  5            3         136
>>
>>
>>
>> And here is the background of what i am trying to do. In case, if any of
>> you have a better suggestion:
>> I run 3 de novo gene prediction algorithm on our newly assembled genome.
>> And I would like to reconcile the gene sets from all of them to genes (a)
>> predicted with exact same boundaries in all three (b) boundaries not same
>> but overlap (c) predicted by only one or two of them.
>>
>> And use reduce in case of (b) to obtain outer most boundaries of
>> overlapping prediction.
>>
>> As a novice users its bit overwhelming for me.
>>
>> Thanks very much,
>> Gowthaman
>>
>>
>


-- 
Gowthaman

Bioinformatics Systems Programmer.
SBRI, 307 West lake Ave N Suite 500
Seattle, WA. 98109-5219
Phone : LAB 206-256-7188 (direct).

	[[alternative HTML version deleted]]

