/
ordination.R
154 lines (140 loc) · 7.17 KB
/
ordination.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
#' @title Ordination Plot
#' @description This function allows for ordination (which helps us to distinguish beta diversity relationships) to be plotted as well as for the corresponding data to be returned.
#' @param obj An object to be converted to a taxmap object with \code{\link[MicrobiomeR]{create_taxmap}}.
#' @param method Choose an ordination method from 'PCoA', 'CCA', 'NMDS' or 'DPCoA', Default: 'PCoA'
#' @param distance Choose a \code{\link[phyloseq]{distance}} method or use a pre-computed \code{\link{dist}}-class object, Default: 'wunifrac'
#' @param color Choose the group or factor of which colors will be mapped to, Default: 'TreatmentGroup'
#' @param title The title of the plot, Default: NULL
#' @param only_data Allows for only ordination data to be generated, Default: FALSE
#' @return By default, it returns an ordination plot.
#' @examples
#' \dontrun{
#' if (interactive()) {
#' # An example ordination plot
#' library(MicrobiomeR)
#' data <- analyzed_silva
#' plot <- ordination_plot(obj = data)
#' plot
#' }
#' }
#' @export
#' @family Visualizations
#' @rdname ordination_plot
#' @seealso View \code{\link{save_ordination_plots}} to save your ordination plot or multiple ordination plots.
#' @importFrom metacoder as_phyloseq
#' @importFrom phyloseq ordinate plot_ordination
#' @importFrom ggplot2 element_text geom_point theme element_blank guide_legend guides ggtitle unit scale_fill_manual labs scale_color_manual
ordination_plot <- function(obj, method = "PCoA", distance = "wunifrac", color = "TreatmentGroup", title = NULL, only_data = FALSE) {
metacoder_object <- validate_MicrobiomeR_format(
obj = create_taxmap(obj),
valid_formats = c("analyzed_format")
)
# Convert taxmap object to a phyloseq object.
phyloseq_object <- metacoder::as_phyloseq(metacoder_object, otu_table = "otu_abundance", phy_tree = "phy_tree")
# Find out how many groups are in the "color" metadata
meta <- phyloseq::sample_data(phyloseq_object)
num_groups <- length(unique(meta[[color]]))
group_names <- unique(meta[[color]])
# Create palette and palette lists
ord_palette <- c("#3288bd", "#d53e4f", "#62954C", "#C59144")
names(ord_palette) <- group_names
ord <- phyloseq::ordinate(phyloseq_object, method, distance)
plot <- phyloseq::plot_ordination(physeq = phyloseq_object, ordination = ord, color = color) +
ggplot2::geom_point(size = 2) + ggplot2::theme(
text = ggplot2::element_text(size = 11, family = "Arial", face = "bold"),
axis.text.y = ggplot2::element_text(margin = ggplot2::unit(c(0.5, 0.5, 0.5, 0.5), "cm")),
axis.text.x = ggplot2::element_text(margin = ggplot2::unit(c(0.5, 0.5, 0.5, 0.5), "cm")), panel.grid.major = ggplot2::element_blank(), panel.grid.minor = element_blank(),
axis.ticks.length = ggplot2::unit(-0.25, "cm"),
strip.background = ggplot2::element_rect(fill = "white"),
strip.text = ggplot2::element_text(colour = "black"),
panel.background = element_blank(),
axis.line = ggplot2::element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill = NA, size = 1)
) + ggplot2::stat_ellipse(
geom = "polygon", alpha = .2,
ggplot2::aes(fill = !!sym(color)),
show.legend = FALSE
)
plot <- plot + ggplot2::scale_fill_manual(values = ord_palette)
plot <- plot + ggplot2::labs(color = color) + ggplot2::scale_color_manual(values = ord_palette)
if (!is.null(title)) {
plot <- plot + ggplot2::ggtitle(title) + ggplot2::annotate("segment", x = Inf, xend = -Inf, y = Inf, yend = Inf, color = "black", lwd = 1) +
ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(fill = NA)))
}
df <- phyloseq::plot_ordination(phyloseq_object, ord, justDF = TRUE, axes = 1:4)
if (only_data == TRUE) {
return(df)
} else {
return(plot)
}
}
#' @title Ordination Plots
#' @description Generate plots for a list of ordination methods and distances.
#' @param obj An object to be converted to a metacoder object with \code{\link[MicrobiomeR]{create_taxmap}}.
#' @param methods A list of ordination methods, Default: 'c("PCoA", "NMDS")'
#' @param distances A list of distance methods, Default: 'c("wunifrac", "unifrac", "bray")'
#' @param color Choose the group or factor of which colors will be mapped to, Default: 'TreatmentGroup'
#' @param select_otu_table The data table to use in the observation data. Default: "otu_proportions"
#' @return Returns a melted dataframe.
#' @export
#' @family Visualizations
#' @rdname ordination_plots
#' @rdname ordination_plot
ordination_plots <- function(obj, methods = c("PCoA", "NMDS"), distances = c("wunifrac", "unifrac", "bray"),
color = "TreatmentGroup", select_otu_table = "otu_proportions") {
if (is.null(methods)) {
methods <- c("PCoA", "NMDS")
} else if (length(methods) < 2) {
stop("Use the ordination_plot function for generating a plot for 1 method or 1 distance.")
}
if (is.null(distances)) {
distances <- c("wunifrac", "unifrac", "bray")
} else if (length(methods) < 2 && length(distances) < 2) {
stop("Use the ordination_plot function for generating a plot for 1 method and 1 distance.")
}
ordination_plots <- list()
for (m in methods) {
for (d in distances) {
ordination_plots[[paste0(m, "_", d)]] <- ordination_plot(obj, method = m, distance = d, color = color, title = NULL, only_data = FALSE)
}
}
return(ordination_plots)
}
#' @title Save Ordination Plots
#' @description This function saves ordination plots stored in a listlike object to an output folder.
#' @param ord_plots An ordination plot list generated by ordination_plots.
#' @param format The format of the output image. Default: 'tiff'
#' @param start_path The starting path of the output directory. Default: 'output'
#' @param ... An optional list of parameters to use in the output_dir function.
#' @return An output directory that contains ordination plots.
#' @pretty_print TRUE
#' @details This function creates an appropriate output directory, where it saves publication ready
#' plots.
#' @examples
#' \dontrun{
#' if (interactive()) {
#' }
#' }
#' @export
#' @family Visualizations
#' @rdname save_ordination_plots
#' @seealso
#' \code{\link[MicrobiomeR]{output_dir}}
#'
#' \code{\link[ggplot2]{ggsave}}
#' @importFrom ggplot2 ggsave
#' @importFrom crayon green
#' @importFrom glue glue
save_ordination_plots <- function(ord_plots, format = "tiff", start_path = "output", ...) {
# Create the relative path to the ordination plots. By default the path will be <pwd>/output/<experiment>/ordination/<format(Sys.time(), "%Y-%m-%d_%s")>
# With the parameters set the full path will be <pwd>/output/<experiment>/ordination/<extra_path>.
full_path <- output_dir(start_path = start_path, plot_type = "ordination", ...)
message(glue::glue(crayon::yellow("Saving Ordination plots to the following directory: \n", "\r\t{full_path}")))
# Iterate the plot list and save them in the proper directory
for (method in names(ord_plots)) {
if (method != "metacoder_object") {
message(crayon::green("Saving the {method} ordination plot."))
ggplot2::ggsave(filename = sprintf("%s_ordination.%s", tolower(method), format), plot = ord_plots[[method]], device = format, path = full_path, dpi = 500)
}
}
}