Skip to content

Commit

Permalink
fix for new version of ALDEx2 1.33
Browse files Browse the repository at this point in the history
  • Loading branch information
@yiluheihei
yiluheihei committed Nov 4, 2023
1 parent e4fea2f commit 577f9da
Show file tree
Hide file tree
Showing 2 changed files with 101 additions and 62 deletions.
94 changes: 52 additions & 42 deletions R/DA-aldex.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,9 +91,9 @@ run_aldex <- function(ps,
denom = c("all", "iqlr", "zero", "lvha"),
paired = FALSE) {
stopifnot(inherits(ps, "phyloseq"))
ps <- check_rank_names(ps) %>%
ps <- check_rank_names(ps) %>%
check_taxa_rank( taxa_rank)

denom <- match.arg(denom, c("all", "iqlr", "zero", "lvha"))
p_adjust <- match.arg(
p_adjust,
Expand All @@ -109,9 +109,9 @@ run_aldex <- function(ps,
c("t.test", "wilcox.test", "kruskal", "glm_anova")
)
if (method %in% c("t.test", "wilcox.test")) {
test <- "t"
test_method <- "t"
} else {
test <- "kw"
test_method <- "kw"
}

# check whether group is valid, write a function
Expand All @@ -137,7 +137,7 @@ run_aldex <- function(ps,
groups <- sample_meta[[group]]
abd <- abundances(ps_summarized, norm = TRUE)

test_fun <- ifelse(test == "t", aldex_t, aldex_kw)
test_fun <- ifelse(test_method == "t", aldex_t, aldex_kw)
test_para <- list(
reads = abd,
conditions = groups,
Expand All @@ -146,7 +146,7 @@ run_aldex <- function(ps,
denom = denom,
p_adjust = p_adjust
)
if (test == "t") {
if (test_method == "t") {
test_para <- c(test_para, paired = paired)
}

Expand Down Expand Up @@ -217,7 +217,7 @@ aldex_t <- function(reads,
if (!inherits(reads, "aldex.clr")) {
reads_clr <- ALDEx2::aldex.clr(
reads = reads,
conds = conditions,
conds = as.character(conditions),
mc.samples = mc_samples,
denom = denom
)
Expand All @@ -234,8 +234,7 @@ aldex_t <- function(reads,
pvalue <- purrr::map_dfc(
mc_instance_ldf,
t_fast,
group = conditions, paired = paired
)
group = conditions, paired = paired)
} else {
pvalue <- purrr::map_dfc(
mc_instance_ldf,
Expand All @@ -244,10 +243,28 @@ aldex_t <- function(reads,
)
}

padj <- purrr::map_dfc(pvalue, p.adjust, method = p_adjust)
# expect value
e_pvalue <- rowMeans(pvalue)
e_padj <- rowMeans(padj)
padj_greater <- purrr::map_dfc(
pvalue,
\(x) p.adjust (2 * x, method = p_adjust)
)
padj_less <- purrr::map_dfc(
pvalue,
\(x) p.adjust (2 * (1 - x), method = p_adjust)
)

# making this into a two-sided test
pvalue_greater <-2 * pvalue
pvalue_less <- 2 * (1 - pvalue)

# making sure the max p-value is 1
pvalue_greater <- apply(pvalue_greater, c(1, 2), \(x) min(x, 1))
pvalue_less <- apply(pvalue_less, c(1, 2), \(x) min(x, 1))

# get the expected value of p value and adjusted p value
e_pvalue <- cbind(rowMeans(pvalue_greater), rowMeans(pvalue_less)) |>
apply(1, min)
e_padj <- cbind(rowMeans(padj_greater), rowMeans(padj_less)) |>
apply(1, min)

# effect size
ef <- ALDEx2::aldex.effect(
Expand Down Expand Up @@ -404,7 +421,7 @@ convert_instance <- function(mc_instance, mc_samples) {
}


# fast test function modified from ALDEx2
# fast test function modified from ALDEx2::t.fast
#' @importFrom stats pt
t_fast <- function(x, group, paired = FALSE) {
grp1 <- group == unique(group)[1]
Expand All @@ -431,21 +448,19 @@ t_fast <- function(x, group, paired = FALSE) {
nonpara = "n"
)
df <- length(idx1) - 1
res <- pt(abs(t), df = df, lower.tail = FALSE) * 2
} else {
t <- multtest::mt.teststat(x,
as.numeric(grp1),
test = "t",
t <- multtest::mt.teststat(x,
as.numeric(grp1),
test = "t",
nonpara = "n"
)
s1 <- apply(x[, grp1], 1, sd)
s2 <- apply(x[, grp2], 1, sd)
df <- ((s1^2 / n1 + s2^2 / n2)^2) / ((s1^2 / n1)^2 / (n1 - 1) +
df <- ((s1^2 / n1 + s2^2 / n2)^2) / ((s1^2 / n1)^2 / (n1 - 1) +
(s2^2 / n2)^2 / (n2 - 1))
res <- pt(abs(t), df = df, lower.tail = FALSE) * 2
}

res
return(pt(t, df = df, lower.tail = FALSE))
}

# wilcox.fast function replaces wilcox.test
Expand All @@ -457,6 +472,7 @@ t_fast <- function(x, group, paired = FALSE) {
# * uses multtest
#' @importFrom stats psignrank pnorm pwilcox wilcox.test
wilcox_fast <- function(x, group, paired = FALSE) {
stopifnot(ncol(x) == length(group))
grp1 <- group == unique(group)[1]
grp2 <- group == unique(group)[2]
n1 <- sum(grp1)
Expand All @@ -466,17 +482,12 @@ wilcox_fast <- function(x, group, paired = FALSE) {
xt <- t(x)
if (paired) {
any_ties <- any(
apply(
xt[grp1, ] - xt[grp2, ], 2,
function(y) length(unique(y))
) != ncol(x) / 2
apply(xt[grp1, ] - xt[grp2, ], 2, function(y) length(unique(y))) !=
ncol(x) / 2
)
} else {
any_ties <- any(
apply(
xt, 2,
function(y) length(unique(y))
) != ncol(x)
apply(xt, 2, function(y) length(unique(y))) != ncol(x)
)
}

Expand All @@ -487,39 +498,38 @@ wilcox_fast <- function(x, group, paired = FALSE) {
function(i) {
wilcox.test(
i[grp1], i[grp2],
paired = paired, correct = FALSE
paired = paired,
alternative = "greater",
correct = FALSE
)$p.value
}
)

return(res)
}

if (paired) {
if (n1 != n2) stop("Cannot pair uneven groups.")
x_diff <- xt[grp1, ] - xt[grp2, ]
v <- apply(x_diff, 2, function(y) sum(rank(abs(y))[y > 0]))
topscore <- (n1 * (n1 + 1)) / 2
v_lower <- ifelse(v > topscore / 2, topscore - v, v)
if (sum(grp1) < 50) {
if (sum(grp1) < 50) {
# as per wilcox test, use exact -- ASSUMES NO TIES!!
v_p <- psignrank(v_lower, n1) * 2
# psignrank returns non-zero for W = mean
res <- ifelse(v_p > 1, 1, v_p)
res <- psignrank(v - 1, n1, lower.tail = FALSE)
} else { # Use normal approximation
v_std <- (topscore / 2 - v_lower) /
v_std <- (v - topscore / 2) /
sqrt(n1 * (n1 + 1) * (2 * n1 + 1) / 24)
res <- pnorm(v_std, lower.tail = FALSE) * 2
res <- pnorm(v_std, lower.tail = FALSE)
}
} else {
w_std <- multtest::mt.teststat(x, as.numeric(grp1), test = "wilcoxon")
if (sum(grp1) < 50 && sum(grp2) < 50) {
# as per wilcox test, use exact -- ASSUMES NO TIES!!
w_var <- sqrt((n1 * n2) * (n1 + n2 + 1) / 12)
w <- abs(w_std) * w_var + (n1 * n2) / 2
w_p <- pwilcox(w - 1, n1, n2, lower.tail = FALSE) * 2
# pwilcox returns non-zero for W = mean
res <- ifelse(w_p > 1, 1, w_p)
w <- w_std * w_var + (n1 * n2) / 2
res <- pwilcox(w - 1, n1, n2, lower.tail = FALSE)
} else { # Use normal approximation
res <- pnorm(abs(w_std), lower.tail = FALSE) * 2
res <- pnorm(w_std, lower.tail = FALSE)
}
}

Expand Down
69 changes: 49 additions & 20 deletions tests/testthat/test-aldex.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,51 +25,80 @@ x <- as(phyloseq::otu_table(kostic_crc), "matrix")[1:20, 1:20]
groups <- phyloseq::sample_data(kostic_crc)[["DIAGNOSIS"]][1:20]
idx1 <- groups == "Tumor"
idx2 <- groups == "Healthy"
x_clr <- ALDEx2::aldex.clr(x, groups, mc.samples = 12)
x_clr <- suppressWarnings(ALDEx2::aldex.clr(x, groups, mc.samples = 12))
instance <- convert_instance(x_clr@analysisData, 12)[[1]]
# keep the same number of samples for two groups
idx2 <- which(idx2)[1:4]
instance <- instance[c(which(idx1), idx2)]
new_groups <- groups[c(which(idx1), idx2)]

test_that("t_fast equal to t.test", {
test_that("t_fast equal to t.test and ALDEx2:::t.fast", {
# unpaired
t_res_up <- apply(
instance, 1,
function(x) t.test(x[1:4], x[5:8], alternative = "greater")$p.value
)
expect_equal(t_res_up, t_fast(instance, new_groups))
expect_equal(
apply(
instance, 1,
function(x) t.test(x[1:4], x[5:8])$p.value
),
t_fast(instance, new_groups)
t_res_up,
ALDEx2:::t.fast(instance, new_groups, paired = FALSE)$p
)

# paired
t_res_p <- apply(
instance, 1,
function(x) t.test(x[1:4], x[5:8],
paired = TRUE,
alternative = "greater")$p.value
)
expect_equal(
apply(
instance, 1,
function(x) t.test(x[1:4], x[5:8], paired = TRUE)$p.value
),
t_res_p,
t_fast(instance, new_groups, paired = TRUE)
)
expect_equal(
t_res_p,
ALDEx2:::t.fast(instance, new_groups, paired = TRUE)$p
)

})

test_that("wilcox_fast equal to wilcox.test", {
# unpaired
wilcox_res_p <- apply(
instance, 1,
function(x) wilcox.test(x[1:4], x[5:8],
alternative = "greater")$p.value
)
expect_equal(
apply(
instance, 1,
function(x) wilcox.test(x[1:4], x[5:8])$p.value
),
wilcox_res_p,
wilcox_fast(instance, new_groups)
)
expect_equal(
wilcox_res_p,
ALDEx2:::wilcox.fast(instance, new_groups, paired = FALSE))

# paired
wilcox_res_p <- apply(
instance, 1,
function(x) wilcox.test(x[1:4], x[5:8],
paired = TRUE,
alternative = "greater")$p.value
)
expect_equal(
apply(
instance, 1,
function(x) wilcox.test(x[1:4], x[5:8], paired = TRUE)$p.value
),
wilcox_res_p,
wilcox_fast(instance, new_groups, paired = TRUE)
)
expect_equal(
wilcox_res_p,
ALDEx2:::wilcox.fast(instance, new_groups, paired = TRUE)
)

expect_equal(
apply(
instance, 1, function(x) {
wilcox.test(x[1:4], x[5:8], paired = TRUE, exact = TRUE)$p.value
wilcox.test(x[1:4], x[5:8],
alternative = "greater",
paired = TRUE, exact = TRUE)$p.value
}
),
wilcox_fast(instance, new_groups, paired = TRUE)
Expand Down

0 comments on commit 577f9da

Please sign in to comment.