\documentclass[a4paper]{article}
%\usepackage{hyperref}
\usepackage[utf8]{inputenc} %unicode support
\usepackage[margin=0.75in]{geometry}
\usepackage{titlesec}
\usepackage{natbib}

%\usepackage{xr}
%\externaldocument{Roediger_OxBio}
\bibliographystyle{natbib}
\newcommand\sectionbreak{\ifnum\value{section}>1\clearpage\fi}
\newcommand\subsectionbreak{\ifnum\value{subsection}>1\clearpage\fi}
\newcommand\subsubsectionbreak{\ifnum\value{subsubsection}>1\clearpage\fi}

\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\thetable}{S\arabic{table}}
\renewcommand{\theequation}{S\arabic{equation}}
\DeclareUnicodeCharacter{00B1}{ }

\makeatletter
\renewcommand{\@maketitle}{
\newpage
\null
\vskip 2em%
\begin{center}%
{\large \@author \par}%
{\LARGE \@title \par}%
\end{center}%
\par} \makeatother

\title{Supplement to: ``chipPCR: an R Package to Pre-Process Raw Data of Amplification Curves''}
\author{Stefan R\"{o}diger*, Micha\l{} Burdukiewicz and Peter Schierack}
\date{}

\begin{document}

\maketitle




%\usepackage{Sweave}
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Supplement to: "chipPCR: an R Package to Pre-Process Raw Data of Amplification Curves"}
%\VignetteDepends{chipPCR}
%\VignetteKeywords{amplification curve pro-processing chipPCR R}
%\VignettePackage{chipPCR}

<<captures,echo=FALSE>>=
library(knitr)
opts_chunk$set(fig.lp="figure:",tidy=TRUE, tidy.opts=list(width.cutoff=60))

fig1_cap <- "The amplification curves were generated with the 
\\textsl{AmpSim} function. All Cqs are unique since random 
values were added to the starting Cq of 25. The parameter $noise = 
0.03$ adds some scatter to the data used for plotting the amplification curves."

fig1_scap <- "Simulation of a qPCR experiment using \\textsl{AmpSim} 
function."


fig2_cap <- "Working principle of the \\textsl{th}.\\textsl{cyc} function.
The function provides two modes (\\textbf{A)} is the linear 
regression. \\textbf{B)} Quadratic regression for the calculation of the 
Cq. In both cases the highest R squared value determines how
many left and right neighbors are used above and the below the defined threshold level ."

fig2_scap <- "Working principle of the \\textsl{th}.\\textsl{cyc} function."

fig3_cap <- "Application of the \\textsl{th}.\\textsl{cyc} function for the analysis of 
ccPCR data. The \\textsl{CPP} function was used to pre-process the data. 
Subsequently, the data were analyzed using the \\textsl{th}.\\textsl{cyc} function 
in the linear regression mode. The threshold level ($r = 50$) was 
identical for all data. The Cq is shown in minutes. The range used 
for the calculation of the Cq is indicated in red. Negative curves are 
automatically excluded from the analysis if the 90\\% percentile is lower 
or equal to the threshold level ($r$)."
fig3_scap <- "Application of the \\textsl{th}.\\textsl{cyc} function for the analysis 
of ccPCR data."

fig4_cap <- "Application of the \\textsl{CPP} and \\textsl{th}.\\textsl{cyc} 
functions. \\textbf{A)} The raw data of the VIMCFX96\\_60 dataset were 
\\textbf{B)} pre-processed with the \\textsl{CPP} function and finally plotted. 
The parameter $trans$ was set to $TRUE$, which leads to a linear trend correction 
and base-lining. By default a Savitzky-Golay filter was used to smooth the data. 
The data were normalized between 0 and 1 ($method.norm = 'minm'$). \\textbf{C)} 
All Cqs were calculated with \\textsl{th}.\\textsl{cyc} function. The Cq for the 
raw data was $17.25 \\pm 0.5$ (at $r = 2575$) and $17.1 \\pm 0.1$ (at $r = 0.1$) 
for the pre-processed data. Our results indicate that the dispersion of the Cq 
values was slightly lower for the pre-processed data."
fig4_scap <- "Application of the \\textsl{CPP} and 
\\textsl{th}.\\textsl{cyc} functions."


fig5_cap <- "Comparison of the normalization methods with the \\textsl{CPP} 
function. The VIMCFX96\\_60 dataset (96-well plate cycler, Bio-Rad CFX96, 
EvaGreen detection) was used. \\emph{(A)} Plot of raw data for all amplification 
curves. The signals are superimposed to circa 2200 RFU and the inter-sample 
baseline and plateau shift is high. Note the positive trend 
(\\textcolor{red}{\\textendash}, fitted with an ordinary least squares method) 
in the background range of cycles 1 to 15. All subsequent plots were processed 
with the \\textsl{CPP} function. By default, the curves are base-lined, smoothed 
(Savitzky-Golay smoother) and the slope corrected by a linear regression ($trans 
= TRUE$). \\emph{(B)} base-lined raw data, \\emph{(C)} \\emph{Min-Max 
normalization}, \\emph{(D)} \\emph{Max normalization}, \\emph{(E)} 
\\emph{lugn-normalization} with a cut off 3\\% 
and \\emph{(F)} \\emph{zscore-normalization}."

fig5_scap <- "Comparison of the normalization methods with the \\textsl{CPP} function."

fig6_cap <- "Amplification standard curve simulation and regression analysis. 
\\emph{(A)} \\textsl{AmpSim} was used to synthesize a qPCR experiment of six 
dilutions (three replicates per dilution) standard samples. The Cqs were 
determined by the $SDM$ method (solid black vertical lines). \\emph{(B)} \\textsl{effcalc} 
was used to automatically perform a linear regression (decadic logarithm 
of input concentration versus the Cq). The 95\\% confidence interval is shown be 
the light-blue solid lines."
fig6_scap <- "Amplification standard curve simulation and regression 
analysis."

fig7_cap <- "Calculation of the amplification efficiency.
Data of a VideoScan HCU dilution experiment (C54 dataset) 
were analyzed. \\emph{(A)} Visualization of the raw data. One of the three 
dilutions contains a missing value due to a sensor error. \\emph{(B, top 
panel)} The \\textsl{CPP} function was used to baseline, remove the missing value 
(\\textcolor{red}{\\textendash}) and smooth 
(\\textcolor{black}{\\textendash}, \\textcolor{red}{\\textendash}, 
\\textcolor{green}{\\textendash}) the raw data. \\emph{(B, bottom panel)}. 
The Cqs ($SDM$ by the \\textsl{diffQ} function) were then determined. An amplification 
efficiency of 87.3~\\% was calculated with the \\textsl{effcalc} function."

fig7_scap <- "Calculation of the amplification efficiency."

fig8_cap <- "Imputation of missing values in the data for generating amplification curves. 
\\emph{(A)} Raw data were generated using the \\textsl{AmpSim} 
simulation function. \\emph{(B)} A missing value was introduced in the 
transition phase. The missing value was imputed either by \\emph{(C)} 
linear approximation or \\emph{(D)} a cubic spline approximation. The 
spline approximation nearly reconstituted the original curve."

fig8_scap <- "Imputation of missing values in the data for generating amplification curves."

fig9_cap <- "Quantification cycle (Cq) by the second derivative maximum 
method. Raw data (\\textbullet) were generated by the \\textsl{AmpSim} 
function. The inflection point is the point where 
the slope is maximum and the curvature is zero. The first derivative of the 
amplification curve has a first derivative maximum ($FDM$) at the 
inflection point. The second derivative maximum method ($SDM$) needs to 
differentiate a curve to the second order prior to quantification. The second 
derivative exhibits a zero-crossing at the $FDM$. The function $y = f(x)$ is 
numerically derived by the five-point stencil. This method is
assumption free regarding the function $f$. \\textsl{inder} calculates the 
approximate $SDM$. The $SDM$ is calculated from a derived cubic spline. 
Similarly, the first approximate derivative maximum ($FDM$), second derivative 
minimum ($SDm$), and approximate second derivative center ($SDC$, geometric mean 
of $SDM$ and $SDm$) are available. $FDM$, $SDm$ and $SDC$ values can be used to 
further characterize the amplification process."

fig9_scap <- "Quantification cycle (Cq) by the second derivative maximum 
method."

fig10_cap <- "Plot of all data from C127EGHP and calculate the $SDM$ (Second 
Derivative Maximum) values with the \\textsl{diffQ2} function.
\\emph{(A)} Plot the samples detected with EvaGreen and 
\\emph{(B)} shows the same samples detected with the Hydrolysis probe for 
\\textsl{MLC-2v}. \\emph{(C)} Stripchart of the Cq values (\\textbullet) with the 
median (\\textcolor{red}{\\textendash}) and the median absolute deviation 
(\\textcolor{blue}{\\textendash~\\textendash}). 
This result indicates, that the variance derived from detection with hydrolysis
probes is higher than for the samples detected with EvaGreen. Note: the $inder$ parameter is set as TRUE."
  
fig10_scap <- "Plot of all data from C127EGHP and calculate the $SDM$ (Second 
Derivative Maximum) values with the \\textsl{diffQ2} function."

fig11_cap <- "Amplification curve profiles from the Bio-Rad iQ5 thermo 
cycler 
for the human gene \\textit{HPRT1}.
\\emph{(A)} The \\textsl{C126EG595} dataset was used with 
96 replicates of equal starting numbers of template molecules. Vertical 
lines represent the Cq ($SDM$ method) determined by the $inder$ method on 
amplification curves fitted with a 5-parameter curve function. Curves with 
Cqs less than 14.5 are indicated in red (\\textcolor{red}{\\textendash}). 
\\emph{(B)} Second derivatives of the amplification curves. Note that 
after differentiation all inter sample baseline and plateau shifts are 
similar. \\emph{(C)} Histogram (class width = 0.05 Cq) of the Cq values 
($SDM$). Cqs were mainly at circa 15.7 (N = 80) while some amplification 
curves had a Cq less than 15.5 (N = 16)."
fig11_scap <- "Amplification curve profiles from the Bio-Rad iQ5 thermo 
cycler 
for the human gene \\textit{HPRT1}."

fig12_cap <- "Signal analysis using the VIMCFX96\\_60 dataset (96-well 
plate cycler (Bio-Rad CFX96)). All cycles (ROI: 1 -- 40) were analyzed by 
the \\textsl{MFIaggr} function. The density plot (right upper panel) and 
quantile-quantile analysis (right lower panel) show no normal 
distribution. Owing to the sigmoidal structure of the curve, the density function can be considered bimodal."

fig12_scap <- "Signal analysis using the VIMCFX96\\_60 dataset (96-well 
plate cycler (Bio-Rad CFX96))."

fig13_cap <- "Helicase-dependent amplification (HDA) of \\textsl{vimentin} (\textit{vim}) in the 
VideoScan Platform. The HDA was performed at 65 degrees Celsius. Three 
concentrations of input DNA (D1, D2, D3) were used. The amplification curves 
were smoothed by a moving average (windowsize 3) and base-lined by a robust 
linear regression by computing MM-type regression estimator. The 
\\textsl{th}.\\textsl{cyc} function was used to determine the time required to 
reach the threshold level of 0.05 (--) arbitrary fluorescence units." 

fig13_scap <- "Helicase-dependent amplification (HDA) of \\textsl{vimentin} (\textit{vim})."

fig14_cap <- "The \\textsl{plotCurves} function. Notice: The function plots many 
curves on one plot in separate cells allowing for quick assessment. Missing 
values were artificially introduced at random positions to selected curves of the 
VIMCFX96\\_60 data set (solid black line). A colored box (topleft of each plot) 
indicates the sample name and if the data contain missing values. The red rug 
indicates the position of the missing values. The red lined shows the 
amplification curve after unsupervised pre-processing (using an instance of 
\\textsl{CPP})."
fig14_scap <- "The \\textsl{plotCurves} function."

fig15_cap <- "Use of \\textsl{MFIaggr} to test for heteroskedacity. The data were 
aggregated with the \\textsl{MFIaggr} function and assigned to the object $res$. 
The standard deviation was transformed to the variance. The plot shows the cycle 
dependent variance measured at 60 degrees Celsius (annealing phase; A, B) and 69 
degrees Celsius (elongation phase, C, D) of 96 qPCR replicate amplification 
curves. The first cycles 1 to 10 and next the cycles 1 to 40 were analyzed. The 
Breusch-Pagan test of the \\textsl{MFIaggr} confirmed the heteroskedasticity in the 
amplification curve data. The VIMCFX96\\_60 and VIMCFX96\\_69 datasets were 
used."

fig15_scap <- "Use of \\textsl{MFIaggr} to test for heteroskedasticity."

fig16_cap <- "The \\textsl{inder} function calculates numeric derivatives on smoothed 
data, which results in data points not observable in reality. The rounder 
function averages such result to the real values of cycle number. An 
amplification curve was simulated with the \\textsl{AmpSim} function."

fig16_scap <- "Application of rounder to average numeric derivatives to the real 
values of cycle number."

fig17_cap <- "\\textsl{lm.coefs} a function to compute coefficients for linear 
models. The function is a wrapper for functions to perform normal (least 
squares) and robust linear regression. If the computation of the robust linear regression fails, 
then \\textsl{lm.coefs} performs a linear regression using the least squares 
method. lmrob, MM-type estimators for linear regression; rq, quantile regression 
fit; least, least squares linear regression; rfit, Rank-based estimates of 
regression coefficients. m, slope; n, asymmetry."

fig17_scap <- "\\textsl{lm.coefs} a function to compute coefficients for linear 
models."

fig18_cap <- "\\textsl{bg}.\\textsl{max} tries to estimate the range between the 
background and the plateau phase of an amplification reaction. \\emph{(A)} in 
absence and \\emph{(B)} presence of noise. The data were simulated with the 
\\textsl{AmpSim} function."

fig18_scap <- "\\textsl{bg}.\\textsl{max} tries to estimate the range between the 
background and the plateau phase of an amplification reaction"

fig19_cap <- "Application of the \\textsl{bg}.\\textsl{max} function. 
Amplification curve data from a capillary convective PCR were used \\emph{(A)} 
as raw data and \\emph{(B)} pre-processed (smoothed (moving average, window size 
3), base-lined and trend corrected (robust MM-estimator)) with the \\textsl{CPP} 
function. The output of the \\textsl{CPP} function was used by 
\\textsl{bg}.\\textsl{max} to detected the start and the end of the 
amplification reaction. The start and end were reliably estimated (range between 
``bg.stop'' and ``amp.stop''). There was no significant difference between raw 
and pre-processed data."

fig19_scap <- "Application of the \\textsl{bg}.\\textsl{max} function to detect 
the start and end of an amplification reaction in a capillary convective PCR."

fig20_cap <- "Inspection of the reps384 dataset. The reps384 dataset was used 
for the analysis of the impact of imputed missing values. Three areas of the 
curve data were defined as ``Linear phase'' (red, cycle 1 -- 10), ``Exponential 
phase'' (blue, cycle 11 -- 33), ``Plateau phase'' (green, cycle 34 -- 40)."

fig20_scap <- "Inspection of the reps384 dataset."

fig21_cap <- "Analysis and interpretation of real-time amplification curves. 
\\emph{(A)} Before procession: The fluorescence values are plotted against the cycle, which results 
in sigmoidal shaped amplification curves 
(\\textcolor{black}{\\textendash},~\\textcolor{red}{\\textendash}). Measurements may 
occasionally contain missing values (``NA'', \\textcolor{red}{\\textendash}) and 
outliers (orange circle, \\textcolor{black}{\\textendash}), due to noise 
introduced by the detection system or sensor errors. Outliers are often 
present in the first cycle due to sensor adjustments. The signal difference 
between the background phase (first cycles) and the plateau phase (last cycles) 
is quantifiable as signal-to-noise ratio (SNR). The SNR between different 
samples (e.g., \\textcolor{black}{\\textendash} and 
\\textcolor{black}{\\textendash}) can vary. For interpretation it is recommended 
to normalize the data. Negative samples (\\textcolor{blue}{\\textendash}) need 
to be identified automatically. \\emph{(B)} Pre-processed raw data, where NAs 
were imputed and the noise slightly removed. The curves were adjusted to have 
the same baseline and plateau level. The quantification cycles (Cq) of the 
positive reactions are determined in the exponential phase (``threshold method'' 
is used in this example). Negative samples are automatically set to zero." 

fig21_scap <- "Analysis and interpretation of real-time amplification curves."

