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 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?

I am using 1.6.9 of GenomicFeatures, R is 2.14.2 and my txdb was built
using the command:

txdb <- makeTranscriptDbFromUCSC(genome = 'hg19', tablename = 'ensGene')

Has anyone encountered this? Any help would be great.

Thanks, Fong

	[[alternative HTML version deleted]]

