Package 'qpcR'

Akaike's second-order corrected Information Criterion

Description

Calculates the second-order corrected Akaike Information Criterion for objects of class pcrfit, nls, lm, glm or any other models from which coefficients and residuals can be extracted. This is a modified version of the original AIC which compensates for bias with small $n$. As qPCR data usually has $\frac{n}{k} < 40$ (see original reference), AICc was implemented to correct for this.

Usage

AICc(object)

Arguments

object

a fitted model.

Details

Extends the AIC such that

$AICc = AIC+\frac{2k(k + 1)}{n - k - 1}$

with $k$ = number of parameters, and $n$ = number of observations. For large $n$, AICc converges to AIC.

Value

The second-order corrected AIC value.

Author(s)

Andrej-Nikolai Spiess

References

Akaike Information Criterion Statistics.
Sakamoto Y, Ishiguro M and Kitagawa G.
D. Reidel Publishing Company (1986).

Regression and Time Series Model Selection in Small Samples.
Hurvich CM & Tsai CL.
Biometrika (1989), 76: 297-307.

See Also

AIC, logLik.

Examples

m1 <- pcrfit(reps, 1, 2, l5)
AICc(m1)

---

Calculation of Akaike weights/relative likelihoods/delta-AICs

Description

Calculates Akaike weights from a vector of AIC values.

Usage

akaike.weights(x)

Arguments

x	a vector containing the AIC values.

Details

Although Akaike's Information Criterion is recognized as a major measure for selecting models, it has one major drawback: The AIC values lack intuitivity despite higher values meaning less goodness-of-fit. For this purpose, Akaike weights come to hand for calculating the weights in a regime of several models. Additional measures can be derived, such as $\Delta(AIC)$ and relative likelihoods that demonstrate the probability of one model being in favor over the other. This is done by using the following formulas:

delta AICs:

$\Delta_i(AIC) = AIC_i - min(AIC)$

relative likelihood:

$L \propto exp\left(-\frac{1}{2}\Delta_i(AIC)\right)$

Akaike weights:

$w_i(AIC) = \frac{exp\left(-\frac{1}{2}\Delta_i(AIC)\right)}{\sum_{k=1}^K exp\left(-\frac{1}{2}\Delta_k(AIC)\right)}$

Value

A list containing the following items:

deltaAIC	the $\Delta(AIC)$ values.
rel.LL	the relative likelihoods.
weights	the Akaike weights.

Author(s)

Andrej-Nikolai Spiess

References

Classical literature:
Akaike Information Criterion Statistics.
Sakamoto Y, Ishiguro M and Kitagawa G.
D. Reidel Publishing Company (1986).

Model selection and inference: a practical information-theoretic approach.
Burnham KP & Anderson DR.
Springer Verlag, New York, USA (2002).

A good summary:
AIC model selection using Akaike weights. Wagenmakers EJ & Farrell S.
Psychonomic Bull Review (2004), 11: 192-196.

See Also

AIC, logLik.

Examples

## Apply a list of different sigmoidal models to data
## and analyze GOF statistics with Akaike weights
## on 8 different sigmoidal models.
modList <- list(l7, l6, l5, l4, b7, b6, b5, b4)
aics <- sapply(modList, function(x) AIC(pcrfit(reps, 1, 2, x))) 
akaike.weights(aics)$weights

---

Calculation of qPCR efficiency using dilution curves and replicate bootstrapping

Description

This function calculates the PCR efficiency from a classical qPCR dilution experiment. The threshold cycles are plotted against the logarithmized concentration (or dilution) values, a linear regression line is fit and the efficiency calculated by $E = 10^{\frac{-1}{slope}}$. A graph is displayed with the raw values plotted with the threshold cycle and the linear regression curve. The threshold cycles are calculated either by some arbitrary fluorescence value (i.e. as given by the qPCR software) or calculated from the second derivative maximum of the dilution curves. If values to be predicted are given, they are calculated from the curve and also displayed within. calib2 uses a bootstrap approach if replicates for the dilutions are supplied. See 'Details'.

Usage

calib(refcurve, predcurve = NULL, thresh = "cpD2", dil = NULL, 
       group = NULL, plot = TRUE, conf = 0.95, B = 200)

Arguments

refcurve	a 'modlist' containing the curves for calibration.
predcurve	an (optional) 'modlist' containing the curves for prediction.
thresh	the fluorescence value from which the threshold cycles are defined. Either "cpD2" or a numeric value.
dil	a vector with the concentration (or dilution) values corresponding to the calibration curves.
group	a factor defining the group membership for the replicates. See 'Examples'.
plot	logical. Should the fitting (bootstrapping) be displayed? If FALSE, only values are returned.
conf	the confidence interval. Defaults to 95%, can be omitted with NULL.
B	the number of bootstraps.

Details

calib2 calculates confidence intervals for efficiency, AICc, adjusted $R^2_{adj}$ and the prediction curve concentrations. If single replicates per dilution are supplied by the user, confidence intervals for the prediction curves are calculated based on asymptotic normality. If multiple replicates are supplied, the regression curves are calculated by randomly sampling one of the replicates from each dilution group. The confidence intervals are then calculated from the bootstraped results.

Value

A list with the following components:

eff	the efficiency.
AICc	the second-order corrected AIC.
Rsq.ad	the adjusted $R^2_{adj}$.
predconc	the (log) concentration of the predicted curves.
conf.boot	a list containing the confidence intervals for the efficiency, the AICc, Rsq.ad and the predicted concentrations.

A plot is also supplied for efficiency, AICc, Rsq.ad and predicted concentrations including confidence intervals in red.

Author(s)

Andrej-Nikolai Spiess

Examples

## Define calibration curves,
## dilutions (or copy numbers) 
## and curves to be predicted.
## Do background subtraction using
## average of first 8 cycles. No replicates.
CAL <- modlist(reps, fluo = c(2, 6, 10, 14, 18, 22), 
               baseline = "mean", basecyc = 1:8)
COPIES <- c(100000, 10000, 1000, 100, 10, 1)
PRED <- modlist(reps, fluo = c(3, 7, 11), 
                baseline = "mean", basecyc = 1:8)

## Conduct normal quantification using
## the second derivative maximum of first curve.
calib(refcurve = CAL, predcurve = PRED, thresh = "cpD2", 
       dil = COPIES, plot = FALSE) 

## Using a defined treshold value.
#calib(refcurve = CAL, predcurve = PRED, thresh = 0.5, dil = COPIES) 

## Using six dilutions with four replicates/dilution.
## Not run: 
#CAL2 <- modlist(reps, fluo = 2:25)
#calib(refcurve = CAL2, predcurve = PRED, thresh = "cpD2", 
#      dil = COPIES, group = gl(6,4)) 

## End(Not run)

---

Cy0 alternative to threshold cycles as in Guescini et al. (2008)

Description

An alternative to the classical crossing point/threshold cycle estimation as described in Guescini et al (2002). A tangent is fit to the first derivative maximum (point of inflection) of the modeled curve and the intersection with the x-axis is calculated.

Usage

Cy0(object, plot = FALSE, add = FALSE, ...)

Arguments

object	a fitted object of class 'pcrfit'.
plot	if TRUE, displays a plot of Cy0.
add	if TRUE, a plot is added to any other existing plot, i.e. as from plot.pcrfit.
...	other parameters to be passed to plot.pcrfit or points.

Details

The function calculates the first derivative maximum (cpD1) of the curve and the slope and fluorescence $F_{cpD2}$ at that point. Cy0 is then calculated by $Cy_0 = cpD1 - \frac{F_{cpD2}}{slope}$.

Value

The $Cy_0$ value.

Author(s)

Andrej-Nikolai Spiess

References

A new real-time PCR method to overcome significant quantitative inaccuracy due to slight amplification inhibition.
Guescini M, Sisti D, Rocchi MB, Stocchi L & Stocchi V.
BMC Bioinformatics (2008), 9: 326.

Examples

## Single curve with plot.
m1 <- pcrfit(reps, 1, 2, l5)
Cy0(m1, plot = TRUE)

## Add to 'efficiency' plot.
efficiency(m1)
Cy0(m1, add = TRUE)

## Compare s.d. of replicates between
## Cy0 and cpD2 method. cpD2 wins!
ml1 <- modlist(reps, model = l4)
cy0 <- sapply(ml1, function(x) Cy0(x))
cpd2 <- sapply(ml1, function(x) efficiency(x, plot = FALSE)$cpD2)
tapply(cy0, gl(7, 4), function(x) sd(x))
tapply(cpd2, gl(7, 4), function(x) sd(x))

---

The amplification efficiency curve of a fitted object

Description

Calculates the efficiency curve from the fitted object by $E_n = \frac{F_n}{F_{n-1}}$, with $E$ = efficiency, $F$ = raw fluorescence, $n$ = Cycle number. Alternatively, a cubic spline interpolation can be used on the raw data as in Shain et al. (2008).

Usage

eff(object, method = c("sigfit", "spline"), sequence = NULL, baseshift = NULL, 
    smooth = FALSE, plot = FALSE)

Arguments

object	an object of class 'pcrfit'.
method	the efficiency curve is either calculated from the sigmoidal fit (default) or a cubic spline interpolation.
sequence	a 3-element vector (from, to, by) defining the sequence for the efficiency curve. Defaults to [min(Cycles), max(Cycles)] with 100 points per cycle.
baseshift	baseline shift value in case of type = "spline". See documentation to maxRatio.
smooth	logical. If TRUE and type = "spline", invokes a 5-point convolution filter (filter). See documentation to maxRatio.
plot	should the efficiency be plotted?

Details

For more information about the curve smoothing, baseline shifting and cubic spline interpolation for the method as in Shain et al. (2008), see 'Details' in maxRatio.

Value

A list with the following components:

eff.x	the cycle points.
eff.y	the efficiency values at eff.x.
effmax.x	the cycle number with the highest efficiency.
effmax.y	the maximum efficiency.

Author(s)

Andrej-Nikolai Spiess

References

A new method for robust quantitative and qualitative analysis of real-time PCR.
Shain EB & Clemens JM.
Nucleic Acids Research (2008), 36, e91.

Examples

## With default 100 points per cycle.
m1 <- pcrfit(reps, 1, 7, l5)
eff(m1, plot = TRUE) 

## Not all data and only 10 points per cycle.
eff(m1, sequence = c(5, 35, 0.1), plot = TRUE) 

## When using cubic splines it is preferred 
## to use the smoothing option.
#eff(m1, method = "spline", plot = TRUE, smooth = TRUE, baseshift = 0.3)

---

Calculation of qPCR efficiency and other important qPCR parameters

Description

