R/plotPOD.R

In markusboenn/POD: Probability of Detection for Qualitative PCR Methods

#'
#' @title Generate Plot to Analyze Single Lab PCR Outcomes
#' @description Show POD curve and LOD value to validate qualitative PCR methods of a single laboratory.
#'
#' @param obj A list returned by \code{\link{analyzeSingleLab}}.
#' @param model Simple or full model
#' @param qLOD The quantile(s) for LOD to be shown in the plot. Multiplied by \eqn{100} if less than one.
#' @param show.ci Show the confidence interval of the LOD in the plot.
#' @param show.warnings Show the warning regarding significant deviation from \eqn{1} in the plot.
#' @param wmark Logical. Show a watermark at the upper right corner of the plot.
#' @param unit A string indicating the unit of the data.
#'
#' @details The graph generated by this function gives the laboratory-specific rates of detection (RODs) as blue diamonds. The blue curve denotes the mean POD curve along with the corresponding \eqn{95 \%} confidence range highlighted as the grey band. The POD curve under ideal conditions is displayed as the black dashed curve.
#'
#' If model is set to "auto", a plausiblity test is applied to determine if the POD curve bases on the simplified or on full parameter estimation. If the corrective parameter determined from the full model significantly differs from \eqn{1}, a message is shown in the plot. Testing for significant deviation is currently done by checking the condition \eqn{1-b>0.2}. The threshould \eqn{0.2} has been determined empirically to agree with the original webtool and might be changed in future versions of the package.
#'
#' Three cases can be distinguished. First, the value for the slope parameter b is significantly less than \eqn{1}. This means the average amplification probability is higher at higher dilution levels than at lower dilution levels. Such a situation can be related to: inhibitory matrix effects, a large variability in the amplification process from the one test to another under repeatability conditions, or accidental problems causing false positives if the number of copies of the target DNA sequence is less than \eqn{1}. Second, the calculated POD curve indicates sensitivity better than achievable according to the theoretical POD curve. Third, the number of positive test results is significantly higher than expected at nominal copies of nominal DNA concentrations in \eqn{[0.5,1.5]}. In this case check the correctness of the serial dilution.
#'
#' Another warning appears if the LOD of interest exceeds the highest number of considered nominal copies.
#'
#' The unit is add to the LOD value, in front of the confidence intervall.
#'
#' @return The passed list 'obj' is returned invisibly.
#'
#' @export
#' @rdname plotPOD
#' @importFrom graphics abline axis lines mtext plot points polygon segments text title
#' @importFrom grDevices adjustcolor
#' @importFrom stats ks.test
#' @importFrom utils packageVersion
#'
#' @examples
#' x <- cbind(
#'  X=c(0.1,1,2,5,10,20),
#'  S=c( 0,5,6,6,6,6 ),
#'  N=c( 6,6,6,6,6,6 )
#' )
#' obj <- analyzeSingleLab(x=x)
#' plotPOD(obj)

#C Parameter 'app' controls the appearance of the figure. The default is "quodata" and aims to generate a plot looking similar to the one generated by the original webtool of Quodata (<http://quodata.de/content/validation-qualitative-pcr-methods-single-laboratory/>).

