Thanks Steve. That clarifies what is going on to me. I suppose the only way
to solve this to take the first gene id from an index that contains > 1
gene ids.

Fong


On Sun, Jan 13, 2013 at 2:36 AM, Steve Lianoglou <
mailinglist.honeypot@gmail.com> wrote:

> Hi Fong,
>
> On Sun, Jan 13, 2013 at 3:09 AM, Fong Chun Chan <fongchun@alumni.ubc.ca>
> wrote:
> > Hi,
> >
> > I was hoping someone could explain this weird situation to me when I am
> > trying to extract exon information using the GenomicFeatures Bioconductor
> > package.  Here is what I've done:
> >
> > txdb <- loadFeatures(txdbFile)
> > allExons <- exons(txdb, columns = c('gene_id', 'exon_id'))
> > annotDf <- values(allExons)
> > dim(annotDf)
> >
> > [1] 541825 2
> >
> > According to this code, I've extracted a total of 541825 exons which are
> > associated with x number of genes. Now when I do the following:
> >
> > length(unlist(annotDf[, 'gene_id']))
> >
> > [1] 546664
> >
> > How is it possible that there are suddenly an additional 4839 genes added
> > to the 'gene_id' column?
>
> The first clue is that the gene_id column is stored as a *List type,
> which you correctly realized because you are `unlisting` it it to get
> its length. The reason that it's a list is so that one could store
> more than one element per index, otherwise they would have just used a
> character (or factor) vector.
>
> So, what you can do is to see what's going on here -- I am using the
> latest version of R and the "TxDb.Hsapiens.UCSC.hg19.knownGene"
> package, so things might not be exactly the same for you, but I said
> up my vars to mimc yours, so:
>
> R> gene.ids <- annotDf$gene_id
> R> counts <- elementLengths(gene.ids) ## number of elements per bin
> R> table(counts)
> counts
>      1      2      3      4      5      6      9     15     22
> 282939   3788    101      7      6      1      4      3      3
>
> So we see that several gene ids can be assigned to one exon. You might
> as why, and I reckon the answer is that the exons are defined as
> "unique" simply by their chr,strand,start,end combo and several genes
> (or transcripts) might share the same exon (in "genomic space"). Let's
> see:
>
> R> allExons[elementLengths(gene.ids) > 1]
> GRanges with 3913 ranges and 2 metadata columns:
>          seqnames               ranges strand   |
> gene_id   exon_id
>             <Rle>            <IRanges>  <Rle>   |
> <CompressedCharacterList> <integer>
>      [1]     chr1 [  324288,   324345]      +   |
> 100132287,100133331        12
>      [2]     chr1 [10490159, 10490625]      +   |
> 100526739,378708       897
>      [3]     chr1 [10490804, 10491458]      +   |
> 100526739,378708       898
>      [4]     chr1 [10493899, 10494022]      +   |
> 100526739,378708       899
>      [5]     chr1 [10494714, 10494747]      +   |
> 100526739,378708       900
>      [6]     chr1 [10500404, 10500470]      +   |
> 100526739,378708       901
>      [7]     chr1 [10511434, 10512060]      +   |
> 100526739,1325       904
>      [8]     chr1 [16355621, 16355794]      +   |
> 1187,1188      1504
>      [9]     chr1 [19923471, 19923603]      +   |
> 100532736,440574      1788
>      ...      ...                  ...    ... ...
>  ...       ...
>   [3905]     chrY [26944570, 26944641]      -   |
> 57054,57135    274737
>   [3906]     chrY [26946960, 26947031]      -   |
> 57054,57135    274738
>   [3907]     chrY [26949403, 26949474]      -   |
> 57054,57135    274739
>   [3908]     chrY [26950875, 26951014]      -   |
> 57054,57135    274740
>   [3909]     chrY [26951104, 26951167]      -   |
> 57054,57135    274741
>   [3910]     chrY [26951604, 26951655]      -   |
> 57054,57135    274742
>   [3911]     chrY [26952216, 26952307]      -   |
> 57054,57135    274743
>   [3912]     chrY [26952582, 26952728]      -   |
> 57054,57135    274744
>   [3913]     chrY [26959330, 26959639]      -   |
> 57054,57135    274745
>
> If you hop to any of these regions in your favorite genome browser,
> you should find overlapping gene annotations there that share exact
> exon boundaries across transcripts.
>
> > The reason this is a problem is that I would like
> > to create a 3 column output where I use the command:
> >
> > data.frame( geneID = unlist(annotDf[, 'gene_id']), exonID = annotDf[,
> > 'exon_id'], exonCounts = exonCounts)
> >
> > Where exonCounts is generated with the countOverlaps() function using
> > allExons and reads from a bam file.  This command is failing because it
> > seems that the number of gene_ids is not pairing up with the exon_ids
> > properly. Is there something I am missing here?
>
> Hopefully this clarification helps?
>
> And also:
>
> > I am using 1.6.9 of GenomicFeatures, R is 2.14.2 and my txdb was built
> > using the command:
>
> You should upgrade to the latest version of R and bioconductor
> packages to get the best functionality (and help).
>
> Hope that helps,
> -steve
>
> --
> Steve Lianoglou
> Graduate Student: Computational Systems Biology
>  | Memorial Sloan-Kettering Cancer Center
>  | Weill Medical College of Cornell University
> Contact Info: http://cbio.mskcc.org/~lianos/contact
>

	[[alternative HTML version deleted]]

