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]]