.Rhistory

q()
library(DESeq2)
q()
5 + 2
## Load required libraries
library(DESeq2)
knitr::opts_chunk$set(echo = T)
knitr::opts_chunk$set(message=F)
knitr::opts_chunk$set(warning=F)
knitr::opts_chunk$set(comment = NA)
knitr::opts_chunk$set(fig.width=9, fig.height=9)
knitr::opts_chunk$set(error = TRUE)
## Load required libraries
library(DESeq2)
library(data.table)
library(tidyverse)
library(vegan)
library(lmPerm)
library(viridis)
library(grid)
library(gridExtra)
library(cowplot)
library(devtools)
install_github("eastmallingresearch/Metabarcoding_pipeline/scripts")
library(metafuncs)
TAXCONF = 0.6   # Sets the taxonomy confidence level to get "rank" in taxonomy files
TOPOTU = 10    # Number of Top OTUs for summary information
# Run constants
Factor1 = "trial"
Factor2 = "cultivar"
design = x ~ trial / block + planting_date * cultivar
# colour blind palette
cbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")
## Custom functions
gfunc <- function(countData,colData,title) {
#### Rarefaction curve plotter ####
colData <- colData[names(countData),]
# descending order each sample
DT <- data.table(apply(countData,2,sort,decreasing=T))
# get cummulative sum of each sample
DT <- cumsum(DT)
# log the count values
DT <- log10(DT)
# set values larger than maximum for each column to NA
DT <- data.table(apply(DT,2,function(x) {x[(which.max(x)+1):length(x)]<- NA;x}))
# remove rows with all NA
DT <- DT[rowSums(is.na(DT)) != ncol(DT), ]
# add a count column to the data table
DT$x <- seq(1,nrow(DT))
# melt the data table for easy plotting
MDT <- melt(DT,id.vars="x")
# create an empty ggplot object from the data table
g <- ggplot(data=MDT,aes(x=x,y=value,colour=variable))
# remove plot background and etc.
g <- g + theme_classic_thin() %+replace%
theme(legend.position="none",axis.title=element_blank())
# plot cumulative reads
g <- g + geom_line(size=1.5) + scale_colour_viridis(discrete=T)
# add axis lables
g <- g + ggtitle(title)
#g <- g + ylab(expression("Log"[10]*" aligned sequenecs"))+xlab("OTU count")
# print the plot
g
}
metadata <- "sample_metadata.txt"
# Load
#ubiome_BAC <- loadData("data/BAC.otu_table.txt",metadata,"data/BAC.sintax.taxa",RHB="BAC")
#ubiome_FUN <- loadData("data/FUN.otu_table.txt",metadata,"data/FUN.sintax.taxa",RHB="FUN")
ubiome_BAC <- loadData("data/BAC.zotu_table.txt",metadata,"data/zBAC.sintax.taxa",RHB="BAC")
ubiome_FUN <- loadData("data/FUN.zotu_table.txt",metadata,"data/zFUN.sintax.taxa",RHB="FUN")
# Filter Plant, Chloroplast, and Eukaryote OTUs
# Fungi: Plantae OTUs
cat("Fungi:", length(grep("Plantae", ubiome_FUN$taxData$kingdom)), "Plantae OTUs\n")
# Bacteria: Chloroplast (Streptophyta) and Eukaryote OTUs
cat(
"Bacteria:", length(grep("Streptophyta",ubiome_BAC$taxData$genus)), "Chloroplast OTUs;",
length(grep("Eukaryota",ubiome_BAC$taxData$kingdom)), "Eukaryote OTUs\n"
)
# Filter Chloroplast and Eukaryote
filt <- rownames(
ubiome_BAC$taxData[
grepl("Streptophyta", ubiome_BAC$taxData$genus) &
as.numeric(ubiome_BAC$taxData$g_conf) >= TAXCONF,
]
)
filt <- c(filt, rownames(ubiome_BAC$taxData[grep("Eukaryota", ubiome_BAC$taxData$kingdom), ]))
cat("Bacteria: removing",length(filt),"OTUs")
ubiome_BAC$taxData <- ubiome_BAC$taxData[!rownames(ubiome_BAC$taxData) %in% filt,]
ubiome_BAC$countData <- ubiome_BAC$countData[!rownames(ubiome_BAC$countData) %in% filt,]
# filtered for minimum of 1000 reads
ubiome_FUN$dds <- ubiom_to_des(ubiome_FUN,filter=expression(colSums(countData)>=1000))
ubiome_BAC$dds <- ubiom_to_des(ubiome_BAC,filter=expression(colSums(countData)>=1000))
invisible(mapply(assign, names(ubiome_BAC), ubiome_BAC, MoreArgs = list(envir = globalenv())))
rare_bac <- gfunc(as.data.frame(counts(dds)), as.data.frame(colData(dds)), "Bacteria ZOTU")
invisible(mapply(assign, names(ubiome_FUN), ubiome_FUN, MoreArgs = list(envir = globalenv())))
rare_fun <- gfunc(as.data.frame(counts(dds)), as.data.frame(colData(dds)), "Fungi ZOTU")
rarefaction_plots <- grid.arrange(
rare_bac, rare_fun,
left = textGrob(label = expression("Log"[10] * " aligned sequenecs"), rot = 90),
bottom = "ZOTU count", nrow = 2
)
ggsave(filename = "rarefaction_plots.png", plot = rarefaction_plots, path = "figures/")
rarefaction_plots
View(rarefaction_plots)
globalenv()
# Fungi
invisible(mapply(assign,names(ubiome_FUN),ubiome_FUN,MoreArgs=list(envir=globalenv())))
print(
"Raw reads","\n\n",
"Total raw reads:\t\t", sum(countData),"\n",
"Mean raw reads per sample:\t", mean(colSums(countData)),"\n",
"Median raw reads per sample:", median(colSums(countData)),"\n",
"Max raw reads per sample:\t", max(colSums(countData)),"\n\n",
"Min raw reads per sample:\t", min(colSums(countData)),"\n"
)
cat(
"Raw reads","\n\n",
"Total raw reads:\t\t", sum(countData),"\n",
"Mean raw reads per sample:\t", mean(colSums(countData)),"\n",
"Median raw reads per sample:", median(colSums(countData)),"\n",
"Max raw reads per sample:\t", max(colSums(countData)),"\n\n",
"Min raw reads per sample:\t", min(colSums(countData)),"\n"
)
cat(
"Raw reads","\n\n",
"Total raw reads:\t\t", sum(countData),"\n",
"Mean raw reads per sample:\t", mean(colSums(countData)),"\n",
"Median raw reads per sample:\t", median(colSums(countData)),"\n",
"Max raw reads per sample:\t", max(colSums(countData)),"\n",
"Min raw reads per sample:\t", min(colSums(countData)),"\n"
)
#colSums(countData)
nct <- counts(dds,normalize=T)
cat("Normalised reads","\n\n",
"Total normalised reads:\t\t", sum(nct),"\n",
"Mean normalised reads per sample:\t", mean(colSums(nct)),"\n",
"Median normalised reads per sample:\t", median(colSums(nct)),"\n",
"Min normalised reads per sample:\t", min(colSums(nct)),"\n",
"Max normalised reads per sample:\t", max(colSums(nct)),"\n\n"
)
colSums(countData)
round(colSums(counts(dds,normalize=T)),0)
cat(
"Total OTUs:\t\t", nrow(taxData),"\n\n",
"Raw reads per OTU summary", "\n\n",
"Mean raw reads per OTU:\t", mean(rowSums(countData)),"\n",
"Median raw per OTU:\t", median(rowSums(countData)),"\n",
"OTU raw Min reads:\t\t", min(rowSums(countData)),"\n",
"OTU raw Max reads:\t\t", max(rowSums(countData)),"\n\n"
)
cat(
"Total OTUs:\t\t", nrow(taxData),"\n\n",
"Raw reads per OTU summary", "\n\n",
"Mean raw reads per OTU:\t", mean(rowSums(countData)),"\n",
"Median raw per OTU:\t\t", median(rowSums(countData)),"\n",
"OTU raw Min reads:\t\t", min(rowSums(countData)),"\n",
"OTU raw Max reads:\t\t", max(rowSums(countData)),"\n\n"
)
cat(
"Normalised reads per OTU summary","\n\n",
"Mean normalised reads per OTU:\t", mean(rowSums(nct)),"\n",
"Median normalised reads per OTU:\t", median(rowSums(nct)),"\n",
"OTU normalised Min reads:\t\t", min(rowSums(nct)),"\n",
"OTU normalised Max reads:\t\t", max(rowSums(nct)),"\n\n"
)
cat(
"Normalised reads per OTU summary","\n\n",
"Mean normalised reads per OTU:\t\t", mean(rowSums(nct)),"\n",
"Median normalised reads per OTU:\t", median(rowSums(nct)),"\n",
"OTU normalised Min reads:\t\t", min(rowSums(nct)),"\n",
"OTU normalised Max reads:\t\t", max(rowSums(nct)),"\n\n"
)
cat(
"Total OTUs:\t\t", nrow(taxData),"\n\n",
"Raw reads per OTU summary", "\n\n",
"Mean raw reads per OTU:\t", mean(rowSums(countData)),"\n",
"Median raw per OTU:\t\t", median(rowSums(countData)),"\n",
"OTU raw Min reads:\t\t", min(rowSums(countData)),"\n",
"OTU raw Max reads:\t\t", max(rowSums(countData)),"\n\n"
)
cat(
"Normalised reads per OTU summary","\n\n",
"Mean normalised reads per OTU:\t\t", mean(rowSums(nct)),"\n",
"Median normalised reads per OTU:\t", median(rowSums(nct)),"\n",
"OTU normalised Min reads:\t\t", min(rowSums(nct)),"\n",
"OTU normalised Max reads:\t\t", max(rowSums(nct)),"\n\n"
)
y <- rowSums(nct)
y <- y[order(y,decreasing = T)]
# proportion
xy <- y/sum(y)
cat("Total " ,TOPOTU,"OTUs:\n")
data.frame(counts=y[1:TOPOTU],Prop=xy[1:TOPOTU],rank=taxData[names(y)[1:TOPOTU],]$rank)
cat(
"Total OTUs:\t\t", nrow(taxData),"\n\n",
"Raw reads per OTU summary", "\n\n",
"Mean raw reads per OTU:\t", mean(rowSums(countData)),"\n",
"Median raw per OTU:\t\t", median(rowSums(countData)),"\n",
"OTU raw Min reads:\t\t", min(rowSums(countData)),"\n",
"OTU raw Max reads:\t\t", max(rowSums(countData)),"\n\n"
)
cat(
"Normalised reads per OTU summary","\n\n",
"Mean normalised reads per OTU:\t\t", mean(rowSums(nct)),"\n",
"Median normalised reads per OTU:\t", median(rowSums(nct)),"\n",
"OTU normalised Min reads:\t\t", min(rowSums(nct)),"\n",
"OTU normalised Max reads:\t\t", max(rowSums(nct)),"\n\n"
)
y <- rowSums(nct)
y <- y[order(y,decreasing = T)]
# proportion
xy <- y/sum(y)
cat("Top " ,TOPOTU, "OTUs:\n")
data.frame(counts=y[1:TOPOTU],Prop=xy[1:TOPOTU],rank=taxData[names(y)[1:TOPOTU],]$rank)
# Proportion of OTUs which can be assigned (with the given confidence) at each taxonomic rank
tx <- copy(taxData)
setDT(tx)
cols <- names(tx)[9:15]
tx[,(cols):=lapply(.SD,as.factor),.SDcols=cols]
data.table(
rank=c("kingdom","phylum","class","order","family","genus","species"),
"0.8"=round(unlist(lapply(cols,function(col) sum(as.number(tx[[col]])>=0.8)/nrow(tx))),2),
"0.65"=round(unlist(lapply(cols,function(col) sum(as.number(tx[[col]])>=0.65)/nrow(tx))),2),
"0.5"=round(unlist(lapply(cols,function(col) sum(as.number(tx[[col]])>=0.5)/nrow(tx))),2)
)
tx <-taxData[rownames(dds),]
nc <- counts(dds,normalize=T)
ac <- sum(nc)
data.table(
rank=c("kingdom","phylum","class","order","family","genus","species"),
"0.8"=round(unlist(lapply(cols ,function(col)(sum(nc[which(as.numeric(tx[[col]])>=0.8),])/ac *100))),2),
"0.65"=round(unlist(lapply(cols,function(col)(sum(nc[which(as.numeric(tx[[col]])>=0.65),])/ac *100))),2),
"0.5"=round(unlist(lapply(cols,function(col)(sum(nc[which(as.numeric(tx[[col]])>=0.5),])/ac *100))),2)
)
design <- x ~ trial / block + planting_date * cultivar # + Error(Plot)
# get the diversity index data
all_alpha_ord <- plot_alpha(counts(dds,normalize=T),colData(dds),design="trial",returnData=T)
# join diversity indices and metadata
all_alpha_ord <- all_alpha_ord[
as.data.table(colData(dds), keep.rownames = "Samples"),
on = "Samples"
]
all_alpha_ord[, trial.block := as.factor(paste0(trial, block))]
design <- x ~ trial + trial.block + planting_date * cultivar
design <- x ~ trial / block + planting_date * cultivar # + Error(Plot)
setkey(all_alpha_ord, S.chao1)
all_alpha_ord[, measure := as.numeric(as.factor(S.chao1))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
all_alpha_ord[, trial.block := as.factor(paste0(trial, block))]
design <- x ~ trial + trial.block + planting_date * cultivar
setkey(all_alpha_ord, S.chao1)
all_alpha_ord[, measure := as.numeric(as.factor(S.chao1))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
design <- x ~ trial / trial.block + planting_date * cultivar
setkey(all_alpha_ord, S.chao1)
all_alpha_ord[, measure := as.numeric(as.factor(S.chao1))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
setkey(all_alpha_ord,shannon)
all_alpha_ord[,measure:=as.numeric(as.factor(shannon))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
design <- x ~ trial / trial.block + trial *planting_date * cultivar
setkey(all_alpha_ord,shannon)
all_alpha_ord[,measure:=as.numeric(as.factor(shannon))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
design <- x ~ trial / block + trial *planting_date * cultivar
setkey(all_alpha_ord,shannon)
all_alpha_ord[,measure:=as.numeric(as.factor(shannon))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
colData(dds)
View(colData(dds))
all_alpha_ord
# get the diversity index data
all_alpha_ord <- plot_alpha(counts(dds,normalize=T),colData(dds),design="trial",returnData=T)
# join diversity indices and metadata
all_alpha_ord <- all_alpha_ord[
as.data.table(colData(dds), keep.rownames = "Samples"),
on = "Samples"
]
all_alpha_ord[, trial.block := as.factor(paste0(trial, block))]
design <- x ~ trial / trial.block + trial *planting_date * cultivar
setkey(all_alpha_ord, S.chao1)
all_alpha_ord[, measure := as.numeric(as.factor(S.chao1))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
setkey(all_alpha_ord,shannon)
all_alpha_ord[,measure:=as.numeric(as.factor(shannon))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
design <- x ~ + trial * planting_date * cultivar + trial:trial.block
setkey(all_alpha_ord,shannon)
all_alpha_ord[,measure:=as.numeric(as.factor(shannon))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
# get the diversity index data
all_alpha_ord <- plot_alpha(counts(dds,normalize=T),colData(dds),design="trial",returnData=T)
# join diversity indices and metadata
all_alpha_ord <- all_alpha_ord[
as.data.table(colData(dds), keep.rownames = "Samples"),
on = "Samples"
]
all_alpha_ord[, trial.block := as.factor(paste0(trial, block))]
design <- x ~ + trial * planting_date * cultivar + trial / trial.block
setkey(all_alpha_ord, S.chao1)
all_alpha_ord[, measure := as.numeric(as.factor(S.chao1))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
setkey(all_alpha_ord,shannon)
all_alpha_ord[,measure:=as.numeric(as.factor(shannon))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
setkey(all_alpha_ord,simpson)
all_alpha_ord[,measure:=as.numeric(as.factor(shannon))]
summary(aovp(update(design,measure~.),all_alpha_ord,seqs=T))
dds <- dds[rowSums(counts(dds, normalize=T))>5,]
### PCA ###
# perform PC decomposition of DES object
mypca <- des_to_pca(dds)
# to get pca plot axis into the same scale create a dataframe of PC scores multiplied by their variance
d <-t(data.frame(t(mypca$x)*mypca$percentVar))
round(mypca$percentVar[1:4],3)
apply(mypca$x[,1:4],2,function(x){
summary(aov(update(design,x~.),data=as.data.frame(cbind(x,colData(dds)))))
})
q()
.vsc.attach()
## Load required libraries
library(DESeq2)
library(data.table)
library(tidyverse)
library(vegan)
library(lmPerm)
library(viridis)
library(grid)
library(gridExtra)
library(cowplot)
library(devtools)
install_github("eastmallingresearch/Metabarcoding_pipeline/scripts")
library(metafuncs)
q()
# library(car)
library(cowplot)
library(data.table)
library(DESeq2)
library(DHARMa)
library(ggpubr)
library(grid)
library(gridExtra)
library(iNEXT)
library(kableExtra)
library(knitr)
library(lmPerm)
library(MASS)
library(pscl)
# library(rcompanion)
library(seqinr)
library(tidyverse)
library(vegan)
library(viridis)
# devtools::install_github("eastmallingresearch/Metabarcoding_pipeline/scripts")
library(metafuncs)
load("FUN.RData")
significant_asvs[, c("Site", "ASV", "Taxonomy", "Abundance", "coef", "var", "p_adjusted")] %>%
kbl("html", digits = 3) %>%
kable_styling("striped", full_width = T) %>%
column_spec(3, extra_css = "word-wrap: break-word;") %>%
kable_save("tables/FUN_canker_site_asvs.html")
significant_asvs[, c("Site", "ASV", "Taxonomy", "Abundance", "coef", "var", "p_adjusted")] %>%
kbl("html", digits = 3) %>%
kable_styling("striped", full_width = T) %>%
column_spec(3, extra_css = "word-wrap: break-word;") %>%
save_kable("tables/FUN_canker_site_asvs.html")
significant_asvs[, c("Site", "ASV", "Taxonomy", "Abundance", "coef", "var", "p_adjusted")] %>%
kbl("html", digits = 3) %>%
kable_styling("striped", full_width = T) %>%
column_spec(3, extra_css = "word-wrap: break-word; width: 20%;") %>%
save_kable("tables/FUN_canker_site_asvs.html")
significant_asvs[, c("Site", "ASV", "Taxonomy", "Abundance", "coef", "var", "p_adjusted")] %>%
kbl("html", digits = 3) %>%
kable_styling("striped", full_width = T) %>%
column_spec(3, extra_css = "word-wrap: break-word; width: 200%;") %>%
save_kable("tables/FUN_canker_site_asvs.html")
significant_asvs[, c("Site", "ASV", "Taxonomy", "Abundance", "coef", "var", "p_adjusted")] %>%
kbl("html", digits = 3) %>%
kable_styling("striped", full_width = T) %>%
column_spec(3, extra_css = "word-wrap: break-all; width: 200%;") %>%
save_kable("tables/FUN_canker_site_asvs.html")
significant_asvs[, c("Site", "ASV", "Taxonomy", "Abundance", "coef", "var", "p_adjusted")] %>%
kbl("html", digits = 3) %>%
kable_styling("striped", full_width = T) %>%
column_spec(3, extra_css = "word-wrap: break-all; width: 10cm;") %>%
save_kable("tables/FUN_canker_site_asvs.html")
significant_asvs[, c("Site", "ASV", "Taxonomy", "Abundance", "coef", "var", "p_adjusted")] %>%
kbl("html", digits = 3) %>%
kable_styling("striped", full_width = T) %>%
column_spec(3, extra_css = "word-wrap: break-all; white-space: normal; width: 10cm;") %>%
save_kable("tables/FUN_canker_site_asvs.html")
significant_asvs[, c("Site", "ASV", "Taxonomy", "Abundance", "coef", "var", "p_adjusted")] %>%
kbl("html", digits = 3) %>%
kable_styling("striped", full_width = F) %>%
column_spec(3, extra_css = "word-wrap: break-all; white-space: normal; width: 10cm;") %>%
save_kable("tables/FUN_canker_site_asvs.html")
significant_asvs$Taxonomy
significant_asvs$Taxonomy %>% unlist
load("BAC.RData")
alpha_box <- ggboxplot(
data = alpha_combined, x = "Site", y = "value",
color = "Storage", add = "jitter", facet.by = c("measure", "kingdom"),
palette = cbPalette, scales = "free_y", legend = "bottom",
ylab = "Mean diversity index", xlab = "Site"
) #+ guides(color = guide_legend(position = "right"))
alpha_combined <- rbind(fun_alpha, bac_alpha) %>%
subset(select = c("Site", "Storage", "Scion", "Target", "S.chao1", "shannon", "simpson")) %>%
mutate(kingdom = ifelse(Target == "ITS", "Fungi", "Bacteria")) %>%
mutate(kingdom = factor(kingdom, levels = c("Fungi", "Bacteria"))) %>%
mutate(Storage = factor(Storage, levels = c("no", "yes"))) %>%
rename(Chao1 = S.chao1, Shannon = shannon, Simpson = simpson) %>%
pivot_longer(
cols = c("Chao1", "Shannon", "Simpson"), names_to = "measure", values_to = "value"
)
alpha_box <- ggboxplot(
data = alpha_combined, x = "Site", y = "value",
color = "Storage", add = "jitter", facet.by = c("measure", "kingdom"),
palette = cbPalette, scales = "free_y", legend = "bottom",
ylab = "Mean diversity index", xlab = "Site"
) #+ guides(color = guide_legend(position = "right"))
ggsave(
filename = "alpha_box.png", plot = alpha_box, path = "figures/",
height = 24, width = 24, units = "cm"
)
alpha_box_fungi <- ggboxplot(
data = alpha_combined[alpha_combined$kingdom == "Fungi", ], x = "Site", y = "value",
color = "Storage", add = "jitter", facet.by = "measure",
palette = cbPalette, scales = "free_y", legend = "bottom",
ylab = "Mean diversity index", xlab = "Site"
)
ggsave(
filename = "alpha_box_fungi.png", plot = alpha_box_fungi, path = "figures/",
height = 12, width = 24, units = "cm"
)
alpha_box_bacteria <- ggboxplot(
data = alpha_combined[alpha_combined$kingdom == "Bacteria", ], x = "Site", y = "value",
color = "Storage", add = "jitter", facet.by = "measure",
palette = cbPalette, scales = "free_y", legend = "bottom",
ylab = "Mean diversity index", xlab = "Site"
)
ggsave(
filename = "alpha_box_bacteria.png", plot = alpha_box_bacteria, path = "figures/",
height = 12, width = 24, units = "cm"
)
alpha_box_combined <- ggarrange(
alpha_box_fungi, alpha_box_bacteria,
ncol = 1, nrow = 2, labels = c("A", "B"),
common.legend = T, legend = "bottom"
)
ggsave(
filename = "alpha_box_combined.png", plot = alpha_box_combined, path = "figures/",
height = 24, width = 24, units = "cm"
)
ggsave(
filename = "alpha_box_combined.png", plot = alpha_box_combined, path = "figures/",
height = 20, width = 24, units = "cm"
)
alpha_box <- ggboxplot(
data = alpha_combined, x = "Site", y = "value",
color = "Storage", add = "jitter", facet.by = c("kingdom", "measure"),
palette = cbPalette, scales = "free_y", legend = "bottom",
ylab = "Mean diversity index", xlab = "Site"
) #+ guides(color = guide_legend(position = "right"))
ggsave(
filename = "alpha_box.png", plot = alpha_box, path = "figures/",
height = 24, width = 24, units = "cm"
)
alpha_box <- ggboxplot(
data = alpha_combined, x = "Site", y = "value",
color = "Storage", add = "jitter", facet.by = c("kingdom", "measure"),
palette = cbPalette, scales = "free", legend = "bottom",
ylab = "Mean diversity index", xlab = "Site"
) #+ guides(color = guide_legend(position = "right"))
ggsave(
filename = "alpha_box.png", plot = alpha_box, path = "figures/",
height = 24, width = 24, units = "cm"
)
alpha_box <- ggboxplot(
data = alpha_combined, x = "Site", y = "value",
color = "Storage", add = "jitter", facet.by = c("measure", "kingdom"),
palette = cbPalette, scales = "free_y", legend = "bottom",
ylab = "Mean diversity index", xlab = "Site"
) #+ guides(color = guide_legend(position = "right"))
ggsave(
filename = "alpha_box.png", plot = alpha_box, path = "figures/",
height = 24, width = 24, units = "cm"
)
q()