This function calculates the PCR efficiency of a model of class 'pcrfit', including several other important values for qPCR quantification like the first and second derivatives and the corresponding maxima thereof (i.e. threshold cycles). These values can subsequently be used for the calculation of PCR kinetics, fold induction etc. All values are included in a graphical output of the fit. Additionally, several measures of goodness-of-fit are calculated, i.e. the Akaike Information Criterion (AIC), the residual variance and the $R^2$ value.

Usage

efficiency(object, plot = TRUE, type = "cpD2", thresh = NULL, 
           shift = 0, amount = NULL, ...)

Arguments

object	an object of class 'pcrfit'.
plot	logical. If TRUE, a graph is displayed. If FALSE, values are printed out.
type	the method of efficiency estimation. See 'Details'.
thresh	an (optional) numeric value for a fluorescence threshold border. Overrides type.
shift	a user defined shift in cycles from the values defined by type. See 'Examples'.
amount	the template amount or molecule number for quantitative calibration.
...	other parameters to be passed to eff or plot.pcrfit.

Details

The efficiency is always (with the exception of type = "maxRatio") calculated from the efficiency curve (in blue), which is calculated according to $E_n = \frac{F_n}{F_{n-1}}$ from the fitted curve, but taken from different points at the curve, as to be defined in type:

"cpD2" taken from the maximum of the second derivative curve,
"cpD1" taken from the maximum of the first derivative curve,
"maxE" taken from the maximum of the efficiency curve,
"expR" taken from the exponential region by $expR = cpD2-(cpD1-cpD2)$,
"CQ" taken from the 20% value of the fluorescence at "cpD2" as developed by Corbett Research (comparative quantification),
"Cy0" the intersection of a tangent on the first derivative maximum with the abscissa as calculated according to Guescini et al. or
a numeric value taken from the threshold cycle output of the PCR software, i.e. 15.24 as defined in type or
a numeric value taken from the fluorescence threshold output of the PCR software as defined in thresh.

The initial fluorescence $F_0$ for relative or absolute quantification is either calculated by setting $x = 0$ in the sigmoidal model of object giving init1 or by calculating an exponential model down (init2) with $F_0 = \frac{F_n}{E^n}$, with $F_n$ = raw fluorescence, $E$ = PCR efficiency and $n$ = the cycle number defined by type. If a template amount is defined, a conversion factor $cf = \frac{amount}{F_0}$ is given. The different measures for goodness-of-fit give an overview for the validity of the efficiency estimation. First and second derivatives are calculated from the fitted function and the maxima of the derivatives curve and the efficiency curve are obtained.

If type = "maxRatio", the maximum efficiency is calculated from the cubic spline interpolated raw fluorescence values and therefore NOT from the sigmoidal fit. This is a different paradigm and will usually result in fairly the same threshold cycles as with type = "cpD2", but the efficiencies are generally lower. See documentation to maxRatio. This method is usually not applied for calculating efficiencies that are to be used for relative quantification, but one might try...

Value

A list with the following components:

eff	the PCR efficiency.
resVar	the residual variance.
AICc	the bias-corrected Akaike Information Criterion.
AIC	the Akaike Information Criterion.
Rsq	the $R^2$ value.
Rsq.ad	the adjusted $R_{adj}^2$ value.
cpD1	the first derivative maximum (point of inflection in 'l4' or 'b4' models, can be used for defining the threshold cycle).
cpD2	the second derivative maximum (turning point of cpD1, more often used for defining the threshold cycle).
cpE	the PCR cycle with the highest efficiency.
cpR	the PCR cycle within the exponential region calculated as under 'Details'.
cpT	the PCR cycle corresponding to the fluorescence threshold as defined in thresh.
Cy0	the PCR threshold cycle 'Cy0' according to Guescini et al. See 'Details'.
cpCQ	the PCR cycle corresponding to the 20% fluorescence value at 'cpD2'.
cpMR	the PCR cycle corresponding to the 'maxRatio', if this was selected.
fluo	the raw fluorescence value at the point defined by type or thresh.
init1	the initial template fluorescence from the sigmoidal model, calculated as under 'Details'.
init2	the initial template fluorescence from an exponential model, calculated as under 'Details'.
cf	the conversion factor between raw fluorescence and template amount, if the latter is defined.

If object was of type 'modlist', the results are given as a matrix, with samples in columns.

Note

In some curves that are fitted with the 'b5'/'l5' models, the 'f' (asymmetry) parameter can be extremely high due to severe asymmetric structure. The efficiency curve deduced from the coefficients of the fit can then be very extreme in the exponential region. It is strongly advised to use efficiency(object, method = "spline") so that eff calculates the curve from a cubic spline of the original data points (see 'Examples').

Three parameter models ('b3' or 'l3') do not work very well in calculating the PCR efficiency. It is advisable not to take too many cycles of the plateau phase prior to fitting the model as this has a strong effect on the validity of the efficiency estimates.

Author(s)

Andrej-Nikolai Spiess

References

Validation of a quantitative method for real time PCR kinetics.
Liu W & Saint DA.
BBRC (2002), 294: 347-353.

A new real-time PCR method to overcome significant quantitative inaccuracy due to slight amplification inhibition.
Guescini M, Sisti D, Rocchi MB, Stocchi L & Stocchi V.
BMC Bioinformatics (2008), 9: 326.

Examples

## Fitting initial model.
m1 <-  pcrfit(reps, 1, 2, l4)
efficiency(m1)
 
## Using one cycle 'downstream'
## of second derivative max.
efficiency(m1, type = "cpD2", shift = -1)

## Using "maxE" method, with calculation of PCR efficiency
## 2 cycles 'upstream' from the cycle of max efficiency.
efficiency(m1, type = "maxE", shift = 2)

## Using the exponential region.
efficiency(m1, type = "expR")

## Using threshold cycle (i.e. 15.32) 
## from PCR software.
efficiency(m1, type = 15.32)

## Using Cy0 method from Guescini et al. (2008)
## add Cy0 tangent. 
efficiency(m1, type = "Cy0")
Cy0(m1, add = TRUE)

## Using a defined fluorescence
## threshold value from PCR software.
efficiency(m1, thresh = 1)
 
## Using the first 30 cycles and a template amount
## (optical calibration).
m2 <-  pcrfit(reps[1:30, ], 1, 2, l5)
efficiency(m2, amount = 1E3)

## Using 'maxRatio' method from Shain et al. (2008)
## baseshifting essential!
efficiency(m1, type = "maxRatio", baseshift = 0.2)

## Using the efficiency from a cubic spline fit
## of the 'eff' function.
efficiency(m1, method = "spline")

## Not run: 
## On a modlist with plotting
## of the efficiencies.
ml1 <- modlist(reps, model = l5)
res <- sapply(ml1, function(x) efficiency(x)$eff)
barplot(as.numeric(res))

## End(Not run)

---

Evidence ratio for model comparisons with AIC, AICc or BIC

Description

The evidence ratio

$\frac{1}{exp(-0.5 \cdot (IC2 - IC1))}$

is calculated for one of the information criteria $IC = AIC, AICc, BIC$ either from two fitted models or two numerical values. Models can be compared that are not nested and where the f-test on residual-sum-of-squares is not applicable.

Usage

evidence(x, y, type = c("AIC", "AICc", "BIC"))

Arguments

x	a fitted object or numerical value.
y	a fitted object or numerical value.
type	any of the three Information Criteria AIC, AICc or BIC.

Details

Small differences in values can mean substantial more 'likelihood' of one model over the other. For example, a model with AIC = -130 is nearly 150 times more likely than a model with AIC = -120.

Value

A value of the first model x being more likely than the second model y. If large, first model is better. If small, second model is better.

Author(s)

Andrej-Nikolai Spiess

Examples

## Compare two four-parameter and five-parameter
## log-logistic models.
m1 <- pcrfit(reps, 1, 2, l4)
m2 <- pcrfit(reps, 1, 2, l5)
evidence(m2, m1)

## Ratio of two AIC's.
evidence(-120, -123)

---

Comparison of all sigmodal models within the exponential region

Description

The exponential region of the qPCR data is identified by the studentized outlier method, as in expfit. The root-mean-squared-error (RMSE) of all available sigmoidal models within this region is then calculated. The result of the fits are plotted and models returned in order of ascending RMSE.

Usage

expcomp(object, ...)

Arguments

object	an object of class 'pcrfit'.
...	other parameters to be passed to expfit.

Details

The following sigmoidal models are fitted: b4, b5, b6, b7, l4, l5, l6, l7

Value

A dataframe with names of the models, in ascending order of RMSE.

Author(s)

Andrej-Nikolai Spiess

Examples

m1 <- pcrfit(reps, 1, 2, l4)
expcomp(m1)

---

Calculation of PCR efficiency by fitting an exponential model

Description

An exponential model is fit to a window of defined size on the qPCR raw data. The window is identified either by the second derivative maximum 'cpD2' (default), 'studentized outlier' method as described in Tichopad et al. (2003), the 'midpoint' method (Peirson et al., 2003) or by subtracting the difference of cpD1 and cpD2 from cpD2 ('ERBCP', unpublished).

Usage

expfit(object, method = c("cpD2", "outlier", "midpoint", "ERBCP"),
       model = c("exp", "linexp"), offset = 0, pval = 0.05, n.outl = 3, 
       n.ground = 1:5, corfact = 1, fix = c("top", "bottom", "middle"), 
       nfit = 5, plot = TRUE, ...)

Arguments

object	an object of class 'pcrfit'.
method	one of the four possible methods to be used for defining the position of the fitting window.
model	which exponential model to use. expGrowth is default, but the linear-exponential model linexp can also be chosen.
offset	for method = "cpD2", the cycle offset from second derivative maximum.
pval	for method = "outlier", the p-value for the outlier test.
n.outl	for method = "outlier", the number of successive outlier cycles.
n.ground	for method = "midpoint", the number of cycles in the noisy ground phase to calculate the standard deviation from.
corfact	for method = "ERBCP", the correction factor for finding the exponential region. See 'Details'.
fix	for methods "midpoint" and "ERBCP", the orientation of the fitting window based on the identified point. See 'Details'.
nfit	the size of the fitting window.
plot	logical. If TRUE, a graphical display of the curve and the fitted region is shown.
...	other parameters to be passed to the plotting function.

Details

The exponential growth function $f(x) = a \cdot exp(b \cdot x) + c$ is fit to a subset of the data. Calls efficiency for calculation of the second derivative maximum, takeoff for calculation of the studentized residuals and 'outlier' cycle, and midpoint for calculation of the exponential phase 'midpoint'. For method 'ERBCP' (Exponential Region By Crossing Points), the exponential region is calculated by $expR = cpD2 - \code{corfact} \cdot (cpD1-cpD2)$. The efficiency is calculated from the exponential fit with $E = exp(b)$ and the inital template fluorescence $F_0 = a$.