plotPOD <- function(obj, model=c("auto", "simple", "full"), qLOD=95, show.ci=TRUE, show.warnings=FALSE, wmark=TRUE, unit=""){
    if(class(obj) != "pod"){
        warning("Variable 'obj' is not of class 'pod'. Did you use 'analyzeSingleLab' to generate 'obj'?")
        return(obj)
    }
    x <- obj$x
    XLIM <- range(x$X);
    YTICK <- seq(0, 1, 0.05)
    YTICKLAB <- as.character(YTICK)
    if( qLOD[1] < 1 ){qLOD <- qLOD*100}
    qLOD <- sort(qLOD)
    
    XLIM <- range(c(XLIM, 100))
    YTICKLAB <- sapply(YTICKLAB, function(v){ switch(nchar(v), paste0(v, ".00"), v, paste0(v, "0"), v) } )

    xsmooth <- .smoothLog( x$X )
    WARN1 <- WARN2 <- WARN3 <- WARN4 <- ""
    model <- model[1]

    
    .b <- obj$b
    bfull <- obj$fit.glm.full$b
    bdiff <- .b-bfull
    tol <- 0.2; # the threshould 0.2 has been determined empirically and might differ from the threshould used in the webtool
    
    FIT <- obj[["fit.glm.simple"]]
    if(model=="full"){ FIT <- obj[["fit.glm.full"]]; }
    if(model=="auto"){
        if( bdiff > tol ){
            cat("Using full GLM\n")
            FIT <- obj[["fit.glm.full"]]
            WARN1 <- paste0("* The value for the slope parameter b is significantly less than ", .b, ". This means the average amplification probability is higher at higher dilution levels than at lower dilution levels.")
        }
        #cat(" ( ", .b, " vs. ", bfull, " [", bdiff, "] )\n", sep="")
    }
    
    pd <- FIT$cip95
    lambda <- FIT$lambda
    b <- FIT$b
    LodList <- FIT$LOD
    if( ! any(as.character(qLOD) %in% rownames(LodList)) ){
        warning( "Provided 'qLOD' values not available. Using first quantile (", rownames(LodList)[1], ")" )
        qLOD <- rownames(LodList)[1]
    }
    if( LodList[as.character(qLOD), nrow(LodList)][1] > max(x$X) ){
        WARN4 <- "* The LOD exceeds the highest number of considered nominal copies. Hence, the validation criterion is not fulfilled."
    }
    
    # compute POD curves
    PodIdeal <- computePOD(xsmooth, 1, 1)
    PodFit <- computePOD(xsmooth, lambda, b)
    
    # compute difference between ideal and fitted curve
    # suppress warnings arrising due to ties
    ks <- suppressWarnings(ks.test(PodIdeal, PodFit, alternative="less"))
    if(ks$statistic == 0){
        WARN2 <- "* The calculated POD curve indicates sensitivity better than achievable according to the theoretical POD curve."
        ks2 <- suppressWarnings(ks.test(PodIdeal, PodFit, alternative="greater"))
        if( ks2$p.value < 0.05 ){ 
            WARN2 <- paste0(WARN2, " It is most likely that the given nominal copies are incorrect. Please check it again.")
        }
    }
    obj$ks <- ks
 
    
    x.warn3 <- which( x$X >= 0.5 & x$X <= 1.5 )
    if( length(x.warn3) ){
        xx.warn3 <- round(x$X[x.warn3], 2)
        candgt <- candlt <- c()
        #DEBUG print(pd[ pd$X %in% xx.warn3, ]); print(x[x.warn3,])
        for( i in 1:length( xx.warn3 ) ){
            if( round(pd[ pd$X == xx.warn3[i], "Upper" ], 2) < x[x.warn3[i],]$Y ){ candgt <- c(candgt, x$X[x.warn3][i]) }
            if( pd[ pd$X == xx.warn3[i], "Lower" ] > x[x.warn3[i],]$Y ){ candlt <- c(candlt, x$X[x.warn3][i]) }
        }
        if(length(candgt)){
            candgt <- paste(collapse=", ", candgt)
            WARN3 <- paste0("* The number of positive test results is significantly higher than expected at nominal copies of ", paste(collapse=", ", candgt), ". Please check the correctness of the serial dilution.")
        }
    }
    WARNALL <- gsub(paste(collapse="\n", c(WARN1, WARN2, WARN3, WARN4)), pattern="[\n]{1,}", replacement="\n")
    WARNALL <- gsub(WARNALL, pattern="^\n|\n$", replacement="")
    
    # actual plot
    plot( x=1,y=1, type="n", ylim=c(0,1), xlim=XLIM, log="x", ylab="POD and ROD", xlab="Number of DNA copies", yaxt="n" )
    axis(2, at=YTICK, labels=YTICKLAB, las=2)
    if( show.warnings && nchar(WARNALL) ){ title( sub="Warnings appeared. Check command line output.", cex.sub=1.5, font.sub=2, col.sub=2 ) }

    # add helper lines
    abline(h=seq(0,1,0.05), lwd=0.5, col="gray80")
    
    # compute and add 95 % confidence range
    lines(x=pd$X, y=pd$Upper, col=.COLOR.SHADE)
    lines(x=pd$X, y=pd$Lower, col=.COLOR.SHADE)
    polygon(c(pd$X, rev(pd$X), pd$X[1]), c(pd$Upper, rev(pd$Lower), pd$Upper[1]), col = adjustcolor(.COLOR.SHADE, alpha.f = 0.5), border = NA)

    # add observed data
    points(x=x$X, y=x$Y, pch=18, col=.COLOR.OBS, cex=2)

    # add ideal POD function
    lines( x=xsmooth, y=PodIdeal, col="black", lty=2 )

    # add fitted POD function
    lines( x=xsmooth, y=PodFit, col=.COLOR.OBS, lty=1 )

    
    # add LOD
    unit <- gsub( unit, pattern="^ *| *$", replacement="" ) # drop leading and trailing spaces
    if(nchar(unit)){ unit <- paste0(" ", unit) }
    for( nLOD in as.character(qLOD) ){
        LOD <- LodList[nLOD, ]
        nLOD <- as.numeric(nLOD)
        LODr <- round(LOD[1], 2)
        LODll <- round(LOD[2], 3)
        LODul <- round(LOD[3], 3)
        PodFit <- computePOD(LOD[1], lambda, b)
        if( LOD[1] > 1e-3 ){
            segments(LOD[1], 0, LOD[1], PodFit,     col=.COLOR.LOD, lty=2, lwd=4)
            segments(min(x$X), nLOD/100, LOD[1], nLOD/100, col=.COLOR.LOD, lty=2, lwd=4)
            segments(LOD[2], 0, LOD[3], 0,   col=.COLOR.LOD, lty=1, lwd=4)
            lab <- bquote(paste(LOD[.(nLOD)*'%']== .(LODr)))
            if( length(qLOD)==1 && show.ci ){
                lab <- bquote( paste( sep="", LOD[.(nLOD)*'%']==.(LODr), .(unit), " [", .(LODll), ";", .(LODul), "]" ) )
            }
            text(x=LOD[1], y=0.05, pos=4, labels=lab, col=.COLOR.LOD)
        }
    }
    
    # digital watermark
    if( wmark ){
        mtext(text=paste0("POD ", packageVersion("POD")), side=3, line=0, adj=1, col="gray", cex=0.75)
    }
    
    obj$warn <- c(WARN1, WARN2, WARN3, WARN4)
    if( nchar(WARNALL) ){
        cat("The following warning(s) appeared:\n")
        cat(WARNALL, "\n", sep="")
    }
    
    
 
return(invisible(obj))
}