Commit 6f783900 authored by Etienne Rifa's avatar Etienne Rifa
Browse files

more outputs + infos for beta + misc changes

parent 56dde1e4
Pipeline #39411 passed with stage
in 15 seconds
......@@ -5,7 +5,7 @@
#' @param psobj a phyloseq object (output from decontam or generate_phyloseq)
#' @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 col A metadata column (among sample_variables(data)).
#' @param cov Covariable names comma separated vector.
#' @param cov Covariable names comma separated vector (last covariable is used to plot samples with different shape).
#' @param dist0 Dissimilarity index, partial match to "unifrac", "wunifrac", "dpcoa", "jsd", "manhattan", "euclidean", "canberra", "bray", "kulczynski", "jaccard", "gower", "altGower", "morisita", "horn", "mountford", "raup" , "binomial", "chao", "cao" or "mahalanobis".
#' @param ord0 Currently supported method options are: c("DCA", "CCA", "RDA", "CAP", "DPCoA", "NMDS", "MDS", "PCoA")
#' @param output The output file directory.
......@@ -34,50 +34,48 @@ diversity_beta_light <- function(psobj, rank = "ASV", col = NULL, cov = NULL, di
flog.debug(cov)
if(!is.null(cov)){
cov1 = unlist(strsplit(cov, ","))
cov1 <- unlist(strsplit(cov, ","))
}
# metrics <- sapply(strsplit(measures,","), '[')
col1 = unlist(strsplit(col, "[+]"))
col1 <- col
if(length(col1)==1){
fun <- paste("psobj = subset_samples(psobj, !is.na(",col,"))",sep="")
eval(parse(text=fun))
}else{
fun <- paste("psobj = subset_samples(psobj, !is.na(",col1[1],"))",sep="")
eval(parse(text=fun))
fun <- paste("psobj = subset_samples(psobj, !is.na(",col1[2],"))",sep="")
eval(parse(text=fun))
}
stop("Only one factor in 'col' argument, add more covariable in 'cov' argument.")
}
if(rank=="ASV"){
flog.info('No glom ...')
data_rank = psobj
data_rank <- psobj
}else{
data_rank = tax_glom(psobj, rank)
data_rank <- tax_glom(psobj, rank)
}
# Figure
resBeta = list()
flog.info('Plot ...')
resBeta <- list()
if(!is.null(cov)){
p1 <- plot_samples(data_rank, ordinate(data_rank, ord0, dist0), color = col, shape = cov1[1] ) +
p1 <- plot_samples(data_rank, ordinate(data_rank, ord0, dist0), color = col, shape = cov1[length(cov1)] ) +
theme_bw() + ggtitle(glue::glue("{ord0} + {dist0}")) + stat_ellipse() + scale_shape_manual(values = 0:10)
}else{
p1 <- plot_samples(data_rank, ordinate(data_rank, ord0, dist0), color = col) +
theme_bw() + ggtitle(glue::glue("{ord0} + {dist0}")) + stat_ellipse()
}
# plot(p1)
resBeta$plot = p1 + theme(axis.text.x = element_text(angle = 45, hjust=1),
flog.info('Plot ok...')
resBeta$plot <- p1 + theme(axis.text.x = element_text(angle = 45, hjust=1),
,axis.text=element_text(size=18),
axis.title=element_text(size=16,face="bold"),
strip.text.x = element_text(size = 18,face="bold"),
title=element_text(size=16,face="bold"))
if(tests){
otable = otu_table(data_rank)
otable <- otu_table(data_rank)
flog.info(glue::glue('Tests on {dist0} ...'))
mdata = data.frame(sample_data(data_rank))
mdata <- data.frame(sample_data(data_rank))
mdata$Depth <- sample_sums(data_rank)
if( any(grepl(dist0, c("unifrac", "wunifrac", "dpcoa", "jsd") )) ){
dist1 <-phyloseq::distance(data_rank, dist0)
......@@ -85,18 +83,21 @@ diversity_beta_light <- function(psobj, rank = "ASV", col = NULL, cov = NULL, di
dist1 <<- vegdist(t(otable), method = dist0)
}
if(!is.null(cov)){
resBC = adonis(as.formula(paste('dist1 ~ Depth +', paste(cov1, collapse="+"), "+", col)), data = mdata, permutations = 1000)
form1 <- as.formula(paste('dist1 ~ Depth +', paste(cov1, collapse="+"), "+", col))
resBC <- adonis(form1, data = mdata, permutations = 1000)
}else{
resBC = adonis(as.formula(paste('dist1 ~ Depth +', col)), data = mdata, permutations = 1000)
form1 <- as.formula(paste('dist1 ~ Depth +', col))
resBC <- adonis(form1, data = mdata, permutations = 1000)
}
#PairwiseAdonis
# resBC2 = pairwise.adonis2(as.formula( paste('BC.dist ~ ', col,sep="") ), data = mdata)
if(length(col1)>1){
fact1 <- apply( mdata[,c(col1)] , 1 , paste , collapse = "-" )
resBC2 = pairwise.adonis(dist1, fact1, p.adjust.m='fdr')
resBC2 <- pairwise.adonis(dist1, fact1, p.adjust.m='fdr')
} else {
resBC2 = pairwise.adonis(dist1, mdata[,c(col1)], p.adjust.m='fdr')
resBC2 <- pairwise.adonis(dist1, mdata[,c(col1)], p.adjust.m='fdr')
}
write.table(resBC$aov.tab, file=paste0(output,'/',col,'_permANOVA.txt'), sep="\t")
......@@ -104,9 +105,12 @@ diversity_beta_light <- function(psobj, rank = "ASV", col = NULL, cov = NULL, di
ggsave(glue::glue("{output}/beta_diversity.eps"), plot=resBeta$plot, height = 20, width = 30, units="cm", dpi = 500, device="eps")
resBeta$permanova = resBC$aov.tab
resBeta$pairwisepermanova = resBC2
resBeta$permanova <- resBC$aov.tab
resBeta$permanova_formula <- format(form1)
resBeta$pairwisepermanova <- resBC2
resBeta$test_table <- mdata
resBeta$dist <- dist1
}
return(resBeta)
}
}
\ No newline at end of file
......@@ -20,7 +20,7 @@ heatmap_fun <- function(data = data, column1 = "", top = 20, output = "./plot_he
output1 <- paste(getwd(),'/',output,'/',sep='')
if(!dir.exists(output1)){
dir.create(output1)
dir.create(output1, recursive = TRUE)
}
flog.info('Done.')
......@@ -33,7 +33,7 @@ heatmap_fun <- function(data = data, column1 = "", top = 20, output = "./plot_he
colnames(data.com)
p.heat <- ggplot(data.com, aes(x = Sample, y = Tax)) + geom_tile(aes(fill = Abundance))
p.heat <- p.heat + scale_fill_distiller("Abundance", palette = "RdYlBu") + theme_bw()
p.heat <- p.heat + scale_fill_distiller("Relative\nabundance", palette = "RdYlBu") + theme_bw()
# Make bacterial names italics
p.heat <- p.heat + theme(axis.text.y = element_text(colour = 'black',
......
......@@ -12,7 +12,9 @@ ASVenn_fun(
subset = "",
lvls = NULL,
krona = "",
shared = TRUE
shared = TRUE,
verbose = 2,
ggplotmode = FALSE
)
}
\arguments{
......@@ -31,6 +33,10 @@ ASVenn_fun(
\item{krona}{Krona of exclusive ASV or shared with informed level and others. Must be among levels of column1 argument.}
\item{shared}{shared [TRUE] or exclusive [FALSE] mode.}
\item{verbose}{Verbose level. (1: quiet, 2: print infos, 3: print infos + debug)}
\item{ggplotmode}{if TRUE plot the Venn diagram using ggplot}
}
\value{
Returns list with venn diagram and table with shared features. Exports a venn diagram with corresponding tabulated file.
......
......@@ -10,7 +10,8 @@ VENNFUN(
TITRE = TITRE,
output = "./",
refseq1 = NULL,
alltax = NULL
alltax = NULL,
ggplotmode = FALSE
)
}
\arguments{
......@@ -25,6 +26,8 @@ VENNFUN(
\item{refseq1}{Reference sequences.}
\item{alltax}{Taxonomy table.}
\item{ggplotmode}{if TRUE plot the Venn diagram using ggplot}
}
\value{
Exports a venn diagram with corresponding tabulated file.
......
......@@ -22,7 +22,7 @@ diversity_beta_light(
\item{col}{A metadata column (among sample_variables(data)).}
\item{cov}{Covariable names comma separated vector.}
\item{cov}{Covariable names comma separated vector (last covariable is used to plot samples with different shape).}
\item{dist0}{Dissimilarity index, partial match to "unifrac", "wunifrac", "dpcoa", "jsd", "manhattan", "euclidean", "canberra", "bray", "kulczynski", "jaccard", "gower", "altGower", "morisita", "horn", "mountford", "raup" , "binomial", "chao", "cao" or "mahalanobis".}
......
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