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

