R/metabolomics_plots.R
In NicWir/VisomX: User-Friendly and Comprehensive Analysis and Visualization of Omics Data
Defines functions
met.plot_PLSImpScatter
Documented in
met.plot_PLSImpScatter
#' Plot a heatmap of compounds with missing values
#'
#' \code{met.plot_missval} generates a heatmap of compounds
#' with missing values to discover whether values are missing at random or not.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container with missing values.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF and PNG file?
#' @param fontsize (Numeric) Font sized used in the plot.
#' @return The input mSet object with added heatmap indicating whether values are missing (0) or not (1). The plot can be retrieved from within R via \code{ComplexHeatmap::draw(mSetObj$imgSet$missval_heatmap.plot)}
#' (generated by \code{\link[ComplexHeatmap]{Heatmap}}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_missval <- function (mSetObj, plot = TRUE, imgName = "Plots/MissValPlot", format = "pdf", export = TRUE, fontsize = 12)
{
assay <- mSetObj[["dataSet"]][["data_orig"]]
if (!any(is.na(assay)) & !any(assay==0)) {
stop("No missing or zero values in '", deparse(substitute(se)),
"'", call. = FALSE)
}
imgName = paste(imgName,".", format,
sep = "")
# Make assay data binary (1 = valid value, 0 = missing value)
df <- assay %>% data.frame(check.names = FALSE)
df <- t(sapply(df, as.numeric))
colnames(df) <- rownames(assay)
missval <-
df[apply(df, 1, function(x)
any(is.na(x))) | apply(df, 1, function(x)
any(x==0)),]
missval <- ifelse(is.na(missval)|missval==0, 0, 1)
# Plot binary heatmap
ht2 = ComplexHeatmap::Heatmap(missval, col = c("white", "black"),
column_names_side = "top", show_row_names = FALSE,
show_column_names = TRUE, name = "Missing values pattern",
column_names_gp = grid::gpar(fontsize = fontsize), heatmap_legend_param = list(at = c(0,
1), labels = c("Missing value", "Valid value")))
if (export == TRUE){
message(paste0("Exporting Missing values pattern to:\n\"", getwd(), "/", imgName, "\""))
if (format == "pdf") {
grDevices::pdf(imgName)
ComplexHeatmap::draw(ht2, heatmap_legend_side = "top")
grDevices::dev.off()
}
if (format == "png") {
grDevices::png(imgName,
width = 6, height = 6, units = 'in', res = 300)
ComplexHeatmap::draw(ht2, heatmap_legend_side = "top")
grDevices::dev.off()
}
mSetObj$imgSet$missval_heatmap <- imgName
}
if (plot == TRUE){
ComplexHeatmap::draw(ht2, heatmap_legend_side = "top")
}
mSetObj$imgSet$missval_heatmap.plot <- ComplexHeatmap::draw(ht2, heatmap_legend_side = "top")
return(mSetObj)
}
#' Visualize intensities of compounds with missing values
#'
#' \code{met.plot_detect} generates density and CumSum plots
#' of compound intensities with and without missing values.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container with missing values.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF and PNG file?
#' @param basesize (Numeric) Font sized used in the plot.
#' @return Density and CumSum plots of intensities of proteins with and without missing values. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$missval_density.plot)}.
#' (generated by \code{\link[ggplot2]{ggplot}}).
#' @import ggplot2
#' @references adapted from package DEP (\url{https://bioconductor.org/packages/DEP/}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_detect <- function (mSetObj, plot = TRUE, imgName = "Plots/DensMissPlot", format = "pdf", export = TRUE, basesize = 10)
{
assay <- mSetObj[["dataSet"]][["data_orig"]] %>% replace(.==0, NA)
if (!any(is.na(assay)) & !any(assay==0)) {
stop("No missing or zero values in '", if_else(!is.null(mSetObj$dataSet$prenorm),
deparse(substitute(mSetObj[["dataSet"]][["prenorm"]])),
if_else(!is.null(mSetObj$dataSet$data_proc),
deparse(substitute(mSetObj[["dataSet"]][["data_proc"]])),deparse(substitute(mSetObj[["dataSet"]][["data_orig"]])))),
"'", call. = FALSE)
}
imgName = paste(imgName,".", format,
sep = "")
df <- assay %>% data.frame(check.names = FALSE)+1
df <- t(sapply(df, as.numeric)) %>% data.frame() %>% log10()
colnames(df) <- rownames(assay)
df <- df %>% tibble::rownames_to_column() %>% tidyr::gather(ID, val, -rowname)
stat <- df %>% group_by(rowname) %>% summarize(mean = mean(val,
na.rm = TRUE), missval = any(is.na(val)|val==0))
cumsum <- stat %>% group_by(missval) %>% arrange(mean) %>%
mutate(num = 1, cs = cumsum(num), cs_frac = cs/n())
p1 <- ggplot2::ggplot(stat, aes(mean, col = missval)) +
geom_density(na.rm = TRUE, show.legend = FALSE) +
stat_density(geom="line", position ="identity") +
labs(x = expression(log[10] ~ "Intensity"), y = "Density") +
guides(col = guide_legend(title = paste0(
"Missing values\n(",
round(
100 * nrow(cumsum[cumsum$missval == TRUE, ]) / nrow(cumsum[cumsum$missval == FALSE, ]),
digits = 1
),
"% missing)"
))) +
theme_minimal(base_size = basesize)
p2 <- ggplot2::ggplot(cumsum, aes(mean, cs_frac, col = missval)) +
geom_line() + labs(x = expression(log[10] ~ "Intensity"),
y = "Cumulative fraction") +
guides(col = guide_legend(title = paste0(
"Missing values\n(",
round(
100 * nrow(cumsum[cumsum$missval == TRUE, ]) / nrow(cumsum[cumsum$missval ==
FALSE, ]),
digits = 1
),
"% missing)"
))) +
theme_minimal(base_size = basesize)
df_per_group <- df
df_per_group$ID <- stringr::str_replace(df_per_group$ID, "_[:digit:]+$", "")
stat_per_group <- df_per_group %>% group_by(rowname, ID) %>% summarize(mean = mean(val,
na.rm = TRUE), missval = any(is.na(val)))
cumsum_per_group <- stat_per_group %>% group_by(missval) %>% arrange(mean) %>%
mutate(num = 1, cs = cumsum(num), cs_frac = cs/n())
p1_per_group <- ggplot2::ggplot(stat_per_group, aes(mean, col = missval)) +
geom_density(na.rm = TRUE, show.legend = FALSE) +
stat_density(geom = "line", position = "identity") +
labs(x = expression(log[2] ~ "Intensity"), y = "Density") +
guides(col = guide_legend(title = paste0(
"Missing values\n(",
round(
100 * nrow(cumsum_per_group[cumsum_per_group$missval == TRUE, ]) / nrow(cumsum_per_group[cumsum_per_group$missval ==
FALSE, ]),
digits = 1
),
"% missing)"
))) +
theme_minimal(base_size = basesize)
p2_per_group <-
ggplot2::ggplot(cumsum_per_group, aes(mean, cs_frac, col = missval)) +
geom_line() + labs(x = expression(log[2] ~ "Intensity"),
y = "Cumulative fraction") +
guides(col = guide_legend(title = paste0(
"Missing values\n(",
round(
100 * nrow(cumsum_per_group[cumsum_per_group$missval == TRUE, ]) / nrow(cumsum_per_group[cumsum_per_group$missval ==
FALSE, ]),
digits = 1
),
"% missing)"
))) +
theme_minimal(base_size = basesize)
if (export == TRUE){
message(paste0("Exporting density distributions and cumulative fraction of proteins with and without missing values to:\n\"", getwd(), "/", imgName, "\""))
if (format == "pdf") {
grDevices::pdf(imgName)
}
if (format == "png") {
grDevices::png(imgName,
width = 6,
height = 6,
units = 'in',
res = 300
)
}
gridExtra::grid.arrange(p1, p2, ncol = 1)
grDevices::dev.off()
mSetObj$imgSet$missval_density <- imgName
}
if (plot == TRUE){
gridExtra::grid.arrange(p1, p2, ncol = 1)
}
mSetObj$imgSet$missval_density.plot <- gridExtra::grid.arrange(p1, p2, ncol = 1)
return(mSetObj)
}
#' Visualize p values determined by ANOVA analysis.
#'
#' \code{met.plot_ANOVA} visualizes the significance of each compound in the dataset.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after ANOVA analysis (\code{\link[VisomX]{met.ANOVA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @return The input mSet object with added scatter plot with the p value for each compound analyzed. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$anova.plot)}.
#' (generated by \code{\link[ggplot2]{ggplot}}).
#' @import ggplot2
#' @references adapted from \code{\link[MetaboAnalystR]{PlotANOVA}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_ANOVA <- function (mSetObj = NA, imgName = "ANOVA_plot", format = "pdf", dpi = NULL,
width = NA, subtitle = FALSE, export = TRUE, plot=TRUE)
{
mSetObj <- mSetObj
lod <- -log10(mSetObj$analSet$aov$p.adj)
df <- data.frame(names=names(lod), values=lod, row.names = NULL)
imgName = paste(imgName,".", format,
sep = "")
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 7
}
else {
w <- width
}
h <- w * 6/9
mSetObj$imgSet$anova <- imgName
g <- ggplot2::ggplot(df, aes(x=as.numeric(row.names(df)), y=values)) +
geom_point(aes(fill=mSetObj$analSet$aov$inx.imp), color="black", shape=21, size=2, alpha=0.6) +
scale_fill_manual(values = c("#1170AA", "#EF6F6A"), labels=c("unsignificant", "significant")) +
geom_text(aes(label = names), size=2, hjust=1, vjust=-1) +
labs(x="Metabolites", y="-Log10(adjusted p-value)", title="One-way ANOVA") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.title = element_blank(), legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)) +
geom_hline(yintercept = -log10(mSetObj[["analSet"]][["aov"]][["raw.thresh"]]), linetype = "longdash",
alpha = 0.4)
if(subtitle==TRUE){
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are followed by [space]letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
g <- g + labs(subtitle=stringr::str_wrap(subtitle.text), 55)
}
if(export == TRUE){
message(paste0("Exporting anova plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
print(g)
grDevices::dev.off()
}
mSetObj$imgSet$anova.plot <- g
if(plot==TRUE){
print(g)
}
return(mSetObj)
}
#' Visualize effect of normalization on sample density distributions.
#'
#' \code{met.plot_SampleNormSummary} visualizes the density distributions of samples before and after normalization, transformation, and scaling.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after normalization (\code{\link[VisomX]{met.normalize}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param show_prenorm (Logical, \code{TRUE} or \code{FALSE}) Shall density distributions before normalization be displayed?
#' @return The input mSet object with added plot: the top is a density plot and the bottom is a box plot.
#' In a boxplot, the bottom and top of the box are always the 25th and 75th percentile (the lower and upper quartiles, or Q1 and Q3, respectively), and the band near the middle of the box is always the 50th percentile (the median or Q2). The upper whisker is located at the smaller of the maximum x value and Q3 + 1.5 x IQR (Interquantile Range), whereas the lower whisker is located at the larger of the smallest x value and Q1 - 1.5 x IQR.
#' The plot can be retrieved from within R via \code{print(mSetObj$imgSet$summary_norm_sample.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotSampleNormSummary}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_SampleNormSummary <- function (mSetObj = NA, imgName = "SampleNormSummary", format = "png", dpi = NULL,
width = NA, show_prenorm = TRUE, export = TRUE, plot=TRUE)
{
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- if_else(show_prenorm, 10, 5.5)
}
else if (width == 0) {
w <- if_else(show_prenorm, 10, 5.5)
}
else if (width > 0) {
w = width
}
h = w * if_else(show_prenorm, 1.25, 2.3)
if(!is.null(mSetObj$dataSet$prenorm)){
proc.data <- data.frame(mSetObj$dataSet$prenorm)
} else if(!is.null(mSetObj$dataSet$data_proc)){
proc.data <- mSetObj$dataSet$data_proc
} else {
proc.data <- qs::qread("data_proc.qs")
}
if(show_prenorm){
p <- function(){
graphics::layout(matrix(c(1, 2, 2, 2, 3, 4, 4, 4), 4, 2, byrow = FALSE))
pre.inx <- if (nrow(proc.data) < 50) {
1:nrow(proc.data)
} else {
sample(1:nrow(proc.data), 50)
}
namesVec <- rownames(proc.data[pre.inx, , drop = FALSE])
nm.inx <- namesVec %in% rownames(mSetObj$dataSet$norm)
namesVec <- namesVec[nm.inx]
pre.inx <- pre.inx[nm.inx]
norm.inx <- match(namesVec, rownames(mSetObj$dataSet$norm))
namesVec <- substr(namesVec, 1, 20)
rangex.pre <- range(proc.data[pre.inx, , drop = FALSE], na.rm = T)
rangex.norm <- range(mSetObj$dataSet$norm[norm.inx, , drop = FALSE],
na.rm = T)
x.label <- if_else(mSetObj$dataSet$type == "conc", "Concentration", "Intensity")
y.label <- "Samples"
op <- graphics::par(mar = c(6.5, 7, 0, 0), xaxt = "s")
plot(stats::density(apply(proc.data, 1, mean, na.rm = TRUE)), col = "darkblue",
las = 2, lwd = 2, main = "", xlab = "", ylab = "")
graphics::mtext("Density", 2, 5)
op <- graphics::par(mar = c(5.75, 8, 4, 0), xaxt = "s")
graphics::boxplot(t(proc.data[pre.inx, , drop = FALSE]), names = namesVec,
ylim = rangex.pre, las = 2, col = "lightgreen",
horizontal = T)
graphics::mtext(x.label, 1, 4.7, font = 2)
graphics::mtext(font = 2, "Before Normalization", 3, 1)
op <- graphics::par(mar = c(6.5, 7, 0, 2), xaxt = "s")
plot(stats::density(apply(mSetObj$dataSet$norm, 1, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
op <- graphics::par(mar = c(5.75, 8, 4, 2), xaxt = "s")
graphics::boxplot(t(mSetObj$dataSet$norm[norm.inx, , drop = FALSE]),
names = namesVec, ylim = rangex.norm, las = 2, col = "lightgreen",
ylab = "", horizontal = T)
subtitle.text <- paste0(
if_else(
mSetObj[["dataSet"]][["rownorm.method"]] == "N/A" &
mSetObj[["dataSet"]][["trans.method"]] == "N/A" &
mSetObj[["dataSet"]][["scale.method"]] == "N/A" ,
"(no normalization/transformation)",
stringr::str_replace_all(
paste0(
"(",
mSetObj[["dataSet"]][["rownorm.method"]],
", ",
mSetObj[["dataSet"]][["trans.method"]],
", ",
mSetObj[["dataSet"]][["scale.method"]],
")"
),
"N/A",
""
)
)
) %>% str_replace_all(., "(?<![:alpha:]), ", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., ",[:blank:](?![:alpha:])", "") # Remove all commas that are followed by [space]letter.
graphics::mtext(font = 2, paste0("After Normalization\n", stringr::str_wrap(subtitle.text, 55)), 3, 1)
graphics::mtext(font = 2, paste("Normalized", x.label), 1, 4.7)
}
} else {
p <- function(){
graphics::layout(matrix(c(1, 2, 2, 2), 4, 1, byrow = FALSE))
pre.inx <- if (nrow(proc.data) < 50) {
1:nrow(proc.data)
} else {
sample(1:nrow(proc.data), 50)
}
namesVec <- rownames(proc.data[pre.inx, , drop = FALSE])
nm.inx <- namesVec %in% rownames(mSetObj$dataSet$norm)
namesVec <- namesVec[nm.inx]
pre.inx <- pre.inx[nm.inx]
norm.inx <- match(namesVec, rownames(mSetObj$dataSet$norm))
namesVec <- substr(namesVec, 1, 20)
rangex.pre <- range(proc.data[pre.inx, , drop = FALSE], na.rm = T)
rangex.norm <- range(mSetObj$dataSet$norm[norm.inx, , drop = FALSE],
na.rm = T)
x.label <- if_else(mSetObj$dataSet$type == "conc", "Concentration", "Intensity")
y.label <- "Samples"
op <- graphics::par(mar = c(6.5, 7, 0, 0), xaxt = "s")
plot(density(apply(mSetObj$dataSet$norm, 1, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
graphics::mtext("Density", 2, 5)
op <- graphics::par(mar = c(5.75, 8, 4, 0), xaxt = "s")
graphics::boxplot(t(mSetObj$dataSet$norm[norm.inx, , drop = FALSE]),
names = namesVec, ylim = rangex.norm, las = 2, col = "lightgreen",
ylab = "", horizontal = T)
subtitle.text <- paste0(
if_else(
mSetObj[["dataSet"]][["rownorm.method"]] == "N/A" &
mSetObj[["dataSet"]][["trans.method"]] == "N/A" &
mSetObj[["dataSet"]][["scale.method"]] == "N/A" ,
"(no normalization/transformation)",
stringr::str_replace_all(
paste0(
"(",
mSetObj[["dataSet"]][["rownorm.method"]],
", ",
mSetObj[["dataSet"]][["trans.method"]],
", ",
mSetObj[["dataSet"]][["scale.method"]],
")"
),
"N/A",
""
)
)
) %>% str_replace_all(., "(?<![:alpha:]), ", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., ",[:blank:](?![:alpha:])", "") # Remove all commas that are followed by [space]letter.
graphics::mtext(font = 2, paste0("After Normalization\n", stringr::str_wrap(subtitle.text, 55)), 3, 1)
graphics::mtext(font = 2, paste("Normalized", x.label), 1, 4.7)
}
}
if(export == TRUE){
message(paste0("Exporting density distributions of samples before and after normalization to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$summary_norm_sample <- imgName
}
temp <- tempfile()
grDevices::pdf(file=temp)
p(); mSetObj$imgSet$summary_norm_sample.plot <- grDevices::recordPlot()
grDevices::dev.off()
file.remove(temp)
if(plot == TRUE){
p()
}
return(mSetObj)
}
#' Visualize effect of normalization on feature density distributions.
#'
#' \code{met.plot_FeatureNormSummary} visualizes the density distributions of features before and after normalization, transformation, and scaling.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after normalization (\code{\link[VisomX]{met.normalize}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param pre.inx (Numeric vector, or \code{NULL}) Provide index of compounds in your dataset that shall be diplayed with box plots. If \code{NULL}, the compounds with minimum and maximum values as well as up to 50 additional features are selected automatically.
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param show_prenorm (Logical, \code{TRUE} or \code{FALSE}) Shall density distributions before normalization be displayed?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @return The input mSet object with added plot: the top is a density plot and the bottom is a box plot.
#' In a boxplot, the bottom and top of the box are always the 25th and 75th percentile (the lower and upper quartiles, or Q1 and Q3, respectively), and the band near the middle of the box is always the 50th percentile (the median or Q2). The upper whisker is located at the smaller of the maximum x value and Q3 + 1.5 x IQR (Interquantile Range), whereas the lower whisker is located at the larger of the smallest x value and Q1 - 1.5 x IQR.
#' The plot can be retrieved from within R via \code{print(mSetObj$imgSet$summary_norm_feature.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotFeatureNormSummary}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_FeatureNormSummary <- function (mSetObj = NA, imgName = "FeatureNormSummary", format = "pdf", dpi = NULL, pre.inx = NULL,
width = NA, show_prenorm = TRUE, export = TRUE, plot = TRUE)
{
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- if_else(show_prenorm, 10.5, 5.75)
}
else if (width == 0) {
w = if_else(show_prenorm, 7.2, 3.6)
}
else if (width > 0) {
w = width
}
h = w * if_else(show_prenorm, 1.1905, 2)
if(!is.null(mSetObj$dataSet$prenorm)){
proc.data <- data.frame(mSetObj$dataSet$prenorm)
} else if(!is.null(mSetObj$dataSet$data_proc)){
proc.data <- mSetObj$dataSet$data_proc
} else {
proc.data <- qs::qread("data_proc.qs")
}
if(show_prenorm){
p <- function(){
graphics::layout(matrix(c(1, 2, 2, 2, 3, 4, 4, 4), 4, 2, byrow = FALSE))
if (is.null(pre.inx)) {
if (ncol(proc.data) > 50) {
min.inx <- which.min(colMeans(proc.data))
max.inx <- which.max(colMeans(proc.data))
order.inx <- order(colMeans(proc.data))
pre.inx <- unname(c(min.inx,
order.inx[!order.inx %in% c(min.inx, max.inx)][round(seq(2, (ncol(proc.data) - 2),
ncol(proc.data) / 50))],
max.inx))
} else {
pre.inx <- 1:nrow(proc.data)
}
}
namesVec <- colnames(proc.data[, pre.inx, drop = FALSE])
nm.inx <- namesVec %in% colnames(mSetObj$dataSet$norm)
namesVec <- namesVec[nm.inx]
pre.inx <- pre.inx[nm.inx]
norm.inx <- match(namesVec, colnames(mSetObj$dataSet$norm))
namesVec <- substr(namesVec, 1, 12)
rangex.pre <- range(proc.data[, pre.inx, drop = FALSE], na.rm = T)
rangex.norm <- range(mSetObj$dataSet$norm[, norm.inx, drop = FALSE],
na.rm = T)
x.label <- if_else(mSetObj$dataSet$type == "conc", "Concentration", "Intensity")
y.label <- if (mSetObj$dataSet$type == "conc") {
"Compounds"
} else if (mSetObj$dataSet$type == "specbin") {
"Spectra Bins"
} else if (mSetObj$dataSet$type == "nmrpeak") {
"Peaks (ppm)"
} else if (mSetObj$dataSet$type == "mspeak") {
"Peaks (mass)"
} else {
"Peaks(mz/rt)"
}
if (anal.type == "roc" & mSetObj$dataSet$roc_cols ==
1) {
op <- graphics::par(mar = c(4, 7, 4, 0), xaxt = "s")
graphics::plot.new()
}
else {
op <- graphics::par(mar = c(4, 7, 4, 0), xaxt = "s")
plot(density(apply(proc.data, 2, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
graphics::mtext("Density", 2, 5)
graphics::mtext(font = 2, "Before Normalization", 3, 1)
}
op <- graphics::par(mar = c(7, 7, 0.5, 0), xaxt = "s")
graphics::boxplot(proc.data[, pre.inx, drop = FALSE], names = namesVec,
ylim = rangex.pre, las = 2, col = "lightgreen",
horizontal = T, show.names = T)
graphics::mtext(x.label, 1, 5, font = 2)
if (anal.type == "roc" & mSetObj$dataSet$roc_cols ==
1) {
op <- graphics::par(mar = c(4, 7, 4, 2), xaxt = "s")
graphics::plot.new()
}
else {
op <- graphics::par(mar = c(4, 7, 4, 2), xaxt = "s")
plot(density(apply(mSetObj$dataSet$norm, 2, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
subtitle.text <- paste0(
if_else(
mSetObj[["dataSet"]][["rownorm.method"]] == "N/A" &
mSetObj[["dataSet"]][["trans.method"]] == "N/A" &
mSetObj[["dataSet"]][["scale.method"]] == "N/A" ,
"(no normalization/transformation)",
stringr::str_replace_all(
paste0(
"(",
mSetObj[["dataSet"]][["rownorm.method"]],
", ",
mSetObj[["dataSet"]][["trans.method"]],
", ",
mSetObj[["dataSet"]][["scale.method"]],
")"
),
"N/A",
""
)
)
) %>% str_replace_all(., "(?<![:alpha:]), ", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., ",[:blank:](?![:alpha:])", "") # Remove all commas that are followed by [space]letter.
graphics::mtext(paste0("After Normalization"), 3, 1, font = 2)
}
op <- graphics::par(mar = c(7, 7, 0.5, 2), xaxt = "s")
graphics::boxplot(mSetObj$dataSet$norm[, norm.inx, drop = FALSE], names = namesVec,
ylim = rangex.norm, las = 2, col = "lightgreen",
horizontal = T, show.names = T)
graphics::mtext(font = 2, paste("Normalized", x.label), 1, 5)
}
} else {
p <- function(){
graphics::layout(matrix(c(1, 2, 2, 2), 4, 1, byrow = FALSE))
if(is.null(pre.inx)){
if (ncol(proc.data) > 50) {
min.inx <- which.min(colMeans(proc.data))
max.inx <- which.max(colMeans(proc.data))
order.inx <- order(colMeans(proc.data))
pre.inx <- unname(c(min.inx,
order.inx[!order.inx %in% c(min.inx, max.inx)][round(seq(2, (ncol(proc.data) - 2),
ncol(proc.data) / 50))],
max.inx))
} else {
pre.inx <- 1:nrow(proc.data)
}
}
namesVec <- colnames(proc.data[, pre.inx, drop = FALSE])
nm.inx <- namesVec %in% colnames(mSetObj$dataSet$norm)
namesVec <- namesVec[nm.inx]
pre.inx <- pre.inx[nm.inx]
norm.inx <- match(namesVec, colnames(mSetObj$dataSet$norm))
namesVec <- substr(namesVec, 1, 12)
rangex.pre <- range(proc.data[, pre.inx, drop = FALSE], na.rm = T)
rangex.norm <- range(mSetObj$dataSet$norm[, norm.inx, drop = FALSE],
na.rm = T)
x.label <- if_else(mSetObj$dataSet$type == "conc", "Concentration", "Intensity")
y.label <- if (mSetObj$dataSet$type == "conc") {
"Compounds"
}
else if (mSetObj$dataSet$type == "specbin") {
"Spectra Bins"
}
else if (mSetObj$dataSet$type == "nmrpeak") {
"Peaks (ppm)"
}
else if (mSetObj$dataSet$type == "mspeak") {
"Peaks (mass)"
}
else {
"Peaks(mz/rt)"
}
if (anal.type == "roc" & mSetObj$dataSet$roc_cols ==
1) {
op <- graphics::par(mar = c(7, 7, 0, 0), xaxt = "s")
graphics::plot.new()
}
else {
op <- graphics::par(mar = c(7, 7, 4, 0), xaxt = "s")
plot(density(apply(mSetObj$dataSet$norm, 2, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
graphics::mtext("Density", 2, 5)
}
op <- graphics::par(mar = c(7, 7, 0, 2), xaxt = "s")
graphics::boxplot(mSetObj$dataSet$norm[, norm.inx, drop = FALSE], names = namesVec,
ylim = rangex.norm, las = 2, col = "lightgreen",
horizontal = T, show.names = T)
subtitle.text <- paste0(
if_else(
mSetObj[["dataSet"]][["rownorm.method"]] == "N/A" &
mSetObj[["dataSet"]][["trans.method"]] == "N/A" &
mSetObj[["dataSet"]][["scale.method"]] == "N/A" ,
"(no normalization/transformation)",
stringr::str_replace_all(
paste0(
"(",
mSetObj[["dataSet"]][["rownorm.method"]],
", ",
mSetObj[["dataSet"]][["trans.method"]],
", ",
mSetObj[["dataSet"]][["scale.method"]],
")"
),
"N/A",
""
)
)
) %>% str_replace_all(., "(?<![:alpha:]), ", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., ",[:blank:](?![:alpha:])", "") # Remove all commas that are followed by [space]letter.
graphics::mtext(font = 2, paste0("After Normalization\n", stringr::str_wrap(subtitle.text, 55)), 3, 1)
graphics::mtext(font = 2, paste("Normalized", x.label), 1, 4.7)
}
}
if(export == TRUE){
message(paste0("Exporting density distributions of features before and after normalization to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$summary_norm_feature <- imgName
}
temp <- tempfile()
grDevices::pdf(file=temp)
p(); mSetObj$imgSet$summary_norm_feature.plot <- grDevices::recordPlot()
grDevices::dev.off()
file.remove(temp)
if(plot == TRUE){
p()
}
return(mSetObj)
}
#' Pattern hunter, plot correlation heatmap
#'
#' \code{met.plot_CorrHeatMap_Samples} visualizes the correlations between samples.
#'
#' @param mSetObj Input name of the created mSet object.
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param cor.method (Character) Indicate the correlation method, \code{'pearson'}, \code{'spearman'}, or \code{'kendall'}.
#' @param colors (Character) Indicate the colors for the heatmap, \code{"bwm"} for default, \code{"gbr"} for red/green, \code{"heat"} for heat colors, \code{"topo"} for topo colors, and \code{"gray"} for gray scale.
#' @param viewOpt (Character) Indicate \code{"overview"} to get an overview of the heatmap, and \code{"detail"} to get a detailed view of the heat map.
#' @param fix.col (Logical) Fix colors (\code{TRUE}) or not (\code{FALSE}).
#' @param no.clst (Logical) Indicate if the correlations should be clustered (\code{TRUE}) or not (\code{FALSE}).
#' @param corrCutoff (Numeric) Set correlation cut off.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @return The input mSet object with added heat map showing the correlations between samples. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$corr.heatmap_samples.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotCorrHeatMap}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_CorrHeatMap_Samples <- function (mSetObj = NA, imgName = "correlation_samples", format = "pdf", dpi = 72,
width = NA, cor.method = "pearson", colors="bwm", viewOpt="overview", fix.col = FALSE,
no.clst = FALSE, corrCutoff = 0, plot = TRUE, export = TRUE)
{
target = "row"
main <- xlab <- ylab <- NULL
data <- mSetObj$dataSet$norm
corrCutoff <- as.numeric(corrCutoff)
if (target == "row") {
data <- t(data)
}
if (ncol(data) > 1000) {
filter.val <- apply(data.matrix(data), 2, stats::IQR, na.rm = T)
rk <- rank(-filter.val, ties.method = "random")
data <- as.data.frame(data[, rk <= 1000])
print("Data is reduced to 1000 vars ..")
}
colnames(data) <- substr(colnames(data), 1, 18)
corr.mat <- cor(data, method = cor.method)
corr.mat[abs(corr.mat) < corrCutoff] <- 0
mSetObj$analSet$pwcor <- list()
mSetObj$analSet$pwcor$data <- data
mSetObj$analSet$pwcor$cor.method <- cor.method
mSetObj$analSet$pwcor$no.clst <- no.clst
if (colors == "gbr") {
colors <- grDevices::colorRampPalette(c("green", "black",
"red"), space = "rgb")(256)
}
else if (colors == "heat") {
colors <- grDevices::heat.colors(256)
}
else if (colors == "topo") {
colors <- grDevices::topo.colors(256)
}
else if (colors == "gray") {
colors <- grDevices::colorRampPalette(c("grey90", "grey10"))(256)
}
else {
colors <- rev(grDevices::colorRampPalette(c(RColorBrewer::brewer.pal(10,
"RdBu"), "#001833"))(256))
}
imgName = paste(imgName, ".", format,
sep = "")
if (viewOpt == "overview") {
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 7.2
}
else {
w <- 7.2
}
h <- w
}
else {
if (ncol(corr.mat) > 50) {
myH <- ncol(corr.mat) * 12 + 40
}
else if (ncol(corr.mat) > 20) {
myH <- ncol(corr.mat) * 12 + 60
}
else {
myH <- ncol(corr.mat) * 12 + 120
}
h <- round(myH/72, 2)
if (is.na(width)) {
w <- h
}
else if (width == 0) {
w <- h <- 7.2
}
else {
w <- h <- 7.2
}
}
if (no.clst) {
rowv = FALSE
colv = FALSE
dendro = "none"
}
else {
rowv = TRUE
colv = TRUE
dendro = "both"
}
if (fix.col) {
breaks <- seq(from = -1,
to = 1,
length = 257)
res <-
pheatmap::pheatmap(
corr.mat,
fontsize = 11,
fontsize_row = 11,
cluster_rows = colv,
cluster_cols = rowv,
color = colors,
breaks = breaks,
silent=T
)
}
else {
res <-
pheatmap::pheatmap(
corr.mat,
fontsize = 11,
fontsize_row = 11,
cluster_rows = colv,
cluster_cols = rowv,
color = colors,
silent=T
)
}
if (export == TRUE) {
message(paste0("Exporting correlation heat map of samples to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
grid::grid.newpage()
grid::grid.draw(res$gtable)
grDevices::dev.off()
mSetObj$imgSet$corr.heatmap_samples <- imgName
}
if (!no.clst) {
new.ord <- res$tree_row$order
corr.mat <- corr.mat[new.ord, new.ord]
mSetObj$analSet$pwcor$new.ord <- new.ord
}
if(plot == TRUE){
grid::grid.newpage()
grid::grid.draw(res$gtable)
}
mSetObj$imgSet$corr.heatmap_samples.plot <- res
message(paste0("Exporting CSV table file with sample correlations to:\n\"", getwd(), "/correlation_table_samples.csv\""))
fast.write.csv(signif(corr.mat, 5), file = "correlation_table_samples.csv")
return(mSetObj)
}
#' Pattern hunter, plot correlation heatmap
#'
#' \code{met.plot_CorrHeatMap_Features} visualizes the correlations between features.
#'
#' @param mSetObj Input name of the created mSet object.
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param cor.method Indicate the correlation method, 'pearson', 'spearman', or 'kendall'.
#' @param colors Indicate the colors for the heatmap, "bwm" for default, "gbr" for red/green, "heat" for heat colors, "topo" for topo colors, and "gray" for gray scale.
#' @param viewOpt Indicate "overview" to get an overview of the heatmap, and "detail" to get a detailed view of the heatmap.
#' @param fix.col (Logical) Fix colors (\code{TRUE}) or not (\code{FALSE}).
#' @param no.clst (Logical) Indicate if the correlations should be clustered (\code{TRUE}) or not (\code{FALSE}).
#' @param corrCutoff set correlation cut off
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @return The input mSet object with added heat map showing the correlations between features. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$corr.heatmap_features.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotCorrHeatMap}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
#' @export
met.plot_CorrHeatMap_Features <- function (mSetObj = NA, imgName = "correlation_features", format = "pdf", dpi = NULL,
width = NA, cor.method = "pearson", colors="bwm", viewOpt="overview", fix.col = TRUE,
no.clst = FALSE, corrCutoff = 0, plot = TRUE, export = TRUE)
{
target = "col"
main <- xlab <- ylab <- NULL
data <- mSetObj$dataSet$norm
corrCutoff <- as.numeric(corrCutoff)
if (target == "row") {
data <- t(data)
}
if (ncol(data) > 1000) {
filter.val <- apply(data.matrix(data), 2, stats::IQR, na.rm = T)
rk <- rank(-filter.val, ties.method = "random")
data <- as.data.frame(data[, rk <= 1000])
print("Data is reduced to 1000 vars ..")
}
colnames(data) <- substr(colnames(data), 1, 18)
corr.mat <- cor(data, method = cor.method)
corr.mat[abs(corr.mat) < corrCutoff] <- 0
mSetObj$analSet$pwcor <- list()
mSetObj$analSet$pwcor$data <- data
mSetObj$analSet$pwcor$cor.method <- cor.method
mSetObj$analSet$pwcor$no.clst <- no.clst
if (colors == "gbr") {
colors <- grDevices::colorRampPalette(c("green", "black",
"red"), space = "rgb")(256)
}
else if (colors == "heat") {
colors <- grDevices::heat.colors(256)
}
else if (colors == "topo") {
colors <- grDevices::topo.colors(256)
}
else if (colors == "gray") {
colors <- grDevices::colorRampPalette(c("grey90", "grey10"))(256)
}
else {
colors <- rev(grDevices::colorRampPalette(c(RColorBrewer::brewer.pal(10,
"RdBu"), "#001833"))(256))
}
imgName = paste(imgName, ".", format,
sep = "")
if (viewOpt == "overview") {
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 7.2
}
else {
w <- 7.2
}
h <- w
}
else {
if (ncol(corr.mat) > 50) {
myH <- ncol(corr.mat) * 12 + 40
}
else if (ncol(corr.mat) > 20) {
myH <- ncol(corr.mat) * 12 + 60
}
else {
myH <- ncol(corr.mat) * 12 + 120
}
h <- round(myH/72, 2)
if (is.na(width)) {
w <- h
}
else if (width == 0) {
w <- h <- 7.2
}
else {
w <- h <- 7.2
}
}
if (no.clst) {
rowv = FALSE
colv = FALSE
dendro = "none"
}
else {
rowv = TRUE
colv = TRUE
dendro = "both"
}
if (fix.col) {
breaks <- seq(from = -1,
to = 1,
length = 257)
res <-
pheatmap::pheatmap(
corr.mat,
fontsize = 4 ^ (2 - (ncol(mSetObj$dataSet$norm) / 110)),
fontsize_row = 4 ^ (2 - (ncol(mSetObj$dataSet$norm) / 110)),
cluster_rows = colv,
cluster_cols = rowv,
color = colors,
border_color = NA,
breaks = breaks,
silent=T
)
}
else {
res <-
pheatmap::pheatmap(
corr.mat,
fontsize = 4 ^ (2 - (ncol(mSetObj$dataSet$norm) / 110)),
fontsize_row = 4 ^ (2 - (ncol(mSetObj$dataSet$norm) / 110)),
cluster_rows = colv,
cluster_cols = rowv,
color = colors,
border_color = NA,
silent=T
)
}
if (export == TRUE) {
message(paste0("Exporting correlation heat map of features to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
grid::grid.newpage()
grid::grid.draw(res$gtable)
grDevices::dev.off()
mSetObj$imgSet$corr.heatmap_features <- res
}
if (!no.clst) {
new.ord <- res$tree_row$order
corr.mat <- corr.mat[new.ord, new.ord]
mSetObj$analSet$pwcor$new.ord <- new.ord
}
if(plot == TRUE){
grid::grid.newpage()
grid::grid.draw(res$gtable)
}
mSetObj$imgSet$corr.heatmap_features.plot <- res
message(paste0("Exporting CSV table file with feature correlations to:\n\"", getwd(), "/correlation_table_features.csv\""))
fast.write.csv(signif(corr.mat, 5), file = "correlation_table_features.csv")
return(mSetObj)
}
#' Plot PCA scree plot
#'
#' \code{met.plot_PCAScree} visualizes the proportion of variance explained for each principal component.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param scree.num Numeric, input a number to indicate the number of principal components to display in the scree plot.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added line plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pca.scree.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPCAScree}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
#' @export
met.plot_PCAScree <- function (mSetObj = NA, imgName = "PCA_ScreePlot", format = "pdf", dpi = NULL,
width = NA, scree.num, plot = TRUE, export = FALSE)
{
stds <- mSetObj$analSet$pca$std[1:scree.num]
pcvars <- mSetObj$analSet$pca$variance[1:scree.num]
cumvars <- mSetObj$analSet$pca$cum.var[1:scree.num]
ylims <- range(c(pcvars, cumvars))
extd <- (ylims[2] - ylims[1])/10
miny <- ifelse(ylims[1] - extd > 0, ylims[1] - extd, 0)
maxy <- ifelse(ylims[2] + extd > 1, 1, ylims[2] + extd)
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 10
}
else if (width == 0) {
w <- 8
}
else {
w <- width
}
h <- w * 2/3
p <- function(){
graphics::par(mar = c(5, 5, 6, 3))
plot(pcvars, type = "l", col = "blue", main = "Scree plot",
xlab = "PC index", ylab = "Variance explained",
ylim = c(miny, maxy), axes = F)
graphics::text(pcvars, labels = paste(100 * round(pcvars, 3), "%"),
adj = c(-0.3, -0.5), srt = 45, xpd = T)
graphics::points(pcvars, col = "red")
graphics::lines(cumvars, type = "l", col = "green")
graphics::text(cumvars, labels = paste(100 * round(cumvars, 3), "%"),
adj = c(-0.3, -0.5), srt = 45, xpd = T)
graphics::points(cumvars, col = "red")
graphics::abline(v = 1:scree.num, lty = 3)
graphics::axis(2)
graphics::axis(1, 1:length(pcvars), 1:length(pcvars))
}
if (export == TRUE) {
message(paste0("Exporting PCA scree plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$pca.scree <- imgName
}
if (plot == TRUE){
p()
}
p(); mSetObj$imgSet$pca.scree.plot <- grDevices::recordPlot()
return(mSetObj)
}
#' Create 2D PCA score plot
#'
#' \code{met.plot_PCA2DScore} visualizes clusters of samples based on their similarity in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param pcx Specify the principal component on the x-axis
#' @param pcy Specify the principal component on the y-axis
#' @param reg Numeric, input a number between 0 and 1, 0.95 will display the 95 percent confidence regions, and 0 will not.
#' @param show Display sample names, \code{1} = show names, \code{0} = do not show names.
#' @param grey.scale Use grey-scale colors, \code{1} = grey-scale, \code{0} = not grey-scale.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pca.score2d_PCx_PCy.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPCA2DScore}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PCA2DScore <- function (mSetObj = NA, imgName = "PCA_2DScores", format = "pdf", dpi = NULL, subtitle = FALSE,
width = NA, pcx, pcy, reg = 0.95, show = 1, grey.scale = 0, plot = TRUE, export = TRUE)
{
cls <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
xlabel = paste("PC", pcx, "(", round(100 * mSetObj$analSet$pca$variance[pcx],
1), "%)")
ylabel = paste("PC", pcy, "(", round(100 * mSetObj$analSet$pca$variance[pcy],
1), "%)")
pc1 = mSetObj$analSet$pca$x[, pcx]
pc2 = mSetObj$analSet$pca$x[, pcy]
text.lbls <- substr(names(pc1), 1, 14)
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 7.2
}
else if (width == 0) {
w <- 7.2
}
else {
w <- width
}
if (exists("group_factor")) {
legend.nm <- unique(as.character(sort(group_factor)))
}
else{
legend.nm <- unique(as.character(sort(cls)))
}
h <- w-1
p <- function(){
graphics::plot.new()
op <- graphics::par(mar=c(5,5,if_else(subtitle==TRUE, 7.9, 3),3))
l <- graphics::legend(0, 0, bty='n', legend = legend.nm,
plot=FALSE, pch=c(1, 2), lty=c(1, 2))
# calculate right margin width in ndc
w <- graphics::grconvertX(l$rect$w, to='ndc') - graphics::grconvertX(0, to='ndc')
graphics::par(omd=c(0, 1-w, 0, 1))
if (cls.type == "disc") {
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.factor(as.numeric(levels(cls))[cls])
}
else {
cls <- cls
}
lvs <- levels(cls)
pts.array <- array(0, dim = c(100, 2, length(lvs)))
for (i in 1:length(lvs)) {
inx <- cls == lvs[i]
groupVar <- stats::var(cbind(pc1[inx], pc2[inx]), na.rm = T)
groupMean <- cbind(mean(pc1[inx], na.rm = T), mean(pc2[inx],
na.rm = T))
pts.array[, , i] <- ellipse::ellipse(groupVar, centre = groupMean,
level = reg, npoints = 100)
}
xrg <- range(pc1, pts.array[, 1, ])
yrg <- range(pc2, pts.array[, 2, ])
x.ext <- (xrg[2] - xrg[1])/12
y.ext <- (yrg[2] - yrg[1])/12
xlims <- c(xrg[1] - x.ext, xrg[2] + x.ext)
ylims <- c(yrg[1] - y.ext, yrg[2] + y.ext)
cols <- GetColorSchema(cls, grayscale = F)
uniq.cols <- unique(cols)
plot(pc1, pc2, xlab = xlabel, xlim = xlims, ylim = ylims,
ylab = ylabel, type = "n",
col = cols, pch = as.numeric(cls) + 1,
cex.lab = 2,
cex.axis = 2,
cex.main = 2,
cex.sub = 2)
graphics::title("PCA Scores Plot", line = if_else(subtitle==TRUE, 4.5, 1), cex.main = 2.5)
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::mtext(line = 3.9, padj = 1, cex=1.7, if_else(subtitle==TRUE, stringr::str_wrap(subtitle.text, 55), NULL))
graphics::grid(col = "lightgray", lty = "dotted", lwd = 1)
if (length(uniq.cols) > 1) {
names(uniq.cols) <- unique(as.character(sort(cls)))
}
for (i in 1:length(lvs)) {
if (length(uniq.cols) > 1) {
graphics::polygon(pts.array[, , i], col = grDevices::adjustcolor(uniq.cols[lvs[i]],
alpha = 0.2), border = NA)
}
else {
graphics::polygon(pts.array[, , i], col = grDevices::adjustcolor(uniq.cols,
alpha = 0.2), border = NA)
}
if (grey.scale) {
graphics::lines(pts.array[, , i], col = grDevices::adjustcolor("black",
alpha = 0.5), lty = 2)
}
}
pchs <- GetShapeSchema(mSetObj, show, grey.scale)
if (grey.scale) {
cols <- rep("black", length(cols))
}
if (show == 1) {
graphics::text(pc1, pc2, label = text.lbls, pos = 4, xpd = T,
cex = 0.75)
graphics::points(pc1, pc2, pch = pchs, col = cols)
}
else {
if (length(uniq.cols) == 1) {
graphics::points(pc1, pc2, pch = pchs, col = cols, cex = 1)
}
else {
if (grey.scale == 1 | (exists("shapeVec") &&
all(shapeVec >= 0))) {
my.cols <- grDevices::adjustcolor(cols, alpha.f = 0.4)
my.cols[pchs == 21] <- "black"
graphics::points(pc1, pc2, pch = pchs, col = my.cols,
bg = grDevices::adjustcolor(cols, alpha.f = 0.4), cex = 1.8)
}
else {
graphics::points(pc1, pc2, pch = 21, bg = grDevices::adjustcolor(cols,
alpha.f = 0.4), cex = 2)
}
}
}
uniq.pchs <- unique(pchs)
if (length(uniq.cols) != length(levels(cls))) {
if (mSetObj$dataSet$type.cls.lbl == "integer") {
names(cols) <- as.numeric(cls)
}
else {
names(cols) <- as.character(cls)
}
match.inx <- match(levels(cls), names(cols))
uniq.cols <- cols[match.inx]
}
if (grey.scale) {
uniq.cols <- "black"
}
graphics::legend(graphics::par('usr')[2], graphics::par('usr')[4], bty='n', legend = legend.nm, pch = uniq.pchs,
col = uniq.cols, cex = 1.3, xpd = NA)
}
else {
plot(pc1, pc2, xlab = xlabel, ylab = ylabel, type = "n",
main = "Scores Plot")
graphics::points(pc1, pc2, pch = 15, col = "magenta")
graphics::text(pc1, pc2, label = text.lbls, pos = 4, col = "blue",
xpd = T, cex = 0.8)
}
graphics::par(op)
}
if (export == TRUE){
message(paste0("Exporting 2D PCA scores plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[['imgSet']][[paste0('pca.score2d_PC', pcx, "_PC", pcy)]] <- imgName
}
p(); mSetObj[['imgSet']][[paste0('pca.score2d_PC', pcx, "_PC", pcy, ".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
if(plot==TRUE){
p()
}
return(mSetObj)
}
#' Plot PCA loadings and also set up the matrix for display
#'
#' \code{met.plot_PCA2DLoading} visualizes which features are responsible for the patterns seen among the samples in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param inx1 Numeric, indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 Numeric, indicate the number of the principal component for the y-axis of the loading plot.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pca.loading_PCx_PCy.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPCALoading}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PCA2DLoading <- function (mSetObj = NA, imgName = "PCA_2DLoadings", format = "pdf", dpi = NULL, subtitle = FALSE,
width = NA, inx1, inx2, export = TRUE, plot=TRUE)
{
loadings <- as.matrix(cbind(mSetObj$analSet$pca$rotation[,
inx1], mSetObj$analSet$pca$rotation[, inx2]))
ord.inx <- order(-abs(loadings[, 1]), -abs(loadings[, 2]))
loadings <- signif(loadings[ord.inx, ], 5)
ldName1 = paste("Loadings PC", inx1, "(", round(100 * mSetObj$analSet$pca$variance[inx1],
1), "%)")
ldName2 = paste("Loadings PC", inx2, "(", round(100 * mSetObj$analSet$pca$variance[inx2],
1), "%)")
colnames(loadings) <- c(ldName1, ldName2)
mSetObj$analSet$pca$imp.loads <- loadings
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 7.2
}
else if (width == 0) {
w <- 7.2
}
else {
w <- width
}
h <- w-0.5
plotType <- mSetObj$analSet$pca$loading.type
p <- function(){
graphics::plot.new()
graphics::par(mar = c(6, 6.2, if_else(subtitle==TRUE, 7.9, 3), 6), mgp=c(4.5, 1, 0))
plot(loadings[, 1], loadings[, 2], las = 1, xlab = ldName1,
ylab = ldName2,
cex.lab = 2,
cex.axis = 2,
cex.main = 2,
cex.sub = 2)
graphics::title("PCA Loading Scatterplot", line = if_else(subtitle==TRUE, 4.5, 1), cex.main = 2.5)
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::mtext(line = 3.9, padj = 1, cex=1.7, if_else(subtitle==TRUE, stringr::str_wrap(subtitle.text, 55), NULL))
mSetObj$pca.axis.lims <- graphics::par("usr")
graphics::grid(col = "lightgray", lty = "dotted", lwd = 1)
graphics::points(loadings[, 1], loadings[, 2], pch = 19, col = grDevices::adjustcolor("#339933",
alpha.f = 0.4))
if (plotType == "all") {
graphics::text(loadings[, 1], loadings[, 2], labels = substr(rownames(loadings),
1, 16), pos = 4, col = "gray30", xpd = T)
}
else if (plotType == "gray30") {
if (length(mSetObj$custom.cmpds) > 0) {
hit.inx <- rownames(loadings) %in% mSetObj$custom.cmpds
graphics::text(loadings[hit.inx, 1], loadings[hit.inx, 2],
labels = rownames(loadings)[hit.inx], pos = 4,
col = "blue", xpd = T)
}
}
}
if (export == TRUE){
message(paste0("Exporting 2D PCA loadings plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[['imgSet']][[paste0('pca.loading_PC', inx1, "_PC", inx2)]] <- imgName
}
p(); mSetObj[['imgSet']][[paste0('pca.loading_PC', inx1, "_PC", inx2, ".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
if (plot == TRUE){
p()
}
return(mSetObj)
}
#' Create 2D PLS-DA score plot
#'
#' \code{met.plot_PLS2DScore} visualizes clusters of samples based on their similarity in partial least squares-discriminant analysis.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the loading plot.
#' @param reg (Numeric) Enter a number between 0 and 1, 0.95 will display the 95 percent confidence regions, and 0 will not.
#' @param show Display sample names, \code{1} = show names, \code{0} = do not show names.
#' @param grey.scale Use grey-scale colors, \code{1} = grey-scale, \code{0} = not grey-scale.
#' @param use.sparse (Logical) Use a sparse algorithm (\code{TRUE}) or not (\code{FALSE}).
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.score2d_PCx_PCy.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLS2DScore}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS2DScore <- function (mSetObj = NA, imgName = "PLSDA_2DScore", format = "pdf", dpi = NULL,
width = NA, inx1=1, inx2=2, reg = 0.95, show = 1, grey.scale = 0, subtitle=FALSE,
use.sparse = FALSE, plot = TRUE, export = TRUE)
{
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 7.2
}
else if (width == 0) {
w <- 7.2
}
else {
w <- width
}
h <- w-1
cls1 <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
lv1 <- mSetObj$analSet$plsr$scores[, inx1]
lv2 <- mSetObj$analSet$plsr$scores[, inx2]
xlabel <- paste("Component", inx1, "(", round(100 *
mSetObj$analSet$plsr$Xvar[inx1]/mSetObj$analSet$plsr$Xtotvar,
1), "%)")
ylabel <- paste("Component", inx2, "(", round(100 *
mSetObj$analSet$plsr$Xvar[inx2]/mSetObj$analSet$plsr$Xtotvar,
1), "%)")
if (exists("group_factor")) {
legend.nm <- unique(as.character(sort(group_factor)))
}
else{
legend.nm <- unique(as.character(sort(cls1)))
}
p<- function(){
graphics::plot.new()
l <- graphics::legend(0, 0, bty='n', legend = legend.nm,
plot=FALSE, pch=c(1, 2), lty=c(1, 2))
# calculate right margin width in ndc
w <- graphics::grconvertX(l$rect$w, to='ndc') - graphics::grconvertX(0, to='ndc')
graphics::par(omd=c(0, 1-w, 0, 1))
graphics::par(mar = c(5, 5, if_else(subtitle==TRUE, 7.9, 3), 3))
text.lbls <- substr(rownames(mSetObj$dataSet$norm), 1, 12)
if (cls.type == "integer") {
cls <- as.factor(as.numeric(levels(cls1))[cls1])
}
else {
cls <- cls1
}
lvs <- levels(cls)
pts.array <- array(0, dim = c(100, 2, length(lvs)))
for (i in 1:length(lvs)) {
inx <- cls1 == lvs[i]
groupVar <- stats::var(cbind(lv1[inx], lv2[inx]), na.rm = T)
groupMean <- cbind(mean(lv1[inx], na.rm = T), mean(lv2[inx],
na.rm = T))
pts.array[, , i] <- ellipse::ellipse(groupVar, centre = groupMean,
level = reg, npoints = 100)
}
xrg <- range(lv1, pts.array[, 1, ])
yrg <- range(lv2, pts.array[, 2, ])
x.ext <- (xrg[2] - xrg[1])/12
y.ext <- (yrg[2] - yrg[1])/12
xlims <- c(xrg[1] - x.ext, xrg[2] + x.ext)
ylims <- c(yrg[1] - y.ext, yrg[2] + y.ext)
cols <- GetColorSchema(cls, grayscale = F)
uniq.cols <- unique(cols)
plot(lv1, lv2, xlab = xlabel, xlim = xlims, ylim = ylims,
ylab = ylabel, type = "n",
cex.lab = 2,
cex.axis = 2,
cex.main = 2,
cex.sub = 2)
graphics::title("PLS-DA Scores Plot", line = if_else(subtitle==TRUE, 4.5, 1), cex.main = 2.5)
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::mtext(line = 3.9, padj = 1, cex=1.7, if_else(subtitle==TRUE, stringr::str_wrap(subtitle.text, 55), NULL))
graphics::grid(col = "lightgray", lty = "dotted", lwd = 1)
if (length(uniq.cols) > 1) {
names(uniq.cols) <- legend.nm
}
for (i in 1:length(lvs)) {
if (length(uniq.cols) > 1) {
graphics::polygon(pts.array[, , i], col = grDevices::adjustcolor(uniq.cols[lvs[i]],
alpha = 0.2), border = NA)
}
else {
graphics::polygon(pts.array[, , i], col = grDevices::adjustcolor(uniq.cols,
alpha = 0.2), border = NA)
}
if (grey.scale) {
graphics::lines(pts.array[, , i], col = grDevices::adjustcolor("black",
alpha = 0.5), lty = 2)
}
}
pchs <- GetShapeSchema(mSetObj, show, grey.scale)
if (grey.scale) {
cols <- rep("black", length(cols))
}
if (show == 1) {
graphics::text(lv1, lv2, label = text.lbls, pos = 4, xpd = T, cex = 0.75)
graphics::points(lv1, lv2, pch = pchs, col = cols)
}
else {
if (length(uniq.cols) == 1) {
graphics::points(lv1, lv2, pch = pchs, col = cols, cex = 1)
}
else {
if (grey.scale == 1 | (exists("shapeVec") &&
all(shapeVec >= 0))) {
my.cols <- grDevices::adjustcolor(cols, alpha.f = 0.4)
my.cols[pchs == 21] <- "black"
graphics::points(lv1, lv2, pch = pchs, col = my.cols, bg = grDevices::adjustcolor(cols,
alpha.f = 0.4), cex = 1.8)
}
else {
graphics::points(lv1, lv2, pch = 21, bg = grDevices::adjustcolor(cols,
alpha.f = 0.4), cex = 2)
}
}
}
uniq.pchs <- unique(pchs)
if (grey.scale) {
uniq.cols <- "black"
}
if (length(uniq.cols) != length(levels(cls))) {
if (cls.type == "integer") {
names(cols) <- as.numeric(cls)
}
else {
names(cols) <- as.character(cls)
}
match.inx <- match(levels(cls), names(cols))
uniq.cols <- cols[match.inx]
}
if (length(uniq.pchs) != length(levels(cls))) {
names(pchs) <- as.character(cls)
match.inx <- match(levels(cls), names(pchs))
uniq.pchs <- pchs[match.inx]
}
graphics::legend(graphics::par('usr')[2], graphics::par('usr')[4], bty='n', legend = legend.nm, pch = uniq.pchs,
col = uniq.cols, cex = 1.3, xpd = NA)
graphics::par(graphics::par(mar=c(5,5,3,3)))
}
if (export == TRUE) {
message(paste0("Exporting 2D PLS-DA scores plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[['imgSet']][[paste0('pls.score2d_PC', inx1, "_PC", inx2)]] <- imgName
}
p(); mSetObj[['imgSet']][[paste0('pls.score2d_PC', inx1, "_PC", inx2, ".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
if (plot == TRUE){
p()
}
return(mSetObj)
}
#' Plot PLS-DA loadings and also set up the matrix for display
#'
#' \code{met.plot_PLS2DLoading} visualizes which features are responsible for the patterns seen among the samples in partial least squares-discriminant analysis.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param inx1 Numeric, indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 Numeric, indicate the number of the principal component for the y-axis of the loading plot.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.loading_PCx_PCy.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLSLoading}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS2DLoading <- function (mSetObj = NA, imgName = "PLSDA_2DLoadings", format = "pdf", dpi = NULL, subtitle=FALSE,
width = NA, inx1 = 1, inx2 = 2, plot = TRUE, export = TRUE)
{
load1 <- mSetObj$analSet$plsr$loadings[, inx1]
load2 <- mSetObj$analSet$plsr$loadings[, inx2]
loadings = as.matrix(cbind(load1, load2))
ord.inx <- order(-abs(loadings[, 1]), -abs(loadings[, 2]))
loadings <- signif(loadings[ord.inx, ], 5)
ldName1 <- paste(
"Loadings Comp. ",
inx1,
"(",
round(
100 *
mSetObj$analSet$plsr$Xvar[inx1] /
mSetObj$analSet$plsr$Xtotvar,
1
),
"%)"
)
ldName2 <- paste(
"Loadings Comp. ",
inx2,
"(",
round(
100 *
mSetObj$analSet$plsr$Xvar[inx2] /
mSetObj$analSet$plsr$Xtotvar,
1
),
"%)"
)
colnames(loadings) <- c(ldName1, ldName2)
mSetObj$analSet$plsr$imp.loads <- loadings
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 10
}
else if (width == 0) {
w <- 10
}
else {
w <- width
}
h <- 0.85*w
plotType <- mSetObj$analSet$plsr$loading.type
p <- function(){
graphics::par(mar = c(6, 6, if_else(subtitle==TRUE, 7.9, 3), 6), mgp=c(4.5, 1, 0))
plot(loadings[, 1], loadings[, 2], las = 2, xlab = ldName1,
ylab = ldName2,
cex.lab = 2,
cex.axis = 2,
cex.main = 2,
cex.sub = 2)
graphics::title("PLS-DA Loading Scatterplot", line = if_else(subtitle==TRUE, 5, 1), cex.main = 2.3)
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::mtext(line = 3.9, padj = 1, cex=1.7, if_else(subtitle==TRUE, stringr::str_wrap(subtitle.text, 55), NULL))
mSetObj$pls.axis.lims <- graphics::par("usr")
graphics::grid(col = "lightgray", lty = "dotted", lwd = 1)
graphics::points(loadings[, 1], loadings[, 2], pch = 19, col = grDevices::adjustcolor("#CC6663",
alpha.f = 0.4))
if (plotType == "all") {
graphics::text(loadings[, 1], loadings[, 2], labels = substr(rownames(loadings),
1, 16), pos = 4, col = "gray30", xpd = T)
}
else if (plotType == "custom") {
if (length(mSetObj$custom.cmpds) > 0) {
hit.inx <- colnames(mSetObj$dataSet$norm) %in% mSetObj$custom.cmpds
graphics::text(loadings[hit.inx, 1], loadings[hit.inx, 2],
labels = rownames(loadings)[hit.inx], pos = 4,
col = "gray30", xpd = T)
}
}
}
if (export == TRUE) {
message(paste0("Exporting 2D PLS-DA loadings plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[['imgSet']][[paste0('pls.loading_PC', inx1, "_PC", inx2)]] <- imgName
}
p(); mSetObj[['imgSet']][[paste0('pls.loading_PC', inx1, "_PC", inx2, ".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
if (plot == TRUE){
p()
}
return(mSetObj)
}
#' Plot PLS important features
#'
#' \code{met.plot_PLS_Imp} visualizes the features' importance in the PLS-DA model.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (\code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param type (Character) Indicate the variables of importance type to use, \code{"vip"} to use VIP scores, or \code{"coef"} for coefficients
#' @param feat.nm (Character) Indicate the name of the feature. If \code{type = "vip"}, choose \code{"Comp. 1"}, \code{"Comp. 2"}, etc., depending on the component for which the VIP scores should be shown. If \code{type = "coef"}, choose \code{"coef.mean"} for average coefficients, or the name of a sample group in your dataset.
#' @param feat.num (Numeric) Indicate the number of features to show in the plot.
#' @param color.BW (Logical) \code{TRUE} to use black and white, or \code{FALSE} to not.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param title (Logical) \code{TRUE} to add a title with the used normalization, transformation, and scaling method, or \code{FALSE} to not add any title.
#' @return The input mSet object with added dot plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.imp_type_feat.nm.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLS.Impg}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS_Imp <- function (mSetObj = NA, imgName, format = "pdf", dpi = NULL,
width = NA, type = "vip", feat.nm = c("Comp. 1", "Comp. 2"), feat.num = 15,
color.BW = FALSE, plot = TRUE, export = TRUE, title=FALSE)
{
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 8
}
else if (width == 0) {
w <- 7
}
else {
w <- width
}
h <- w
p <- function(){
if (type == "vip") {
mSetObj$analSet$plsda$imp.type <- "vip"
vips <- mSetObj$analSet$plsda$vip.mat[, feat.nm]
met.plot_ImpVar(mSetObj, vips, paste0("VIP scores - ", feat.nm), feat.num,
color.BW, title = title)
}
else {
mSetObj$analSet$plsda$imp.type <- "coef"
data <- mSetObj$analSet$plsda$coef.mat[, feat.nm]
met.plot_ImpVar(mSetObj, data, paste0("Coefficients - ", feat.nm), feat.num,
color.BW, title = title)
}
}
if (export == TRUE) {
message(paste0("Exporting PLS important features plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[["imgSet"]][[paste0("pls.imp_", type, "_", feat.nm)]] <- imgName
}
p(); mSetObj[["imgSet"]][[paste0("pls.imp_", type, "_", feat.nm,".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
return(mSetObj)
}
#' Plot PLS important features
#'
#' \code{met.plot_ImpVar} is an internal function used by \code{met.plot_PLS_Imp} to visualize the features' importance in the PLS-DA model.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imp.vec A vector of importan variables.
#' @param xlbl The x-label
#' @param feat.num (Numeric) Indicate the number of features to show in the plot.
#' @param color.BW (Logical) \code{TRUE} to use black and white, or \code{FALSE} to not.
#' @param title (Logical) \code{TRUE} to add a title with the used normalization, transformation, and scaling method, or \code{FALSE} to not add any title.
#' @return The input mSet object with added plot object.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotImpVar}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
met.plot_ImpVar <- function (mSetObj = NA, imp.vec, xlbl, feat.num = 15, color.BW = FALSE, title=FALSE)
{
cls.len <- length(levels(mSetObj$dataSet$cls))
if (cls.len == 2) {
rt.mrg <- 5
}
else if (cls.len == 3) {
rt.mrg <- 6
}
else if (cls.len == 4) {
rt.mrg <- 7
}
else if (cls.len == 5) {
rt.mrg <- 8
}
else if (cls.len == 6) {
rt.mrg <- 9
}
else {
rt.mrg <- 11
}
op <- graphics::par(mar = c(5, 7, if_else(title==TRUE, 4, 3), rt.mrg))
if (feat.num <= 0) {
feat.num = 15
}
if (feat.num > length(imp.vec)) {
feat.num <- length(imp.vec)
}
imp.vec <- rev(sort(imp.vec))[1:feat.num]
imp.vec <- sort(imp.vec)
mns <- by(mSetObj$dataSet$norm[, names(imp.vec)], mSetObj$dataSet$cls,
function(x) {
apply(x, 2, mean, trim = 0.1)
})
mns <- t(matrix(unlist(mns), ncol = feat.num, byrow = TRUE))
vip.nms <- substr(names(imp.vec), 1, 14)
names(imp.vec) <- NULL
dotcolor <- ifelse(color.BW, "darkgrey", "#585855")
graphics::dotchart(imp.vec, bg = dotcolor, xlab = xlbl, cex = 1.3)
graphics::mtext(side = 2, at = 1:feat.num, vip.nms, las = 2, line = 1)
if(title==TRUE){
title.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::title(stringr::str_wrap(title.text, 55), line = 1, cex.main = 1.3)
}
axis.lims <- graphics::par("usr")
shift <- 2 * graphics::par("cxy")[1]
lgd.x <- axis.lims[2] + shift
x <- rep(lgd.x, feat.num)
y <- 1:feat.num
graphics::par(xpd = T)
nc <- ncol(mns)
colorpalette <- ifelse(color.BW, "Greys", "RdYlBu")
col <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(7, colorpalette))(nc)
if (color.BW)
col <- rev(col)
bg <- matrix("", nrow(mns), nc)
for (m in 1:nrow(mns)) {
bg[m, ] <- (col[nc:1])[rank(mns[m, ])]
}
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.factor(as.numeric(levels(mSetObj$dataSet$cls))[mSetObj$dataSet$cls])
}
else {
cls <- mSetObj$dataSet$cls
}
cls.lbl <- levels(cls)
for (n in 1:ncol(mns)) {
graphics::points(x, y, bty = "n", pch = 22, bg = bg[, n],
cex = 3)
graphics::text(x[1], axis.lims[4], cls.lbl[n], srt = 45, adj = c(0.2,
0.5))
x <- x + shift/1.25
}
col <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(7, colorpalette))(50)
if (color.BW)
col <- rev(col)
nc <- length(col)
x <- rep(x[1] + shift, nc)
shifty <- (axis.lims[4] - axis.lims[3])/3
starty <- axis.lims[3] + shifty
endy <- axis.lims[3] + 2 * shifty
y <- seq(from = starty, to = endy, length = nc)
graphics::points(x, y, bty = "n", pch = 15, col = rev(col), cex = 2)
graphics::text(x[1], endy + shifty/8, "High")
graphics::text(x[1], starty - shifty/8, "Low")
graphics::par(op)
}
#' Plot PLS important features
#'
#' \code{met.plot_PLSImpScatter} visualizes the features' importance in the PLS-DA model by comparing VIP scores and coefficients.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (\code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param feat.nm (Character) Indicate the name of the feature. Choose \code{"coef.mean"} for average coefficients, or the name of a sample group in your dataset.
#' @param vip.nm (Character or a character vector) Choose \code{"Comp. 1"}, \code{"Comp. 2"}, etc., depending on the component for which the VIP scores should be shown. If more than one component are specified (the default is \code{c("Comp. 1", "Comp. 2")}), an average VIP score is calculated.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param vip.thresh (Numeric) Draw a vertical line in the plot indicating a chosen VIP relevance threshold (the default is \code{1}).
#' @param show.title (Logical) \code{TRUE} to add a title to the plot, or \code{FALSE} to not.
#' @param title (Character) Define the title if \code{show.title = TRUE}
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.ImpScatter_plot_feat.nm_.plot)}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
#' @import ggplot2
met.plot_PLSImpScatter <- function(mSetObj, imgName = "PLS_ImpScatter", format = "png", dpi = 300, width = NA,
feat.nm="coef.mean", vip.nm = c("Comp. 1", "Comp. 2"), plot = TRUE, export = FALSE, vip.thresh = 1,
show.title = FALSE, title){
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 8
}
else {
w <- width
}
h <- w * if_else(show.title==TRUE, 1.1, 1.0)
df <- data.frame(mSetObj$analSet$plsda$vip.mat, mSetObj$analSet$plsda$coef.mat[rownames(mSetObj$analSet$plsda$vip.mat),], check.names = TRUE)
feat <- make.names(feat.nm)
if(length(vip.nm)>1){
df$vip_avg <- rowMeans(df[,grep(paste0(make.names(vip.nm), collapse="|"), colnames(df))])
vip <- "vip_avg"
} else {
vip <- make.names(vip.nm)
}
p <-
ggplot(df, aes_string(x = vip, y = feat)) +
geom_point(color = "Black", fill="Gray", shape=21, size = 3.5, alpha = 0.6) +
ggrepel::geom_text_repel(
aes(label = rownames(df)),
box.padding = unit(0.25,
"lines"),
point.padding = unit(0.5, "lines")
) +
geom_vline(xintercept = vip.thresh,
linetype = "longdash", alpha = 0.4) +
ylab(if_else(feat.nm=="coef.mean", feat.nm, paste0("Coef. ", feat.nm))) +
theme_bw(base_size = 18) +
theme(axis.text=element_text(size = 20),
axis.title = element_text(size=24,face="bold"))
if(length(vip.nm)>1){
p <- p + xlab(paste0("average VIP (", paste(vip.nm, collapse = ", "), ")"))
} else {
p <- p + xlab(paste0("VIP - ", vip.nm))
}
if(show.title == TRUE){
p <- p + ggtitle(title) + theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
if (export == TRUE){
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
print(p)
grDevices::dev.off()
mSetObj[["imgSet"]][[paste0("pls.ImpScatter_", feat.nm)]] <- imgName
}
mSetObj[["imgSet"]][[paste0("pls.ImpScatter_plot_",feat.nm, ".plot")]] <- p
if (plot == TRUE){
print(p)
}
return(mSetObj)
}
#' Create heat map of hierarchically clustered features and samples
#'
#' \code{met.plot_heatmap} plots a heat map with optional filtering for top features based on results from t-tests/ANOVA, VIP or randomforest.
#'
#' @param mSetObj Input name of the created mSet object.
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (\code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param dataOpt (Character) Set data options. Choose \code{"norm"} to use data after normalization or any other string to use original values.
#' @param scaleOpt (Character) Set value standardization option. \code{"row"} for autoscaling of features, \code{"column"} for autoscaling of samples, \code{"none"} for no autoscaling.
#' @param smplDist (Character) Enter the sample distance method: \code{"euclidean"} for Euclidean distance, \code{"correlation"} for Pearson distance, or \code{"minkowski"} for Minkowski method.
#' @param clstDist (Character) Enter the feature clustering method: \code{"ward.d"}, \code{"average"}, \code{"complete"}, or \code{"single"}.
#' @param palette Input color palette choice: \code{"bwm"} for Blue-White-Red, \code{"gbr"} for Green-Black-Red, \code{"heat"} for Red-Yellow, \code{"topo"} for topological colors, \code{"gray"} for gray-scale.
#' @param viewOpt Set view options, choose \code{"detail"} or \code{"overview"}.
#' @param rowV (Logical) Shall clustering be applied to samples?
#' @param colV (Logical) Shall clustering be applied to features?
#' @param var.inx (Numeric vector) \emph{optional}: Provide indeces of selected features to be plotted.
#' @param border (Logical) Indicate whether or not to show cell-borders.
#' @param grp.ave (Logical) Shall only group averages be shown (as opposed to displaying every sample)?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added heat map. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$heatmap.avg.plot)} or \code{print(mSetObj$imgSet$heatmap.all.plot)}, depending on \code{grp.ave}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotHeatMap}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_heatmap <- function (mSetObj = NA, imgName = "Heatmap_features", format = "pdf", dpi = NULL,
width = NA, dataOpt = "norm", scaleOpt = "row",
smplDist = "euclidean", clstDist = "ward.D", palette = "bwm",
viewOpt = "detail", rowV = T, colV = T, var.inx = NULL,
border = T, grp.ave = F, plot = TRUE, export = TRUE)
{
imgName = paste(imgName, ".", format,
sep = "")
cls <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
cls.class <- mSetObj$dataSet$type.cls.lbl
mSetObj$analSet$htmap <- list(dist.par = smplDist, clust.par = clstDist)
if (dataOpt == "norm") {
my.data <- mSetObj$dataSet$norm
}
else {
if(!is.null(mSetObj$dataSet$prenorm)){
row.norm <- mSetObj$dataSet$prenorm
} else {
my.data <- qs::qread("prenorm.qs")
}
}
if (is.null(var.inx)) {
hc.dat <- as.matrix(my.data)
}
else {
hc.dat <- as.matrix(my.data[, var.inx])
}
colnames(hc.dat) <- substr(colnames(hc.dat), 1, 18)
if (cls.class == "integer") {
hc.cls <- as.factor(as.numeric(levels(cls))[cls])
}
else {
hc.cls <- cls
}
if (grp.ave) {
lvs <- levels(hc.cls)
my.mns <- matrix(ncol = ncol(hc.dat), nrow = length(lvs))
for (i in 1:length(lvs)) {
inx <- hc.cls == lvs[i]
my.mns[i, ] <- apply(hc.dat[inx, ], 2, mean)
}
rownames(my.mns) <- lvs
colnames(my.mns) <- colnames(hc.dat)
hc.dat <- my.mns
hc.cls <- as.factor(lvs)
}
if (palette == "gbr") {
colors <- grDevices::colorRampPalette(c("green", "black",
"red"), space = "rgb")(256)
}
else if (palette == "heat") {
colors <- grDevices::heat.colors(256)
}
else if (palette == "topo") {
colors <- grDevices::topo.colors(256)
}
else if (palette == "gray") {
colors <- grDevices::colorRampPalette(c("grey90", "grey10"),
space = "rgb")(256)
}
else {
colors <- rev(grDevices::colorRampPalette(RColorBrewer::brewer.pal(10,
"RdBu"))(256))
}
if (cls.type == "disc") {
annotation <- data.frame(class = hc.cls)
rownames(annotation) <- rownames(hc.dat)
}
else {
annotation <- NA
}
plot_dims <- get_pheatmap_dims(t(hc.dat), annotation, viewOpt,
width)
h <- plot_dims$height
w <- plot_dims$width
if (grp.ave) {
w <- nrow(hc.dat) * 25 + 300
w <- round(w/72, 2)
}
if (border) {
border.col <- "grey60"
}
else {
border.col <- NA
}
p <- function(){
if (cls.type == "disc") {
cols <- GetColorSchema(cls, palette == "gray")
uniq.cols <- unique(cols)
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.factor(as.numeric(levels(cls))[cls])
}
else {
cls <- cls
}
names(uniq.cols) <- unique(as.character(sort(cls)))
ann_colors <- list(class = uniq.cols)
pheatmap::pheatmap(t(hc.dat), annotation = annotation,
annotation_colors = ann_colors, fontsize = 8, fontsize_row = 8,
clustering_distance_rows = smplDist, clustering_distance_cols = smplDist,
clustering_method = clstDist, border_color = border.col,
cluster_rows = colV, cluster_cols = rowV, scale = scaleOpt,
color = colors)
}
else {
stats::heatmap(hc.dat, Rowv = rowTree, Colv = colTree, col = colors,
scale = "column")
}
}
if (export == TRUE) {
message(paste0("Exporting heat map of features to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$heatmap <- imgName
}
if (grp.ave == TRUE) {
p(); mSetObj$imgSet$heatmap.avg.plot <- grDevices::recordPlot()
} else {
p(); mSetObj$imgSet$heatmap.all.plot <- grDevices::recordPlot()
}
if (plot == TRUE) {
p()
}
return(mSetObj)
}
#' Volcano Plot
#'
#' \code{met.plot_volcano} generates a volcano plot for a specified contrast.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after ANOVA analysis (\code{\link[VisomX]{met.ANOVA.Anal}}) if \code{test = "anova"}.
#' @param grp1 (Character) Enter name of the first group for the contrast \code{grp1 vs. grp2}. If both group arguments are empty, the first two names in the list of groups are selected.
#' @param grp2 (Character) Enter name of the second group for the contrast \code{grp1 vs. grp2}. If both group arguments are empty, the first two names in the list of groups are selected.
#' @param test (Character) Choose a statistical tests. For \code{test = "ttest"}, \code{met.plot_volcano} runs \code{\link[VisomX]{met.Ttests.Anal}} with the chosen test parameters.
#' For \code{test = "anova"}, \code{\link[VisomX]{met.ANOVA.Anal}} must have been applied previously on the mSetObj.
#' @param paired (Logical) Is the data paired (\code{TRUE}) or not (\code{FALSE}). Only applicable for \code{test = "ttest"}.
#' @param nonpar (Logical) Use a non-parametric test (\code{TRUE}) or not (\code{FALSE}). Only applicable for \code{test = "ttest"}.
#' @param equal.var (Logical) The two groups have equal variance (\code{TRUE}) or not (\code{FALSE}). Only applicable for \code{test = "ttest"}.
#' @param log2fc.thresh (Numeric) Enter a relevance threshold for log2 fold changes, highlighted in the plot by vertical lines and colored compounds
#' @param threshp (Numeric) Enter a significance threshold for features based on T-test or ANOVA test results, highlighted in the plot by an horizontal line and colored compounds
#' @param pval.type (Character) Display and apply significance threshold to \code{"raw"} p values or adjusted p values (\code{"fdr"}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}). If \code{NULL}, the name \code{Plots/Volcano_grp1_vs_grp2} is assigned.
#' @param format (\code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param add_names (Logical) Display labels of significant features (\code{TRUE}) or not (\code{FALSE}).
#' @param label_size (Numeric) Font size for feature labels (if \code{add_names = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param silent (Logical) Shall the results of \code{\link[VisomX]{met.Ttests.Anal}} be printed in the console (\code{TRUE}) or not (\code{FALSE})?
#' @param test_condition (Logical) Add a subtitle with the applied data transformation and scaling methods be displayed below the plot title (\code{TRUE}) or not (\code{FALSE}).
#' @return The input mSet object with added volcano plot (generated by \code{\link[ggplot2]{ggplot}}). The plot can be retrieved from within R via \code{print(mSetObj$imgSet$volcano$grp1_vs_grp2.plot)}.
#' @references adapted from \code{\link[DEP]{plot_volcano}} (\url{https://bioconductor.org/packages/DEP/}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @import ggplot2
#' @export
met.plot_volcano <- function (mSetObj = NA, grp1, grp2, test = "ttest", paired = FALSE, nonpar = FALSE, equal.var = TRUE,
log2fc.thresh = 1, threshp = 0.05, pval.type = "fdr", imgName = NULL, format = "pdf",
add_names = TRUE, label_size = 3, dpi = NULL, width = NA, plot = TRUE, export = TRUE,
silent = FALSE, test_condition = FALSE)
{
if(!is.null(grp1) && !is.null(grp2)){
assertthat::assert_that(grp1 %in% levels(mSetObj$dataSet$cls),
grp2 %in% levels(mSetObj$dataSet$cls),
msg = paste0("'",grp1, "' or '", grp2, "' are not valid conditions in the dataset. Valid conditions are:\n",
paste(levels(mSetObj$dataSet$cls), collapse = ", ")))
}
if(is.null(grp1) && is.null(grp2)){
grp1 <- unique(levels(mSetObj$dataSet$cls))[1]
grp2 <- unique(levels(mSetObj$dataSet$cls))[2]
}
if(test=="ttest"){
mSetObj <- met.Ttests.Anal(mSetObj, grp1 = grp1, grp2 = grp2, nonpar, threshp, paired,
equal.var, pval.type, silent = silent)
p.value <- mSetObj$analSet$tt[[paste0(grp1,"_vs_",grp2)]][["p.value"]]
if (pval.type == "fdr") {
p.value <- p.adjust(p.value, "fdr")
}
} else if(test=="anova"&&!is.null(mSetObj$analSet$aov)){
p.value <- stats::setNames(mSetObj$analSet$aov$post.hoc.all[,'p.adj'], rownames(mSetObj$analSet$aov$post.hoc.all))
p.value <- sapply(1:length(p.value), function (x) str_split(p.value[x], "; "), simplify = TRUE)
p.value <-
stats::setNames(as.numeric(sapply(p.value, "[", match(
x = paste0(grp2, "-", grp1), table = unlist(str_split(
mSetObj$analSet$aov$post.hoc.all$contrasts[1], "; "
))
))),
rownames(mSetObj$analSet$aov$post.hoc.all))
} else if(test=="anova"&&is.null(mSetObj$analSet$aov)){
stop("Please perform 'met.ANOVA.Anal()' on your mSet before running 'met.plot_volcano' with 'test = \"anova\"'.")
}
inx.p <- p.value <= threshp
p.log <- -log10(p.value)
mSetObj <- met.FC.Anal(mSetObj, log2fc.thresh = log2fc.thresh, grp1 = grp1, grp2 = grp2)
max.xthresh <- log2fc.thresh
min.xthresh <- -log2fc.thresh
fc.log <- mSetObj$analSet$fc[[paste0(grp1,"_vs_",grp2)]][["fc.log"]]
fc.all <- mSetObj$analSet$fc[[paste0(grp1,"_vs_",grp2)]][["fc.all"]]
inx.up <- mSetObj$analSet$fc[[paste0(grp1,"_vs_",grp2)]][["inx.up"]]
inx.down <- mSetObj$analSet$fc[[paste0(grp1,"_vs_",grp2)]][["inx.down"]]
keep.inx <- names(inx.p) %in% names(inx.up)
inx.p <- inx.p[keep.inx]
p.log <- p.log[keep.inx]
inx.imp <- (inx.up | inx.down) & inx.p
sig.var <- cbind(fc.all[inx.imp, drop = F], fc.log[inx.imp,
drop = F], p.value[inx.imp, drop = F], p.log[inx.imp,
drop = F])
if (pval.type == "fdr") {
colnames(sig.var) <- c("FC", "log2(FC)",
"p.adjusted", "-log10(p)")
}
else {
colnames(sig.var) <- c("FC", "log2(FC)",
"raw.pval", "-log10(p)")
}
ord.inx <- order(sig.var[, 4], abs(sig.var[, 2]), decreasing = T)
sig.var <- sig.var[ord.inx, , drop = F]
sig.var <- signif(sig.var, 5)
fileName <- "volcano.csv"
fast.write.csv(signif(sig.var, 5), file = fileName)
volcano <- list(pval.type = pval.type, raw.threshx = log2fc.thresh,
raw.threshy = threshp, paired = paired, max.xthresh = max.xthresh,
min.xthresh = min.xthresh, thresh.y = -log10(threshp),
fc.all = fc.all, fc.log = fc.log, inx.up = inx.up, inx.down = inx.down,
p.log = p.log, inx.p = inx.p, sig.mat = sig.var)
if(!exists("volcano", mSetObj$analSet)){
mSetObj$analSet$volcano <- list()
}
mSetObj$analSet$volcano[[paste0(grp1, "_vs_", grp2)]] <- volcano
if(!is.character(imgName)){
imgName <- paste0("Plots/Volcano_", grp1, "_vs_", grp2)
}
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 8
}
else {
w <- width
}
h <- w * 5.5/6.5
if(!exists("volcano", mSetObj$imgSet)){
mSetObj$imgSet$volcano <- list()
}
vcn <- mSetObj$analSet$volcano[[paste0(grp1, "_vs_", grp2)]]
imp.inx <- (vcn$inx.up | vcn$inx.down) & vcn$inx.p
de <- data.frame(cbind(vcn$fc.log, vcn$p.log))
de$Status <- "Non-SIG"
de$Status[vcn$inx.p & vcn$inx.up] <- "UP"
de$Status[vcn$inx.p & vcn$inx.down] <- "DOWN"
de$Status <- as.factor(de$Status)
mycols <- levels(de$Status)
mycols[mycols == "UP"] <- "#e31f26"
mycols[mycols == "DOWN"] <- "#387fb9"
mycols[mycols == "Non-SIG"] <- "#525352"
de$delabel <- NA
de$delabel[imp.inx] <- rownames(de)[imp.inx]
if (pval.type == "fdr") {
de$shape <- ifelse(de[,2] > 5.9, "triangle", "circle")
} else {
de$shape <- "circle"
}
if (pval.type == "fdr") {
if (max(de[,2]) > 6.0) {
de[,2][de[,2] > 5.9] <- 5.9
}
}
title_text <- paste0(grp1, " vs. ", grp2, if_else(
test_condition == TRUE,
paste0(
"\n(",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", "none"),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", "none"),
")"
),
"")
)
p <- ggplot(data = de, aes(x = de[, 1], y = de[, 2],
col = Status, label = delabel)) +
geom_point(aes(shape = shape), size = 2.5, alpha = 0.6) +
scale_shape(guide = "none") + # Hide shapes from legend
scale_color_manual(labels = c("Increased", "Not significant", "Decreased"),
values = mycols) +
geom_vline(xintercept = c(-vcn$raw.threshx, vcn$raw.threshx),
linetype = "longdash", alpha = 0.4) +
geom_hline(yintercept = -log10(vcn$raw.threshy), linetype = "longdash",
alpha = 0.4) +
labs(title = title_text) +
labs(x = "log2(fold change)", y = "-log10(p value)") +
theme_bw(base_size = 20) +
theme(legend.position = "bottom",title = element_text(size = exp(-nchar(paste0(grp1, " vs. ", grp2))/40)*if_else(
test_condition == TRUE, 25, 30)),
axis.title = element_text(size = 22), legend.title = element_text(size = 22), plot.title = element_text(hjust = 0.5)) +
ggh4x::force_panelsizes(rows = unit(if_else(test_condition==TRUE, 0.57, 0.61)*w, "in"),
cols = unit(if_else(test_condition==TRUE, 0.67, 0.72)*w, "in")) +
coord_cartesian(xlim = c(-max(abs(de[, 1])) - 0.05*max(abs(de[, 1])),
max(abs(de[, 1])) + 0.05*max(abs(de[, 1]))),
expand = TRUE)
# scale_x_continuous(breaks = scales::pretty_breaks(n = max(abs(de[,1])) + 1))
if (add_names) {
p <- p + ggrepel::geom_text_repel(
aes(label = delabel),
size = label_size,
box.padding = unit(0.25,
"lines"),
point.padding = unit(0.1, "lines"),
segment.size = 0.5,
show.legend = FALSE,
max.overlaps = 16
)
}
if (pval.type == "fdr") {
p <- p + labs(y = expression(-log[10] ~ "(adj. p-value)"))
}
if (export == TRUE){
message(paste0("Exporting volcano plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
print(p)
grDevices::dev.off()
mSetObj$imgSet$volcano[[paste0(grp1, "_vs_", grp2)]] <- imgName
}
mSetObj$imgSet$volcano[[paste0(grp1, "_vs_", grp2, ".plot")]] <- p
if (plot == TRUE){
print(p)
}
return(mSetObj)
}
#' Create 3D PCA loading plot
#'
#' \code{met.plot_PCA3DLoading} visualizes which features are responsible for the patterns seen among the samples in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format format.
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the loading plot.
#' @param inx3 (Numeric) Indicate the number of the principal component for the z-axis of the loading plot.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as file with the chosen format?
#' @return The input mSet object with added 3D scatter plot. The plot can be retrieved from within R by executing \code{met.print_PCA3DLoading(mSetObj$imgSet$pca.loading3d.plot)}.
#' @export
met.plot_PCA3DLoading <- function (mSetObj = NA, imgName = "PCA3DLoading", format = "json", inx1 = 1,
inx2 = 2, inx3 = 3, export = F)
{
cls <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
cls.class <- mSetObj$dataSet$type.cls.lbl
pca <- mSetObj$analSet$pca
pca3d <- list()
pca3d$loading$axis <- paste("Loading PC", c(inx1, inx2,
inx3), sep = "")
coords <- data.frame(t(signif(pca$rotation[, 1:3], 5)))
dists <- apply(coords, 2, function(x) dist(rbind(x, c(0, 0, 0))))
pca3d$loading$cols <- GetRGBColorGradient(dists)
cols_vectors <- stringr::str_replace_all(pca3d$loading$cols, "rgba\\(", "") %>% stringr::str_replace_all(., "\\)", "")
cols_vectors <- sapply(1:length(cols_vectors), function (x) str_split(cols_vectors[x], ","))
pca3d$loading$cols_hex <-
sapply(1:length(cols_vectors), function (x)
grDevices::rgb(cols_vectors[[x]][1], cols_vectors[[x]][2], cols_vectors[[x]][3],
alpha = cols_vectors[[x]][4],
maxColorValue = 255
))
colnames(coords) <- NULL
pca3d$loading$xyz <- coords
pca3d$loading$name <- rownames(pca$rotation)
pca3d$loading$entrez <- rownames(pca$rotation)
if (cls.class == "integer") {
clss <- as.character(sort(as.factor(as.numeric(levels(cls))[cls])))
}
else {
clss <- as.character(cls)
}
if (all.numeric(clss)) {
clss <- paste("Group", clss)
}
pca3d$cls = clss
imgName = paste(imgName, ".", format, sep = "")
if(export==TRUE){
message(paste0("Exporting 3D PCA loading plot to:\n\"", getwd(), imgName, "\""))
json.mat <- jsonlite::toJSON(pca3d)
sink(imgName)
cat(json.mat)
sink()
mSetObj$imgSet$pca.loading3d <- imgName
}
mSetObj$imgSet$pca.loading3d.plot <- pca3d
return(mSetObj)
}
#' Plot a generated 3D PCA loading plot in RStudio
#'
#' \code{met.print_PCA3DLoading} opens a new graphics device to show the generated 3D PCA loading plot.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after \code{\link[VisomX]{met.plot_PCA3DLoading}}.
#' @return A 3D scatter plot. To close the graphics device, execute \code{rgl::close3d()}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.print_PCA3DLoading <- function (mSetObj = NA)
{
colors <- mSetObj$imgSet$pca.loading3d.plot$loading$cols_hex
rgl::plot3d(x= mSetObj$imgSet$pca.loading3d.plot$loading$xyz[1,],
y= mSetObj$imgSet$pca.loading3d.plot$loading$xyz[3,],
z= mSetObj$imgSet$pca.loading3d.plot$loading$xyz[2,],
col = colors,
xlab = mSetObj$imgSet$pca.loading3d.plot$loading$axis[1],
ylab = mSetObj$imgSet$pca.loading3d.plot$loading$axis[3],
zlab = mSetObj$imgSet$pca.loading3d.plot$loading$axis[2],
type="s" , size = 1.5)
rgl::text3d(x=mSetObj$imgSet$pca.loading3d.plot$loading$xyz[1,],
y=mSetObj$imgSet$pca.loading3d.plot$loading$xyz[3,],
z=mSetObj$imgSet$pca.loading3d.plot$loading$xyz[2,],
mSetObj$imgSet$pca.loading3d.plot$loading$name, cex=0.6, adj = c(1.3, 1.1))
}
#' Create 3D PCA Score plot
#'
#' \code{met.plot_PCA3DScore} visualizes clusters of samples based on their similarity in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format format.
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the score plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the score plot.
#' @param inx3 (Numeric) Indicate the number of the principal component for the z-axis of the score plot.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as file with the chosen format?
#' @return The input mSet object with added 3D scatter plot. The plot can be retrieved from within R by executing \code{met.print_PCA3DScore(mSetObj$imgSet$pca.score3d.plot)}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PCA3DScore <- function (mSetObj = NA, imgName = "PCA3DScore", format = "json", inx1 = 1,
inx2 = 2, inx3 = 3, export = F)
{
cls <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
cls.class <- mSetObj$dataSet$type.cls.lbl
pca <- mSetObj$analSet$pca
pca3d <- list()
pca3d$score$axis <- paste("PC", c(inx1, inx2, inx3),
" (", 100 * round(mSetObj$analSet$pca$variance[c(inx1,
inx2, inx3)], 3), "%)", sep = "")
coords <- data.frame(t(signif(pca$x[, c(inx1, inx2, inx3)],
5)))
colnames(coords) <- NULL
pca3d$score$xyz <- coords
pca3d$score$name <- rownames(mSetObj$dataSet$norm)
if (cls.class == "integer") {
cls <- as.character(sort(as.factor(as.numeric(levels(cls))[cls])))
}
else {
cls <- as.character(cls)
}
if (all.numeric(cls)) {
cls <- paste("Group", cls)
}
pca3d$score$facA <- cls
cols <- unique(GetColorSchema(as.factor(cls)))
pca3d$score$colors <- apply(col2rgb(cols), 2, function(x) {
paste("rgb(", paste(x, collapse = ","), ")", sep = "")
})
imgName = paste(imgName, ".", format, sep = "")
if(export==TRUE){
message(paste0("Exporting 3D PCA score plot to:\n\"", getwd(), imgName, "\""))
json.obj <- jsonlite::toJSON(pca3d)
sink(imgName)
cat(json.obj)
sink()
mSetObj$imgSet$pca.score3d <- imgName
}
mSetObj$imgSet$pca.score3d.plot <- pca3d
return(mSetObj)
}
#' Plot a generated 3D PCA scores plot in RStudio
#'
#' \code{met.print_PCA3DScore} opens a new graphics device to show the generated 3D PCA score plot.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after \code{\link[VisomX]{met.plot_PCA3DScore}}.
#' @return A 3D scatter plot. To close the graphics device, execute \code{rgl::close3d()}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.print_PCA3DScore <- function (mSetObj = NA)
{
colors <- GetColorSchema(as.factor(cls))
rgl::plot3d(x= mSetObj$imgSet$pca.score3d.plot$score$xyz[1,],
y= mSetObj$imgSet$pca.score3d.plot$score$xyz[3,],
z= mSetObj$imgSet$pca.score3d.plot$score$xyz[2,],
col = colors,
xlab = mSetObj$imgSet$pca.score3d.plot$score$axis[1],
ylab = mSetObj$imgSet$pca.score3d.plot$score$axis[3],
zlab = mSetObj$imgSet$pca.score3d.plot$score$axis[2],
type="s" , size = 1.5)
rgl::par3d(windowRect = c(0, 0, 512, 512))
rgl::legend3d(
"right",
legend = levels(factor(
mSetObj$imgSet$pca.score3d.plot$score$facA
)),
col = unique(colors),
pch = 16
)
}
#' Create 3D PLS-DA loading plot
#'
#' \code{met.plot_PLS3DLoading} visualizes which features are responsible for the patterns seen among the samples in partial least squares-discriminant analysis.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format format.
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the loading plot.
#' @param inx3 (Numeric) Indicate the number of the principal component for the z-axis of the loading plot.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as file with the chosen format?
#' @return The input mSet object with added 3D scatter plot. The plot can be retrieved from within R by executing \code{met.print_PLS3DLoading(mSetObj$imgSet$pls.loading3d.plot)}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS3DLoading <- function (mSetObj = NA, imgName = "PLS3DLoading", format = "json", inx1 = 1,
inx2 = 2, inx3 = 3, export=F)
{
pls = mSetObj$analSet$plsr
coords <- signif(as.matrix(cbind(pls$loadings[, inx1], pls$loadings[,
inx2], pls$loadings[, inx3])), 5)
pls3d <- list()
pls3d$loading$axis <- paste("Loading Comp.", c(inx1, inx2,
inx3), sep = "")
coords0 <- coords <- data.frame(t(signif(pls$loadings[, c(inx1,
inx2, inx3)], 5)))
colnames(coords) <- NULL
pls3d$loading$xyz <- coords
pls3d$loading$name <- rownames(pls$loadings)
pls3d$loading$entrez <- rownames(pls$loadings)
dists <- apply(coords0, 2, function(x) dist(rbind(x, c(0, 0, 0))))
cols <- GetRGBColorGradient(dists)
pls3d$loading$cols <- cols
cols_vectors <- stringr::str_replace_all(cols, "rgba\\(", "") %>% stringr::str_replace_all(., "\\)", "")
cols_vectors <- sapply(1:length(cols_vectors), function (x) str_split(cols_vectors[x], ","))
pls3d$loading$cols_hex <-
sapply(1:length(cols_vectors), function (x)
grDevices::rgb(cols_vectors[[x]][1], cols_vectors[[x]][2], cols_vectors[[x]][3],
alpha = cols_vectors[[x]][4],
maxColorValue = 255
))
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.character(sort(as.factor(as.numeric(levels(mSetObj$dataSet$cls))[mSetObj$dataSet$cls])))
}
else {
cls <- as.character(mSetObj$dataSet$cls)
}
if (all.numeric(cls)) {
cls <- paste("Group", cls)
}
pls3d$cls = cls
imgName = paste(imgName, ".", format, sep = "")
if(export==TRUE){
message(paste0("Exporting 3D PLS-DA loading plot to:\n\"", getwd(), imgName, "\""))
json.mat <- jsonlite::toJSON(pls3d)
sink(imgName)
cat(json.mat)
sink()
mSetObj$imgSet$pls.loading3d <- imgName
}
mSetObj$imgSet$pls.loading3d.plot <- pls3d
return(mSetObj)
}
#' Plot a generated 3D PLS-DA loading plot in RStudio
#'
#' \code{met.print_PLS3DLoading} opens a new graphics device to show the generated 3D PLS-DA loading plot.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after \code{\link[VisomX]{met.plot_PLS3DLoading}}.
#' @return A 3D scatter plot. To close the graphics device, execute \code{rgl::close3d()}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.print_PLS3DLoading <- function (mSetObj = NA)
{
colors <- mSetObj$imgSet$pls.loading3d.plot$loading$cols_hex
rgl::plot3d(x= mSetObj$imgSet$pls.loading3d.plot$loading$xyz[1,],
y= mSetObj$imgSet$pls.loading3d.plot$loading$xyz[3,],
z= mSetObj$imgSet$pls.loading3d.plot$loading$xyz[2,],
col = colors,
xlab = mSetObj$imgSet$pls.loading3d.plot$loading$axis[1],
ylab = mSetObj$imgSet$pls.loading3d.plot$loading$axis[3],
zlab = mSetObj$imgSet$pls.loading3d.plot$loading$axis[2],
type="s" , size = 1.5)
rgl::text3d(x=mSetObj$imgSet$pls.loading3d.plot$loading$xyz[1,],
y=mSetObj$imgSet$pls.loading3d.plot$loading$xyz[3,],
z=mSetObj$imgSet$pls.loading3d.plot$loading$xyz[2,],
mSetObj$imgSet$pls.loading3d.plot$loading$name, cex=0.6, adj = c(1.3, 1.1))
}
#' Create 3D PLS-DA Score plot
#'
#' \code{met.plot_PLS3DScore} visualizes clusters of samples based on their similarity in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format format.
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the score plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the score plot.
#' @param inx3 (Numeric) Indicate the number of the principal component for the z-axis of the score plot.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as file with the chosen format?
#' @return The input mSet object with added 3D scatter plot. The plot can be retrieved from within R by executing \code{met.print_PLS3DScore(mSetObj$imgSet$pls.score3d.plot)}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS3DScore <- function (mSetObj = NA, imgName = "PLS3DScore", format = "json", inx1 = 1,
inx2 = 2, inx3 = 3, export=F)
{
cls1 <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
cls.class <- mSetObj$dataSet$type.cls.lbl
pls3d <- list()
pls3d$score$axis <- paste("Component", c(inx1, inx2,
inx3), " (", round(100 * mSetObj$analSet$plsr$Xvar[c(inx1,
inx2, inx3)]/mSetObj$analSet$plsr$Xtotvar, 1), "%)",
sep = "")
coords <- data.frame(t(signif(mSetObj$analSet$plsr$score[,
c(inx1, inx2, inx3)], 5)))
colnames(coords) <- NULL
pls3d$score$xyz <- coords
pls3d$score$name <- rownames(mSetObj$dataSet$norm)
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.character(sort(as.factor(as.numeric(levels(cls1))[cls1])))
}
else {
cls <- as.character(cls1)
}
if (all.numeric(cls)) {
cls <- paste("Group", cls)
}
pls3d$score$facA <- cls
cols <- unique(GetColorSchema(cls1))
pls3d$score$colors <- apply(col2rgb(cols), 2, function(x) {
paste("rgb(", paste(x, collapse = ","), ")", sep = "")
})
imgName = paste(imgName, ".", format, sep = "")
if(export==TRUE){
message(paste0("Exporting 3D PLS-DA score plot to:\n\"", getwd(), imgName, "\""))
json.obj <- jsonlite::toJSON(pls3d)
sink(imgName)
cat(json.obj)
sink()
mSetObj$imgSet$pls.score3d <- imgName
}
mSetObj$imgSet$pls.score3d.plot <- pls3d
return(mSetObj)
}
#' Plot a generated 3D PLS-DA scores plot in RStudio
#'
#' \code{met.print_PLS3DScore} opens a new graphics device to show the generated 3D PLS-DA score plot.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after first after \code{\link[VisomX]{met.plot_PLS3DScore}}.
#' @return A 3D scatter plot. To close the graphics device, execute \code{rgl::close3d()}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.print_PLS3DScore <- function (mSetObj = NA)
{
colors <- GetColorSchema(as.factor(cls))
rgl::plot3d(x= mSetObj$imgSet$pls.score3d.plot$score$xyz[1,],
y= mSetObj$imgSet$pls.score3d.plot$score$xyz[3,],
z= mSetObj$imgSet$pls.score3d.plot$score$xyz[2,],
col = colors,
xlab = mSetObj$imgSet$pls.score3d.plot$score$axis[1],
ylab = mSetObj$imgSet$pls.score3d.plot$score$axis[3],
zlab = mSetObj$imgSet$pls.score3d.plot$score$axis[2],
type="s" , size = 1.5)
rgl::par3d(windowRect = c(0, 0, 512, 512))
rgl::legend3d(
"right",
legend = levels(factor(
mSetObj$imgSet$pls.score3d.plot$score$facA
)),
col = unique(colors),
pch = 16
)
}
#' Plot PLS-DA classification performance using different components, permutation
#'
#' Permutation test is a technique for testing a hypothesis of no effect, when the distribution of the test statistic is unknown. The objective of this test is to confirm that the initial model is superior to other models obtained by permuting the class labels and randomly assigning them to different individuals or, in other words, to answer the question "what is the model's performance if the groups are formed randomly".
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param title (Logical) \code{TRUE} to add a title with the used normalization, transformation, and scaling method, or \code{FALSE} to not add any title.
#' @return The input mSet object with added permutation test plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.permut.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLS.Permutation}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS.Permutation <- function (mSetObj = NA, imgName = "PLSDA-Permutation", format = "pdf", dpi = NULL,
width = NA, plot = TRUE, export = TRUE, title = FALSE)
{
bw.vec <- mSetObj$analSet$plsda$permut
len <- length(bw.vec)
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 8
}
else if (width == 0) {
w <- 7
}
else {
w <- width
}
h <- w * 6/8
p <- function(){
graphics::par(mar = c(5, 5, if_else(title==TRUE, 3, 2), 4))
hst <- graphics::hist(
bw.vec,
breaks = "FD",
freq = T,
ylab = "Frequency",
xlab = paste0(
"Permutation test statistics",
if_else(
mSetObj[["analSet"]][["plsda"]][["permut.type"]] == "separation distance",
" (separation distance)",
" (prediction accuracy)"
)
),
col = "#abd9e9",
main = ""
)
h <- max(hst$counts)
graphics::arrows(x0=bw.vec[1]+max(bw.vec)/20, y0=if_else(title==TRUE, h/1.45, h/5),
x1=bw.vec[1], y1=h/35, col ="#d73027", lwd=1.5, length=0.15, angle=25, xpd = T)
if(title==TRUE){
title.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::title(stringr::str_wrap(title.text, 55), line = 1, cex.main = 1.2)
}
graphics::text(x=bw.vec[1]+max(bw.vec)/20, y=if_else(title==TRUE, h/1.15, h/3.5), labels=paste("Observed \n statistic \n",
mSetObj$analSet$plsda$permut.p), xpd = T, cex=0.95)
}
if(export == TRUE){
message(paste0("Exporting plots with PLS-DA model permutation results to:\n\"", getwd(), "/", imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$pls.permut <- imgName
}
p(); mSetObj$imgSet$pls.permut.plot <- grDevices::recordPlot()
grDevices::dev.off()
if (plot == TRUE){
p()
}
return(mSetObj)
}
#' Plot PLS-DA classification performance using different components
#'
#' \code{met.plot_PLS.Crossvalidation} performs cross-validation (CV) on the generated PLS-DA model, where a fraction of data is held back, and the model trained on the rest. In each CV, the predicted data are compared with the original data, and the sum of squared errors is calculated. The prediction error is then summed over all samples (Predicted Residual Sum of Squares or PRESS). For convenience, the PRESS is divided by the initial sum of squares and subtracted from 1 to resemble the scale of the R2. Good predictions will have low PRESS or high Q2. Generally speaking, a model with an R2 (and Q2) value above 0.7 can be considered predictive. It is possible to have negative Q2, which means that your model is not at all predictive or is overfitted.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param title (Logical) \code{TRUE} to add a title with the used normalization, transformation, and scaling method, or \code{FALSE} to not add any title.
#' @return The input mSet object with added cross validation test plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.crossvalidation.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLS.Classification}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS.Crossvalidation <- function (mSetObj = NA, imgName = "PLSDA-CrossValidation", format = "pdf", dpi = NULL,
width = NA, plot = TRUE, export = TRUE, title = FALSE)
{
mSetObj <- mSetObj
res <- mSetObj$analSet$plsda$fit.info
colnames(res) <- 1:ncol(res)
best.num <- mSetObj$analSet$plsda$best.num
choice <- mSetObj$analSet$plsda$choice
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 7
}
else if (width == 0) {
w <- 7
}
else {
w <- width
}
h <- w * 5/6.7
p <- function(){
graphics::par(mar = c(5, 5, if_else(title==TRUE, 3, 2), 7))
graphics::barplot(res, beside = TRUE, col = c("lightblue", "mistyrose",
"lightcyan"), ylim = c(0, 1.05), xlab = "Number of components",
ylab = "Performance")
if (choice == "Q2") {
graphics::text((best.num - 1) * 3 + best.num + 2.5, res[3, best.num] +
0.02, labels = "*", cex = 2.5, col = "red")
}
else if (choice == "R2") {
graphics::text((best.num - 1) * 3 + best.num + 1.5, res[2, best.num] +
0.02, labels = "*", cex = 2.5, col = "red")
}
else {
graphics::text((best.num - 1) * 3 + best.num + 0.5, res[1, best.num] +
0.02, labels = "*", cex = 2.5, col = "red")
}
xpos <- ncol(res) * 3 + ncol(res) + 1
graphics::legend(xpos, 1, rownames(res), fill = c("lightblue",
"mistyrose", "lightcyan"), xpd = T)
if(title==TRUE){
title.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all commas that are not followed by [space]letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::title(stringr::str_wrap(title.text, 55), line = 1, cex.main = 1.2)
}
}
if(export == TRUE){
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$pls.crossvalidation <- imgName
}
p(); mSetObj$imgSet$pls.crossvalidation.plot <- grDevices::recordPlot()
grDevices::dev.off()
if(plot==TRUE){
p()
}
return(mSetObj)
}
NicWir/VisomX documentation built on Dec. 8, 2024, 1:27 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions met.plot_PLSImpScatter
Documented in met.plot_PLSImpScatter
#' Plot a heatmap of compounds with missing values
#'
#' \code{met.plot_missval} generates a heatmap of compounds
#' with missing values to discover whether values are missing at random or not.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container with missing values.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF and PNG file?
#' @param fontsize (Numeric) Font sized used in the plot.
#' @return The input mSet object with added heatmap indicating whether values are missing (0) or not (1). The plot can be retrieved from within R via \code{ComplexHeatmap::draw(mSetObj$imgSet$missval_heatmap.plot)}
#' (generated by \code{\link[ComplexHeatmap]{Heatmap}}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_missval <- function (mSetObj, plot = TRUE, imgName = "Plots/MissValPlot", format = "pdf", export = TRUE, fontsize = 12)
{
assay <- mSetObj[["dataSet"]][["data_orig"]]
if (!any(is.na(assay)) & !any(assay==0)) {
stop("No missing or zero values in '", deparse(substitute(se)),
"'", call. = FALSE)
}
imgName = paste(imgName,".", format,
sep = "")
# Make assay data binary (1 = valid value, 0 = missing value)
df <- assay %>% data.frame(check.names = FALSE)
df <- t(sapply(df, as.numeric))
colnames(df) <- rownames(assay)
missval <-
df[apply(df, 1, function(x)
any(is.na(x))) | apply(df, 1, function(x)
any(x==0)),]
missval <- ifelse(is.na(missval)|missval==0, 0, 1)
# Plot binary heatmap
ht2 = ComplexHeatmap::Heatmap(missval, col = c("white", "black"),
column_names_side = "top", show_row_names = FALSE,
show_column_names = TRUE, name = "Missing values pattern",
column_names_gp = grid::gpar(fontsize = fontsize), heatmap_legend_param = list(at = c(0,
1), labels = c("Missing value", "Valid value")))
if (export == TRUE){
message(paste0("Exporting Missing values pattern to:\n\"", getwd(), "/", imgName, "\""))
if (format == "pdf") {
grDevices::pdf(imgName)
ComplexHeatmap::draw(ht2, heatmap_legend_side = "top")
grDevices::dev.off()
}
if (format == "png") {
grDevices::png(imgName,
width = 6, height = 6, units = 'in', res = 300)
ComplexHeatmap::draw(ht2, heatmap_legend_side = "top")
grDevices::dev.off()
}
mSetObj$imgSet$missval_heatmap <- imgName
}
if (plot == TRUE){
ComplexHeatmap::draw(ht2, heatmap_legend_side = "top")
}
mSetObj$imgSet$missval_heatmap.plot <- ComplexHeatmap::draw(ht2, heatmap_legend_side = "top")
return(mSetObj)
}
#' Visualize intensities of compounds with missing values
#'
#' \code{met.plot_detect} generates density and CumSum plots
#' of compound intensities with and without missing values.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container with missing values.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF and PNG file?
#' @param basesize (Numeric) Font sized used in the plot.
#' @return Density and CumSum plots of intensities of proteins with and without missing values. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$missval_density.plot)}.
#' (generated by \code{\link[ggplot2]{ggplot}}).
#' @import ggplot2
#' @references adapted from package DEP (\url{https://bioconductor.org/packages/DEP/}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_detect <- function (mSetObj, plot = TRUE, imgName = "Plots/DensMissPlot", format = "pdf", export = TRUE, basesize = 10)
{
assay <- mSetObj[["dataSet"]][["data_orig"]] %>% replace(.==0, NA)
if (!any(is.na(assay)) & !any(assay==0)) {
stop("No missing or zero values in '", if_else(!is.null(mSetObj$dataSet$prenorm),
deparse(substitute(mSetObj[["dataSet"]][["prenorm"]])),
if_else(!is.null(mSetObj$dataSet$data_proc),
deparse(substitute(mSetObj[["dataSet"]][["data_proc"]])),deparse(substitute(mSetObj[["dataSet"]][["data_orig"]])))),
"'", call. = FALSE)
}
imgName = paste(imgName,".", format,
sep = "")
df <- assay %>% data.frame(check.names = FALSE)+1
df <- t(sapply(df, as.numeric)) %>% data.frame() %>% log10()
colnames(df) <- rownames(assay)
df <- df %>% tibble::rownames_to_column() %>% tidyr::gather(ID, val, -rowname)
stat <- df %>% group_by(rowname) %>% summarize(mean = mean(val,
na.rm = TRUE), missval = any(is.na(val)|val==0))
cumsum <- stat %>% group_by(missval) %>% arrange(mean) %>%
mutate(num = 1, cs = cumsum(num), cs_frac = cs/n())
p1 <- ggplot2::ggplot(stat, aes(mean, col = missval)) +
geom_density(na.rm = TRUE, show.legend = FALSE) +
stat_density(geom="line", position ="identity") +
labs(x = expression(log[10] ~ "Intensity"), y = "Density") +
guides(col = guide_legend(title = paste0(
"Missing values\n(",
round(
100 * nrow(cumsum[cumsum$missval == TRUE, ]) / nrow(cumsum[cumsum$missval == FALSE, ]),
digits = 1
),
"% missing)"
))) +
theme_minimal(base_size = basesize)
p2 <- ggplot2::ggplot(cumsum, aes(mean, cs_frac, col = missval)) +
geom_line() + labs(x = expression(log[10] ~ "Intensity"),
y = "Cumulative fraction") +
guides(col = guide_legend(title = paste0(
"Missing values\n(",
round(
100 * nrow(cumsum[cumsum$missval == TRUE, ]) / nrow(cumsum[cumsum$missval ==
FALSE, ]),
digits = 1
),
"% missing)"
))) +
theme_minimal(base_size = basesize)
df_per_group <- df
df_per_group$ID <- stringr::str_replace(df_per_group$ID, "_[:digit:]+$", "")
stat_per_group <- df_per_group %>% group_by(rowname, ID) %>% summarize(mean = mean(val,
na.rm = TRUE), missval = any(is.na(val)))
cumsum_per_group <- stat_per_group %>% group_by(missval) %>% arrange(mean) %>%
mutate(num = 1, cs = cumsum(num), cs_frac = cs/n())
p1_per_group <- ggplot2::ggplot(stat_per_group, aes(mean, col = missval)) +
geom_density(na.rm = TRUE, show.legend = FALSE) +
stat_density(geom = "line", position = "identity") +
labs(x = expression(log[2] ~ "Intensity"), y = "Density") +
guides(col = guide_legend(title = paste0(
"Missing values\n(",
round(
100 * nrow(cumsum_per_group[cumsum_per_group$missval == TRUE, ]) / nrow(cumsum_per_group[cumsum_per_group$missval ==
FALSE, ]),
digits = 1
),
"% missing)"
))) +
theme_minimal(base_size = basesize)
p2_per_group <-
ggplot2::ggplot(cumsum_per_group, aes(mean, cs_frac, col = missval)) +
geom_line() + labs(x = expression(log[2] ~ "Intensity"),
y = "Cumulative fraction") +
guides(col = guide_legend(title = paste0(
"Missing values\n(",
round(
100 * nrow(cumsum_per_group[cumsum_per_group$missval == TRUE, ]) / nrow(cumsum_per_group[cumsum_per_group$missval ==
FALSE, ]),
digits = 1
),
"% missing)"
))) +
theme_minimal(base_size = basesize)
if (export == TRUE){
message(paste0("Exporting density distributions and cumulative fraction of proteins with and without missing values to:\n\"", getwd(), "/", imgName, "\""))
if (format == "pdf") {
grDevices::pdf(imgName)
}
if (format == "png") {
grDevices::png(imgName,
width = 6,
height = 6,
units = 'in',
res = 300
)
}
gridExtra::grid.arrange(p1, p2, ncol = 1)
grDevices::dev.off()
mSetObj$imgSet$missval_density <- imgName
}
if (plot == TRUE){
gridExtra::grid.arrange(p1, p2, ncol = 1)
}
mSetObj$imgSet$missval_density.plot <- gridExtra::grid.arrange(p1, p2, ncol = 1)
return(mSetObj)
}
#' Visualize p values determined by ANOVA analysis.
#'
#' \code{met.plot_ANOVA} visualizes the significance of each compound in the dataset.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after ANOVA analysis (\code{\link[VisomX]{met.ANOVA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @return The input mSet object with added scatter plot with the p value for each compound analyzed. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$anova.plot)}.
#' (generated by \code{\link[ggplot2]{ggplot}}).
#' @import ggplot2
#' @references adapted from \code{\link[MetaboAnalystR]{PlotANOVA}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_ANOVA <- function (mSetObj = NA, imgName = "ANOVA_plot", format = "pdf", dpi = NULL,
width = NA, subtitle = FALSE, export = TRUE, plot=TRUE)
{
mSetObj <- mSetObj
lod <- -log10(mSetObj$analSet$aov$p.adj)
df <- data.frame(names=names(lod), values=lod, row.names = NULL)
imgName = paste(imgName,".", format,
sep = "")
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 7
}
else {
w <- width
}
h <- w * 6/9
mSetObj$imgSet$anova <- imgName
g <- ggplot2::ggplot(df, aes(x=as.numeric(row.names(df)), y=values)) +
geom_point(aes(fill=mSetObj$analSet$aov$inx.imp), color="black", shape=21, size=2, alpha=0.6) +
scale_fill_manual(values = c("#1170AA", "#EF6F6A"), labels=c("unsignificant", "significant")) +
geom_text(aes(label = names), size=2, hjust=1, vjust=-1) +
labs(x="Metabolites", y="-Log10(adjusted p-value)", title="One-way ANOVA") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.title = element_blank(), legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)) +
geom_hline(yintercept = -log10(mSetObj[["analSet"]][["aov"]][["raw.thresh"]]), linetype = "longdash",
alpha = 0.4)
if(subtitle==TRUE){
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are followed by [space]letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
g <- g + labs(subtitle=stringr::str_wrap(subtitle.text), 55)
}
if(export == TRUE){
message(paste0("Exporting anova plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
print(g)
grDevices::dev.off()
}
mSetObj$imgSet$anova.plot <- g
if(plot==TRUE){
print(g)
}
return(mSetObj)
}
#' Visualize effect of normalization on sample density distributions.
#'
#' \code{met.plot_SampleNormSummary} visualizes the density distributions of samples before and after normalization, transformation, and scaling.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after normalization (\code{\link[VisomX]{met.normalize}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param show_prenorm (Logical, \code{TRUE} or \code{FALSE}) Shall density distributions before normalization be displayed?
#' @return The input mSet object with added plot: the top is a density plot and the bottom is a box plot.
#' In a boxplot, the bottom and top of the box are always the 25th and 75th percentile (the lower and upper quartiles, or Q1 and Q3, respectively), and the band near the middle of the box is always the 50th percentile (the median or Q2). The upper whisker is located at the smaller of the maximum x value and Q3 + 1.5 x IQR (Interquantile Range), whereas the lower whisker is located at the larger of the smallest x value and Q1 - 1.5 x IQR.
#' The plot can be retrieved from within R via \code{print(mSetObj$imgSet$summary_norm_sample.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotSampleNormSummary}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_SampleNormSummary <- function (mSetObj = NA, imgName = "SampleNormSummary", format = "png", dpi = NULL,
width = NA, show_prenorm = TRUE, export = TRUE, plot=TRUE)
{
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- if_else(show_prenorm, 10, 5.5)
}
else if (width == 0) {
w <- if_else(show_prenorm, 10, 5.5)
}
else if (width > 0) {
w = width
}
h = w * if_else(show_prenorm, 1.25, 2.3)
if(!is.null(mSetObj$dataSet$prenorm)){
proc.data <- data.frame(mSetObj$dataSet$prenorm)
} else if(!is.null(mSetObj$dataSet$data_proc)){
proc.data <- mSetObj$dataSet$data_proc
} else {
proc.data <- qs::qread("data_proc.qs")
}
if(show_prenorm){
p <- function(){
graphics::layout(matrix(c(1, 2, 2, 2, 3, 4, 4, 4), 4, 2, byrow = FALSE))
pre.inx <- if (nrow(proc.data) < 50) {
1:nrow(proc.data)
} else {
sample(1:nrow(proc.data), 50)
}
namesVec <- rownames(proc.data[pre.inx, , drop = FALSE])
nm.inx <- namesVec %in% rownames(mSetObj$dataSet$norm)
namesVec <- namesVec[nm.inx]
pre.inx <- pre.inx[nm.inx]
norm.inx <- match(namesVec, rownames(mSetObj$dataSet$norm))
namesVec <- substr(namesVec, 1, 20)
rangex.pre <- range(proc.data[pre.inx, , drop = FALSE], na.rm = T)
rangex.norm <- range(mSetObj$dataSet$norm[norm.inx, , drop = FALSE],
na.rm = T)
x.label <- if_else(mSetObj$dataSet$type == "conc", "Concentration", "Intensity")
y.label <- "Samples"
op <- graphics::par(mar = c(6.5, 7, 0, 0), xaxt = "s")
plot(stats::density(apply(proc.data, 1, mean, na.rm = TRUE)), col = "darkblue",
las = 2, lwd = 2, main = "", xlab = "", ylab = "")
graphics::mtext("Density", 2, 5)
op <- graphics::par(mar = c(5.75, 8, 4, 0), xaxt = "s")
graphics::boxplot(t(proc.data[pre.inx, , drop = FALSE]), names = namesVec,
ylim = rangex.pre, las = 2, col = "lightgreen",
horizontal = T)
graphics::mtext(x.label, 1, 4.7, font = 2)
graphics::mtext(font = 2, "Before Normalization", 3, 1)
op <- graphics::par(mar = c(6.5, 7, 0, 2), xaxt = "s")
plot(stats::density(apply(mSetObj$dataSet$norm, 1, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
op <- graphics::par(mar = c(5.75, 8, 4, 2), xaxt = "s")
graphics::boxplot(t(mSetObj$dataSet$norm[norm.inx, , drop = FALSE]),
names = namesVec, ylim = rangex.norm, las = 2, col = "lightgreen",
ylab = "", horizontal = T)
subtitle.text <- paste0(
if_else(
mSetObj[["dataSet"]][["rownorm.method"]] == "N/A" &
mSetObj[["dataSet"]][["trans.method"]] == "N/A" &
mSetObj[["dataSet"]][["scale.method"]] == "N/A" ,
"(no normalization/transformation)",
stringr::str_replace_all(
paste0(
"(",
mSetObj[["dataSet"]][["rownorm.method"]],
", ",
mSetObj[["dataSet"]][["trans.method"]],
", ",
mSetObj[["dataSet"]][["scale.method"]],
")"
),
"N/A",
""
)
)
) %>% str_replace_all(., "(?<![:alpha:]), ", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., ",[:blank:](?![:alpha:])", "") # Remove all commas that are followed by [space]letter.
graphics::mtext(font = 2, paste0("After Normalization\n", stringr::str_wrap(subtitle.text, 55)), 3, 1)
graphics::mtext(font = 2, paste("Normalized", x.label), 1, 4.7)
}
} else {
p <- function(){
graphics::layout(matrix(c(1, 2, 2, 2), 4, 1, byrow = FALSE))
pre.inx <- if (nrow(proc.data) < 50) {
1:nrow(proc.data)
} else {
sample(1:nrow(proc.data), 50)
}
namesVec <- rownames(proc.data[pre.inx, , drop = FALSE])
nm.inx <- namesVec %in% rownames(mSetObj$dataSet$norm)
namesVec <- namesVec[nm.inx]
pre.inx <- pre.inx[nm.inx]
norm.inx <- match(namesVec, rownames(mSetObj$dataSet$norm))
namesVec <- substr(namesVec, 1, 20)
rangex.pre <- range(proc.data[pre.inx, , drop = FALSE], na.rm = T)
rangex.norm <- range(mSetObj$dataSet$norm[norm.inx, , drop = FALSE],
na.rm = T)
x.label <- if_else(mSetObj$dataSet$type == "conc", "Concentration", "Intensity")
y.label <- "Samples"
op <- graphics::par(mar = c(6.5, 7, 0, 0), xaxt = "s")
plot(density(apply(mSetObj$dataSet$norm, 1, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
graphics::mtext("Density", 2, 5)
op <- graphics::par(mar = c(5.75, 8, 4, 0), xaxt = "s")
graphics::boxplot(t(mSetObj$dataSet$norm[norm.inx, , drop = FALSE]),
names = namesVec, ylim = rangex.norm, las = 2, col = "lightgreen",
ylab = "", horizontal = T)
subtitle.text <- paste0(
if_else(
mSetObj[["dataSet"]][["rownorm.method"]] == "N/A" &
mSetObj[["dataSet"]][["trans.method"]] == "N/A" &
mSetObj[["dataSet"]][["scale.method"]] == "N/A" ,
"(no normalization/transformation)",
stringr::str_replace_all(
paste0(
"(",
mSetObj[["dataSet"]][["rownorm.method"]],
", ",
mSetObj[["dataSet"]][["trans.method"]],
", ",
mSetObj[["dataSet"]][["scale.method"]],
")"
),
"N/A",
""
)
)
) %>% str_replace_all(., "(?<![:alpha:]), ", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., ",[:blank:](?![:alpha:])", "") # Remove all commas that are followed by [space]letter.
graphics::mtext(font = 2, paste0("After Normalization\n", stringr::str_wrap(subtitle.text, 55)), 3, 1)
graphics::mtext(font = 2, paste("Normalized", x.label), 1, 4.7)
}
}
if(export == TRUE){
message(paste0("Exporting density distributions of samples before and after normalization to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$summary_norm_sample <- imgName
}
temp <- tempfile()
grDevices::pdf(file=temp)
p(); mSetObj$imgSet$summary_norm_sample.plot <- grDevices::recordPlot()
grDevices::dev.off()
file.remove(temp)
if(plot == TRUE){
p()
}
return(mSetObj)
}
#' Visualize effect of normalization on feature density distributions.
#'
#' \code{met.plot_FeatureNormSummary} visualizes the density distributions of features before and after normalization, transformation, and scaling.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after normalization (\code{\link[VisomX]{met.normalize}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param pre.inx (Numeric vector, or \code{NULL}) Provide index of compounds in your dataset that shall be diplayed with box plots. If \code{NULL}, the compounds with minimum and maximum values as well as up to 50 additional features are selected automatically.
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param show_prenorm (Logical, \code{TRUE} or \code{FALSE}) Shall density distributions before normalization be displayed?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @return The input mSet object with added plot: the top is a density plot and the bottom is a box plot.
#' In a boxplot, the bottom and top of the box are always the 25th and 75th percentile (the lower and upper quartiles, or Q1 and Q3, respectively), and the band near the middle of the box is always the 50th percentile (the median or Q2). The upper whisker is located at the smaller of the maximum x value and Q3 + 1.5 x IQR (Interquantile Range), whereas the lower whisker is located at the larger of the smallest x value and Q1 - 1.5 x IQR.
#' The plot can be retrieved from within R via \code{print(mSetObj$imgSet$summary_norm_feature.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotFeatureNormSummary}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_FeatureNormSummary <- function (mSetObj = NA, imgName = "FeatureNormSummary", format = "pdf", dpi = NULL, pre.inx = NULL,
width = NA, show_prenorm = TRUE, export = TRUE, plot = TRUE)
{
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- if_else(show_prenorm, 10.5, 5.75)
}
else if (width == 0) {
w = if_else(show_prenorm, 7.2, 3.6)
}
else if (width > 0) {
w = width
}
h = w * if_else(show_prenorm, 1.1905, 2)
if(!is.null(mSetObj$dataSet$prenorm)){
proc.data <- data.frame(mSetObj$dataSet$prenorm)
} else if(!is.null(mSetObj$dataSet$data_proc)){
proc.data <- mSetObj$dataSet$data_proc
} else {
proc.data <- qs::qread("data_proc.qs")
}
if(show_prenorm){
p <- function(){
graphics::layout(matrix(c(1, 2, 2, 2, 3, 4, 4, 4), 4, 2, byrow = FALSE))
if (is.null(pre.inx)) {
if (ncol(proc.data) > 50) {
min.inx <- which.min(colMeans(proc.data))
max.inx <- which.max(colMeans(proc.data))
order.inx <- order(colMeans(proc.data))
pre.inx <- unname(c(min.inx,
order.inx[!order.inx %in% c(min.inx, max.inx)][round(seq(2, (ncol(proc.data) - 2),
ncol(proc.data) / 50))],
max.inx))
} else {
pre.inx <- 1:nrow(proc.data)
}
}
namesVec <- colnames(proc.data[, pre.inx, drop = FALSE])
nm.inx <- namesVec %in% colnames(mSetObj$dataSet$norm)
namesVec <- namesVec[nm.inx]
pre.inx <- pre.inx[nm.inx]
norm.inx <- match(namesVec, colnames(mSetObj$dataSet$norm))
namesVec <- substr(namesVec, 1, 12)
rangex.pre <- range(proc.data[, pre.inx, drop = FALSE], na.rm = T)
rangex.norm <- range(mSetObj$dataSet$norm[, norm.inx, drop = FALSE],
na.rm = T)
x.label <- if_else(mSetObj$dataSet$type == "conc", "Concentration", "Intensity")
y.label <- if (mSetObj$dataSet$type == "conc") {
"Compounds"
} else if (mSetObj$dataSet$type == "specbin") {
"Spectra Bins"
} else if (mSetObj$dataSet$type == "nmrpeak") {
"Peaks (ppm)"
} else if (mSetObj$dataSet$type == "mspeak") {
"Peaks (mass)"
} else {
"Peaks(mz/rt)"
}
if (anal.type == "roc" & mSetObj$dataSet$roc_cols ==
1) {
op <- graphics::par(mar = c(4, 7, 4, 0), xaxt = "s")
graphics::plot.new()
}
else {
op <- graphics::par(mar = c(4, 7, 4, 0), xaxt = "s")
plot(density(apply(proc.data, 2, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
graphics::mtext("Density", 2, 5)
graphics::mtext(font = 2, "Before Normalization", 3, 1)
}
op <- graphics::par(mar = c(7, 7, 0.5, 0), xaxt = "s")
graphics::boxplot(proc.data[, pre.inx, drop = FALSE], names = namesVec,
ylim = rangex.pre, las = 2, col = "lightgreen",
horizontal = T, show.names = T)
graphics::mtext(x.label, 1, 5, font = 2)
if (anal.type == "roc" & mSetObj$dataSet$roc_cols ==
1) {
op <- graphics::par(mar = c(4, 7, 4, 2), xaxt = "s")
graphics::plot.new()
}
else {
op <- graphics::par(mar = c(4, 7, 4, 2), xaxt = "s")
plot(density(apply(mSetObj$dataSet$norm, 2, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
subtitle.text <- paste0(
if_else(
mSetObj[["dataSet"]][["rownorm.method"]] == "N/A" &
mSetObj[["dataSet"]][["trans.method"]] == "N/A" &
mSetObj[["dataSet"]][["scale.method"]] == "N/A" ,
"(no normalization/transformation)",
stringr::str_replace_all(
paste0(
"(",
mSetObj[["dataSet"]][["rownorm.method"]],
", ",
mSetObj[["dataSet"]][["trans.method"]],
", ",
mSetObj[["dataSet"]][["scale.method"]],
")"
),
"N/A",
""
)
)
) %>% str_replace_all(., "(?<![:alpha:]), ", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., ",[:blank:](?![:alpha:])", "") # Remove all commas that are followed by [space]letter.
graphics::mtext(paste0("After Normalization"), 3, 1, font = 2)
}
op <- graphics::par(mar = c(7, 7, 0.5, 2), xaxt = "s")
graphics::boxplot(mSetObj$dataSet$norm[, norm.inx, drop = FALSE], names = namesVec,
ylim = rangex.norm, las = 2, col = "lightgreen",
horizontal = T, show.names = T)
graphics::mtext(font = 2, paste("Normalized", x.label), 1, 5)
}
} else {
p <- function(){
graphics::layout(matrix(c(1, 2, 2, 2), 4, 1, byrow = FALSE))
if(is.null(pre.inx)){
if (ncol(proc.data) > 50) {
min.inx <- which.min(colMeans(proc.data))
max.inx <- which.max(colMeans(proc.data))
order.inx <- order(colMeans(proc.data))
pre.inx <- unname(c(min.inx,
order.inx[!order.inx %in% c(min.inx, max.inx)][round(seq(2, (ncol(proc.data) - 2),
ncol(proc.data) / 50))],
max.inx))
} else {
pre.inx <- 1:nrow(proc.data)
}
}
namesVec <- colnames(proc.data[, pre.inx, drop = FALSE])
nm.inx <- namesVec %in% colnames(mSetObj$dataSet$norm)
namesVec <- namesVec[nm.inx]
pre.inx <- pre.inx[nm.inx]
norm.inx <- match(namesVec, colnames(mSetObj$dataSet$norm))
namesVec <- substr(namesVec, 1, 12)
rangex.pre <- range(proc.data[, pre.inx, drop = FALSE], na.rm = T)
rangex.norm <- range(mSetObj$dataSet$norm[, norm.inx, drop = FALSE],
na.rm = T)
x.label <- if_else(mSetObj$dataSet$type == "conc", "Concentration", "Intensity")
y.label <- if (mSetObj$dataSet$type == "conc") {
"Compounds"
}
else if (mSetObj$dataSet$type == "specbin") {
"Spectra Bins"
}
else if (mSetObj$dataSet$type == "nmrpeak") {
"Peaks (ppm)"
}
else if (mSetObj$dataSet$type == "mspeak") {
"Peaks (mass)"
}
else {
"Peaks(mz/rt)"
}
if (anal.type == "roc" & mSetObj$dataSet$roc_cols ==
1) {
op <- graphics::par(mar = c(7, 7, 0, 0), xaxt = "s")
graphics::plot.new()
}
else {
op <- graphics::par(mar = c(7, 7, 4, 0), xaxt = "s")
plot(density(apply(mSetObj$dataSet$norm, 2, mean, na.rm = TRUE)),
col = "darkblue", las = 2, lwd = 2, main = "",
xlab = "", ylab = "")
graphics::mtext("Density", 2, 5)
}
op <- graphics::par(mar = c(7, 7, 0, 2), xaxt = "s")
graphics::boxplot(mSetObj$dataSet$norm[, norm.inx, drop = FALSE], names = namesVec,
ylim = rangex.norm, las = 2, col = "lightgreen",
horizontal = T, show.names = T)
subtitle.text <- paste0(
if_else(
mSetObj[["dataSet"]][["rownorm.method"]] == "N/A" &
mSetObj[["dataSet"]][["trans.method"]] == "N/A" &
mSetObj[["dataSet"]][["scale.method"]] == "N/A" ,
"(no normalization/transformation)",
stringr::str_replace_all(
paste0(
"(",
mSetObj[["dataSet"]][["rownorm.method"]],
", ",
mSetObj[["dataSet"]][["trans.method"]],
", ",
mSetObj[["dataSet"]][["scale.method"]],
")"
),
"N/A",
""
)
)
) %>% str_replace_all(., "(?<![:alpha:]), ", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., ",[:blank:](?![:alpha:])", "") # Remove all commas that are followed by [space]letter.
graphics::mtext(font = 2, paste0("After Normalization\n", stringr::str_wrap(subtitle.text, 55)), 3, 1)
graphics::mtext(font = 2, paste("Normalized", x.label), 1, 4.7)
}
}
if(export == TRUE){
message(paste0("Exporting density distributions of features before and after normalization to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$summary_norm_feature <- imgName
}
temp <- tempfile()
grDevices::pdf(file=temp)
p(); mSetObj$imgSet$summary_norm_feature.plot <- grDevices::recordPlot()
grDevices::dev.off()
file.remove(temp)
if(plot == TRUE){
p()
}
return(mSetObj)
}
#' Pattern hunter, plot correlation heatmap
#'
#' \code{met.plot_CorrHeatMap_Samples} visualizes the correlations between samples.
#'
#' @param mSetObj Input name of the created mSet object.
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param cor.method (Character) Indicate the correlation method, \code{'pearson'}, \code{'spearman'}, or \code{'kendall'}.
#' @param colors (Character) Indicate the colors for the heatmap, \code{"bwm"} for default, \code{"gbr"} for red/green, \code{"heat"} for heat colors, \code{"topo"} for topo colors, and \code{"gray"} for gray scale.
#' @param viewOpt (Character) Indicate \code{"overview"} to get an overview of the heatmap, and \code{"detail"} to get a detailed view of the heat map.
#' @param fix.col (Logical) Fix colors (\code{TRUE}) or not (\code{FALSE}).
#' @param no.clst (Logical) Indicate if the correlations should be clustered (\code{TRUE}) or not (\code{FALSE}).
#' @param corrCutoff (Numeric) Set correlation cut off.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @return The input mSet object with added heat map showing the correlations between samples. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$corr.heatmap_samples.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotCorrHeatMap}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_CorrHeatMap_Samples <- function (mSetObj = NA, imgName = "correlation_samples", format = "pdf", dpi = 72,
width = NA, cor.method = "pearson", colors="bwm", viewOpt="overview", fix.col = FALSE,
no.clst = FALSE, corrCutoff = 0, plot = TRUE, export = TRUE)
{
target = "row"
main <- xlab <- ylab <- NULL
data <- mSetObj$dataSet$norm
corrCutoff <- as.numeric(corrCutoff)
if (target == "row") {
data <- t(data)
}
if (ncol(data) > 1000) {
filter.val <- apply(data.matrix(data), 2, stats::IQR, na.rm = T)
rk <- rank(-filter.val, ties.method = "random")
data <- as.data.frame(data[, rk <= 1000])
print("Data is reduced to 1000 vars ..")
}
colnames(data) <- substr(colnames(data), 1, 18)
corr.mat <- cor(data, method = cor.method)
corr.mat[abs(corr.mat) < corrCutoff] <- 0
mSetObj$analSet$pwcor <- list()
mSetObj$analSet$pwcor$data <- data
mSetObj$analSet$pwcor$cor.method <- cor.method
mSetObj$analSet$pwcor$no.clst <- no.clst
if (colors == "gbr") {
colors <- grDevices::colorRampPalette(c("green", "black",
"red"), space = "rgb")(256)
}
else if (colors == "heat") {
colors <- grDevices::heat.colors(256)
}
else if (colors == "topo") {
colors <- grDevices::topo.colors(256)
}
else if (colors == "gray") {
colors <- grDevices::colorRampPalette(c("grey90", "grey10"))(256)
}
else {
colors <- rev(grDevices::colorRampPalette(c(RColorBrewer::brewer.pal(10,
"RdBu"), "#001833"))(256))
}
imgName = paste(imgName, ".", format,
sep = "")
if (viewOpt == "overview") {
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 7.2
}
else {
w <- 7.2
}
h <- w
}
else {
if (ncol(corr.mat) > 50) {
myH <- ncol(corr.mat) * 12 + 40
}
else if (ncol(corr.mat) > 20) {
myH <- ncol(corr.mat) * 12 + 60
}
else {
myH <- ncol(corr.mat) * 12 + 120
}
h <- round(myH/72, 2)
if (is.na(width)) {
w <- h
}
else if (width == 0) {
w <- h <- 7.2
}
else {
w <- h <- 7.2
}
}
if (no.clst) {
rowv = FALSE
colv = FALSE
dendro = "none"
}
else {
rowv = TRUE
colv = TRUE
dendro = "both"
}
if (fix.col) {
breaks <- seq(from = -1,
to = 1,
length = 257)
res <-
pheatmap::pheatmap(
corr.mat,
fontsize = 11,
fontsize_row = 11,
cluster_rows = colv,
cluster_cols = rowv,
color = colors,
breaks = breaks,
silent=T
)
}
else {
res <-
pheatmap::pheatmap(
corr.mat,
fontsize = 11,
fontsize_row = 11,
cluster_rows = colv,
cluster_cols = rowv,
color = colors,
silent=T
)
}
if (export == TRUE) {
message(paste0("Exporting correlation heat map of samples to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
grid::grid.newpage()
grid::grid.draw(res$gtable)
grDevices::dev.off()
mSetObj$imgSet$corr.heatmap_samples <- imgName
}
if (!no.clst) {
new.ord <- res$tree_row$order
corr.mat <- corr.mat[new.ord, new.ord]
mSetObj$analSet$pwcor$new.ord <- new.ord
}
if(plot == TRUE){
grid::grid.newpage()
grid::grid.draw(res$gtable)
}
mSetObj$imgSet$corr.heatmap_samples.plot <- res
message(paste0("Exporting CSV table file with sample correlations to:\n\"", getwd(), "/correlation_table_samples.csv\""))
fast.write.csv(signif(corr.mat, 5), file = "correlation_table_samples.csv")
return(mSetObj)
}
#' Pattern hunter, plot correlation heatmap
#'
#' \code{met.plot_CorrHeatMap_Features} visualizes the correlations between features.
#'
#' @param mSetObj Input name of the created mSet object.
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param cor.method Indicate the correlation method, 'pearson', 'spearman', or 'kendall'.
#' @param colors Indicate the colors for the heatmap, "bwm" for default, "gbr" for red/green, "heat" for heat colors, "topo" for topo colors, and "gray" for gray scale.
#' @param viewOpt Indicate "overview" to get an overview of the heatmap, and "detail" to get a detailed view of the heatmap.
#' @param fix.col (Logical) Fix colors (\code{TRUE}) or not (\code{FALSE}).
#' @param no.clst (Logical) Indicate if the correlations should be clustered (\code{TRUE}) or not (\code{FALSE}).
#' @param corrCutoff set correlation cut off
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @return The input mSet object with added heat map showing the correlations between features. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$corr.heatmap_features.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotCorrHeatMap}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
#' @export
met.plot_CorrHeatMap_Features <- function (mSetObj = NA, imgName = "correlation_features", format = "pdf", dpi = NULL,
width = NA, cor.method = "pearson", colors="bwm", viewOpt="overview", fix.col = TRUE,
no.clst = FALSE, corrCutoff = 0, plot = TRUE, export = TRUE)
{
target = "col"
main <- xlab <- ylab <- NULL
data <- mSetObj$dataSet$norm
corrCutoff <- as.numeric(corrCutoff)
if (target == "row") {
data <- t(data)
}
if (ncol(data) > 1000) {
filter.val <- apply(data.matrix(data), 2, stats::IQR, na.rm = T)
rk <- rank(-filter.val, ties.method = "random")
data <- as.data.frame(data[, rk <= 1000])
print("Data is reduced to 1000 vars ..")
}
colnames(data) <- substr(colnames(data), 1, 18)
corr.mat <- cor(data, method = cor.method)
corr.mat[abs(corr.mat) < corrCutoff] <- 0
mSetObj$analSet$pwcor <- list()
mSetObj$analSet$pwcor$data <- data
mSetObj$analSet$pwcor$cor.method <- cor.method
mSetObj$analSet$pwcor$no.clst <- no.clst
if (colors == "gbr") {
colors <- grDevices::colorRampPalette(c("green", "black",
"red"), space = "rgb")(256)
}
else if (colors == "heat") {
colors <- grDevices::heat.colors(256)
}
else if (colors == "topo") {
colors <- grDevices::topo.colors(256)
}
else if (colors == "gray") {
colors <- grDevices::colorRampPalette(c("grey90", "grey10"))(256)
}
else {
colors <- rev(grDevices::colorRampPalette(c(RColorBrewer::brewer.pal(10,
"RdBu"), "#001833"))(256))
}
imgName = paste(imgName, ".", format,
sep = "")
if (viewOpt == "overview") {
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 7.2
}
else {
w <- 7.2
}
h <- w
}
else {
if (ncol(corr.mat) > 50) {
myH <- ncol(corr.mat) * 12 + 40
}
else if (ncol(corr.mat) > 20) {
myH <- ncol(corr.mat) * 12 + 60
}
else {
myH <- ncol(corr.mat) * 12 + 120
}
h <- round(myH/72, 2)
if (is.na(width)) {
w <- h
}
else if (width == 0) {
w <- h <- 7.2
}
else {
w <- h <- 7.2
}
}
if (no.clst) {
rowv = FALSE
colv = FALSE
dendro = "none"
}
else {
rowv = TRUE
colv = TRUE
dendro = "both"
}
if (fix.col) {
breaks <- seq(from = -1,
to = 1,
length = 257)
res <-
pheatmap::pheatmap(
corr.mat,
fontsize = 4 ^ (2 - (ncol(mSetObj$dataSet$norm) / 110)),
fontsize_row = 4 ^ (2 - (ncol(mSetObj$dataSet$norm) / 110)),
cluster_rows = colv,
cluster_cols = rowv,
color = colors,
border_color = NA,
breaks = breaks,
silent=T
)
}
else {
res <-
pheatmap::pheatmap(
corr.mat,
fontsize = 4 ^ (2 - (ncol(mSetObj$dataSet$norm) / 110)),
fontsize_row = 4 ^ (2 - (ncol(mSetObj$dataSet$norm) / 110)),
cluster_rows = colv,
cluster_cols = rowv,
color = colors,
border_color = NA,
silent=T
)
}
if (export == TRUE) {
message(paste0("Exporting correlation heat map of features to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
grid::grid.newpage()
grid::grid.draw(res$gtable)
grDevices::dev.off()
mSetObj$imgSet$corr.heatmap_features <- res
}
if (!no.clst) {
new.ord <- res$tree_row$order
corr.mat <- corr.mat[new.ord, new.ord]
mSetObj$analSet$pwcor$new.ord <- new.ord
}
if(plot == TRUE){
grid::grid.newpage()
grid::grid.draw(res$gtable)
}
mSetObj$imgSet$corr.heatmap_features.plot <- res
message(paste0("Exporting CSV table file with feature correlations to:\n\"", getwd(), "/correlation_table_features.csv\""))
fast.write.csv(signif(corr.mat, 5), file = "correlation_table_features.csv")
return(mSetObj)
}
#' Plot PCA scree plot
#'
#' \code{met.plot_PCAScree} visualizes the proportion of variance explained for each principal component.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param scree.num Numeric, input a number to indicate the number of principal components to display in the scree plot.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added line plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pca.scree.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPCAScree}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
#' @export
met.plot_PCAScree <- function (mSetObj = NA, imgName = "PCA_ScreePlot", format = "pdf", dpi = NULL,
width = NA, scree.num, plot = TRUE, export = FALSE)
{
stds <- mSetObj$analSet$pca$std[1:scree.num]
pcvars <- mSetObj$analSet$pca$variance[1:scree.num]
cumvars <- mSetObj$analSet$pca$cum.var[1:scree.num]
ylims <- range(c(pcvars, cumvars))
extd <- (ylims[2] - ylims[1])/10
miny <- ifelse(ylims[1] - extd > 0, ylims[1] - extd, 0)
maxy <- ifelse(ylims[2] + extd > 1, 1, ylims[2] + extd)
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 10
}
else if (width == 0) {
w <- 8
}
else {
w <- width
}
h <- w * 2/3
p <- function(){
graphics::par(mar = c(5, 5, 6, 3))
plot(pcvars, type = "l", col = "blue", main = "Scree plot",
xlab = "PC index", ylab = "Variance explained",
ylim = c(miny, maxy), axes = F)
graphics::text(pcvars, labels = paste(100 * round(pcvars, 3), "%"),
adj = c(-0.3, -0.5), srt = 45, xpd = T)
graphics::points(pcvars, col = "red")
graphics::lines(cumvars, type = "l", col = "green")
graphics::text(cumvars, labels = paste(100 * round(cumvars, 3), "%"),
adj = c(-0.3, -0.5), srt = 45, xpd = T)
graphics::points(cumvars, col = "red")
graphics::abline(v = 1:scree.num, lty = 3)
graphics::axis(2)
graphics::axis(1, 1:length(pcvars), 1:length(pcvars))
}
if (export == TRUE) {
message(paste0("Exporting PCA scree plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$pca.scree <- imgName
}
if (plot == TRUE){
p()
}
p(); mSetObj$imgSet$pca.scree.plot <- grDevices::recordPlot()
return(mSetObj)
}
#' Create 2D PCA score plot
#'
#' \code{met.plot_PCA2DScore} visualizes clusters of samples based on their similarity in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param pcx Specify the principal component on the x-axis
#' @param pcy Specify the principal component on the y-axis
#' @param reg Numeric, input a number between 0 and 1, 0.95 will display the 95 percent confidence regions, and 0 will not.
#' @param show Display sample names, \code{1} = show names, \code{0} = do not show names.
#' @param grey.scale Use grey-scale colors, \code{1} = grey-scale, \code{0} = not grey-scale.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pca.score2d_PCx_PCy.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPCA2DScore}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PCA2DScore <- function (mSetObj = NA, imgName = "PCA_2DScores", format = "pdf", dpi = NULL, subtitle = FALSE,
width = NA, pcx, pcy, reg = 0.95, show = 1, grey.scale = 0, plot = TRUE, export = TRUE)
{
cls <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
xlabel = paste("PC", pcx, "(", round(100 * mSetObj$analSet$pca$variance[pcx],
1), "%)")
ylabel = paste("PC", pcy, "(", round(100 * mSetObj$analSet$pca$variance[pcy],
1), "%)")
pc1 = mSetObj$analSet$pca$x[, pcx]
pc2 = mSetObj$analSet$pca$x[, pcy]
text.lbls <- substr(names(pc1), 1, 14)
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 7.2
}
else if (width == 0) {
w <- 7.2
}
else {
w <- width
}
if (exists("group_factor")) {
legend.nm <- unique(as.character(sort(group_factor)))
}
else{
legend.nm <- unique(as.character(sort(cls)))
}
h <- w-1
p <- function(){
graphics::plot.new()
op <- graphics::par(mar=c(5,5,if_else(subtitle==TRUE, 7.9, 3),3))
l <- graphics::legend(0, 0, bty='n', legend = legend.nm,
plot=FALSE, pch=c(1, 2), lty=c(1, 2))
# calculate right margin width in ndc
w <- graphics::grconvertX(l$rect$w, to='ndc') - graphics::grconvertX(0, to='ndc')
graphics::par(omd=c(0, 1-w, 0, 1))
if (cls.type == "disc") {
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.factor(as.numeric(levels(cls))[cls])
}
else {
cls <- cls
}
lvs <- levels(cls)
pts.array <- array(0, dim = c(100, 2, length(lvs)))
for (i in 1:length(lvs)) {
inx <- cls == lvs[i]
groupVar <- stats::var(cbind(pc1[inx], pc2[inx]), na.rm = T)
groupMean <- cbind(mean(pc1[inx], na.rm = T), mean(pc2[inx],
na.rm = T))
pts.array[, , i] <- ellipse::ellipse(groupVar, centre = groupMean,
level = reg, npoints = 100)
}
xrg <- range(pc1, pts.array[, 1, ])
yrg <- range(pc2, pts.array[, 2, ])
x.ext <- (xrg[2] - xrg[1])/12
y.ext <- (yrg[2] - yrg[1])/12
xlims <- c(xrg[1] - x.ext, xrg[2] + x.ext)
ylims <- c(yrg[1] - y.ext, yrg[2] + y.ext)
cols <- GetColorSchema(cls, grayscale = F)
uniq.cols <- unique(cols)
plot(pc1, pc2, xlab = xlabel, xlim = xlims, ylim = ylims,
ylab = ylabel, type = "n",
col = cols, pch = as.numeric(cls) + 1,
cex.lab = 2,
cex.axis = 2,
cex.main = 2,
cex.sub = 2)
graphics::title("PCA Scores Plot", line = if_else(subtitle==TRUE, 4.5, 1), cex.main = 2.5)
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::mtext(line = 3.9, padj = 1, cex=1.7, if_else(subtitle==TRUE, stringr::str_wrap(subtitle.text, 55), NULL))
graphics::grid(col = "lightgray", lty = "dotted", lwd = 1)
if (length(uniq.cols) > 1) {
names(uniq.cols) <- unique(as.character(sort(cls)))
}
for (i in 1:length(lvs)) {
if (length(uniq.cols) > 1) {
graphics::polygon(pts.array[, , i], col = grDevices::adjustcolor(uniq.cols[lvs[i]],
alpha = 0.2), border = NA)
}
else {
graphics::polygon(pts.array[, , i], col = grDevices::adjustcolor(uniq.cols,
alpha = 0.2), border = NA)
}
if (grey.scale) {
graphics::lines(pts.array[, , i], col = grDevices::adjustcolor("black",
alpha = 0.5), lty = 2)
}
}
pchs <- GetShapeSchema(mSetObj, show, grey.scale)
if (grey.scale) {
cols <- rep("black", length(cols))
}
if (show == 1) {
graphics::text(pc1, pc2, label = text.lbls, pos = 4, xpd = T,
cex = 0.75)
graphics::points(pc1, pc2, pch = pchs, col = cols)
}
else {
if (length(uniq.cols) == 1) {
graphics::points(pc1, pc2, pch = pchs, col = cols, cex = 1)
}
else {
if (grey.scale == 1 | (exists("shapeVec") &&
all(shapeVec >= 0))) {
my.cols <- grDevices::adjustcolor(cols, alpha.f = 0.4)
my.cols[pchs == 21] <- "black"
graphics::points(pc1, pc2, pch = pchs, col = my.cols,
bg = grDevices::adjustcolor(cols, alpha.f = 0.4), cex = 1.8)
}
else {
graphics::points(pc1, pc2, pch = 21, bg = grDevices::adjustcolor(cols,
alpha.f = 0.4), cex = 2)
}
}
}
uniq.pchs <- unique(pchs)
if (length(uniq.cols) != length(levels(cls))) {
if (mSetObj$dataSet$type.cls.lbl == "integer") {
names(cols) <- as.numeric(cls)
}
else {
names(cols) <- as.character(cls)
}
match.inx <- match(levels(cls), names(cols))
uniq.cols <- cols[match.inx]
}
if (grey.scale) {
uniq.cols <- "black"
}
graphics::legend(graphics::par('usr')[2], graphics::par('usr')[4], bty='n', legend = legend.nm, pch = uniq.pchs,
col = uniq.cols, cex = 1.3, xpd = NA)
}
else {
plot(pc1, pc2, xlab = xlabel, ylab = ylabel, type = "n",
main = "Scores Plot")
graphics::points(pc1, pc2, pch = 15, col = "magenta")
graphics::text(pc1, pc2, label = text.lbls, pos = 4, col = "blue",
xpd = T, cex = 0.8)
}
graphics::par(op)
}
if (export == TRUE){
message(paste0("Exporting 2D PCA scores plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[['imgSet']][[paste0('pca.score2d_PC', pcx, "_PC", pcy)]] <- imgName
}
p(); mSetObj[['imgSet']][[paste0('pca.score2d_PC', pcx, "_PC", pcy, ".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
if(plot==TRUE){
p()
}
return(mSetObj)
}
#' Plot PCA loadings and also set up the matrix for display
#'
#' \code{met.plot_PCA2DLoading} visualizes which features are responsible for the patterns seen among the samples in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param inx1 Numeric, indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 Numeric, indicate the number of the principal component for the y-axis of the loading plot.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pca.loading_PCx_PCy.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPCALoading}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PCA2DLoading <- function (mSetObj = NA, imgName = "PCA_2DLoadings", format = "pdf", dpi = NULL, subtitle = FALSE,
width = NA, inx1, inx2, export = TRUE, plot=TRUE)
{
loadings <- as.matrix(cbind(mSetObj$analSet$pca$rotation[,
inx1], mSetObj$analSet$pca$rotation[, inx2]))
ord.inx <- order(-abs(loadings[, 1]), -abs(loadings[, 2]))
loadings <- signif(loadings[ord.inx, ], 5)
ldName1 = paste("Loadings PC", inx1, "(", round(100 * mSetObj$analSet$pca$variance[inx1],
1), "%)")
ldName2 = paste("Loadings PC", inx2, "(", round(100 * mSetObj$analSet$pca$variance[inx2],
1), "%)")
colnames(loadings) <- c(ldName1, ldName2)
mSetObj$analSet$pca$imp.loads <- loadings
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 7.2
}
else if (width == 0) {
w <- 7.2
}
else {
w <- width
}
h <- w-0.5
plotType <- mSetObj$analSet$pca$loading.type
p <- function(){
graphics::plot.new()
graphics::par(mar = c(6, 6.2, if_else(subtitle==TRUE, 7.9, 3), 6), mgp=c(4.5, 1, 0))
plot(loadings[, 1], loadings[, 2], las = 1, xlab = ldName1,
ylab = ldName2,
cex.lab = 2,
cex.axis = 2,
cex.main = 2,
cex.sub = 2)
graphics::title("PCA Loading Scatterplot", line = if_else(subtitle==TRUE, 4.5, 1), cex.main = 2.5)
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::mtext(line = 3.9, padj = 1, cex=1.7, if_else(subtitle==TRUE, stringr::str_wrap(subtitle.text, 55), NULL))
mSetObj$pca.axis.lims <- graphics::par("usr")
graphics::grid(col = "lightgray", lty = "dotted", lwd = 1)
graphics::points(loadings[, 1], loadings[, 2], pch = 19, col = grDevices::adjustcolor("#339933",
alpha.f = 0.4))
if (plotType == "all") {
graphics::text(loadings[, 1], loadings[, 2], labels = substr(rownames(loadings),
1, 16), pos = 4, col = "gray30", xpd = T)
}
else if (plotType == "gray30") {
if (length(mSetObj$custom.cmpds) > 0) {
hit.inx <- rownames(loadings) %in% mSetObj$custom.cmpds
graphics::text(loadings[hit.inx, 1], loadings[hit.inx, 2],
labels = rownames(loadings)[hit.inx], pos = 4,
col = "blue", xpd = T)
}
}
}
if (export == TRUE){
message(paste0("Exporting 2D PCA loadings plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[['imgSet']][[paste0('pca.loading_PC', inx1, "_PC", inx2)]] <- imgName
}
p(); mSetObj[['imgSet']][[paste0('pca.loading_PC', inx1, "_PC", inx2, ".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
if (plot == TRUE){
p()
}
return(mSetObj)
}
#' Create 2D PLS-DA score plot
#'
#' \code{met.plot_PLS2DScore} visualizes clusters of samples based on their similarity in partial least squares-discriminant analysis.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the loading plot.
#' @param reg (Numeric) Enter a number between 0 and 1, 0.95 will display the 95 percent confidence regions, and 0 will not.
#' @param show Display sample names, \code{1} = show names, \code{0} = do not show names.
#' @param grey.scale Use grey-scale colors, \code{1} = grey-scale, \code{0} = not grey-scale.
#' @param use.sparse (Logical) Use a sparse algorithm (\code{TRUE}) or not (\code{FALSE}).
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.score2d_PCx_PCy.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLS2DScore}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS2DScore <- function (mSetObj = NA, imgName = "PLSDA_2DScore", format = "pdf", dpi = NULL,
width = NA, inx1=1, inx2=2, reg = 0.95, show = 1, grey.scale = 0, subtitle=FALSE,
use.sparse = FALSE, plot = TRUE, export = TRUE)
{
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 7.2
}
else if (width == 0) {
w <- 7.2
}
else {
w <- width
}
h <- w-1
cls1 <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
lv1 <- mSetObj$analSet$plsr$scores[, inx1]
lv2 <- mSetObj$analSet$plsr$scores[, inx2]
xlabel <- paste("Component", inx1, "(", round(100 *
mSetObj$analSet$plsr$Xvar[inx1]/mSetObj$analSet$plsr$Xtotvar,
1), "%)")
ylabel <- paste("Component", inx2, "(", round(100 *
mSetObj$analSet$plsr$Xvar[inx2]/mSetObj$analSet$plsr$Xtotvar,
1), "%)")
if (exists("group_factor")) {
legend.nm <- unique(as.character(sort(group_factor)))
}
else{
legend.nm <- unique(as.character(sort(cls1)))
}
p<- function(){
graphics::plot.new()
l <- graphics::legend(0, 0, bty='n', legend = legend.nm,
plot=FALSE, pch=c(1, 2), lty=c(1, 2))
# calculate right margin width in ndc
w <- graphics::grconvertX(l$rect$w, to='ndc') - graphics::grconvertX(0, to='ndc')
graphics::par(omd=c(0, 1-w, 0, 1))
graphics::par(mar = c(5, 5, if_else(subtitle==TRUE, 7.9, 3), 3))
text.lbls <- substr(rownames(mSetObj$dataSet$norm), 1, 12)
if (cls.type == "integer") {
cls <- as.factor(as.numeric(levels(cls1))[cls1])
}
else {
cls <- cls1
}
lvs <- levels(cls)
pts.array <- array(0, dim = c(100, 2, length(lvs)))
for (i in 1:length(lvs)) {
inx <- cls1 == lvs[i]
groupVar <- stats::var(cbind(lv1[inx], lv2[inx]), na.rm = T)
groupMean <- cbind(mean(lv1[inx], na.rm = T), mean(lv2[inx],
na.rm = T))
pts.array[, , i] <- ellipse::ellipse(groupVar, centre = groupMean,
level = reg, npoints = 100)
}
xrg <- range(lv1, pts.array[, 1, ])
yrg <- range(lv2, pts.array[, 2, ])
x.ext <- (xrg[2] - xrg[1])/12
y.ext <- (yrg[2] - yrg[1])/12
xlims <- c(xrg[1] - x.ext, xrg[2] + x.ext)
ylims <- c(yrg[1] - y.ext, yrg[2] + y.ext)
cols <- GetColorSchema(cls, grayscale = F)
uniq.cols <- unique(cols)
plot(lv1, lv2, xlab = xlabel, xlim = xlims, ylim = ylims,
ylab = ylabel, type = "n",
cex.lab = 2,
cex.axis = 2,
cex.main = 2,
cex.sub = 2)
graphics::title("PLS-DA Scores Plot", line = if_else(subtitle==TRUE, 4.5, 1), cex.main = 2.5)
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::mtext(line = 3.9, padj = 1, cex=1.7, if_else(subtitle==TRUE, stringr::str_wrap(subtitle.text, 55), NULL))
graphics::grid(col = "lightgray", lty = "dotted", lwd = 1)
if (length(uniq.cols) > 1) {
names(uniq.cols) <- legend.nm
}
for (i in 1:length(lvs)) {
if (length(uniq.cols) > 1) {
graphics::polygon(pts.array[, , i], col = grDevices::adjustcolor(uniq.cols[lvs[i]],
alpha = 0.2), border = NA)
}
else {
graphics::polygon(pts.array[, , i], col = grDevices::adjustcolor(uniq.cols,
alpha = 0.2), border = NA)
}
if (grey.scale) {
graphics::lines(pts.array[, , i], col = grDevices::adjustcolor("black",
alpha = 0.5), lty = 2)
}
}
pchs <- GetShapeSchema(mSetObj, show, grey.scale)
if (grey.scale) {
cols <- rep("black", length(cols))
}
if (show == 1) {
graphics::text(lv1, lv2, label = text.lbls, pos = 4, xpd = T, cex = 0.75)
graphics::points(lv1, lv2, pch = pchs, col = cols)
}
else {
if (length(uniq.cols) == 1) {
graphics::points(lv1, lv2, pch = pchs, col = cols, cex = 1)
}
else {
if (grey.scale == 1 | (exists("shapeVec") &&
all(shapeVec >= 0))) {
my.cols <- grDevices::adjustcolor(cols, alpha.f = 0.4)
my.cols[pchs == 21] <- "black"
graphics::points(lv1, lv2, pch = pchs, col = my.cols, bg = grDevices::adjustcolor(cols,
alpha.f = 0.4), cex = 1.8)
}
else {
graphics::points(lv1, lv2, pch = 21, bg = grDevices::adjustcolor(cols,
alpha.f = 0.4), cex = 2)
}
}
}
uniq.pchs <- unique(pchs)
if (grey.scale) {
uniq.cols <- "black"
}
if (length(uniq.cols) != length(levels(cls))) {
if (cls.type == "integer") {
names(cols) <- as.numeric(cls)
}
else {
names(cols) <- as.character(cls)
}
match.inx <- match(levels(cls), names(cols))
uniq.cols <- cols[match.inx]
}
if (length(uniq.pchs) != length(levels(cls))) {
names(pchs) <- as.character(cls)
match.inx <- match(levels(cls), names(pchs))
uniq.pchs <- pchs[match.inx]
}
graphics::legend(graphics::par('usr')[2], graphics::par('usr')[4], bty='n', legend = legend.nm, pch = uniq.pchs,
col = uniq.cols, cex = 1.3, xpd = NA)
graphics::par(graphics::par(mar=c(5,5,3,3)))
}
if (export == TRUE) {
message(paste0("Exporting 2D PLS-DA scores plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[['imgSet']][[paste0('pls.score2d_PC', inx1, "_PC", inx2)]] <- imgName
}
p(); mSetObj[['imgSet']][[paste0('pls.score2d_PC', inx1, "_PC", inx2, ".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
if (plot == TRUE){
p()
}
return(mSetObj)
}
#' Plot PLS-DA loadings and also set up the matrix for display
#'
#' \code{met.plot_PLS2DLoading} visualizes which features are responsible for the patterns seen among the samples in partial least squares-discriminant analysis.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param subtitle (Logical, \code{TRUE} or \code{FALSE}) Shall the applied data transformation and scaling methods be displayed below the plot title?
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param inx1 Numeric, indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 Numeric, indicate the number of the principal component for the y-axis of the loading plot.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.loading_PCx_PCy.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLSLoading}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS2DLoading <- function (mSetObj = NA, imgName = "PLSDA_2DLoadings", format = "pdf", dpi = NULL, subtitle=FALSE,
width = NA, inx1 = 1, inx2 = 2, plot = TRUE, export = TRUE)
{
load1 <- mSetObj$analSet$plsr$loadings[, inx1]
load2 <- mSetObj$analSet$plsr$loadings[, inx2]
loadings = as.matrix(cbind(load1, load2))
ord.inx <- order(-abs(loadings[, 1]), -abs(loadings[, 2]))
loadings <- signif(loadings[ord.inx, ], 5)
ldName1 <- paste(
"Loadings Comp. ",
inx1,
"(",
round(
100 *
mSetObj$analSet$plsr$Xvar[inx1] /
mSetObj$analSet$plsr$Xtotvar,
1
),
"%)"
)
ldName2 <- paste(
"Loadings Comp. ",
inx2,
"(",
round(
100 *
mSetObj$analSet$plsr$Xvar[inx2] /
mSetObj$analSet$plsr$Xtotvar,
1
),
"%)"
)
colnames(loadings) <- c(ldName1, ldName2)
mSetObj$analSet$plsr$imp.loads <- loadings
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 10
}
else if (width == 0) {
w <- 10
}
else {
w <- width
}
h <- 0.85*w
plotType <- mSetObj$analSet$plsr$loading.type
p <- function(){
graphics::par(mar = c(6, 6, if_else(subtitle==TRUE, 7.9, 3), 6), mgp=c(4.5, 1, 0))
plot(loadings[, 1], loadings[, 2], las = 2, xlab = ldName1,
ylab = ldName2,
cex.lab = 2,
cex.axis = 2,
cex.main = 2,
cex.sub = 2)
graphics::title("PLS-DA Loading Scatterplot", line = if_else(subtitle==TRUE, 5, 1), cex.main = 2.3)
subtitle.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::mtext(line = 3.9, padj = 1, cex=1.7, if_else(subtitle==TRUE, stringr::str_wrap(subtitle.text, 55), NULL))
mSetObj$pls.axis.lims <- graphics::par("usr")
graphics::grid(col = "lightgray", lty = "dotted", lwd = 1)
graphics::points(loadings[, 1], loadings[, 2], pch = 19, col = grDevices::adjustcolor("#CC6663",
alpha.f = 0.4))
if (plotType == "all") {
graphics::text(loadings[, 1], loadings[, 2], labels = substr(rownames(loadings),
1, 16), pos = 4, col = "gray30", xpd = T)
}
else if (plotType == "custom") {
if (length(mSetObj$custom.cmpds) > 0) {
hit.inx <- colnames(mSetObj$dataSet$norm) %in% mSetObj$custom.cmpds
graphics::text(loadings[hit.inx, 1], loadings[hit.inx, 2],
labels = rownames(loadings)[hit.inx], pos = 4,
col = "gray30", xpd = T)
}
}
}
if (export == TRUE) {
message(paste0("Exporting 2D PLS-DA loadings plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[['imgSet']][[paste0('pls.loading_PC', inx1, "_PC", inx2)]] <- imgName
}
p(); mSetObj[['imgSet']][[paste0('pls.loading_PC', inx1, "_PC", inx2, ".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
if (plot == TRUE){
p()
}
return(mSetObj)
}
#' Plot PLS important features
#'
#' \code{met.plot_PLS_Imp} visualizes the features' importance in the PLS-DA model.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (\code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param type (Character) Indicate the variables of importance type to use, \code{"vip"} to use VIP scores, or \code{"coef"} for coefficients
#' @param feat.nm (Character) Indicate the name of the feature. If \code{type = "vip"}, choose \code{"Comp. 1"}, \code{"Comp. 2"}, etc., depending on the component for which the VIP scores should be shown. If \code{type = "coef"}, choose \code{"coef.mean"} for average coefficients, or the name of a sample group in your dataset.
#' @param feat.num (Numeric) Indicate the number of features to show in the plot.
#' @param color.BW (Logical) \code{TRUE} to use black and white, or \code{FALSE} to not.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param title (Logical) \code{TRUE} to add a title with the used normalization, transformation, and scaling method, or \code{FALSE} to not add any title.
#' @return The input mSet object with added dot plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.imp_type_feat.nm.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLS.Impg}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS_Imp <- function (mSetObj = NA, imgName, format = "pdf", dpi = NULL,
width = NA, type = "vip", feat.nm = c("Comp. 1", "Comp. 2"), feat.num = 15,
color.BW = FALSE, plot = TRUE, export = TRUE, title=FALSE)
{
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 8
}
else if (width == 0) {
w <- 7
}
else {
w <- width
}
h <- w
p <- function(){
if (type == "vip") {
mSetObj$analSet$plsda$imp.type <- "vip"
vips <- mSetObj$analSet$plsda$vip.mat[, feat.nm]
met.plot_ImpVar(mSetObj, vips, paste0("VIP scores - ", feat.nm), feat.num,
color.BW, title = title)
}
else {
mSetObj$analSet$plsda$imp.type <- "coef"
data <- mSetObj$analSet$plsda$coef.mat[, feat.nm]
met.plot_ImpVar(mSetObj, data, paste0("Coefficients - ", feat.nm), feat.num,
color.BW, title = title)
}
}
if (export == TRUE) {
message(paste0("Exporting PLS important features plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj[["imgSet"]][[paste0("pls.imp_", type, "_", feat.nm)]] <- imgName
}
p(); mSetObj[["imgSet"]][[paste0("pls.imp_", type, "_", feat.nm,".plot")]] <- grDevices::recordPlot()
grDevices::dev.off()
return(mSetObj)
}
#' Plot PLS important features
#'
#' \code{met.plot_ImpVar} is an internal function used by \code{met.plot_PLS_Imp} to visualize the features' importance in the PLS-DA model.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imp.vec A vector of importan variables.
#' @param xlbl The x-label
#' @param feat.num (Numeric) Indicate the number of features to show in the plot.
#' @param color.BW (Logical) \code{TRUE} to use black and white, or \code{FALSE} to not.
#' @param title (Logical) \code{TRUE} to add a title with the used normalization, transformation, and scaling method, or \code{FALSE} to not add any title.
#' @return The input mSet object with added plot object.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotImpVar}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
met.plot_ImpVar <- function (mSetObj = NA, imp.vec, xlbl, feat.num = 15, color.BW = FALSE, title=FALSE)
{
cls.len <- length(levels(mSetObj$dataSet$cls))
if (cls.len == 2) {
rt.mrg <- 5
}
else if (cls.len == 3) {
rt.mrg <- 6
}
else if (cls.len == 4) {
rt.mrg <- 7
}
else if (cls.len == 5) {
rt.mrg <- 8
}
else if (cls.len == 6) {
rt.mrg <- 9
}
else {
rt.mrg <- 11
}
op <- graphics::par(mar = c(5, 7, if_else(title==TRUE, 4, 3), rt.mrg))
if (feat.num <= 0) {
feat.num = 15
}
if (feat.num > length(imp.vec)) {
feat.num <- length(imp.vec)
}
imp.vec <- rev(sort(imp.vec))[1:feat.num]
imp.vec <- sort(imp.vec)
mns <- by(mSetObj$dataSet$norm[, names(imp.vec)], mSetObj$dataSet$cls,
function(x) {
apply(x, 2, mean, trim = 0.1)
})
mns <- t(matrix(unlist(mns), ncol = feat.num, byrow = TRUE))
vip.nms <- substr(names(imp.vec), 1, 14)
names(imp.vec) <- NULL
dotcolor <- ifelse(color.BW, "darkgrey", "#585855")
graphics::dotchart(imp.vec, bg = dotcolor, xlab = xlbl, cex = 1.3)
graphics::mtext(side = 2, at = 1:feat.num, vip.nms, las = 2, line = 1)
if(title==TRUE){
title.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::title(stringr::str_wrap(title.text, 55), line = 1, cex.main = 1.3)
}
axis.lims <- graphics::par("usr")
shift <- 2 * graphics::par("cxy")[1]
lgd.x <- axis.lims[2] + shift
x <- rep(lgd.x, feat.num)
y <- 1:feat.num
graphics::par(xpd = T)
nc <- ncol(mns)
colorpalette <- ifelse(color.BW, "Greys", "RdYlBu")
col <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(7, colorpalette))(nc)
if (color.BW)
col <- rev(col)
bg <- matrix("", nrow(mns), nc)
for (m in 1:nrow(mns)) {
bg[m, ] <- (col[nc:1])[rank(mns[m, ])]
}
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.factor(as.numeric(levels(mSetObj$dataSet$cls))[mSetObj$dataSet$cls])
}
else {
cls <- mSetObj$dataSet$cls
}
cls.lbl <- levels(cls)
for (n in 1:ncol(mns)) {
graphics::points(x, y, bty = "n", pch = 22, bg = bg[, n],
cex = 3)
graphics::text(x[1], axis.lims[4], cls.lbl[n], srt = 45, adj = c(0.2,
0.5))
x <- x + shift/1.25
}
col <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(7, colorpalette))(50)
if (color.BW)
col <- rev(col)
nc <- length(col)
x <- rep(x[1] + shift, nc)
shifty <- (axis.lims[4] - axis.lims[3])/3
starty <- axis.lims[3] + shifty
endy <- axis.lims[3] + 2 * shifty
y <- seq(from = starty, to = endy, length = nc)
graphics::points(x, y, bty = "n", pch = 15, col = rev(col), cex = 2)
graphics::text(x[1], endy + shifty/8, "High")
graphics::text(x[1], starty - shifty/8, "Low")
graphics::par(op)
}
#' Plot PLS important features
#'
#' \code{met.plot_PLSImpScatter} visualizes the features' importance in the PLS-DA model by comparing VIP scores and coefficients.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (\code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param feat.nm (Character) Indicate the name of the feature. Choose \code{"coef.mean"} for average coefficients, or the name of a sample group in your dataset.
#' @param vip.nm (Character or a character vector) Choose \code{"Comp. 1"}, \code{"Comp. 2"}, etc., depending on the component for which the VIP scores should be shown. If more than one component are specified (the default is \code{c("Comp. 1", "Comp. 2")}), an average VIP score is calculated.
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param vip.thresh (Numeric) Draw a vertical line in the plot indicating a chosen VIP relevance threshold (the default is \code{1}).
#' @param show.title (Logical) \code{TRUE} to add a title to the plot, or \code{FALSE} to not.
#' @param title (Character) Define the title if \code{show.title = TRUE}
#' @return The input mSet object with added scatter plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.ImpScatter_plot_feat.nm_.plot)}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
#' @import ggplot2
met.plot_PLSImpScatter <- function(mSetObj, imgName = "PLS_ImpScatter", format = "png", dpi = 300, width = NA,
feat.nm="coef.mean", vip.nm = c("Comp. 1", "Comp. 2"), plot = TRUE, export = FALSE, vip.thresh = 1,
show.title = FALSE, title){
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 8
}
else {
w <- width
}
h <- w * if_else(show.title==TRUE, 1.1, 1.0)
df <- data.frame(mSetObj$analSet$plsda$vip.mat, mSetObj$analSet$plsda$coef.mat[rownames(mSetObj$analSet$plsda$vip.mat),], check.names = TRUE)
feat <- make.names(feat.nm)
if(length(vip.nm)>1){
df$vip_avg <- rowMeans(df[,grep(paste0(make.names(vip.nm), collapse="|"), colnames(df))])
vip <- "vip_avg"
} else {
vip <- make.names(vip.nm)
}
p <-
ggplot(df, aes_string(x = vip, y = feat)) +
geom_point(color = "Black", fill="Gray", shape=21, size = 3.5, alpha = 0.6) +
ggrepel::geom_text_repel(
aes(label = rownames(df)),
box.padding = unit(0.25,
"lines"),
point.padding = unit(0.5, "lines")
) +
geom_vline(xintercept = vip.thresh,
linetype = "longdash", alpha = 0.4) +
ylab(if_else(feat.nm=="coef.mean", feat.nm, paste0("Coef. ", feat.nm))) +
theme_bw(base_size = 18) +
theme(axis.text=element_text(size = 20),
axis.title = element_text(size=24,face="bold"))
if(length(vip.nm)>1){
p <- p + xlab(paste0("average VIP (", paste(vip.nm, collapse = ", "), ")"))
} else {
p <- p + xlab(paste0("VIP - ", vip.nm))
}
if(show.title == TRUE){
p <- p + ggtitle(title) + theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
if (export == TRUE){
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
print(p)
grDevices::dev.off()
mSetObj[["imgSet"]][[paste0("pls.ImpScatter_", feat.nm)]] <- imgName
}
mSetObj[["imgSet"]][[paste0("pls.ImpScatter_plot_",feat.nm, ".plot")]] <- p
if (plot == TRUE){
print(p)
}
return(mSetObj)
}
#' Create heat map of hierarchically clustered features and samples
#'
#' \code{met.plot_heatmap} plots a heat map with optional filtering for top features based on results from t-tests/ANOVA, VIP or randomforest.
#'
#' @param mSetObj Input name of the created mSet object.
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (\code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param dataOpt (Character) Set data options. Choose \code{"norm"} to use data after normalization or any other string to use original values.
#' @param scaleOpt (Character) Set value standardization option. \code{"row"} for autoscaling of features, \code{"column"} for autoscaling of samples, \code{"none"} for no autoscaling.
#' @param smplDist (Character) Enter the sample distance method: \code{"euclidean"} for Euclidean distance, \code{"correlation"} for Pearson distance, or \code{"minkowski"} for Minkowski method.
#' @param clstDist (Character) Enter the feature clustering method: \code{"ward.d"}, \code{"average"}, \code{"complete"}, or \code{"single"}.
#' @param palette Input color palette choice: \code{"bwm"} for Blue-White-Red, \code{"gbr"} for Green-Black-Red, \code{"heat"} for Red-Yellow, \code{"topo"} for topological colors, \code{"gray"} for gray-scale.
#' @param viewOpt Set view options, choose \code{"detail"} or \code{"overview"}.
#' @param rowV (Logical) Shall clustering be applied to samples?
#' @param colV (Logical) Shall clustering be applied to features?
#' @param var.inx (Numeric vector) \emph{optional}: Provide indeces of selected features to be plotted.
#' @param border (Logical) Indicate whether or not to show cell-borders.
#' @param grp.ave (Logical) Shall only group averages be shown (as opposed to displaying every sample)?
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @return The input mSet object with added heat map. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$heatmap.avg.plot)} or \code{print(mSetObj$imgSet$heatmap.all.plot)}, depending on \code{grp.ave}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotHeatMap}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_heatmap <- function (mSetObj = NA, imgName = "Heatmap_features", format = "pdf", dpi = NULL,
width = NA, dataOpt = "norm", scaleOpt = "row",
smplDist = "euclidean", clstDist = "ward.D", palette = "bwm",
viewOpt = "detail", rowV = T, colV = T, var.inx = NULL,
border = T, grp.ave = F, plot = TRUE, export = TRUE)
{
imgName = paste(imgName, ".", format,
sep = "")
cls <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
cls.class <- mSetObj$dataSet$type.cls.lbl
mSetObj$analSet$htmap <- list(dist.par = smplDist, clust.par = clstDist)
if (dataOpt == "norm") {
my.data <- mSetObj$dataSet$norm
}
else {
if(!is.null(mSetObj$dataSet$prenorm)){
row.norm <- mSetObj$dataSet$prenorm
} else {
my.data <- qs::qread("prenorm.qs")
}
}
if (is.null(var.inx)) {
hc.dat <- as.matrix(my.data)
}
else {
hc.dat <- as.matrix(my.data[, var.inx])
}
colnames(hc.dat) <- substr(colnames(hc.dat), 1, 18)
if (cls.class == "integer") {
hc.cls <- as.factor(as.numeric(levels(cls))[cls])
}
else {
hc.cls <- cls
}
if (grp.ave) {
lvs <- levels(hc.cls)
my.mns <- matrix(ncol = ncol(hc.dat), nrow = length(lvs))
for (i in 1:length(lvs)) {
inx <- hc.cls == lvs[i]
my.mns[i, ] <- apply(hc.dat[inx, ], 2, mean)
}
rownames(my.mns) <- lvs
colnames(my.mns) <- colnames(hc.dat)
hc.dat <- my.mns
hc.cls <- as.factor(lvs)
}
if (palette == "gbr") {
colors <- grDevices::colorRampPalette(c("green", "black",
"red"), space = "rgb")(256)
}
else if (palette == "heat") {
colors <- grDevices::heat.colors(256)
}
else if (palette == "topo") {
colors <- grDevices::topo.colors(256)
}
else if (palette == "gray") {
colors <- grDevices::colorRampPalette(c("grey90", "grey10"),
space = "rgb")(256)
}
else {
colors <- rev(grDevices::colorRampPalette(RColorBrewer::brewer.pal(10,
"RdBu"))(256))
}
if (cls.type == "disc") {
annotation <- data.frame(class = hc.cls)
rownames(annotation) <- rownames(hc.dat)
}
else {
annotation <- NA
}
plot_dims <- get_pheatmap_dims(t(hc.dat), annotation, viewOpt,
width)
h <- plot_dims$height
w <- plot_dims$width
if (grp.ave) {
w <- nrow(hc.dat) * 25 + 300
w <- round(w/72, 2)
}
if (border) {
border.col <- "grey60"
}
else {
border.col <- NA
}
p <- function(){
if (cls.type == "disc") {
cols <- GetColorSchema(cls, palette == "gray")
uniq.cols <- unique(cols)
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.factor(as.numeric(levels(cls))[cls])
}
else {
cls <- cls
}
names(uniq.cols) <- unique(as.character(sort(cls)))
ann_colors <- list(class = uniq.cols)
pheatmap::pheatmap(t(hc.dat), annotation = annotation,
annotation_colors = ann_colors, fontsize = 8, fontsize_row = 8,
clustering_distance_rows = smplDist, clustering_distance_cols = smplDist,
clustering_method = clstDist, border_color = border.col,
cluster_rows = colV, cluster_cols = rowV, scale = scaleOpt,
color = colors)
}
else {
stats::heatmap(hc.dat, Rowv = rowTree, Colv = colTree, col = colors,
scale = "column")
}
}
if (export == TRUE) {
message(paste0("Exporting heat map of features to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$heatmap <- imgName
}
if (grp.ave == TRUE) {
p(); mSetObj$imgSet$heatmap.avg.plot <- grDevices::recordPlot()
} else {
p(); mSetObj$imgSet$heatmap.all.plot <- grDevices::recordPlot()
}
if (plot == TRUE) {
p()
}
return(mSetObj)
}
#' Volcano Plot
#'
#' \code{met.plot_volcano} generates a volcano plot for a specified contrast.
#'
#' @param mSetObj Input name of the created mSet object,
#' Data container after ANOVA analysis (\code{\link[VisomX]{met.ANOVA.Anal}}) if \code{test = "anova"}.
#' @param grp1 (Character) Enter name of the first group for the contrast \code{grp1 vs. grp2}. If both group arguments are empty, the first two names in the list of groups are selected.
#' @param grp2 (Character) Enter name of the second group for the contrast \code{grp1 vs. grp2}. If both group arguments are empty, the first two names in the list of groups are selected.
#' @param test (Character) Choose a statistical tests. For \code{test = "ttest"}, \code{met.plot_volcano} runs \code{\link[VisomX]{met.Ttests.Anal}} with the chosen test parameters.
#' For \code{test = "anova"}, \code{\link[VisomX]{met.ANOVA.Anal}} must have been applied previously on the mSetObj.
#' @param paired (Logical) Is the data paired (\code{TRUE}) or not (\code{FALSE}). Only applicable for \code{test = "ttest"}.
#' @param nonpar (Logical) Use a non-parametric test (\code{TRUE}) or not (\code{FALSE}). Only applicable for \code{test = "ttest"}.
#' @param equal.var (Logical) The two groups have equal variance (\code{TRUE}) or not (\code{FALSE}). Only applicable for \code{test = "ttest"}.
#' @param log2fc.thresh (Numeric) Enter a relevance threshold for log2 fold changes, highlighted in the plot by vertical lines and colored compounds
#' @param threshp (Numeric) Enter a significance threshold for features based on T-test or ANOVA test results, highlighted in the plot by an horizontal line and colored compounds
#' @param pval.type (Character) Display and apply significance threshold to \code{"raw"} p values or adjusted p values (\code{"fdr"}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}). If \code{NULL}, the name \code{Plots/Volcano_grp1_vs_grp2} is assigned.
#' @param format (\code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param add_names (Logical) Display labels of significant features (\code{TRUE}) or not (\code{FALSE}).
#' @param label_size (Numeric) Font size for feature labels (if \code{add_names = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param silent (Logical) Shall the results of \code{\link[VisomX]{met.Ttests.Anal}} be printed in the console (\code{TRUE}) or not (\code{FALSE})?
#' @param test_condition (Logical) Add a subtitle with the applied data transformation and scaling methods be displayed below the plot title (\code{TRUE}) or not (\code{FALSE}).
#' @return The input mSet object with added volcano plot (generated by \code{\link[ggplot2]{ggplot}}). The plot can be retrieved from within R via \code{print(mSetObj$imgSet$volcano$grp1_vs_grp2.plot)}.
#' @references adapted from \code{\link[DEP]{plot_volcano}} (\url{https://bioconductor.org/packages/DEP/}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @import ggplot2
#' @export
met.plot_volcano <- function (mSetObj = NA, grp1, grp2, test = "ttest", paired = FALSE, nonpar = FALSE, equal.var = TRUE,
log2fc.thresh = 1, threshp = 0.05, pval.type = "fdr", imgName = NULL, format = "pdf",
add_names = TRUE, label_size = 3, dpi = NULL, width = NA, plot = TRUE, export = TRUE,
silent = FALSE, test_condition = FALSE)
{
if(!is.null(grp1) && !is.null(grp2)){
assertthat::assert_that(grp1 %in% levels(mSetObj$dataSet$cls),
grp2 %in% levels(mSetObj$dataSet$cls),
msg = paste0("'",grp1, "' or '", grp2, "' are not valid conditions in the dataset. Valid conditions are:\n",
paste(levels(mSetObj$dataSet$cls), collapse = ", ")))
}
if(is.null(grp1) && is.null(grp2)){
grp1 <- unique(levels(mSetObj$dataSet$cls))[1]
grp2 <- unique(levels(mSetObj$dataSet$cls))[2]
}
if(test=="ttest"){
mSetObj <- met.Ttests.Anal(mSetObj, grp1 = grp1, grp2 = grp2, nonpar, threshp, paired,
equal.var, pval.type, silent = silent)
p.value <- mSetObj$analSet$tt[[paste0(grp1,"_vs_",grp2)]][["p.value"]]
if (pval.type == "fdr") {
p.value <- p.adjust(p.value, "fdr")
}
} else if(test=="anova"&&!is.null(mSetObj$analSet$aov)){
p.value <- stats::setNames(mSetObj$analSet$aov$post.hoc.all[,'p.adj'], rownames(mSetObj$analSet$aov$post.hoc.all))
p.value <- sapply(1:length(p.value), function (x) str_split(p.value[x], "; "), simplify = TRUE)
p.value <-
stats::setNames(as.numeric(sapply(p.value, "[", match(
x = paste0(grp2, "-", grp1), table = unlist(str_split(
mSetObj$analSet$aov$post.hoc.all$contrasts[1], "; "
))
))),
rownames(mSetObj$analSet$aov$post.hoc.all))
} else if(test=="anova"&&is.null(mSetObj$analSet$aov)){
stop("Please perform 'met.ANOVA.Anal()' on your mSet before running 'met.plot_volcano' with 'test = \"anova\"'.")
}
inx.p <- p.value <= threshp
p.log <- -log10(p.value)
mSetObj <- met.FC.Anal(mSetObj, log2fc.thresh = log2fc.thresh, grp1 = grp1, grp2 = grp2)
max.xthresh <- log2fc.thresh
min.xthresh <- -log2fc.thresh
fc.log <- mSetObj$analSet$fc[[paste0(grp1,"_vs_",grp2)]][["fc.log"]]
fc.all <- mSetObj$analSet$fc[[paste0(grp1,"_vs_",grp2)]][["fc.all"]]
inx.up <- mSetObj$analSet$fc[[paste0(grp1,"_vs_",grp2)]][["inx.up"]]
inx.down <- mSetObj$analSet$fc[[paste0(grp1,"_vs_",grp2)]][["inx.down"]]
keep.inx <- names(inx.p) %in% names(inx.up)
inx.p <- inx.p[keep.inx]
p.log <- p.log[keep.inx]
inx.imp <- (inx.up | inx.down) & inx.p
sig.var <- cbind(fc.all[inx.imp, drop = F], fc.log[inx.imp,
drop = F], p.value[inx.imp, drop = F], p.log[inx.imp,
drop = F])
if (pval.type == "fdr") {
colnames(sig.var) <- c("FC", "log2(FC)",
"p.adjusted", "-log10(p)")
}
else {
colnames(sig.var) <- c("FC", "log2(FC)",
"raw.pval", "-log10(p)")
}
ord.inx <- order(sig.var[, 4], abs(sig.var[, 2]), decreasing = T)
sig.var <- sig.var[ord.inx, , drop = F]
sig.var <- signif(sig.var, 5)
fileName <- "volcano.csv"
fast.write.csv(signif(sig.var, 5), file = fileName)
volcano <- list(pval.type = pval.type, raw.threshx = log2fc.thresh,
raw.threshy = threshp, paired = paired, max.xthresh = max.xthresh,
min.xthresh = min.xthresh, thresh.y = -log10(threshp),
fc.all = fc.all, fc.log = fc.log, inx.up = inx.up, inx.down = inx.down,
p.log = p.log, inx.p = inx.p, sig.mat = sig.var)
if(!exists("volcano", mSetObj$analSet)){
mSetObj$analSet$volcano <- list()
}
mSetObj$analSet$volcano[[paste0(grp1, "_vs_", grp2)]] <- volcano
if(!is.character(imgName)){
imgName <- paste0("Plots/Volcano_", grp1, "_vs_", grp2)
}
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 9
}
else if (width == 0) {
w <- 8
}
else {
w <- width
}
h <- w * 5.5/6.5
if(!exists("volcano", mSetObj$imgSet)){
mSetObj$imgSet$volcano <- list()
}
vcn <- mSetObj$analSet$volcano[[paste0(grp1, "_vs_", grp2)]]
imp.inx <- (vcn$inx.up | vcn$inx.down) & vcn$inx.p
de <- data.frame(cbind(vcn$fc.log, vcn$p.log))
de$Status <- "Non-SIG"
de$Status[vcn$inx.p & vcn$inx.up] <- "UP"
de$Status[vcn$inx.p & vcn$inx.down] <- "DOWN"
de$Status <- as.factor(de$Status)
mycols <- levels(de$Status)
mycols[mycols == "UP"] <- "#e31f26"
mycols[mycols == "DOWN"] <- "#387fb9"
mycols[mycols == "Non-SIG"] <- "#525352"
de$delabel <- NA
de$delabel[imp.inx] <- rownames(de)[imp.inx]
if (pval.type == "fdr") {
de$shape <- ifelse(de[,2] > 5.9, "triangle", "circle")
} else {
de$shape <- "circle"
}
if (pval.type == "fdr") {
if (max(de[,2]) > 6.0) {
de[,2][de[,2] > 5.9] <- 5.9
}
}
title_text <- paste0(grp1, " vs. ", grp2, if_else(
test_condition == TRUE,
paste0(
"\n(",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", "none"),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", "none"),
")"
),
"")
)
p <- ggplot(data = de, aes(x = de[, 1], y = de[, 2],
col = Status, label = delabel)) +
geom_point(aes(shape = shape), size = 2.5, alpha = 0.6) +
scale_shape(guide = "none") + # Hide shapes from legend
scale_color_manual(labels = c("Increased", "Not significant", "Decreased"),
values = mycols) +
geom_vline(xintercept = c(-vcn$raw.threshx, vcn$raw.threshx),
linetype = "longdash", alpha = 0.4) +
geom_hline(yintercept = -log10(vcn$raw.threshy), linetype = "longdash",
alpha = 0.4) +
labs(title = title_text) +
labs(x = "log2(fold change)", y = "-log10(p value)") +
theme_bw(base_size = 20) +
theme(legend.position = "bottom",title = element_text(size = exp(-nchar(paste0(grp1, " vs. ", grp2))/40)*if_else(
test_condition == TRUE, 25, 30)),
axis.title = element_text(size = 22), legend.title = element_text(size = 22), plot.title = element_text(hjust = 0.5)) +
ggh4x::force_panelsizes(rows = unit(if_else(test_condition==TRUE, 0.57, 0.61)*w, "in"),
cols = unit(if_else(test_condition==TRUE, 0.67, 0.72)*w, "in")) +
coord_cartesian(xlim = c(-max(abs(de[, 1])) - 0.05*max(abs(de[, 1])),
max(abs(de[, 1])) + 0.05*max(abs(de[, 1]))),
expand = TRUE)
# scale_x_continuous(breaks = scales::pretty_breaks(n = max(abs(de[,1])) + 1))
if (add_names) {
p <- p + ggrepel::geom_text_repel(
aes(label = delabel),
size = label_size,
box.padding = unit(0.25,
"lines"),
point.padding = unit(0.1, "lines"),
segment.size = 0.5,
show.legend = FALSE,
max.overlaps = 16
)
}
if (pval.type == "fdr") {
p <- p + labs(y = expression(-log[10] ~ "(adj. p-value)"))
}
if (export == TRUE){
message(paste0("Exporting volcano plot to:\n\"", getwd(), imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
print(p)
grDevices::dev.off()
mSetObj$imgSet$volcano[[paste0(grp1, "_vs_", grp2)]] <- imgName
}
mSetObj$imgSet$volcano[[paste0(grp1, "_vs_", grp2, ".plot")]] <- p
if (plot == TRUE){
print(p)
}
return(mSetObj)
}
#' Create 3D PCA loading plot
#'
#' \code{met.plot_PCA3DLoading} visualizes which features are responsible for the patterns seen among the samples in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format format.
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the loading plot.
#' @param inx3 (Numeric) Indicate the number of the principal component for the z-axis of the loading plot.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as file with the chosen format?
#' @return The input mSet object with added 3D scatter plot. The plot can be retrieved from within R by executing \code{met.print_PCA3DLoading(mSetObj$imgSet$pca.loading3d.plot)}.
#' @export
met.plot_PCA3DLoading <- function (mSetObj = NA, imgName = "PCA3DLoading", format = "json", inx1 = 1,
inx2 = 2, inx3 = 3, export = F)
{
cls <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
cls.class <- mSetObj$dataSet$type.cls.lbl
pca <- mSetObj$analSet$pca
pca3d <- list()
pca3d$loading$axis <- paste("Loading PC", c(inx1, inx2,
inx3), sep = "")
coords <- data.frame(t(signif(pca$rotation[, 1:3], 5)))
dists <- apply(coords, 2, function(x) dist(rbind(x, c(0, 0, 0))))
pca3d$loading$cols <- GetRGBColorGradient(dists)
cols_vectors <- stringr::str_replace_all(pca3d$loading$cols, "rgba\\(", "") %>% stringr::str_replace_all(., "\\)", "")
cols_vectors <- sapply(1:length(cols_vectors), function (x) str_split(cols_vectors[x], ","))
pca3d$loading$cols_hex <-
sapply(1:length(cols_vectors), function (x)
grDevices::rgb(cols_vectors[[x]][1], cols_vectors[[x]][2], cols_vectors[[x]][3],
alpha = cols_vectors[[x]][4],
maxColorValue = 255
))
colnames(coords) <- NULL
pca3d$loading$xyz <- coords
pca3d$loading$name <- rownames(pca$rotation)
pca3d$loading$entrez <- rownames(pca$rotation)
if (cls.class == "integer") {
clss <- as.character(sort(as.factor(as.numeric(levels(cls))[cls])))
}
else {
clss <- as.character(cls)
}
if (all.numeric(clss)) {
clss <- paste("Group", clss)
}
pca3d$cls = clss
imgName = paste(imgName, ".", format, sep = "")
if(export==TRUE){
message(paste0("Exporting 3D PCA loading plot to:\n\"", getwd(), imgName, "\""))
json.mat <- jsonlite::toJSON(pca3d)
sink(imgName)
cat(json.mat)
sink()
mSetObj$imgSet$pca.loading3d <- imgName
}
mSetObj$imgSet$pca.loading3d.plot <- pca3d
return(mSetObj)
}
#' Plot a generated 3D PCA loading plot in RStudio
#'
#' \code{met.print_PCA3DLoading} opens a new graphics device to show the generated 3D PCA loading plot.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after \code{\link[VisomX]{met.plot_PCA3DLoading}}.
#' @return A 3D scatter plot. To close the graphics device, execute \code{rgl::close3d()}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.print_PCA3DLoading <- function (mSetObj = NA)
{
colors <- mSetObj$imgSet$pca.loading3d.plot$loading$cols_hex
rgl::plot3d(x= mSetObj$imgSet$pca.loading3d.plot$loading$xyz[1,],
y= mSetObj$imgSet$pca.loading3d.plot$loading$xyz[3,],
z= mSetObj$imgSet$pca.loading3d.plot$loading$xyz[2,],
col = colors,
xlab = mSetObj$imgSet$pca.loading3d.plot$loading$axis[1],
ylab = mSetObj$imgSet$pca.loading3d.plot$loading$axis[3],
zlab = mSetObj$imgSet$pca.loading3d.plot$loading$axis[2],
type="s" , size = 1.5)
rgl::text3d(x=mSetObj$imgSet$pca.loading3d.plot$loading$xyz[1,],
y=mSetObj$imgSet$pca.loading3d.plot$loading$xyz[3,],
z=mSetObj$imgSet$pca.loading3d.plot$loading$xyz[2,],
mSetObj$imgSet$pca.loading3d.plot$loading$name, cex=0.6, adj = c(1.3, 1.1))
}
#' Create 3D PCA Score plot
#'
#' \code{met.plot_PCA3DScore} visualizes clusters of samples based on their similarity in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after principal component analysis (\code{\link[VisomX]{met.PCA.Anal}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format format.
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the score plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the score plot.
#' @param inx3 (Numeric) Indicate the number of the principal component for the z-axis of the score plot.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as file with the chosen format?
#' @return The input mSet object with added 3D scatter plot. The plot can be retrieved from within R by executing \code{met.print_PCA3DScore(mSetObj$imgSet$pca.score3d.plot)}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PCA3DScore <- function (mSetObj = NA, imgName = "PCA3DScore", format = "json", inx1 = 1,
inx2 = 2, inx3 = 3, export = F)
{
cls <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
cls.class <- mSetObj$dataSet$type.cls.lbl
pca <- mSetObj$analSet$pca
pca3d <- list()
pca3d$score$axis <- paste("PC", c(inx1, inx2, inx3),
" (", 100 * round(mSetObj$analSet$pca$variance[c(inx1,
inx2, inx3)], 3), "%)", sep = "")
coords <- data.frame(t(signif(pca$x[, c(inx1, inx2, inx3)],
5)))
colnames(coords) <- NULL
pca3d$score$xyz <- coords
pca3d$score$name <- rownames(mSetObj$dataSet$norm)
if (cls.class == "integer") {
cls <- as.character(sort(as.factor(as.numeric(levels(cls))[cls])))
}
else {
cls <- as.character(cls)
}
if (all.numeric(cls)) {
cls <- paste("Group", cls)
}
pca3d$score$facA <- cls
cols <- unique(GetColorSchema(as.factor(cls)))
pca3d$score$colors <- apply(col2rgb(cols), 2, function(x) {
paste("rgb(", paste(x, collapse = ","), ")", sep = "")
})
imgName = paste(imgName, ".", format, sep = "")
if(export==TRUE){
message(paste0("Exporting 3D PCA score plot to:\n\"", getwd(), imgName, "\""))
json.obj <- jsonlite::toJSON(pca3d)
sink(imgName)
cat(json.obj)
sink()
mSetObj$imgSet$pca.score3d <- imgName
}
mSetObj$imgSet$pca.score3d.plot <- pca3d
return(mSetObj)
}
#' Plot a generated 3D PCA scores plot in RStudio
#'
#' \code{met.print_PCA3DScore} opens a new graphics device to show the generated 3D PCA score plot.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after \code{\link[VisomX]{met.plot_PCA3DScore}}.
#' @return A 3D scatter plot. To close the graphics device, execute \code{rgl::close3d()}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.print_PCA3DScore <- function (mSetObj = NA)
{
colors <- GetColorSchema(as.factor(cls))
rgl::plot3d(x= mSetObj$imgSet$pca.score3d.plot$score$xyz[1,],
y= mSetObj$imgSet$pca.score3d.plot$score$xyz[3,],
z= mSetObj$imgSet$pca.score3d.plot$score$xyz[2,],
col = colors,
xlab = mSetObj$imgSet$pca.score3d.plot$score$axis[1],
ylab = mSetObj$imgSet$pca.score3d.plot$score$axis[3],
zlab = mSetObj$imgSet$pca.score3d.plot$score$axis[2],
type="s" , size = 1.5)
rgl::par3d(windowRect = c(0, 0, 512, 512))
rgl::legend3d(
"right",
legend = levels(factor(
mSetObj$imgSet$pca.score3d.plot$score$facA
)),
col = unique(colors),
pch = 16
)
}
#' Create 3D PLS-DA loading plot
#'
#' \code{met.plot_PLS3DLoading} visualizes which features are responsible for the patterns seen among the samples in partial least squares-discriminant analysis.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format format.
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the loading plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the loading plot.
#' @param inx3 (Numeric) Indicate the number of the principal component for the z-axis of the loading plot.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as file with the chosen format?
#' @return The input mSet object with added 3D scatter plot. The plot can be retrieved from within R by executing \code{met.print_PLS3DLoading(mSetObj$imgSet$pls.loading3d.plot)}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS3DLoading <- function (mSetObj = NA, imgName = "PLS3DLoading", format = "json", inx1 = 1,
inx2 = 2, inx3 = 3, export=F)
{
pls = mSetObj$analSet$plsr
coords <- signif(as.matrix(cbind(pls$loadings[, inx1], pls$loadings[,
inx2], pls$loadings[, inx3])), 5)
pls3d <- list()
pls3d$loading$axis <- paste("Loading Comp.", c(inx1, inx2,
inx3), sep = "")
coords0 <- coords <- data.frame(t(signif(pls$loadings[, c(inx1,
inx2, inx3)], 5)))
colnames(coords) <- NULL
pls3d$loading$xyz <- coords
pls3d$loading$name <- rownames(pls$loadings)
pls3d$loading$entrez <- rownames(pls$loadings)
dists <- apply(coords0, 2, function(x) dist(rbind(x, c(0, 0, 0))))
cols <- GetRGBColorGradient(dists)
pls3d$loading$cols <- cols
cols_vectors <- stringr::str_replace_all(cols, "rgba\\(", "") %>% stringr::str_replace_all(., "\\)", "")
cols_vectors <- sapply(1:length(cols_vectors), function (x) str_split(cols_vectors[x], ","))
pls3d$loading$cols_hex <-
sapply(1:length(cols_vectors), function (x)
grDevices::rgb(cols_vectors[[x]][1], cols_vectors[[x]][2], cols_vectors[[x]][3],
alpha = cols_vectors[[x]][4],
maxColorValue = 255
))
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.character(sort(as.factor(as.numeric(levels(mSetObj$dataSet$cls))[mSetObj$dataSet$cls])))
}
else {
cls <- as.character(mSetObj$dataSet$cls)
}
if (all.numeric(cls)) {
cls <- paste("Group", cls)
}
pls3d$cls = cls
imgName = paste(imgName, ".", format, sep = "")
if(export==TRUE){
message(paste0("Exporting 3D PLS-DA loading plot to:\n\"", getwd(), imgName, "\""))
json.mat <- jsonlite::toJSON(pls3d)
sink(imgName)
cat(json.mat)
sink()
mSetObj$imgSet$pls.loading3d <- imgName
}
mSetObj$imgSet$pls.loading3d.plot <- pls3d
return(mSetObj)
}
#' Plot a generated 3D PLS-DA loading plot in RStudio
#'
#' \code{met.print_PLS3DLoading} opens a new graphics device to show the generated 3D PLS-DA loading plot.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after \code{\link[VisomX]{met.plot_PLS3DLoading}}.
#' @return A 3D scatter plot. To close the graphics device, execute \code{rgl::close3d()}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.print_PLS3DLoading <- function (mSetObj = NA)
{
colors <- mSetObj$imgSet$pls.loading3d.plot$loading$cols_hex
rgl::plot3d(x= mSetObj$imgSet$pls.loading3d.plot$loading$xyz[1,],
y= mSetObj$imgSet$pls.loading3d.plot$loading$xyz[3,],
z= mSetObj$imgSet$pls.loading3d.plot$loading$xyz[2,],
col = colors,
xlab = mSetObj$imgSet$pls.loading3d.plot$loading$axis[1],
ylab = mSetObj$imgSet$pls.loading3d.plot$loading$axis[3],
zlab = mSetObj$imgSet$pls.loading3d.plot$loading$axis[2],
type="s" , size = 1.5)
rgl::text3d(x=mSetObj$imgSet$pls.loading3d.plot$loading$xyz[1,],
y=mSetObj$imgSet$pls.loading3d.plot$loading$xyz[3,],
z=mSetObj$imgSet$pls.loading3d.plot$loading$xyz[2,],
mSetObj$imgSet$pls.loading3d.plot$loading$name, cex=0.6, adj = c(1.3, 1.1))
}
#' Create 3D PLS-DA Score plot
#'
#' \code{met.plot_PLS3DScore} visualizes clusters of samples based on their similarity in principal component analysis.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format format.
#' @param inx1 (Numeric) Indicate the number of the principal component for the x-axis of the score plot.
#' @param inx2 (Numeric) Indicate the number of the principal component for the y-axis of the score plot.
#' @param inx3 (Numeric) Indicate the number of the principal component for the z-axis of the score plot.
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as file with the chosen format?
#' @return The input mSet object with added 3D scatter plot. The plot can be retrieved from within R by executing \code{met.print_PLS3DScore(mSetObj$imgSet$pls.score3d.plot)}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS3DScore <- function (mSetObj = NA, imgName = "PLS3DScore", format = "json", inx1 = 1,
inx2 = 2, inx3 = 3, export=F)
{
cls1 <- mSetObj$dataSet$cls
cls.type <- mSetObj$dataSet$cls.type
cls.class <- mSetObj$dataSet$type.cls.lbl
pls3d <- list()
pls3d$score$axis <- paste("Component", c(inx1, inx2,
inx3), " (", round(100 * mSetObj$analSet$plsr$Xvar[c(inx1,
inx2, inx3)]/mSetObj$analSet$plsr$Xtotvar, 1), "%)",
sep = "")
coords <- data.frame(t(signif(mSetObj$analSet$plsr$score[,
c(inx1, inx2, inx3)], 5)))
colnames(coords) <- NULL
pls3d$score$xyz <- coords
pls3d$score$name <- rownames(mSetObj$dataSet$norm)
if (mSetObj$dataSet$type.cls.lbl == "integer") {
cls <- as.character(sort(as.factor(as.numeric(levels(cls1))[cls1])))
}
else {
cls <- as.character(cls1)
}
if (all.numeric(cls)) {
cls <- paste("Group", cls)
}
pls3d$score$facA <- cls
cols <- unique(GetColorSchema(cls1))
pls3d$score$colors <- apply(col2rgb(cols), 2, function(x) {
paste("rgb(", paste(x, collapse = ","), ")", sep = "")
})
imgName = paste(imgName, ".", format, sep = "")
if(export==TRUE){
message(paste0("Exporting 3D PLS-DA score plot to:\n\"", getwd(), imgName, "\""))
json.obj <- jsonlite::toJSON(pls3d)
sink(imgName)
cat(json.obj)
sink()
mSetObj$imgSet$pls.score3d <- imgName
}
mSetObj$imgSet$pls.score3d.plot <- pls3d
return(mSetObj)
}
#' Plot a generated 3D PLS-DA scores plot in RStudio
#'
#' \code{met.print_PLS3DScore} opens a new graphics device to show the generated 3D PLS-DA score plot.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after first after \code{\link[VisomX]{met.plot_PLS3DScore}}.
#' @return A 3D scatter plot. To close the graphics device, execute \code{rgl::close3d()}.
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.print_PLS3DScore <- function (mSetObj = NA)
{
colors <- GetColorSchema(as.factor(cls))
rgl::plot3d(x= mSetObj$imgSet$pls.score3d.plot$score$xyz[1,],
y= mSetObj$imgSet$pls.score3d.plot$score$xyz[3,],
z= mSetObj$imgSet$pls.score3d.plot$score$xyz[2,],
col = colors,
xlab = mSetObj$imgSet$pls.score3d.plot$score$axis[1],
ylab = mSetObj$imgSet$pls.score3d.plot$score$axis[3],
zlab = mSetObj$imgSet$pls.score3d.plot$score$axis[2],
type="s" , size = 1.5)
rgl::par3d(windowRect = c(0, 0, 512, 512))
rgl::legend3d(
"right",
legend = levels(factor(
mSetObj$imgSet$pls.score3d.plot$score$facA
)),
col = unique(colors),
pch = 16
)
}
#' Plot PLS-DA classification performance using different components, permutation
#'
#' Permutation test is a technique for testing a hypothesis of no effect, when the distribution of the test statistic is unknown. The objective of this test is to confirm that the initial model is superior to other models obtained by permuting the class labels and randomly assigning them to different individuals or, in other words, to answer the question "what is the model's performance if the groups are formed randomly".
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param title (Logical) \code{TRUE} to add a title with the used normalization, transformation, and scaling method, or \code{FALSE} to not add any title.
#' @return The input mSet object with added permutation test plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.permut.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLS.Permutation}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS.Permutation <- function (mSetObj = NA, imgName = "PLSDA-Permutation", format = "pdf", dpi = NULL,
width = NA, plot = TRUE, export = TRUE, title = FALSE)
{
bw.vec <- mSetObj$analSet$plsda$permut
len <- length(bw.vec)
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 8
}
else if (width == 0) {
w <- 7
}
else {
w <- width
}
h <- w * 6/8
p <- function(){
graphics::par(mar = c(5, 5, if_else(title==TRUE, 3, 2), 4))
hst <- graphics::hist(
bw.vec,
breaks = "FD",
freq = T,
ylab = "Frequency",
xlab = paste0(
"Permutation test statistics",
if_else(
mSetObj[["analSet"]][["plsda"]][["permut.type"]] == "separation distance",
" (separation distance)",
" (prediction accuracy)"
)
),
col = "#abd9e9",
main = ""
)
h <- max(hst$counts)
graphics::arrows(x0=bw.vec[1]+max(bw.vec)/20, y0=if_else(title==TRUE, h/1.45, h/5),
x1=bw.vec[1], y1=h/35, col ="#d73027", lwd=1.5, length=0.15, angle=25, xpd = T)
if(title==TRUE){
title.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all / that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all / that are not followed by a letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::title(stringr::str_wrap(title.text, 55), line = 1, cex.main = 1.2)
}
graphics::text(x=bw.vec[1]+max(bw.vec)/20, y=if_else(title==TRUE, h/1.15, h/3.5), labels=paste("Observed \n statistic \n",
mSetObj$analSet$plsda$permut.p), xpd = T, cex=0.95)
}
if(export == TRUE){
message(paste0("Exporting plots with PLS-DA model permutation results to:\n\"", getwd(), "/", imgName, "\""))
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$pls.permut <- imgName
}
p(); mSetObj$imgSet$pls.permut.plot <- grDevices::recordPlot()
grDevices::dev.off()
if (plot == TRUE){
p()
}
return(mSetObj)
}
#' Plot PLS-DA classification performance using different components
#'
#' \code{met.plot_PLS.Crossvalidation} performs cross-validation (CV) on the generated PLS-DA model, where a fraction of data is held back, and the model trained on the rest. In each CV, the predicted data are compared with the original data, and the sum of squared errors is calculated. The prediction error is then summed over all samples (Predicted Residual Sum of Squares or PRESS). For convenience, the PRESS is divided by the initial sum of squares and subtracted from 1 to resemble the scale of the R2. Good predictions will have low PRESS or high Q2. Generally speaking, a model with an R2 (and Q2) value above 0.7 can be considered predictive. It is possible to have negative Q2, which means that your model is not at all predictive or is overfitted.
#'
#' @param mSetObj Input name of the created mSet object.
#' Data container after partial least squares-discriminant analysis (\code{\link[VisomX]{met.PLSR.Anal}} and \code{\link[VisomX]{met.PLSDA.CV}}).
#' @param imgName (Character) Enter a name for the image file (if \code{export = TRUE}).
#' @param format (Character, \code{"png"} or \code{"pdf"}) image file format (if \code{export = TRUE}).
#' @param dpi (Numeric) resolution of the image file (if \code{export = TRUE}). If \code{NULL}, the resolution will be chosen automatically based on the chosen file format (300 dpi for PNG, 72 dpi for PDF)
#' @param width (Numeric) width of the the image file in inches (if \code{export = TRUE}).
#' @param plot (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be returned in the RStudio 'Plots' pane?
#' @param export (Logical, \code{TRUE} or \code{FALSE}) Shall the plot be exported as PDF or PNG file?
#' @param title (Logical) \code{TRUE} to add a title with the used normalization, transformation, and scaling method, or \code{FALSE} to not add any title.
#' @return The input mSet object with added cross validation test plot. The plot can be retrieved from within R via \code{print(mSetObj$imgSet$pls.crossvalidation.plot)}.
#' @references adapted from \code{\link[MetaboAnalystR]{PlotPLS.Classification}} (\url{https://github.com/xia-lab/MetaboAnalystR}).
#' @author Nicolas T. Wirth \email{mail.nicowirth@gmail.com}
#' Technical University of Denmark
#' License: GNU GPL (>= 2)
#' @export
met.plot_PLS.Crossvalidation <- function (mSetObj = NA, imgName = "PLSDA-CrossValidation", format = "pdf", dpi = NULL,
width = NA, plot = TRUE, export = TRUE, title = FALSE)
{
mSetObj <- mSetObj
res <- mSetObj$analSet$plsda$fit.info
colnames(res) <- 1:ncol(res)
best.num <- mSetObj$analSet$plsda$best.num
choice <- mSetObj$analSet$plsda$choice
imgName = paste(imgName, ".", format,
sep = "")
if (is.na(width)) {
w <- 7
}
else if (width == 0) {
w <- 7
}
else {
w <- width
}
h <- w * 5/6.7
p <- function(){
graphics::par(mar = c(5, 5, if_else(title==TRUE, 3, 2), 7))
graphics::barplot(res, beside = TRUE, col = c("lightblue", "mistyrose",
"lightcyan"), ylim = c(0, 1.05), xlab = "Number of components",
ylab = "Performance")
if (choice == "Q2") {
graphics::text((best.num - 1) * 3 + best.num + 2.5, res[3, best.num] +
0.02, labels = "*", cex = 2.5, col = "red")
}
else if (choice == "R2") {
graphics::text((best.num - 1) * 3 + best.num + 1.5, res[2, best.num] +
0.02, labels = "*", cex = 2.5, col = "red")
}
else {
graphics::text((best.num - 1) * 3 + best.num + 0.5, res[1, best.num] +
0.02, labels = "*", cex = 2.5, col = "red")
}
xpos <- ncol(res) * 3 + ncol(res) + 1
graphics::legend(xpos, 1, rownames(res), fill = c("lightblue",
"mistyrose", "lightcyan"), xpd = T)
if(title==TRUE){
title.text <- paste0(
"(",
stringr::str_replace(mSetObj[["dataSet"]][["rownorm.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["trans.method"]], "N/A", ""),
"/",
stringr::str_replace(mSetObj[["dataSet"]][["scale.method"]], "N/A", ""),
")"
) %>% str_replace_all(., "(?<![:alpha:])/", "") %>% # Remove all commas that are not preceded by a letter.
str_replace_all(., "/(?![:alpha:])", "") %>% # Remove all commas that are not followed by [space]letter.
str_replace_all(., "\\(\\)", "(no normalization/transformation)") # Remove '()'.
graphics::title(stringr::str_wrap(title.text, 55), line = 1, cex.main = 1.2)
}
}
if(export == TRUE){
if (format == "pdf") {
grDevices::pdf(
file = imgName,
width = w,
height = h,
bg = "white",
onefile = FALSE
)
}
else {
grDevices::png(imgName,
width = w,
height = h,
units = 'in',
res = if_else(is.null(dpi), if_else(format=="pdf", 72, 300), dpi)
)
}
p()
grDevices::dev.off()
mSetObj$imgSet$pls.crossvalidation <- imgName
}
p(); mSetObj$imgSet$pls.crossvalidation.plot <- grDevices::recordPlot()
grDevices::dev.off()
if(plot==TRUE){
p()
}
return(mSetObj)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.