Value

A list with the following components:

point	the point within the exponential region as identified by one of the three methods.
cycles	the cycles of the identified region.
eff	the efficiency calculated from the exponential fit.
AIC	the Akaike Information Criterion of the fit.
resVar	the residual variance of the fit.
RMSE	the root-mean-squared-error of the fit.
init	the initial template fluorescence.
mod	the exponential model of class 'nls'.

Author(s)

Andrej-Nikolai Spiess

References

Standardized determination of real-time PCR efficiency from a single reaction set-up.
Tichopad A, Dilger M, Schwarz G & Pfaffl MW.
Nucleic Acids Research (2003), 31:e122.

Comprehensive algorithm for quantitative real-time polymerase chain reaction.
Zhao S & Fernald RD.
J Comput Biol (2005), 12:1047-64.

Examples

## Using default SDM method.
m1 <- pcrfit(reps, 1, 2, l5)
expfit(m1)

## Using 'outlier' method.
expfit(m1, method = "outlier")

## Linear exponential model.
expfit(m1, model = "linexp")

---

The chi-square goodness-of-fit

Description

Calculates $\chi^2$, reduced $\chi_{\nu}^2$ and the $\chi^2$ fit probability for objects of class pcrfit, lm, glm, nls or any other object with a call component that includes formula and data. The function checks for replicated data (i.e. multiple same predictor values). If replicates are not given, the function needs error values, otherwise NA's are returned.

Usage

fitchisq(object, error = NULL)

Arguments

object	a single model of class 'pcrfit', a 'replist' or any fitted model of the above.
error	in case of a model without replicates, a single error for all response values or a vector of errors for each response value.

Details

The variance of a fit $s^2$ is also characterized by the statistic $\chi^2$ defined as followed:

$\chi^2 \equiv \sum_{i=1}^n \frac{(y_i - f(x_i))^2}{\sigma_i^2}$

The relationship between $s^2$ and $\chi^2$ can be seen most easily by comparison with the reduced $\chi^2$:

$\chi_\nu^2 = \frac{\chi^2}{\nu} = \frac{s^2}{\langle \sigma_i^2 \rangle}$

whereas $\nu$ = degrees of freedom (N - p), and $\langle \sigma_i^2 \rangle$ is the weighted average of the individual variances. If the fitting function is a good approximation to the parent function, the value of the reduced chi-square should be approximately unity, $\chi_\nu^2 = 1$. If the fitting function is not appropriate for describing the data, the deviations will be larger and the estimated variance will be too large, yielding a value greater than 1. A value less than 1 can be a consequence of the fact that there exists an uncertainty in the determination of $s^2$, and the observed values of $\chi_\nu^2$ will fluctuate from experiment to experiment. To assign significance to the $\chi^2$ value, we can use the integral probability

$P_\chi(\chi^2;\nu) = \int_{\chi^2}^\infty P_\chi(x^2, \nu)dx^2$

which describes the probability that a random set of n data points sampled from the parent distribution would yield a value of $\chi^2$ equal to or greater than the calculated one. This is calculated by $1 - pchisq(\chi^2, \nu)$.

Value

A list with the following items:

chi2	the $\chi^2$ value.
chi2.red	the reduced $\chi_\nu^2$.
p.value	the fit probability as described above.

Author(s)

Andrej-Nikolai Spiess

References

Data Reduction and Error Analysis for the Physical Sciences.
Bevington PR & Robinson DK.
McGraw-Hill, New York (2003).

Applied Regression Analysis.
Draper NR & Smith H.
Wiley, New York, 1998.

Examples

## Using replicates by making a 'replist'.
ml1 <- modlist(reps, fluo = 2:5)
rl1 <- replist(ml1, group = c(1, 1, 1, 1))
fitchisq(rl1[[1]])

## Using single model with added error.
m1 <- pcrfit(reps, 1, 2, l5)
fitchisq(m1, 0.1)

---

Batch calculation of qPCR fit parameters/efficiencies/threshold cycles with simple output, especially tailored to high-throughput data

Description

This is a cut-down version of pcrbatch, starting with data of class 'modlist', which delivers a simple dataframe output, with either the parameters of the fit or calculated threshold cycles/efficiencies. The column names are deduced from the run names. All calculations have been error-protected through tryCatch, so whenever there is any kind of error (parameter extraction, efficiency estimation etc), NA is returned. This function can be used with high throughput data quite conveniently. All methods as in pcrbatch are available. The results are automatically copied to the clipboard.

Usage

getPar(x, type = c("fit", "curve"), cp = "cpD2", eff = "sigfit", ...)

Arguments

x	an object of class 'pcrfit' or 'modlist'.
type	fit will extract the fit parameters, curve will invoke efficiency and return threshold cycles/efficiencies.
cp	which method for threshold cycle estimation. Any of the methods in efficiency, i.e. "cpD2" (default), "cpD1", "maxE", "expR", "Cy0", "CQ", "maxRatio".
eff	which method for efficiency estimation. Either "sigfit" (default), "sliwin" or "expfit".
...	other parameters to be passed to efficiency, sliwin or expfit.

Details

Takes about 4 sec for 100 runs on a Pentium 4 Quad-Core (3 Ghz) when using type = "curve". When using type = "fit", the fitted model parameters are returned. If type = "curve", threshold cycles and efficiencies are calculated by efficiency based on the parameters supplied in ... (default cpD2).

Value

A dataframe, which is automatically copied to the clipboard.

Author(s)

Andrej-Nikolai Spiess.

Examples

## Simple example with fit parameters.
ml1 <- modlist(rutledge, model = l5)
getPar(ml1, type = "fit")

## Using a mechanistic model such as
## 'mak3' and extracting D0 values
## => initial template fluorescence.
ml2 <- modlist(rutledge, 1, 2:41, model = mak3)
res <- getPar(ml2, type = "fit")
barplot(log10(res[1, ]), las = 2)

---

Outlier summary for objects of class 'modlist' or 'replist'

Description

For model lists of class 'modlist' or 'replist', is.outlier returns a vector of logicals for each run if they are outliers (i.e. sigmoidal or kinetic) or not.

Usage

is.outlier(object)

Arguments

object

an object of class 'modlist' or 'replist'.

Value

A vector of logicals with run names.

Author(s)

Andrej-Nikolai Spiess

See Also

KOD.

Examples

## Analyze in respect to amplification
## efficiency outliers.
ml1 <- modlist(reps, 1, 2:5)
res1 <- KOD(ml1, check = "uni2")

## Which runs are outliers?
outl <- is.outlier(res1)
outl
which(outl)

## Not run: 
## Test for sigmoidal outliers
## with the 'testdat' dataset.
ml2 <- modlist(testdat, model = l5, check = "uni2")
is.outlier(ml2)    

## End(Not run)

---

(K)inetic (O)utlier (D)etection using several methods

Description

Identifies and/or removes qPCR runs according to several published methods or own ideas. The univariate measures are based on efficiency or difference in first/second derivative maxima. Multivariate methods are implemented that describe the structure of the curves according to several fixpoints such as first/second derivative maximum, slope at first derivative maximum or plateau fluorescence. These measures are compared with a set of curves using the mahalanobis distance with a robust covariance matrix and calculation of statistics by a $\chi^2$ distribution. See 'Details'.

Usage

KOD(object, method = c("uni1", "uni2", "multi1", "multi2", "multi3"),
    par = parKOD(), remove = FALSE, verbose = TRUE, plot = TRUE,  ...)

Arguments

object	an object of class 'modlist' or 'replist'.
method	which method to use for kinetic outlier identification. Method "uni1" is default. See 'Details' for all methods.
par	parameters for the different methods. See parKOD.
remove	logical. If TRUE, outlier runs are removed and the object is updated. If FALSE, the individual qPCR runs are tagged as 'outliers' or not. See 'Details'.
verbose	logical. If TRUE, all calculation steps and results are displayed on the console.
plot	logical. If TRUE, a multivariate plot is displayed.
...	any other parameters to be passed to sliwin, efficiency or expfit.

Details

The following methods for the detection of kinetic outliers are implemented
uni1: KOD method according to Bar et al. (2003). Outliers are defined by removing the sample efficiency from the replicate group and testing it against the remaining samples' efficiencies using a Z-test:

$P = 2 \cdot \left[1 - \Phi\left(\frac{e_i - \mu_{train}}{\sigma_{train}}\right)\right] < 0.05$

uni2: This method from the package author is more or less a test on sigmoidal structure for the individual curves. It is different in that there is no comparison against other curves from a replicate set. The test is simple: The difference between first and second derivative maxima should be less than 10 cycles:

$\left(\frac{\partial^3 F(x;a,b,...)}{\partial x^3} = 0\right) - \left(\frac{\partial^2 F(x;a,b...)}{\partial x^2} = 0\right) < 10$

Sounds astonishingly simple, but works: Runs are defines as 'outliers' that really failed to amplify, i.e. have no sigmoidal structure or are very shallow. It is the default setting in modlist.

multi1: KOD method according to Tichopad et al. (2010). Assuming two vectors with first and second derivative maxima $t_1$ and $t_2$ from a 4-parameter sigmoidal fit within a window of points around the first derivative maximum, a linear model $t_2 = t_1 \cdot b + a + \tau$ is made. Both $t_1$ and the residuals from the fit $\tau = t_2 - \hat{t_2}$ are Z-transformed:

$t_1(norm) = \frac{t_1 - \bar{t}_1}{{\sigma_t}_1}, \; {\tau_1}_{norm} = \frac{\tau_1 - \bar{\tau}_1}{{\sigma_\tau}_1}$

Both $t_1$ and $\tau$ are used for making a robust covariance matrix. The outcome is plugged into a mahalanobis distance analysis using the 'adaptive reweighted estimator' from package 'mvoutlier' and p-values for significance of being an 'outlier' are deduced from a $\chi^2$ distribution. If more than two parameters are supplied, princomp is used instead.

multi2: Second KOD method according to Tichopad et al. (2010), mentioned in the paper. Uses the same pipeline as multi1, but with the slope at the first derivative maximum and maximum fluorescence as parameters:

$\frac{\partial F(x;a,b,...)}{\partial x}, F_{max}$

multi3: KOD method according to Sisti et al. (2010). Similar to multi2, but uses maximum fluorescence, slope at first derivative maximum and y-value at first derivative maximum as fixpoints:

$\frac{\partial F(x;a,b,...)}{\partial x}, F\left(\frac{\partial^2 F(x;a,b,...)}{\partial x^2} = 0\right), F_{max}$

All essential parameters for the methods can be tweaked by parKOD. See there and in 'Examples'.