fig22_cap <- "Smoother and filter methods of the \\emph{chipPCR} package. 
\\emph{(A)} Raw data were generated using the \\textsl{AmpSim} simulation 
function. \\emph{(B)} The difference of the raw data to the smoothed data was 
plotted. ``savgol'' (Savitzky-Golay smoothing), ``lowess'' (locally-weighted 
polynomial regression), ``mova3'' (moving average with window size of 3), 
``smooth'' (cubic smoothing spline), ``spline'' (Interpolating cubic spline), 
``supsmu'' (Friedman's SuperSmoother), ``whit1'' (weighted Whittaker smoothing 
with a finite difference penalty of order 1), ``whit2'' (weighted Whittaker 
smoothing with a finite difference penalty of order 2). The ``savgol'', 
``smooth'', ``spline'' ``whit1'' , and ``whit2'' nearly preserved the original 
curve. The other functions resulted in alterations in the transition phases of 
the amplification curve. Optimized time series smoother, like the Kalman filter 
\\citep{Tusell_2010}, are not yet integrated in this package."

fig22_scap <- "Smoother and filter methods of the \\emph{chipPCR} package."

fig23_cap <- "Sample for analysis using the \\textsl{MFIaggr} function. 
The VIMCFX96\\_60 dataset (96-well plate cycler (Bio--Rad CFX96)) was used. 
Either all or a subset of the cycles (ROI: 1 -- 10) or all cycles (ROI: 1 -- 40) 
(Figure~\\ref{figure:MFIaggr_all}) were analyzed. The density plot (right upper 
panel) and quantile-quantile analysis (right lower panel) show no normal 
distribution. Due to the sigmoidal structure of the curve, the density function 
can be considered bimodal."

fig23_scap <- "Signal analysis by the \\textsl{MFIaggr} function."

fig24_cap <- "Analysis of two PCR conditions with the \\textsl{MFIaggr} function. 
The VIMCFX96\\_60 and VIMCFX96\\_69 datasets (96-well plate cycler (Bio--Rad 
CFX96)) were used. A subset of cycles (ROI: 2 -- 10) was analyzed for a qPCR 
measured during the annealing phase and elongation phase. The density 
plot (right upper panel) and quantile-quantile analysis (right lower panel) show 
no normal distribution. The measurements during the annealing phase and 
elongation phase show a similar distribution but differ in their mean fluorescence intensity (MFI) levels."

fig24_scap <- "Analysis of two PCR condition with the \\textsl{MFIaggr} function."

fig25_cap <- "Calculation of the amplification efficiency (AE) from a qPCR 
experiment. DNA of \\textsl{vimentin} (\\emph{A}) and \\textsl{MLC-2v} (\\emph{B}) was diluted and 
amplified in a Roche Light Cycler 1.5 (C60.amp dataset). Cqs ($SDM$) were 
calculated by the \\textsl{inder} function and analyzed with \\textsl{effcalc}. 
The AE was (\\emph{C}) 95.1~\\% for \\textsl{vimentin} and  
(\\emph{D}) 94.1~\\%  for \\textsl{MLC-2v}."

fig25_scap <- "Calculation of the amplification efficiency for the C60.amp dataset."

@

\begin{figure}[ht]
\centering
\scalebox{0.6}{
\includegraphics[clip=true,trim=7cm 7cm 7cm 7cm]{logo.pdf}
}
\end{figure}

%create TOC
\tableofcontents

\pagebreak

\begin{abstract} 
\textbf{Background: } 
Both the quantitative real-time polymerase chain reaction (qPCR) and isothermal 
amplification are standard methods used for the quantification of nucleic acids (DNA, 
RNA). Numerous real-time read-out technologies with different technical 
platforms have been developed so far. However, analysis of amplification curves consists of cascaded steps implemented in a similar manner across all existing technological platforms. Despite the growing interest in amplification-based techniques, there are 
only few open source tools for pre-processing real-time amplification data. The availability of a software for pre-processing raw amplification data is mandatory in different scenarios, for example during the development, optimization and improvement of the functionality of instruments.

\textbf{Results and Conclusion: } %if any 
$\emph{chipPCR}$ is a versatile \textbf{R} package for pre-processing 
(e.g., imputation of missing values, smoothing) and quality analysis of raw data for amplification curves coming from conventional quantitative polymerase chain reactions (qPCR), and quantitative isothermal amplification (qIA) reactions. The package 
contains datasets, which were generated by helicase-dependent amplification 
(HDA) or polymerase chain reaction (PCR) under various temperature conditions 
and detection systems, such as hydrolysis probes and intercalating dyes. The 
structure of the packages is amenable for integration to Web-based and standalone 
\emph{shiny} applications.

\end{abstract}

\section{Availability, requirements and setting up a working environment}
\label{sec:work_environment}

\emph{Convention: According to the MIQE guidelines (``Minimum Information for 
publication of Quantitative real-time PCR Experiments'' 
\citep{bustin_miqe_2009}) is the threshold cycle (Ct) referred to as 
quantification cycle (Cq). We use the expression Cq exclusively regardless of 
the amplification method and mathematical principle. The abbreviations MFI, RFU 
and refMFI refer to arbitrary units of mean fluorescence intensities.}

The software is available in an \textbf{R} environment or through web browser applications:
\begin{itemize}
\item Project name: chipPCR, 
\item Project homepage (development): https://github.com/michbur/chipPCR, 
\item Project homepage at CRAN: http://cran.r-project.org/web/packages/chipPCR/index.html, 
\item Operating System: Platform independent, 
\item Other requirement: R 3.1.0 or higher, 
\item License: GPL-3
\end{itemize}

We use \textbf{R}'s \emph{S4} class system (see \emph{methods} 
package) to separate the interface and the implementation because, unlike the 
\textbf{R}'s \emph{S3} class system, it requires the explicit declaration of classes and the inheritance and relationship for each class or method. Therefore, the number and types of objects in slots in an instance of a class 
have to be established at the time of the class definition. Objects belonging to the 
class are validated against this definition and have to fulfill it at any time. 
\textsl{setGeneric} declares generic functions. \emph{S4} methods are declared by calls 
to \textsl{setMethod} together with the name of generics and signatures of the arguments.
Signatures are used for identification of classes of one or more arguments of the methods. 

\emph{S4} classes therefore require a higher development effort than \emph{S3} 
classes, but offer more stringent control over the contents of created objects. Additional 
information (e.g., results, parameters) can also be included. \emph{S4} assures better control 
on the object structure and chosen method dispatch \citep{Karatzoglou_2004}. 

For high-throughput capability, we avoided loops in the core structures of the 
$\emph{chipPCR}$ package and partially used parallel computing (\textsl{smoother} 
function) to keep the code fast. This package supports the use of most popular R packages, thus providing a
communication layer required for parallel computing: \emph{parallel} and \emph{snow}.

This vignette can be viewed from \textbf{R} using command: vignette("chipPCR"). 
All experimental details for the datasets are described in the $\emph{chipPCR}$ 
package manual and in the citations. To start an analysis it is required to 
choose a dataset (as shown below).

<<load_data,message=FALSE,results='asis'>>=
# Load chipPCR 
require(chipPCR)
# Load package for table formatting
require(xtable)
# Print table
print(xtable(head(C60.amp[, 1L:5]), caption = "First five cycles of imported data."))
@

All datasets used in the following examples can be loaded similarly. $\emph{chipPCR}$ 
relies on the \textbf{R} environment, and dedicated \textbf{R} packages 
(e.g., $\emph{RDML}$) as default data format and standard import and export formats \citep{perkins_2012, blagodatskikh_2014, RDCT2014c}.

Graphical user interfaces (GUIs) are important to make software usable for researchers not fluent in \textbf{R}. Several 
\textbf{R} GUI projects have been proposed \citep{rodiger_rkward_2012}. Selected functionality of 
$\emph{chipPCR}$ originates from \textbf{RKWard} GUI plugins 
\citep{pabinger_2014}. The \emph{shiny} framework \citep{shiny_2014} can be used 
to build and deploy GUIs for the desktop or as services for 
interactive web applications on servers. Prerequisites are an installation of 
\textbf{R}, with installed $\emph{shiny}$ and $\emph{chipPCR}$ packages and a 
modern web browser. It is possible to run a GUI in a web browser on a local 
machine without an Internet connection. Ad-blocking software may cause malfunctions 
and should be turned off. \emph{shiny} enables the building of plugin-like applications with highly customizable widgets (e.g., sliders) for an 
efficient extension. \emph{shiny} applications can be updated live and in an interactive manner. 
The user interfaces can be built entirely using \textbf{R} and operates in any 
\textbf{R} environment (cross-platform). The functions \textsl{AmpSim}, 
\textsl{th.cyc}, \textsl{bg.max} and \textsl{amptester} are a part of 
\emph{shiny} GUIs. Examples and case studies for the mentioned functions are presented in 
the following sections.

\section{Pre-processing raw data for DNA amplification plots}
\label{sec:problems}

Quantitative polymerase chain reaction (qPCR) and quantitative isothermal 
amplification (qIA) are standard methods used for amplifying nucleic acids (e.g., genomic 
DNA, copy DNA) \citep{pabinger_2014}. The Taguchi methods provide general optimization frameworks used in engineering 
optimization processes and other related disciplines. However, applications include 
PCR application too \citep{cobb_1994, thanakiatkrai_2012}. Basically, it is possible 
to determine the optimal conditions for PCR-based assays. A software implementation for 
\textbf{R} is available in the $\emph{qualityTools}$ package \citep{qualityTools}.
Recently amplification 
methods with continuous temperature gradients (e.g., microfluidics, capillary 
convective PCR (ccPCR)) emerged 
\citep{chou_rapid_2011,rodiger_nucleic_2014,spiess_impact_2014}. Isothermal 
amplification is a monocyclic reaction at a constant temperature. Conversely, 
PCR is a polycyclic reaction with repeated thermal cycling condition steps 
(denaturation, annealing, elongation) and measurements at discrete cycle steps. 
The curve shape of isothermal amplification reactions do not necessarily follow 
an S-shaped structure and the measurement is time-based (continuous, not 
mandatory equidistant) in contrast to a cycle-based (discontinuous) measurement 
of qPCRs. The number of measure points is usually higher than in qPCR experiments. All these 
amplification methods are used in different real-time monitoring technologies, 
such as our previously reported VideoScan technology 
\citep{roediger_highly_2013}, microfluidic systems, point-of-care devices, 
microbead-chip technologies and commercial real-time thermo cyclers 
\citep{chang_2012, roediger_bead_qPCR_2013, rodiger_nucleic_2014}. Real-time 
technologies enable the quantification of nucleic acids by calculation of 
specific curve parameters like the quantification cycle (Cq) and the 
amplification efficiency (AE) \citep{ruijter_2013, tellinghuisen_2014}.

The data quality of experimental instruments is often not suitable for end-user 
analysis and presentation. Moreover, analysis of raw data may lead to 
misinterpretations and false performance estimates under certain conditions. 
Therefore, novel technologies usually depend on software for pre-procession of 
raw data. Pre-processing specifically addresses raw data inspection, steps to 
transform raw data in a compatible format for successive analysis steps (e.g., 
smoothing, imputation of missing values), data reduction (e.g., removal of 
invalid sets), noise reduction and data quality management. Noise is challenging 
because derivative processes as used for ``cycle of quantification'' methods 
(e.g., Second Derivative Maximum method) leads to an amplification of noise 
\citep{larionov_2005, tuomi_2010, roediger_RJ_2013, ruijter_2013, 
tellinghuisen_2014}. Pre-processing algorithms remove stochastic errors and 
artefacts (e.g., noise, photo-bleaching effects, degassing effects, different 
signal levels) as illustrated in Figure~\ref{figure:problems}. However, 
misinterpretations are more likely if arbitrary manual corrections are not 
performed. A manual alteration is in contradiction to reproducible research. In 
particular, open source scientific software and the associated input and output 
data are central structural elements to enable recomputability and 
reproducibility of results \citep{Blanton_2014, Jacobs_2014, Stodden_2014}. 

\textbf{R} is widely used for the analysis of qPCR data. Most \textbf{R} 
packages focus on the read-in, processing and post-processing of datasets 
originating from commercial qPCR systems. The fundamental steps of amplification curve 
analysis are: (1) raw data read-in, (2) amplification curve pre-processing 
(e.g., noise reduction, outlier removal), (3) amplification curve processing 
(e.g., Cq and AE calculation), (4) post-processing and quantification of 
secondary parameters (e.g., Delta-Delta-Ct for gene expression analysis), (5) 
data export, (6) visualization and (7) report generation. Sophisticated 
\textbf{R} packages for the steps 1 and 3--7 are available from Bioconductor and 
CRAN \citep{Dvinge_2009, Zhang_2010, heckmann_2011, perkins_2012, 
gehlenborg_2013, huntley_2013, zhang_2013, mccall_2014, pabinger_2014}. However, 
there is no \textbf{R} package for pre-processing and quality analysis of 
raw data for amplification curves. This need also applies to other existing software solutions (compare 
\citep{pabinger_2014}). Pre-processing in most commercial cyclers is a ``black 
box'' model, where inner subroutines are not available for inspection. This approach limits
reproduction of analysis on other platforms and introduces difficulties for transfer of experimental design and setup or for the adequate use of statistical tools. Moreover, it is desirable to set up 
work-flows in an open environment, which enables downstream analyses and offers 
powerful tools for data visualizations and automatic report generation.

The $\emph{chipPCR}$ package was developed to automatize pre-processing, ease 
data analysis/visualization and offer a quality control for the statistical data 
analysis of qPCR or qIA experiments (see Section~\ref{sec:functions}). 
$\emph{chipPCR}$ is primarily targeted at developers of novel systems. 
Nonetheless, users who process raw data of commercial systems can also utilize 
its functionalities.

<<problems,warning=FALSE,message=FALSE,fig.show='hold',fig.cap=fig21_cap,fig.scap=fig21_scap,warning=FALSE,fig.width = 11, fig.height = 8>>=
# Use AmpSim to generate amplification curves with 40 cycles
# and different Cq's.
res.pos <- AmpSim(cyc = 1:40, noise = TRUE, b.eff = -12, nnl = 0.02)
res.pos[5, 2] <- res.pos[5, 2] * 6

res.low <- AmpSim(cyc = 1:40, noise = TRUE, b.eff = -20, bl = 0.5, 
		  ampl = 0.58, Cq = 33)
# Add missing value to res.low at cycle 31
res.low[31, 2] <- NA

res.neg <- AmpSim(cyc = 1:40, b.eff = -0.1, bl = 0.05, ampl = 0.4, Cq = 1, 
		  noise = FALSE, nnl = 0.5)
		      
res.pos.CPP <- cbind(1:40, CPP(res.pos[, 1], res.pos[, 2], 
		     bg.outliers = TRUE, smoother = TRUE, method = "smooth", 
		      method.norm = "minm", method.reg = "lmrob")$y)
		      
res.low.NA <- cbind(1:40, CPP(res.low[, 1], res.low[, 2], smoother = TRUE, 
		    method = "smooth", bg.outliers = TRUE, method.norm = "minm", 
		    method.reg = "lmrob")$y)
		      
res.neg.exc <- cbind(1:40, amptester(res.neg[, 2]))

par(mfrow = c(1,2), las = 0, bty = "n", cex.axis = 1.5, cex.lab = 1.5, 
    font = 2, cex.main = 1.8, oma = c(1,1,1,1))
plot(NA, NA, xlim = c(1,40), ylim = c(0, max(res.pos[, 2])), xlab = "Cycle", 
     ylab = "Raw fluorescence")
mtext("A", cex = 2, side = 3, adj = 0, font = 2)

lines(res.pos, lwd = 2)
lines(res.low, col = 2, lwd = 2)
arrows(38, min(res.low[, 2], na.rm = TRUE), 38, max(res.low[, 2], 
       na.rm = TRUE), code=3, lwd=3, angle=90, col="grey")
text(38, max(res.low[, 2], na.rm = TRUE) * 0.7,"SNR", cex=1.2)

arrows(29,0.42,31,0.51, lwd=2)
text(29,0.38, "NA", cex=1.2)

points(res.pos[5, 1], res.pos[5, 2], pch=21, cex=4, lwd=5, col="orange")
text(res.pos[5, 1], res.pos[5, 2] * 1.2, "Outlier", cex=1.2)

lines(res.neg, col = 4, lwd = 2)
text(20, mean(res.neg[, 2]) *0.9, "No amplification", cex=1.2, col = 
"blue")


plot(NA, NA, xlim = c(1,40), ylim = c(0, max(res.pos[, 2])), xlab = "Cycle", 
     ylab = "Pre-processed fluorescence")
abline(h = 0.03, lty = 2, lwd = 2)
mtext("B", cex = 2, side = 3, adj = 0, font = 2)

lines(res.pos.CPP, lwd = 2)
lines(res.low.NA, col = 2, lwd = 2)
lines(res.neg.exc, col = 4, lwd = 2)

legend(1, 1, c("Positive (outlier removed)", "Positive (scaled)", 
       "Negative", "Threshold line nof Cq"), 
       col = c("black", "red", "blue", "black"), lty = c(1,1,1,2), 
       lwd = 2, bty = "n")
