genorm calculations

DESCRIPTION

@@ -16,14 +16,15 @@ BugReports: https://github.com/MahShaaban/pcr/issues
 Depends: R (>= 3.4.0)
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 6.0.1
+RoxygenNote: 6.1.1
 Imports: devtools,
   magrittr,
   dplyr,
   tidyr,
   purrr,
   readr,
-  ggplot2
+  ggplot2,
+  FSA
 Suggests: testthat,
     knitr,
     rmarkdown,

NAMESPACE

@@ -1,16 +1,31 @@
 # Generated by roxygen2: do not edit by hand
 
+export(.pcr_amount)
+export(.pcr_average)
+export(.pcr_calibrate)
+export(.pcr_cv)
+export(.pcr_error)
+export(.pcr_normalize)
+export(.pcr_plot_analyze)
+export(.pcr_plot_assess)
+export(.pcr_sd)
+export(.pcr_trend)
 export(pcr_analyze)
 export(pcr_assess)
 export(pcr_curve)
 export(pcr_dct)
 export(pcr_ddct)
 export(pcr_efficiency)
+export(pcr_genorm)
 export(pcr_lm)
+export(pcr_m)
+export(pcr_nf)
+export(pcr_stability)
 export(pcr_standard)
 export(pcr_test)
 export(pcr_ttest)
 export(pcr_wilcox)
+importFrom(FSA,geomean)
 importFrom(devtools,install_github)
 importFrom(dplyr,bind_cols)
 importFrom(dplyr,bind_rows)
@@ -36,6 +51,7 @@ importFrom(ggplot2,ggplot)
 importFrom(magrittr,"%>%")
 importFrom(purrr,map)
 importFrom(readr,read_csv)
+importFrom(stats,aggregate)
 importFrom(stats,coefficients)
 importFrom(stats,confint)
 importFrom(stats,cor)

R/analyses_fun.R

