Bug in 'gsea' function ('phenoTest' package)
0
1
Entering edit mode
@damien-croteau-chonka-13392
Last seen 4.5 years ago
Boston, MA, USA

Bug Description

I've found a bug in the "gsea" function in the "phenoTest" Bioconductor package that I've been able to reproduce below:

library("phenoTest")
library("GSEABase")
library("devtools")

# setup test data as in package vignette
set.seed(777)

data(eset)

Tumor.size <- rnorm(ncol(eset), 50, 2)
pData(eset) <- cbind(pData(eset), Tumor.size)
pData(eset)[1:20, 'lymph.node.status'] <- 'positive'

survival <- matrix(c("Relapse", "Months2Relapse"),
                   ncol = 2,
                   byrow = TRUE)
colnames(survival) <- c('event', 'time')
vars2test <-
    list(survival = survival,
         categorical = 'lymph.node.status',
         continuous = 'Tumor.size')

epheno <- ExpressionPhenoTest(eset,vars2test,p.adjust.method='none')

sign1 <- sample(featureNames(eset))[1:20]
sign2 <- sample(featureNames(eset))[1:50]
mySignature <- list(sign1, sign2)
names(mySignature) <- c('My first signature', 'Another signature')

myGeneSetA <- GeneSet(geneIds = sign1, setName = 'My first signature')
myGeneSetB <- GeneSet(geneIds = sign2, setName = 'Another signature')
mySignature <- GeneSetCollection(myGeneSetA, myGeneSetB)

The example gene set enrichment analysis (GSEA) from the vignette runs fine:

my.gsea <-
    gsea(
        x = epheno,
        gsets = mySignature,
        B = 1000,
        p.adjust = 'BH',
        center = TRUE,
        pval.comp.method = "original"
    )

A customized GSEA does not run fine:

my.gsea2 <-
    gsea(
        x = epheno,
        gsets = mySignature,
        B = 1000,
        p.adjust = 'BH',
        pval.comp.method = "signed" # [KEY CHANGE]
    )

# Error in getNesPval(nes, pval.comp.method, pval.smooth.tail) : 
#     object 'escore' not found

I've commented the test code above to highlight the key parameter change that reveals the bug. If I follow the traceback, I land here:

traceback()

# 13: getNesPval(nes, pval.comp.method, pval.smooth.tail)
# 12: gseaSignificance(x, y, pval.comp.method, pval.smooth.tail)
# 11: (function (x, y, z) 
#     gseaSignificance(x, y, pval.comp.method, pval.smooth.tail))(x = dots[[1L]][[1L]], 
#         y = dots[[2L]][[1L]])
# 10: mapply(function(x, y, z) gseaSignificance(x, y, pval.comp.method, 
#         pval.smooth.tail), x = es, y = es.sim)
# 9: getSummary(es, es.sim, fchr, p.adjust.method = p.adjust.method, 
#        pval.comp.method, pval.smooth.tail, signatures, test, fewGsets)
# 8: cbind(n = unlist(lapply(x$s, length)), getSummary(es, es.sim, 
#        fchr, p.adjust.method = p.adjust.method, pval.comp.method, 
#        pval.smooth.tail, signatures, test, fewGsets))
# 7: gseaSignificance(x, p.adjust.method, pval.comp.method, pval.smooth.tail)
# 6: gseaSignificance(x, p.adjust.method, pval.comp.method, pval.smooth.tail)
# 5: FUN(X[[i]], ...)
# 4: lapply(x, function(x) gseaSignificance(x, p.adjust.method, pval.comp.method, 
#        pval.smooth.tail))
# 3: gseaSignificance(sim, p.adjust.method = p.adjust.method, pval.comp.method = pval.comp.method, 
#        pval.smooth.tail = pval.smooth.tail)
# 2: gseaSignificance(sim, p.adjust.method = p.adjust.method, pval.comp.method = pval.comp.method, 
#        pval.smooth.tail = pval.smooth.tail)
# 1: gsea(x = epheno, gsets = mySignature, B = 1000, p.adjust = "BH", 
#        pval.comp.method = "signed")

