Entering edit mode
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
