Commit a699c70a authored by Etienne Rifa's avatar Etienne Rifa
Browse files

bars alpha beta

bars: allow to keep sample order given by metadata
alpha: suppress fdr for pairwise wilcox and legend display
beta: plotly output
parent 2cac0388
......@@ -82,7 +82,7 @@ aggregate_top_taxa <- function (x, top, level){
#' @export
bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, sample_labels = FALSE, split = FALSE,
bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, sample_labels = FALSE, split = FALSE, split_sid_order = FALSE,
relative = TRUE, autoorder = TRUE, ylab = "Abundance", outfile="plot_compo.html", verbose = TRUE){
if(verbose){
......@@ -135,6 +135,8 @@ if( all(Ord1 != sample_variables(data))){
orderedOrd1 <- meltdat[,Ord1]
}
if(sample_labels){
lab1 = "sample.id"
}else{
......@@ -159,7 +161,7 @@ if( all(Ord1 != sample_variables(data))){
fun = glue( "df1$g <- factor(df1$g, levels = as.character(unique(orderedOrd1)))")
eval(parse(text=fun))
fun = glue( "meltdat${Ord1} <- factor(meltdat${Ord1}, levels = as.character(levels(orderedOrd1)))")
fun = glue( "meltdat${Ord1} <- factor(meltdat${Ord1}, levels = as.character(unique(orderedOrd1)))") # keep sample ids order from metadata.
eval(parse(text=fun))
subp1 <- df1 %>% plot_ly(
......@@ -184,7 +186,7 @@ if( all(Ord1 != sample_variables(data))){
tt <- levels(meltdat$variable)
meltdat$variable <- factor(meltdat$variable, levels= c("Other", tt[tt!="Other"]))
fun = glue( "meltdat${Ord1} <- factor(meltdat${Ord1}, levels = as.character(levels(orderedOrd1)))")
fun = glue( "meltdat${Ord1} <- factor(meltdat${Ord1}, levels = as.character(unique(orderedOrd1)))")
eval(parse(text=fun))
p1=plot_ly(meltdat, x = ~sample.id, y = ~value, type = 'bar', name = ~variable, color = ~variable) %>% #, color = ~variable
......@@ -215,12 +217,17 @@ if( all(Ord1 != sample_variables(data))){
return(p1)
} else {
flog.info('Splitted plot...')
p1 = meltdat %>% group_by(across({Ord1})) %>%
dplyr::group_map(~ plot_ly(data=., x = ~sample.id, y = ~value, type = 'bar',
name = ~variable,
color = ~variable, legendgroup = ~variable,
showlegend = (.y == levels(meltdat[, Ord1])[1])),
keep = TRUE) %>%
if(split_sid_order){
meltdat$sample.id = factor(meltdat$sample.id, levels = unique(meltdat$sample.id)) # keep sample ids order from metadata.
}
p1 = meltdat %>% arrange(across({Ord1})) %>% group_by(across({Ord1})) %>%
mutate(across(where(is.character), as.factor)) %>%
dplyr::group_map(~ plot_ly(data=., x = ~sample.id, y = ~value, type = 'bar',
name = ~variable,
color = ~variable, legendgroup = ~variable,
showlegend = (.y == levels(meltdat[, Ord1])[1])),
keep = TRUE) %>%
subplot(nrows = 1, shareX = TRUE, shareY=TRUE, titleX = FALSE) %>%
layout(title="",
xaxis = list(title = glue("{Ord1} = {levels(meltdat[, Ord1])[1]}")),
......
......@@ -15,10 +15,11 @@ alphaPlot <- function(data = data, col1 = "", col2 = "", measures = c("Shannon")
p <- plot_richness(data,x=col1, color=col2, measures=measures)
}
p$layers <- p$layers[-1]
if(col2 != ""){legpos = "right"}else{legpos = "none"}
p <- p + ggtitle('Alpha diversity indexes') + geom_boxplot(alpha = 1, outlier.shape = NA) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust=1),
legend.position = "none",axis.text=element_text(size=15),
legend.position = legpos, axis.text=element_text(size=15),
axis.title=element_text(size=16,face="bold"),
strip.text.x = element_text(size = 18,face="bold"),
title=element_text(size=16,face="bold"))
......@@ -44,7 +45,7 @@ alphaPlotly <- function(data=data, alpha=alpha, col1='', col2='', measures=c("Sh
#'
#' @param data a phyloseq object (output from decontam or generate_phyloseq)
#' @param output Output directory
#' @param column1 Column name of first factor to test (covariable in ANOVA).
#' @param column1 Column name of first factor to test (covariable in ANOVA), or principal factor if only one factor to test.
#' @param column2 Column name of second factor to test (last factor in ANOVA).
#' @param column3 Column name of subjects in case of repeated mesures (mixed model)
#' @param supcovs One or more supplementary covariables to test in anova (provided as comma separated vector)
......@@ -104,8 +105,8 @@ diversity_alpha_fun <- function(data = data, output = "./plot_div_alpha/", colum
#alpha diversité tableau
resAlpha = list()
flog.info('Alpha diversity tab ...')
resAlpha$alphatable <- estimate_richness(data, measures = measures )
# row.names(resAlpha$alphatable) <- gsub("X","",row.names(resAlpha$alphatable))
write.table(resAlpha$alphatable,paste(output,'/alphaDiversity_table.csv',sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
flog.info('Done.')
......@@ -115,9 +116,10 @@ diversity_alpha_fun <- function(data = data, output = "./plot_div_alpha/", colum
resAlpha$plot = p
ggsave(paste(output,'/alpha_diversity.eps',sep=''), plot=p, height = 15, width = 30, units="cm", dpi = 500, device="eps")
anova_data <- cbind(sample_data(data), resAlpha$alphatable)
alphatable <- estimate_richness(data, measures = measures )
anova_data <- cbind(sample_data(data), alphatable)
anova_data$Depth <- sample_sums(data)
resAlpha$alphatable <- anova_data
if(length(levels(as.factor(anova_data[,column1])))>1 & mean(table(anova_data[,column1])) != 1 ){
flog.info('ANOVA ...')
# variables <- paste(sep=" + ", "Depth", var1)
......@@ -163,7 +165,7 @@ diversity_alpha_fun <- function(data = data, output = "./plot_div_alpha/", colum
# }
flog.info(paste("\n##pvalues of pairwise wilcox test on ", m, "with FDR correction \n"), sep="")
fun <- paste("wilcox_res1 <- pairwise.wilcox.test(anova_data$",m,", anova_data[,column1], p.adjust.method='fdr')", sep="")
fun <- paste("wilcox_res1 <- pairwise.wilcox.test(anova_data$",m,", anova_data[,column1], p.adjust.method='none')", sep="")
eval(parse(text = fun))
# print(round(wilcox_res1$p.value,3))
wilcox_col1 <- round(wilcox_res1$p.value,3)
......@@ -175,7 +177,7 @@ diversity_alpha_fun <- function(data = data, output = "./plot_div_alpha/", colum
if(column2 != ''){
flog.info(paste("\n##pvalues of pairwise wilcox test on ", m, "with FDR correction \n"), sep="")
fun <- paste("wilcox_res1 <- pairwise.wilcox.test(anova_data$",m,", anova_data[,column2], p.adjust.method='fdr')", sep="")
fun <- paste("wilcox_res1 <- pairwise.wilcox.test(anova_data$",m,", anova_data[,column2], p.adjust.method='none')", sep="")
eval(parse(text = fun))
wilcox_col2_fdr <- round(wilcox_res1$p.value,3)
......
......@@ -94,9 +94,18 @@ diversity_beta_light <- function(psobj, rank = "ASV", col = NULL, cov = NULL, di
}else{
p1 <- plot_samples(data_rank, ordinate(data_rank, ord0, dist0), color = col, axes = axes ) +
theme_bw() + ggtitle(glue::glue("{ord0} + {dist0}"))
if(ellipse){p1 <- p1 + stat_ellipse()}
p1 <- phyloseq::plot_ordination(physeq = data_rank, ordination = ordinate(data_rank, ord0, dist0), axes = axes)
p1$layers[[1]] <- NULL
sample.id = sample_names(data_rank)
sdata <- sample_data(data_rank)
fact <- as.matrix(sdata[,col])
p1 <- p1 + aes(color = fact, sample.id = sample.id)
p1 <- p1 + stat_ellipse(aes(group = fact))
p1 <- p1 + geom_point() + theme_bw()
resBeta$plotly1 <- ggplotly(p1, tooltip=c("x", "y", "sample.id")) %>% config(toImageButtonOptions = list(format = "svg"))
}
# plot(p1)
flog.info('Plot ok...')
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment