Skip to content
Snippets Groups Projects
Commit df064dd4 authored by Simeon's avatar Simeon
Browse files

Formatting improvements for pdf knitting (FEMSEc submission)

Additional descriptions and formatting improvements to submit a knitted PDF instead of the Rmd as supplementary material.
parent 979f4475
No related branches found
No related tags found
No related merge requests found
colors_oomyc_genera <- colors_oomyc_genera[order(names(colors_oomyc_genera))]
colors_oomyc_genera["Phytophthora"] <- viridis_reds[1]
colors_oomyc_genera["Pythium"] <- viridis_reds[4]
#define ggplot theme (shared between genus and species level plots)
shared_theme <- theme(legend.key.size = unit(3, "mm"),
legend.position = "none",
axis.title.x = element_blank(),
plot.margin = margin(20,10,10,10),
legend.text.align = 0)
shared_theme_abu <- theme(legend.key.size = unit(3, "mm"),
legend.position = "none",
axis.title.x = element_blank(),
plot.margin = margin(20,10,10,10),
legend.text.align = 0)
#Format legend labels
if(customLabels){
labelling <- c(expression(paste("uncultured ", italic("Apodachlya"))),
expression(paste("uncultured ", italic("Lagenidium"))),
expression(paste("uncultured ", italic("Oomycetes"))),
expression(paste("uncultured ", italic("Phytophthora"))),
expression(paste("uncultured ", italic("Pythium"))),
expression(italic("Achlya"), italic("Albugo"), italic("Aphanomyces"),
italic("Apodachlya"), italic("Atkinsiella"), italic("Bremia"),
italic("Brevilegnia"), italic("Dictyuchus"), italic("Eurychasma"),
italic("Geolegnia"), italic("Globisporangium"), italic("Halocrusticida"),
italic("Haptoglossa"), italic("Hyaloperonospora"), italic("Lagena"),
italic("Lagenidium"), italic("Leptolegnia"), italic("Myzocytiopsis"), italic("Lagenidium"), italic("Leptolegnia"), italic("Myzocytiopsis"),
italic("Paralagenidium"), italic("Peronospora"), italic("Phragmosporangium"), italic("Paralagenidium"), italic("Peronospora"), italic("Phragmosporangium"),
italic("Phytophthora"), italic("Phytopythium"), italic("Pilasporangium"), italic("Phytophthora"), italic("Phytopythium"), italic("Pilasporangium"),
...@@ -7,506 +33,480 @@ italic("Saprolegnia"), italic("Thraustotheca"), italic("Wilsoniana")), NA) ...@@ -7,506 +33,480 @@ italic("Saprolegnia"), italic("Thraustotheca"), italic("Wilsoniana")), NA)
}else{ }else{
labelling <- unique(ASVs_funguild$Genus) labelling <- unique(ASVs_funguild$Genus)
} }
guild_g <- ggplot(ASVs_funguild, aes(x = guild_genus, fill = Genus)) + geom_bar() + shared_theme + guild_g <- ggplot(ASVs_funguild, aes(x = guild_genus, fill = Genus)) +
scale_fill_manual(values = colors_oomyc_genera, na.value = "grey", labels = labelling) + geom_bar() + shared_theme +
scale_fill_manual(values = colors_oomyc_genera,
na.value = "grey", labels = labelling) +
ylab("Number of ASVs") + ylim(0,1500) ylab("Number of ASVs") + ylim(0,1500)
guild_s <- ggplot(ASVs_funguild, aes(x = guild_species, fill = Genus)) + geom_bar() + shared_theme + guild_s <- ggplot(ASVs_funguild, aes(x = guild_species, fill = Genus)) +
scale_fill_manual(values = colors_oomyc_genera, na.value = "grey", labels = labelling) + geom_bar() + shared_theme +
scale_fill_manual(values = colors_oomyc_genera,
na.value = "grey", labels = labelling) +
ylab("Number of ASVs") + ylim(0,1500) ylab("Number of ASVs") + ylim(0,1500)
guild_g_abu <- ggplot(ASVs_funguild, aes(x = guild_genus, y = Abundance_sum/1000000, fill = Genus)) + guild_g_abu <- ggplot(ASVs_funguild, aes(x = guild_genus,
y = Abundance_sum/1000000,
fill = Genus)) +
geom_col() + shared_theme_abu + geom_col() + shared_theme_abu +
scale_fill_manual(values = colors_oomyc_genera, na.value = "grey", labels = labelling) + scale_fill_manual(values = colors_oomyc_genera,
na.value = "grey", labels = labelling) +
ylab("Abundance in M reads") + ylim(0,2.5) ylab("Abundance in M reads") + ylim(0,2.5)
guild_s_abu <- ggplot(ASVs_funguild, aes(x = guild_species, y = Abundance_sum/1000000, fill = Genus)) + guild_s_abu <- ggplot(ASVs_funguild, aes(x = guild_species,
y = Abundance_sum/1000000,
fill = Genus)) +
geom_col() + shared_theme_abu + geom_col() + shared_theme_abu +
scale_fill_manual(values = colors_oomyc_genera, na.value = "grey", labels = labelling) + scale_fill_manual(values = colors_oomyc_genera,
na.value = "grey", labels = labelling) +
ylab("Abundance in M reads") + ylim(0,2.5) ylab("Abundance in M reads") + ylim(0,2.5)
guilds <- ggarrange(guild_g, guild_g_abu, guild_s, guild_s_abu, nrow = 2, ncol = 2, guilds <- ggarrange(guild_g, guild_g_abu, guild_s,
guild_s_abu, nrow = 2, ncol = 2,
common.legend = TRUE, legend = "right" ,align = "hv", common.legend = TRUE, legend = "right" ,align = "hv",
labels = c("Genus level", "", "Species level", ""), font.label = list(size = 14)) labels = c("Genus level", "", "Species level", ""),
font.label = list(size = 14))
guilds guilds
save_plot(guilds ,plot_name = "FunGuild_combo") save_plot(guilds ,plot_name = "FunGuild_combo_Fig2")
cat("Chunk successfully run")
## CHANGE ME to change the width (in cm) of the output.
wp <- 17.5
# CHANGE ME to change the height (in cm) of the output.
hp <- 15
# CHANGE ME to change the resolution (in dpi) of the output.
res <- 150
# calculate averages for display in graphs
avgNames <- c("PPUntr", "PPBait","ChrUntr", "ChrBait",
"O-PUntr", "O-PBait", "OthUntr", "OthBait")
avgs <- c(mean(ranksUntrPP$Abundance),
mean(ranksBaitPP$Abundance),
mean(ranksUntrAlgae$Abundance),
mean(ranksBaitAlgae$Abundance),
mean(ranksUntrOwoP$Abundance),
mean(ranksBaitOwoP$Abundance),
mean(rankedOthersUntr$OtherASVcounts),
mean(rankedOthersBait$OtherASVcounts))
names(avgs) <- avgNames
# multiply to transform relative abundances into avg percentages and round to
# one decimal
avgs <- avgs*100
avgs <- format(round(avgs, 1), nsmall = 1)
# Define fixed color schemes for consistency between plots
plotvars <- c("Chrysophyceae", "All oomycetes",
"Other oomycetes", "Other classes", "PP")
# Just to catch others/exceptions, make a general palette that will be
# overwritten below
plotPalette <- viridis(length(plotvars))
names(plotPalette) <- plotvars
plotPalette["Other classes"] <- "#C0C0C0"
plotPalette["Chrysophyceae"] <- viridis_yellows[2]
plotPalette["PP"] <- viridis_reds[4]
plotPalette["Other oomycetes"] <- viridis_reds[1]
plotPalette["All oomycetes"] <- viridis_reds[1]
# plot total ASV counts, sorted by rank and color by class.
AllBait <- ggplot(data = OthersRankPP,
aes(x = Rank, y = Abundance, fill = Set))
AllUntr <- ggplot(data = UntrOthersRankPP,
aes(x = Rank, y = Abundance, fill = Set))
customTheme <- theme(legend.position = "none",
axis.title = element_blank(),
title = element_text(size = 6))
customThemeA <- theme(legend.position = "none",
legend.title = element_blank(),
legend.key.size = unit(5,"mm"),
axis.title.x = element_blank())
xUpperLim <- nrow(OthersRankPP) + 1
ytextl = 0.05
xtextl = 15
ytexts = 0.75
xtexts = 45
middleline <- geom_blank() #geom_hline(yintercept = 0.5,
# linetype="dashed", color="black")
pAllBait <- AllBait + geom_col() +
geom_col(data = AlgaeInvRankOom) +
geom_col(data = OomRankPP) +
geom_col(data = ranksBaitPP) +
scale_fill_manual(values = plotPalette) +
ggtitle("After enrichment") + xlab("Sample rank") +
labs(fill = "") + ylim(0,1) + xlim(0,xUpperLim) +
customThemeA + middleline +
annotate(geom = 'text', label = paste("Mean PP", avgs["PPBait"], "%"),
x = xtextl, y = ytextl) +
ylab("rel. Abundance")
pAllUntr <- AllUntr + geom_col() +
geom_col(data = UntrAlgaeInvRankOom) +
geom_col(data = UntrOomRankPP) +
geom_col(data = ranksUntrPP) +
scale_fill_manual(values = plotPalette) +
ggtitle("Before enrichment") + xlab("Sample rank") +
labs(fill = "") + ylim(0,1) + xlim(0,xUpperLim) +
customThemeA + middleline +
annotate(geom = 'text', label = paste("Mean PP", avgs["PPUntr"], "%"),
x = xtextl, y = ytextl) +
ylab("rel. Abundance")
AlgBait <- ggplot(data = ranksBaitAlgae,
aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Chrysophyceae - after enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["ChrBait"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
OomBait <- ggplot(data = ranksBaitOwoP,
aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Other oomycetes - after enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["O-PBait"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
OthBait <- ggplot(data = rankedOthersBait,
aes(x = Rank, y = OtherASVcounts, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Other classes - after enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["OthBait"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
AlgUntr <- ggplot(data = ranksUntrAlgae,
aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Chrysophyceae - before enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["ChrUntr"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
OomUntr <- ggplot(data = ranksUntrOwoP,
aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Other oomycetes - before enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["O-PUntr"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
OthUntr <- ggplot(data = rankedOthersUntr,
aes(x = Rank, y = OtherASVcounts, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Other classes - before enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["OthUntr"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
#A.combo <- ggarrange(pAllUntr, pAllBait, nrow=2,
# labels = c("A"), vjust = 1, common.legend = TRUE, legend = "bottom")
untr_smallplots <- ggarrange(AlgUntr, OomUntr, OthUntr, nrow = 3)
untr_smallplots_anot <- annotate_figure(untr_smallplots,
left = text_grob("Rel. abundance",
rot = 90, size = 11))
bait_smallplots <- ggarrange(AlgBait, OomBait, OthBait, nrow = 3)
bait_smallplots_anot <- annotate_figure(bait_smallplots,
left = text_grob("Rel. abundance",
rot = 90, size = 11))
combo <- ggarrange(pAllUntr, untr_smallplots_anot, pAllBait,
bait_smallplots_anot, nrow = 2, ncol = 2, widths = c(2,1),
labels = c("A", "", "B", ""), common.legend = TRUE,
legend = "bottom") + theme(plot.margin = margin(10,10,10,10))
#Save plot
save_plot(combo, plot_name = "Class_ranking_Fig3")
combo
cat("Chunk successfully run")
library(rstatix)
pairwise_test <- function(ASV_list){
ps.temp <- prune_taxa(ASV_list, ps.trans)
rankTemp <- cuphyr::make_ranked_sums(ps.temp)
rankTemp <- as.data.frame(rankTemp)
teststats <- pairwise_t_test(rankTemp, formula = Abundance ~ Bait)
return(teststats)}
ChrysUntrVsBait <- pairwise_test(subASVsAlg)
subASVsO_P <- setdiff(subASVsOom, subASVsPP)
OomUntrVsBait <- pairwise_test(subASVsO_P)
PPUntrVsBait <- pairwise_test(subASVsPP)
subASVsOth <- setdiff(taxa_names(ps), subASVsOom)
subASVsOth <- setdiff(subASVsOth, subASVsAlg)
OthUntrVsBait <- pairwise_test(subASVsOth)
all_stats <- list(PPUntrVsBait, ChrysUntrVsBait, OomUntrVsBait, OthUntrVsBait)
names(all_stats) <- c("PP", "Chrysophyceae", "Other oomycetes", "Other classes")
all_stats <- bind_rows(all_stats,.id = "Genera from")
all_stats
cat("Chunk successfully run") cat("Chunk successfully run")
topnt2$labels
topnt2$plot_env$fill
topnt2$plot_env$mdf
topnt2$plot_env$fill
topnt2$plot_env$mdf$Genus %>% unique()
topnt2$plot_env$mdf$Genus %>% unique() %>% sort()
names(colors_lvl2)
names(colors_lvl2) %>% sort()
#CHANGE ME to the amount of top taxlvl to plot (e.g. 'numt = 20' plots Top20 Taxa at taxlvl)
numt <- 5
#CHANGE ME to the amount of top taxa you want to plot (e.g. 'numt2 = 25' plots Top25 Taxa at taxlvl2)
numt2 <- 10
#CHANGE ME to the primary parameter of interest (x-axis). Accepted values are the column headers in the descriptor file.
primaryPar <- "Soil"
##CHANGE ME to the secondary parameter of interest (panels). Accepted values are the column headers in the descriptor file.
secondaryPar <- "Bait"
##CHANGE ME to the taxonomic level of interest (color coding). Accepted values are available taxonomic levels. 'taxlvl' should be higher than 'taxlvl2', e.g. Class and Genus
taxlvl <- "Class"
taxlvl2 <- "Genus"
#Adds formatted legend labels for the specific output in the manuscript. WARNING! Will overwrite actual values and needs to be turned off/adjusted manually when data is updated
customLabels <- TRUE
##CHANGE ME to change the width (in cm) of the output.
wp <- 20 wp <- 20
#CHANGE ME to change the height (in cm) of the output.
hp <- 15 hp <- 15
#CHANGE ME to change the resolution (in dpi) of the output. seqs <- readDNAStringSet(file = file.path(path,"iso_ASVs_refs_v3.fasta"),
res <- 300 format = "fasta")
#Get total number of reads and number of Pythium/Phytophthora reads (without contro samples!) to find out percentage of Pyhium and Phytophthora align <- AlignSeqs(DNAStringSet(seqs), iterations = 100, refinements = 5)
total_reads_no_contr <- ps_tbl %>% seqs_phang <- phangorn::phyDat(as(align, "matrix"), type = "DNA")
dplyr::filter(Pass == "Yes") %>% seqs_dm <- phangorn::dist.ml(seqs_phang)
dplyr::mutate(sumAbu = sum(Abundance)) %>% seqs_treeNJ <- NJ(seqs_dm)
select(sumAbu) %>% seqs_fit = pml(seqs_treeNJ, data = seqs_phang)
unique() %>% fitGTRseqs <- update(seqs_fit, k = 4, inv = 0.2)
unlist() %>% fitGTRseqs <- optim.pml(fitGTRseqs, model = "GTR", optInv = TRUE,
unname() optGamma = TRUE, rearrangement = "stochastic",
pythium_reads_no_contr <- ps_tbl %>% control = pml.control(trace = 0))
dplyr::filter(Pass == "Yes") %>% GTR_tree <- treeio::as.treedata(fitGTRseqs$tree, type = "ml")
dplyr::group_by(Genus) %>% iso_ASV_tree <- ggtree(GTR_tree) + geom_tree() +
summarise(sumAbu = sum(Abundance)) %>% geom_treescale(x = 0.5) +
dplyr::filter(Genus == "Pythium") %>% geom_tiplab() + xlim(-0.1,3)
select(sumAbu) %>% save_plot(iso_ASV_tree, plot_name = "IsolateSeqs_ASVs_ref_tree_Supfig4")
unlist() %>% iso_ASV_tree
unname()
phytoph_reads_no_contr <- ps_tbl %>%
dplyr::filter(Pass == "Yes") %>%
dplyr::group_by(Genus) %>%
summarise(sumAbu = sum(Abundance)) %>%
dplyr::filter(Genus == "Phytophthora") %>%
select(sumAbu) %>%
unlist() %>%
unname()
prctg_pyth_phyt_str <- paste0("\n\nTotal number of reads in the samples: ",
total_reads_no_contr ,
", number of Pythium reads: ",
pythium_reads_no_contr,
", number of Phytophthora reads: ",
phytoph_reads_no_contr,
"\nPercentage of Pythium total: ",
pythium_reads_no_contr/total_reads_no_contr*100,
"%\nPercentage of Phytophthora total: ",
phytoph_reads_no_contr/total_reads_no_contr*100, "%")
#Most abundant at specified levels, phyloseq objects
ps.topnt <- cuphyr::abundant_tax_physeq(
physeq = ps.trans, lvl = taxlvl, top = numt,
ignore_na = TRUE, silent = FALSE)
ps.topnt2 <- cuphyr::abundant_tax_physeq(
physeq = ps.trans,lvl = taxlvl2, top = numt2,
ignore_na = TRUE, silent = FALSE)
#Abundance lists (sorted by abundance)
taxlist_topnt1 <- cuphyr::abundant_tax_physeq(
physeq = ps.trans, lvl = taxlvl, top = numt,
output_format = "tops", ignore_na = TRUE)
taxlist_topnt2 <- cuphyr::abundant_tax_physeq(
physeq = ps.trans, lvl = taxlvl2, top = numt2,
output_format = "tops", ignore_na = TRUE)
#tidy evaluation translation magic (needed so the dplyr functions that extract 'tax_lookup' below can parse the strings from the variable)
lvl1 <- sym(taxlvl)
lvl2 <- sym(taxlvl2)
lvl1 <- enquo(lvl1)
lvl2 <- enquo(lvl2)
#make lookup table to guide coloring
tax_lookup <- ps_tbl %>%
group_by(!!lvl1, !!lvl2) %>%
summarise() %>%
filter(!!lvl2 %in% taxlist_topnt2)
#Define fixed color schemes for classes for consistency between plots. This is a color scheme specified for this particular published figure. Can be used for orientation in general cases.
#Background colors to be overwritten for special cases
colors_lvl1 <- viridis(length(taxlist_topnt1), option = "D", begin = 1 , end = 0.5)
names(colors_lvl1) <- taxlist_topnt1
colors_lvl2 <- viridis(length(taxlist_topnt2), option = "D", begin = 1 , end = 0.5)
names(colors_lvl2) <- taxlist_topnt2
colors_combo <- c(colors_lvl1, colors_lvl2)
#Special cases
colors_combo["Chrysophyceae"] <- viridis_yellows[2]
colors_combo["Pedospumella"] <- viridis_yellows[1]
colors_combo["'Spumella-like'"] <- viridis_light[1]
colors_combo["Oomycetes"] <- viridis_reds[4]
colors_combo["Phytophthora"] <- viridis_reds[1]
colors_combo["Pythium"] <- viridis_reds[4]
colors_combo["Aphanomyces"] <- viridis_blues[1]
colors_combo["Globisporangium"] <- "#4C0000FF"
colors_combo["Pythiogeton"] <- viridis_blues[2]
colors_combo["Saprolegnia"] <- viridis_blues[4]
#Format legend labels
if(customLabels){
labelling_lvl2 <- c(expression(paste(italic("Spumella"), "-like")),
"uncultured fungus", expression(italic("Aphanomyces"),
italic("Globisporangium"), italic("Haptoglossa"), italic("Pedospumella"),
italic("Phytophthora"), italic("Pythiogeton"), italic("Pythiopsis"),
italic("Pythium")))
}else{
labelling_lvl2 <- names(colors_lvl2) %>% sort()
}
#Re-label the facets
enrich_labs <- c("Before enrichment", "After enrichment")
names(enrich_labs) <- c("Before enrichment", "Enriched")
theme_bar <- theme(legend.position = "bottom",
legend.key.size = unit(0.4, "cm"),
legend.spacing.x = unit(0.3, 'cm'),
text = element_text(size = 10),
strip.text.x = element_text(size = 10),
strip.background = element_blank(),
axis.text.x = element_text(angle = 0, hjust = 0.5))
topnt <- plot_bar(ps.topnt, x = primaryPar, fill = taxlvl,
title = paste("Top", numt, "Classes")) +
facet_wrap(paste0("~", secondaryPar),
labeller = labeller(Bait = enrich_labs)) +
scale_x_discrete(breaks = c("S05" ,"S15", "S25", "S35", "S45" ,"S55")) +
scale_fill_manual(values = colors_combo) +
ylab("Relative abundance") +
theme_bar +
xlab("Soil samples S01-S64")
topnt_summary <- cuphyr::summarise_physeq(physeq = ps,
ASV_sublist = taxa_names(ps.topnt),
sublist_id =
paste0("top ", numt, " ", taxlvl),
samp_names = FALSE)
topnt2 <- plot_bar(ps.topnt2, x = primaryPar, fill = taxlvl2,
title = paste("Top", numt2, "Genera")) +
facet_wrap(paste0("~", secondaryPar),
labeller = labeller(Bait = enrich_labs)) +
scale_x_discrete(breaks = c("S05" ,"S15", "S25", "S35", "S45" ,"S55")) +
scale_fill_manual(values = colors_combo) +
ylab("Relative abundance") +
theme_bar + xlab("Soil samples S01-S64")
topnt2_summary <- cuphyr::summarise_physeq(physeq = ps,
ASV_sublist = taxa_names(ps.topnt2),
sublist_id =
paste0("top ", numt2, " ", taxlvl2),
samp_names = FALSE)
#Combine the plots
combobar1 <- ggarrange(topnt, topnt2, nrow = 2, labels = c("A", "B"), align = "hv")
#Save
save_plot(combobar1, plot_name = "Topn_barplot_Fig2")
#Print to standard out
cat(topnt_summary, "\n\n", topnt2_summary)
cat(prctg_pyth_phyt_str)
combobar1
cat("Chunk successfully run") cat("Chunk successfully run")
#CHANGE ME to the amount of top taxlvl to plot (e.g. 'numt = 20' plots Top20 Taxa at taxlvl)
numt <- 5
#CHANGE ME to the amount of top taxa you want to plot (e.g. 'numt2 = 25' plots Top25 Taxa at taxlvl2)
numt2 <- 10
#CHANGE ME to the primary parameter of interest (x-axis). Accepted values are the column headers in the descriptor file.
primaryPar <- "Soil"
##CHANGE ME to the secondary parameter of interest (panels). Accepted values are the column headers in the descriptor file.
secondaryPar <- "Bait"
##CHANGE ME to the taxonomic level of interest (color coding). Accepted values are available taxonomic levels. 'taxlvl' should be higher than 'taxlvl2', e.g. Class and Genus
taxlvl <- "Class"
taxlvl2 <- "Genus"
#Adds formatted legend labels for the specific output in the manuscript. WARNING! Will overwrite actual values and needs to be turned off/adjusted manually when data is updated
customLabels <- TRUE
##CHANGE ME to change the width (in cm) of the output. ##CHANGE ME to change the width (in cm) of the output.
wp <- 20 wp <- 12
#CHANGE ME to change the height (in cm) of the output. #CHANGE ME to change the height (in cm) of the output.
hp <- 15 hp <- 10
#CHANGE ME to change the resolution (in dpi) of the output. # ASVs (before LULU!) that exactly match isolates seqs (identified by vsearch),
res <- 300 # read in from fasta
#Get total number of reads and number of Pythium/Phytophthora reads (without contro samples!) to find out percentage of Pyhium and Phytophthora exASVs_seqs <- readDNAStringSet(paste0(path,"/exact_ASVs_v3.fasta"))
total_reads_no_contr <- ps_tbl %>% iso_seqs <- readDNAStringSet(paste0(path,"/trim_isos_v3.fasta"))
dplyr::filter(Pass == "Yes") %>% iso_ASV_match <- read_delim(paste0(path,"/match_trim-isolates_ASVs_v3.txt"),
dplyr::mutate(sumAbu = sum(Abundance)) %>% delim = "\t",
select(sumAbu) %>% col_names = c("Iso_Soil" ,"OTU",
unique() %>% "Match", "Mism", "Gap"))
unlist() %>% #Parse headers to get ASV IDs listed in vector
unname() exASVs <- names(exASVs_seqs) %>% str_remove("\\|\\|.*")
pythium_reads_no_contr <- ps_tbl %>% isos <- names(iso_seqs)
dplyr::filter(Pass == "Yes") %>%
dplyr::group_by(Genus) %>%
summarise(sumAbu = sum(Abundance)) %>%
dplyr::filter(Genus == "Pythium") %>%
select(sumAbu) %>%
unlist() %>%
unname()
phytoph_reads_no_contr <- ps_tbl %>%
dplyr::filter(Pass == "Yes") %>%
dplyr::group_by(Genus) %>%
summarise(sumAbu = sum(Abundance)) %>%
dplyr::filter(Genus == "Phytophthora") %>%
select(sumAbu) %>%
unlist() %>%
unname()
prctg_pyth_phyt_str <- paste0("\n\nTotal number of reads in the samples: ",
total_reads_no_contr ,
", number of Pythium reads: ",
pythium_reads_no_contr,
", number of Phytophthora reads: ",
phytoph_reads_no_contr,
"\nPercentage of Pythium total: ",
pythium_reads_no_contr/total_reads_no_contr*100,
"%\nPercentage of Phytophthora total: ",
phytoph_reads_no_contr/total_reads_no_contr*100, "%")
#Most abundant at specified levels, phyloseq objects
ps.topnt <- cuphyr::abundant_tax_physeq(
physeq = ps.trans, lvl = taxlvl, top = numt,
ignore_na = TRUE, silent = FALSE)
ps.topnt2 <- cuphyr::abundant_tax_physeq(
physeq = ps.trans,lvl = taxlvl2, top = numt2,
ignore_na = TRUE, silent = FALSE)
#Abundance lists (sorted by abundance)
taxlist_topnt1 <- cuphyr::abundant_tax_physeq(
physeq = ps.trans, lvl = taxlvl, top = numt,
output_format = "tops", ignore_na = TRUE)
taxlist_topnt2 <- cuphyr::abundant_tax_physeq(
physeq = ps.trans, lvl = taxlvl2, top = numt2,
output_format = "tops", ignore_na = TRUE)
#tidy evaluation translation magic (needed so the dplyr functions that extract 'tax_lookup' below can parse the strings from the variable)
lvl1 <- sym(taxlvl)
lvl2 <- sym(taxlvl2)
lvl1 <- enquo(lvl1)
lvl2 <- enquo(lvl2)
#make lookup table to guide coloring
tax_lookup <- ps_tbl %>%
group_by(!!lvl1, !!lvl2) %>%
summarise() %>%
filter(!!lvl2 %in% taxlist_topnt2)
#Define fixed color schemes for classes for consistency between plots. This is a color scheme specified for this particular published figure. Can be used for orientation in general cases.
#Background colors to be overwritten for special cases
colors_lvl1 <- viridis(length(taxlist_topnt1), option = "D", begin = 1 , end = 0.5)
names(colors_lvl1) <- taxlist_topnt1
colors_lvl2 <- viridis(length(taxlist_topnt2), option = "D", begin = 1 , end = 0.5)
names(colors_lvl2) <- taxlist_topnt2
colors_combo <- c(colors_lvl1, colors_lvl2)
#Special cases
colors_combo["Chrysophyceae"] <- viridis_yellows[2]
colors_combo["Pedospumella"] <- viridis_yellows[1]
colors_combo["'Spumella-like'"] <- viridis_light[1]
colors_combo["Oomycetes"] <- viridis_reds[4]
colors_combo["Phytophthora"] <- viridis_reds[1]
colors_combo["Pythium"] <- viridis_reds[4]
colors_combo["Aphanomyces"] <- viridis_blues[1]
colors_combo["Globisporangium"] <- "#4C0000FF"
colors_combo["Pythiogeton"] <- viridis_blues[2]
colors_combo["Saprolegnia"] <- viridis_blues[4]
#Format legend labels
if(customLabels){
labelling_lvl2 <- c(expression(paste(italic("Spumella"), "-like")),
"uncultured fungus", expression(italic("Aphanomyces"),
italic("Globisporangium"), italic("Haptoglossa"), italic("Pedospumella"),
italic("Phytophthora"), italic("Pythiogeton"), italic("Pythiopsis"),
italic("Pythium")))
}else{
labelling_lvl2 <- names(colors_lvl2) %>% sort()
}
#Re-label the facets #Re-label the facets
enrich_labs <- c("Before enrichment", "After enrichment") enrich_labs <- c("Before enrichment", "After enrichment")
names(enrich_labs) <- c("Before enrichment", "Enriched") names(enrich_labs) <- c("Before enrichment", "Enriched")
theme_bar <- theme(legend.position = "bottom", samp_sums <- ps_tbl %>%
legend.key.size = unit(0.4, "cm"), dplyr::group_by(SampleIDs) %>%
legend.spacing.x = unit(0.3, 'cm'), summarise(sum_Abundance = sum(Abundance))
text = element_text(size = 10), ps_tr_tbl <- left_join(ps_tbl, samp_sums, by = "SampleIDs") %>%
strip.text.x = element_text(size = 10), dplyr::mutate(Rel_Abundance = Abundance/sum_Abundance) %>%
strip.background = element_blank(), dplyr::mutate(Prcnt_Abundance = Rel_Abundance*100)
axis.text.x = element_text(angle = 0, hjust = 0.5)) exASV_counts <- ps_tr_tbl %>%
topnt <- plot_bar(ps.topnt, x = primaryPar, fill = taxlvl, dplyr::select(OTU:Bait, Sanger, Rel_Abundance, Prcnt_Abundance) %>%
title = paste("Top", numt, "Classes")) + dplyr::filter(OTU %in% exASVs) %>%
facet_wrap(paste0("~", secondaryPar), dplyr::filter(Abundance >= 1)
labeller = labeller(Bait = enrich_labs)) + iso_ASVs <- ggplot(exASV_counts, aes(x = OTU,
scale_x_discrete(breaks = c("S05" ,"S15", "S25", "S35", "S45" ,"S55")) + y = reorder(Soil, dplyr::desc(Soil)),
scale_fill_manual(values = colors_combo) + fill = Prcnt_Abundance)) +
ylab("Relative abundance") + geom_tile() +
theme_bar + scale_fill_viridis() +
xlab("Soil samples S01-S64") facet_wrap(~Bait, labeller = labeller(Bait = enrich_labs)) +
topnt_summary <- cuphyr::summarise_physeq(physeq = ps, theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
ASV_sublist = taxa_names(ps.topnt), panel.grid = element_blank(),
sublist_id = axis.title = element_blank(),
paste0("top ", numt, " ", taxlvl), strip.background = element_rect(fill = "white", color = NA)) +
samp_names = FALSE) labs(fill = "Abund. (%)") +
topnt2 <- plot_bar(ps.topnt2, x = primaryPar, fill = taxlvl2, scale_x_discrete(limits = exASVs)
title = paste("Top", numt2, "Genera")) + save_plot(iso_ASVs, plot_name = "IsolateSeqs_ASVs_heatmap_Fig4")
facet_wrap(paste0("~", secondaryPar), iso_ASVs
labeller = labeller(Bait = enrich_labs)) +
scale_x_discrete(breaks = c("S05" ,"S15", "S25", "S35", "S45" ,"S55")) +
scale_fill_manual(values = colors_combo, labels = labelling_lvl2) +
ylab("Relative abundance") +
theme_bar + xlab("Soil samples S01-S64")
topnt2_summary <- cuphyr::summarise_physeq(physeq = ps,
ASV_sublist = taxa_names(ps.topnt2),
sublist_id =
paste0("top ", numt2, " ", taxlvl2),
samp_names = FALSE)
#Combine the plots
combobar1 <- ggarrange(topnt, topnt2, nrow = 2, labels = c("A", "B"), align = "hv")
#Save
save_plot(combobar1, plot_name = "Topn_barplot_Fig2")
#Print to standard out
cat(topnt_summary, "\n\n", topnt2_summary)
cat(prctg_pyth_phyt_str)
combobar1
cat("Chunk successfully run") cat("Chunk successfully run")
#CHANGE ME to the amount of top taxlvl to plot (e.g. 'numt = 20' plots Top20 Taxa at taxlvl)
numt <- 5
#CHANGE ME to the amount of top taxa you want to plot (e.g. 'numt2 = 25' plots Top25 Taxa at taxlvl2)
numt2 <- 10
#CHANGE ME to the primary parameter of interest (x-axis). Accepted values are the column headers in the descriptor file.
primaryPar <- "Soil"
##CHANGE ME to the secondary parameter of interest (panels). Accepted values are the column headers in the descriptor file.
secondaryPar <- "Bait"
##CHANGE ME to the taxonomic level of interest (color coding). Accepted values are available taxonomic levels. 'taxlvl' should be higher than 'taxlvl2', e.g. Class and Genus
taxlvl <- "Class"
taxlvl2 <- "Genus"
#Adds formatted legend labels for the specific output in the manuscript. WARNING! Will overwrite actual values and needs to be turned off/adjusted manually when data is updated
customLabels <- TRUE
##CHANGE ME to change the width (in cm) of the output. ##CHANGE ME to change the width (in cm) of the output.
wp <- 20 wp <- 17.5
#CHANGE ME to change the height (in cm) of the output. #CHANGE ME to change the height (in cm) of the output.
hp <- 15 hp <- 10
#CHANGE ME to change the resolution (in dpi) of the output. #CHANGE ME to change the resolution (in dpi) of the output.
res <- 300 res <- 150
#Get total number of reads and number of Pythium/Phytophthora reads (without contro samples!) to find out percentage of Pyhium and Phytophthora #calculate averages for display in graphs
total_reads_no_contr <- ps_tbl %>% avgNames <- c("PPUntr", "PPBait","ChrUntr", "ChrBait")
dplyr::filter(Pass == "Yes") %>% avgs <- c(mean(ranksUntrPP_tot$Abundance),mean(ranksBaitPP_tot$Abundance),
dplyr::mutate(sumAbu = sum(Abundance)) %>% mean(ranksUntrAlgae_tot$Abundance),mean(ranksBaitAlgae_tot$Abundance))
select(sumAbu) %>% names(avgs) <- avgNames
unique() %>% avgs <- format(round(avgs, 1), nsmall = 1)
unlist() %>% #Define fixed color schemes for consistency between plots
unname() plotvars <- c("Chrysophyceae_tot", "All oomycetes",
pythium_reads_no_contr <- ps_tbl %>% "Other oomycetes", "Other classes", "PP_tot")
dplyr::filter(Pass == "Yes") %>% getPalette <- colorRampPalette(viridis(length(plotvars)))
dplyr::group_by(Genus) %>% plotPalette <- getPalette(length(plotvars))
summarise(sumAbu = sum(Abundance)) %>% names(plotPalette) <- plotvars
dplyr::filter(Genus == "Pythium") %>% plotPalette["Chrysophyceae_tot"] <- viridis_yellows[2]
select(sumAbu) %>% plotPalette["PP_tot"] <- viridis_reds[4]
unlist() %>% customTheme <- theme(legend.position = "none",
unname() axis.title = element_blank(),
phytoph_reads_no_contr <- ps_tbl %>% title = element_text(size = 7))
dplyr::filter(Pass == "Yes") %>% customThemeA <- theme(legend.position = "none",
dplyr::group_by(Genus) %>% legend.title = element_blank(),
summarise(sumAbu = sum(Abundance)) %>% legend.key.size = unit(5,"mm"),
dplyr::filter(Genus == "Phytophthora") %>% axis.title = element_blank(),
select(sumAbu) %>% title = element_text(size = 10))
unlist() %>% xUpperLim <- nrow(ranksBaitAlgae_tot) + 1
unname() ytextl = 0.05
prctg_pyth_phyt_str <- paste0("\n\nTotal number of reads in the samples: ", xtextl = 15
total_reads_no_contr , ytexts = 90000
", number of Pythium reads: ", xtexts = 45
pythium_reads_no_contr, pAllBait <- ggplot(data = ranksBaitPP_tot,
", number of Phytophthora reads: ", aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
phytoph_reads_no_contr, scale_fill_manual(values = plotPalette) +
"\nPercentage of Pythium total: ", ggtitle("PP - after enrichment") +
pythium_reads_no_contr/total_reads_no_contr*100, xlim(0,xUpperLim) +
"%\nPercentage of Phytophthora total: ", customTheme +
phytoph_reads_no_contr/total_reads_no_contr*100, "%") annotate(geom = 'text', label = paste("Mean", avgs["PPBait"]),
#Most abundant at specified levels, phyloseq objects x = xtexts, y = ytexts) +
ps.topnt <- cuphyr::abundant_tax_physeq( scale_y_continuous(limits = c(0,105000))
physeq = ps.trans, lvl = taxlvl, top = numt, pAllUntr <- ggplot(data = ranksUntrPP_tot,
ignore_na = TRUE, silent = FALSE) aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
ps.topnt2 <- cuphyr::abundant_tax_physeq( scale_fill_manual(values = plotPalette) +
physeq = ps.trans,lvl = taxlvl2, top = numt2, ggtitle("PP - before enrichment") +
ignore_na = TRUE, silent = FALSE) xlim(0,xUpperLim) +
#Abundance lists (sorted by abundance) customTheme +
taxlist_topnt1 <- cuphyr::abundant_tax_physeq( annotate(geom = 'text', label = paste("Mean", avgs["PPUntr"]),
physeq = ps.trans, lvl = taxlvl, top = numt, x = xtexts, y = ytexts) +
output_format = "tops", ignore_na = TRUE) scale_y_continuous(limits = c(0,105000))
taxlist_topnt2 <- cuphyr::abundant_tax_physeq( AlgBait <- ggplot(data = ranksBaitAlgae_tot,
physeq = ps.trans, lvl = taxlvl2, top = numt2, aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
output_format = "tops", ignore_na = TRUE) scale_fill_manual(values = plotPalette) +
#tidy evaluation translation magic (needed so the dplyr functions that extract 'tax_lookup' below can parse the strings from the variable) ggtitle("Chrysophyceae - after enrichment") +
lvl1 <- sym(taxlvl) xlim(0,xUpperLim) +
lvl2 <- sym(taxlvl2) customTheme +
lvl1 <- enquo(lvl1) annotate(geom = 'text', label = paste("Mean", avgs["ChrBait"]),
lvl2 <- enquo(lvl2) x = xtexts, y = ytexts) +
#make lookup table to guide coloring scale_y_continuous(limits = c(0,105000))
tax_lookup <- ps_tbl %>% AlgUntr <- ggplot(data = ranksUntrAlgae_tot,
group_by(!!lvl1, !!lvl2) %>% aes(x = Rank, y = Abundance, fill = Set)) +
summarise() %>% geom_col() +
filter(!!lvl2 %in% taxlist_topnt2) scale_fill_manual(values = plotPalette) +
#Define fixed color schemes for classes for consistency between plots. This is a color scheme specified for this particular published figure. Can be used for orientation in general cases. ggtitle("Chrysophyceae - before enrichment") +
#Background colors to be overwritten for special cases xlim(0,xUpperLim) +
colors_lvl1 <- viridis(length(taxlist_topnt1), option = "D", begin = 1 , end = 0.5) customTheme +
names(colors_lvl1) <- taxlist_topnt1 annotate(geom = 'text', label = paste("Mean", avgs["ChrUntr"]),
colors_lvl2 <- viridis(length(taxlist_topnt2), option = "D", begin = 1 , end = 0.5) x = xtexts, y = ytexts) +
names(colors_lvl2) <- taxlist_topnt2 scale_y_continuous(limits = c(0,105000))
colors_combo <- c(colors_lvl1, colors_lvl2) A.combo <- ggarrange(pAllUntr, pAllBait, nrow = 2, labels = c("A"), vjust = 1)
#Special cases B.combo <- ggarrange(AlgUntr, AlgBait, nrow = 2, labels = c("B"), vjust = 1)
colors_combo["Chrysophyceae"] <- viridis_yellows[2] A.combo.anot <- annotate_figure(A.combo, left = text_grob("Abundance", rot = 90))
colors_combo["Pedospumella"] <- viridis_yellows[1] combo.AB <- ggarrange(A.combo.anot, B.combo, ncol = 2) +
colors_combo["'Spumella-like'"] <- viridis_light[1] theme(plot.margin = margin(10,10,10,10)) +
colors_combo["Oomycetes"] <- viridis_reds[4] scale_fill_manual(values = plotPalette)
colors_combo["Phytophthora"] <- viridis_reds[1] #Save plot
colors_combo["Pythium"] <- viridis_reds[4] save_plot(combo.AB, plot_name = "Class_absolute_ranking_Supfig3")
colors_combo["Aphanomyces"] <- viridis_blues[1] combo.AB
colors_combo["Globisporangium"] <- "#4C0000FF" cat("Chunk successfully run")
colors_combo["Pythiogeton"] <- viridis_blues[2] sessionInfo()
colors_combo["Saprolegnia"] <- viridis_blues[4] ## CHANGE ME to change the width (in cm) of the output.
#Format legend labels wp <- 17.5
if(customLabels){ # CHANGE ME to change the height (in cm) of the output.
labelling_lvl2 <- c(expression(paste(italic("Spumella"), "-like")), hp <- 15
"uncultured fungus", expression(italic("Aphanomyces"), # CHANGE ME to change the resolution (in dpi) of the output.
italic("Globisporangium"), italic("Haptoglossa"), italic("Pedospumella"), res <- 150
italic("Phytophthora"), italic("Pythiogeton"), italic("Pythiopsis"), # calculate averages for display in graphs
italic("Pythium"))) avgNames <- c("PPUntr", "PPBait","ChrUntr", "ChrBait",
}else{ "O-PUntr", "O-PBait", "OthUntr", "OthBait")
labelling_lvl2 <- names(colors_lvl2) %>% sort() avgs <- c(mean(ranksUntrPP$Abundance),
} mean(ranksBaitPP$Abundance),
#Re-label the facets mean(ranksUntrAlgae$Abundance),
enrich_labs <- c("Before enrichment", "After enrichment") mean(ranksBaitAlgae$Abundance),
names(enrich_labs) <- c("Before enrichment", "Enriched") mean(ranksUntrOwoP$Abundance),
theme_bar <- theme(legend.position = "bottom", mean(ranksBaitOwoP$Abundance),
legend.key.size = unit(0.4, "cm"), mean(rankedOthersUntr$OtherASVcounts),
legend.spacing.x = unit(0.3, 'cm'), mean(rankedOthersBait$OtherASVcounts))
text = element_text(size = 10), names(avgs) <- avgNames
strip.text.x = element_text(size = 10), # multiply to transform relative abundances into avg percentages and round to
strip.background = element_blank(), # one decimal
axis.text.x = element_text(angle = 0, hjust = 0.5), avgs <- avgs*100
legend.text.align = 0) avgs <- format(round(avgs, 1), nsmall = 1)
topnt <- plot_bar(ps.topnt, x = primaryPar, fill = taxlvl, # Define fixed color schemes for consistency between plots
title = paste("Top", numt, "Classes")) + plotvars <- c("Chrysophyceae", "All oomycetes",
facet_wrap(paste0("~", secondaryPar), "Other oomycetes", "Other classes", "PP")
labeller = labeller(Bait = enrich_labs)) + # Just to catch others/exceptions, make a general palette that will be
scale_x_discrete(breaks = c("S05" ,"S15", "S25", "S35", "S45" ,"S55")) + # overwritten below
scale_fill_manual(values = colors_combo) + plotPalette <- viridis(length(plotvars))
ylab("Relative abundance") + names(plotPalette) <- plotvars
theme_bar + plotPalette["Other classes"] <- "#C0C0C0"
xlab("Soil samples S01-S64") plotPalette["Chrysophyceae"] <- viridis_yellows[2]
topnt_summary <- cuphyr::summarise_physeq(physeq = ps, plotPalette["PP"] <- viridis_reds[4]
ASV_sublist = taxa_names(ps.topnt), plotPalette["Other oomycetes"] <- viridis_reds[1]
sublist_id = plotPalette["All oomycetes"] <- viridis_reds[1]
paste0("top ", numt, " ", taxlvl), # plot total ASV counts, sorted by rank and color by class.
samp_names = FALSE) AllBait <- ggplot(data = OthersRankPP,
topnt2 <- plot_bar(ps.topnt2, x = primaryPar, fill = taxlvl2, aes(x = Rank, y = Abundance, fill = Set))
title = paste("Top", numt2, "Genera")) + AllUntr <- ggplot(data = UntrOthersRankPP,
facet_wrap(paste0("~", secondaryPar), aes(x = Rank, y = Abundance, fill = Set))
labeller = labeller(Bait = enrich_labs)) + customTheme <- theme(legend.position = "none",
scale_x_discrete(breaks = c("S05" ,"S15", "S25", "S35", "S45" ,"S55")) + axis.title = element_blank(),
scale_fill_manual(values = colors_combo, labels = labelling_lvl2) + title = element_text(size = 6))
ylab("Relative abundance") + customThemeA <- theme(legend.position = "none",
theme_bar + xlab("Soil samples S01-S64") legend.title = element_blank(),
topnt2_summary <- cuphyr::summarise_physeq(physeq = ps, legend.key.size = unit(5,"mm"),
ASV_sublist = taxa_names(ps.topnt2), axis.title.x = element_blank())
sublist_id = xUpperLim <- nrow(OthersRankPP) + 1
paste0("top ", numt2, " ", taxlvl2), ytextl = 0.05
samp_names = FALSE) xtextl = 15
#Combine the plots ytexts = 0.75
combobar1 <- ggarrange(topnt, topnt2, nrow = 2, labels = c("A", "B"), align = "hv") xtexts = 45
#Save middleline <- geom_blank() #geom_hline(yintercept = 0.5,
save_plot(combobar1, plot_name = "Topn_barplot_Fig2") # linetype="dashed", color="black")
#Print to standard out pAllBait <- AllBait + geom_col() +
cat(topnt_summary, "\n\n", topnt2_summary) geom_col(data = AlgaeInvRankOom) +
cat(prctg_pyth_phyt_str) geom_col(data = OomRankPP) +
combobar1 geom_col(data = ranksBaitPP) +
scale_fill_manual(values = plotPalette) +
ggtitle("After enrichment") + xlab("Sample rank") +
labs(fill = "") + ylim(0,1) + xlim(0,xUpperLim) +
customThemeA + middleline +
annotate(geom = 'text', label = paste("Mean PP", avgs["PPBait"], "%"),
x = xtextl, y = ytextl) +
ylab("Rel. abundance")
pAllUntr <- AllUntr + geom_col() +
geom_col(data = UntrAlgaeInvRankOom) +
geom_col(data = UntrOomRankPP) +
geom_col(data = ranksUntrPP) +
scale_fill_manual(values = plotPalette) +
ggtitle("Before enrichment") + xlab("Sample rank") +
labs(fill = "") + ylim(0,1) + xlim(0,xUpperLim) +
customThemeA + middleline +
annotate(geom = 'text', label = paste("Mean PP", avgs["PPUntr"], "%"),
x = xtextl, y = ytextl) +
ylab("Rel. abundance")
AlgBait <- ggplot(data = ranksBaitAlgae,
aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Chrysophyceae - after enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["ChrBait"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
OomBait <- ggplot(data = ranksBaitOwoP,
aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Other oomycetes - after enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["O-PBait"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
OthBait <- ggplot(data = rankedOthersBait,
aes(x = Rank, y = OtherASVcounts, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Other classes - after enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["OthBait"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
AlgUntr <- ggplot(data = ranksUntrAlgae,
aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Chrysophyceae - before enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["ChrUntr"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
OomUntr <- ggplot(data = ranksUntrOwoP,
aes(x = Rank, y = Abundance, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Other oomycetes - before enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["O-PUntr"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
OthUntr <- ggplot(data = rankedOthersUntr,
aes(x = Rank, y = OtherASVcounts, fill = Set)) + geom_col() +
scale_fill_manual(values = plotPalette) +
ggtitle("Other classes - before enrichment") +
xlim(0,xUpperLim) + customTheme + middleline +
annotate(geom = 'text', label = paste("Mean", avgs["OthUntr"], "%"),
x = xtexts, y = ytexts) +
scale_y_continuous(breaks = c(0,0.5,1), limits = c(0,1))
#A.combo <- ggarrange(pAllUntr, pAllBait, nrow=2,
# labels = c("A"), vjust = 1, common.legend = TRUE, legend = "bottom")
untr_smallplots <- ggarrange(AlgUntr, OomUntr, OthUntr, nrow = 3)
untr_smallplots_anot <- annotate_figure(untr_smallplots,
left = text_grob("Rel. abundance",
rot = 90, size = 11))
bait_smallplots <- ggarrange(AlgBait, OomBait, OthBait, nrow = 3)
bait_smallplots_anot <- annotate_figure(bait_smallplots,
left = text_grob("Rel. abundance",
rot = 90, size = 11))
combo <- ggarrange(pAllUntr, untr_smallplots_anot, pAllBait,
bait_smallplots_anot, nrow = 2, ncol = 2, widths = c(2,1),
labels = c("A", "", "B", ""), common.legend = TRUE,
legend = "bottom") + theme(plot.margin = margin(10,10,10,10))
#Save plot
save_plot(combo, plot_name = "Class_ranking_Fig3")
combo
cat("Chunk successfully run") cat("Chunk successfully run")
# CHANGE ME according to the location containing 'seqtab_nochim.rds' and other relevant files.
path = "Data_files/"
# CHANGE ME to 'TRUE' to list all samples and generate an empty metadata file
optional_sample_check <- "FALSE"
# CHANGE ME to TRUE to update the cuphyr package. This is a package the author published to separate out general purpose functions from this script. If the package is not installed, it will be installed, regardless of the value chosen here.
update_cuphyr <- TRUE
################# NO CHANGES NECESSARY BELOW #################
#global knitr settings to determine the output of the markdown document upon 'knitting'
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(root.dir = paste0(path))
#Messages and warnings are turned off globally for a cleaner output. Change to TRUE for complete output
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(tidy = TRUE)
#CuPhyR is the author's own package that contains several functions used to process phyloseq objects. It is hosted on GitHub and installed/updated using this command if update_cuphyr is TRUE
if (update_cuphyr | !"cuphyr" %in% rownames(installed.packages())) {
devtools::install_github("simeross/cuphyr")
}
#Loading required packages
library(phyloseq)
library(dada2)
.Rhistory
This diff is collapsed.
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment