---
title: "Using RDML with other R packages"
author: "Konstantin A. Blagodatskikh, Micha&#322; Burdukiewicz, Stefan R&ouml;diger, Andrej-Nikolai Spiess"
date: "`r Sys.Date()`"
output: 
rmarkdown::html_vignette:
toc: true
pandoc_args: [
"--number-sections"
]
bibliography: "RDML.bib"
vignette: >
  %\VignetteIndexEntry{Using RDML with other R packages}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
```{r, echo=FALSE, message=FALSE}
knitr::opts_chunk$set(collapse=T, comment="#>")
```

# A Practical Example for Usage of the **RMDL** package

A key benefit of the **RMDL** package is that it enables further statistical 
analysis of RDML data. This section describes **R** packages 
that link with the **RMDL** package. It will provide examples on how 
other packages allow the manipulation of RDML objects.

In this example we will employ the **qpcR** package [@ritz_qpcr_2008] to calculate Cq values after the selection of an optimal sigmoidal model, as suggested by the Akaike's An Information Criterion. The Cq values will then be used to calculate the amplification efficiency from a calibration curve, based on the `effcalc` function of the **chipPCR** package [@roediger2015chippcr].

Note: The data used here serve as example only. Overall, the quality of the 
measurement is not appropriate for further application in a study. Here, we show 
that the estimated amplification efficiency from the StepOne system differs 
from that of the proposed pipeline (**RDML** $\rightarrow$ **qpcR** $\rightarrow$ **chipPCR**).

In principle, the following code examples can be combined with report generating 
toolkits for data analysis pipelines, such as **Nozzle** (see @gehlenborg_nozzle:_2013).

# Preparation of the data

In the **RDML** vignette, it was demonstrated how to read-in RDML data. In this section we will 
continue with the built-in RDML example file `stepone_std.rdml`. This file was 
obtained during the measurement of human DNA concentration by a *LightCycler 
96* instrument (Roche) and the *XY-Detect* kit (Syntol, Russia). The file is opened as 
described before:

```{r, results="hide"}
# Load the RDML package and use its functions to `extract` the required data
library(RDML)
filename <- system.file("extdata/stepone_std.rdml", package="RDML")
raw_data <- RDML$new(filename=filename)
```

# Working with metadata of the example file

At this stage, all data from the `stepone_std.rdml` RDML file are available. For 
example, the amplification efficiency as provided by the StepOne™ Real-Time 
PCR System software can be fetched from the package via `raw_data$target[["RNase 
P"]]$amplificationEfficiency`. In another example, we extract the information 
about the target in this experiment.

```{r}
raw_data$target
```
To get some information about the run we can use:

```{r}
raw_data$experiment[["Standard Curve Example"]]$run
```

# Extraction of raw amplification curve data

For convenience, we use the pipe function `%>%` from the **magrittr** package for 
further analysis. In the following, we fetch the amplification curve data from the RDML file.

```{r}
raw_data_tab <- raw_data$AsTable(
    # Custom name pattern 'position~sample~sample.type~target~run.id'
    name.pattern=paste(
        react$position,
        react$sample$id,
        private$.sample[[react$sample$id]]$type$value,
        data$tar$id,
        run$id$id, # run id added to names
        sep="~"))
# Get all fluorescence data and assign them to the object fdata
fdata <- as.data.frame(raw_data$GetFData(raw_data_tab, long.table=FALSE))
```

# Plotting of raw amplification curve data

The plotting of the raw data is an important step to visually inspect the data. 
In this example, we use the **ggplot2** package [@wickham_ggplot2_2009] instead 
of the default **R** graphics functions. 

```{r}
# Load the ggplot2 package for plotting
library(ggplot2)
# Load the reshape2 package to rearrange the data
library(reshape2)
# Rearrange and plot the raw data
fdata_gg <- melt(fdata, id.vars="cyc")
ggplot(data=fdata_gg, aes(x=cyc, y=value, color=variable)) + 
    geom_line() + labs(x="Cycle", y="RFU") + theme_light() +
    theme(legend.position="top",
    legend.direction="horizontal")

```

# Analysis of the qPCR amplification curve data by a custom made function