lines(c(15.1,15.1), c(-1,0.03), lwd = 2, col = "black")
text(14, 0.06, "Cq")
lines(c(28.5,28.5), c(-1,0.03), lwd = 2, col = "red")
text(27, 0.06, "Cq", col = "red")
@

\section{Functions of the chipPCR package}
\label{sec:functions}

The main functions of the $\emph{chipPCR}$ package are:

\begin{itemize}
\item \textsl{AmpSim}: simulatea S-shaped amplification plots based on a 5-parameter model accompanied by \textsl{AmpSim.gui} a \emph{shiny} GUI, for \textsl{AmpSim},
\item \textsl{bg.max}: detects the start and end of an 
amplification reaction,
\item \textsl{CPP}: easily accesses several pre-processing functions,
\item \textsl{fixNA}: imputes missing values in a data column,
\item \textsl{inder}: interpolates first and second derivatives 
interpolation using the five-point stencil (accompanied by 
\textsl{rounder} function),
\item \textsl{MFIaggr}: analyzes a bulk of replicates of an amplification 
reaction, 
\item \textsl{smoother}: smoothens the data for curve plotting by different methods 
(e.g., moving average, Savitzky-Golay smoothing).
\end{itemize}

Additionally, further auxiliary functions are:

\begin{itemize}
\item \textsl{amptester}: detects the start and end of an amplification curve, 
\item \textsl{effcalc}: calculates the amplification efficiency, 
\item \textsl{humanrater}: rates curves using a graphical interface, 
\item \textsl{lm.coefs}: computes linear model coefficients, 
\item \textsl{normalizer}: normalizes data between a user defined range, 
\item \textsl{plotCurves}: plots many curves on one plot in separate cells allowing for quick assessment, 
\item \textsl{th.cyc}: calculates the number of cycles for which the fluorescence reporter signal exceeds a defined threshold, called the threshold cycle (Cq),
\end{itemize}

These auxillary functions are used for post-processing (e.g., Cq calculation, 
AE calculation) and quality analysis. Here we provide more information on the functionality of the chipPCR package. Selected functionality is used in the 
\emph{RDML} \citep{blagodatskikh_2014} package. Most of the functions are 
central elements of other $\emph{chipPCR}$ functions. For example, 
\textsl{fixNA} is embedded in most functions to prevent errors due to missing 
values.

\section{\textsl{AmpSim} - a function for simulating amplification curves}
\label{sec:AmpSim}

The \textsl{AmpSim} function simulates amplification reactions. Use 
cases include teaching, algorithm testing or the comparison of an experimental 
system to the predicted (``optimal'') model. \textsl{AmpSim} uses a 5-parameter 
model (Equation~\ref{eq:richards}), which is commonly used for the 
simulation of amplification curves \citep{ritz_2008, spiess_2008, gerhard_2014}.

\begin{equation} \label{eq:richards}
fluo = bl + \frac{ampl - bl}{1 + \exp{(b.eff * (\log{cyc} - \log{Cq}))}}
\end{equation}

\textsl{AmpSim} has the intrinsic property to generate unique results if the 
$noise$ parameter is set to $TRUE$. This is due to the addition of normally distributed noise (as
per \textsl{rnorm} function from \emph{stats} package), for identically replicated
noisy dataset random seed (for example \textsl{set.seed(123)}). The 
amplification curves of Figure~\ref{figure:AmpSim_effcalc}~\emph{A} were 
generated with the same starting parameter of \textsl{AmpSim} with some noise 
added. \textsl{AmpSim}.\textsl{gui} is a \emph{shiny} GUI \citep{shiny_2014} 
implementation for \textsl{AmpSim}. A GUI for the simulation and analysis of 
amplification reactions can be invoked by pasting the following code snippet 
in an \textbf{R} console.

<<AmpSim_GUI,eval = FALSE>>=
# Load the shiny package (chipPCR should already be loaded).
# Run from a R console following commands.
require(shiny)

# Invoke the shiny AmpSim app in the default web browser.
runApp(paste(find.package("chipPCR")[1],"/AmpSim.gui", sep = ""))

# Alternatively call shiny app AmpSim from gist
runGist('https://gist.github.com/michbur/e1def41598f1d0c1e2e6')
@

The function opens a $\emph{chipPCR}$ webpage in a default web browser 
(Figure~\ref{figure:browser}). All parameters of the \textsl{AmpSim} 
function may be set in the left part of the interface. 
In addition, the GUI shows some information calculated by the 
\textsl{bg}.\textsl{max} function in a summary field and a plot below the simulated 
amplification curve.

\begin{figure}[ht]
\centering
\scalebox{0.43}{
\includegraphics{browser.png}
}
\caption{Locally running \emph{shiny} \textsl{AmpSim}.\textsl{gui} app. 
\emph{(Top)} The plot of the \textsl{AmpSim}.\textsl{gui} is shown in a web 
browser (\textbf{Iceweasel}, v. 29.0.1) along with the parameters (left panel) 
and the estimation by the \textsl{th}.\textsl{cyc} function. The code 
(``server.R'', ``ui.R'') of the \emph{shiny} app is shown in the right panel. All parameters 
(e.g., Cq value, baseline) of the \textsl{AmpSim} function are accessible. 
\emph{(Bottom)} Additionally, \textsl{AmpSim}.\textsl{gui} shows the plot output 
and the textual results of the \textsl{bg}.\textsl{max} function.}
\label{figure:browser}
\end{figure}

\textsl{AmpSim} has several parameters controlling simulation of 
amplification curves. $b.eff$ and $Cq$ are strongly connected. 
Thus changing one of them changes both values. $Cq$ is used to 
define an approximate Cq value. The expression ``approximate Cq'' value is 
used because the calculated Cq value varies depending on the preferred 
Cq quantification method (e.g., second derivative maximum ($SDM$) method, 
threshold method). \textsl{AmpSim} is used to simulate data with noise 
(based on \textsl{rnorm}, \emph{stats}), signal-to-noise ratios, 
photo-bleaching and other influences on a qPCR reaction. The following 
example illustrates the use of \textsl{AmpSim} 
(Figure~\ref{figure:AmpSim_random}).

<<AmpSim_random,warning=FALSE,message=FALSE,fig.show='hold',fig.cap=fig1_cap,fig.scap=fig1_scap,fig.width = 11, fig.height = 8>>=
# Draw an empty plot for 40 cycles with user defined parameters.

par(las = 0, bty = "n", oma = c(.5,.5,.5,.5))
plot(NA, NA, xlim = c(1,40), ylim = c(0,1.1), xlab = "Cycle", ylab = "RFU")
colors <- rainbow(8)

# Create eight amplification curves. The approximate Cqs are synthesized 
# as temporary Cqs by adding a random value to a starting Cq of 25. Note: 
# ``noise'' is set TRUE with a level of nnl = 0.03. This adds some scatter 
# to the amplification curves.

sim <- sapply(1L:8, function(i) {
  Cq.tmp <- 25 + rnorm(1) * 5
  
  tmp <- AmpSim(1:40, Cq = Cq.tmp, noise = TRUE, nnl = 0.03)
  lines(tmp, col = colors[i], lwd = 2)
  
  # Add the approximate Cq values to the plot
  text(3, 1 - i / 10, paste("Cq ", round(Cq.tmp, 2)), col = colors[i])
})
@

\textsl{AmpSim} was used to 
illustrate the use of \textsl{inder} (Figure~\ref{figure:SDM}), the 
\textsl{fixNA} (Figure~\ref{figure:fixNA}) and the 
\textsl{smoother} (Figure~\ref{figure:smoothing_intro}) functions.

\section{Single-blinded, randomized judging of amplification curves}
\label{sec:humanrater}

Humans show bias towards interpreting data for a particular outcome. A single-blinded 
and randomized experiment aims to reduce bias in the results. We developed the 
\textsl{humanrater} function to evaluate the quality of curve data (e.g., 
amplification or melting curve data) in a randomized, half-blinded manner. The 
function allows interactive rating of a curve for a certain characteristic. 
\textsl{humanrater} draws individual graphs of a curve and prompts an input 
field for the user. The application of this function are numerous (e.g., comparing the 
human rating and the rating of a machine or the rating of several individual 
experts). A list of designations to characterize the amplification curve can be 
specified. The names of elements can be specified (e.g., by short designations 
used during rating). Defaults are $y$ for ``yes'', $a$ for ``ambiguous'' and $n$ 
for ``no''. It is possible to supply longer or shorter designations for lists. In 
our example, we used \textsl{humanrater} in the \textbf{RKWard} GUI 
(Figure~\ref{figure:humanrater}). We aimed to characterize amplification curves 
which were randomly drawn from our simulated ``testdata'' dataset.


<<humanrater_gui,eval = FALSE>>=
# Create a set of data to be analyzed by humanrater.
# The function AmpSim creates amplification curves which follow a nearly
# optimal sigmoidal curve shape or just noise.

testdata <- data.frame(1:35,
                       AmpSim(Cq = 15, noise = TRUE)[, 2], 
                       AmpSim(Cq = 25, noise = TRUE)[, 2],
                       rnorm(35),
                       AmpSim(Cq = 35, noise = TRUE)[, 2],
                       rnorm(35),
                       AmpSim(Cq = 45, noise = TRUE)[, 2])

# Use testdata as input for humanrater and assign the results to the
# object human.test.
# check testdata for significance of amplification in two repeats.

human.test1 <- humanrater(testdata, repeats = 2)
@

\begin{figure}[ht]
\centering
\scalebox{0.43}{
\includegraphics{humanrater.png}
}
\caption{Application of \textsl{humanrater} in a working instance of \textbf{RKWard}. 
\textsl{humanrater} was used to analyze a row of amplification curves. 
\emph{(A)} All data are anonymous and can be randomized during the rating. The 
number of repeats for the rating and the categories (e.g., $y$ for ``yes'', $a$ 
for ``ambiguous'' and $n$ for ``no'') can be defined by the user. The function 
has an option to present the curves at random (default). \emph{(B)} The user 
gets as result a tabular output, including the result of each run and the 
conformity of the runs (see table in console).} 
\label{figure:humanrater}
\end{figure}

\section{Inspection and analysis of data for amplification curves}
\label{sec:MFIaggr}

The following section briefly describes $\emph{chipPCR}$ functions used for 
visualizing and analyzing data for amplification curves. The functions 
\textsl{MFIaggr} and \textsl{plotCurves} (Section~\ref{sec:plotCurves}) were 
developed for a rapid and convenient inspection of raw data. \textsl{MFIaggr} is 
a powerful analytical and graphical tool for fast multiple comparison of 
cycle-dependent signal dispersion and distribution. The continuous response 
variable $y$ is used to describe the relationships to $n$ continuous predictor 
variables $x_i$, where $i \in \{1, ..., n\}$. Use cases include the comparison 
of independent reaction vessels or the analysis of replicate experiments.

The idea is to analyze only a region of interest (ROI) from a dataset as defined 
by the parameter $llul$ (\textbf{l}ower \textbf{l}imit and \textbf{u}pper 
\textbf{l}imit). A ROI can be cycles or a time frame. \textsl{MFIaggr} is a 
relative of the \textsl{MFIerror} function from the \emph{MBmca} 
\citep{roediger_RJ_2013} package but allows a finer grained data analysis for specific 
parts of a plotted curve. \textsl{MFIaggr} returns an object of the class list 
with the columns ``Cycle'', ``Location'' (Mean, Median), ``Deviation'' (Standard 
Deviation, Median Absolute Deviation) and ``Coefficient of Variation''. If the 
option $rob$ is $TRUE$ the function calculates out the median and the median 
absolute deviation (MAD) instead of the mean and standard deviation. 
\textsl{MFIaggr} has The results for the ROI can be invoke by $@stats$. The 
output includes the mean, median, standard deviation (sd), median absolute 
deviation (mad), inter quartile range (IQR), medcouple (robust measure of 
skewness), skewness (Pearson's second skewness coefficient; 
$skewness~=~3~(mean(x)~-~median(x))~/~sd(x)$), signal-to-noise ratio (SNR), 
variance-to-mean ratio (VRM), number of missing values (NAs) and results from a 
linear fit of the ROI (intercept, slope, r.squared). We included the 
Breusch-Pagan test to test for heteroskedasticity in a linear regression model 
(see Section~\ref{sec:MFIaggr}). In our example we analyzed the raw fluorescence 
from 96 replicates of a qPCR experiment for the human gene \textit{vimentin}. 
The \textsl{MFIaggr} plot shows that the first ten cycles (noise) follow a 
normal distribution (Figure~\ref{figure:MFIaggr_intro}). In contrast, the 
analysis of all cycles expectedly shows a distribution, which significantly 
differs from a normal distribution (Figure~\ref{figure:MFIaggr_all}). Setting 
the option $CV = FALSE$ shows the relative standard deviation (RSD,~\%). The 
variance between the amplification curves of replicates should be low. Other 
results of \textsl{MFIaggr} include the density analysis ($@density$), the 
quantile ($@qqnorm.data$), and the results of the linear regression ($@lm.roi$) 
from the ROI. In particular, this function might be useful for quality 
management during the development of high-throughput technologies.

<<MFIaggr_intro,warning=FALSE,message=FALSE,fig.show='hold', fig.cap=fig23_cap,fig.scap=fig23_scap,fig.width = 11, fig.height = 8>>=
par(las = 0, bty = "n", cex.axis = 1.2, cex.lab = 1.2, 
    font = 2, cex.main = 1.2, oma = c(1,1,1,1))

plot(MFIaggr(VIMCFX96_60[, 1], VIMCFX96_60[, 2:ncol(VIMCFX96_60)], 
     llul = c(1,10)), CV = FALSE)

# plot(MFIaggr(VIMCFX96_60[, 1], VIMCFX96_60[, 2:ncol(VIMCFX96_60)], 
#      llul = c(1,40)), CV = FALSE)
@

The function can be used to compare two conditions of a qPCR experiment. In 
our example we tried to spot differences between a measurement during the 
annealing phase and the elongation phase 
(Figure~\ref{figure:MFIaggr_conditions}).

<<MFIaggr_conditions,warning=FALSE,message=FALSE,fig.show='hold', fig.cap=fig24_cap,fig.scap=fig24_scap,fig.width = 11, fig.height = 8>>=
par(las = 0, bty = "n", cex.axis = 1.2, cex.lab = 1.2, 
    font = 2, cex.main = 1.2, oma = c(1,1,1,1))

plot(x = MFIaggr(VIMCFX96_60, fluo = c(2L:10)), y = MFIaggr(VIMCFX96_69, fluo = c(2L:10)))
@

An analysis via the 
\emph{shiny} \textsl{MFIaggr.gui} app is shown in Figure~\ref{figure:MFIaggr2}.

\begin{figure}[ht]
\centering
\scalebox{0.43}{
\includegraphics{MFIaggr_gui.png}
}
\caption{Example of \emph{shiny} \textsl{MFIaggr}.\textsl{gui} app. \emph{(Top)} 
The plot of the \textsl{AmpSim}.\textsl{gui} is shown in a web browser 
(\textbf{Iceweasel}, v. 32.0) along with the parameters (left panel) and the code 
(``server.R'', ``ui.R'').} 
\label{figure:MFIaggr2}
\end{figure}



In our example we analyzed the raw fluorescence from 96 replicates 
($\emph{VIMCFX96\_60}$ dataset) of a qPCR experiment for the human gene 
\textit{vimentin}. The \textsl{MFIaggr} plot shows that the analysis of 
all cycles is non-normally distributed (Figure~\ref{figure:MFIaggr_all}). 

<<MFIaggr_all,fig.show='hold',fig.cap=fig12_cap,fig.scap=fig12_scap>>=
plot(MFIaggr(VIMCFX96_60[, 1], VIMCFX96_60[, 2:ncol(VIMCFX96_60)], 
             llul = c(1,40)), CV = FALSE)
@

\textsl{MFIaggr} analyzes the heteroskedasticity. Heteroskedasticity (“hetero” = 
different, “skedasis” = dispersion) is present if the variance (error term) is 
not constant. If the error terms do not have a constant variance, they are 
homoskedastic. Analysis of the heteroskedasticity gives insight into the 
characteristics of a system. In the following example we compared the 
$\emph{VIMCFX96\_60}$ and $\emph{VIMCFX96\_69}$ datasets both obtained from the 
same qPCR run in a Bio-Rad CFX96 during an annealing phase of $60\,^{\circ}\mathrm{C}$ and an elongation phase of $69\,^{\circ}\mathrm{C}$, respectively. The heteroskedasticity increased 
expectedly during the amplification reaction. The variance in the 
elongation phase (Figure~\ref{figure:MFIaggr_all}C and D) was lower than 
in the annealing phase (Figure~\ref{figure:MFIaggr_all}A and B). The 
heteroskedasticity was significant during the 
first 15 cycles at $60\,^{\circ}\mathrm{C}$ (Figure~\ref{figure:MFIaggr_all}A).

