Commit 99aa1bf7 authored by Etienne Rifa's avatar Etienne Rifa
Browse files

bars_fun, heatmap_fun

bars_fun: add ylab argument
heatmap: add feature to use preformated phyloseq object, add table output that can be used to create boxplot with specific taxa with normalized abundances.
parent 54820cf2
Pipeline #39731 passed with stage
in 15 seconds
...@@ -82,8 +82,8 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV", ...@@ -82,8 +82,8 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV",
if(length(lvls)==0){ if(length(lvls)==0){
flog.error('You must provide levels...') flog.error('You must provide levels...')
stop() stop()
} else if(length(lvls)>5){ } else if(length(lvls)>7){
flog.error('Venn diagram is limited to 5 levels') flog.error('Venn diagram is limited to 7 levels')
stop() stop()
} else{ } else{
if(!all(lvls %in% na.omit(levels(as.factor(sample_data(data)[,column1]@.Data[[1]])) ))){ if(!all(lvls %in% na.omit(levels(as.factor(sample_data(data)[,column1]@.Data[[1]])) ))){
...@@ -93,19 +93,6 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV", ...@@ -93,19 +93,6 @@ ASVenn_fun <- function(data = data, output = "./ASVenn/", rank = "ASV",
} }
if(length(lvls)==0){
flog.error('You must provide levels...')
stop()
} else if(length(lvls)>5){
flog.error('Venn diagram is limited to 5 levels')
stop()
} else{
if(!all(lvls %in% na.omit(levels(as.factor(sample_data(data)[,column1]@.Data[[1]])) ))){
flog.error('Your levels are not present in metadata...')
stop()
}
}
#Nombre d'espèce par matrice #Nombre d'espèce par matrice
flog.info('Parsing factor ...') flog.info('Parsing factor ...')
......
...@@ -41,6 +41,7 @@ rarefaction <- function(data = data, col = NULL, step = 100, ggplotly = TRUE){ ...@@ -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 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 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 relative Plot relative (TRUE, default) or raw abundance plot (FALSE)
#' @param ylab Y axis title ("Abundance")
#' @param outfile Output html file. #' @param outfile Output html file.
#' #'
#' @return Returns barplots in an interactive plotly community plot #' @return Returns barplots in an interactive plotly community plot
...@@ -56,8 +57,7 @@ rarefaction <- function(data = data, col = NULL, step = 100, ggplotly = TRUE){ ...@@ -56,8 +57,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, bars_fun <- function(data = data, rank = "Genus", top = 10, Ord1 = NULL, Fact1 = NULL, split = FALSE,
relative = TRUE, relative = TRUE, ylab = "Abundance", outfile="plot_compo.html", verbose = TRUE){
outfile="plot_compo.html", verbose = TRUE){
if(verbose){ if(verbose){
invisible(flog.threshold(INFO)) invisible(flog.threshold(INFO))
...@@ -125,13 +125,12 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){ ...@@ -125,13 +125,12 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
legendgroup = ~g, legendgroup = ~g,
showlegend = FALSE showlegend = FALSE
) %>% layout(xaxis = list(zeroline = FALSE,showline = FALSE, showgrid = FALSE), ) %>% layout(xaxis = list(zeroline = FALSE,showline = FALSE, showgrid = FALSE),
yaxis=list(showticklabels = FALSE,title = "",showgrid = FALSE)) yaxis=list(showticklabels = FALSE,title = ylab, showgrid = FALSE))
if(relative){ if(relative){
flog.info('Plotting relative...') flog.info('Plotting relative...')
#relative abondance #relative abondance
plottitle = "Relative abundance"
otable=apply(otable,2, function(x){Tot=sum(x); x/Tot}) otable=apply(otable,2, function(x){Tot=sum(x); x/Tot})
dat= as.data.frame(t(otable)) dat= as.data.frame(t(otable))
dat <- cbind.data.frame(sdata, dat) dat <- cbind.data.frame(sdata, dat)
...@@ -143,7 +142,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){ ...@@ -143,7 +142,7 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
eval(parse(text=fun)) eval(parse(text=fun))
p1=plot_ly(meltdat, x = ~sample.id, y = ~value, type = 'bar', name = ~variable, color = ~variable) %>% #, color = ~variable p1=plot_ly(meltdat, x = ~sample.id, y = ~value, type = 'bar', name = ~variable, color = ~variable) %>% #, color = ~variable
layout(title="Relative abundance", yaxis = list(title = 'Relative abundance'), xaxis = xform, barmode = 'stack') layout(title="", yaxis = list(title = ylab), xaxis = xform, barmode = 'stack')
if(length(df1$x) != length(unique(df1$g))){ if(length(df1$x) != length(unique(df1$g))){
p1 <- subplot(p1, subp1, nrows = 2, shareX = T, heights=c(0.95,0.05)) %>% p1 <- subplot(p1, subp1, nrows = 2, shareX = T, heights=c(0.95,0.05)) %>%
...@@ -152,9 +151,8 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){ ...@@ -152,9 +151,8 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
}else{ }else{
flog.info('Plotting raw...') flog.info('Plotting raw...')
#raw abundance #raw abundance
plottitle = "Raw abundance"
p1=plot_ly(meltdat, x = ~sample.id, y = ~value, type = 'bar', name = ~variable, color = ~variable) %>% #, color = ~variable p1=plot_ly(meltdat, x = ~sample.id, y = ~value, type = 'bar', name = ~variable, color = ~variable) %>% #, color = ~variable
layout(title="Raw abundance", yaxis = list(title = 'Raw abundance'), xaxis = xform, barmode = 'stack') layout(title="", yaxis = list(title = ylab), xaxis = xform, barmode = 'stack')
if(length(df1$x) != length(unique(df1$g))){ if(length(df1$x) != length(unique(df1$g))){
p1 <- subplot(p1, subp1, nrows = 2, shareX = T, heights=c(0.95,0.05)) %>% p1 <- subplot(p1, subp1, nrows = 2, shareX = T, heights=c(0.95,0.05)) %>%
...@@ -178,9 +176,9 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){ ...@@ -178,9 +176,9 @@ if( all(Ord1 != sample_variables(data)) | all(Fact1 != sample_variables(data))){
showlegend = (.y == levels(meltdat[, Ord1])[1])), showlegend = (.y == levels(meltdat[, Ord1])[1])),
keep = TRUE) %>% keep = TRUE) %>%
subplot(nrows = 1, shareX = TRUE, shareY=TRUE, titleX = FALSE) %>% subplot(nrows = 1, shareX = TRUE, shareY=TRUE, titleX = FALSE) %>%
layout(title=plottitle, layout(title="",
xaxis = list(title = glue("{Ord1} = {unique(meltdat[, Ord1])[1]}")), xaxis = list(title = glue("{Ord1} = {unique(meltdat[, Ord1])[1]}")),
yaxis = list(title = 'Relative abundance'), yaxis = list(title = ylab),
barmode = 'stack') barmode = 'stack')
for (i in 2:length(unique(meltdat[, Ord1]))) { for (i in 2:length(unique(meltdat[, Ord1]))) {
......
...@@ -6,7 +6,10 @@ ...@@ -6,7 +6,10 @@
#' @param output Output directory #' @param output Output directory
#' @param column1 Column name of factor to plot with (among sample_variables(data)). #' @param column1 Column name of factor to plot with (among sample_variables(data)).
#' @param top Number of top features to plot. #' @param top Number of top features to plot.
#' @param relative If TRUE, calculate relative abundance, needs a phyloseq object with raw abundance.
#' @param aggregate If TRUE, aggregate abundance of non "top" taxa in a new taxa named "Other".
#' @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 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 legend Legend title ("Abundance")
#' #'
#' @return Generate an heatmap with top taxa (html plotly version in output directory) #' @return Generate an heatmap with top taxa (html plotly version in output directory)
#' @importFrom htmlwidgets saveWidget #' @importFrom htmlwidgets saveWidget
...@@ -16,7 +19,8 @@ ...@@ -16,7 +19,8 @@
# Decontam Function # Decontam Function
heatmap_fun <- function(data = data, column1 = "", top = 20, output = "./plot_heatmap/", rank = "Species"){ heatmap_fun <- function(data = data, column1 = "", top = 20, output = "./plot_heatmap/", rank = "Species", relative = TRUE, aggregate = TRUE, legend = "Abundance"){
LL = list()
output1 <- paste(getwd(),'/',output,'/',sep='') output1 <- paste(getwd(),'/',output,'/',sep='')
if(!dir.exists(output1)){ if(!dir.exists(output1)){
...@@ -24,16 +28,37 @@ heatmap_fun <- function(data = data, column1 = "", top = 20, output = "./plot_he ...@@ -24,16 +28,37 @@ heatmap_fun <- function(data = data, column1 = "", top = 20, output = "./plot_he
} }
flog.info('Done.') flog.info('Done.')
psobj.top <- microbiome::aggregate_top_taxa(data, rank, top = top) ps.glom.rel <- psobj.top <- dataglom <- data
if(!is.null(rank)){
dataglom <- tax_glom(data, rank)
tt <- tax_table(dataglom)
taxa <- tt[,rank]
taxa_names(dataglom) <- taxa
}
otable <- otu_table(dataglom)
sdata <- as.data.frame(as.matrix(sample_data(dataglom)))
data.com <- reshape2::melt(otable)
data.com$xlabel <- as.factor(sdata[as.character(data.com$Var2),match(column1,names(sdata))])
names(data.com) <- c("Tax", "Sample", "Abundance", "xlabel")
data.com$Tax = factor(data.com$Tax, levels = sort(unique(as.character(data.com$Tax))))
if(aggregate){psobj.top <- microbiome::aggregate_top_taxa(data, rank, top = top)}
if(relative){
ps.glom.rel <- microbiome::transform(psobj.top, "compositional")
plot.composition.relAbun <- microbiome::plot_composition(ps.glom.rel, x.label = column1)
data.com <- plot.composition.relAbun$data
colnames(data.com)
}
ps.glom.rel <- microbiome::transform(psobj.top, "compositional")
plot.composition.relAbun <- microbiome::plot_composition(ps.glom.rel, x.label = column1)
data.com <- plot.composition.relAbun$data
colnames(data.com)
p.heat <- ggplot(data.com, aes(x = Sample, y = Tax)) + geom_tile(aes(fill = Abundance)) p.heat <- ggplot(data.com, aes(x = Sample, y = Tax)) + geom_tile(aes(fill = Abundance))
p.heat <- p.heat + scale_fill_distiller("Relative\nabundance", palette = "RdYlBu") + theme_bw() p.heat <- p.heat + scale_fill_distiller(legend, palette = "RdYlBu") + theme_bw()
# Make bacterial names italics # Make bacterial names italics
p.heat <- p.heat + theme(axis.text.y = element_text(colour = 'black', p.heat <- p.heat + theme(axis.text.y = element_text(colour = 'black',
...@@ -57,5 +82,8 @@ heatmap_fun <- function(data = data, column1 = "", top = 20, output = "./plot_he ...@@ -57,5 +82,8 @@ heatmap_fun <- function(data = data, column1 = "", top = 20, output = "./plot_he
flog.info(paste(output1,"heatmap_",column1,".html",sep='')) flog.info(paste(output1,"heatmap_",column1,".html",sep=''))
saveWidget(pltly.heat, file= paste(output1,"heatmap_",column1,".html",sep='')) saveWidget(pltly.heat, file= paste(output1,"heatmap_",column1,".html",sep=''))
return(p.heat)
LL$plot <- p.heat
LL$table <- data.com
return(LL)
} }
...@@ -12,7 +12,9 @@ bars_fun( ...@@ -12,7 +12,9 @@ bars_fun(
Fact1 = NULL, Fact1 = NULL,
split = FALSE, split = FALSE,
relative = TRUE, relative = TRUE,
outfile = "plot_compo.html" ylab = "Abundance",
outfile = "plot_compo.html",
verbose = TRUE
) )
} }
\arguments{ \arguments{
...@@ -30,6 +32,8 @@ bars_fun( ...@@ -30,6 +32,8 @@ bars_fun(
\item{relative}{Plot relative (TRUE, default) or raw abundance plot (FALSE)} \item{relative}{Plot relative (TRUE, default) or raw abundance plot (FALSE)}
\item{ylab}{Y axis title ("Abundance")}
\item{outfile}{Output html file.} \item{outfile}{Output html file.}
} }
\value{ \value{
......
...@@ -9,7 +9,10 @@ heatmap_fun( ...@@ -9,7 +9,10 @@ heatmap_fun(
column1 = "", column1 = "",
top = 20, top = 20,
output = "./plot_heatmap/", output = "./plot_heatmap/",
rank = "Species" rank = "Species",
relative = TRUE,
aggregate = TRUE,
legend = "Abundance"
) )
} }
\arguments{ \arguments{
...@@ -22,6 +25,12 @@ heatmap_fun( ...@@ -22,6 +25,12 @@ heatmap_fun(
\item{output}{Output directory} \item{output}{Output directory}
\item{rank}{Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)} \item{rank}{Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)}
\item{relative}{If TRUE, calculate relative abundance, needs a phyloseq object with raw abundance.}
\item{aggregate}{If TRUE, aggregate abundance of non "top" taxa in a new taxa named "Other".}
\item{legend}{Legend title ("Abundance")}
} }
\value{ \value{
Generate an heatmap with top taxa (html plotly version in output directory) Generate an heatmap with top taxa (html plotly version in output directory)
......
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