Value

An object of the same class as in object that is 'tagged' in its name (**name**) if it is an outlier and also with an item $isOutlier with outlier information (see is.outlier). If remove = TRUE, the outlier runs are removed (and the fitting updated in case of a 'replist').

Author(s)

Andrej-Nikolai Spiess

References

Kinetic Outlier Detection (KOD) in real-time PCR.
Bar T, Stahlberg A, Muszta A & Kubista M.
Nucl Acid Res (2003), 31: e105.

Quality control for quantitative PCR based on amplification compatibility test.
Tichopad A, Bar T, Pecen L, Kitchen RR, Kubista M &, Pfaffl MW.
Methods (2010), 50: 308-312.

Shape based kinetic outlier detection in real-time PCR.
Sisti D, Guescini M, Rocchi MBL, Tibollo P, D'Atri M & Stocchi V.
BMC Bioinformatics (2010), 11: 186.

See Also

Function is.outlier to get an outlier summary.

Examples

## kinetic outliers:
## on a 'modlist', using efficiency from sigmoidal fit
## and alpha = 0.01. 
## F7.3 detected as outlier (shallower => low efficiency)
ml1 <- modlist(reps, 1, c(2:5, 28), model = l5)
res1 <- KOD(ml1, method = "uni1", par = parKOD(eff = "sliwin", alpha = 0.01))
plot(res1)

## Sigmoidal outliers:
## remove runs without sigmoidal structure.
ml2 <- modlist(testdat, model = l5)
res2 <- KOD(ml2, method = "uni2", remove = TRUE)
plot(res2, which = "single")

## Not run: 
## Multivariate outliers:
## a few runs are identified.
ml3 <- modlist(reps, model = l5)
res3 <- KOD(ml3, method = "multi1")

## On a 'replist', several outliers identified.
rl3 <- replist(ml3, group = gl(7, 4))
res4 <- KOD(rl3, method = "uni1")

## End(Not run)

---

Calculation of likelihood ratios for nested models

Description

Calculates the likelihood ratio and p-value from a chi-square distribution for two nested models.

Usage

llratio(objX, objY)

Arguments

objX	Either a value of class logLik or a model for which logLik can be applied.
objY	Either a value of class logLik or a model for which logLik can be applied.

Details

The likelihood ratio statistic is

$LR = \frac{f(X, \hat{\phi}, \hat{\psi})}{f(X, \phi, \hat{\psi_0})}$

The usual test statistic is

$\Lambda = 2 \cdot (l(\hat{\phi}, \hat{\psi}) - l(\phi, \hat{\psi_0}))$

Following the large sample theory, if $H_0$ is true, then

$\Lambda \sim \chi_p^2$

Value

A list containing the following items:

ratio	the likelihood ratio statistic.
df	the change in parameters.
p.value	the p-value from a $\chi^2$ distribution. See Details.

Author(s)

Andrej-Nikolai Spiess

See Also

AIC, logLik.

Examples

## Compare l5 and l4 model.
m1 <- pcrfit(reps, 1, 2, l5)
m2 <- pcrfit(reps, 1, 2, l4)
llratio(m1, m2)

---

Formal lack-Of-Fit test of a nonlinear model against a one-way ANOVA model

Description

Tests the nonlinear model against a more general one-way ANOVA model and from a likelihood ratio test. P-values are derived from the F- and $\chi^2$ distribution, respectively.

Usage

LOF.test(object)

Arguments

object

an object of class 'replist', 'pcrfit' or 'nls', which was fit with replicate response values.

Details

The one-way ANOVA model is constructed from the data component of the nonlinear model by factorizing each of the predictor values. Hence, the nonlinear model becomes a submodel of the one-way ANOVA model and we test both models with the null hypothesis that the ANOVA model can be simplified to the nonlinear model (Lack-of-fit test). This is done by two approaches:

1) an F-test (Bates & Watts, 1988).
2) a likelihood ratio test (Huet et al, 2004).

P-values are derived from an F-distribution (1) and a $\chi^2$ distribution (2).

Value

A list with the following components:

pF	the p-value from the F-test against the one-way ANOVA model.
pLR	the p-value from the likelihood ratio test against the one-way ANOVA model.

Author(s)

Andrej-Nikolai Spiess

References

Nonlinear Regression Analysis and its Applications.
Bates DM & Watts DG.
John Wiley & Sons (1988), New York.

Statistical Tools for Nonlinear Regression: A Practical Guide with S-PLUS and R Examples.
Huet S, Bouvier A, Poursat MA & Jolivet E.
Springer Verlag (2004), New York, 2nd Ed.

Examples

## Example with a 'replist'
## no lack-of-fit.
ml1 <- modlist(reps, fluo = 2:5, model = l5)
rl1 <- replist(ml1, group = c(1, 1, 1, 1))
LOF.test(rl1)

## Example with a 'nls' fit
## => there is a lack-of-fit.
DNase1 <- subset(DNase, Run == 1)
fm1DNase1 <- nls(density ~ SSlogis(log(conc), Asym, xmid, scal), DNase1) 
LOF.test(fm1DNase1)

---

Calculation of qPCR efficiency by the 'linear regression of efficiency' method

Description

The LRE method is based on a linear regression of raw fluorescence versus efficiency, with the final aim to obtain cycle dependent individual efficiencies $E_n$. A linear model is then fit to a sliding window of defined size(s) and within a defined border. Regression coefficients are calculated for each window, and from the window of maximum regression, parameters such as PCR efficiency and initial template fluorescence are calculated. See 'Details' for more information. This approach is quite similar to the one in sliwin, but while sliwin regresses cycle number versus log(fluorescence), LRE regresses raw fluorescence versus efficiency. Hence, the former is based on assuming a constant efficiency for all cycles while the latter is based on a per-cycle individual efficiency.

Usage

LRE(object, wsize = 6, basecyc = 1:6, base = 0, border = NULL, 
    plot = TRUE, verbose = TRUE, ...)

Arguments

object	an object of class 'pcrfit'.
wsize	the size(s) of the sliding window(s), default is 6. A sequence such as 4:6 can be used to optimize the window size.
basecyc	if base != 0, which cycles to use for an initial baseline estimation based on the averaged fluorescence values.
base	either 0 for no baseline optimization, or a scalar defining multiples of the standard deviation of all baseline points obtained from basecyc. These are iteratively subtracted from the raw data. See 'Details' and 'Examples'.
border	either NULL (default) or a two-element vector which defines the border from the take-off point to points nearby the upper asymptote (saturation phase). See 'Details'.
plot	if TRUE, the result is plotted with the fluorescence/efficiency curve, sliding window, regression line and baseline.
verbose	logical. If TRUE, more information is displayed in the console window.
...	only used internally for passing the parameter matrix.

Details

To avoid fits with a high $R^2$ in the baseline region, some border in the data must be defined. In LRE, this is by default (base = NULL) the region in the curve starting at the take-off cycle ($top$) as calculated from takeoff and ending at the transition region to the upper asymptote (saturation region). The latter is calculated from the first and second derivative maxima: $asympt = cpD1 + (cpD1 - cpD2)$. If the border is to be set by the user, border values such as c(-2, 4) extend these values by $top + border[1]$ and $asympt + border[2]$. The efficiency is calculated by $E_n = \frac{F_n}{F_{n-1}}$ and regressed against the raw fluorescence values $F$: $E = F\beta + \epsilon$. For the baseline optimization, 100 baseline values $Fb_i$ are interpolated in the range of the data:

$F_{min} \le Fb_i \le base \cdot \sigma(F_{basecyc[1]}...F_{basecyc[2]})$

and subtracted from $F_n$. For all iterations, the best regression window in terms of $R^2$ is found and its parameters returned. Two different initial template fluorescence values $F_0$ are calculated in LRE:

init1: Using the single maximum efficiency $E_{max}$ (the intercept of the best fit) and the fluorescence at second derivative maximum $F_{cpD2}$, by

$F_0 = \frac{F_{cpD2}}{E_{max}^{cpD2}}$

init2: Using the cycle dependent efficiencies $E_n$ from $n = 1$ to the near-lowest integer (floor) cycle of the second derivative maximum $n = \lfloor cpD2 \rfloor$, and the fluorescence at the floor of the second derivative maximum $F_{\lfloor cpD2 \rfloor}$, by

$F_0 = \frac{F_{\lfloor cpD2 \rfloor}}{\prod E_n}$

This approach corresponds to the paradigm described in Rutledge & Stewart (2008), by using cycle-dependent and decreasing efficiencies $\Delta_E$ to calculate $F_0$.

Value

A list with the following components:

eff	the maximum PCR efficiency $E_{max}$ calculated from the best window.
rsq	the maximum $R^2$.
base	the optimized baseline value.
window	the best window found within the borders.
parMat	a matrix containing the parameters as above for each iteration.
init1	the initial template fluorescence $F_0$ assuming constant efficiency $E_{max}$ as described under 'Details'.
init2	the initial template fluorescence $F_0$, assuming cycle-dependent efficiency $E_n$ as described under 'Details'.

Author(s)

Andrej-Nikolai Spiess

References

A kinetic-based sigmoidal model for the polymerase chain reaction and its application to high-capacity absolute quantitative real-time PCR.
Rutledge RG & Stewart D.
BMC Biotech (2008), 8: 47.

Examples

## Not run: 
## Sliding window of size 5 between take-off point 
## and 3 cycles upstream of the upper asymptote 
## turning point, one standard deviation baseline optimization.
m1 <- pcrfit(reps, 1, 2, l4)
LRE(m1, wsize = 5, border = c(0, 3), base = 1)

## End(Not run)

---

The maxRatio method as in Shain et al. (2008)

Description

The maximum ratio (MR) is determined along the cubic spline interpolated curve of $\frac{F_n}{F_{n-1}}$ and the corresponding cycle numbers FCN and its adjusted version FCNA are then calculated for MR.

Usage

maxRatio(x, method = c("spline", "sigfit"), baseshift = NULL, 
         smooth = TRUE, plot = TRUE, ...)

Arguments

x	an object of class 'pcrfit' (single run) or 'modlist' (multiple runs).
method	the parameters are either calculated from the cubic spline interpolation (default) or a sigmoidal fit.
baseshift	numerical. Shift value in case of type = "spline". See 'Details'.
smooth	logical. If TRUE and type = "spline", invokes a 5-point convolution filter (filter). See 'Details'.
plot	Should diagnostic plots be displayed?
...	other parameters to be passed to eff or plot.

Details

