On Tue, Sep 14, 2010 at 5:43 PM, Elizabeth Purdom <epurdom@stat.berkeley.edu
> wrote:

>  Hello,
>
> I would like to get the range of genes from a TranscriptDb object. Ideally,
> I would like to do:
>
> txdb <- loadFeatures(system.file("extdata", "UCSC_knownGene_sample.sqlite",
>                                   package="GenomicFeatures"))
> genes(txdb)
>
> And have this return the entire genomic range of the gene (much like
> transcripts() returns the genomic range of the transcript, including the
> introns). However, this does not exist. So I have made a little function to
> do this:
>
> genes<-function(txdb){
>    txByGene<-transcriptsBy(txdb,"gene")
>    geneStarts<-sapply(start(txByGene),min)
>    geneEnds<-sapply(end(txByGene),max)
>    geneChr<-sapply(seqnames(txByGene),unique)
>    geneStrand<-sapply(strand(txByGene),unique)
>  GRanges(seqnames=geneChr,strand=geneStrand,ranges=IRanges(geneStarts,geneEnds),gene_name=names(txByGene))
>
> }
>
> However, this has run me into all kinds of problems because the same gene
> name can be used for transcripts on different chromosomes, namely the X and
> Y chromosome (which makes me wonder about the safety of the assumption that
> gene names are unique -- does this count as unique?).
>

http://en.wikipedia.org/wiki/Pseudoautosomal_region


>
> Before I go further down this route, is there not a simple way to get this
> information?
>
>
Stating the obvious, probably, but this is not a question with a simple
answer.  You may have to settle on some sort of heuristic (choose a random
region, the longest region, the one containing the most transcripts, etc.)
to sort out such cases or simply ignore them altogether (throw them away).
 Unfortunately, these oddities are present in the human genome (and other
expanded genomes) and it is not possible to totally avoid them.

Sean

	[[alternative HTML version deleted]]

