[BioC] Problems in running easyRNASeq

pepap jpasulka at gmail.com
Fri Mar 29 12:06:33 CET 2013


Hello,

I have same problem as reported previously (https://stat.ethz.ch/pipermail/bioconductor/2013-January/050555.html), "Error in .doBasicCount(obj)". How and what to fix ?

>count.table <- easyRNASeq(
+ filenames=c("siCTRL_5314.bam",
+  "siCTRL_5315.bam",
+  "siZ9_5310.bam",
+  "siZ9_5312.bam"),
+ filesDirectory="/Volumes/Cellar/DNASeq_20121016/02.human/Orig_tophat_mapping",
+ organism="Hsapiens", 
+ chr.sizes=as.list(seqlengths(Hsapiens)), 
+ readLength=101L, 
+ annotationMethod="rda", 
+ annotationFile="./exon.annotation.biomart.rda",
+ format="bam", 
+ count="exons")
Checking arguments... 
Fetching annotations... 
Summarizing counts... 
Processing siCTRL_5314.bam
Error in .doBasicCount(obj) : The genomicAnnotation slot is empty

>sessionInfo()
R version 2.15.2 (2012-10-26)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
 [1] BSgenome.Hsapiens.UCSC.hg19_1.3.19 easyRNASeq_1.4.2                  
 [3] ShortRead_1.16.4                   latticeExtra_0.6-24               
 [5] RColorBrewer_1.0-5                 Rsamtools_1.10.2                  
 [7] DESeq_1.10.1                       lattice_0.20-14                   
 [9] locfit_1.5-8                       BSgenome_1.26.1                   
[11] GenomicRanges_1.10.7               Biostrings_2.26.3                 
[13] edgeR_3.0.8                        limma_3.14.4                      
[15] Biobase_2.18.0                     genomeIntervals_1.14.0            
[17] intervals_0.14.0                   IRanges_1.16.6                    
[19] BiocGenerics_0.4.0                 biomaRt_2.14.0                    

loaded via a namespace (and not attached):
 [1] annotate_1.36.0      AnnotationDbi_1.20.7 bitops_1.0-4.2      
 [4] DBI_0.2-5            genefilter_1.40.0    geneplotter_1.36.0  
 [7] grid_2.15.2          hwriter_1.3          RCurl_1.95-4.1      
[10] RSQLite_0.11.2       splines_2.15.2       stats4_2.15.2       
[13] survival_2.37-4      XML_3.95-0.2         xtable_1.7-1        
[16] zlibbioc_1.4.0

>samtools view -H siCTRL_5314.bam
@HD	VN:1.0	SO:coordinate
@SQ	SN:chr1	LN:249250621
@SQ	SN:chr10	LN:135534747
@SQ	SN:chr11	LN:135006516
@SQ	SN:chr12	LN:133851895
@SQ	SN:chr13	LN:115169878
@SQ	SN:chr14	LN:107349540
@SQ	SN:chr15	LN:102531392
@SQ	SN:chr16	LN:90354753
@SQ	SN:chr17	LN:81195210
@SQ	SN:chr18	LN:78077248
@SQ	SN:chr19	LN:59128983
@SQ	SN:chr2	LN:243199373
@SQ	SN:chr20	LN:63025520
@SQ	SN:chr21	LN:48129895
@SQ	SN:chr22	LN:51304566
@SQ	SN:chr3	LN:198022430
@SQ	SN:chr4	LN:191154276
@SQ	SN:chr5	LN:180915260
@SQ	SN:chr6	LN:171115067
@SQ	SN:chr7	LN:159138663
@SQ	SN:chr8	LN:146364022
@SQ	SN:chr9	LN:141213431
@SQ	SN:chrM	LN:16571
@SQ	SN:chrX	LN:155270560
@SQ	SN:chrY	LN:59373566
@PG	ID:TopHat	VN:1.3.2	CL:tophat -p 8 genome  siCTRL_5314_R01.txt siCTRL_5314_R02.txt

I did not use GTF file for annotation, the ".rda" file was created as follows:
#################################################################
library(biomaRt)

ensembl<-useMart("ensembl", dataset="hsapiens_gene_ensembl")
exon.annotation<-getBM(
 c("external_gene_id",
  "ensembl_gene_id",
  "strand",
  "ensembl_transcript_id",
  "chromosome_name",
  "ensembl_exon_id",
  "exon_chrom_start",
  "exon_chrom_end"
  ),
 mart=ensembl
)
exon.annotation$chromosome <- paste( 
 "chr", 
 exon.annotation$chromosome_name, 
 sep="")

library(IRanges)

gAnnot <- RangedData(
 IRanges(start=exon.annotation$exon_chrom_start,
  end=exon.annotation$exon_chrom_end
        ),
  space=exon.annotation$chromosome,
  strand=exon.annotation$strand,
  transcript=exon.annotation$ensembl_transcript_id,
  gene=exon.annotation$ensembl_gene_id,
  exon=exon.annotation$ensembl_exon_id,
  gene_name=exon.annotation$external_gene_id,
  universe="Hsap"
                    )
save(gAnnot, file="exon.annotation.biomart.rda")
#################################################################

Thanks in advance.

Josef Pasulka
CEITEC muni
pepap at chemi.muni.cz


More information about the Bioconductor mailing list