In a first step, the raw fluorescence data can be smoothed by a 5-point convolution filter. This is optional but feasible for many qPCR setups with significant noise in the baseline region, and therefore set to TRUE as default. If baseshift is a numeric value, this is added to each response value $F_n = F_n + baseshift$ (baseline shifting). Finally, a cubic spline is fit with a resolution of 0.01 cycles and the maximum ratio (efficiency) is calculated by $MR = max(\frac{F_n}{F_{n-1}}-1)$. $FCN$ is then calculated as the cycle number at $MR$ and from this the adjusted $FCNA = FCN -log_2(MR)$. Sometimes problems are encountered in which, due to high noise in the background region, randomly high efficiency ratios are calculated. This must be resolved by tweaking the baseshift value.

Value

A list with the following components:

eff	the maximum efficiency. Equals to mr + 1.
mr	the maximum ratio.
fcn	the cycle number at mr.
fcna	an adjusted fcn, as described in Shain et al.
names	the names of the runs as taken from the original dataframe.

Note

This function has been approved by the original author (Eric Shain).

Author(s)

Andrej-Nikolai Spiess

References

A new method for robust quantitative and qualitative analysis of real-time PCR.
Shain EB & Clemens JM.
Nucleic Acids Research (2008), 36: e91.

Examples

## On single curve using baseline shifting.
m1 <- pcrfit(reps, 1, 2, l5)
maxRatio(m1, baseshift = 0.3)     

## On a 'modlist' using baseline shifting.
## Not run: 
ml1 <- modlist(reps, model = l5) 
maxRatio(ml1, baseshift = 0.5)

## End(Not run)

---

Melting curve analysis with (iterative) Tm identification and peak area calculation/cutoff

Description

This function conducts a melting curve analysis from the melting curve data of a real-time qPCR instrument. The data has to be preformatted in a way that for each column of temperature values there exists a corresponding fluorescence value column. See edit(dyemelt) for a proper format. The output is a graph displaying the raw fluorescence curve (black), the first derivative curve (red) and the identified melting peaks. The original data together with the results ($-\frac{\partial F}{\partial T}$ values, $T_m$ values) are returned as a list. An automatic optimization procedure is also implemented which iterates over span.smooth and span.peaks values and finds the optimal parameter combination that delivers minimum residual sum-of-squares of the identified $T_m$ values to known $T_m$ values. For all peaks, the areas can be calculated and only those included which have areas higher than a given cutoff (cut.Area). If no peak was identified meeting the cutoff values, the melting curves are flagged with a 'bad' attribute. See 'Details'.

Usage

meltcurve(data, temps = NULL, fluos = NULL, window = NULL, 
          norm = FALSE, span.smooth = 0.05, span.peaks = 51, 
          is.deriv = FALSE, Tm.opt = NULL, Tm.border = c(1, 1), 
          plot = TRUE, peaklines = TRUE, calc.Area = TRUE, 
          plot.Area = TRUE, cut.Area = 0,...)

Arguments

data	a dataframe containing the temperature and fluorescence data.
temps	a vector of column numbers reflecting the temperature values. If NULL, they are assumed to be 1, 3, 5, ... .
fluos	a vector of column numbers reflecting the fluorescence values. If NULL, they are assumed to be 2, 4, 6, ... .
window	a user-defined window for the temperature region to be analyzed. See 'Details'.
norm	logical. If TRUE, the fluorescence values are scaled between [0, 1].
span.smooth	the window span for curve smoothing. Can be tweaked to optimize $T_m$ identification.
span.peaks	the window span for peak identification. Can be tweaked to optimize $T_m$ identification. Must be an odd number.
is.deriv	logical. Use TRUE, if data is already in first derivative transformed format.
Tm.opt	a possible vector of known $T_m$ values to optimize span.smooth and span.peaks against. See 'Details' and 'Examples'.
Tm.border	for peak area calculation, a vector containing left and right border temperature values from the $T_m$ values. Default is -1/+1 ?C.
plot	logical. If TRUE, a plot with the raw melting curve, derivative curve and identified $T_m$ values is displayed for each sample.
peaklines	logical. If TRUE, lines that show the identified peaks are plotted.
calc.Area	logical. If TRUE, all peak areas are calculated.
plot.Area	logical. If TRUE, the baselined area identified for the peaks is plotted by filling the peaks in red.
cut.Area	a peak area value to identify only those peaks with a higher area.
...	other parameters to be passed to plot.

Details

The melting curve analysis is conducted with the following steps:

1a) Temperature and fluorescence values are selected in a region according to window.
1b) If norm = TRUE, the fluorescence data is scaled into [0, 1] by qpcR:::rescale.
Then, the function qpcR:::TmFind conducts the following steps:
2a) A cubic spline function (splinefun) is fit to the raw fluorescence melt values.
2b) The first derivative values are calculated from the spline function for each of the temperature values.
2c) Friedman's supersmoother (supsmu) is applied to the first derivative values.
2d) Melting peaks ($T_m$) values are identified by qpcR:::peaks.
2e) Raw melt data, first derivative data, best parameters, residual sum-of-squares and identified $T_m$ values are returned.
Peak areas are then calculated by qpcR:::peakArea:
3a) A linear regression curve is fit from the leftmost temperature value ($T_m$ - Tm.border[1]) to the rightmost temperature value ($T_m$ + Tm.border[2]) by lm.
3b) A baseline curve is calculated from the regression coefficients by predict.lm.
3c) The baseline data is subtracted from the first derivative melt data (baselining).
3d) A splinefun is fit to the baselined data.
3e) The area of this spline function is integrated from the leftmost to rightmost temperature value.
4) If calculated peak areas were below cut.Area, the corresponding $T_m$ values are removed.
Finally,
5) A matrix of xyy-plots is displayed using qpcR:::xyy.plot.

is.deriv must be set to TRUE if the exported data was already transformed to $-\frac{\partial F}{\partial T}$ by the PCR system (i.e. Stratagene MX3000P).

If values are given to Tm.opt (see 'Examples'), then meltcurve is iterated over all combinations of span.smooth = seq(0, 0.2, by = 0.01) and span.peaks = seq(11, 201, by = 10). For each iteration, $T_m$ values are calculated and compared to those given by measuring the residual sum-of-squares between the given values Tm.opt and the $Tm$ values obtained during the iteration:

$RSS = \sum_{i=1}^n{(Tm_i - Tm.opt_i)^2}$

The returned list items containing the resulting data frame each has an attribute "quality" which is set to "bad" if none of the peaks met the cut.Area criterion (or "good" otherwise).

Value

A list with as many items as melting curves, named as in data, each containing a data.frame with the temperature (Temp), fluorescence values (Fluo), first derivative (dF.dT) values, (optimized) parameters of span.smooth/span.peaks, residual sum-of-squares (if Tm.opt != NULL), identified melting points (Tm), calculated peak areas (Area) and peak baseline values (baseline).

Note

The peaks function is derived from a R-Help mailing list entry in Nov 2005 by Martin Maechler.

Author(s)

Andrej-Nikolai Spiess

Examples

## Default columns.
data(dyemelt)
res1 <- meltcurve(dyemelt, window = c(75, 86))
res1

## Selected columns and normalized fluo values.
res2 <- meltcurve(dyemelt, temps = c(1, 3), fluos = c(2, 4), 
                  window = c(75, 86), norm = TRUE)  

## Removing peaks based on peak area
## => two peaks have smaller areas and are not included.
res3 <- meltcurve(dyemelt, temps = 1, fluos = 2, window = c(75, 86),  
                  cut.Area = 0.2) 
attr(res3[[1]], "quality")
                 
## If all peak areas do not meet the cutoff value, meltcurve is
## flagged as 'bad'.
res4 <- meltcurve(dyemelt, temps = 1, fluos = 2, window = c(75, 86),  
                  cut.Area = 0.5) 
attr(res4[[1]], "quality")

## Optimizing span and peaks values.
## Not run: 
res5 <- meltcurve(dyemelt[, 1:6], window = c(74, 88), 
                  Tm.opt = c(77.2, 80.1, 82.4, 84.8))

## End(Not run)

---

Calculation of the 'midpoint' region

Description

Calculates the exponential region midpoint using the algorithm described in Peirson et al. (2003).

Usage

midpoint(object, noise.cyc = 1:5)

Arguments

object	a fitted object of class 'pcrfit'.
noise.cyc	the cycles defining the background noise.

Details

The 'midpoint' region is calculated by

$F_{noise} \cdot \sqrt{\frac{F_{max}}{F_{noise}}}$

with $F_{noise}$ = the standard deviation of the background cycles and $F_{max}$ = the maximal fluorescence.

Value

A list with the following components:

f.mp	the 'midpoint' fluorescence.
cyc.mp	the 'midpoint' cycle, as predicted from f.mp.

Author(s)

Andrej-Nikolai Spiess

References

Experimental validation of novel and conventional approaches to quantitative real-time PCR data analysis.
Peirson SN, Butler JN & Foster RG.
Nucleic Acids Research (2003), 31: e73.

Examples

m1 <- pcrfit(reps, 1, 2, l5)
mp <- midpoint(m1) 
plot(m1)
abline(h = mp$f.mp, col = 2)
abline(v = mp$mp, col = 2)

---

Create nonlinear models from a dataframe and coerce them into a list

Description

Essential function to create a list of nonlinear models from the columns (runs) of a qPCR dataframe. This function houses different methods for curve transformation prior to fitting, such as normalization in [0, 1], smoothing, baseline subtraction etc. Runs that failed to fit or that have been identified as kinetic outliers (by default: lack of sigmoidal structure) can be removed automatically as well as their entries in an optionally supplied label vector.

Usage

modlist(x, cyc = 1, fluo = NULL, model = l4, check = "uni2", 
        checkPAR = parKOD(), remove = c("none", "fit", "KOD"), 
        exclude = NULL, labels = NULL, norm = FALSE, 
        baseline = c("none", "mean", "median", "lin", "quad", "parm"),
        basecyc = 1:8, basefac = 1, smooth = NULL, smoothPAR = NULL, 
        factor = 1, opt = FALSE, 
        optPAR = list(sig.level = 0.05, crit = "ftest"), verbose = TRUE, ...)

Arguments

