R/efficiency.r

In rtpcr: qPCR Data Analysis

#' @title Amplification efficiency statistics and standard curves
#'
#' @description
#' The \code{efficiency} function calculates amplification efficiency (E)
#' and related statistics, including slope and coefficient of determination
#' (R\eqn{^2}), and generates standard curves for qPCR assays.
#'
#' @details
#' Amplification efficiency is estimated from standard curves generated by
#' regressing Ct values against the logarithm of template dilution.
#' For each gene, the function reports the slope of the standard curve,
#' amplification efficiency (E), and R\eqn{^2} as a measure of goodness of fit.
#' The function also provides graphical visualization of the standard curves.
#'
#' @author Ghader Mirzaghaderi
#'
#' @export
#'
#' @import tidyr
#' @import dplyr
#' @import reshape2
#' @import ggplot2
#' @import purrr
#'
#' @param df
#' A data frame containing dilution series and corresponding Ct values.
#' The first column should represent dilution levels, and the remaining
#' columns should contain Ct values for different genes.
#'
#' @return
#' A list with the following components:
#' \describe{
#'   \item{efficiency}{A data frame containing slope, amplification efficiency (E),
#'   and R\eqn{^2} statistics for each gene.}
#'   \item{Slope_compare}{A table comparing slopes between genes.}
#'   \item{plot}{A \code{ggplot2} object showing standard curves for all genes.}
#' }
#'
#' @examples
#'
#' # Load example efficiency data
#' data_efficiency
#'
#' # Calculate amplification efficiency and generate standard curves
#' efficiency(data_efficiency)




efficiency <- function(df) {

  
  # renaming the first column
  colnames(df)[1] <- "dilutions"
  dilutions <- df$dilutions
  # Fit the linear regressions and extract the slope and R-squared
  results <- df %>%
    select(-dilutions) %>%
    map_df(~{
      model <- lm(. ~ log10(dilutions))
      Slope <- coef(model)[2]
      R2 <- summary(model)$r.squared
      E <- 10^(-1/coef(model)[2])
      data.frame(Slope, R2, E)
    })
  
  # Add the column names to the results data frame
  results <- cbind(gene = colnames(df)[2:ncol(df)], results)
  colnames(results) <- c("Gene", "Slope", "R2", "E")
  rownames(results) <- NULL
  
  Ct <- df[,2]
  p <- ggplot(data = df) + 
    geom_point(aes(y = Ct, x = log10(dilutions))) +
    geom_smooth(data = df,aes(x = log10(dilutions), y = Ct), formula = y ~ x,
                method = "lm", se = F)

  

  
  # COMPAIRING SLOPES
  # making a long format data
  if(ncol(df) > 2){
    e <- melt(df, id = "dilutions")
    
    dilutions <- e$dilutions
    value <- e$value
    
    lm <- lm(value ~ log10(dilutions) * variable, data = e)
    slopes <- emtrends(lm, pairwise ~ variable, var = "log10(dilutions)")
    
    
    
    fits <- lapply(df[,-1], function(x) lm(x ~ log10(df[,1])))
    mdat <- melt(df, id = "dilutions")
    Ct <- mdat$value
    Gene <- mdat$variable
    variable <- mdat$Gene
    p <- ggplot(data = mdat) + 
      geom_point(aes(y = Ct, x = log10(dilutions), color = Gene)) +
      geom_smooth(data = mdat,aes(x = log10(dilutions), y = Ct, color = Gene), formula = y ~ x,
                  method = "lm", se = F)
  }
  


  if(ncol(df) == 2){
    res <- list(Efficiency = results, plot = p)
  } else {
    res <- list(Efficiency = results, Slope_compare = slopes, plot = p)
  }
   return(res)
}