leekgroup / derfinder

R package for DER Finder, a method for differential expression analysis of RNA-seq data
http://biostatistics.oxfordjournals.org/content/early/2014/01/06/biostatistics.kxt053.short?rss=1
MIT License
22 stars 11 forks source link

exons length error from stopifnot() call in getFlags.R #12

Closed brianhigh closed 9 years ago

brianhigh commented 9 years ago

When trying to reproduce results in the 2013 paper, from this code in analysis_code.R ...

# get the flags:
exons = getAnnotation("hg19","knownGene")
myflags = getFlags(regions = regions.merged.y, exons, "chrY", pctcut = 0.8)

... I am getting this fatal error ...

Error: length(unique(exons$chr)) == 1 is not TRUE

Some tests:

> length(unique(subset(x=exons, seqnames == "chrY", c(chr))))
 [1] 1
> length(unique(unlist(subset(x=exons, seqnames == "chrY", c(chr)))))
 [1] 101
> str(unique(unlist(subset(x=exons, seqnames == "chrY", c(chr)))))
 chr [1:101] "100101116" "100101120" "100132596" "100133941" "100289150" ...
> summary(subset(x=exons, seqnames == "chrY", c(chr)))
     chr           
 Length:1190       
 Class :character  
 Mode  :character  
> head(subset(x=exons, seqnames == "chrY"))
    gene       chr seqnames   start     end width strand
392   85 100101116     chrY 6258442 6258716   275      +
393   85 100101116     chrY 6262141 6262300   160      +
394   85 100101116     chrY 6269164 6269272   109      +
395   85 100101116     chrY 6271629 6271766   138      +
396   85 100101116     chrY 6279348 6279605   258      +
397   85 100101116     chrY 9590765 9591022   258      -
> tail(subset(x=exons, seqnames == "chrY"))
        gene  chr seqnames    start      end width strand
258853 22527 9189     chrY  2354455  2358810  4356      -
258854 22527 9189     chrY  2354455  2358813  4359      -
258855 22527 9189     chrY  2368352  2368580   229      -
258856 22527 9189     chrY  2368858  2369015   158      -
263606 22923 9426     chrY 20137667 20139627  1961      +
263607 22923 9426     chrY 19990140 19992100  1961      -

See: getFlags.R ...

stopifnot(length(unique(exons$chr))==1,
length(unique(regions$chr)) == 1,
unique(exons$chr) == unique(regions$chr))

This is my sessionInfo():

> sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C         LC_TIME=C            LC_COLLATE=C        
 [5] LC_MONETARY=C        LC_MESSAGES=C        LC_PAPER=C           LC_NAME=C           
 [9] LC_ADDRESS=C         LC_TELEPHONE=C       LC_MEASUREMENT=C     LC_IDENTIFICATION=C 

attached base packages:
 [1] splines   grid      parallel  stats4    stats     graphics  grDevices utils     datasets 
[10] methods   base     

other attached packages:
 [1] rtracklayer_1.26.2     GenomicFeatures_1.18.3 AnnotationDbi_1.28.1  
 [4] derfinder_1.0.2        locfdr_1.1-7           devtools_1.7.0        
 [7] HiddenMarkov_1.8-1     Genominator_1.20.0     GenomeGraphs_1.26.0   
[10] DESeq_1.18.0           lattice_0.20-29        locfit_1.5-9.1        
[13] Biobase_2.26.0         edgeR_3.8.5            limma_3.22.4          
[16] biomaRt_2.22.0         Rsamtools_1.18.2       Biostrings_2.34.1     
[19] XVector_0.6.0          GenomicRanges_1.18.4   GenomeInfoDb_1.2.4    
[22] IRanges_2.0.1          S4Vectors_0.4.0        BiocGenerics_0.12.1   
[25] BiocInstaller_1.16.1   sqldf_0.4-10           gsubfn_0.6-6          
[28] proto_0.3-10           RSQLite_1.0.0          DBI_0.3.1             

loaded via a namespace (and not attached):
 [1] BBmisc_1.9              BatchJobs_1.5           BiocParallel_1.0.3     
 [4] GenomicAlignments_1.2.1 RColorBrewer_1.1-2      RCurl_1.95-4.5         
 [7] XML_3.98-1.1            annotate_1.44.0         base64enc_0.1-2        
[10] bitops_1.0-6            brew_1.0-6              checkmate_1.5.1        
[13] chron_2.3-45            codetools_0.2-10        digest_0.6.8           
[16] fail_1.2                foreach_1.4.2           genefilter_1.48.1      
[19] geneplotter_1.44.0      httr_0.6.1              iterators_1.0.7        
[22] sendmailR_1.2-1         stringr_0.6.2           survival_2.37-7        
[25] tcltk_3.1.2             tools_3.1.2             xtable_1.7-4           
[28] zlibbioc_1.12.0
lcolladotor commented 9 years ago

Have you tried this with R 2.15.x?

brianhigh commented 9 years ago

I tried today with R 2.15.x and this issue does not come up with R 2.15.3. So, this issue can be closed.

lcolladotor commented 9 years ago

Awesome =)