[BioC] discrepancy between gene models in easyRNASeq and provided exon annotation

[Johanna Schott [guest]]

Dear list (and dear Nico),

I encountered the following problem while annotating RNASeq reads with the UCSC table refGene (hg19):
The exon coordinates of the gene model that is computed by easyRNAseq are quite different from the exon coordinates that I provide in my (self-made) annotation object.

Here is one gene plus my code as an example:

> # load required packages
> library("easyRNASeq")
> library(BSgenome.Hsapiens.UCSC.hg19)
> 
> # get chromosome sizes
> chr.sizes <- seqlengths(Hsapiens)
> 
> # get annotation object
> annot <- load("gAnnot_refGene_hg19.rda")
> 
> # compare gene model for MYC with exons in annotation object
> 
> object_RNASeq <- easyRNASeq(getwd(), filename = c("sample1.sorted.bam"),
+ readLength = 51L,
+ organism = "Hsapiens",
+ chr.sizes = chr.sizes,
+ format = "bam",
+ annotationMethod = "env",
+ annotationObject = exon_range,
+ count = "genes",
+ outputFormat = "RNAseq",
+ summarization = "geneModels"
+ )

> model <- geneModel(object_RNASeq)
> model[which( model$gene == "MYC" ), ] # here you see the exon coordinates of the gene model

RangedData with 3 rows and 4 value columns across 24 spaces
     space                 ranges |        gene      strand     transcript
  <factor>              <IRanges> | <character> <character>    <character>
1     chr8 [126085394, 126085566] |         MYC           + MYC_transcript
2     chr8 [126087239, 126087353] |         MYC           + MYC_transcript
3     chr8 [126088589, 126088742] |         MYC           + MYC_transcript
         exon
  <character>
1       MYC_1
2       MYC_2
3       MYC_3

> 
> 
> exon_range[which( exon_range$gene == "MYC" ), ] # here you see the exon coordinates that I provided
RangedData with 3 rows and 4 value columns across 24 spaces
     space                 ranges |      strand  transcript        gene
  <factor>              <IRanges> | <character> <character> <character>
1     chr8 [128748314, 128748869] |           +   NM_002467         MYC
2     chr8 [128750493, 128751265] |           +   NM_002467         MYC
3     chr8 [128752641, 128753680] |           +   NM_002467         MYC
                                 exon
                          <character>
1 NM_002467_exon_0_0_chr8_128748315_f
2 NM_002467_exon_1_0_chr8_128750494_f
3 NM_002467_exon_2_0_chr8_128752642_f


I checked genes on different chromosomes, genes with one or several exons. None of this seems to matter. 

What could be the problem? I am completely without a clue...

Thanks a lot for any reply!

Johanna


 -- output of sessionInfo(): 

R version 2.15.1 (2012-06-22)
Platform: x86_64-pc-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] BSgenome.Hsapiens.UCSC.hg19_1.3.19 easyRNASeq_1.4.2                   ShortRead_1.16.2                   latticeExtra_0.6-24               
 [5] RColorBrewer_1.0-5                 Rsamtools_1.10.2                   DESeq_1.10.1                       lattice_0.20-6                    
 [9] locfit_1.5-8                       BSgenome_1.26.1                    GenomicRanges_1.10.5               Biostrings_2.26.2                 
[13] IRanges_1.16.4                     edgeR_3.0.4                        limma_3.14.1                       biomaRt_2.14.0                    
[17] Biobase_2.18.0                     genomeIntervals_1.14.0             BiocGenerics_0.4.0                 intervals_0.13.3                  

loaded via a namespace (and not attached):
 [1] annotate_1.36.0      AnnotationDbi_1.20.3 bitops_1.0-5         DBI_0.2-5            genefilter_1.40.0    geneplotter_1.36.0  
 [7] grid_2.15.1          hwriter_1.3          RCurl_1.95-3         RSQLite_0.11.2       splines_2.15.1       stats4_2.15.1       
[13] survival_2.36-14     XML_3.95-0.1         xtable_1.7-0         zlibbioc_1.4.0      
> 


--
Sent via the guest posting facility at bioconductor.org.