bug (?) in selectHKgenes() related to geomMean() in SLqPCR package (with fix)
0
1
Entering edit mode
Davis, Wade ▴ 350
@davis-wade-2803
Last seen 7.4 years ago
Dear Dr. Kohl: First, thank you for your work on the SLqPCR package. I've found it very useful over the years. For me, one of the most useful functions in SLqPCR is selectHKgenes(). Recently when calling selectHKgenes(), I came across an error that occurs when NA values are found in the relData matrix, which traces back to a bug (I think) in geomMean. Even though the option na.rm=T is used, it produces an error when encountering NAs. Below is a reproducible example with a proposed fix. Sincerely, Wade Davis relData<-matrix(data= rnorm(n=20,mean=1,sd=0.2),nr=4,ncol=5) relData [,1] [,2] [,3] [,4] [,5] [1,] 0.8083619 1.172125 0.8412847 1.0004588 0.5717653 [2,] 1.0717140 1.123366 1.3911562 0.9947363 0.7826146 [3,] 1.1598977 1.161934 1.0233950 0.9287549 1.1451843 [4,] 1.0506174 0.982740 1.0509895 0.8127318 1.0421654 # behaves as expected selectHKgenes(relData=relData,geneSymbol=letters[1:5],na.rm=F) selectHKgenes(relData=relData,geneSymbol=letters[1:5],na.rm=T) # introduce NA relData[1,1]<-NA relData [,1] [,2] [,3] [,4] [,5] [1,] NA 1.029276 1.0139564 1.0834577 0.4496888 [2,] 0.7939039 1.245437 0.9915528 0.9748007 1.5730212 [3,] 0.7362098 1.400378 1.0634846 1.1232321 1.1641223 [4,] 1.0944912 1.017443 0.8877103 0.7744538 0.8976991 # error selectHKgenes(relData=relData,geneSymbol=letters[1:5],na.rm=T) Error in if (any(x < 0)) stop("'x' contains negative value(s)") : missing value where TRUE/FALSE needed #look at geomMean function # compute geometric mean > geomMean(c(NA,1:10),na.rm=TRUE) Error in if (any(x < 0)) stop("'x' contains negative value(s)") : missing value where TRUE/FALSE needed ################################ # fix broken function ################################ geomMean <- function(x, na.rm = FALSE){ if(!is.numeric(x) && !is.complex(x) && !is.logical(x)){ warning("argument is not numeric or logical: returning NA") return(as.numeric(NA)) } if(na.rm) x <- x[!is.na(x)] #### fix is this line ######## if(any(x < 0)) stop("'x' contains negative value(s)") return(prod(x)^(1/length(x))) } assignInNamespace("geomMean", geomMean, "SLqPCR") # try again with fixed function, both functions work as expected geomMean(c(NA,1:10),na.rm=TRUE) [1] 4.528729 selectHKgenes(relData=relData,geneSymbol=letters[1:5],na.rm=T) ############################################################### Step 1 : gene expression stability values M: c b d a e 0.3770762 0.3884963 0.4402875 0.5500149 0.7280952 average expression stability M: 0.496794 gene with lowest stability (largest M value): e Pairwise variation, ( 4 / 5 ): 0.1870483 ############################################################### Step 2 : gene expression stability values M: c b d a 0.2442868 0.3055493 0.3075325 0.5130044 average expression stability M: 0.3425932 gene with lowest stability (largest M value): a Pairwise variation, ( 3 / 4 ): 0.1100028 ############################################################### Step 3 : gene expression stability values M: c d b 0.1498885 0.1755862 0.1910716 average expression stability M: 0.1721821 gene with lowest stability (largest M value): b Pairwise variation, ( 2 / 3 ): 0.06023353 ############################################################### Step 4 : gene expression stability values M: d c 0.1344031 0.1344031 average expression stability M: 0.1344031 $ranking 1 1 3 4 5 "c" "d" "b" "a" "e" $variation 4/5 3/4 2/3 0.18704830 0.11000284 0.06023353 $meanM 5 4 3 2 0.4967940 0.3425932 0.1721821 0.1344031 [[alternative HTML version deleted]]
SLqPCR SLqPCR • 623 views
ADD COMMENT

Login before adding your answer.

Traffic: 293 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