x	a dataframe containing the qPCR data or a single qPCR run of class 'pcrfit'.
cyc	the column containing the cycle data. Defaults to first column.
fluo	the column(s) (runs) to be analyzed. If NULL, all runs will be considered.
model	the model to be used for all runs.
check	the method for kinetic outlier detection. Default is check for sigmoidal structure, see KOD. To turn off, use NULL.
checkPAR	parameters to be supplied to the check method, see KOD.
remove	which runs to remove. Either "none", those which failed to "fit" or from the "KOD" outlier method.
exclude	either "" for samples with missing column names or a regular expression defining columns (samples) to be excluded from modlist. See 'Details'.
labels	a vector containing labels, i.e. for defining replicate groups prior to ratiobatch.
norm	logical. Should the raw data be normalized within [0, 1] before model fitting?
baseline	type of baseline subtraction. See 'Details'.
basecyc	cycle range to be used for baseline subtraction, i.e. 1:5.
basefac	a factor for the baseline value, such as 0.95.
smooth	which curve smoothing method to use. See 'Details'.
smoothPAR	parameters to be supplied to the smoothing functions, supplied as a list. See 'Details'.
factor	a multiplication factor for the fluorescence response values (barely useful, but who knows...).
opt	logical. Should model selection be applied to each model?
optPAR	parameters to be supplied to mselect.
verbose	logical. If TRUE, fitting and tagging results will be displayed in the console.
...	other parameters to be passed to pcrfit.

Details

From version 1.4-0, the following baselining methods are available for the fluorescence values:
baseline = numeric: a numeric value such as baseline = 0.2 for subtracting from each $F_i$.
"mean": subtracts the mean of all basecyc cycles from each $F_i$.
"median": subtracts the median of all basecyc cycles from each $F_i$.
"lin": creates a linear model of all basecyc cycles, predicts $P_i$ over all cycles $i$ from this model, and subtracts $F_i - P_i$.
"quad": creates a quadratic model of all basecyc cycles, predicts $P_i$ over all cycles $i$ from this model, and subtracts $F_i - P_i$.
"parm": extracts the $c$ parameter from the fitted sigmoidal model and subtracts this value from all $F_i$.
It is switched off by default, but in case of data with a high baseline (such as in TaqMan PCR), it should be turned on as otherwise this will give highly underestimated efficiencies and hence wrong init2 values.

From version 1.3-8, the following smoothing methods are available for the fluorescence values:
"lowess": Lowess smoothing, see lowess, parameter in smoothPAR: f.
"supsmu": Friedman's SuperSmoother, see supsmu, parameter in smoothPAR: span.
"spline": Smoothing spline, see smooth.spline, parameter in smoothPAR: spar.
"savgol": Savitzky-Golay smoother, qpcR:::savgol, parameter in smoothPAR: none.
"kalman": Kalman smoother, see arima, parameter in smoothPAR: none.
"runmean": Running mean, see qpcR:::runmean, parameter in smoothPAR: wsize.
"whit": Whittaker smoother, see qpcR:::whittaker, parameter in smoothPAR: lambda.
"ema": Exponential moving average, see qpcR:::EMA, parameter in smoothPAR: alpha.
The author of this package advocates the use of "spline", "savgol" or "whit" because these three smoothers have the least influence on overall curve structure.

In case of unsuccessful model fitting and if remove = "none" (default), the original data is included in the output, albeit with no fitting information. This is useful since using plot.pcrfit on the 'modlist' shows the non-fitted runs. If remove = "fit", the non-fitted runs are automatically removed and will thus not be displayed. If remove = "KOD", by default all runs without sigmoidal structure are removed likewise. If a labels vector lab is supplied, the labels from the failed fits are removed and a new label vector lab_mod is written to the global environment. This way, an initial labeling vector for all samples can be supplied, bad runs and their labels automatically removed and these transfered to downstream analysis (i.e. to ratiobatch) without giving errors. exclude offers an option to exclude samples from the modlist by some regular expression or by using "" for samples with empty column names. See 'Examples'.

Value

A list with each item containing the model from each column. A names item (which is tagged by *NAME*, if fitting failed) containing the column name is attached to each model as well as an item isFitted with either TRUE (fitting converged) or FALSE (a fitting error occured). This information is useful when ratiocalc is to be applied and unsuccessful fits should automatically removed from the given group definition. If kinetic outlier detection is selected, an item isOutlier is attached, defining the run as an outlier (TRUE) or not (FALSE).

Author(s)

Andrej-Nikolai Spiess

See Also

pcrbatch for batch analysis using different methods.

Examples

## Calculate efficiencies and ct values 
## for each run in the 'reps' data, 
## subtract baseline using mean of 
## first 8 cycles.
ml1 <- modlist(reps, model = l5, baseline = "mean")
getPar(ml1, type = "curve")

## 'Crossing points' for the first 3 runs (normalized)
##  and using best model from Akaike weights.
ml2 <- modlist(reps, 1, 2:5, model = l5, norm = TRUE, 
               opt = TRUE, optPAR = list(crit = "weights"))
sapply(ml2, function(x) efficiency(x, plot = FALSE)$cpD2)

## Convert a single run to a 'modlist'.
m <- pcrfit(reps, 1, 2, l5)
ml3 <- modlist(m)

## Using the 'testdat' set
## include failed fits.
ml4 <- modlist(testdat, 1, 2:9,  model = l5)
plot(ml4, which = "single")

## Remove failed fits and update a label vector.
GROUP <- c("g1s1", "g1s2", "g1s3", "g1s4", "g1c1", "g1c2", "g1c3", "g1c4") 
ml5 <- modlist(testdat, 1, 2:9,  model = l5, labels = GROUP, remove = "KOD")
plot(ml5, which = "single")

## Smoothing by EMA and alpha = 0.8.
ml6 <- modlist(reps, model = l5, smooth = "ema",
               smoothPAR = list(alpha = 0.5))
plot(ml6)

## Not run: 
## Use one of the mechanistic models
## get D0 values.
ml7 <- modlist(reps, model = mak3)
sapply(ml7, function(x) coef(x)[1])

## Exclude first sample in each 
## replicate group of dataset 'reps'.
ml8 <- modlist(reps, exclude = ".1")
plot(ml8, which = "single")

## Using weighted fitting:
## weighted by inverse residuals.
ml9 <- modlist(reps, weights = "1/abs(resid)")
plot(ml9, which = "single")

## Use linear model of the first 10
## cycles for baselining.
ml10 <- modlist(reps, basecyc = 1:10, baseline = "lin")
plot(ml10)

## Use a single value for baselining.
ml11 <- modlist(reps, basecyc = 1:10, baseline = 0.5)
plot(ml11)

## End(Not run)

---

Sigmoidal model selection by different criteria

Description

Model selection by comparison of different models using

1) the maximum log likelihood value,
2) Akaike's Information Criterion (AIC),
3) bias-corrected Akaike's Information Criterion (AICc),
4) the estimated residual variance,
5) the p-value from a nested F-test on the residual variance,
6) the p-value from the likelihood ratio,
7) the Akaike weights based on AIC,
8) the Akaike weights based on AICc, and
9) the reduced chi-square, $\chi_\nu^2$, if replicates exist.

The best model is chosen by 5), 6), 8) or 9) and returned as a new model.

Usage

mselect(object, fctList = NULL, sig.level = 0.05, verbose = TRUE, 
        crit = c("ftest", "ratio", "weights", "chisq"), do.all = FALSE, ...)

Arguments

object	an object of class 'pcrfit' or 'replist'.
fctList	a list of functions to be analyzed, i.e. for a non-nested regime. Should also contain the original model.
sig.level	the significance level for the nested F-test.
verbose	logical. If TRUE, the result matrix is displayed in the console.
crit	the criterium for model selection. Either "ftest"/"ratio" for nested models or "weights"/"fitprob" for nested and non-nested models.
do.all	if TRUE, all available sigmoidal models are tested and the best one is selected based on AICc weights.
...	other parameters to be passed to fitchisq.

Details

Criteria 5) and 6) cannot be used for comparison unless the models are nested. Criterion 8), Akaike weights, can be used for nested and non-nested regimes, which also accounts for the reduced $\chi_\nu^2$. For criterion 1) the larger the better. For criteria 2), 3) and 4): the smaller the better. The best model is chosen either from the nested F-test (anova), likelihood ratio (llratio), corrected Akaike weights (akaike.weights) or reduced $\chi_\nu^2$ (fitchisq) and returned as a new model. When using "ftest"/"ratio" the corresponding nested functions are analyzed automatically, i.e. b3/b4/b5/b6/b7; l3/l4/l5/l6/l7. If supplying nested models, please do this with ascending number of parameters.

Value

A model of the best fit selected by one of the criteria above. The new model has an additional list item 'retMat' with a result matrix of the criterion tests.

Author(s)

Andrej-Nikolai Spiess

See Also

llratio, akaike.weights and fitchisq.

Examples

## Choose best model based on F-tests 
## on the corresponding nested models.
m1 <- pcrfit(reps, 1, 2, l4)
m2 <- mselect(m1)
summary(m2)  ## Converted to l7 model!

## Use Akaike weights on non-nested models
## compare to original model.
m2 <- mselect(m1, fctList = list(l4, b5, cm3), crit = "weights")
summary(m2) ## Converted to b5 model!

## Try all sigmoidal models.
m3 <- pcrfit(reps, 1, 20, l4)
mselect(m3, do.all = TRUE) ## l7 wins by far!

## On replicated data using reduced chi-square.
ml1 <- modlist(reps, fluo = 2:5, model = l4)
rl1 <- replist(ml1, group = c(1, 1, 1, 1))
mselect(rl1, crit = "chisq")  ## converted to l6!

---

Parameters that can be changed to tweak the kinetic outlier methods

Description

A control function with different list items that change the performance of the different (kinetic) outlier functions as defined in KOD.

Usage

parKOD(eff = c("sliwin", "sigfit", "expfit"), train = TRUE, 
       alpha = 0.05, cp.crit = 10, cut = c(-6, 2))

Arguments

eff	uni1. The efficiency method to be used. Either sliwin, sigfit or expfit.
train	uni1. If TRUE, the sample's efficiency is NOT included in the calculation of the average efficiency (default), if FALSE it is.
cp.crit	uni2. The cycle difference between first and second derivative maxima, default is 10.
cut	multi1. A 2-element vector defining the lower and upper border from the first derivative maximum from where to cut the complete curve.
alpha	the p-value cutoff value for all implemented statistical tests.

Details

For more details on the function of the parameters within the different kinetic and sigmoidal outlier methods, see KOD.

Value

If called, returns a list with the parameters as items.

Author(s)

Andrej-Nikolai Spiess

Examples

## Multivariate outliers,
## adjusting the 'cut' parameter.
ml1 <- modlist(reps, 1, 2:5, model = l5)
res1 <- KOD(ml1, method = "multi1", par = parKOD(cut = c(-5, 2)))

---

Batch calculation of qPCR efficiency and other qPCR parameters

Description

This function batch calculates the results obtained from efficiency, sliwin, expfit, LRE or the coefficients from any of the makX/cm3 models on a dataframe containing many qPCR runs. The input can also be a list obtained from modlist, which simplifies things in many cases. The output is a dataframe with the estimated parameters and model descriptions. Very easy to use on datasheets containing many qPCR runs, i.e. as can be imported from Excel. The result is automatically copied to the clipboard.

