From bee7c53b8b32c8c2d650e34b410624663ee295a2 Mon Sep 17 00:00:00 2001 From: Jakob Russel Date: Mon, 3 May 2021 12:28:11 +0200 Subject: [PATCH] Fix bug in euler/venn functions which falsely omited some taxa with plot=FALSE --- DESCRIPTION | 3 ++- R/ps_euler.R | 49 +++++++++++++++++++++++++++++++++---------------- R/ps_venn.R | 49 +++++++++++++++++++++++++++++++++---------------- 3 files changed, 68 insertions(+), 33 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 383e759..f5cb460 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,8 @@ Package: MicEco Title: Various functions for microbial community data -Version: 0.9.14 +Version: 0.9.15 Authors@R: person("Jakob", "Russel", email = "russel2620@gmail.com", role = c("aut", "cre")) +Maintainer: Jakob Russel Description: Collection of functions for microbiome analyses. E.g. fitting neutral models and standardized effect sizes of phylogenetic beta diversities, and much more. Depends: R (>= 3.2.5) Imports: stats, utils, phyloseq, foreach, doSNOW, picante, vegan, snow, bbmle, Hmisc, abind, reshape2, eulerr, pheatmap diff --git a/R/ps_euler.R b/R/ps_euler.R index d8520ac..7606c71 100644 --- a/R/ps_euler.R +++ b/R/ps_euler.R @@ -37,26 +37,43 @@ ps_euler <- function(ps, group, fraction = 0, weight = FALSE, type = "percent", ps_mat <- ps_mat[, -1] ps_mat_bin <- (ps_mat>0)*1 - if(weight){ - df <- eulerr::euler(ps_mat_bin, weights = rowMeans(ps_mat)) - } else { - df <- eulerr::euler(ps_mat_bin) - } - if(plot){ + if(weight){ + df <- eulerr::euler(ps_mat_bin, weights = rowMeans(ps_mat)) + } else { + df <- eulerr::euler(ps_mat_bin) + } plot(df, quantities = list(type=type), ...) } else { + # Find taxa in all combinations + combis <- lapply(2:ncol(ps_mat), function(k) lapply(lapply(1:(ncol(combn(1:ncol(ps_mat_bin), m = k))), + function(y) ps_mat_bin[, combn(1:ncol(ps_mat_bin), m = k)[, y]]), + function(x) rownames(x[rowSums(x) >= k, , drop=FALSE]))) + + # Find taxa in singles singles <- apply(ps_mat_bin, 2, function(x) names(x[x > 0])) - combis <- do.call(c, lapply(2:ncol(ps_mat), - function(k) lapply(lapply(1:(ncol(combn(1:ncol(ps_mat_bin), m = k))), - function(y) ps_mat_bin[, combn(1:ncol(ps_mat_bin), m = k)[, y]]), - function(x) rownames(x[rowSums(x) >= k, ])))) - names(combis) <- do.call(c, lapply(2:ncol(ps_mat), function(k) apply(combn(colnames(ps_mat_bin), m = k), 2, function(x) paste(x, collapse = " & ")))) - combined <- c(lapply(seq_along(singles), function(x) setdiff(singles[[x]], do.call(c, singles[-x]))), - lapply(seq_along(combis)[1:(length(combis)-1)], function(x) setdiff(combis[[x]], do.call(c, combis[-x]))), - combis[length(combis)]) - names(combined) <- c(names(singles), names(combis)) - return(combined) + # Keep only those NOT in the same combination space + singles <- lapply(seq_along(singles), function(x) setdiff(singles[[x]], do.call(c, singles[-x]))) + combis <- lapply(combis, function(cc) lapply(seq_along(cc), function(x) setdiff(cc[[x]], do.call(c, cc[-x])))) + + # Names + names(singles) <- colnames(ps_mat_bin) + for(i in 2:ncol(ps_mat)){ + names(combis[[i-1]]) <- apply(combn(colnames(ps_mat_bin), m = i), 2, function(x) paste(x, collapse = "__")) + } + + # Recursively go through combination space from complex to simple to keep only those in unique combinations + combis <- rev(combis) + combis_new <- list() + for(i in seq_along(combis)){ + if(i == 1) { + combis_new[[i]] <- combis[[i]] + } else { + combis_new[[i]] <- lapply(combis[[i]], function(x) setdiff(x, unlist(combis_new))) + } + } + combis_new <- c(singles, unlist(combis_new, recursive = FALSE)) + return(combis_new[sapply(combis_new, function(x) length(x)>0)]) } } diff --git a/R/ps_venn.R b/R/ps_venn.R index 6d372a5..503af4d 100644 --- a/R/ps_venn.R +++ b/R/ps_venn.R @@ -37,26 +37,43 @@ ps_venn <- function(ps, group, fraction = 0, weight = FALSE, type = "percent", r ps_mat <- ps_mat[, -1] ps_mat_bin <- (ps_mat>0)*1 - if(weight){ - df <- eulerr::venn(ps_mat_bin, weights = rowMeans(ps_mat)) - } else { - df <- eulerr::venn(ps_mat_bin) - } - if(plot){ + if(weight){ + df <- eulerr::venn(ps_mat_bin, weights = rowMeans(ps_mat)) + } else { + df <- eulerr::venn(ps_mat_bin) + } plot(df, quantities = list(type=type), ...) } else { + # Find taxa in all combinations + combis <- lapply(2:ncol(ps_mat), function(k) lapply(lapply(1:(ncol(combn(1:ncol(ps_mat_bin), m = k))), + function(y) ps_mat_bin[, combn(1:ncol(ps_mat_bin), m = k)[, y]]), + function(x) rownames(x[rowSums(x) >= k, , drop=FALSE]))) + + # Find taxa in singles singles <- apply(ps_mat_bin, 2, function(x) names(x[x > 0])) - combis <- do.call(c, lapply(2:ncol(ps_mat), - function(k) lapply(lapply(1:(ncol(combn(1:ncol(ps_mat_bin), m = k))), - function(y) ps_mat_bin[, combn(1:ncol(ps_mat_bin), m = k)[, y]]), - function(x) rownames(x[rowSums(x) >= k, ])))) - names(combis) <- do.call(c, lapply(2:ncol(ps_mat), function(k) apply(combn(colnames(ps_mat_bin), m = k), 2, function(x) paste(x, collapse = " & ")))) - combined <- c(lapply(seq_along(singles), function(x) setdiff(singles[[x]], do.call(c, singles[-x]))), - lapply(seq_along(combis)[1:(length(combis)-1)], function(x) setdiff(combis[[x]], do.call(c, combis[-x]))), - combis[length(combis)]) - names(combined) <- c(names(singles), names(combis)) - return(combined) + # Keep only those NOT in the same combination space + singles <- lapply(seq_along(singles), function(x) setdiff(singles[[x]], do.call(c, singles[-x]))) + combis <- lapply(combis, function(cc) lapply(seq_along(cc), function(x) setdiff(cc[[x]], do.call(c, cc[-x])))) + + # Names + names(singles) <- colnames(ps_mat_bin) + for(i in 2:ncol(ps_mat)){ + names(combis[[i-1]]) <- apply(combn(colnames(ps_mat_bin), m = i), 2, function(x) paste(x, collapse = "__")) + } + + # Recursively go through combination space from complex to simple to keep only those in unique combinations + combis <- rev(combis) + combis_new <- list() + for(i in seq_along(combis)){ + if(i == 1) { + combis_new[[i]] <- combis[[i]] + } else { + combis_new[[i]] <- lapply(combis[[i]], function(x) setdiff(x, unlist(combis_new))) + } + } + combis_new <- c(singles, unlist(combis_new, recursive = FALSE)) + return(combis_new[sapply(combis_new, function(x) length(x)>0)]) } }