<<MFIaggr_heteroskedasticity,fig.show='hold',fig.cap=fig15_cap,fig.scap=fig15_scap,fig.width = 11, fig.height = 8>>=
par(mfrow = c(2,2), bty = "n")
# Create a helper function "hsk.test" to analyze the heteroskedasticity
# and the variance.
hsk.test <- function(x, y, llul = c(1,15), main = "") {
  res <- MFIaggr(x, y, llul = llul)
  head(res)
  plot(res[, 1], res[, 3]^2, xlab = "Cycle", 
       ylab = "Variance of refMFI", 
       xlim = llul, ylim = c(min(res[llul[1]:llul[2], 3]^2), 
                             max(res[llul[1]:llul[2], 3]^2)), main = main, 
       pch = 19, type = "b")
  abline(v = llul, col = "grey", lty = 2, lwd = 2)
  legend("top", paste0("Breusch-Pagan test p-value: \n", 
                       format(summary(res, print = FALSE)[14], digits = 2)), bty = "n")
}

hsk.test(VIMCFX96_60[, 1], VIMCFX96_60[, 2:ncol(VIMCFX96_60)], 
         llul = c(1,15),
         main = "ROI Cycle 1 to 15\nAnnealing phase") 
mtext("A", cex = 2, side = 3, adj = 0)

hsk.test(VIMCFX96_60[, 1], VIMCFX96_60[, 2:ncol(VIMCFX96_60)], llul = 
           c(1,40), main = "ROI Cycle 1 to 40\nAnnealing phase") 
mtext("B", cex = 2, side = 3, adj = 0)

hsk.test(VIMCFX96_69[, 1], VIMCFX96_69[, 2:ncol(VIMCFX96_69)], llul = 
           c(1,15), main = "ROI Cycle 1 to 15\nElongation phase") 
mtext("C", cex = 2, side = 3, adj = 0)

hsk.test(VIMCFX96_69[, 1], VIMCFX96_69[, 2:ncol(VIMCFX96_69)], llul = 
           c(1,40), main = "ROI Cycle 1 to 40\nElongation phase") 
mtext("D", cex = 2, side = 3, adj = 0)
@

\subsection{Data overview - \textsl{plotCurves}}
\label{sec:plotCurves}

\textsl{plotCurves} visualizes many curves on one plot in separate cells
allowing quick assessment of experiments (Figure~\ref{figure:plotCurves}). In 
addition to this, \textsl{plotCurves} has an option to run an unsupervised 
\textsl{CPP} pre-processing step on the raw data. This smooths the 
data (Savitzky-Golay smoothing), removes missing values (spline 
interpolation by default) and performs a background subtraction (base-lining 
to zero). \textsl{plotCurves} has a colored indicator for rapid visualization of 
dataset with potentially problematic amplification curves. The plot output is arranged 
in a table-like fashion, where each curve is presented in a different cell.


<<plotCurves,fig.show='hold',fig.cap=fig14_cap,fig.scap=fig14_scap,warning=FALSE,fig.width = 11, fig.height = 8>>=
y <- VIMCFX96_60[, 2L:9]
# Introduce some missing values.
y[c(10, 22, 3, 25, 26, 15, 27, 23, 4), c(5, 7, 4, 2, 1)] <- NA

# Show plot with raw data and missing values (black line) and show 
# plots with pre-processed data and imputed missing values (red line).
plotCurves(VIMCFX96_60[, 1], y, nrow = 2, type = "l", CPP = TRUE)
@


\section{Proposed workflow}

In the previous section we showed different methods for investigating specific 
properties of the measured data. Next, we focused on pre-processing methods of the 
$\emph{chipPCR}$ package. Here, we wish to state only
function names and give some information on their working principle. Details will 
be explained in the subsequent sections in order to avoid confusion of the reader. We show the application of the \textsl{CPP} function, as a proposed workflow for 
customized pre-processor functions. Data were taken from the $\emph{VIMCFX96\_60}$ data 
set. This dataset was measured with a Bio-Rad CFX96 thermo-cycler with 96 
replicates (see $\emph{chipPCR}$ manual for experimental details).

The \textsl{CPP} function is a wrapper for pre-processing algorithms, which 
includes normalization, background subtraction, outlier removal using the (\textsl{fixNA} 
function) in the background range and tests for positive amplification 
reactions. \textsl{CPP} uses the \textsl{bg}.\textsl{max} function to automatically estimate the start of the amplification process. The background range is 
often noisy, which makes it hard to determine a meaningful background value. 
Therefore, \textsl{CPP} can optionally remove outliers by finding the value with 
the smallest and largest difference from the mean as provided by the 
$rm.outlier$ function from the $\emph{outlier}$ package \citep{Komsta_2011}. 
$rm.outlier$ detects these outliers by a simple rule without statistical testing 
and replaces them by the sample mean. Outliers herein refer to the smallest 
and largest value, which has maximum difference from the sample mean. The slope 
of the background range is often unequal to zero. By setting the parameter 
$trans$ it is possible to apply a simple correction of the slope. This slope correction can either be done by 
a robust linear regression that computes MM-type regression estimators, or by a 
nonparametric rank-based estimator or a standard linear regression model. 
\textsl{CPP} uses by default a robust linear regression (MM-type estimator) as 
integrated in the \textsl{lm}.\textsl{coefs} function. A defined range of the 
amplification curve (typically the background range) is used to extrapolate the 
linear trend over the entire dataset. However, this step has to be performed 
with caution since this operation affects the AE. The background is assumed to 
be constant for the entire measurement. Caution is needed when using $trans$ 
with time series (see \textsl{lm} from the $\emph{stats}$ package for details). 
Additionally, all data are normalized between the minimum and maximum. This is 
taken care of by the \textsl{normalizer} function. Smoothing of the data is 
finally done based on an instance of the \textsl{smoother} function. By default, 
a Savitsky-Golay filter was used to smooth the data. The following code is a 
representative example for the use of \textsl{CPP} 
(Figure~\ref{figure:workflow}). Note that warnings in following code chunks were 
suppressed.

<<workflow,fig.show='hold',fig.cap=fig4_cap,fig.scap=fig4_scap,warning=FALSE,fig.width = 11, fig.height = 8>>=
layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE), respect = TRUE)

par(las = 0, bty = "n", oma = c(.5,.5,.5,.5))

th.cyc.raw <- apply(VIMCFX96_60[, -1], 2, function(i) {
  th.cyc(VIMCFX96_60[, 1], i, r = 2575)[1,1]})

res.CPP <- apply(VIMCFX96_60[, -1], 2, function(i) {
  CPP(VIMCFX96_60[, 1], i, trans = TRUE, 
      method.norm = "minm")[["y.norm"]]})

th.cyc.CPP <- apply(res.CPP, 2, function(i) {
  th.cyc(VIMCFX96_60[, 1], i, r = 0.1)[1,1]})

matplot(VIMCFX96_60[, -1], type = "l", pch = 19, col = 1, lty = 1, 
        xlab = "Cycle", ylab = "Raw fluorescence", main = "Raw")
abline(h = 2575, lty = 2)
mtext("A", cex = 1.2, side = 3, adj = 0, font = 2)

matplot(res.CPP, type = "l", pch = 19, col = 1, lty = 1, xlab = "Cycle", 
        ylab = "Fluorescence", main = "CPP")
abline(h = 0.1, lty = 2)
mtext("B", cex = 1.2, side = 3, adj = 0, font = 2)

boxplot(data.frame(Raw = th.cyc.raw, CPP = th.cyc.CPP), ylab = "Cq (Ct)", 
        notch = TRUE)
mtext("C", cex = 1.2, side = 3, adj = 0, font = 2)
@

\section{Imputation of missing values in amplification curve data - \textsl{fixNA}}
\label{sec:fixNA}

Experimental technologies may produce missing values (NA) at 
random due to sensor drop-outs or other technical difficulties. Upon ecnountering a missing value, many analytical functions stop to progress or discard entire datasets. Such a conduct is reasonable 
in cases where the data structure is unknown. However, in the case of 
data used for plotting amplification curves, it is justified to impute NAs because the structure 
generally resembles an S-shaped curve. Standard approaches include substitution 
with most frequent values, mean value imputation, last value carried forward, 
bootstrapping, or substitution by correlation with replicate measurements 
\citep{Harrell_2001}. In the case of amplification curves other approaches are 
favorable. In particular, the transition phases (e.g., background phase to 
exponential phase) are potentially prone to bias.

The NAs may be caused by detector problems, acquisition error or 
other assorted problems. There are different ways to handle missing 
values. One approach is to ignore NAs, which is generally acceptable. 
However, in case of further calculation it is often necessary to handle 
cases of missing values in a way that the next calculation steps can be 
performed. Missing values can be eliminated by a imputation, which 
encompasses various approaches. This includes calculating a location parameter 
(e.g., mean, median) or other significant values (e.g., minimum, maximum, modus) 
of a data column. However, in non-linear processes such as amplification 
processes its is reasonable to estimate the missing values from a trend. 

The function \textsl{fixNA} imputes missing values in a single column of data 
(response) either by linear approximation or an 
approximation by cubic splines (default) (Figure~\ref{figure:fixNA}). Other 
smoothing functions such as the Savitzky-Golay smoothing filter have the 
intrinsic capability to remove missing values \citep{Savitzky_1964, 
Eilers_2003}. However, such functionality is not yet implemented in this package. This linear 
approach is useful but may be problematic on the phases other than background or 
plateau phases of an amplification reaction. The parameter $spline$ of 
\textsl{fixNA} enables a trend estimation on cubic splines and may be more 
appropriate in most scenarios.

The $\emph{reps384}$ dataset from the $\emph{qpcR}$ package 
\citep{Spiess_qpcR_2014} was used to compare the influence of imputation on 
real-world data. The experimental details have been described in 
\citep{ruijter_2013}. The dataset consists of 379 replicate amplification curves 
(see documentation of the $\emph{qpcR}$ package for details). Our \textit{in-silico} 
experiment was designed as followed: Either one or three missing values were 
artificially added to each amplification curve at random positions. We separated 
the amplification curve into three different regions (``Linear phase'' (cycle 1 
-- 10), ``Exponential phase'' (cycle 11 -- 34) and ``Plateau phase'' (cycle 34 
-- 40), Figure~\ref{figure:fixNA_data}) and investigated the impact on the qPCR 
parameters ``Cq (SDM, Cy0)'' and curve background parameter $bg$. The background 
was calculated as mean and standard deviation. The performance of the imputed 
models was analyzed with the goodness-of-fit in this region by the commonly used 
normalized root-mean-squared-error (NRMSE). We compared the imputation by 
``linear approximation'' ($fixNA(x, y, spline = FALSE)$) and ``spline 
approximation'' ($fixNA(x, y, spline = TRUE)$). Our results show that imputation 
with the spline method worked reliably and introduced no significant bias to all 
the investigated parameters. The 
Kruskal-Wallis rank sum test (\textsl{kruskal.test}, $\emph{stats}$ package 
\citep{RCT_2013}) was used to compare the linear and spline-based imputation.

Our \textit{in-silico} experiments showed that cubic spline interpolation yielded the most probable 
values and therefore led to the least effect on tested statistical parameters 
(Cq, background signal, Pearson correlation coefficient) on the exponential phase 
and is therefore the recommended approach to remove missing values 
(Figure~\ref{figure:fixNA_data}). We observed no significant bias using cubic spline 
interpolation (Table~\ref{table:fixNA}). The performance of \textsl{fixNA} using 
cubic splines was better than a linear interpolation (Figure~\ref{figure:fixNA}). 
However, the linear approximation might be applicable in measurements with high 
sampling rates (e.g., isothermal amplification) (not shown). Any method requires a 
minimum number of data points as foundation for a meaningful imputation. 
\textsl{fixNA} attempts to take care of such pitfalls. By rule of thumbs we 
determined that the number of missing elements in relation to the total number 
of elements should not exceed 30~\%. In a case where more than 30~\% of all values are NAs, the  \textsl{fixNA} function gives a warning.

Our results for the given experimental setting support the following statements. 
(I) The imputation of missing values by spline interpolation and linear methods 
introduce no significant bias on the tested parameters ``Cq (SDM, Cy0)'', 
``bg'' and the accompanied quality measure NRMSE (Table~\ref{table:fixNA}). (II) 
We found that the difference between the linear and spline imputation method is 
negligible (p~1).

<<fixNA_data,fig.show='hold',fig.cap=fig20_cap,fig.scap=fig20_scap,message=FALSE,fig.width = 11, fig.height = 8>>=
library(qpcR)
library(chipPCR)
cols <- adjustcolor(2:4, 0.6)
plot(NA, NA, xlim = c(1,45), ylim = c(min(reps384[, -1]), max(reps384[, -1])), 
     col = 1, pch = 19, type = "b", xlab = "Cycle", ylab = "Fluorescence")
rect(0.8, min(reps384[, -1]), 10.2, max(reps384[, -1]), border = NA, col = cols[1])
rect(10.8, min(reps384[, -1]), 33.2, max(reps384[, -1]), border = NA, col = cols[2])
rect(33.8, min(reps384[, -1]), 45, max(reps384[, -1]), border = NA, col = cols[3])

apply(reps384[, -1], 2, function(i) lines(reps384[, 1], i))
@


<<fixNA,fig.show='hold',fig.cap=fig8_cap,fig.scap=fig8_scap>>=
# Simulation of an ideal amplification curve with 40 cycles
# The other parameter of the AmpSim function are identical to
# the default.

res <- AmpSim(cyc = 1:40)

# Introduce a missing value (cycle 18) in the transition between 
# the background and the exponential phase.

res.NA <- res
res.NA[18, 2] <- NA

# Helper function to highlight the position of the missing value.
abliner <- function(x1 = 17.5, x2 = 18.5, y1 = 0.09, y2 = 0.14) {
  abline(v = c(x1, x2), col = "red")
  abline(h = c(y1, y2), col = "red")
}

par(las = 0, mfrow = c(2,2), bty = "n")
plot(res, xlab = "Cycles", ylab = "refMFI", type = "b", pch = 20, 
     main = "Without NA")
abliner()
mtext("A", cex = 1.2, side = 3, adj = 0, font = 2)
res.NA.linear <- fixNA(res.NA[, 1], res.NA[, 2], spline = FALSE, 
                       verbose = FALSE)

plot(res.NA, xlab = "Cycles", ylab = "refMFI", type = "b", pch = 20, 
     main = "With NA during transition")
abliner()
mtext("B", cex = 1.2, side = 3, adj = 0, font = 2)

res.NA.spline <- fixNA(res.NA[, 1], res.NA[, 2], spline = TRUE, 
                       verbose = FALSE)

plot(res.NA.linear, xlab = "Cycles", ylab = "refMFI", type = "b", 
     pch = 20, main = "Linear imputed\n NA")
abliner()
mtext("C", cex = 1.2, side = 3, adj = 0, font = 2)

plot(res.NA.spline, xlab = "Cycles", ylab = "refMFI", type = "b", 
     pch = 20, main = "Spline imputed\n NA")
abliner()
mtext("D", cex = 1.2, side = 3, adj = 0, font = 2)
par(mfrow = c(1,1))
@

\clearpage


<<fixNA_simulation,echo=TRUE,eval=FALSE>>=
# Evaluation of fixNA by introducing missing values and comparing efficiency measures
# for imputed and raw data

library(qpcR)
library(chipPCR)

linear_range <- 1L:10
exponential_range <- 11L:33
plateau_range <- 34L:45

raw.eff <- t(sapply(2L:ncol(reps384), function(i) {
  fit.raw <- pcrfit(reps384, cyc = 1, fluo = i)
  c(cpD2.raw = efficiency(fit.raw, type = "cpD2", plot = FALSE)[["cpD2"]],
    Cy0.raw = efficiency(fit.raw, type = "Cy0", plot = FALSE)[["Cy0"]],
    background = mean(reps384[1L:5, i]))
}))

#calculate raw efficiency measures
raw.df <- as.vector(apply(raw.eff, 2, function(i)
  c(mean(i), sd(i))))

#reformat raw data
raw.df <- cbind(data.frame("none", 0, "-"), t(raw.df))
colnames(raw.df) <- c("Imputation", "NA.number", "Phase", 
                      "cpD2.mean", "cpD2.sd",
                      "Cy0.mean", "Cy0.sd",
                      "background.mean", "background.sd")

#helper function introducing NA into data, fixing them and calculating efficiency 
gen.fix <- function(number_points, phase, spline){
  #copy raw data
  temp_dat <- reps384[, -1]
  
  #introduce NA(s)
  for(i in 1L:ncol(temp_dat))
    temp_dat[sample(phase, 1), i] <- NA
  
  #impute missing values using fixNA
  fixed_dat <- sapply(1L:ncol(temp_dat), function(i)
    fixNA(reps384[, 1], temp_dat[, i], spline = spline))
  
  #compute efficiency measures
  t(sapply(2L:ncol(fixed_dat), function(i) {
    fit.fix <- pcrfit(cbind(reps384[, 1], fixed_dat), cyc = 1, fluo = i)
    c(cpD2.fix = efficiency(fit.fix, type = "cpD2", plot = FALSE)[["cpD2"]],
      Cy0.fix = efficiency(fit.fix, type = "Cy0", plot = FALSE)[["Cy0"]],
      background = mean(fixed_dat[1L:5, i]))
  }))
}