Usage

pcrbatch(x, cyc = 1, fluo = NULL, 
         methods = c("sigfit", "sliwin", "expfit", "LRE"),
         model = l4, check = "uni2", checkPAR = parKOD(), 
         remove = c("none", "fit", "KOD"), exclude = NULL, 
         type = "cpD2", labels = NULL, norm = FALSE, 
         baseline = c("none", "mean", "median", "lin", "quad", "parm"),
         basecyc = 1:8, basefac = 1, smooth = NULL, smoothPAR = NULL, 
         factor = 1, opt = FALSE, optPAR = list(sig.level = 0.05, crit = "ftest"), 
         group = NULL, names = c("group", "first"), plot = TRUE, 
         verbose = TRUE, ...)

Arguments

x	a dataframe containing the qPCR raw data from the different runs or a list obtained from modlist.
cyc	the column containing the cycle data. Defaults to first column.
fluo	the column(s) (runs) to be analyzed. If NULL, all runs will be considered.
methods	a character vector defining the methods to use. See 'Details'.
model	the model to be used for all runs.
check	the method for outlier detection in KOD. Default is check for sigmoidal structure.
checkPAR	parameters to be supplied to the check method.
remove	which runs to remove. Either none, those which failed to fit or from the outlier methods.
exclude	either "" for samples with missing column names or a regular expression defining columns (samples) to be excluded from pcrbatch. See 'Details' and 'Examples' in modlist.
type	the point on the amplification curve from which the efficiency is estimated. See efficiency.
labels	a vector containing labels, i.e. for defining replicate groups prior to ratiobatch.
norm	logical. Should the raw data be normalized within [0, 1] before model fitting?
baseline	type of baseline subtraction. See 'Details' in modlist.
basecyc	cycle range to be used for baseline subtraction, i.e. 1:5.
basefac	a factor when using averaged baseline cycles, such as 0.95.
smooth	which curve smoothing method to use. See modlist.
smoothPAR	parameters to be supplied to the smoothing functions, supplied as a list. See modlist.
factor	a multiplication factor for the fluorescence response values (barely useful, but who knows...).
opt	logical. Should model selection be applied to each model?
optPAR	parameters to be supplied to mselect.
group	a vector containing the grouping for possible replicates.
names	how to name the grouped fit. Either 'group_1, ...' or the first name of the replicates.
plot	logical. If TRUE, the single runs are plotted from the internal 'modlist' for diagnostics.
verbose	logical. If TRUE, fitting and tagging results will be displayed in the console.
...	other parameters to be passed to downstream methods.

Details

The methods vector is used for defining the different methods from which pcrbatch will concatenate the results. The mechanistic models (mak2, mak2i, mak3, mak3i, cm3) are omitted by default, because fitting is time-expensive. If they should be included, just add "mak3" to methods. See 'Examples'. The qPCR raw data should be arranged with the cycle numbers in the first column with the name "Cycles". All subsequent columns must be plain raw data with sensible column descriptions. If replicates are defined by group, the output will contain a numbering of groups (i.e. "group_1" for the first replicate group). The model selection process is optional, but we advocate using this for obtaining better parameter estimates. Normalization has been described to improve certain qPCR analyses, but this has still to be independently evaluated. Background subtraction is done as in modlist and efficiency. In case of unsuccessful model fitting or lack of sigmoidal structure, the names are tagged by *NAME* or **NAME**, respectively (if remove = "none"). However, if remove = "fit" or remove = "KOD", the failed runs are excluded from the output. Similar to modlist, if a labels vector lab is supplied, the labels from the failed fits are removed and a new label vector lab_mod is written to the global environment.

Value

A dataframe with the results in columns containing the calculated values, fit parameters and (tagged) model name together with the different methods used as the name prefix. A plot shows a plot matrix of all amplification curves/sigmoidal fits and failed amplifications marked with asterisks.

Note

IMPORTANT: When subsequent use of ratiocalc is desired, use pcrbatch on the single run level with group = NULL and remove = "none", so that ratiocalc can automatically delete the failed runs from its group definition. Otherwise error propagation will fail.

Author(s)

Andrej-Nikolai Spiess

See Also

The function modlist for creating a list of models, which is used internally by pcrbatch.

Examples

## First 4 runs and return parameters of fit
## do background subtraction using mean the first 5 cycles.
pcrbatch(reps, fluo = 2:5, baseline = "mean", basecyc = 1:5)

## Not run: 
##  First 8 runs, with 4 replicates each, l5 model.
pcrbatch(reps, fluo = 2:9, model = l5, group = c(1,1,1,1,2,2,2,2))

## Using model selection (Akaike weights) 
## on the first 4 runs, runs 1 and 2 are replicates.
pcrbatch(reps, fluo = 2:5, group = c(1,1,2,3), 
         opt = TRUE, optPAR = list(crit = "weights"))

## Fitting a sigmoidal and 'mak3' mechanistic model.
pcrbatch(reps, methods = c("sigfit", "mak3"))

## Converting a 'modlist' to 'pcrbatch'.
ml5 <- modlist(reps, 1, 2:5, b5)
res5 <- pcrbatch(ml5)

## Using Whittaker smoothing.
pcrbatch(reps, smooth = "whit")

## End(Not run)

---

Bootstrapping and jackknifing qPCR data

Description

Confidence intervals for the estimated parameters and goodness-of-fit measures are calculated for a nonlinear qPCR data fit by either
a) boostrapping the residuals of the fit or
b) jackknifing and refitting the data.

Confidence intervals can also be calculated for all parameters obtained from the efficiency analysis.

Usage

pcrboot(object, type = c("boot", "jack"), B = 100, njack = 1,
        plot = TRUE, do.eff = TRUE, conf = 0.95, verbose = TRUE, ...)

Arguments

object	an object of class 'pcrfit'.
type	either bootstrapping or jackknifing.
B	numeric. The number of iterations.
njack	numeric. In case of type = "jack", how many datapoints to exclude. Defaults to leave-one-out.
plot	should the fitting and final results be displayed as a plot?
do.eff	logical. If TRUE, efficiency analysis will be performed.
conf	the confidence level.
verbose	logical. If TRUE, the iterations will be printed on the console.
...	other parameters to be passed on to the plotting functions.

Details

Non-parametric bootstrapping is applied using the centered residuals.
1) Obtain the residuals from the fit:

$\hat{\varepsilon}_t = y_t - f(x_t, \hat{\theta})$

2) Draw bootstrap pseudodata:

$y_{t}^{\ast} = f(x_t, \hat{\theta}) + \epsilon_{t}^{\ast}$

where $\epsilon_{t}^{\ast}$ are i.i.d. from distribution $\hat{F}$, where the residuals from the original fit are centered at zero.
3) Fit $\hat\theta^\ast$ by nonlinear least-squares.
4) Repeat B times, yielding bootstrap replications

$\hat\theta^{\ast 1}, \hat\theta^{\ast 2}, \ldots, \hat\theta^{\ast B}$

One can then characterize the EDF and calculate confidence intervals for each parameter:

$\theta \in [EDF^{-1}(\alpha/2), EDF^{-1}(1-\alpha/2)]$

The jackknife alternative is to perform the bootstrap on the data-predictor vector, i.e. eliminating a certain number of datapoints.
If the residuals are correlated or have non-constant variance the latter is recommended. This may be the case in qPCR data, as the variance in the low fluorescence region (ground phase) is usually much higher than in the rest of the curve.

Value

A list containing the following items:

ITER	a list containing each of the results from the iterations.
CONF	a list containing the confidence intervals for each item in ITER.

Each item contains subitems for the coefficients (coef), root-mean-squared error (rmse), residual sum-of-squares (rss), goodness-of-fit measures (gof) and the efficiency analysis (eff). If plot = TRUE, all data is plotted as boxplots including confidence intervals.

Author(s)

Andrej-Nikolai Spiess

References

Nonlinear regression analysis and its applications.
Bates DM & Watts DG.
Wiley, Chichester, UK, 1988.

Nonlinear regression.
Seber GAF & Wild CJ.
Wiley, New York, 1989.

Boostrap accuracy for non-linear regression models.
Roy T.
J Chemometics (1994), 8: 37-44.

Examples

## Simple bootstrapping with
## too less iterations...
par(ask = FALSE)
m1 <- pcrfit(reps, 1, 2, l4)
pcrboot(m1, B = 20)

## Jackknifing with leaving
## 5 datapoints out.
m2 <- pcrfit(reps, 1, 2, l4)
pcrboot(m2, type = "jack", njack = 5, B = 20)

---

Workhorse function for qPCR model fitting

Description

This is the workhorse function of the qpcR package that fits one of the available models to qPCR data using (weighted) nonlinear least-squares (Levenberg-Marquardt) fitting from nlsLM of the 'minpack.lm' package.

Usage

pcrfit(data, cyc = 1, fluo, model = l4, start = NULL, 
       offset = 0, weights = NULL, verbose = TRUE, ...)

Arguments

data	the name of the dataframe containing the qPCR runs.
cyc	the column containing the cycle data. Defaults to 1.
fluo	the column(s) containing the raw fluorescence data of the run(s). If more than one column is given, the model will be built with the replicates. See 'Details' and 'Examples'.
model	the model to be used for the analysis. Defaults to 'l4'.
start	a vector of starting values that can be supplied externally.
offset	an offset cycle number from the second derivative cut-off cycle for all MAK methods. See 'Details' and 'Example'.
weights	a vector with same length as data containing possible weights for the nonlinear fit, or an expression to calculate weights from. See 'Details'.
verbose	logical. If TRUE, fitting and convergence results will be displayed in the console.
...	other parameters to be passed to nlsLM.

Details

This is a newer (from qpcR 1.3-7 upwards) version of pcrfit. It is a much simpler implementation containing only the LM-Algorithm for fitting, but this fitting routine has proven to be so robust that other optimization routines (such as in optim) could safely be removed. The fitting is done with the new nlsLM function of the 'minpack.lm' package, which gives a model of class 'nls' as output.

This function is to be used at the single run level or on replicates (by giving several columns). The latter will build a single model based on the replicate values. If many models should be built on a cohort of replicates, use modlist and replist.

The offset value defines the offset cycles from the second derivative maximum that is used as a cut-off criterion in the MAK methods. See 'Examples'.

