Open suprarai opened 1 month ago
Hello Surya,
Thanks for your message! The plot_otc_curves
function in the OTC
package currently does not allow plot customizations. For now, you can try the following function plot_otc_curves_flexible
, which I adapted to rename the x and y axis labels (function parameters name_xlab
and name_ylab
). The function now outputs the ggplot2
object itself, which can be saved and customized more flexibly with ggplot2
functions. To use the function, copy it and run it in your R session:
plot_otc_curves_flexible <- function(otc_curves,
load_data = FALSE,
CB = TRUE,
alpha = 0.05,
upper_t = 250,
lower_t = 30,
name_xlab = "threshold in nm",
name_ylab = "Optimal Transport Colocalization"){
if (load_data) {
env <- new.env()
nm <- load(otc_curves, env)[1]
otc_curves <- env[[nm]]
}
otc_means <- aggregate(otc_curves$OTC, by = list(otc_curves$class,
otc_curves$t), mean)
colnames(otc_means) <- c("class", "t", "OTC")
if (CB) {
cb_name = paste("CB_", toString(alpha), "_", sep = "")
otc_means$quantile = 0
for (class in unique(otc_curves$class)) {
data <- otc_curves[otc_curves$class == class, ]
n <- length(unique(data$no_pic))
Ts <- unique(data$t)
p <- length(Ts)
df <- data.frame(matrix(ncol = p, nrow = n))
colnames(df) <- Ts
for (j in 1:n) {
df[j, ] <- data[data$no_pic == j, ]$OTC
}
df_center <- scale(df, scale = FALSE)
cov_matrix = matrix(0, p, p)
for (i in 1:nrow(df_center)) {
cov_matrix = cov_matrix + as.numeric(df_center[i,
]) %*% t(as.numeric(df_center[i, ]))
}
cov_matrix = 1/(n - 1) * cov_matrix
sample = MASS::mvrnorm(n = 1000, mu = rep(0, p),
Sigma = cov_matrix)
sample = abs(sample)
sups = apply(sample, 1, max)
alpha = 0.05
u_alpha = quantile(sups, probs = alpha/2)
otc_means[otc_means$class == class, ]$quantile <- u_alpha
}
otc_means$upper <- otc_means$OTC + otc_means$quantile/sqrt(n)
otc_means$lower <- otc_means$OTC - otc_means$quantile/sqrt(n)
otc_means$upper_thresh <- 1
otc_means$lower_thresh <- 0
otc_means$upper <- apply(otc_means[, c("upper", "upper_thresh")],
1, min)
otc_means$lower <- apply(otc_means[, c("lower", "lower_thresh")],
1, max)
otc_means <- subset(otc_means, select = -c(lower_thresh,
upper_thresh))
}
else {
cb_name = "_"
}
otc_means <- otc_means[otc_means$t >= lower_t & otc_means$t <=
upper_t, ]
p <- ggplot(data = otc_means, aes(x = t, y = OTC, color = class)) +
geom_line() +
geom_ribbon(aes(x = t, ymax = upper, ymin = lower, fill = class),
alpha = 0.3) +
scale_x_continuous(name = name_xlab) +
scale_y_continuous(name = name_ylab, lim = c(0, 1)) +
guides(color = guide_legend(title = ""),
linetype = guide_legend(title = ""),
fill = guide_legend(title = "")) +
theme(legend.key.height = unit(2, "lines"),
legend.key.width = unit(2, "lines"))
return(p)
}
For example, if you want to add a vertical line at 50 on the x-axis, you can do this as follows:
p <- plot_otc_curves_flexible(otc_curves = otc_curves,
name_ylab = "new_ylab",
name_xlab = "new_xlab")
p + ggplot2::geom_vline(xintercept = 50)
Let me know if that helps :)
Best, Julia
Dear Julia, Thank you very much for the code. It means a lot to me.
I found a easier way. I first saved csv file and then recreated 95% OTC confidence band using ggplot2 and stat_summary function. This way, I can plot multiple OTC bands for multiple objects (lets say A, B, C, and D) with another object (lets say X) in a same figure.
But, I have another issue. I tried stat_summary function to make 95% confidence band. However, the width of the band (upper and lower threshold) is "bigger" as compared to the one created using OTC::plot_otc_curves. I guess is you might have used different formula to calculate the upper and lower threshold for the confidence band.
Could you please help me so that I can recreate exact similar plot from csv file as well ?. The csv file only spits "t" and "OTC" values. Thanks in advance! Best, Surya
Hi, I want to change some of the parameters on the graph like changing x and y axis labels, adding the line on the X and Y axis.
I try to do it on this step but could not able to make it.
plot otc curves
OTC::plot_otc_curves(otc_curves = otc_curves, output_path ="../results", output_name = "demo")
Could you please help me if it is possible to do it ? Best, Surya