Skip to content

Commit

Permalink
simplify differential_tests code
Browse files Browse the repository at this point in the history
  • Loading branch information
@AliciaSchep
AliciaSchep committed Oct 12, 2017
1 parent 8a5a34c commit 243d192
Show file tree
Hide file tree
Showing 3 changed files with 58 additions and 78 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: chromVAR
Type: Package
Title: Chromatin Variation Across Regions
Version: 0.99.0
Version: 0.99.0.9999
Date: 2016-12-20
Authors@R: c(
person("Alicia", "Schep", email = "[email protected]", role = c("aut","cre")),
Expand Down
1 change: 1 addition & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,7 @@ importFrom(graphics,text)
importFrom(plotly,ggplotly)
importFrom(plotly,plotlyOutput)
importFrom(plotly,renderPlotly)
importFrom(stats,aggregate)
importFrom(stats,anova)
importFrom(stats,approx)
importFrom(stats,as.dist)
Expand Down
133 changes: 56 additions & 77 deletions R/differential_tests.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,87 +10,78 @@
#' @return data.frame with p value and adjusted p value
#' @export
#' @author Alicia Schep
#' @importFrom stats aggregate
#' @examples
#' # Load very small example results from computeDeviations
#' data(mini_dev, package = "chromVAR")
#' difdev <- differentialDeviations(mini_dev, "Cell_Type")
differentialDeviations <- function(object,
groups,
alternative = c("two.sided", "less",
"greater"),
parametric = TRUE) {
groups,
alternative = c("two.sided", "less",
"greater"),
parametric = TRUE) {
stopifnot(is(object,"chromVARDeviations"))
if (length(groups) == 1 && groups %in% colnames(colData(object))){
groups <- colData(object)[groups]
} else if (length(groups) != ncol(object)){
if (length(groups) == 1 && groups %in% colnames(colData(object))) {
groups <- colData(object)[[groups]]
} else if (length(groups) != ncol(object)) {
stop("invalid groups input, must be vector of lench ncol(object) or column",
" name from colData(object)")
}

inputs <- lapply(split(seq_len(ncol(object)),groups),
function(x) deviations(object)[,x])
groups <- as.factor(groups)

alternative <- match.arg(alternative)
inputs <- deviations(object)

if (parametric) {
if (length(inputs) == 2) {
if (nlevels(groups) == 2) {
# t-test
p_val <- vapply(seq_len(nrow(inputs[[1]])),
function(x) t_helper(inputs[[1]][x,],
inputs[[2]][x, ],
alternative),
0)
p_val <- apply(inputs, 1, t_helper, groups, alternative)
} else {
# anova
reshaped <- lapply(seq_len(nrow(inputs[[1]])),
function(x) lapply(inputs, function(y) y[x,]))
p_val <- vapply(reshaped, function(x) do.call(anova_helper, x), 0)
p_val <- apply(inputs, 1, anova_helper, groups)
}
} else {
if (length(inputs) == 2) {
if (nlevels(groups) == 2) {
# wilcoxon
p_val <- vapply(seq_len(nrow(inputs[[1]])),
function(x) wilcoxon_helper(inputs[[1]][x, ],
inputs[[2]][x, ],
alternative),
0)
p_val <- apply(inputs, 1, wilcoxon_helper, groups, alternative)
} else {
# kruskal-wallis
reshaped <- lapply(seq_len(nrow(inputs[[1]])),
function(x) lapply(inputs, function(y) y[x,]))
p_val <- vapply(reshaped, function(x) do.call(kw_helper, x), 0)
p_val <- apply(inputs, 1, kw_helper, groups)
}
}

p_adj <- p.adjust(p_val, method = "BH")
return(data.frame(p_value = p_val, p_value_adjusted = p_adj))
}

kw_helper <- function(...) {
inputs <- list(...)
vals <- do.call(c, inputs)
group <- factor(do.call(c, lapply(seq_along(inputs),
function(x) rep(x, length(inputs[[x]])))))
res <- kruskal.test(vals ~ group)
return(res$p.value)
}

wilcoxon_helper <- function(x, y, alternative) {
return(wilcox.test(x, y, alternative = alternative, paired = FALSE)$p.value)
t_helper <- function(x, groups, alternative) {
splitx <- split(x, groups)
return(t.test(splitx[[1]],splitx[[2]],
alternative = alternative,
paired = FALSE,
var.equal = FALSE)$p.value)
}


t_helper <- function(x, y, alternative) {
return(t.test(x, y, alternative = alternative, paired = FALSE,
var.equal = FALSE)$p.value)
anova_helper <- function(x, groups) {
tmpdf <- data.frame(groups = groups, devs = x)
res <- oneway.test(devs ~ groups, tmpdf, var.equal = FALSE)
return(res$p.value)
}

anova_helper <- function(...) {
inputs <- list(...)
vals <- do.call(c, inputs)
group <- factor(do.call(c, lapply(seq_along(inputs),
function(x) rep(x, length(inputs[[x]])))))
anova_res <- oneway.test(vals ~ group, var.equal = FALSE)
return(anova_res$p.value)
kw_helper <- function(x, groups) {
tmpdf <- data.frame(groups = groups, devs = x)
res <- kruskal.test(devs ~ groups, tmpdf)
return(res$p.value)
}

wilcoxon_helper <- function(x, groups, alternative) {
splitx <- split(x, groups)
return(wilcox.test(splitx[[1]], splitx[[2]],
alternative = alternative,
paired = FALSE)$p.value)
}

#' differentialVariability
#'
Expand All @@ -108,45 +99,33 @@ anova_helper <- function(...) {
#' difvar <- differentialVariability(mini_dev, "Cell_Type")
differentialVariability <- function(object, groups, parametric = TRUE) {
stopifnot(is(object,"chromVARDeviations"))
if (length(groups) == 1 && groups %in% colnames(colData(object))){
groups <- colData(object)[groups]
} else if (length(groups) != ncol(object)){
if (length(groups) == 1 && groups %in% colnames(colData(object))) {
groups <- colData(object)[[groups]]
} else if (length(groups) != ncol(object)) {
stop("invalid groups input, must be vector of lench ncol(object) or column",
" name from colData(object)")
}

groups <- as.factor(groups)
inputs <- deviationScores(object)
# Brown-Forsythe test
inputs <- lapply(split(seq_len(ncol(object)),groups),
function(x) deviationScores(object)[,x])
reshaped <- lapply(seq_len(nrow(inputs[[1]])),
function(x) lapply(inputs, function(y) y[x, ]))
if (parametric) {
p_val <- vapply(reshaped, function(x) do.call(bf_var_test, x), 0)
p_val <- apply(inputs, 1, bf_var_test, groups)
} else {
p_val <- vapply(reshaped, function(x) do.call(bf_kw_var_test, x), 0)
p_val <- apply(inputs, 1, bf_kw_var_test, groups)
}
p_adj <- p.adjust(p_val, method = "BH")
return(data.frame(p_value = p_val, p_value_adjusted = p_adj))
}

bf_var_test <- function(...) {
inputs <- list(...)
medians <- vapply(inputs, median, 0, na.rm = TRUE)
median_diff <- do.call(c, lapply(seq_along(inputs),
function(x) abs(inputs[[x]] -
medians[x])))
group <- factor(do.call(c, lapply(seq_along(inputs),
function(x) rep(x, length(inputs[[x]])))))
return(anova(lm(median_diff ~ group))[1, 5])
bf_var_test <- function(x, groups) {
medians <- aggregate(x, list(groups), median)$x
median_diff <- abs(x - unsplit(medians, groups))
return(anova(lm(median_diff ~ groups))[1, 5])
}

bf_kw_var_test <- function(...) {
inputs <- list(...)
medians <- vapply(inputs, median, 0, na.rm = TRUE)
median_diff <- do.call(c, lapply(seq_along(inputs),
function(x) abs(inputs[[x]] -
medians[x])))
group <- factor(do.call(c, lapply(seq_along(inputs),
function(x) rep(x, length(inputs[[x]])))))
return(res <- kruskal.test(median_diff ~ group)$p.value)
bf_kw_var_test <- function(x, groups) {
medians <- aggregate(x, list(groups), median)$x
median_diff <- abs(x - unsplit(medians, groups))
return(kruskal.test(median_diff ~ groups)$p.value)
}

0 comments on commit 243d192

Please sign in to comment.