# Load the RDML package and use its functions to `extract` the required data
library(RDML7)
filename <- system.file("extdata/stepone_std.rdml", package="RDML7")
raw_data <- readRDML(filename)

raw_data$target[["RNase P"]]$amplificationEfficiency
#> [1] 93.91181

names(raw_data$experiment[["Standard Curve Example"]]$run$Run001)
#>  [1] "id"                            "description"                  
#>  [3] "documentation"                 "experimenter"                 
#>  [5] "instrument"                    "dataCollectionSoftware"       
#>  [7] "backgroundDeterminationMethod" "cqDetectionMethod"            
#>  [9] "thermalCyclingConditions"      "pcrFormat"                    
#> [11] "runDate"                       "react"
raw_data$experiment[["Standard Curve Example"]]$run$Run001$instrument
#> [1] "Applied Biosystems StepOne™ Instrument"
raw_data$experiment[["Standard Curve Example"]]$run$Run001$pcrFormat
#> <RDML7::pcrFormatType>
#>  @ rows       : int 8
#>  @ columns    : int 12
#>  @ rowLabel   : <RDML7::labelFormatType>
#>  .. @ value   : chr "ABC"
#>  .. @ variants: chr [1:3] "ABC" "123" "A1a1"
#>  @ columnLabel: <RDML7::labelFormatType>
#>  .. @ value   : chr "123"
#>  .. @ variants: chr [1:3] "ABC" "123" "A1a1"

fdata <- 
  asTable(raw_data,
          # Custom name pattern
          namePattern = "{position}:{sample}~{sampleType}~{target}") |> 
  # Get all fluorescence data
  getFData(raw_data, request = _, long.table=FALSE)

# 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")

# 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(class(res_model)[1] == "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"), as.factor("unk"),
              10000, 5000, 2500, 1250, 625)
Cq_values <- matrix(Cq_SDM[, "Cq"], nrow=length(dilution), ncol=3, byrow=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")

# 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", "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.")

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.")