#' Aberration calculation functions
#'
#' @param data Count or case data.
#' @param aberr_prefix Prefix of the aberrations in the data.
#' @param power Power of aberration.
#' @param X Sum of detected aberrations.
#' @param X2 Quadratic sum of detected aberrations.
#' @param N Number of cells analysed.
#' @param mean Mean.
#' @param var Variance.
#' @param assessment_u Expected \eqn{u}-value of the assessment. For a Poisson distribution this should be unity.
#' @param type Type of input data. Either "count" and "case".
#' @param aberr_module Aberration module.
#' @name calculate_aberr
NULL
# > NULL
#' @rdname calculate_aberr
calculate_aberr_power <- function(data, aberr_prefix = "C", power = 1) {
# Prepare data
aberr_data <- data %>%
.[grep(aberr_prefix, names(.))] %>%
t() %>%
as.data.frame()
powers <- aberr_data %>%
rownames() %>%
regmatches(., regexpr("[0-9]+", .)) %>%
as.numeric() %>%
`^`(power)
# Calculate aberration powers
aberr <- numeric(length = ncol(aberr_data))
for (i in seq_len(ncol(aberr_data))) {
aberr[i] <- sum(aberr_data[, i] * powers)
}
return(aberr)
}
#' @rdname calculate_aberr
calculate_aberr_mean <- function(X, N) {
mean <- X / N
return(mean)
}
#' @rdname calculate_aberr
calculate_aberr_var <- function(X, X2, N) {
var <- (X2 - X^2 / N) / (N - 1)
return(var)
}
#' @rdname calculate_aberr
calculate_aberr_disp_index <- function(mean, var) {
disp_index <- var / mean
return(disp_index)
}
#' @rdname calculate_aberr
calculate_aberr_u_value <- function(X, N, mean, var, assessment_u = 1) {
u_value <- (var / mean - assessment_u) * sqrt((N - 1) / (2 * (1 - 1 / X)))
return(u_value)
}
#' @rdname calculate_aberr
#' @importFrom rlang .data
init_aberr_table <- function(data, type = c("count", "case"), aberr_module = c("dicentrics", "translocations", "micronuclei")) {
# Validate parameters
type <- match.arg(type)
aberr_module <- match.arg(aberr_module)
if (type == "count") {
data <- data %>%
dplyr::mutate(
N = 0,
X = 0,
mean = 0,
var = 0,
DI = 0,
u = 0
) %>%
dplyr::select("D", "N", "X", dplyr::everything()) %>%
dplyr::mutate(
D = as.numeric(.data$D)
) %>%
dplyr::mutate(
dplyr::across(
.cols = c("N", "X", grep("C", names(.), value = TRUE)),
.fns = as.integer
)
)
} else if (type == "case") {
if (aberr_module %in% c("dicentrics", "micronuclei")) {
data <- data %>%
dplyr::mutate(
N = 0,
X = 0,
y = 0,
y_err = 0,
DI = 0,
u = 0
) %>%
dplyr::select("N", "X", dplyr::everything()) %>%
dplyr::mutate(
dplyr::across(
.cols = c("N", "X", grep("C", names(.), value = TRUE)),
.fns = as.integer
)
)
} else if (aberr_module == "translocations") {
data <- data %>%
dplyr::mutate(
N = 0,
X = 0,
Fp = 0,
Fp_err = 0,
DI = 0,
u = 0,
Xc = 0,
Fg = 0,
Fg_err = 0
) %>%
dplyr::select("N", "X", dplyr::everything()) %>%
dplyr::mutate(
dplyr::across(
.cols = c("N", "X", grep("C", names(.), value = TRUE)),
.fns = as.integer
)
)
}
}
return(data)
}
#' Calculate aberrations table
#'
#' @param data Count or case data.
#' @param type Type of input data. Either "count" and "case".
#' @param assessment_u Expected \eqn{u}-value of the assessment. For a Poisson distribution this should be unity.
#' @param aberr_module Aberration module, required for \code{type = "case"}.
#'
#' @return Data frame containing cell count (\eqn{N}), aberrations (\eqn{X}),
#' and other coefficients (dispersion index, \eqn{u}-value, ...), as well as
#' raw count or case \code{data}.
#' @export
#' @importFrom rlang .data
calculate_aberr_table <- function(data, type = c("count", "case"), aberr_module = c("dicentrics", "translocations", "micronuclei"), assessment_u = 1) {
# Validate parameters
type <- match.arg(type)
aberr_module <- match.arg(aberr_module)
if (type == "count") {
data <- data %>%
dplyr::as_tibble() %>%
dplyr::mutate(
N = as.integer(rowSums(.[grep("C", names(.))])),
X = calculate_aberr_power(., power = 1),
X2 = calculate_aberr_power(., power = 2),
mean = calculate_aberr_mean(.data$X, .data$N),
var = calculate_aberr_var(.data$X, .data$X2, .data$N),
DI = calculate_aberr_disp_index(.data$mean, .data$var),
u = calculate_aberr_u_value(.data$X, .data$N, .data$mean, .data$var, assessment_u = assessment_u)
) %>%
dplyr::mutate(
dplyr::across(
.cols = c(grep("C", names(.), value = TRUE)),
.fns = as.integer
)
) %>%
dplyr::select(-"X2") %>%
dplyr::select("D", "N", "X", dplyr::everything())
} else if (type == "case") {
data <- data %>%
dplyr::as_tibble() %>%
dplyr::mutate(
N = as.integer(rowSums(.[grep("C", names(.))])),
X = calculate_aberr_power(., power = 1),
X2 = calculate_aberr_power(., power = 2),
var = calculate_aberr_var(.data$X, .data$X2, .data$N),
mean = calculate_aberr_mean(.data$X, .data$N),
std_err = sqrt(.data$var / .data$N),
DI = calculate_aberr_disp_index(.data$mean, .data$var),
u = calculate_aberr_u_value(.data$X, .data$N, .data$mean, .data$var, assessment_u = assessment_u)
) %>%
dplyr::mutate(
dplyr::across(
.cols = c("N", "X", grep("C", names(.), value = TRUE)),
.fns = as.integer
)
) %>%
dplyr::select(-"X2", -"var") %>%
dplyr::select("N", "X", dplyr::everything())
# Rename mean and std_err columns
if (aberr_module %in% c("dicentrics", "micronuclei")) {
data <- data %>%
dplyr::select(
"N", "X", dplyr::matches("^C[0-9]+$"),
"y" = "mean", "y_err" = "std_err",
"DI", "u"
)
} else if (aberr_module == "translocations") {
data <- data %>%
dplyr::select(
"N", "X", dplyr::matches("^C[0-9]+$"),
"Fp" = "mean", "Fp_err" = "std_err",
"DI", "u"
)
}
}
return(data)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.