---
title: "Algorithms for Automatized Detection of Hook Effect-bearing Amplification Curves"
author: "The PCRedux package authors"
date: "`r Sys.Date()`"
output: 
  rmarkdown::pdf_document:
    number_sections: true
    toc: true
bibliography: "literature.bib"
vignette: >
  %\VignetteIndexEntry{Algorithms for Automatized Detection of Hook Effect-bearing Amplification Curves}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

  \begin{figure}[ht]
  \centering
  \scalebox{0.6}{
  \includegraphics[clip=true,trim=1cm 1cm 1cm 1cm]{Logo.pdf}
  }
  \end{figure}
  
  The PDF version of this document is available [online](https://github.com/devSJR/PCRedux/raw/master/docs/articles/SI1.pdf).
  
  **The document is part of the publication *M. Burdukiewicz, A.-N. Spiess, K.A. Blagodatskikh, W. Lehmann, P. Schierack, S. Rödiger, Algorithms for automated detection of hook effect-bearing amplification curves, Biomolecular Detection and Quantification. (2018). doi:10.1016/j.bdq.2018.08.001.***

  
# Abstract
  This is a supplemental document for the study *Algorithms for Automatized Detection of Hook Effect-bearing Amplification Curves*.
  Quantitative real-time PCR (qPCR) is a widely used method for gene expression analysis, forensics and medical diagnostics [@dvinge_htqpcr:_2009; @martins_dna_2015; @sauer_differentiation_2016].
  
  Numerous algorithms have been developed to extract features from amplification curves such as the cycle of quantification and the amplification efficiency [@ruijter_evaluation_2013]. There is an agreement, that these algorithms need to be evaluated and benchmarked for their performance [@kemperman_mircomp-shiny:_2017]. But at an earlier level it is important to have a solid foundation for the data preprocessing [@spiess_impact_2015; @spiess_system-specific_2016; @ronde_practical_2017].
  Digitalization of processes holds the promise that potential human mistakes can be spotted and that diagnostic processes can be automatized.
  
  The aim of the study is to provide software tools and algorithms, which assists qPCR users during the analysis and quality management of their data. In particular, this study shows how it is possible to automatically detect hook effects (see @barratt_improving_2002) or hook effect-like curvatures. 
  
  
# Introduction
  
  The functions and data 
  presented in the paper are available from <https://github.com/devSJR/PCRedux>. 
  The data, including the RDML file, are part of the \texttt{PCRedux} package and are 
  made available in the CSV or RDML format [@roediger_enabling_2017] for vendors independent analysis.
  
  All analyses were implemented and conducted with the **R** statistical computing 
  language [@R_language; @roediger2015r] and dedicated integrated development 
  environments such as **RKWard** [@roediger_rkward_2012]. Further documentation 
  can be found in the help files of the **R** packages.
  
  
  
# Installation
  
  The **hookreg()** and **hookregNL()** functions are part of the \texttt{PCRedux} package for the **R** statistical computing language.  Download from CRAN <http://cran.r-project.org/> the **R** version for the required operating system and install **R**. Then start **R** and type in the prompt:
  
  ```{R, eval=FALSE, echo=TRUE}
# Select your local mirror
install.packages("PCRedux")
```

The \texttt{PCRedux} package should just install. If this fails make sure you have write access to the destination directory and follow the instructions of the **R** documentation:

```{R, eval=FALSE, echo=TRUE}
# The following command points to the help for download and install of packages 
# from CRAN-like repositories or from local files.
?install.packages()
```

The package can be installed as the latest development version using the [devtools](https://cran.r-project.org/package=devtools) R package.

```{R, eval=FALSE, echo=TRUE}
# Install devtools, if you haven't already.
install.packages("devtools")

library(devtools)
install_github("devSJR/PCRedux")
```

It is recommended to use software with an integrated development environment such as \texttt{RKWard} [@roediger_rkward_2012]. To work with RDML data it is recommend to use the \textit{RDML} package ($\geq$~v.~0.9-9) by invoking the **rdmlEdit()** function (for details see @roediger_enabling_2017) or the rdmlEdit GUI web server (section \ref{rdmlEditserver}). The RDML file `hookreg.rdml` contains the amplification curve data. However, other software package (e.g., [@lefever_rdml_2009; @ruijter_rdml-ninja_2015]) can also be used to work with the RMDL data file format.

# Results for the analysis of the `hookreg.rdml` data set by **humanrater()**

All calculations in the following sections were employed on the `hookreg.rdml` data 
**R** environment by the \texttt{RDML} package [@roediger_enabling_2017]. An overview of 
the used samples and the qPCR detection chemistries and the classification by two humans ("Hook effect-like Rater 1", "Hook effect-like Rater 2") is shown in Table \ref{Table_human_rated}.

```{r eval=TRUE, echo=FALSE, results="asis"}
# Load the RDML package for reading of the hookreg.rdml-file from the PCRedux
# package.

library(RDML)
library(PCRedux)

# A comprehensive description of the RDML-file import can be found in Roediger 
# et al. (2017) Bioinformatics

raw_data <- RDML$new(filename=system.file("hookreg.rdml", package="PCRedux"))

raw_data_tab <- raw_data$AsTable(name.pattern=paste(
  react$position,
  react$sample$id,
  # run id added to names
  sep="~"))
data <- as.data.frame(raw_data$GetFData(raw_data_tab, long.table=FALSE))
```

All amplification curves were plotted according to their experiment conditions. They differed in the target molecules (e.g., \textit{MLC-2v}, \textit{BRCA1}) and the detection chemistries (e.g., EvaGreen, SybrGreen, hydrolysis probes). Figure \ref{ampcurves} shows seven plots for the corresponding experiments. The amplification curves were not preprocessed to preserve the curvature. Selected amplification curves were noisy (e.g., Figure \ref{ampcurves}F), had overshots or undershot in the background phase (e.g., Figure \ref{ampcurves}E-G), a short hook phase (e.g., Figure \ref{ampcurves}D). Amplification curves of Figure \ref{ampcurves}A, D, F and F exhibited a clearly visible hook effect or a hook like effect.

```{r ampcurves,fig.width=7.25, fig.height=9.25, echo=TRUE, message=FALSE, warning=FALSE, fig.cap="Amplification curves. A) Synthetic template, detected with Syto-13. B) Human \\textit{MLC-2v}, detected with a hydrolysis probe. C) \\textit{S27a} housekeeping gene, detected with SybrGreen I. D) Whole genome amplification, detected with EvaGreen. E) Human \\textit{BRCA1} gene, detected with a hydrolysis probe. F) Human \\textit{NRAS} gene, detected with a hydrolysis probe. G) Water control, detected with a hydrolysis probe. See Table \\ref{Table_human_rated} for details. RFU, relative fluorescence units. \\label{ampcurves}"}
par(mfrow=c(4,2))

# Plot all data of the hookreg.rdml-file according to their type.
# Synthetic template, detected with Syto-13
matplot(data[, 1], data[, 2:13], type="l", lty=1, lwd=2, ylab="RFU", xlab="Cycle")
mtext("A", cex = 1.8, side = 3, adj = 0, font = 2)

# Human MLC-2v, detected with a hydrolysis probe.
matplot(data[, 1], data[, 14:45], type="l", lty=1, lwd=2, ylab="RFU", xlab="Cycle")
mtext("B", cex = 1.8, side = 3, adj = 0, font = 2)

# S27a housekeeping gene, detected with SybrGreen I.
matplot(data[, 1], data[, 46:69], type="l", lty=1, lwd=2, ylab="RFU", xlab="Cycle")
mtext("C", cex = 1.8, side = 3, adj = 0, font = 2)

# Whole genome amplification, detected with EvaGreen.
matplot(data[, 1], data[, 70:71], type="l", lty=1, lwd=2, ylab="RFU", xlab="Cycle")
mtext("D", cex = 1.8, side = 3, adj = 0, font = 2)

# Human BRCA1 gene, detected with a hydrolysis probe.
matplot(data[, 1], data[, 72:87], type="l", lty=1, lwd=2, ylab="RFU", xlab="Cycle")
mtext("E", cex = 1.8, side = 3, adj = 0, font = 2)

# Human NRAS gene, detected with a hydrolysis probe.
matplot(data[, 1], data[, 88:95], type="l", lty=1, lwd=2, ylab="RFU", xlab="Cycle")
mtext("F", cex = 1.8, side = 3, adj = 0, font = 2)

# Water control, detected with a hydrolysis probe.
matplot(data[, 1], data[, 96:97], type="l", lty=1, lwd=2, ylab="RFU", xlab="Cycle")
mtext("G", cex = 1.8, side = 3, adj = 0, font = 2)
```

Printout of all measured samples, their rating by two humans (rater 1 and rater 
2) with their dichotomous ratings (0, no hook; 1, hook) and their sources.

- The boggy data (qpcR::boggy) set was taken from the \texttt{qpcR} package 
[@Ritz2008; @spiess_highly_2008].
- The C127EGHP data (chipPCR::C127EGHP) set was taken from the \texttt{chipPCR} package 
[@roediger2015chippcr].
- The testdat data (qpcR::testdat) set was taken from the \texttt{qpcR} package 
[@Ritz2008; @spiess_highly_2008].
- Other data were prepared by Evrogen laboratory experiments.


```{r eval=TRUE, echo=FALSE, results="asis"}
library(readxl)
library(xtable)
options(xtable.comment=FALSE)

Table_human_rated <- read_xlsx(path=system.file("Table_human_rated.xlsx", 
                                                package="PCRedux"))

print(xtable(Table_human_rated, digits=0, 
             caption = "Overview of the used amplification curve data. The 
samples names, data source (origin of data either from an existing data set or  
prepared for this study), the detection chemistries (intercalator (Syto-13, 
SyberGreenI, EvaGreen), hydrolysis probes (TaqMan (Cy5/BHQ2) , TaqMan 
(HEX/BHQ1))) and calculations by tow humans.", 
label='Table_human_rated'),
size = "\\tiny",
include.rownames = FALSE,
include.colnames = TRUE,
caption.placement = "top",
comment=FALSE,
table.placement = "!ht", scalebox='0.65'
)
```

# Results for the analysis with **hookreg()** and **hookregNL()**

This section contains the results of the analysis of the amplification curve 
data with the **hookreg()** function and the **hookregNL()** function. As in the 
previous sections, all code was commented to make it reproducible. Some rows in 
Table \ref{res_hookreg_table} and Table \ref{res_hookregNL_table} appear to be 
empty. This expected behavior may occur in cases where the corresponding 
functions were not able to calculate the coefficients due to a failed model fit 
or violation of the truncation criterion.

## Results for the analysis of the `hookreg.rdml` data set with **hookreg()** \label{sectionHookreg}

The following code was used to analyze the `hookreg.rdml` data set with **hookreg()** function. The **hookreg()** function fits a linear model to a region of interest. The linear model is used to decide if the amplification curve as a hook effect or hook effect-like curvature.

```{r, echo=TRUE, message=FALSE, warning=FALSE}
# Load PCRedux package to obtain the data and make the hookreg() function
# available.
library(PCRedux)

# `data` is a temporary data frame of the hook.rdml amplification curve data file.
# Apply the hookreg() function over the amplification curves and arrange the 
# results in the data frame `res_hookreg`.

res_hookreg <-  data.frame(sample=colnames(data)[-1], 
                           t(sapply(2L:ncol(data), function(i) {
                             hookreg(x=data[, 1], y=data[, i])
                           })))

# Fetch the calculated parameters from the calculations with the hookreg() 
# function as a table `res_hookreg_table`.

res_hookreg_table <- data.frame(sample=as.character(res_hookreg[["sample"]]),
                                intercept=signif(res_hookreg[["intercept"]], 2),
                                slope=signif(res_hookreg[["slope"]], 1),
                                hook.start=signif(res_hookreg[["hook.start"]], 0),
                                hook.delta=signif(res_hookreg[["hook.delta"]], 0),
                                p.value=signif(res_hookreg[["p.value"]], 4),
                                CI.low=signif(res_hookreg[["CI.low"]], 2),
                                CI.up=signif(res_hookreg[["CI.up"]], 2),
                                hook.fit=res_hookreg[["hook.fit"]],
                                hook.CI=res_hookreg[["hook.CI"]],
                                hook=res_hookreg[["hook"]]
)
```

Finally a pretty printout (Table \ref{res_hookreg_table}) of the 
results from the **hookreg()** function for the `hookreg.rdml` data set with the 
following code was prepared.

```{r results='asis', echo=TRUE, message=FALSE, warning=FALSE}
# Load the xtable to create a LaTeX table from the `res_hookreg_table`.
library(xtable)
options(xtable.comment=FALSE)
print(xtable(res_hookreg_table, 
             caption = "Results from the hookreg() function for the hookreg.rdml 
             data set.", 
             label='res_hookreg_table'),
      size = "\\tiny",
      include.rownames = FALSE,
      include.colnames = TRUE,
      caption.placement = "top",
      comment=FALSE,
      table.placement = "!ht", scalebox='0.65'
)
```
The results of the **hookreg()** function are fairly comprehensive. The meaning of the columns is as followed:

\begin{itemize}
\item \emph{intercept}, is the intercept from the start of the potential hook to the end of the amplification curve.
\item \emph{slope} is the slope from the start of the potential hook to the end of the amplification curve. A negative slope is indicative for a hook effect.
\item \emph{hook.start} is the estimated starting cycle of the hook region. 
\item \emph{hook.delta} is the number of cycles from the \emph{hook.start} to the end of the amplification curve. 
\item \emph{p.value} describes the significant relationship between the variables in the linear regression model.
\item \emph{CI.low} and \emph{CI.up} is the confidence interval (low and up) for the slope parameters in the fitted linear model.
\item \emph{hook.fit} is a logical parameter indicating if the fit is significant at a default threshold of 0.005.
\item \emph{hook.CI} is a logical parameter indicating if the slope of fitted linear model is within the confidence interval (0.995).
\item \emph{ hook} is a logical parameter, which combines the significance test and confidence interval test (negative slope).
\end{itemize}


## Results for the analysis of the `hookreg.rdml` data set with **hookregNL()** \label{sectionHookregNL}

The following code was used to analyze the `hookreg.rdml` data set with **hookregNL()** function. The procedure is similar to the analysis with the **hookreg()** function.

The **hookreg()** function fits a six parameter sigmoidal model to amplification curve. The non-linear model 

$$
f(x) = c + k \cdot x + \frac{d - c}{(1 + exp(b(log(x) - log(e))))^f}
\newline
$$

is used to decide, based on the *k* parameter, if the amplification curve as a hook effect or hook effect-like curvature.

```{r, echo=TRUE, message=FALSE, warning=FALSE}
# Note that the PCRedux package needs to be loaded (see above).
# Load the qpcR package to prevent messages during the start.
suppressMessages(library(qpcR))

# `data` is a temporary data frame of the hook.rdml amplification curve data file.
# Apply the hookregNL() function over the amplification curves and arrange the 
# results in the data frame `res_hookregNL`.
# Not that `suppressMessages()` to prevent warning messages from the qpcR package.

res_hookregNL <- data.frame(sample=colnames(data)[-1],
                            t(suppressMessages(sapply(2L:ncol(data), function(i) {
                              hookregNL(x=data[, 1], y=data[, i])
                            }))))

res_hookregNL_table <- data.frame(sample=as.character(res_hookregNL[["sample"]]),
                                  slope=signif(as.numeric(res_hookregNL[["slope"]]), 1),
                                  CI.low=signif(as.numeric(res_hookregNL[["CI.low"]]), 2),
                                  CI.up=signif(as.numeric(res_hookregNL[["CI.up"]]), 2),
                                  hook.CI=unlist(res_hookregNL[["hook"]])
)
```

Finally we prepare a pretty printout (Table \ref{res_hookregNL_table}) of the 
results from the **hookregNL()** function for the `hookreg.rdml` data set with 
the following code with the code shown next.

The results of the **hookregNL()** function are less comprehensive then from the **hookreg()** function . The meaning of the columns is as followed:

\begin{itemize}
\item \emph{slope} is the slope from the start of the potential hook to the end of the amplification curve that was fitted by a six parameter model. A negative slope is indicative for a hook effect.
\item \emph{CI.low} and \emph{CI.up} is the confidence interval (low and up) for the slope parameters in the fitted linear model.
\item \emph{ hook} is a logical parameter, which combines the significance test and confidence interval test (negative slope).
\end{itemize}

```{r results='asis', echo=TRUE, message=FALSE, warning=FALSE}
library(xtable)
options(xtable.comment=FALSE)

print(xtable(res_hookregNL_table, 
             caption = "Results from the hookregNL() function for the 
             hookreg.rdml data set.", 
             label='res_hookregNL_table'),
      size = "\\tiny",
      include.rownames = FALSE,
      include.colnames = TRUE,
      caption.placement = "top",
      comment=FALSE,
      table.placement = "!ht", scalebox='0.65'
)
```


# Comparison of the **hookreg()** and **hookregNL()** methods

The decisions from the human classification (see Table \ref{Table_human_rated}) and the results from the machine decision (section \ref{sectionHookreg} and section \ref{sectionHookregNL}) were aggregated in Table \ref{method_comparison}. 

```{r, echo=FALSE, message=FALSE, warning=FALSE}
# Load the readxl package to obtain the classifications by the human.
# The classification data are stored in an EXCEL file, which is contained in the 
# PCRedux package.

library(readxl)

# Read the data as EXCEL file.
Table_human_rated <- read_xlsx(path=system.file("Table_human_rated.xlsx", 
                                                package="PCRedux"))

# Aggregate the results from the statistical analysis and the human classifications.

res <- data.frame(sample=Table_human_rated[, "Sample"],
                  hr=Table_human_rated[, "Hook effect-like\nRater 2"],
                  hookreg=res_hookreg_table[, "hook"],
                  hookregNL=data.frame(unlist(res_hookregNL[, "hook"]))
)
colnames(res) <- c("Sample",
                   "Human rater",
                   "hookreg",
                   "hookregNL"
)
```

Finally a pretty printout (Table \ref{method_comparison}) of the 
aggregated data set with the following code was prepared:


```{r, echo=TRUE, message=FALSE, warning=FALSE}
# A simple logic was applied to improve the classification result. In this case
# the assumption was, that an amplification curve has an hook effect or hook effect-like
# curvature, if either the hookreg() or hookregNL() function are positive.

meta_hookreg <- sapply(1:nrow(res), function(i){
  ifelse(res[i, "hookreg"] == 1 || res[i, "hookregNL"] == 1, 1, 0)
})

res_out <- data.frame(Sample=res[["Sample"]], res[["Human rater"]], 
                      res_hookreg[["hook"]], res_hookregNL_table[["hook.CI"]], 
                      meta_hookreg)

colnames(res_out) <- c("Sample",
                       "Human rater",
                       "hookreg",
                       "hookregNL",
                       "hookreg and hoohkreNL combined"
)
```

```{r results='asis', echo=TRUE, message=FALSE, warning=FALSE}
library(xtable)
options(xtable.comment=FALSE)

print(xtable(res_out, digits=0, 
             caption = "Aggregated decisions from the human classification and 
the results from the machine decision of the hookreg() and hookregNL() 
functions.", label='method_comparison'), 
caption.placement = "top",
scalebox='0.65')
```

The performance of the **hookreg()** and **hookregNL()** functions was analyzed 
with the **performeR()** function of the \texttt{PCRedux} package (Table 
\ref{res_performeR}). The methods were adopted from @brenner_variation_1997 and 
@kuhn_building_2008. Note that the formula for the calculations of the 
sensitivity, specificity, precision, Negative predictive value, fall-out, false 
negative rate, false discovery rate, Accuracy, F1 score, Matthews correlation 
coefficient and kappa by Cohen are described in the documentation of the 
\texttt{PCRedux} package.

```{r echo=TRUE, message=FALSE, warning=FALSE}
res_performeR <- signif(t(rbind(
  hookreg=performeR(res_out[["hookreg"]], res_out[["Human rater"]]),
  hookregNL=performeR(res_out[["hookregNL"]], res_out[["Human rater"]]),
  combined_hookreg=performeR(res_out[["hookreg and hoohkreNL combined"]], 
                             res_out[["Human rater"]])
)), 4)

colnames(res_performeR) <- c("hookreg", "hookregNL", "hookreg and hookregNL")
```

```{r results='asis', echo=TRUE, message=FALSE, warning=FALSE}
library(xtable)
options(xtable.comment=FALSE)

print(xtable(res_performeR, digits=4, 
             caption = "Analysis of the performance of both algorithms. The 
performance of the individual test and the combination of the tests is shown. 
Note that the classification improved if the hookreg() and hookregNL() function 
were combined by a logical statement. The measure were determined with the 
\\textit{performeR()} function from the \\texttt{PCRedux} package. Sensitivity, 
TPR; Specificity, SPC; Precision, PPV; Negative predictive value, NPV; Fall-out, 
FPR; False negative rate, FNR; False discovery rate, FDR; Accuracy, ACC; F1 
score, F1; Matthews correlation coefficient, MCC, Cohen's kappa (binary 
classification), $\\kappa$", label='res_performeR'),
size = "normalsize",
include.rownames = TRUE,
include.colnames = TRUE,
caption.placement = "top",
comment=FALSE,
table.placement = "!ht", scalebox='0.75'
)
```

# Funding
This work was funded by the Federal Ministry of Education and Research (BMBF) InnoProfile-Transfer-Project 03IPT611X and in part by "digilog: Digitale und analoge Begleiter für eine alternde Bevölkerung" (Gesundheitscampus Brandenburg, Brandenburg Ministry for Science, Research and Culture).
We thank Franziska Dinter (BTU) for reevaluation of the amplification curve data and Maria Tokarenko (Evrogen JSC) for wet lab experiments conduction.


# References