Possible Solution

Here is a link to the line in the current package source where the traceback stops (https://github.com/Bioconductor-mirror/phenoTest/blob/release-3.5/R/gsea.R#L286).

The internal "getNesPval" function itself is defined above (https://github.com/Bioconductor-mirror/phenoTest/blob/release-3.5/R/gsea.R#L194)

It appears that the "escore" variable needed deep inside the "getNesPval" function is not defined there (https://github.com/Bioconductor-mirror/phenoTest/blob/release-3.5/R/gsea.R#L213).

While in debug mode, if I manually redefined the "getNesPval" function to include "escore" as a parameter, the top-level "gsea" command was able to finish successfully.

Session Info

session_info()

# Session info -------------------------------------------------------------------------------------------
#  setting  value                       
#  version  R version 3.4.0 (2017-04-21)
#  system   x86_64, linux-gnu           
#  ui       RStudio (1.0.143)           
#  language (EN)                        
#  collate  en_US.UTF-8                 
#  tz       <NA>                        
#  date     2017-08-31                  
# 
# Packages -----------------------------------------------------------------------------------------------
#  package              * version  date       source        
#  acepack                1.4.1    2016-10-29 CRAN (R 3.4.0)
#  affy                   1.54.0   2017-05-02 Bioconductor  
#  affyio                 1.46.0   2017-05-02 Bioconductor  
#  annotate             * 1.54.0   2017-05-02 Bioconductor  
#  AnnotationDbi        * 1.38.2   2017-08-16 Bioconductor  
#  backports              1.1.0    2017-05-22 CRAN (R 3.4.0)
#  base                 * 3.4.0    2017-05-02 local         
#  base64enc              0.1-3    2015-07-28 CRAN (R 3.4.0)
#  Biobase              * 2.36.2   2017-07-19 Bioconductor  
#  BiocGenerics         * 0.22.0   2017-05-02 Bioconductor  
#  BiocInstaller          1.26.0   2017-05-02 Bioconductor  
#  biomaRt                2.32.1   2017-07-19 Bioconductor  
#  BioNet                 1.36.0   2017-07-18 Bioconductor  
#  Biostrings             2.44.2   2017-07-25 Bioconductor  
#  bit                    1.1-12   2014-04-09 CRAN (R 3.4.0)
#  bit64                  0.9-7    2017-05-08 CRAN (R 3.4.0)
#  bitops                 1.0-6    2013-08-17 CRAN (R 3.4.0)
#  blob                   1.1.0    2017-06-17 CRAN (R 3.4.0)
#  BMA                  * 3.18.7   2017-05-12 CRAN (R 3.4.0)
#  Category               2.42.1   2017-07-18 Bioconductor  
#  caTools                1.17.1   2014-09-10 CRAN (R 3.4.0)
#  cellHTS2               2.40.0   2017-07-18 Bioconductor  
#  checkmate              1.8.3    2017-07-03 CRAN (R 3.4.0)
#  cluster                2.0.6    2017-03-16 CRAN (R 3.4.0)
#  codetools              0.2-15   2016-10-05 CRAN (R 3.4.0)
#  colorspace             1.3-2    2016-12-14 CRAN (R 3.4.0)
#  compiler               3.4.0    2017-05-02 local         
#  data.table             1.10.4   2017-02-01 CRAN (R 3.4.0)
#  datasets             * 3.4.0    2017-05-02 local         
#  DBI                    0.7      2017-06-18 CRAN (R 3.4.0)
#  DelayedArray           0.2.7    2017-07-19 Bioconductor  
#  DEoptimR               1.0-8    2016-11-19 CRAN (R 3.4.0)
#  devtools             * 1.13.3   2017-08-02 CRAN (R 3.4.0)
#  digest                 0.6.12   2017-01-27 CRAN (R 3.4.0)
#  ellipse                0.3-8    2013-04-13 CRAN (R 3.4.0)
#  ff                     2.2-13   2014-04-09 CRAN (R 3.4.0)
#  foreach                1.4.3    2015-10-13 CRAN (R 3.4.0)
#  foreign                0.8-69   2017-06-21 CRAN (R 3.4.0)
#  Formula                1.2-2    2017-07-10 CRAN (R 3.4.0)
#  gdata                  2.18.0   2017-06-06 CRAN (R 3.4.0)
#  genefilter             1.58.1   2017-07-19 Bioconductor  
#  GenomeInfoDb           1.12.2   2017-07-19 Bioconductor  
#  GenomeInfoDbData       0.99.0   2017-05-02 Bioconductor  
#  GenomicRanges          1.28.4   2017-07-19 Bioconductor  
#  ggplot2              * 2.2.1    2016-12-30 CRAN (R 3.4.0)
#  gmp                    0.5-13.1 2017-03-10 CRAN (R 3.4.0)
#  gplots                 3.0.1    2016-03-30 CRAN (R 3.4.0)
#  graph                * 1.54.0   2017-05-02 Bioconductor  
#  graphics             * 3.4.0    2017-05-02 local         
#  grDevices            * 3.4.0    2017-05-02 local         
#  grid                   3.4.0    2017-05-02 local         
#  gridExtra              2.2.1    2016-02-29 CRAN (R 3.4.0)
#  GSEABase             * 1.38.0   2017-07-18 Bioconductor  
#  gtable                 0.2.0    2016-02-26 CRAN (R 3.4.0)
#  gtools                 3.5.0    2015-05-29 CRAN (R 3.4.0)
#  Heatplus             * 2.22.0   2017-05-02 Bioconductor  
#  hgu133a.db             3.2.3    2017-07-18 Bioconductor  
#  Hmisc                  4.0-3    2017-05-02 CRAN (R 3.4.0)
#  hopach                 2.36.0   2017-07-18 Bioconductor  
#  htmlTable              1.9      2017-01-26 CRAN (R 3.4.0)
#  htmltools              0.3.6    2017-04-28 CRAN (R 3.4.0)
#  htmlwidgets            0.9      2017-07-10 CRAN (R 3.4.0)
#  HTSanalyzeR            2.28.0   2017-08-30 Bioconductor  
#  hwriter                1.3.2    2014-09-10 CRAN (R 3.4.0)
#  igraph                 1.1.2    2017-07-21 CRAN (R 3.4.0)
#  inline               * 0.3.14   2015-04-13 CRAN (R 3.4.0)
#  IRanges              * 2.10.2   2017-07-19 Bioconductor  
#  iterators              1.0.8    2015-10-13 CRAN (R 3.4.0)
#  KernSmooth             2.23-15  2015-06-29 CRAN (R 3.4.0)
#  knitr                  1.17     2017-08-10 CRAN (R 3.4.0)
#  lattice                0.20-35  2017-03-25 CRAN (R 3.4.0)
#  latticeExtra           0.6-28   2016-02-09 CRAN (R 3.4.0)
#  lazyeval               0.2.0    2016-06-12 CRAN (R 3.4.0)
#  leaps                * 3.0      2017-01-10 CRAN (R 3.4.0)
#  limma                  3.32.5   2017-08-16 Bioconductor  
#  locfit                 1.5-9.1  2013-04-20 CRAN (R 3.4.0)
#  magrittr               1.5      2014-11-22 CRAN (R 3.4.0)
#  MASS                   7.3-47   2017-04-21 CRAN (R 3.4.0)
#  Matrix                 1.2-11   2017-08-16 CRAN (R 3.4.0)
#  matrixStats            0.52.2   2017-04-14 CRAN (R 3.4.0)
#  memoise                1.1.0    2017-04-21 CRAN (R 3.4.0)
#  methods              * 3.4.0    2017-05-02 local         
#  mgcv                   1.8-18   2017-07-28 CRAN (R 3.4.0)
#  munsell                0.4.3    2016-02-13 CRAN (R 3.4.0)
#  mvtnorm                1.0-6    2017-03-02 CRAN (R 3.4.0)
#  nlme                   3.1-131  2017-02-06 CRAN (R 3.4.0)
#  nnet                   7.3-12   2016-02-02 CRAN (R 3.4.0)
#  oligoClasses           1.38.0   2017-05-03 Bioconductor  
#  parallel             * 3.4.0    2017-05-02 local         
#  pcaPP                  1.9-72   2017-06-27 CRAN (R 3.4.0)
#  phenoTest            * 1.24.0   2017-08-30 Bioconductor  
#  pkgconfig              2.0.1    2017-03-21 CRAN (R 3.4.0)
#  plyr                   1.8.4    2016-06-08 CRAN (R 3.4.0)
#  prada                  1.52.0   2017-07-18 Bioconductor  
#  preprocessCore         1.38.1   2017-07-19 Bioconductor  
#  RankProd               3.2.0    2017-08-30 Bioconductor  
#  RBGL                   1.52.0   2017-05-03 Bioconductor  
#  RColorBrewer           1.1-2    2014-12-07 CRAN (R 3.4.0)
#  Rcpp                   0.12.12  2017-07-15 CRAN (R 3.4.0)
#  RCurl                  1.95-4.8 2016-03-01 CRAN (R 3.4.0)
#  rlang                  0.1.2    2017-08-09 CRAN (R 3.4.0)
#  Rmpfr                  0.6-1    2016-11-15 CRAN (R 3.4.0)
#  robustbase           * 0.92-7   2016-12-09 CRAN (R 3.4.0)
#  rpart                  4.1-11   2017-04-21 CRAN (R 3.4.0)
#  rrcov                * 1.4-3    2016-09-06 CRAN (R 3.4.0)
#  RSQLite                2.0      2017-06-19 CRAN (R 3.4.0)
#  rstudioapi             0.6      2016-06-27 CRAN (R 3.4.0)
#  S4Vectors            * 0.14.3   2017-07-19 Bioconductor  
#  scales                 0.4.1    2016-11-09 CRAN (R 3.4.0)
#  SNPchip                2.22.0   2017-07-18 Bioconductor  
#  splines                3.4.0    2017-05-02 local         
#  splots                 1.42.0   2017-07-18 Bioconductor  
#  stats                * 3.4.0    2017-05-02 local         
#  stats4               * 3.4.0    2017-05-02 local         
#  stringi                1.1.5    2017-04-07 CRAN (R 3.4.0)
#  stringr                1.2.0    2017-02-18 CRAN (R 3.4.0)
#  SummarizedExperiment   1.6.3    2017-07-19 Bioconductor  
#  survival             * 2.41-3   2017-04-04 CRAN (R 3.4.0)
#  tibble                 1.3.3    2017-05-28 CRAN (R 3.4.0)
#  tools                  3.4.0    2017-05-02 local         
#  utils                * 3.4.0    2017-05-02 local         
#  vsn                    3.44.0   2017-07-18 Bioconductor  
#  withr                  2.0.0    2017-07-28 CRAN (R 3.4.0)
#  XML                  * 3.98-1.9 2017-06-19 CRAN (R 3.4.0)
#  xtable                 1.8-2    2016-02-05 CRAN (R 3.4.0)
#  XVector                0.16.0   2017-05-03 Bioconductor  
#  yaml                   2.1.14   2016-11-12 CRAN (R 3.4.0)
#  zlibbioc               1.22.0   2017-05-03 Bioconductor
phenoTest • 1.2k views
ADD COMMENT

Login before adding your answer.

Traffic: 427 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6