#calculate results for linear imputation
res.linear <- lapply(c(1, 3), function(number_of_points)
  lapply(list(linear_range, exponential_range, plateau_range), function(phase)
    apply(gen.fix(number_of_points, phase, FALSE), 2, function(i)
      c(mean(i), sd(i)))))
linear.df <- data.frame(num = unlist(lapply(c(1, 3), rep, 3)), 
                        region = rep(c("linear", "exponential", "plateau"), 2),
                        do.call(rbind, lapply(unlist(res.linear, recursive = FALSE), as.vector))) 


#calculate results for spline imputation
res.spline <- lapply(c(1, 3), function(number_of_points)
  lapply(list(linear_range, exponential_range, plateau_range), function(phase)
    apply(gen.fix(number_of_points, phase, TRUE), 2, function(i)
      c(mean(i), sd(i)))))
spline.df <- data.frame(num = unlist(lapply(c(1, 3), rep, 3)), 
                        region = rep(c("linear", "exponential", "plateau"), 2),
                        do.call(rbind, lapply(unlist(res.linear, recursive = FALSE), as.vector)))

#join and format results
sim.res <- cbind(unlist(lapply(c("linear", "spline"), rep, 6)), 
                 rbind(linear.df, spline.df))
colnames(sim.res) <- c("Imputation", "NA.number", "Phase", 
                       "cpD2.mean", "cpD2.sd",
                       "Cy0.mean", "Cy0.sd",
                       "background.mean", "background.sd")

#join simulation results with results for raw data
fixNA.evaluation <- rbind(raw.df, sim.res)
xtable(fixNA.evaluation[, c(1L:7)], digit = 4)

xtable(fixNA.evaluation[, c(1L:3, 8L:9)], digit = 4)
@


\begin{table}[ht]
\centering
\begin{tabular}{rlrlrrrr}
  \hline
 & Imputation & NA.number & phase & cpD2.mean & cpD2.sd & Cy0.mean & Cy0.sd \\ 
  \hline
1 & none & 0.0000 & - & 19.2931 & 0.1464 & 11.1806 & 1.1326 \\ 
  2 & linear & 1.0000 & linear & 19.2932 & 0.1466 & 11.1801 & 1.1337 \\ 
  3 & linear & 1.0000 & exponential & 19.2875 & 0.1472 & 11.1631 & 1.1351 \\ 
  4 & linear & 1.0000 & plateau & 19.2932 & 0.1469 & 11.1814 & 1.1336 \\ 
  5 & linear & 3.0000 & linear & 19.2929 & 0.1468 & 11.1803 & 1.1339 \\ 
  6 & linear & 3.0000 & exponential & 19.2878 & 0.1465 & 11.1636 & 1.1351 \\ 
  7 & linear & 3.0000 & plateau & 19.2932 & 0.1466 & 11.1815 & 1.1335 \\ 
  8 & spline & 1.0000 & linear & 19.2932 & 0.1466 & 11.1801 & 1.1337 \\ 
  9 & spline & 1.0000 & exponential & 19.2875 & 0.1472 & 11.1631 & 1.1351 \\ 
  10 & spline & 1.0000 & plateau & 19.2932 & 0.1469 & 11.1814 & 1.1336 \\ 
  11 & spline & 3.0000 & linear & 19.2929 & 0.1468 & 11.1803 & 1.1339 \\ 
  12 & spline & 3.0000 & exponential & 19.2878 & 0.1465 & 11.1636 & 1.1351 \\ 
  13 & spline & 3.0000 & plateau & 19.2932 & 0.1466 & 11.1815 & 1.1335 \\ 
   \hline
\end{tabular}
\caption{Results of fixNA data
       imputation, part 1: cpD2 and Cy0. NA column determines how many NA values were 
       introduced in a given phase.} 
\end{table}

\begin{table}[ht]
\centering
\begin{tabular}{rlrlrr}
  \hline
 & Imputation & NA.number & phase & background.mean & background.sd \\ 
  \hline
1 & none & 0.0000 & - & 4567.5648 & 184.7082 \\ 
  2 & linear & 1.0000 & linear & 4568.2296 & 184.7895 \\ 
  3 & linear & 1.0000 & exponential & 4568.0775 & 184.6828 \\ 
  4 & linear & 1.0000 & plateau & 4568.0775 & 184.6828 \\ 
  5 & linear & 3.0000 & linear & 4568.2859 & 184.5738 \\ 
  6 & linear & 3.0000 & exponential & 4568.0775 & 184.6828 \\ 
  7 & linear & 3.0000 & plateau & 4568.0775 & 184.6828 \\ 
  8 & spline & 1.0000 & linear & 4568.2296 & 184.7895 \\ 
  9 & spline & 1.0000 & exponential & 4568.0775 & 184.6828 \\ 
  10 & spline & 1.0000 & plateau & 4568.0775 & 184.6828 \\ 
  11 & spline & 3.0000 & linear & 4568.2859 & 184.5738 \\ 
  12 & spline & 3.0000 & exponential & 4568.0775 & 184.6828 \\ 
  13 & spline & 3.0000 & plateau & 4568.0775 & 184.6828 \\ 
   \hline
\end{tabular}
\caption{Results of fixNA data
       imputation, part 2: background. NA column determines how many NA values were 
       introduced in a given phase.} 
\end{table}


\begin{figure}[ht]
\centering
\scalebox{0.6}{
\includegraphics{fixNA_test1.png}
}
\caption{Imputation of missing values in data for plotting amplification curves by splines and 
linear trends. Missing values were artificially introduced into the 
$\emph{reps384}$ dataset from the $\emph{qpcR}$ package and imputed by the 
\textsl{fixNA} function. We found no significant difference between raw data and 
data with imputed missing values. $R^2$ and correlation coefficients of curves 
were close to 1 with p-value $< 10^{-6}$.}
\label{figure:test_fixNA}
\end{figure}


\section{Smoothing and filtering}
\label{sec:smoothing_filterring}

For data presentation it is often 
useful to smooth or filter the data. Smoothing and filtering are different 
approaches with similar goals. Both pre-process an input signal as output 
for subsequent analysis steps. Filtering uses methods of signal processing and 
takes a data input and applies a function to form an output. Smoothing in contrast 
uses statistical approaches, like local regression models (e.g., least square 
estimate) or cubic splines. We developed the 
\textsl{smoother} function, which is a wrapper for smoother functions and 
filters commonly used to process data for amplification curves. \textsl{smoother} 
inherited traits (Table~\ref{table:tab1}) of the parent functions 
\citep{spiess_impact_2014}. However, the functionality of \textsl{smoother} 
greatly outgrows applications only in amplification curve analysis. 
Incorporating most of the best proven algorithms, we offer the user a powerful 
tool to access the methods while minimizing the drawback of learning the syntax of 
specific functions. \textsl{smoother} was enhanced by the functionality of 
\textsl{fixNA} and \textsl{CPP}. Figure~\ref{figure:smoothing_intro} shows 
results of the \textsl{smoother} function on data used for amplification curves.

<<smoothing_intro,fig.show='hold',fig.cap=fig22_cap,fig.scap=fig22_scap,warning=FALSE,fig.width = 11, fig.height = 8>>=
# Simulate and amplification curve with the AmpSim function
tmp <- AmpSim(cyc = 1:35, bl = 0)

par(las = 0, bty = "n", cex.axis = 1.5, cex.lab = 1.5, 
    font = 2, cex.main = 1.8, oma = c(1,1,1,1), fig = c(0,1,0.55,1))
plot(tmp, type = "b", col = 1, pch = 20, xlab = "", ylab = "RFU", 
      main = "Raw data")
mtext("A", cex = 2, side = 3, adj = 0, font = 2)

# Apply all (parameter method = "all") smoothers/filter with the default 
# setting to the amplification curve of the object tmp. Smoothers / Filters:
# Savitzky-Golay smoothing filter
# locally-weighted polynomial regression
# moving average, windowsize 3
# cubic spline smooth
# standard cubic spline smooth
# Friedman's SuperSmoother
# weighted Whittaker smoothing with first order finite difference penalty
# weighted Whittaker smoothing with a second order finite difference penalty
res <- smoother(tmp[, 1], tmp[, 2], method = "all", CPP = FALSE)

# Calculate the difference between the ideal curve (tmp) and the smoothed curves
# (res) and assign the results to the object res.out
res.out <- cbind(cycle = tmp[, 1], tmp[, 2] - res)

# Plot the smoothed curves
par(fig = c(0,1,0,0.65), new = TRUE)
plot(NA, NA, type = "b", col = 2, pch = 15, xlim = c(1,35), ylim = c(-0.1,0.1),
     xlab = "Cycle", ylab = "delta refMFI (raw - smoothed)",
     main = "Smoothed / Filtered data")
     
mtext("B", cex = 2, side = 3, adj = 0, font = 2) 
legend(1.5, 0.1, ncol = 2, colnames(res.out[, 2:9]), pch = 15:22, 
       lwd = 2, col = c(2:9))

# Plot the results.
tmp <- sapply(2:9, function(i) {
      lines(res.out[, 1], res.out[, i], type = "b", col = i, pch = i + 13)
     }
)
@

The presence of noise may cause many false estimates for the $FDM$ (first 
derivative maximum) and $SDM$ (second derivative maximum). To reduce noise, it 
is possible to smooth the first derivative of the amplification curve. Many 
methods integrated the moving average as a first pre-processing step (e.g., 
\citep{shain_2008}). The moving average filter is a linear filter, which 
sequentially replaces data points with the average of the neighbor data points. 
The average is calculated from a defined span (``window'') of odd count (e.g., 
3, 5). The ``average'' herein may refer to the arithmetic mean, the median, the 
geometric or the exponential mean. The \textsl{smoother} function uses 
exclusively the arithmetic mean. The moving average is intuitive and easy to 
implement. However, such processed data lags behind a trend and ignores rapid 
changes, leading to a forerun of few cycles \citep{spiess_impact_2014}. This is 
in particular problematic during the exponential phase. Splines apply non-parametric 
regression by local cubic polynomials between knot points \citep{Nie_2012}. 
Other examples for smoothers include the Savitzky-Golay smoothing filter, Friedman's 
SuperSmoother, and the Weighted Whittaker smoother (see the \textsl{smoother} 
function for details).

Selected smoothing and filtering functions (e.g., Savitsky-Golay filter) 
assume uniform (equally spaced) sampling. The function 
\textsl{smoother} and \textsl{CPP} (inherited from \textsl{smoother}) give a 
warning in such cases. It is recommended to pre-process the data to have equally 
spaced values. The \textsl{smoother} function enables users to tune the behavior of 
the chosen smoothing algorithm by using nearly all parameters available in 
called subroutines and at the same time uniformizes input and output. It should be 
noted that smoothing may alter the curve shape and thus lead to artificial 
results. Smoothed data are easier to read but introduce a bias to the 
pre-processed data. Therefore, the prime use of smoothers is in processing data for 
visualization purposes. However, it is not recommended to smooth unsupervised signals prior to statistical procedures (e.g., least-squares curve 
fitting). All smoothing algorithms are ``lossy'' to a certain extent and may 
change the curve shape significantly. In particular, the residual evaluation of 
a fit may lead to false prediction, because noise after smoothing may be 
mistaken for signal. Signals obtained after curve smoothing can be used to locate peaks, but such a procedure should be used cautiously for measuring peaks \citep{spiess_impact_2014}.

\begin{table*}[b]
\label{table:tab1}
\caption{ Smoothing and filtering methods of the $\emph{chipPCR}$ package. The 
parameter $lowess$ (locally-weighted polynomial regression, LOWESS) can be 
adjusted by the parameters $f$ and $iter$. The parameter $mova$ (moving average) 
can be adjusted by the windowsize parameter $movaww$. The parameter $savgol$ 
(Savitzky-Golay smoothing filter) can be adjusted by the parameter $p$. The 
parameter $smooth$ (cubic spline smooth) can be adjusted by the parameter 
$df.fact$. A $df.fact$ value of 1 smooths the data least while a value of 0.5 
smooths the curve most. The parameter $spline$ (standard cubic spline 
smooth) has no additional parameter. The parameter $supsmu$ (Friedman's 
SuperSmoother) can be adjusted by the parameter \textsl{span}. The parameter 
$whit1$ (first order finite difference penalty) and $whit2$ (second order finite 
difference penalty) for Weighted Whittaker smoothing filter can be adjusted by 
the parameter $lambda$. For further details on the smoothers refers to the 
documentation of their parent functions.
}
\begin{center}
\begin{tabular}{lccc}
  \hline
  Method & Parameter & value & Parent\\ \hline
  LOWESS & $lowess$ & \textit{f} & \textsl{lowess}, \emph{stats} \\
  Cubic spline & $smooth$ & \textit{df.fact} & \textsl{smooth.spline}, 
\emph{stats} \\
  Interpolating Splines & $spline$ & - & \textsl{spline}, \emph{stats} \\
  Friedman's ``super smoother'' & $supsmu$ & \textit{span} & \textsl{supsmu}, 
\emph{stats}\\
  Savitsky-Golay & $savgol$ & - & \textsl{sgolayfilt}, \emph{signal} \\
  Moving Average & $mova$ & \textit{movaww} (3, 5, ...) & \textsl{filter}, 
\emph{stats} \\
  Whittaker & $whit1$, $whit2$ & \textit{lambda} & \textsl{whit1}, 
\textsl{whit2},
\emph{ptw}  \\
  All smoother & $all$ & defaults &
  \\ \hline
\end{tabular}
\end{center}
\end{table*}


\section{\textsl{bg}.\textsl{max} - a function to estimate the start and end of an amplification reaction}
\label{sec:bgmax}

The following paragraphs describe methods used in literature to detect the 
background range of amplification curves. Background range herein refers to a 
level of reporter fluorescence signal measured before any specific amplification is detectable. 
The raw data (e.g., fluorescence intensity) measured after each step (cycle or 
time point) follow a non-linear progress. Currently none of them is implemented 
in \textbf{R} function. The easiest way to classify them is by the extend of 
assumptions made before applying of a method.

The simplest approach is to treat 
the background fluorescence as a value constant during the whole amplification 
reaction. In this case the noise could be approximated as the mean or median of 
fluorescence values in the lag phase \citep{frank_2009} or their standard deviations 
\citep{peirson_2003}. The more sophisticated way of approximating constant 
background fluorescence requires optimizing its value to achieve linearity of 
the model fit on the semi-logarithmic plot in the log-linear phase 
\citep{frank_2009}. The later procedure is greatly enhanced by performing 
further computations only on a subset of consecutive measurements for which 
calculated efficiencies have the lowest variance. Other methods invoke the 
assumption that background fluorescence is a constant value and instead describe 
it as a function of the cycle number. The algorithm in \textbf{SoFar} 
\citep{wilhelm_2003} fits a nonlinear saturation function to points measured
before the start of the exponential growth phase. Parameters of the saturation 
function are chosen to minimize the sum of squared residuals of the fitted 
function. Then the value of saturation function is calculated for all data 
points and subtracted from measured values giving corrected values of 
fluorescence, which are used in next calculations. 

Some approaches make even less assumptions regarding the form of the background 
noise. The taking-difference linear regression method uses the premise that 
changes of fluorescence between subsequent cycles could exclusively be caused by the 
amplification of the product \citep{rao_2013}. Thus, the corrected values are 
calculated by simply subtracting the fluorescence value in the former cycle 
from fluorescence in the latter. The real fluorescence 
value in the first cycle is unknown since the number of cycles used for 
the computation is reduced by one. 

The \textbf{Real-Time PCR Miner} algorithm is nearly assumption-free 
\citep{zhao_2005}. The principle is that background 
fluorescence is similar in the small groups of subsequent measurements. So the 
first step of the algorithm is division of subsequent measurement points 
belonging to the exponential phase of amplification in at least four-elemental 
groups. For each set of points a pair of the estimate of the 
efficiency and the significance of model representing the relation between the 
fluorescence value and the cycle number is calculated. The estimates with the highest 
significance are the most influential in the computation of the final 
efficiency. 

Finding the beginning of the lag phase and end of plateau phase is important for 
the goodness-of-fit for both exponential-phase-only and S-shaped models. There 
are two strategies. The first narrows the area of the search to the neighborhood 
of their theoretical values determined by a fitted model for the amplification 
reaction. To this group belongs \textbf{SoFar} \citep{wilhelm_2003}. The algorithm looks 
for the start and the end of the exponential phase near the second derivatives 
of the function representing the relation between logarithm of the fluorescence 
and the cycle number. The available correction guarantees that the start of 
amplification has a higher value than background noise. A very similar procedure 
is implemented in \textbf{Real-Time PCR Miner} \citep{zhao_2005}, where background noise 
is used as parameter in implemented models to calculate a theoretical start of the amplification process. The end of amplification process is detected 
by calculating the third derivative of an implemented S-shaped model. The second 
approach does not require theoretical values. A very intuitive solution, 
designated take-off point, by \citet{tichopad_2003} 
describes the lag phase using a linear function. Random deviations are taken 
into account as standardized residuals. The method starts with a fitting of a 
linear function to the first three measurement points. If none of the residuals is 
statistically identified as an outlier, the algorithm fits a new linear 
model to the first four measurement points and so on. The procedure stops when 
two last points are designated as outliers. The first of aforementioned outliers 
defines the end of lag phase. It is worth noting that this algorithm is 
versatile enough to detect the beginning of the plateau phase.

The algorithm of \textsl{bg}.\textsl{max} is based on the assumption that the 
signal difference of successive cycles in the linear ground phase is 
approximately constant. The signal drastically changes during transition into the early exponential phase. First data are smoothed by the Friedman's 'super 
smoother' as found in ``supsmu''. Thereof the approximate first and second 
derivative are calculated by a five-point stencil \textsl{inder}. The difference 
of cycles at the maxima of the first and second approximate derivative and a 
correction factor are used to estimate the range before the exponential phase. 
This simple function finds the background range without modeling the function. 
The start of the background range is defined be a ``fixed'' value. Since many 
signals tend to overshoot in the first cycles, a default value of 2 (for qPCR) is 
chosen. \textsl{bg}.\textsl{max} tries also to estimate the end of an 
amplification reaction (Figure~\ref{figure:bgmax}). See section 
\textsl{bg}.\textsl{max} ``Details'' of the $\emph{chipPCR}$ manual for further 
details. This function is a rational basis for trimming of unwanted data.

<<bgmax,fig.cap=fig18_cap,fig.scap=fig18_scap,fig.show='hold',message=FALSE,results='hide',fig.width = 11, fig.height = 8>>=
par(las = 0, mfrow = c(2,1), bty = "n", oma = c(.5,.5,.5,.5))

res <- AmpSim(cyc = 1:40, Cq = 25)
plot(res, xlim = c(1,40), ylim = c(-0.1,1), xlab = "Cycles", 
     ylab = "refMFI", 
     main = "Background Range Estimation\n in the Absence of Noise", 
     type = "b", pch = 20)
background <- bg.max(res[, 1], res[, 2])
mtext("A", cex = 2, side = 3, adj = 0, font = 2)

points(background[, 3], col = "red", type = "b", pch = 20)
points(background[, 4], col = "blue", type = "b", pch = 20)
abline(v = background@bg.start)
text(background@bg.start, 0.2, "Background start", pos = 4)
abline(v = background@bg.stop, col = "blue")
text(background@bg.stop, 0.25, "Background stop", pos = 4, 
     col = "blue")
abline(v = background@amp.stop, col = "green")
text(background@amp.stop, 0.3, "Plateau transition", pos = 4, 
     col = "green")
legend(4, 1, c("Raw data", "First derivative", "Second derivative"), 
       pch = rep(20,3), col = c(1,2,4), bty = "n")

res <- AmpSim(cyc = 1:40, Cq = 25, noise = TRUE)
plot(res, xlim = c(1,40), ylim = c(-0.1,1), xlab = "Cycles", 
     ylab = "refMFI", 
     main = "Background Range Estimation\n in the Presence of Noise", 
     type = "b", pch = 20)
mtext("B", cex = 2, side = 3, adj = 0, font = 2)
background <- bg.max(res[, 1], res[, 2])

points(background[, 3], col = "red", type = "b", pch = 20)
points(background[, 4], col = "blue", type = "b", pch = 20)
abline(v = background@bg.start)
text(background@bg.start, 0.2, "Background start", pos = 4)
abline(v = background@bg.stop, col = "blue")
text(background@bg.stop, 0.25, "Background stop", pos = 4, col = "blue")
abline(v = background@amp.stop, col = "green")
text(background@amp.stop, 0.3, "Plateau transition", pos = 4, col = 
       "green")
legend(4, 1, c("Raw data", "First derivative", "Second derivative"), 
       pch = rep(20,3), col = c(1,2,4), bty = "n")
par(mfrow = c(1,1))
@

The \textsl{bg}.\textsl{max} algorithm was used to analyze amplification 
data from a capillary convective PCR ($\emph{capillaryPCR}$ dataset, 
$\emph{chipPCR}$ package). The raw data (Figure~\ref{figure:bgmax_ccPCR}~A) and 
pre-processed data (Figure~\ref{figure:bgmax_ccPCR}~B) using the \textsl{CPP} 
showed comparable curvatures. We observed no significant difference between the raw 
and pre-processed data.

<<bgmax_ccPCR,fig.cap=fig19_cap,fig.scap=fig19_scap,fig.show='hold',message=FALSE,results='hide',warning=FALSE,fig.width = 11, fig.height = 8>>=
# Set parameter for the plot.
par(mfrow = c(2,1), las = 0, bty = "n")

# Use of bg.max for time-dependent measurements. Amplification curves 
# from the capillaryPCR dataset were processed in a loop. The results of 
#  bg.max are added to the plot. 

colors <- rainbow(8)

plot(NA, NA, xlim = c(0,75), ylim = c(-200,1300), xlab = "Time (min)", 
     ylab = "Voltage (micro V)", main = "ccPCR - Raw data")
mtext("A", cex = 1.5, side = 3, adj = 0)
for (i in c(1,3,5,7)) {
  x <- capillaryPCR[1L:750, i]
  y <- capillaryPCR[1:750, i + 1]
  res.bg <- summary(bg.max(x, y))
  lines(x, y, type = "b", pch = 20, col = colors[i], cex = 0.5)
  lines(c(res.bg[2], res.bg[2], res.bg[4], res.bg[4]), 
        c(-150, -50, -150, -50), col = colors[i], lwd = 1.5)
  text(10, 1200 - i * 50, 
       paste("bg.start: ", res.bg[1], ", bg.stop: ", res.bg[2], 
             ", amp.stop: ", res.bg[4]), col = colors[i], cex = 0.6)
}

plot(NA, NA, xlim = c(0,75), ylim = c(-200,1300), xlab = "Time (min)", 
     ylab = "Voltage (micro V)", main = "ccPCR - Pre-processed")
mtext("B", cex = 1.5, side = 3, adj = 0)
for (i in c(1,3,5,7)) {
  x <- capillaryPCR[1L:750, i]
  y <- CPP(capillaryPCR[1L:750, i], capillaryPCR[1:750, i + 1], 
           method = "mova", trans = TRUE, bg.range = c(1,105), 
           bg.outliers = TRUE)[["y.norm"]]
  res.bg <- summary(bg.max(x, y))
  lines(x, y, type = "b", pch = 20, col = colors[i], cex = 0.5)
  lines(c(res.bg[2], res.bg[2], res.bg[4], res.bg[4]), 
        c(-150, -50, -150, -50), col = colors[i], lwd = 1.5)
  text(10, 1200 - i * 50, 
       paste("bg.start: ", res.bg[1], ", bg.stop: ", res.bg[2], 
             ", amp.stop: ", res.bg[4]), col = colors[i], cex = 0.6)
}
@

\section{Normalization of amplification curve data}

As illustrated in Figure~\ref{figure:normalization}~\emph{A}, it is a common 
characteristic of data used for plotting amplification curves, that the fluorescence values in the baseline and plateau region 
vary between samples (e.g., due to high background, variances in dye 
quantities). Minimum and maximum values differ within an experiment. For 
visualization it is recommended to scale the data. This helps to grasp the data 
faster and alleviates the comparison of data from different measurements and/or 
scaling. \textsl{normalizer} is a function used to normalize any data by different 
methods (see details). To scale the fluorescence between 0 and 1 a \emph{Min-Max 
normalization} (Equation~\ref{eq:normalization}) is recommended 
\citep{roediger_RJ_2013}. We propose a quantile based normalization as 
alternative (Equation~\ref{eq:quantile_normalization}) since quantiles are less 
affected by outliers. The method can be invoked by the parameter $norm = 
"luqn"$. Although this does not scale all values between zero and one, we found 
it to be useful for noisy data. The parameter $qnL$ is symmetrically used to set 
the level for quantiles. By default, the 3~\% and 97~\% quantiles are used 
for the normalization. In addition, a maximum normalization
(Equation~\ref{eq:max_normalization}, 
Figure~\ref{figure:normalization}~\emph{D}) and a standard score normalization
(Equation~\ref{eq:zscore_normalization}, 
Figure~\ref{figure:normalization}~\emph{F}) is implemented.

\begin{equation} \label{eq:normalization}
RFU_{minmax} = \frac{RFU - \min(RFU)}{\max(RFU) - \min(RFU)}
\end{equation}

\begin{equation} \label{eq:max_normalization}
RFU_{max} = \frac{RFU}{\max(RFU)}
\end{equation}

\begin{equation} \label{eq:quantile_normalization}
RFU_{luqn} = \frac{RFU - Q_{p}(RFU)}{Q_{1 - p}(RFU) - 
Q_{p}(RFU)}
\end{equation}

\begin{equation} \label{eq:zscore_normalization}
RFU_{zscore} = \frac{RFU - \bar{x}_{RFU}}{s_{RFU}}
\end{equation}

The parameter $qnL$ is a user defined quantile, which is used for the 
quantile 
normalization. 
\begin{itemize}
\item A quantile normalization herein refers to an approach, which is less 
prone to 
outliers than a normalization to the minimum and the maximum of an 
amplification plot.
\item minm does a min-max normalization between 0 and 1 (see 
\citep{roediger_RJ_2013} for explanation).
\item max does a normalization to the maximum value (MFI/max(MFI)).
\item lugn does a quantile normalization based on a symmetric proportion 
as defined by the $qnl$ parameter (e.g., $qnl$ = 0.03 equals 3 and 97 
percent quantiles).
\item zscore performs a z-score normalization with a mean of 0 and a 
standard deviation of 1.
\end{itemize}

<<normalization,fig.show='hold',fig.cap=fig5_cap,fig.scap=fig5_scap,fig.width = 11, fig.height = 8>>=
par(mfrow = c(2,3), las = 0, bty = "n", oma = c(.5,.5,.5,.5))
tmp <- VIMCFX96_60

plot(NA, NA, xlim = c(1,40), ylim = c(0, 6000), xlab = "Cycle", 
     ylab = "RFU", main = "Raw data")
mtext("A", cex = 1.2, side = 3, adj = 0, font = 2) 
lin <- apply(tmp[, -1], 2, function(x) lines(tmp[, 1], x))
abline(lm(rowMeans(tmp[2:10, 2L:ncol(tmp)]) ~ tmp[2:10, 1]), col = 2)

plot(NA, NA, xlim = c(1,40), ylim = c(0, 3300), xlab = "Cycle", 
     ylab = "RFU", main = "Baselined data")
mtext("B", cex = 1.2, side = 3, adj = 0, font = 2) 
lin <- apply(tmp[, -1], 2, function(x) lines(tmp[, 1], CPP(tmp[, 1], x, 
                                                           method.norm = "none")$y))


plot(NA, NA, xlim = c(1,40), ylim = c(0, 1.15), xlab = "Cycle", 
     ylab = "RFU", main = "MinMax-Normalization")
mtext("C", cex = 1.2, side = 3, adj = 0, font = 2) 
lin <- apply(tmp[, -1], 2, function(x) lines(tmp[, 1], CPP(tmp[, 1], x, 
                                                           method.norm = "minm")$y))

plot(NA, NA, xlim = c(1,40), ylim = c(0, 1.15), xlab = "Cycle", 
     ylab = "RFU", main = "Max-Normalization")
mtext("D", cex = 1.2, side = 3, adj = 0, font = 2) 
lin <- apply(tmp[, -1], 2, function(x) lines(tmp[, 1], CPP(tmp[, 1], x,, 
                                                           method.norm = "max")$y))

plot(NA, NA, xlim = c(1,40), ylim = c(0, 1.15), xlab = "Cycle", 
     ylab = "RFU", main = "luqn-Normalization")
mtext("E", cex = 1.2, side = 3, adj = 0, font = 2) 
lin <- apply(tmp[, -1], 2, function(x) lines(tmp[, 1], CPP(tmp[, 1], x, 
                                                           method.norm = "luqn", qnL = 0.03)$y))

plot(NA, NA, xlim = c(1,40), ylim = c(-1.5, 1.5), xlab = "Cycle", 
     ylab = "RFU", main = "zscore-Normalization")
mtext("F", cex = 1.2, side = 3, adj = 0, font = 2) 
lin <- apply(tmp[, -1], 2, function(x) lines(tmp[, 1], CPP(tmp[, 1], x, 
                                                           method.norm = "zscore")$y))
@


\subsection{Computing linear model coefficients - Background subtraction based 
on linear models}

The slope of the background range is often unequal to zero and in most cases and is 
accompanied by a positive or negative trend. It is however possible to correct the slope by a linear trend extrapolation. The functions \textsl{lm.coefs} and 
\textsl{CPP} are wrappers used by functions to perform normal (linear least squares) 
and robust linear regression. \textsl{lm.coefs} calculates the estimated 
background of an amplification curve. This includes a ordinary least squares 
method (\textsl{lm}, \emph{stats}) and three other robust methods. Robust 
regression methods are less vulnerable to outliers. This feature is especially 
useful, when the background range contains noise. These robust methods are (I) a 
nonparametric rank-based estimator \citep{Kloke_2012}, (II) quantile regression 
\citep{Koenker_2008} and (III) a MM-type estimators for linear regression 
\citep{Todorov_2009}. By default, the MM-type estimator is used. Under the 
assumption that the background is constant, \textsl{CPP} or 
\textsl{lm.coefs} uses a defined range of the amplification curve (e.g., background 
range) to extrapolate a linear trend over the entire data. The coefficients of 
the analysis can be used for a trend-based correction of the entire dataset 
(Figure~\ref{figure:lmcoef}). If the robust linear regression is impossible, 
 \textsl{lm.coefs} performs a linear regression using the least squares method. 
However, this operation effects the AE. Caution should be exercised 
when using trans with time series (see \textsl{lm} from the $\emph{stats}$ 
package for details).

<<lmcoef,fig.show='hold',fig.cap=fig17_cap,fig.scap=fig17_scap,fig.width = 11, fig.height = 8>>=
par(bty = "n")
plot(VIMCFX96_69[, 1], VIMCFX96_69[, 2], type = "l", xlab = "Cycle", 
     ylab = "Fluorescence")
rect(1,0,10,5000)
method <- c("lmrob", "rq", "least", "rfit")
for (i in 1:4) {
  tmp <- lm.coefs(VIMCFX96_69[1:10, 1], VIMCFX96_69[1:10, 2], 
                  method.reg = method[i])
  text(9, 3200 - i * 100, paste(method[i], ":", "m: ", 
                                round(tmp[1,1], 4), "n: ", round(tmp[2,1], 3)))
  abline(a = tmp[1, 1], b = tmp[2, 1], col = i + 1, lwd = 1.5)
}
legend("right", c("Data", "lmrob", "rq", "least", "rfit"), lty = 1, 
       col = 1:5, cex = 0.95)
@

\section{The \textsl{inder} function - an interpolating five-point stencil}
\label{sec:inder}

Many methods for curve analysis require the calculation of derivatives. It is 
possible to solve this by fitting a curve to a function and performing symbolic 
derivation. Unfortunately, this approach causes information loss through the fitting process and unnecessarily adds additional assumptions regarding the relation between cycle 
number and fluorescence level. Hence, we integrated the \textsl{inder} function. 
\textsl{inder} (``in'' and ``der'' = interpolate derivatives) finds numeric 
derivatives by a five-point stencil, a commonly used finite difference method. 
These methods approximate derivative in a given point by adding up products of 
nearby values of function and their weights~\citep{Dahlquist_2008}. This function 
can estimate the approximate quantification cycles (Cq). 
Differentiation is a method for background suppression and reduction of the 
inter sample background amplitude variations 
(Figure~\ref{figure:inder_fit}~A~and~B). Smoothing may enhance the calculation 
of derivatives and optimize the signal-to-noise ratio. Therefore, we 
implemented spline interpolation and Friedman's SuperSmoother. 
However, the use of this smoother is limited in use to other functions such 
as \textsl{bg.max}. The parameter $Nip$ (default $Nip = 4$) is used to define 
how often an interpolation takes place at n equidistant points within the first 
and the last cycle. A high $Nip$ may improve the precision. However, a $Nip$ less 
than 2 and higher than 20 are not meaningful for conventional qPCR with 30 to 50 
cycles. In the context of qIA, a higher $Nip$ might be appropriate.


\subsection{Quantitative description of amplification reactions}

The Cq is a relative value, which depends on the template copy number, 
instrument, reagents, AE and probe technology. Low Cqs 
correlate with high quantities of template copy numbers. Real-time technologies 
enable the quantification of nucleic acids by calculation of specific curve 
parameters like the Cq and the AE based on the kinetics of the amplification curve. The Cq represents the 
number of cycles (time for qIA) needed to reach a defined fluorescence signal 
level in the exponential phase of the amplification curve. The Cq can be 
determined from a fixed threshold value or by various analytical algorithms as 
described elsewhere \citep{ruijter_2009, ruijter_2013, tellinghuisen_2014}. The 
output of \textsl{inder} includes the first derivative maximum ($FDM$) and 
second derivative maximum ($SDM$), which are commonly used in qPCR experiments 
but might be useful for isothermal amplification processes, too. 
Figure~\ref{figure:SDM} shows a typical result of the \textsl{inder} function. 
Following we show three examples that explain properties of \textsl{inder} and 
illustrate applications of the function in combination with other functions.

The \textsl{inder} function calculates numeric derivatives on 
smoothed data, which results in data points not observable in reality. The 
\textsl{rounder} function averages such result to the real values of cycle 
number

<<rounder_show,warning=FALSE,message=FALSE,fig.show='hold',fig.cap=fig16_cap,fig.scap=fig16_scap,fig.width = 11, fig.height = 8>>=
# Simulate an amplification curve with 40 cycles using the AmpSim 
# function.
isPCR <- AmpSim(cyc = 1:40)

# Use inder to calculate the derivatives and assign the results to the 
# object res
res <- inder(isPCR)

# Process res by rounder and assign the results to the object rd
rd <- rounder(res)

# Print details of res and rd. Due to the internal use of interpolating 
# splines in inder are the number of elements in the object res the n-th 
# time of the raw data. In this case 200 virtual instead of 40 real cycles.
head(res)
summary(res)

head(rd)
# summary(rd)
@

Figure~\ref{figure:SDM} illustrates the most important parameters of the 
\textsl{inder} function. We used the \textsl{AmpSim} function to simulate 
an ideal ``noise-free'' amplification curve with the default set to 
calculate the second derivative maximum ($SDM$) with \textsl{inder}. If 
$logy$ is $TRUE$ then a semi-decadic log scale graph (corresponds to the 
linear phase) to illustrate the exponential dynamic of the qPCR 
amplification is used. The parameter $logy$ is $FALSE$ by default. To the 
best of our knowledge, \textsl{inder} is the first tool in \textbf{R}, which 
allows a user to numerically derive his data without fitting them to any 
function or a combination of functions. The universality of this stencil approach 
can find an application even in problems not related to the analysis of 
amplification curves.

<<SDM,fig.show='hold',fig.cap=fig9_cap,fig.scap=fig9_scap,fig.width = 11,fig.height = 8>>=
# Use AmpSim to generate an amplification curve with 40 cycles
# and an approximate Cq of 20 and assign it to the object isPCR.
# isPCR is an object of the class "data.frame".
isPCR <- AmpSim(cyc = 1:40, Cq = 20)

# Invoke the inder function for the object isPCR to interpolate 
# the derivatives of the simulated data as object res. The Nip 
# parameter was set to 5. This leads to smoother curves. res is
# an object of the class "der".
res <- inder(isPCR, Nip = 5)

# Plot the object res and add descriptions to the elements.

par(las = 0, bty = "n", oma = c(.5,.5,.5,.5))

plot(isPCR, xlab = "Cycle", ylab = "RFU", ylim = c(-0.15,1),
     main = "", type = "b", pch = 20, lwd = 2)
colors <- rainbow(4)
# Add graphical elements for the derivatives and the calculated
# Cq values FDM, SDM, SDm and SDC.

lines(res[, "x"], res[, "d1y"], col = "blue", lwd = 2)
lines(res[, "x"], res[, "d2y"], col = "red", lwd = 2)

# Fetch the Cq values from res with the summary function
summ <- summary(res, print = FALSE)

abline(v = summ, col = colors, lwd = 2)
text(15, 0.3, paste("FDM ~ ", round(summ["FDM"], 2)), 
     cex = 1.1, col = colors[1])
text(15, 0.2, paste("SDM ~ ", round(summ["SDM"], 2)), 
     cex = 1.1, col = colors[2])
text(15, - 0.1, paste("SDm ~ ", round(summ["SDm"], 2)), 
     cex = 1.1, col = colors[3])
text(15, 0.7, paste("SDC ~ ", round(summ["SDC"], 2)), 
     cex = 1.1, col = colors[4])

legend(1.1, 0.9, c("Raw data", "First derivative", "Second derivative"), 
       col = c(1,4,2), lty = c(2,1,1), bty = "n")

# Summary of the object res.
summ
@

\textsl{inder} is a helper function, which can be part of other routines. 
Recently, we added this approach to the \textsl{diffQ} function of the 
\emph{MBmca} for improved predictions. The \textsl{diffQ} function is part of a 
routine to calculate the melting points of nucleic acids 
\citep{roediger_RJ_2013}. The $FDM$ and $SDM$ are peak values to determine the 
Cq. We used the \textsl{inder} function in \textsl{diffQ} to compare the Cq 
values between a quantification experiment where the samples were either 
detected with a gene specific hydrolysis probe or the intercalating dye 
EvaGreen. For the analysis we focused on the $SDM$. We found that the samples 
detected with EvaGreen had a slightly lower Cq (Figure~\ref{figure:inder}~A) 
than samples detected with the hydrolysis probe (Figure~\ref{figure:inder}~B). 
The mean variation of the Cq was less in samples where EvaGreen was used for monitoring.

<<inder,fig.show='hold',fig.cap=fig10_cap,fig.scap=fig10_scap,message=FALSE,results='hide'>>=
# Plot all data from C127EGHP and calculate the SDM (Second Derivative 
# Maximum) values with the diffQ2() function (Note: the inder parameter
# is set as TRUE)
# first plot the samples detected with EvaGreen and next the samples 
# detected with the Hydrolysis probe
require(MBmca)

pointer <- function (x, pos = 1, w = 5, stat = TRUE){
  xx <- pos + rep(seq(-0.1, 0.1, length.out = w), ceiling(length(x)/w))
  yy <- sort(x)
  points(xx[1:length(yy)], yy, pch = 19)
  
  if (stat == TRUE)
    x.median <- median(x, na.rm = T)
  x.mad <- mad(x, na.rm = T) * 2
  param <- c(length= 0, code = 3, pch = 15, cex = 1.2)
  arrows(xx[1] * 0.98, x.median, tail(xx, 1) * 1.02, 
         x.median, param, lwd = 3, col = 2)
  arrows(xx[1] * 1.01, x.median + x.mad, tail(xx, 1) * 0.99, 
         x.median + x.mad, param, lwd = 2, lty = 2, col = 4)
  arrows(xx[1] * 1.01, x.median - x.mad, tail(xx, 1) * 0.99, 
         x.median - x.mad, param, lwd = 2, lty = 2, col = 4)
}

amp.liner <- function(range, input, colors = "black") {
  sapply(range, function(i) {
    lines(input[, 2], input[, i], col = colors, pch = 19)
    tmpP <- mcaSmoother(input[, 2], input[, i])
    SDM <- diffQ2(tmpP, inder = TRUE)[["xTm1.2.D2"]][1]
    abline(v = SDM)
    SDM
  }
  )
}

layout(matrix(c(1,3,2,3), 2, 2, byrow = TRUE), respect = TRUE)
par(las = 0, bty = "n")
plot(NA, NA, xlim = c(1,40), ylim = c(0,10), xlab = "Cycle", 
     ylab = "Fluorescence", main = "EvaGreen")
mtext("A", cex = 1.1, side = 3, adj = 0, font = 2)

EG <- amp.liner(range = 3L:34, input = C127EGHP)

plot(NA, NA, xlim = c(1,40), ylim = c(0,10), xlab = "Cycle", 
     ylab = "Fluorescence", main = "Hydrolysis probe")
mtext("B", cex = 1.1, side = 3, adj = 0, font = 2)

HP <- amp.liner(range = 35L:66, input = C127EGHP)

plot(NA, NA, xlim = c(0.8,2.2), ylim = c(13,14), xaxt = "n", 
     xlab = "", ylab = "Cq (SDM, diffQ2)")
text(c(1.05,2), c(13.05,13.05), c("EG", "HP"), cex = 1.2)
mtext("C", cex = 1.1, side = 3, adj = 0, font = 2)
pointer(EG, pos = 1, w = 8)
pointer(HP, pos = 2, w = 8)
@

\section{Quantification cycle calculation by the \textsl{inder} function}

\subsection{The \textsl{inder} function in combination with a 
5-parameter curve fit function}

In the previous example we used smoothing and the \textsl{inder} method to 
calculate the $SDM$. Smoothing alters the peak signal (e.g., peak height 
reduction and peak width increase are commonly encountered problems). An alternative 
technique to determine the $FDM$ or $SDM$ is by fitting the raw data. In the 
next example we used the \textsl{drm} function from the \emph{drc} package 
\citep{Ritz_2005} to fit a five-parameter log-logistic function (S-shaped). The 
\textsl{inder} function was used to calculate the $SDM$ of the predicted models 
(Figure~\ref{figure:inder_fit}).

<<inder_fit,fig.cap=fig11_cap,fig.scap=fig11_scap,fig.show='hold',message=FALSE,fig.width = 11,fig.height = 8,results='hide'>>=
fit.amp <- function(cyc, fluo, plot = FALSE) {
  
  ampl <- quantile(fluo, 0.999)
  bl <- quantile(fluo, 0.001)
  Cq <- round(mean(cyc))
  b.eff <- 1
  
  fit <- nls(fluo ~ bl + ampl / (1 + exp(- (cyc - Cq) / b.eff)), 
             start = list(Cq = Cq, b.eff = b.eff, ampl = ampl, 
                          bl = bl)
  )
  
  res.pred <- data.frame(cyc, predict(fit))
  res <- inder(res.pred[, 1], res.pred[, 2])
  if (plot) {
    lines(res[, 1], res[, 4])
  }
  # SDM
  summary(res)[2]
}

tmp <- C126EG595

out <- apply(tmp[, -1], 2, function(x) fit.amp(tmp[, 1], x))

layout(matrix(c(1,2,1,3), 2, 2, byrow = TRUE))

plot(NA, NA, xlim = c(1,40), ylim = c(min(tmp[, 2L:97]), 
                                      max(tmp[, 2L:97])), 
                                      xlab = "Cycle", 
                                      ylab = "Raw fluorescence")
mtext("A", cex = 1.2, side = 3, adj = 0, font = 2)
for (i in 2L:97) {
  lines(tmp[, 1], tmp[, i], col = ifelse(out[i - 1] < 15.5, "red", 
                                         "black"), lwd = 2)
}
abline(v = out)

plot(NA, NA, xlab = "Cycle", ylab = "RFU''(Cycle)", main = "", 
     xlim = c(0,40), ylim = c(-850, 850))
abline(v = 15.5, lty = 2)
invisible(apply(tmp[, -1], 2, function(x) {
  fit.amp(tmp[, 1], x, plot = TRUE)
}
))
mtext("B", cex = 1.2, side = 3, adj = 0, font = 2)

hist(out, xlab = "Cq (SDM)", main = "", 
     breaks = seq(14.8, 15.8, 0.05), col = rainbow(96))
abline(v = 15.5, lty = 2)
mtext("C", cex = 1.2, side = 3, adj = 0, font = 2)
@

\section{Threshold cycle method}

The \textsl{th}.\textsl{cyc} function calculates the Cq in qPCRs or cycle time in 
quantitative isothermal amplification (qIA) experiments. This method requires proper 
baselining. Note: more sophisticated quantification methods exist in qPCR analysis   
\citep{ruijter_2009, ruijter_2013, tellinghuisen_2014}. This function was 
implemented primarily for the analysis qIA experiments but is usable 
for qPCR too. We implemented a symmetric approximation algorithm 
based on linear and quadratic least squares regression. The Threshold Cycle (Ct) 
is the cycle number at which the reporter fluorescence signal significantly exceeds  a point above the baseline and defined threshold in particular samples. The 
\textsl{th}.\textsl{cyc} calculates the intersection of the user defined Ct 
value (r) and a linear regression or quadratic polynomial in the range of the 
user defined Ct value. Therefore, \textsl{th}.\textsl{cyc} has no requirement to 
fit a ``complex'' non linear model to the entire dataset but rather focuses on 
a specific area. The polynomial is calculated from four neighbor values of the 
fluorescence threshold.

<<thcyc,warning=FALSE,message=FALSE,fig.show='hold',fig.cap=fig2_cap,fig.scap=fig2_scap,fig.width = 11, fig.height = 8>>=
# Raw data from the VIMCFX96_69 dataset.
# Cycles x and Fluoresce values y
x <- VIMCFX96_69[, 1]
y <- VIMCFX96_69[, 2]

par(mfrow = c(2,1), las = 0, bty = "n")

# Plot the raw data
plot(x, y, xlab = "Cycle", ylab = "Fluo", main = "Linear regression", 
     pch = 19)
mtext("A", cex = 1.3, side = 3, adj = 0) 
# Calculate the Cq (Ct) value
res <- th.cyc(x, y, r = 2400, linear = TRUE)
lines(res@input, col = 2, lwd = 2)

# Threshold fluorescence value
abline(h = res[2], col = 3)

# Calculated Ct value
abline(v = res[1], col = 4)
legend("topleft", paste("Cq (Ct) = ", round(res[1], 3)))

plot(x, y, xlab = "Cycle", ylab = "Fluo", main = "Quadratic regression", 
     pch = 19)
mtext("B", cex = 1.3, side = 3, adj = 0) 

# Calculate the Ct value
res <- th.cyc(x, y, r = 2400, linear = FALSE)
lines(res@input, col = 2, lwd = 2)

# Threshold fluorescence value
abline(h = res[2], col = 3)

# Calculated Ct value
abline(v = res[1], col = 4)
legend("topleft", paste("Cq (Ct) = ", round(res[1], 3)))
@


\subsection{Application of the \textsl{th.cyc} function on ccPCR data}
\label{sec:ccPCR}

In the following example we analyzed a continuous amplification reaction. The 
data were taken from the $\emph{capillaryPCR}$ dataset from the 
$\emph{chipPCR}$ package measured with capillary convective PCR (ccPCR) 
technology. We used a modified device of the ccPCR system as described by 
\citep{chou_rapid_2011}. For technical details see the supplement of 
\citep{spiess_impact_2014}. The PCR was performed with SYBR(R) Green fluorescent 
intercalating dye. Thereof the ESE-Log has an excitation of 470 nm and a detection wavelenth of 520 nm in one channel. The data was recorded by the FL Digital Software 
(QIAGEN Lake Constance) and exported as text-based raw data. The raw data were 
noisy and showed an off-set of circa 150 micro Volt with a slightly negative trend 
(Figure~\ref{figure:bgmax_ccPCR}~A). The data were first pre-processed 
(baselined and slightly smoothed) by the \textsl{CPP} function. The \textsl{th.cyc} function 
was used to determine the time required to reach a certain threshold level above 
the defined value (Threshold level ``r'' (80 $\mu$Volts)) 
(Figure~\ref{figure:thcyc_ccPCR}~B).

<<thcyc_ccPCR,fig.show='hold',fig.cap=fig3_cap,fig.scap=fig3_scap,fig.width = 11, fig.height = 8>>=
# Application of the th.cyc method to determine the Cq from a continuous
# amplification reaction.
par(las = 0, bty = "n", oma = c(.5,.5,.5,.5))
       
plot(NA, NA, xlim = c(0,80), ylim = c(0,1200), xlab = "Time (min)", 
     ylab = "Voltage [micro V]", main = "ccPCR - Pre-processed Data")
mtext("B", cex = 2, side = 3, adj = 0)
# Threshold level "r" (50 micro Volts)
for (i in c(1,3,5,7)) {
  y.tmp <- CPP(capillaryPCR[, i], capillaryPCR[, i + 1], trans = TRUE, bg.range = c(1,150))$y.norm
  Ct.tmp <- th.cyc(capillaryPCR[, i], y.tmp, r = 80, linear = FALSE)
  abline(v = Ct.tmp[1])
  text(Ct.tmp[1] * 1.125, 1200, paste(round(Ct.tmp[1], 1), "\nmin"), cex = 0.8)
  lines(capillaryPCR[, i], y.tmp, type = "b", pch = 20 - i) 
  points(Ct.tmp@input, col = "red", pch = 19)
}
abline(h = 80)
legend("topleft", c("Run 1", "Run 2", "Run 3", "Control"), 
       pch = c(19, 17, 15, 13), lwd = 1.3, bty = "n")