During the next steps we will employ the **qpcR** package. The function 
`mselect` performs sigmoid model selection by different criteria (e.g., 
bias-corrected Akaike's Information Criterion). Note: In most cases a sigmoidal 
model with seven parameters was selected.

The function `efficiency` calculates the qPCR Cq values, amplification efficiency and other important qPCR parameters. In this example, we set the parameter type of the `efficiency` 
function to `Cy0`,  which will calculate and report the Cy0 value. According to 
@guescini_new_2008, the Cy0 value is the intersection of a tangent on the first 
derivative maximum with the abscissa. However, for all further analyses we will 
use the second derivative maximum cycle as Cq value. This 
supplementary material will not focus on the selection of a certain Cq method. For an 
objective decision, we would like to guide the reader to the study by 
@ruijter_evaluation_2013.


```{r, echo=TRUE, include=TRUE}
# Use the magrittr package to create pipes
library(magrittr)

# Write a custom function that calculates the Cq values and other curve parameters
library(qpcR)
res_fit <- do.call(cbind, lapply(2L:ncol(fdata), function(block_column) {
    
    res_model <- try(mselect(pcrfit(data=fdata, cyc=1, fluo=block_column), verbose=FALSE, do.all=TRUE), silent=TRUE)
    if(res_model %>% class=="try-error") {
        res_model <- NA
    }
    else{
        try(efficiency(res_model, plot=FALSE, type="Cy0"), silent=TRUE)
    }
        }
    )
)
# Assign column names
colnames(res_fit) <- colnames(fdata)[-1]

# Fetch only the Cq values (second derivative maximum) and combine them in a
# data.frame

Cq_SDM <- res_fit[rownames(res_fit)=="cpD2", ] %>% unlist %>% as.data.frame
colnames(Cq_SDM) <- c("Cq")

# Prepare the dilutions and calculated Cq values for further usage in the effcalc
# function from the chipPCR package

dilution <- c(as.factor("ntc"), as.factor("unk"), 10000, 5000, 2500, 1250, 625)
Cq_values <- matrix(Cq_SDM[, "Cq"], nrow=length(dilution), ncol=3, byrow=TRUE)
```

# Inspection and interpretation of the analyzed data

Below, we arbitarily selected a non-template control (negative) and an 
unknown sample (positive) (two out of `r ncol(fdata)-1` amplification curves) 
for the presentation of the coefficients that were calculated from the raw amplication 
curve data.

```{r, results="asis", echo=FALSE}
knitr::kable(res_fit[, c(1, 5)], caption="Selected data of calculated values 
from the fitted models. The values reported in the table were calculated by a 
chain of functions from the **qpcR** package. In detail, this were the *mselect* 
function to calculate the optimal multiparametric sigmoid model and the 
*efficiency* function to calculate the Cq values and other curve parameters. 
More information about these functions is available from the manual of the 
**qpcR** package. The samples `A01` (ntc, non-template cntrol) and 
`A05` (unk, unknown) were arbitrarily selected for 
presentation. The custom made function assigned a NA (not available) to sample 
`A01` because no Cq value or other curve 
parameters could be calculated from the non-template control.")
```

# Calculation of the amplication efficiency from calibration data

We now use our calculated Cq values from the dilutions steps (`r 625*c(1, 2, 4, 8, 
16)`) to determine the efficiency from the calibration curve. According to the StepOne software, the amplification efficiency is approx. 
`r round(raw_data$target[["RNase P"]]$amplificationEfficiency, 1)`%. 
To confirm these results, the `effcalc` function was employed to determine the coefficients of the 
calibration curve.

The following few lines are needed to calculate the amplication efficiency based on the previously calculated Cq values.

```{r, echo=TRUE}
# Load the chipPCR package
library(chipPCR)

# Use the effcalc function from the chipPCR package to calculate the amplification
# efficiency and store the results in the object res_efficiency

res_efficiency <- effcalc(dilution[-c(1:3)], Cq_values[-c(1:3), ], logx=TRUE)

# Use the %>% function from the magrittr package to plot the results (res_efficiency) 
# from the effcalc function 

res_efficiency %>% plot(., CI=TRUE, ylab="Cq (SDM)", 
                        main="Second Derivative Maximum Method")
```

The amplification efficiency estimated with the customized function was `r res_efficiency@amplification.efficiency`%, 
which is comparable to the value reported in the `stepone_std.rdml` file. 

```{r, results="asis"}
# Combine the sample labels and the Cq values as calculate by the Second 
# Derivative Maximum Method (cpD2).

sample_Cq <- data.frame(sample=c("ntc", "unk", 
                                 10000, 5000, 2500, 1250, 625), 
                        Cq_values)

# Print table of all Cq values
# Use the kable function from the knitr package to print a table

knitr::kable(sample_Cq, caption="Cq values as calculate by the Second Derivative 
Maximum Method (cpD2). ntc, non-template control. unk, unknown sample. X1, X2 
and X3 are the Cq values from a triplicate measurement.")
```

```{r, results="asis"} 
knitr::kable(res_efficiency, caption="Analysis of the amplification efficiency. 
The table reports the concentration-depentent average Cq values from three 
replicates per dilution step. In addition, the standard deviation (SD) and the 
Coefficient of Variation (RSD [%]) are presented. The results indicate that the 
data basis for the calibration curve is valid.")
```

The Cq values (`r sample_Cq[2, -1]`) from the unknown sample `unk` had an 
average Cq of `r paste(round(mean(as.numeric(sample_Cq[2, -1])), 2), "±", 
round(sd(sample_Cq[2, -1]), 2))`.

# References