Since version 1.3-7, an expression given as a character string can be supplied to the weights argument. This expression, which is transferred to qpcR:::wfct, defines how the vector of weights is calculated from the data. In principle, five different parameters can be used to define weights:
"x" relates to the cycles $x_i$,
"y" relates to the raw fluorescence values $y_i$,
"error" relates to the error $\sigma(y_{i, j})$ of replicates $j$,
"fitted" relates to the fitted values $\hat{y}_i$ of the fit,
"resid" relates to the residuals $y_i - \hat{y}_i$ of the fit.
For "fitted" and "resid", the model is fit unweighted by pcrfit, the fitted/residual values extracted and these subsequently used for refitting the model with weights. These parameters can be used solely or combined to create a weights vector for different regimes. The most commonly used are (see also 'Examples'):
Inverse of response (raw fluorescence) $\frac{1}{y_i}$: "1/y"
Square root of predictor (Cycles) $\sqrt{x_i}$: "sqrt(x)"
Inverse square of fitted values: $\frac{1}{\hat{y}^2_i}$: "1/fitted^2"
Inverse variance $\frac{1}{\sigma^2(y_{i, j})}$: "1/error^2"

Value

A model of class 'nls' and 'pcrfit' with the following items attached:

DATA	the initial data used for fitting.
MODEL	the model used for fitting.
call2	the call to pcrfit.
parMat	the trace of the parameter values. Can be used to track problems.
opt.method	defaults to "LM".

Author(s)

Andrej-Nikolai Spiess

References

Bioassay analysis using R.
Ritz C & Streibig JC.
J Stat Soft (2005), 12: 1-22.

A Method for the Solution of Certain Problems in Least Squares.
K. Levenberg.
Quart Appl Math (1944), 2: 164-168.

An Algorithm for Least-Squares Estimation of Nonlinear Parameters.
D. Marquardt.
SIAM J Appl Math (1963), 11: 431-441.

Examples

## Simple l4 fit of F1.1 of the 'reps' dataset.
m1 <- pcrfit(reps, 1, 2, l4) 
plot(m1)

## Supply own starting values.
pcrfit(reps, 1, 2, l4, start = c(-5, -0.05, 11, 16)) 

## Make a replicate model,
## use inverse variance as weights.
m2 <- pcrfit(reps, 1, 2:5, l5, weights = "1/error^2")
plot(m2)

## Fit a mechanistic 'mak2' model
## to -1 cycle from SDM.
m3 <- pcrfit(reps, 1, 2, mak2, offset = -1)
plot(m3)

---

Summarize measures for the goodness-of-fit

Description

Calculates all implemented measures for the goodness-of-fit and returns them as a list. Works for objects of class pcrfit, lm, glm, nls, drc and many others...

Usage

pcrGOF(object, PRESS = FALSE)

Arguments

object	a fitted object.
PRESS	logical. If TRUE, the more calculation intensive $P^2$ is also returned.

Value

A list with all implemented Information criteria (AIC, AICc, BIC), residual variance, root-mean-squared-error, the reduced $\chi_{\nu}^2$ from fitchisq (if replicates) and the PRESS $P^2$ value (if PRESS = TRUE).

Author(s)

Andrej-Nikolai Spiess

Examples

## Single fit without replicates
## including PRESS statistic.
m1 <- pcrfit(reps, 1, 2, l5)
pcrGOF(m1, PRESS = TRUE)

## Fit containing replicates:
## calculation of reduced 
## chi-square included!
m2 <- pcrfit(reps, 1, 2:5, l5)
pcrGOF(m2)

---

Advanced qPCR data import function

Description

Advanced function to easily import/preformat qPCR data from delimited text files, the clipboard or the workspace. The data files can be located in a directory which is automatically browsed for all files. In a series of steps, the data can be imported and transformed to the appropriate format of the 'qpcR' package (such as in dataset reps, with 'Cycles' in the first column and named runs with raw fluorescence data in remaining columns). A dataset can function as a transformation template, and the remaining files in the directory are then formatted according to the established parameters. See 'Details' and tutorial video in http://www.dr-spiess.de/qpcR/tutorials.html.

Usage

pcrimport(file = NA, sep = NA, dec = NA, delCol = NA, delRow = NA,
          format = c(NA, "col", "row"), sampleDat = NA, refDat = NA,
          names = NA, sampleLen = NA, refLen = NA, check = TRUE,
          usePars = TRUE, dirPars = NULL, needFirst = TRUE, ...)

Arguments

file	either a directory such as "c:/temp" containing the data file(s), the Windows "clipboard" or an object in the workspace such as "reps".
sep	the field separator character, i.e. "\t" for tabs.
dec	the decimal seperator, i.e. ".".
delCol	unneeded columns to delete after successful import, i.e. 2, 1:3, seq(1, 5, by = 2), etc....
delRow	unneeded rows to delete after successful import, i.e. 2, 1:3, seq(1, 5, by = 2), etc....
format	how the data is organized, i.e. in columns or rows.
sampleDat	the columns with the raw fluorescence reporter dye data.
refDat	optional columns with the raw fluorescence reference dye data.
names	the row(s) that should be used for naming the runs.
sampleLen	the rows with the reporter dye cycles.
refLen	the rows with the reference dye cycles.
check	logical. If TRUE, a window displaying the transformed data after each step is displayed. This assists in choosing the right parameters.
usePars	logical. If TRUE, then all files in the directory are batch analysed using the stored parameters. See 'Details'.
dirPars	an optional directory such as "c:/pars" where the formatting parameters can be stored. If NULL, the 'qpcR' directory is used.
needFirst	logical. If TRUE, then the (alphabetically) first file in the directory is used for an initial definition of transformation parameters.
...	other parameters to be passed to read.delim.

Details

This function has been designed to offer maximal flexibility in importing qPCR data from all kinds of systems. This is accomplished by asking the user for many formatting options in single steps, with the final goal of obtaining a dataset that is transformed in a way suitable for pcrfit, as in all datasets in this package (i.e. 'reps'): it must be a dataframe with the first column containing the cycle numbers ("Cycles") and all subsequent columns with sensible sample names, such as "S1_1". In detail, the following steps are queried:

1) Location of the file. Either a directory containing the file(s), the (Windows) clipboard or a dataframe in the workspace.
2) How are the fields separated, i.e. by tabs?
3) What is the decimal separator?
4) Which columns can be deleted? For analysis, we only need the raw fluorescence values and sample names. Everything else should be deleted.
5) Which rows can be deleted? Same as above.
6) Are the runs organized in rows or in columns?
After these steps, the unwanted rows/columns are deleted and the data transformed into vertical format (if it was in rows).
7) In which columns are the runs with reporter dye data (i.e. SybrGreen)?
8) If a reference dye (i.e. ROX) was used, in which columns are the runs?
9) How should the runs be named (automatically or from a row/rows containing names)? If more than one row is supplied, the names in the rows are pasted together, i.e. "A4.GAPDH".
10) Which are the rows containing the raw fluorescence data from cycle to cycle for the reporter dye?
11) If a reference dye was used, which are the rows with the cycle to cycle data?
After these steps, a 'Cycles' column is prepended to the data which should then be in the right format for downstream analysis.
ATTENTION: Because of this step, if the imported data also initially had a column containing cycle numbers, these should be removed in steps 2) or 3)!

One major advantage of this function is that the formatting parameters are stored in a file and can be reused with new data, most conveniently when doing a batch analysis of several files in a directory. When needFirst = TRUE, the alphabetically first run in the directory is used for defining the formatting parameters, and if usePars = TRUE these are applied on all remaining datasets. If the initial definition of formatting parameters is not needed, then setting needFirst = FALSE will apply the last stored parameters on all datasets. By using different dirPars, one can establish different formatting options for different qPCR systems.

The function will query (if needFirst = TRUE) all parameters that are defined as NA. For example, using pcrimport(file = "c:/temp", sep = "\t", dec = ".", delCol = c(1, 3), ...) will result in these parameters not being queried.

If reference dye data was supplied, the function checks if the data is of same dimensions than the reporter dye data. The output is then the normalized fuorescence data $\frac{F_{rep}}{F_{ref}}$.

The 'Examples' feature internal datasets, but this function is best understood by the tutorial under http://www.dr-spiess.de/qpcR/tutorials.html.

Value

A list with the transformed data as data.frame list items, suitable for downstream analysis.

Author(s)

Andrej-Nikolai Spiess

Examples

## Not run: 
## EXAMPLE 1:
## Internal dataset format01.txt (in 'add01' directory)
## with 384 runs.
## Tab delimited, 30 cycles, only reporter dye,
## data in rows, and some unneeded columns and rows.
## This is the example data path, but could be any path
## with data such as c:/temp.
PATH <- path.package("qpcR")
PATHall <- paste(PATH, "/add01/", sep = "")
res <- pcrimport(PATHall)

## Answer queries with the following parameters and
## verify the effects in the 'View' windows:
## 1 => data is tab delimited
## 1 => decimal separator is "."
## c(1, 3) => remove columns 1 + 3
## 1:2 => remove rows 1 + 2
## 2 => data is in rows
## 0 => all data is from reporter dye
## 1 => sample names are in row #1
## 0 => reporter data goes until end of table
## Data is stored as dataframe list items and can
## then be analyzed:
ml <- modlist(res[[1]], model = l5)
plot(ml, which = "single")

## Alternative without query:
res <- pcrimport(PATHall, sep = "\t", dec = ".",
                 delCol = c(1, 3), delRow = 1:2,
                 format = "row", sampleDat = 0,
                 names = 1, sampleLen = 0,
                 check = FALSE)

## Do something useful with the data:
ml <- modlist(res[[1]], model = l5)
plot(ml, which = "single")
                 
## EXAMPLE 2:
## Internal datasets format02a.txt - format02d.txt
## (in 'add02' directory) with 96 runs.
## Tab delimited, 40 cycles, reporter dye + reference dye,
## data in columns, and some unneeded columns and rows.
PATH <- path.package("qpcR")
PATHall <- paste(PATH, "/add02/", sep = "")
res <- pcrimport(PATHall)

## Answer queries with the following parameters and
## verify the effects in the 'View' windows:
## 1 => data is tab delimited
## 1 => decimal separator is "."
## 1 => remove column 1 with cycle data
## c(1, 43, 44) => remove rows 1, 43, 44
## 1 => data is in columns
## 1:96 => data columns for reporter dye
## -2 => reference dye stacked under reporter dye
## 1 => sample names are in row #1
## 1:40 => reporter data is in rows 1-40
## -1 => reference data is stacked under samples
## Data is stored as dataframe list items and can
## then be analyzed.

## Alternative without query:
res2 <- pcrimport(PATHall, sep = "\t", dec = ".",
                 delCol = 1, delRow = c(1, 43, 44),
                 format = "col", sampleDat = 1:96,
                 refDat = -2, names = 1,
                 sampleLen = 1:40, refLen = -1,
                 check = FALSE)

## Do something useful with the data:
ml2 <- modlist(res2[[1]], model = l5)
plot(ml2)

## End(Not run)

---