From 1ea939dd78eabda4f7ba1aef76a12a0000b937ab Mon Sep 17 00:00:00 2001
From: Etienne Rifa <etienne.rifa[at]insa-toulouse.fr>
Date: Wed, 8 Sep 2021 17:07:36 +0200
Subject: [PATCH 1/9] update
---
R/bars_fun.R | 30 ++++++++++++++++++++++++------
R/generate_phyloseq_fun.R | 4 ++--
2 files changed, 26 insertions(+), 8 deletions(-)
diff --git a/R/bars_fun.R b/R/bars_fun.R
index d131ce3..0fe3171 100644
--- a/R/bars_fun.R
+++ b/R/bars_fun.R
@@ -41,6 +41,7 @@ rarefaction <- function(data = data, col = NULL, step = 100, ggplotly = TRUE){
#' @param Fact1 Variable used to change X axis tick labels and color (when split = FALSE)
#' @param split if TRUE make a facet_wrap like grouped by Ord1 (default FALSE)
#' @param relative Plot relative (TRUE, default) or raw abundance plot (FALSE)
+#' @param autoorder Order xaxis with gtools::mixedorder function (TRUE).
#' @param ylab Y axis title ("Abundance")
#' @param outfile Output html file.
#'
@@ -57,7 +58,7 @@ rarefaction <- function(data = data, col = NULL, step = 100, ggplotly = TRUE){
bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 = NULL, split = FALSE,
- relative = TRUE, ylab = "Abundance", outfile="plot_compo.html", verbose = TRUE){
+ relative = TRUE, autoorder = TRUE, ylab = "Abundance", outfile="plot_compo.html", verbose = TRUE){
if(verbose){
invisible(flog.threshold(INFO))
@@ -76,13 +77,17 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
# print("get data")
sdata <- as.data.frame(sample_data(psobj.top), stringsAsFactors = TRUE)
- sdata$sample.id = sample_names(psobj.top)
+ # sdata$sample.id = sample_names(psobj.top)
otable = as.data.frame(otu_table(psobj.top))
row.names(otable) = tax_table(psobj.top)[,rank]
# print("melt data")
dat <- as.data.frame(t(otable))
dat <- cbind.data.frame(sdata, dat)
+
+ fun = glue( "dat <- dat[levels(sdata$sample.id), ]")
+ eval(parse(text=fun))
+
meltdat <- reshape2::melt(dat, id.vars=1:ncol(sdata))
tt <- levels(meltdat$variable)
meltdat$variable <- factor(meltdat$variable, levels= c("Other", tt[tt!="Other"]))
@@ -92,16 +97,26 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
# print(levels(meltdat$sample.id))
# save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
- # TODO: Message d'erreur si factor n'est pas dans les sample_data
+
+
+ if(autoorder){
+ fun = glue( "labs = gtools::mixedorder(as.character(meltdat${Ord1}))" )
+ eval(parse(text=fun))
+ }else{
+ labs = 1:nrow(meltdat)
+ }
+
+# print(unique(meltdat$sample.id[labs]))
fun = glue( "xform <- list(categoryorder = 'array',
- categoryarray = unique(meltdat$sample.id[gtools::mixedorder(as.character(meltdat${Ord1}))]),
+ categoryarray = unique(meltdat$sample.id[labs]),
title = 'Samples',
tickmode = 'array',
tickvals = 0:nrow(sdata),
- ticktext = sdata[unique(meltdat$sample.id[gtools::mixedorder(as.character(meltdat${Ord1}))]), '{Fact1}']@.Data[[1]],
+ ticktext = sdata[as.character(unique(meltdat$sample.id[labs])), '{Fact1}']@.Data[[1]],
tickangle = -90)")
eval(parse(text=fun))
+ # print(xform)
# subplot to vizualize groups
@@ -134,6 +149,9 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
otable=apply(otable,2, function(x){Tot=sum(x); x/Tot})
dat= as.data.frame(t(otable))
dat <- cbind.data.frame(sdata, dat)
+ fun = glue( "dat <- dat[levels(sdata$sample.id), ]")
+ eval(parse(text=fun))
+
meltdat <- reshape2::melt(dat, id.vars=1:ncol(sdata))
tt <- levels(meltdat$variable)
meltdat$variable <- factor(meltdat$variable, levels= c("Other", tt[tt!="Other"]))
@@ -177,7 +195,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
keep = TRUE) %>%
subplot(nrows = 1, shareX = TRUE, shareY=TRUE, titleX = FALSE) %>%
layout(title="",
- xaxis = list(title = glue("{Ord1} = {unique(meltdat[, Ord1])[1]}")),
+ xaxis = list(title = glue("{Ord1} = {levels(meltdat[, Ord1])[1]}")),
yaxis = list(title = ylab),
barmode = 'stack')
diff --git a/R/generate_phyloseq_fun.R b/R/generate_phyloseq_fun.R
index fce0935..7f84e6b 100644
--- a/R/generate_phyloseq_fun.R
+++ b/R/generate_phyloseq_fun.R
@@ -53,10 +53,10 @@ generate_phyloseq_fun <- function(dada_res = dada_res, tax.table = tax.table, tr
if(is.null(tree)){
flog.info('Building phyloseq object without tree...')
- data <- phyloseq(dada_res$otu.table, tax_table(as.matrix(tax.table)), sample.metadata, DNAStringSet(sequences))
+ data <- phyloseq(dada_res$otu.table[,row.names(sample.metadata)], tax_table(as.matrix(tax.table)), sample.metadata, DNAStringSet(sequences))
} else{
flog.info('Building phyloseq object with tree...')
- data <- phyloseq(dada_res$otu.table, tax_table(as.matrix(tax.table)), sample.metadata, phy_tree(tree), DNAStringSet(sequences))
+ data <- phyloseq(dada_res$otu.table[,row.names(sample.metadata)], tax_table(as.matrix(tax.table)), sample.metadata, phy_tree(tree), DNAStringSet(sequences))
}
flog.info('Done.')
data_rel <- transform_sample_counts(data, function(x) x / sum(x) )
--
GitLab
From 390a5af8b07b4ce6adfa6fd5d4c0cdff0ebab027 Mon Sep 17 00:00:00 2001
From: Etienne Rifa <etienne.rifa[at]insa-toulouse.fr>
Date: Wed, 15 Sep 2021 16:05:57 +0200
Subject: [PATCH 2/9] msgs
---
R/bars_fun.R | 7 ++++++-
1 file changed, 6 insertions(+), 1 deletion(-)
diff --git a/R/bars_fun.R b/R/bars_fun.R
index 0fe3171..d058821 100644
--- a/R/bars_fun.R
+++ b/R/bars_fun.R
@@ -88,6 +88,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
fun = glue( "dat <- dat[levels(sdata$sample.id), ]")
eval(parse(text=fun))
+ flog.info(' Melting table...')
meltdat <- reshape2::melt(dat, id.vars=1:ncol(sdata))
tt <- levels(meltdat$variable)
meltdat$variable <- factor(meltdat$variable, levels= c("Other", tt[tt!="Other"]))
@@ -98,6 +99,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
# save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
+ flog.info(' Ordering samples...')
if(autoorder){
fun = glue( "labs = gtools::mixedorder(as.character(meltdat${Ord1}))" )
@@ -107,7 +109,8 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
}
# print(unique(meltdat$sample.id[labs]))
-
+ save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
+ flog.info(' Some treatment 1...')
fun = glue( "xform <- list(categoryorder = 'array',
categoryarray = unique(meltdat$sample.id[labs]),
title = 'Samples',
@@ -120,6 +123,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
# subplot to vizualize groups
+ flog.info(' Some treatment 2...')
orderedIDS <- unique(meltdat$sample.id[gtools::mixedorder(as.character(meltdat[,Ord1]))])
orderedOrd1 <- meltdat[,Ord1][gtools::mixedorder(as.character(meltdat[,Ord1]))]
@@ -132,6 +136,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
fun = glue( "meltdat${Ord1} <- factor(meltdat${Ord1}, levels = as.character(unique(orderedOrd1)))")
eval(parse(text=fun))
+ flog.info(' Plot 1 ...')
subp1 <- df1 %>% plot_ly(
type = 'bar',
x = ~x,
--
GitLab
From 4daaa3f5c9a3fcb0142b8b97c9bf384a1de52808 Mon Sep 17 00:00:00 2001
From: Etienne Rifa <etienne.rifa[at]insa-toulouse.fr>
Date: Wed, 15 Sep 2021 17:42:29 +0200
Subject: [PATCH 3/9] update add aggregate_top_taxa
---
NAMESPACE | 4 +++-
R/bars_fun.R | 29 ++++++++++++++++++++++++++++-
man/aggregate_top_taxa.Rd | 18 ++++++++++++++++++
man/bars_fun.Rd | 3 +++
man/diversity_beta_light.Rd | 8 +++++++-
5 files changed, 59 insertions(+), 3 deletions(-)
create mode 100644 man/aggregate_top_taxa.Rd
diff --git a/NAMESPACE b/NAMESPACE
index 5ce1a3d..825937d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,6 +3,7 @@
export(ASVenn_fun)
export(VENNFUN)
export(aggregate_fun)
+export(aggregate_top_taxa)
export(assign_taxo_fun)
export(bars_fun)
export(check_tax_fun)
@@ -84,8 +85,9 @@ importFrom(metacoder,heat_tree)
importFrom(metacoder,parse_phyloseq)
importFrom(metacoder,taxon_names)
importFrom(metacoder,zero_low_counts)
-importFrom(microbiome,aggregate_top_taxa)
+importFrom(microbiome,aggregate_taxa)
importFrom(microbiome,core_members)
+importFrom(microbiome,top_taxa)
importFrom(microbiome,transform)
importFrom(nlme,lme)
importFrom(phangorn,NJ)
diff --git a/R/bars_fun.R b/R/bars_fun.R
index d058821..cb5ea50 100644
--- a/R/bars_fun.R
+++ b/R/bars_fun.R
@@ -31,6 +31,32 @@ rarefaction <- function(data = data, col = NULL, step = 100, ggplotly = TRUE){
}
+#' aggregate_top_taxa from microbiome package
+#'
+#'
+#' @param x phyloseq object
+#' @param top Keep the top-n taxa, and merge the rest under the category 'Other'. Instead of top-n numeric this can also be a character vector listing the groups to combine.
+#' @param level Summarization level (from ‘rank_names(pseq)’)
+#'
+#' @importFrom microbiome aggregate_taxa
+#' @importFrom microbiome top_taxa
+#'
+#' @export
+
+aggregate_top_taxa <- function (x, top, level){
+ x <- aggregate_taxa(x, level)
+ tops <- top_taxa(x, top)
+ tax <- tax_table(x)
+ inds <- which(!rownames(tax) %in% tops)
+ tax[inds, level] <- "Other"
+ tax_table(x) <- tax
+ tt <- tax_table(x)[, level]
+ tax_table(x) <- tax_table(tt)
+ aggregate_taxa(x, level)
+}
+
+
+
#' Barplots plotly
#'
#'
@@ -48,7 +74,6 @@ rarefaction <- function(data = data, col = NULL, step = 100, ggplotly = TRUE){
#' @return Returns barplots in an interactive plotly community plot
#'
#' @import plotly
-#' @importFrom microbiome aggregate_top_taxa
#' @importFrom reshape2 melt
#' @importFrom gtools mixedsort
#' @importFrom dplyr group_map group_by across
@@ -75,6 +100,8 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
Fdata = data
psobj.top <- aggregate_top_taxa(Fdata, rank, top = top)
+
+
# print("get data")
sdata <- as.data.frame(sample_data(psobj.top), stringsAsFactors = TRUE)
# sdata$sample.id = sample_names(psobj.top)
diff --git a/man/aggregate_top_taxa.Rd b/man/aggregate_top_taxa.Rd
new file mode 100644
index 0000000..82707c5
--- /dev/null
+++ b/man/aggregate_top_taxa.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/bars_fun.R
+\name{aggregate_top_taxa}
+\alias{aggregate_top_taxa}
+\title{aggregate_top_taxa from microbiome package}
+\usage{
+aggregate_top_taxa(x, top, level)
+}
+\arguments{
+\item{x}{phyloseq object}
+
+\item{top}{Keep the top-n taxa, and merge the rest under the category 'Other'. Instead of top-n numeric this can also be a character vector listing the groups to combine.}
+
+\item{level}{Summarization level (from ‘rank_names(pseq)’)}
+}
+\description{
+aggregate_top_taxa from microbiome package
+}
diff --git a/man/bars_fun.Rd b/man/bars_fun.Rd
index 3719792..7ff357c 100644
--- a/man/bars_fun.Rd
+++ b/man/bars_fun.Rd
@@ -12,6 +12,7 @@ bars_fun(
Fact1 = NULL,
split = FALSE,
relative = TRUE,
+ autoorder = TRUE,
ylab = "Abundance",
outfile = "plot_compo.html",
verbose = TRUE
@@ -32,6 +33,8 @@ bars_fun(
\item{relative}{Plot relative (TRUE, default) or raw abundance plot (FALSE)}
+\item{autoorder}{Order xaxis with gtools::mixedorder function (TRUE).}
+
\item{ylab}{Y axis title ("Abundance")}
\item{outfile}{Output html file.}
diff --git a/man/diversity_beta_light.Rd b/man/diversity_beta_light.Rd
index 87425ad..b238096 100644
--- a/man/diversity_beta_light.Rd
+++ b/man/diversity_beta_light.Rd
@@ -12,7 +12,9 @@ diversity_beta_light(
dist0 = "bray",
ord0 = "MDS",
output = "./plot_div_beta/",
- tests = TRUE
+ axes = c(1, 2),
+ tests = TRUE,
+ verbose = 2
)
}
\arguments{
@@ -30,7 +32,11 @@ diversity_beta_light(
\item{output}{The output file directory.}
+\item{axes}{Axes to plot (c(1,2))}
+
\item{tests}{Whether to compute tests or not (TRUE/FALSE)}
+
+\item{verbose}{Verbose level. (1: quiet, 2: print infos, 3: print infos + debug)}
}
\value{
Return specific plots and tests in list and output them in the output directory.
--
GitLab
From 84c13bb80c3be087ea156aa6877b4bdeff61c9f9 Mon Sep 17 00:00:00 2001
From: Etienne Rifa <etienne.rifa[at]insa-toulouse.fr>
Date: Fri, 17 Sep 2021 09:58:10 +0200
Subject: [PATCH 4/9] update bars
---
R/bars_fun.R | 49 +++++++++++++++++++++++++++----------------------
1 file changed, 27 insertions(+), 22 deletions(-)
diff --git a/R/bars_fun.R b/R/bars_fun.R
index cb5ea50..76a07c8 100644
--- a/R/bars_fun.R
+++ b/R/bars_fun.R
@@ -64,7 +64,7 @@ aggregate_top_taxa <- function (x, top, level){
#' @param rank Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)
#' @param top Number of top taxa to plot
#' @param Ord1 Variable used to order sample (X axis) or split the barplot if split = TRUE
-#' @param Fact1 Variable used to change X axis tick labels and color (when split = FALSE)
+#' @param sample_labels If true, x axis labels are sample IDS, if false labels displayed are levels from Ord1 argument. (FALSE)
#' @param split if TRUE make a facet_wrap like grouped by Ord1 (default FALSE)
#' @param relative Plot relative (TRUE, default) or raw abundance plot (FALSE)
#' @param autoorder Order xaxis with gtools::mixedorder function (TRUE).
@@ -82,7 +82,7 @@ aggregate_top_taxa <- function (x, top, level){
#' @export
-bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 = NULL, split = FALSE,
+bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, sample_labels = FALSE, split = FALSE,
relative = TRUE, autoorder = TRUE, ylab = "Abundance", outfile="plot_compo.html", verbose = TRUE){
if(verbose){
@@ -92,8 +92,8 @@ bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 =
}
-if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
- stop(paste("Wrong value in Ord1 or Fact1 arguments, please use variables existing in the phyloseq object:", toString(sample_variables(data))))
+if( all(Ord1 != sample_variables(data))){
+ stop(paste("Wrong value in Ord1, please use variables existing in the phyloseq object:", toString(sample_variables(data))))
}
#{paste(sample_variables(data), collapse = ",")}
flog.info('Preprocess...')
@@ -104,7 +104,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
# print("get data")
sdata <- as.data.frame(sample_data(psobj.top), stringsAsFactors = TRUE)
- # sdata$sample.id = sample_names(psobj.top)
+ sdata$sample.id = sample_names(psobj.top)
otable = as.data.frame(otu_table(psobj.top))
row.names(otable) = tax_table(psobj.top)[,rank]
@@ -112,8 +112,8 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
dat <- as.data.frame(t(otable))
dat <- cbind.data.frame(sdata, dat)
- fun = glue( "dat <- dat[levels(sdata$sample.id), ]")
- eval(parse(text=fun))
+ # fun = glue( "dat <- dat[levels(sdata$sample.id), ]")
+ # eval(parse(text=fun))
flog.info(' Melting table...')
meltdat <- reshape2::melt(dat, id.vars=1:ncol(sdata))
@@ -123,7 +123,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
LL=list()
# print(head(meltdat))
# print(levels(meltdat$sample.id))
- # save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
+ save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
flog.info(' Ordering samples...')
@@ -136,17 +136,21 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
}
# print(unique(meltdat$sample.id[labs]))
- save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
- flog.info(' Some treatment 1...')
- fun = glue( "xform <- list(categoryorder = 'array',
- categoryarray = unique(meltdat$sample.id[labs]),
- title = 'Samples',
- tickmode = 'array',
- tickvals = 0:nrow(sdata),
- ticktext = sdata[as.character(unique(meltdat$sample.id[labs])), '{Fact1}']@.Data[[1]],
- tickangle = -90)")
- eval(parse(text=fun))
- # print(xform)
+ if(sample_labels){
+ lab1 = "sample.id"
+ }else{
+ lab1 = Ord1
+ }
+ print(lab1)
+
+ flog.info(' Set labels...')
+ xform <- list(categoryorder = 'array',
+ categoryarray = unique(meltdat$sample.id[labs]),
+ title = 'Samples',
+ tickmode = 'array',
+ tickvals = 0:nrow(sdata),
+ ticktext = sdata[as.character(unique(meltdat$sample.id[labs])), lab1]@.Data[[1]],
+ tickangle = -90)
# subplot to vizualize groups
@@ -155,7 +159,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
orderedOrd1 <- meltdat[,Ord1][gtools::mixedorder(as.character(meltdat[,Ord1]))]
df1 <- cbind.data.frame(x=sdata[orderedIDS, "sample.id"]@.Data[[1]],
- g=sdata[orderedIDS, Fact1]@.Data[[1]],
+ g=sdata[orderedIDS, Ord1]@.Data[[1]],
y=1)
fun = glue( "df1$g <- factor(df1$g, levels = as.character(unique(orderedOrd1)))")
@@ -177,12 +181,13 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
if(relative){
flog.info('Plotting relative...')
+ save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
#relative abondance
otable=apply(otable,2, function(x){Tot=sum(x); x/Tot})
dat= as.data.frame(t(otable))
dat <- cbind.data.frame(sdata, dat)
- fun = glue( "dat <- dat[levels(sdata$sample.id), ]")
- eval(parse(text=fun))
+ # fun = glue( "dat <- dat[levels(sdata$sample.id), ]")
+ # eval(parse(text=fun))
meltdat <- reshape2::melt(dat, id.vars=1:ncol(sdata))
tt <- levels(meltdat$variable)
--
GitLab
From 68f24906b711ce303be9cbd886f335dc9625250e Mon Sep 17 00:00:00 2001
From: Etienne Rifa <etienne.rifa[at]insa-toulouse.fr>
Date: Fri, 17 Sep 2021 14:47:05 +0200
Subject: [PATCH 5/9] bars_fun update replace Fact argument by sample_labels
boolean, allow user to display sample ids rather than group they belong for
non splitted plot. autoorder argument to apply or not automatic ordering of
groups.
---
R/bars_fun.R | 41 ++++++++++++++++-------------------------
man/bars_fun.Rd | 6 +++---
2 files changed, 19 insertions(+), 28 deletions(-)
diff --git a/R/bars_fun.R b/R/bars_fun.R
index 76a07c8..2c794a9 100644
--- a/R/bars_fun.R
+++ b/R/bars_fun.R
@@ -64,10 +64,10 @@ aggregate_top_taxa <- function (x, top, level){
#' @param rank Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)
#' @param top Number of top taxa to plot
#' @param Ord1 Variable used to order sample (X axis) or split the barplot if split = TRUE
-#' @param sample_labels If true, x axis labels are sample IDS, if false labels displayed are levels from Ord1 argument. (FALSE)
+#' @param sample_labels If true, x axis labels are sample IDS, if false labels displayed are levels from Ord1 argument. Ignored if split = TRUE (FALSE)
#' @param split if TRUE make a facet_wrap like grouped by Ord1 (default FALSE)
#' @param relative Plot relative (TRUE, default) or raw abundance plot (FALSE)
-#' @param autoorder Order xaxis with gtools::mixedorder function (TRUE).
+#' @param autoorder Automatic ordering xaxis labels based on Ord1 factor levels with gtools::mixedorder function (TRUE).
#' @param ylab Y axis title ("Abundance")
#' @param outfile Output html file.
#'
@@ -95,14 +95,11 @@ bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, sample_
if( all(Ord1 != sample_variables(data))){
stop(paste("Wrong value in Ord1, please use variables existing in the phyloseq object:", toString(sample_variables(data))))
}
-#{paste(sample_variables(data), collapse = ",")}
+
flog.info('Preprocess...')
Fdata = data
psobj.top <- aggregate_top_taxa(Fdata, rank, top = top)
-
-
- # print("get data")
sdata <- as.data.frame(sample_data(psobj.top), stringsAsFactors = TRUE)
sdata$sample.id = sample_names(psobj.top)
otable = as.data.frame(otu_table(psobj.top))
@@ -112,9 +109,6 @@ if( all(Ord1 != sample_variables(data))){
dat <- as.data.frame(t(otable))
dat <- cbind.data.frame(sdata, dat)
- # fun = glue( "dat <- dat[levels(sdata$sample.id), ]")
- # eval(parse(text=fun))
-
flog.info(' Melting table...')
meltdat <- reshape2::melt(dat, id.vars=1:ncol(sdata))
tt <- levels(meltdat$variable)
@@ -123,25 +117,29 @@ if( all(Ord1 != sample_variables(data))){
LL=list()
# print(head(meltdat))
# print(levels(meltdat$sample.id))
- save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
-
+ # save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
- flog.info(' Ordering samples...')
if(autoorder){
+ flog.info(' Ordering samples...')
fun = glue( "labs = gtools::mixedorder(as.character(meltdat${Ord1}))" )
eval(parse(text=fun))
+
+ orderedIDS <- unique(meltdat$sample.id[gtools::mixedorder(as.character(meltdat[,Ord1]))])
+ orderedOrd1 <- meltdat[,Ord1][gtools::mixedorder(as.character(meltdat[,Ord1]))]
+ orderedOrd1 <- factor(orderedOrd1, levels = gtools::mixedsort(levels(orderedOrd1)))
}else{
labs = 1:nrow(meltdat)
+
+ orderedIDS <- unique(meltdat$sample.id)
+ orderedOrd1 <- meltdat[,Ord1]
}
-# print(unique(meltdat$sample.id[labs]))
if(sample_labels){
lab1 = "sample.id"
}else{
lab1 = Ord1
}
- print(lab1)
flog.info(' Set labels...')
xform <- list(categoryorder = 'array',
@@ -154,20 +152,16 @@ if( all(Ord1 != sample_variables(data))){
# subplot to vizualize groups
- flog.info(' Some treatment 2...')
- orderedIDS <- unique(meltdat$sample.id[gtools::mixedorder(as.character(meltdat[,Ord1]))])
- orderedOrd1 <- meltdat[,Ord1][gtools::mixedorder(as.character(meltdat[,Ord1]))]
-
+ flog.info(' Subplot...')
df1 <- cbind.data.frame(x=sdata[orderedIDS, "sample.id"]@.Data[[1]],
g=sdata[orderedIDS, Ord1]@.Data[[1]],
y=1)
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(unique(orderedOrd1)))")
+ fun = glue( "meltdat${Ord1} <- factor(meltdat${Ord1}, levels = as.character(levels(orderedOrd1)))")
eval(parse(text=fun))
- flog.info(' Plot 1 ...')
subp1 <- df1 %>% plot_ly(
type = 'bar',
x = ~x,
@@ -181,19 +175,16 @@ if( all(Ord1 != sample_variables(data))){
if(relative){
flog.info('Plotting relative...')
- save(list = ls(all.names = TRUE), file = "debug.rdata", envir = environment())
#relative abondance
otable=apply(otable,2, function(x){Tot=sum(x); x/Tot})
dat= as.data.frame(t(otable))
dat <- cbind.data.frame(sdata, dat)
- # fun = glue( "dat <- dat[levels(sdata$sample.id), ]")
- # eval(parse(text=fun))
meltdat <- reshape2::melt(dat, id.vars=1:ncol(sdata))
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(unique(orderedOrd1)))")
+ fun = glue( "meltdat${Ord1} <- factor(meltdat${Ord1}, levels = as.character(levels(orderedOrd1)))")
eval(parse(text=fun))
p1=plot_ly(meltdat, x = ~sample.id, y = ~value, type = 'bar', name = ~variable, color = ~variable) %>% #, color = ~variable
@@ -215,7 +206,7 @@ if( all(Ord1 != sample_variables(data))){
}
}
- # facet_wrap output
+ # Splitted plot output
if(!split) {
if(!is.null(outfile)){
htmlwidgets::saveWidget(p1, outfile)
diff --git a/man/bars_fun.Rd b/man/bars_fun.Rd
index 7ff357c..d4dc2dc 100644
--- a/man/bars_fun.Rd
+++ b/man/bars_fun.Rd
@@ -9,7 +9,7 @@ bars_fun(
rank = "Genus",
top = 10,
Ord1 = NULL,
- Fact1 = NULL,
+ sample_labels = FALSE,
split = FALSE,
relative = TRUE,
autoorder = TRUE,
@@ -27,13 +27,13 @@ bars_fun(
\item{Ord1}{Variable used to order sample (X axis) or split the barplot if split = TRUE}
-\item{Fact1}{Variable used to change X axis tick labels and color (when split = FALSE)}
+\item{sample_labels}{If true, x axis labels are sample IDS, if false labels displayed are levels from Ord1 argument. Ignored if split = TRUE (FALSE)}
\item{split}{if TRUE make a facet_wrap like grouped by Ord1 (default FALSE)}
\item{relative}{Plot relative (TRUE, default) or raw abundance plot (FALSE)}
-\item{autoorder}{Order xaxis with gtools::mixedorder function (TRUE).}
+\item{autoorder}{Automatic ordering xaxis labels based on Ord1 factor levels with gtools::mixedorder function (TRUE).}
\item{ylab}{Y axis title ("Abundance")}
--
GitLab
From 426505d9239e60f07b111c664b4339127a37daab Mon Sep 17 00:00:00 2001
From: Etienne Rifa <etienne.rifa[at]insa-toulouse.fr>
Date: Fri, 15 Oct 2021 11:24:44 +0200
Subject: [PATCH 6/9] manual+namespace update
---
NAMESPACE | 2 ++
man/bars_fun.Rd | 1 +
man/diversity_alpha_fun.Rd | 2 +-
man/diversity_beta_light.Rd | 5 ++++-
4 files changed, 8 insertions(+), 2 deletions(-)
diff --git a/NAMESPACE b/NAMESPACE
index 825937d..30872a7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -66,6 +66,7 @@ importFrom(DESeq2,resultsNames)
importFrom(dplyr,across)
importFrom(dplyr,group_by)
importFrom(dplyr,group_map)
+importFrom(dplyr,left_join)
importFrom(ggpubr,ggtexttable)
importFrom(ggpubr,text_grob)
importFrom(ggpubr,ttheme)
@@ -97,4 +98,5 @@ importFrom(phangorn,pml)
importFrom(plotly,ggplotly)
importFrom(ranacapa,ggrare)
importFrom(reshape2,melt)
+importFrom(tibble,rownames_to_column)
importFrom(venn,venn)
diff --git a/man/bars_fun.Rd b/man/bars_fun.Rd
index d4dc2dc..44373c9 100644
--- a/man/bars_fun.Rd
+++ b/man/bars_fun.Rd
@@ -11,6 +11,7 @@ bars_fun(
Ord1 = NULL,
sample_labels = FALSE,
split = FALSE,
+ split_sid_order = FALSE,
relative = TRUE,
autoorder = TRUE,
ylab = "Abundance",
diff --git a/man/diversity_alpha_fun.Rd b/man/diversity_alpha_fun.Rd
index 4f60f8d..714642e 100644
--- a/man/diversity_alpha_fun.Rd
+++ b/man/diversity_alpha_fun.Rd
@@ -20,7 +20,7 @@ diversity_alpha_fun(
\item{output}{Output directory}
-\item{column1}{Column name of first factor to test (covariable in ANOVA).}
+\item{column1}{Column name of first factor to test (covariable in ANOVA), or principal factor if only one factor to test.}
\item{column2}{Column name of second factor to test (last factor in ANOVA).}
diff --git a/man/diversity_beta_light.Rd b/man/diversity_beta_light.Rd
index b238096..357f044 100644
--- a/man/diversity_beta_light.Rd
+++ b/man/diversity_beta_light.Rd
@@ -14,7 +14,8 @@ diversity_beta_light(
output = "./plot_div_beta/",
axes = c(1, 2),
tests = TRUE,
- verbose = 2
+ verbose = 2,
+ ellipse = TRUE
)
}
\arguments{
@@ -37,6 +38,8 @@ diversity_beta_light(
\item{tests}{Whether to compute tests or not (TRUE/FALSE)}
\item{verbose}{Verbose level. (1: quiet, 2: print infos, 3: print infos + debug)}
+
+\item{ellipse}{Plot ellipse (TRUE)}
}
\value{
Return specific plots and tests in list and output them in the output directory.
--
GitLab
From 15f8df1574a69df21719bd893aa71b5c17016058 Mon Sep 17 00:00:00 2001
From: Etienne Rifa <etienne.rifa[at]insa-toulouse.fr>
Date: Fri, 15 Oct 2021 11:26:29 +0200
Subject: [PATCH 7/9] allow to generate and work without taxonomy, for
filtering before assignation
---
R/decontam_fun.R | 25 ++++++++++++++++++++---
R/generate_phyloseq_fun.R | 42 +++++++++++++++++++++++++++++++--------
2 files changed, 56 insertions(+), 11 deletions(-)
diff --git a/R/decontam_fun.R b/R/decontam_fun.R
index 0f36706..446d927 100644
--- a/R/decontam_fun.R
+++ b/R/decontam_fun.R
@@ -272,7 +272,12 @@ decontam_fun <- function(data = data, domain = "Bacteria", output = "./decontam_
TABf=cbind.data.frame( TABtest_filt, TABf )
names(TABf)[1] = names(TF)[j]
}
- TABff <- cbind(as.matrix(TABf), as.matrix(data_no_filtering@tax_table[row.names(TABf),]))
+
+ if(!is.null(data@tax_table)){
+ TABff <- cbind(as.matrix(TABf), as.matrix(data_no_filtering@tax_table[row.names(TABf),]))
+ }else{
+ TABff <- as.matrix(TABf)
+ }
write.table(TABff, file = paste(output,'/Exclu_out.csv',sep=''), sep = "\t", col.names=NA)
@@ -312,11 +317,25 @@ decontam_fun <- function(data = data, domain = "Bacteria", output = "./decontam_
flog.info(paste("AFTER FILTERING: ",nsamples(data), "samples and", ntaxa(data),"ASVs in otu table") )
flog.info('Writing raw tables.')
- write.table(cbind(otu_table(data),"Consensus Lineage" = apply(tax_table(data), 1, paste, collapse = ";"), "sequences"=as.data.frame(refseq(data)) ), paste(output,"/raw_otu-table.csv",sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
+ if(!is.null(data@tax_table)){
+ write.table(cbind(otu_table(data),"Consensus Lineage" = apply(tax_table(data), 1, paste, collapse = ";"), "sequences"=as.data.frame(refseq(data)) ), paste(output,"/raw_otu-table.csv",sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
+
+ }else{
+ write.table(cbind(otu_table(data), "sequences"=as.data.frame(refseq(data)) ), paste(output,"/raw_otu-table.csv",sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
+
+ }
flog.info('Writing relative tables.')
data_rel <- transform_sample_counts(data, function(x) x / sum(x) )
- write.table(cbind(otu_table(data_rel),"Consensus Lineage" = apply(tax_table(data_rel), 1, paste, collapse = ";"), "sequences"=as.data.frame(refseq(data_rel))),paste(output,"/relative_otu-table.csv",sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
+ if(!is.null(data@tax_table)){
+ write.table(cbind(otu_table(data_rel),"Consensus Lineage" = apply(tax_table(data_rel), 1, paste, collapse = ";"), "sequences"=as.data.frame(refseq(data_rel))),paste(output,"/relative_otu-table.csv",sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
+
+
+ }else{
+ write.table(cbind(otu_table(data_rel), "sequences"=as.data.frame(refseq(data_rel))),paste(output,"/relative_otu-table.csv",sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
+
+ }
+
flog.info('Saving R objects.')
save(data, file=paste(output,'/robjects.Rdata',sep=''))
diff --git a/R/generate_phyloseq_fun.R b/R/generate_phyloseq_fun.R
index 7f84e6b..02cb34e 100644
--- a/R/generate_phyloseq_fun.R
+++ b/R/generate_phyloseq_fun.R
@@ -42,8 +42,16 @@ generate_phyloseq_fun <- function(dada_res = dada_res, tax.table = tax.table, tr
flog.info('Done.')
if(length(setdiff(colnames(dada_res$otu.table), sampledata$sample.id)) > 0){
- flog.info(setdiff(colnames(dada_res$otu.table), sampledata$sample.id))
- stop("ERROR: number of samples in metadata differ from otu table.")
+ message("ERROR: sample names in otu table differs from sample metadata.")
+ print(setdiff(colnames(dada_res$otu.table), sampledata$sample.id))
+ stop("Please check metadata table.",call. = FALSE)
+ }
+
+ if(length(setdiff(sampledata$sample.id, colnames(dada_res$otu.table))) > 0){
+ message("WARNING: One or more samples in metadata are not present in otu table:")
+ print(setdiff(sampledata$sample.id, colnames(dada_res$otu.table)))
+
+ sample.metadata <- sample.metadata[colnames(dada_res$otu.table),]
}
flog.info("Sequences..")
@@ -51,13 +59,31 @@ generate_phyloseq_fun <- function(dada_res = dada_res, tax.table = tax.table, tr
flog.debug('Generate MD5 ids...')
names(sequences) <- sapply(sequences,digest,algo='md5')
- if(is.null(tree)){
- flog.info('Building phyloseq object without tree...')
- data <- phyloseq(dada_res$otu.table[,row.names(sample.metadata)], tax_table(as.matrix(tax.table)), sample.metadata, DNAStringSet(sequences))
- } else{
- flog.info('Building phyloseq object with tree...')
- data <- phyloseq(dada_res$otu.table[,row.names(sample.metadata)], tax_table(as.matrix(tax.table)), sample.metadata, phy_tree(tree), DNAStringSet(sequences))
+ if(!is.null(tax.table)){
+
+ if(is.null(tree)){
+ flog.info('Building phyloseq object without tree...')
+ data <- phyloseq(dada_res$otu.table[,row.names(sample.metadata)], tax_table(as.matrix(tax.table)), sample.metadata, DNAStringSet(sequences))
+ } else{
+ flog.info('Building phyloseq object with tree...')
+ data <- phyloseq(dada_res$otu.table[,row.names(sample.metadata)], tax_table(as.matrix(tax.table)), sample.metadata, phy_tree(tree), DNAStringSet(sequences))
+ }
+
+ }else{
+
+ if(is.null(tree)){
+ flog.info('Building phyloseq object without tree...')
+ data <- phyloseq(dada_res$otu.table[,row.names(sample.metadata)], sample.metadata, DNAStringSet(sequences))
+ } else{
+ flog.info('Building phyloseq object with tree...')
+ data <- phyloseq(dada_res$otu.table[,row.names(sample.metadata)], sample.metadata, phy_tree(tree), DNAStringSet(sequences))
+ }
+
+
}
+
+
+
flog.info('Done.')
data_rel <- transform_sample_counts(data, function(x) x / sum(x) )
--
GitLab
From 2cac0388266ffd98f8bff6d4a82487a65c1d913c Mon Sep 17 00:00:00 2001
From: Etienne Rifa <etienne.rifa[at]insa-toulouse.fr>
Date: Fri, 15 Oct 2021 11:27:35 +0200
Subject: [PATCH 8/9] illumina paired: manage sample with no reads after
filtering
---
R/dada2_fun.R | 42 +++++++++++++++++++++++++++++-------------
1 file changed, 29 insertions(+), 13 deletions(-)
diff --git a/R/dada2_fun.R b/R/dada2_fun.R
index b25e5c8..471dc97 100644
--- a/R/dada2_fun.R
+++ b/R/dada2_fun.R
@@ -30,6 +30,8 @@
#' @import futile.logger
#' @import digest
#' @import phyloseq
+#' @importFrom tibble rownames_to_column
+#' @importFrom dplyr left_join
#' @export
# DADA2 function
@@ -235,19 +237,27 @@ dada2_fun <- function(amplicon = "16S", path = "", outpath = "./dada2_out/", f_t
flog.info('Filtering reads...')
-
- out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(f_trunclen,r_trunclen),
+ out0 <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(f_trunclen,r_trunclen),
maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE, trimLeft=trim_l,
compress=compress, multithread=TRUE)
+ row.names(out0) = sample.names
+ out <- as.data.frame(out0) %>%tibble::rownames_to_column(var = "sample.id")
+
+
+ filtFs_out <- grep("F_filt",list.files(glue::glue("{path}/filtered/"), full.names = TRUE), value = TRUE)
+ filtRs_out <- grep("R_filt",list.files(glue::glue("{path}/filtered/"), full.names = TRUE), value = TRUE)
+
+ samples_names <- grep(".fastq",list.files("./data_arch_links/", full.names = FALSE), value = TRUE)
+
flog.info('Done.')
}
-
+
#COMMON
flog.info('Learning error model...')
- errF <- learnErrors(filtFs, multithread=TRUE)
- errR <- learnErrors(filtRs, multithread=TRUE)
+ errF <- learnErrors(filtFs_out, multithread=TRUE)
+ errR <- learnErrors(filtRs_out, multithread=TRUE)
flog.info('Done.')
@@ -269,8 +279,8 @@ dada2_fun <- function(amplicon = "16S", path = "", outpath = "./dada2_out/", f_t
flog.info('Dereplicating fastq...')
- derepFs <- derepFastq(filtFs, verbose=TRUE)
- derepRs <- derepFastq(filtRs, verbose=TRUE)
+ derepFs <- derepFastq(filtFs_out, verbose=TRUE)
+ derepRs <- derepFastq(filtRs_out, verbose=TRUE)
flog.info('Done.')
flog.info('dada2...')
@@ -285,7 +295,6 @@ dada2_fun <- function(amplicon = "16S", path = "", outpath = "./dada2_out/", f_t
seqtab <- makeSequenceTable(mergers)
-
flog.info('Removing chimeras...')
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
flog.debug(sum(seqtab.nochim))
@@ -295,19 +304,26 @@ dada2_fun <- function(amplicon = "16S", path = "", outpath = "./dada2_out/", f_t
flog.info('Done.')
- track <- cbind.data.frame(out, stockFs, stockRs, sapply(mergers, getN), rowSums(seqtab.nochim))
+ track0 <- cbind.data.frame(stockFs, stockRs, sapply(mergers, getN), rowSums(seqtab.nochim))
+ rownames(track0) <- stringr::str_remove(rownames(track0), "_F_filt.fastq")
+
# If processing a single sample, remove the sapply calls: e.g. replace sapply(dadaFs, getN) with getN(dadaFs)
- colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
- rownames(track) <- sample.names
- head(track)
- write.table(track, paste(outpath,"/read_tracking.csv",sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
+ track <- track0 %>% tibble::rownames_to_column(var="sample.id")
+
+
+ final_track <- out %>% dplyr::left_join(y = track, by = "sample.id")
+ colnames(final_track) <- c("sample.id", "input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
+ head(final_track)
+
+ write.table(final_track, paste(outpath,"/read_tracking.csv",sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
seqtab.export <- seqtab.nochim
colnames(seqtab.export) <- sapply(colnames(seqtab.export), digest::digest, algo="md5")
otu.table <- phyloseq::otu_table(t(seqtab.export), taxa_are_rows = TRUE)
+ colnames(otu.table) <- stringr::str_remove(colnames(otu.table), "_F_filt.fastq")
flog.info('Writing raw tables.')
write.table(cbind(t(seqtab.export), "Sequence" = colnames(seqtab.nochim)), paste(outpath,"/raw_otu-table.csv",sep=''), sep="\t", row.names=TRUE, col.names=NA, quote=FALSE)
--
GitLab
From a699c70a59c70f4a58d52eebb15d10476fe2eeca Mon Sep 17 00:00:00 2001
From: Etienne Rifa <etienne.rifa[at]insa-toulouse.fr>
Date: Fri, 15 Oct 2021 11:28:18 +0200
Subject: [PATCH 9/9] bars alpha beta bars: allow to keep sample order given by
metadata alpha: suppress fdr for pairwise wilcox and legend display beta:
plotly output
---
R/bars_fun.R | 25 ++++++++++++++++---------
R/diversity_alpha_fun.R | 16 +++++++++-------
R/diversity_beta_light.R | 15 ++++++++++++---
3 files changed, 37 insertions(+), 19 deletions(-)
diff --git a/R/bars_fun.R b/R/bars_fun.R
index 2c794a9..4b10019 100644
--- a/R/bars_fun.R
+++ b/R/bars_fun.R
@@ -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]}")),
diff --git a/R/diversity_alpha_fun.R b/R/diversity_alpha_fun.R
index f3b64a3..a775bb2 100644
--- a/R/diversity_alpha_fun.R
+++ b/R/diversity_alpha_fun.R
@@ -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)
diff --git a/R/diversity_beta_light.R b/R/diversity_beta_light.R
index 9f9e530..dd59873 100644
--- a/R/diversity_beta_light.R
+++ b/R/diversity_beta_light.R
@@ -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...')
--
GitLab