@


\subsection{Application of the \textsl{th.cyc} and \textsl{CPP} functions for helicase 
dependent Amplification}
\label{sec:HDA}

In this example for time-dependent reactions, we examined a helicase-dependent 
amplification (HDA). The target was human \textit{vimentin} (\textit{vim}). The VideoScan Platform 
\citep{roediger_highly_2013} was used to monitor the amplification. The HDA was 
performed at $65\,^{\circ}\mathrm{C}$ \citep{spiess_impact_2014}. Three 
concentrations of input DNA (D1, D2, D3) were used (Figure~\ref{figure:HDA}). In 
detail, the IsoAmp(R) III Universal tHDA Kit (Biohelix) was used. Primers and 
templates are described in \citep{roediger_highly_2013}. Reaction mix A was: 10 $\mu$L A. bidest, 1.25 $\mu$L 10xbuffer, 0.75 $\mu$L primer(150~nM final), 0.5 $\mu$L template plasmid. Preincubation: This mixture was incubated for 2 min at $95\,^{\circ}\mathrm{C}$ and immediately placed on ice. Reaction mix B: 5 $\mu$L A. bidest., 1.25 $\mu$L 10x buffer, 2 $\mu$L NaCl, 1.25 $\mu$L MgSO$_{4}$, 1.75 $\mu$L dNTPs, 0.25 $\mu$L EvaGreen, 1 $\mu$L enzyme mix. The mix was covered with 50 $\mu$L mineral oil. Measurement started immediately after after adding buffer B. Three dilutions of input DNA were used (1x (D1), a 1:10  (D2) and a 1:100 (D3) dilution). 
Temperature profile (after Preincubation):

\begin{itemize}
 \item 60 seconds at $65\,^{\circ}\mathrm{C}$,
 \item 11 seconds at $55\,^{\circ}\mathrm{C}$ \& Measurement
\end{itemize}


<<HDA,fig.show='hold',fig.cap=fig13_cap,fig.scap=fig13_scap,message=FALSE,warning=FALSE,fig.width = 11, fig.height = 8>>=
par(mfrow = c(2,1), bty = "n")
plot(NA, NA, xlim = c(0,5000), ylim = c(0.45,0.8), xlab = "Time (sec)", 
     ylab = "Fluorescence", main = "HDA - Raw data")
mtext("A", cex = 2, side = 3, adj = 0)
lines(C85[, 2], C85[, 3], type = "b", col = 2, pch = 20)
lines(C85[, 4], C85[, 5], type = "b", col = 4, pch = 20)
lines(C85[, 6], C85[, 7], type = "b", col = 6, pch = 20)
legend("bottomright", c("D1, 1x", "D2, 1:10", "D3, 1:100"), col = c(2,4,6), 
       pch = rep(20,3), bty = "n")

plot(NA, NA, xlim = c(0,2000), ylim = c(0,0.4), xlab = "Time (sec)", 
     ylab = "Fluorescence", main = "HDA - Pre-processed data")
mtext("B", cex = 2, side = 3, adj = 0)
legend("topleft", c("D1, 1x", "D2, 1:10", "D3, 1:100"), col = c(2,4,6), 
       pch = rep(20,3), bty = "n")

# Define the parameters for the pre-processing by CPP and the th.cyc 
# function.
# smoothing method
sm <- "mova"

# manual range for background
br <- c(2,10)

# time range for analysis
xr <- 3L:200

# method for baseline normalization
lrg <- "least"

# threshold level for the th.cyc function
r <- 0.025
# Calculate in a loop the Cq values (Cycle threshold method) and add the
# calculated time (in minutes) to the plot.
for (i in c(2,4,6)) {
  y.tmp <- CPP(C85[xr, i], C85[xr, i + 1], method = sm, bg.range = br, 
               trans = TRUE)$y.norm
  Ct.tmp <- th.cyc(C85[xr, i], y.tmp, r = r, linear = FALSE)
  abline(v = Ct.tmp[1], col = "grey")
  lines(C85[xr, i], y.tmp, col = i, lwd = 2)
  points(Ct.tmp@input, col = "red", pch = 19)
  text(Ct.tmp[1] * 1.1, 0.36, paste(round(Ct.tmp[1]/60, 1), "\nmin"))
}

# Show the fluorescence value, which defines the threshold.
abline(h = r, lty = 2)
@


\section{Amplification efficiency}

Various factors influence amplification reactions. The amplification efficiency 
(AE) is controlled by a complex interaction of the intrinsic and extrinsic factors 
like reaction conditions, substrate consumption, primer dimmer formation and 
molecule specific reaction rates \citep{mehra_2005}. Some probe systems are 
perceived as bias-introducing. Therefore, qPCR reactions should be corrected 
based on the AE \citep{tuomi_2010, ruijter_2014}. The 
AE can be estimated from individual samples or a set 
of samples to compensate the presence of inhibitors and noise. Indirect methods 
use fitted mathematical models or estimate the AE from absolute fluorescence 
values \citep{liu_2002, tichopad_2003, alvarez_2007, smith_2007, batsch_2008, 
mallona_2011}. The \emph{qpcR} function has many functions for indirect estimation of the AE. However, the most commonly used is the ``direct method'' \citep{liu_2002, 
stahlberg_2003}, whereby the AE is estimated from dilution series of a template. 
The AE of a qPCR reaction is calculated from the slope of the standard curve 
(Equation~\ref{eq:AE}).

\begin{equation} \label{eq:AE}
AE = \frac{10^{(-1/m)}}{2} * 100
\end{equation}

\textsl{effcalc} is used for the automatic calculation of the AE of 
a dilution series (Figure~\ref{figure:AmpSim_effcalc}). An object of the class 
list contains the ``Concentration'', Cqs, deviation of the Cqs, "Coefficient of 
Variance" sequentially in the columns, the AE (\%) 
according to Equation~\ref{eq:AE}, the results of the linear regression and the 
correlation test (Pearson's) (Table~\ref{table:effcalc_output}). The 
\textsl{effcalc} has several options to enhance the plot. For example, it is 
possible to indicate the confidence interval (default CI~=~95~\%). Further 
options are described in the chipPCR manual.

Providing that the smoother is properly adjusted, it is possible to detect 
only the significant peaks while small or to narrow peaks are ignored. 
\textsl{smoother} is used by other functions of $\emph{chipPCR}$ like 
\textsl{CPP}. The example for Figure~\ref{figure:inder} illustrates the 
use of the \textsl{diffQ} and \textsl{diffQ2} functions from the 
\emph{MBmca} package \citep{roediger_RJ_2013} and the integration of the \textsl{inder} function. The 
\textsl{inder} function is used in \textsl{diffQ} and \textsl{diffQ2} for a 
precise peak location, while the approximate $SDM$ is calculated from the 
derivative of a quadratic function at the approximate $SDM$. 

<<AmpSim_effcalc,fig.show='hold',fig.cap=fig6_cap,fig.scap=fig6_scap,fig.width = 11,fig.height = 8,message=FALSE>>=

# Load MBmca package (v. 0.0.3-3 or later)
require(MBmca)

# Create an graphic device for two empty plots.
par(mfrow = c(1,2))
plot(NA, NA, xlim = c(1,45), ylim = c(0.01,1.1), xlab = "Cycles", 
     ylab = "Fluorescence", main = "")
mtext("A", cex = 1.1, side = 3, adj = 0, font = 2)

# Create a sequence of "targeted" Cq values (Cq.t) between 15 and 34 
# cycles.

Cq.t <- rep(seq(15, 34, 3.5), 3)

# In-silico experiment set up: Define the levels for the decadic dilutions
# with concentrations from 100 to 0.001 (six steps) as three replicates.

dilution <- rep(10^(2:-4), 3)

# Create an empty matrix for the results of the concentration
# dependent Cq values.

ma.out <- matrix(data = NA, nrow = 45, ncol = length(Cq.t))

# Use AmpSim to simulate amplification curves at different concentrations. 
# The simulation is performed with the addition of some noise. This 
# generates unique (non-reproducible) amplification curves, even under 
# identical parameter settings.

Cq.out <- vector()

# Simulate a qPCR reaction with AmpSim for 45 cycles and some noise.

for (i in 1L:18) {
  ma.out[1:45, i] <- AmpSim(cyc = c(1:45), b.eff = -50, bl = 0.001, 
                            ampl = 1, Cq = Cq.t[i], noise = TRUE, 
                            nnl = 0.02)[, 2]
  lines(1:45, ma.out[, i])
  tmpP <- mcaSmoother(1:45, ma.out[, i])
  # Calculate the pseudo Second Derivative Maximum (SDM) (Cq) using 
  # the diffQ2 function from the MBmca package.
  Cq.tmp <- diffQ2(tmpP, inder = TRUE)$xTm1.2.D2[1]
  abline(v = Cq.tmp)
  Cq.out <- c(Cq.out, Cq.tmp)
}

# Assign the calculated Cqs to the corresponding concentrations.
tmp <- data.frame(dilution[1:6], Cq.out[1:6], Cq.out[7:12], Cq.out[13:18])

# Determine the amplification efficiency by using the effcalc function.
plot(effcalc(tmp[, 1], tmp[, 2:4]), CI = TRUE)
mtext("B", cex = 1.1, side = 3, adj = 0, font = 2) 
@

Then we analyzed the $\emph{C54}$ dataset from the $\emph{chipPCR}$ package with the \textsl{effcalc} 
function. Herein, a qPCR experiment for the amplification of 
\textit{MLC-2v} was performed using the VideoScan heating/cooling-unit. Calculation of the Cq required pre-processing of the data. One amplification curve contained a missing value 
(Figure~\ref{figure:CPP_C54}~A) which was removed by the spline method of \textsl{CPP}. 
Additionally, the data were baselined (linear model, robust MM-estimator) and smoothed by 
Savitzky-Golay smoothing (Figure~\ref{figure:CPP_C54}~B). A final analysis 
with the \textsl{effcalc} function showed that the AE was 
circa 87.3~\% for the gene \textit{MLC-2v} using the VideoScan HCU 
(Figure~\ref{figure:CPP_C54}~C). However, since only few measurement points were 
tested for this experiment, it just safe to say that the hardware of the HCU works 
reliably under the experimental conditions.

<<CPP_C54,fig.show='hold',fig.cap=fig7_cap,fig.scap=fig7_scap,fig.width = 11,fig.height = 8,message=FALSE>>=
require(MBmca)
par(las = 0, bty = "n", oma = c(.5,.5,.5,.5))
par(fig = c(0,0.5,0,1), new = TRUE)
plot(NA, NA, xlim = c(1,55), ylim = c(0, 0.7), xlab = "Cycle", 
     ylab = "refMFI", main = "Raw data")
just_line <- apply(C54[, c(2:4)], 2, function(y) lines(C54[, 1], y))
mtext("A", cex = 1.2, side = 3, adj = 0, font = 2) 

par(fig = c(0.5,1,0.5,1), new = TRUE)
plot(NA, NA, xlim = c(1,55), ylim = c(0, 0.55), xlab = "Cycle", 
     ylab = "refMFI", main = "pre-processed data")
mtext("B", cex = 1.2, side = 3, adj = 0, font = 2) 

D1 <- cbind(C54[1:35, 1], CPP(C54[1:35, 1], C54[1:35, 2], trans = TRUE, 
                              bg.range = c(1,8))[["y.norm"]])
D2 <- cbind(C54[1:45, 1], CPP(C54[1:45, 1], C54[1:45, 3], trans = 
                                TRUE)[["y.norm"]])
D3 <- cbind(C54[1:55, 1], CPP(C54[1:55, 1], C54[1:55, 4], trans = 
                                TRUE)[["y.norm"]])

lines(D1, col = 1)
lines(D2, col = 2)
lines(D3, col = 3)

dilution <- c(1E0, 1E-3, 1E-6)
Cq.D1 <- diffQ2(D1, inder = TRUE)[["xTm1.2.D2"]][1]
Cq.D2 <- diffQ2(D2, inder = TRUE)[["xTm1.2.D2"]][1]
Cq.D3 <- diffQ2(D3, inder = TRUE)[["xTm1.2.D2"]][1]

res.dil <- data.frame(dilution, rbind(Cq.D1, Cq.D2, Cq.D3))
par(fig = c(0.5,1,0,0.5), new = TRUE)
plot(effcalc(res.dil[, 1], res.dil[, 2]), res.fit = NULL)
@


<<effcalc_output,echo=FALSE,results='asis',message=FALSE>>=
print(xtable(effcalc(res.dil[, 1], res.dil[, 2]), 
      caption = "Output of the effcalc function.", label = "table:effcalc_output"))
@

In another example we used \textsl{effcalc} function to analyze the $\emph{C60.amp}$ data 
set from the $\emph{chipPCR}$ package. All data of the human genes \textit{vimentin} 
(Figure~\ref{figure:effcalc_VIM_MLC}~A) and \textit{MLC-2v} 
(Figure~\ref{figure:effcalc_VIM_MLC}~B) were amplified in a Roche Light Cycler 
1.5 and detected by the HRM dye EvaGreen in independent experiments. As shown in 
the code and Figure~\ref{figure:effcalc_VIM_MLC} it is possible to obtain a 
complete analysis with few commands. The amplification efficiencies for both 
qPCRs was higher than 94~\% (Figure~\ref{figure:effcalc_VIM_MLC}~C~and~D, 
Table~\ref{table:effcalc_output}).

<<effcalc_VIM_MLC,fig.show='hold',fig.cap=fig25_cap,fig.scap=fig25_scap,fig.width = 11, fig.height = 8>>=
colors <- rep(rainbow(7), each = 2)
par(mfrow = c(2,2))

plot(NA, NA, xlim = c(0,44), ylim = c(0, 6), xlab = "Cycles", ylab = "RFU")
legend(0, 6, colnames(C60.amp[, 4L:17]), ncol = 2, col = colors[1:14], 
       pch = 19, bty = "n")
mtext("A", cex = 1.2, side = 3, adj = 0, font = 2)
SDM.vim <- sapply(4L:17, function(i) {
  lines(C60.amp[, 1], C60.amp[, i], col = colors[i - 3])
  SDM <- summary(inder(C60.amp[, 1], C60.amp[, i]), print = FALSE)[2]
}
)

plot(NA, NA, xlim = c(0,44), ylim = c(0, 4), xlab = "Cycles", ylab = "RFU")
legend(0, 4, colnames(C60.amp[, 18L:31]), ncol = 2, col = colors[1:14], 
       pch = 19, bty = "n")
mtext("B", cex = 1.2, side = 3, adj = 0, font = 2)
SDM.mlc2v <- sapply(18L:31, function(i) {
  lines(C60.amp[, 1], C60.amp[, i], col = colors[i - 17])
  SDM <- summary(inder(C60.amp[, 1], C60.amp[, i]), print = FALSE)[2]
}
)

#create vector of dilutions
dil <- sort(rep(10^(0L:-6), 2), TRUE)

res <- cbind(dil, SDM.vim, SDM.mlc2v)

plot(effcalc(res[, 1], res[, 2]))
mtext("C", cex = 1.2, side = 3, adj = 0, font = 2)

plot(effcalc(res[, 1], res[, 3]))
mtext("D", cex = 1.2, side = 3, adj = 0, font = 2)
@

\section{Datasets}

\begin{enumerate}
<<summ.datasets,echo=FALSE,results='asis'>>=
load("datdf.RData")

datdf <- cbind(datdf, t(sapply(datdf[["dat.names"]], function(i) 
  dim(get(substr(i, 1, nchar(i) - 1))))))
tmp <- apply(datdf, 1, function(i) {
  safe.name <- i[1]
  if(grepl("_", safe.name, fixed = TRUE)) {
    safe.name <- sub("_", ".", i, fixed = TRUE)
  }
  cat("\\item Dataset: ", safe.name, "\n\\begin{itemize}\n\\item Dataset type: ", 
      i[3], "\n\\item Description: ", i[2], "\n\\item Number of variables: ", i[4], 
      "\n\\item Number of measurements: ", i[5], "\n\\end{itemize}\n")
})
@
\end{enumerate}

\section{Acknowledgment}

Part of this work was funded by the BMBF InnoProfileTransfer-Projekt 03 IPT 611X. 
Grateful thanks belong to all authors of the cited \textbf{R} packages, the \textbf{R} 
community and \textbf{RKWard} developers. We would like to thank 
Claudia Deutschmann (Brandenburg University of Technology Cottbus - Senftenberg, 
Germany), Ralf Himmelreich (Fraunhofer ICT-IMM, Germany) and Katharina Klat 
(Darmstadt University of Applied Sciences, Germany) for the contribution of 
data.

\textbf{Authors' contributions:} SR conceived the study, and participated in its design and coordination and 
wrote the manuscript. SR and MB jointly developed the software. All authors read 
and approved the final version of the manuscript.

\listoffigures

\bibliography{roediger-burdukiewicz}

\end{document}