@@ -381,19 +381,20 @@ pcr_curve <- function(df, group_var, reference_gene, reference_group,
 #'
 #' @inheritParams pcr_ddct
 #' @inheritParams pcr_curve
-#' @param method A character string; 'delta_delta_ct' default, 'delta_ct' or
-#' 'relative_curve' for invoking a certain analysis model
+#' @param method A character string; 'delta_delta_ct' default, 'delta_ct',
+#' 'relative_curve' or 'genorm' for invoking a certain analysis model.
 #' @param ... Arguments passed to the methods
 #'
 #' @return A data.frame by default, when \code{plot} is TRUE returns a plot.
-#' For details; \link{pcr_ddct}, \link{pcr_dct} and \link{pcr_curve}.
+#' For details; \link{pcr_ddct}, \link{pcr_dct},  \link{pcr_curve} and
+#' \link{pcr_genorm}.
 #'
 #' @details The different analysis methods can be invoked using the
-#' argument method with 'delta_delta_ct' default, 'delta_ct' or
-#' 'relative_curve' for the double delta \eqn{C_T}, delta ct or the standard curve
-#' model respectively. Alternatively, the same methods can be applied by using
-#' the corresponding functions directly: \link{pcr_ddct}, \link{pcr_dct} or
-#' \link{pcr_curve}
+#' argument method with 'delta_delta_ct' default, 'delta_ct', 'relative_curve'
+#' or 'genorm' for the double delta \eqn{C_T}, delta \eqn{C_T}, the standard
+#' curve or the geNorm model respectively. Alternatively, the same methods can
+#' be applied by using the corresponding functions directly: \link{pcr_ddct},
+#' \link{pcr_dct}, \link{pcr_curve} or \link{pcr_genorm}.
 #'
 #' @references Livak, Kenneth J, and Thomas D Schmittgen. 2001. “Analysis of
 #' Relative Gene Expression Data Using Real-Time Quantitative PCR and the
@@ -476,11 +477,141 @@ pcr_curve <- function(df, group_var, reference_gene, reference_group,
 #'            method = 'relative_curve',
 #'            plot = TRUE)
 #'
+#' # applying geNorm method
+#' ## locate and read raw ct data
+#' fl <- system.file('extdata', 'ct5.csv', package = 'pcr')
+#' ct5 <- readr::read_csv(fl)
+#'
+#' # add grouping variable
+#' group_var <- rep(c('control', 'treatment'), each = 2)
+#'
+#' # calculate relative expression change
+#' pcr_genorm(ct5,
+#'            group_var = group_var,
+#'            reference_group = 'control',
+#'            reference_gene = paste0('REF', 1:3))
+#'
 #' @export
 pcr_analyze <- function(df, method = 'delta_delta_ct', ...) {
   switch(method,
          'delta_delta_ct' = pcr_ddct(df, ...),
          'delta_ct' = pcr_dct(df, ...),
-         'relative_curve' = pcr_curve(df, ...))
+         'relative_curve' = pcr_curve(df, ...),
+         'genorm' = pcr_genorm(df, ...))
 }
 
+#' Calculate the geNorm method
+#'
+#' Applies the geNorm method by first calculating a normalization facor from
+#' multiple reference genes and then apply the \link{pcr_dct} using the factor.
+#'
+#' @inheritParams pcr_ddct
+#'
+#' @return A data.frame of 5 columns:
+#' \itemize{
+#'   \item group The unique entries in group_var
+#'   \item gene The column names of df. reference_gene is dropped
+#'   \item factor The name of the normalization factor 'NF'
+#'   \item relative_expression The expression of target genes normalized by
+#'   a reference_gene and calibrated by a reference_group
+#'   \item error The standard deviation of the relative_expression
+#' }
+#'
+#' @examples
+#' ## locate and read raw ct data
+#' fl <- system.file('extdata', 'ct5.csv', package = 'pcr')
+#' ct5 <- readr::read_csv(fl)
+#'
+#' # add grouping variable
+#' group_var <- rep(c('control', 'treatment'), each = 2)
+#'
+#' # calculate relative expression change
+#' pcr_genorm(ct5,
+#'            group_var = group_var,
+#'            reference_group = 'control',
+#'            reference_gene = paste0('REF', 1:3))
+#'
+#' @importFrom magrittr %>%
+#' @importFrom dplyr select mutate group_by summarise
+#'
+#' @export
+pcr_genorm <- function(df, group_var, reference_group, reference_gene) {
+  # caluclate a normalization factor
+  nf <- pcr_nf(df,
+               group_var = group_var,
+               reference_gene = reference_gene,
+               reference_group = reference_group)
+
+  # calculate dct for genes of interest and
+  # transform error term
+  goi <- select(df, -reference_gene) %>%
+    pcr_dct(group_var = group_var, reference_group = reference_group) %>%
+    mutate(error = fold_change*error*log(2)) %>%
+    select(group, gene, fold_change, error)
+
+  # normalize by normalization factors
+  res <- full_join(goi, nf) %>%
+    group_by(group, gene, factor) %>%
+    summarise(relative_expression = fold_change/mean,
+              error = relative_expression * sum((sd/mean)^2 + (error/fold_change)^2)^.5 )
+
+  # return resulsts
+  return(res)
+}
+
+#' Calculate normalization factors
+#'
+#' Calculates a normalization factor from multiple reference genes.
+#'
+#' @inheritParams pcr_ddct
+#'
+#' @return A data.frame of 5 columns
+#' \itemize{
+#'   \item group The unique entries in group_var
+#'   \item gene The column names of df. reference_gene is dropped
+#'   \item mean The average value of the normalization factor
+#'   \item sd The standard deviation of the normalization factor
+#'   \item factor The name of the normalization factor 'NF'
+#' }
+#'
+#' @examples
+#' ## locate and read raw ct data
+#' fl <- system.file('extdata', 'ct5.csv', package = 'pcr')
+#' ct5 <- readr::read_csv(fl)
+#'
+#' # add grouping variable
+#' group_var <- rep(c('control', 'treatment'), each = 2)
+#'
+#' # calculate a normalization factor
+#' pcr_nf(ct5,
+#'        group_var = group_var,
+#'        reference_group = 'control',
+#'        reference_gene = paste0('REF', 1:3))
+#'
+#' @importFrom magrittr %>%
+#' @importFrom dplyr mutate select group_by summarise
+#' @importFrom FSA geomean
+#'
+#' @export
+pcr_nf <- function(df, group_var, reference_group, reference_gene) {
+  if(length(reference_gene) < 2) {
+    stop('reference gene must be a vector of two or more elements.')
+  }
+
+  # select reference gene ct values
+  df_ref <- select(df, reference_gene)
+
+  # calculate the dct
+  dct <- pcr_dct(df_ref, group_var = group_var, reference_group = reference_group)
+
+  # caluculate the mean and sd
+  nf <- dct %>%
+    mutate(error = fold_change*error*log(2)) %>%
+    group_by(group) %>%
+    summarise(mean = geomean(fold_change),
+              sd = mean * sum((error/ (fold_change*length(reference_gene)))^2)^.5) %>%
+    mutate(factor = 'NF')
+
+  # return normalization factor
+  return(nf)
+}

R/helper_fun.R

@@ -49,7 +49,8 @@
 #' @importFrom tidyr gather
 #'
 #' @keywords internal
-
+#'
+#' @export
 .pcr_average <- function(df, group_var, amount, tidy = FALSE) {
   # when group_var is provided
   # calculate the averages using group_var in group_by
@@ -127,7 +128,8 @@
 #' @importFrom dplyr select mutate_if starts_with
 #'
 #' @keywords internal
-
+#'
+#' @export
 .pcr_normalize <- function(df, reference_gene, mode = 'subtract',
                            tidy = FALSE) {
   # get the reference_gene column and unlist
@@ -197,7 +199,8 @@
 #' @importFrom tidyr gather
 #'
 #' @keywords internal
-
+#'
+#' @export
 .pcr_calibrate <- function(df, reference_group, mode = 'subtract',
                            tidy = FALSE) {
   # get the row index of the reference group
@@ -255,7 +258,8 @@
 #' @importFrom stats sd
 #'
 #' @keywords internal
-
+#'
+#' @export
 .pcr_sd <- function(df, group_var, tidy = FALSE) {
   # group_by the group_var
   # calculate standard deviation
@@ -308,7 +312,8 @@
 #' @importFrom tidyr gather
 #'
 #' @keywords internal
-
+#'
+#' @export
 .pcr_error <- function(df, reference_gene, tidy = FALSE) {
   # get the reference_gene column and unlist
   ref <- select(df, reference_gene) %>%
@@ -377,8 +382,8 @@
 #' @importFrom dplyr group_by summarise ungroup
 #'
 #' @keywords internal
-
-
+#'
+#' @export
 .pcr_cv <- function(amounts, group_var, tidy = FALSE) {
   # group_by group_var and calculate cv
   cv <- mutate(amounts, group = group_var) %>%
@@ -435,7 +440,8 @@
 #'            slope = slope)
 #'
 #' @keywords internal
-
+#'
+#' @export
 .pcr_amount <- function(df, intercept, slope) {
   # make a data.frame of intercept, slope ana gene names
   curve <- data_frame(intercept,
@@ -492,7 +498,8 @@
 #' .pcr_trend(ct3, amount = amount)
 #'
 #' @keywords internal
-
+#'
+#' @export
 .pcr_trend <- function(df, amount) {
   # make a trend line using linear regression
   trend_line <- map(df, function(x) {

R/package_data.R

@@ -70,3 +70,25 @@
 #'
 #' @source \url{https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1395339/}
 "ct4"
+
+#' \eqn{C_T} values from qPCR with multiple reference genes
+#'
+#' A dataset containing the \eqn{C_T} values of one gene of interest and three
+#' reference genes in two control and two treatment samples.
+#'
+#' @format A data.frame with 4 rows and 4 variables:
+#' \describe{
+#'   \item{GOT1}{\eqn{C_T} values of the gene of interest}
+#'   \item{REF1}{\eqn{C_T} values of the first reference gene}
+#'   \item{REF2}{\eqn{C_T} values of the second reference gene}
+#'   \item{REF3}{\eqn{C_T} values of the third reference gene}
+#' }
+"ct5"
+
+#' Expression values reference genes
+#'
+#' A dataset containing the relative expression values of ten reference genes
+#' in five different tissues.
+#'
+#' @format A \code{data.frame} with 85 rows and 11 variables
+"hk_tissue"

R/pcr.R

@@ -22,4 +22,5 @@ if(getRversion() >= "2.15.1")  utils::globalVariables(c('group',
                                                         'upper',
                                                         'lower',
                                                         'log_amount',
-                                                        'ct3'))
+                                                        'relative_expression',
+                                                        'fold_change'))

R/plotting_fun.R

@@ -66,6 +66,7 @@
 #' # make a plot
 #' .pcr_plot_analyze(df, method = 'relative_curve')
 #'
+#' @export
 .pcr_plot_analyze <- function(df, method, facets = FALSE) {
   y <- switch (method,
                'delta_delta_ct' = 'relative_expression',
@@ -125,6 +126,7 @@
 #' @importFrom tidyr gather
 #' @importFrom ggplot2 ggplot aes geom_point geom_errorbar facet_wrap geom_smooth
 #'
+#' @export
 .pcr_plot_assess <- function(df, amount, reference_gene, method) {
   switch (method,
     'efficiency' = {