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

>
> >
> >
> >     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.
> >
> I understand the possible pitfalls (and the biology for why they would
> happen), but at the same time, it is still a pretty standard question
> that many databases give some answer to. Generally, I understand gene
> boundaries as the largest region that contains the transcripts annotated
> to be in that gene. This is the case for Ensembl, I believe (where you
> can definitely query the database and get out gene start and stop points
> for a gene id). If there are multiple chromosomes for the same gene id,
> these could be returned as separate entries. I would like to have a
> database query that accomplishes this,


Hi, Elizabeth.  Just to answer your question directly....

It is possible to do this type of thing in SQL and connecting to the public
mysql server at UCSC is possible using RMySQL or some other mysql client
(see here: http://genome.ucsc.edu/FAQ/FAQdownloads.html#download29).  I
provide a few (very lightly tested) examples.  Each runs in less than a
second on our local UCSC mysql mirror.

This little SQL grabs the min(start) and max(end) for each locuslink
directly from the UCSC database.

select value as locuslink,min(txStart) as s,max(txEnd) as e,chrom
  from knownGene kg
  join knownToLocusLink kl
  on kg.name=kl.name
  group by value,chrom;

The next just counts the number of such records returned from above that end
up  with records on more than one chromosome (281 LocusLinks, hg18,
knownGene)

select locuslink,count(*) from
  (select value as locuslink,min(txStart) as s,max(txEnd) as e,
       max(txEnd)-min(txStart) as d,chrom
       from knownGene kg
       join knownToLocusLink kl
       on kg.name=kl.name group by value,chrom) t
       group by t.locuslink having count(*)>1;

This final one pulls out the single longest LocusLink record for each
LocusLink.  Results will then be unique.

select t.* from
  (select value as locuslink,min(txStart) as s,max(txEnd) as e,
       max(txEnd)-min(txStart) as d,chrom
       from knownGene kg
       join knownToLocusLink kl
       on kg.name=kl.name group by value,chrom) t
       group by t.locuslink having d=max(d);

I use LocusLink here to match the ancient UCSC terminology, but, of course,
we are talking about Entrez Gene IDs.

Hope that helps a bit.

Sean



> like the transcripts and exons
> functions, because my naive implementation (in addition to the problems
> with the multiple chrX and chrY mappings) is extremely slow.
>
> The reason I want this is that I want to pull up all reads from a bam
> file associated with this gene, regardless of the transcript. If I give
> per transcript information, I get reads coming up multiple times,
> because the transcripts from the same gene substantially overlap and the
> reads are pulled up per transcript. This might be useful in some
> setting, but generally I'm going to want to then post-process the reads
> per gene to deal with the multiple transcript problem (yes I know genes
> can overlap too, but this is a problem of smaller magnitude, and I think
> I can deal with that issue separately). Moreover, the number of genes is
> much smaller than the number of transcripts, and I want to reduce the
> number of queries of the bam file, as well as the number of reads I pull
> into memory at a time.
>
> Thanks,
> Elizabeth
>
>
>
>
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>

	[[alternative HTML version deleted]]

