#' Visualize Longitudinal Feature
#'
#' Visualize Longitudinal Feature
#'
#' @param df dataframe has the Count, Group, ID, Time
#' @param text feature name
#' @param group.levels The two level's name
#' @param unit time interval unit
#' @param col two color to be used for the two groups (eg., c("red", "blue")).
#' @param ylabel text to be shown on the y-axis of all generated figures (default: "Normalized Count")
#' @param prefix prefix to be used to create directory for the analysis results
#' @import ggplot2
#' @import grDevices
#' @import graphics
#' @references
#' Ahmed Metwally (ametwall@stanford.edu)
#' @examples
#' data(metalonda_test_data)
#' pfx = tempfile()
#' dir.create(file.path(pfx))
#' n.sample = 5
#' n.timepoints = 10
#' n.group = 2
#' Group = factor(c(rep(0, n.sample*n.timepoints), rep(1, n.sample*n.timepoints)))
#' Time = rep(rep(1:n.timepoints, times = n.sample), 2)
#' ID = factor(rep(1:(2*n.sample), each = n.timepoints))
#' points = seq(1, 10, length.out = 10)
#' aggregate.df = data.frame(Count = metalonda_test_data[1,], Time = Time, Group = Group, ID = ID)
#' visualizeFeature(df = aggregate.df, text = rownames(metalonda_test_data)[1],
#' group.levels = Group, prefix = pfx)
#' @export
visualizeFeature = function (df, text, group.levels, unit = "days", ylabel = "Normalized Count",
col = c("blue", "firebrick"), prefix = "Test")
{
cat("Visualizing Feature = ", text, "\n")
Count=0; Time=0; ID=0; Group=0 ## This line is just to pass the CRAN checks for the aes in ggplot2
p = ggplot(df, aes(Time, Count, colour = Group, group = interaction(Group, ID)))
p = p + geom_point(size = 1, alpha = 0.5) + geom_line(size = 1, alpha = 0.7) + theme_bw() +
ggtitle(paste("Feature = ", text, sep = "")) + labs(y = ylabel, x = sprintf("Time (%s)", unit)) +
scale_colour_manual(values = col, breaks = c("0", "1"),
labels = c(group.levels[1], group.levels[2])) +
theme(axis.text.x = element_text(colour="black", size=12, angle=0, hjust=0.5, vjust=0.5, face="bold"),
axis.text.y = element_text(colour="black", size=12, angle=0, hjust=0.5, vjust=0.5, face="bold"),
axis.title.x = element_text(colour="black", size=15, angle=0, hjust=.5, vjust=0.5, face="bold"),
axis.title.y = element_text(colour="black", size=15, angle=90, hjust=.5, vjust=.5, face="bold"),
legend.text=element_text(size=15, face="plain"), legend.title = element_blank(),
plot.title = element_text(hjust = 0.5)) +
theme(legend.position="top") + scale_x_continuous(breaks = waiver())
#print("Prefix = ", prefix)
#ggsave(filename=paste("Feature_", text, ".jpg", sep=""), dpi = 1200, height = 10, width = 15, units = 'cm')
ggsave(filename=paste(prefix, "/", "Feature_", text, ".jpg", sep=""), dpi = 1200, height = 10, width = 15, units = 'cm')
}
#' Visualize the feature trajectory with the fitted Splines
#'
#' Plot the longitudinal features along with the fitted splines
#'
#' @param df dataframe has the Count , Group, ID, Time
#' @param model the fitted model
#' @param method The fitting method (negative binomial, LOWESS)
#' @param group.levels The two level's name
#' @param text feature name
#' @param unit time unit used in the Time vector (hours, days, weeks, months, etc.)
#' @param col two color to be used for the two groups (eg., c("red", "blue")).
#' @param ylabel text to be shown on the y-axis of all generated figures (default: "Normalized Count")
#' @param prefix prefix to be used to create directory for the analysis results
#' @import ggplot2
#' @import grDevices
#' @import graphics
#' @references
#' Ahmed Metwally (ametwall@stanford.edu)
#' @export
visualizeFeatureSpline = function (df, model, method, text, group.levels, unit = "days", ylabel = "Normalized Count",
col = c("blue", "firebrick"), prefix = "Test")
{
cat("Visualizing Splines of Feature = ", text, "\n")
Count=0;Time=0;ID=0;Group=0;lnn=0 ## This line is just to pass the CRAN checks for the aes in ggplot2
dd.null = model$dd.null
dd.0 = model$dd.0
dd.1 = model$dd.1
ln = factor(c(rep("longdash", nrow(df)), rep("longdash", nrow(dd.0)), rep("longdash", nrow(dd.1))))
size = c(rep(1, nrow(df)), rep(1, nrow(dd.0)), rep(1, nrow(dd.1)))
dm = rbind(df[,c("Time", "Count", "Group", "ID")], dd.0, dd.1)
dm$lnn=ln
dm$sz= size
p = ggplot(dm, aes(Time, Count, colour = Group, group = interaction(Group, ID)))
p = p + theme_bw() + geom_point(size=1, alpha=0.5) + geom_line(aes(linetype=lnn), size=1, alpha=0.5) +
ggtitle(paste("Feature = ", text, sep = "")) + labs(y = ylabel, x = sprintf("Time (%s)", unit)) +
scale_colour_manual(values = c(col, col),
breaks = c("0", "1", "fit.0", "fit.1"),
labels = c(group.levels[1], group.levels[2], paste(group.levels[1], ".fit", sep=""), paste(group.levels[2], ".fit", sep="")))+
theme(axis.text.x = element_text(colour="black", size=12, angle=0, hjust=0.5, vjust=0.5, face="bold"),
axis.text.y = element_text(colour="black", size=12, angle=0, hjust=0.5, vjust=0.5, face="bold"),
axis.title.x = element_text(colour="black", size=15, angle=0, hjust=.5, vjust=0.5, face="bold"),
axis.title.y = element_text(colour="black", size=15, angle=90, hjust=.5, vjust=.5, face="bold"),
legend.text=element_text(size=15, face="plain"), legend.title = element_blank(),
plot.title = element_text(hjust = 0.5)) +
theme(legend.position="top") + scale_x_continuous(breaks = waiver()) + guides(linetype=FALSE, size =FALSE)
ggsave(filename=paste(prefix, "/", "Feature_", text, "_CurveFitting_", method, ".jpg", sep=""), dpi = 1200, height = 10, width = 15, units = 'cm')
}
#' Visualize significant time interval
#'
#' Visualize significant time interval
#'
#' @param aggregate.df Dataframe has the Count, Group, ID, Time
#' @param model.ss The fitted model
#' @param method Fitting method (negative binomial or LOWESS)
#' @param start Vector of the start points of the time intervals
#' @param end Vector of the end points of the time intervals
#' @param text Feature name
#' @param group.levels Level's name
#' @param unit time unit used in the Time vector (hours, days, weeks, months, etc.)
#' @param col two color to be used for the two groups (eg., c("red", "blue")).
#' @param ylabel text to be shown on the y-axis of all generated figures (default: "Normalized Count")
#' @param prefix prefix to be used to create directory for the analysis results
#' @import ggplot2
#' @import grDevices
#' @import graphics
#' @references
#' Ahmed Metwally (ametwall@stanford.edu)
#' @export
visualizeArea = function(aggregate.df, model.ss, method, start, end, text, group.levels, unit = "days",
ylabel = "Normalized Count", col = c("blue", "firebrick"), prefix = "Test")
{
cat("Visualizing Significant Intervals of Feature = ", text, "\n")
Time = 0 ## This line is just to pass the CRAN checks for the aes in ggplot2
sub.11 = list()
sub.10 = list()
xx = NULL
for(i in 1:length(start))
{
sub.11[[i]] = subset(model.ss$dd.1, Time >= start[i] & Time <= end[i])
sub.10[[i]] = subset(model.ss$dd.0, Time >= start[i] & Time <= end[i])
cmd = sprintf('geom_ribbon(data=sub.10[[%d]], aes(ymin = sub.11[[%d]]$Count, ymax = Count), colour= "grey3", fill="grey69",
alpha = 0.6)', i, i)
if (i != 1)
{
xx = paste(xx, cmd, sep = "+")
} else
{
xx = cmd
}
}
# ddNULL = model_ss$ddNULL
dd.0 = model.ss$dd.0
dd.1 = model.ss$dd.1
dm = rbind(dd.0, dd.1)
p1 = 'ggplot(dm, aes(Time, Count, colour = Group, group = interaction(Group, ID))) +
theme_bw() + geom_point(size = 1, alpha = 0.5) + geom_line(size = 1, alpha = 0.5) +
ggtitle(paste("Feature = ", text, sep = "")) + labs(y = ylabel, x = sprintf("Time (%s)", unit)) +
scale_colour_manual(values = col,
breaks = c("fit.0", "fit.1"),
labels = c(paste(group.levels[1], ".fit", sep = ""), paste(group.levels[2], ".fit", sep = ""))) +
theme(axis.text.x = element_text(colour = "black", size = 12, angle = 0, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.text.y = element_text(colour = "black", size = 12, angle = 0, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.title.x = element_text(colour = "black", size = 15, angle = 0, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.title.y = element_text(colour = "black", size = 15, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
legend.text = element_text(size = 15, face="plain"), legend.title = element_blank(),
plot.title = element_text(hjust = 0.5)) +
theme(legend.position = "top") + scale_x_continuous(breaks = waiver())'
p2 = xx
p3 = paste(p1, p2, sep="+")
p = eval(parse(text = p3))
ggsave(filename=paste(prefix, "/", "Feature_", text, "_SignificantInterval_", method, ".jpg", sep=""), dpi = 1200, height = 10, width = 15, units = 'cm')
}
#' Visualize all significant time intervals for all tested features
#'
#' Visualize all significant time intervals for all tested features
#'
#' @param interval.details Dataframe has infomation about significant interval (feature name, start, end, dominant, p-value)
#' @param prefix prefix for the output figure
#' @param unit time unit used in the Time vector (hours, days, weeks, months, etc.)
#' @param col two color to be used for the two groups (eg., c("red", "blue")).
#' @import ggplot2
#' @import grDevices
#' @import graphics
#' @references
#' Ahmed Metwally (ametwall@stanford.edu)
#' @export
visualizeTimeIntervals = function(interval.details, prefix = "Test", unit = "days",
col = c("blue", "firebrick"))
{
feature=0;dominant=0;ID=0;Group=0;lnn=0 ## This line is just to pass the CRAN checks for the aes in ggplot2
interval.details$dominant = as.factor(interval.details$dominant)
interval.details$pvalue = as.numeric((interval.details$pvalue))
interval.details = interval.details[order(interval.details$feature), ]
ggplot(interval.details, aes(ymin = start , ymax = end, x = feature, xend = feature)) +
geom_linerange(aes(color = dominant), size = 1) +
coord_flip() + scale_colour_manual(values = col) +
labs(x = "Feature", y = sprintf("Time (%s)", unit), colour="Dominant") +
theme(axis.text.x = element_text(colour = "black", size = 10, angle = 0, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.text.y = element_text(colour = "black", size = 8, angle = 0, vjust = 0.5, face = "bold"),
axis.title.x = element_text(colour = "black", size = 15, angle = 0, hjust = 0.5, vjust = 0.5, face = "bold"),
axis.title.y = element_text(colour = "black", size = 15, angle = 90, hjust = 0.5, vjust = 0.5, face = "bold"),
legend.text = element_text(size = 15, face = "plain")) +
theme(panel.grid.minor = element_blank(),
panel.grid.major.y = element_line(colour = "white", size = 6),
panel.grid.major.x = element_line(colour = "white",size = 0.75)) +
theme(legend.position="top", panel.border = element_rect(colour = "black", fill = NA, size = 2))
ggsave(filename = paste(prefix, "/", prefix, "_MetaLonDA_TimeIntervals.jpg", sep=""), dpi = 1200, height = 30, width = 20, units = 'cm')
}
#' Visualize Area Ratio (AR) empirical distribution
#'
#' Visualize Area Ratio (AR) empirical distribution for each time interval
#'
#' @param permuted Permutation of the permuted data
#' @param text Feature name
#' @param method fitting method
#' @import ggplot2
#' @import grDevices
#' @import graphics
#' @param prefix prefix to be used to create directory for the analysis results
#' @references
#' Ahmed Metwally (ametwall@stanford.edu)
#' @export
visualizeARHistogram = function(permuted, text, method, prefix = "Test"){
cat("Visualizing AR Distribution for Feature = ", text, "\n")
n = ncol(permuted)
r = ceiling(sqrt(n))
c = ceiling(sqrt(n))
xx = paste(prefix, "/", "Feature_", text, "_AR_distribution_", method, ".jpg", sep = "")
jpeg(filename = xx, res = 1200, height = r*5, width = c*5, units = 'cm')
par(mfrow=c(r,c))
for( i in 1:ncol(permuted)){
hist(permuted[,i], xlab = "AR Ratio", ylab = "Frequency",
breaks = 10, col = "yellow", border = "red",
main = paste("Interval # ", i, sep=""), xlim = c(0,1))
}
dev.off()
}
#' Visualize log2 fold-change and significance of each interval as volcano plot
#'
#' Visualize log2 fold-change and significance of each interval as volcano plot
#'
#' @param df Dataframe has a detailed summary about feature's significant intervals
#' @param text Feature name
#' @import ggplot2
#' @import grDevices
#' @import graphics
#' @param prefix prefix to be used to create directory for the analysis results
#' @references
#' Ahmed Metwally (ametwall@stanford.edu)
#' @export
visualizeVolcanoPlot = function(df, text, prefix = "Test"){
adjusted.pvalue_pseudo=0; Significance=0; log2FoldChange=0 ## This line is just to pass the CRAN checks
cat("Visualizing Volcano Plot of Feature = ", text, "\n")
# Highlight features that have an absolute log2 fold change > 1 and p-value < 0.05
df$adjusted.pvalue_pseudo = df$adjusted.pvalue
df$adjusted.pvalue_pseudo[which(df$adjusted.pvalue == 0)] = 0.00001
df$Significance <- "NS"
df$Significance[(abs(df$log2FoldChange) > 1)] <- "FC"
df$Significance[(df$adjusted.pvalue_pseudo<0.05)] <- "FDR"
df$Significance[(df$adjusted.pvalue_pseudo<0.05) & (abs(df$log2FoldChange)>1)] <- "FC_FDR"
table(df$Significance)
df$Significance <- factor(df$Significance, levels=c("NS", "FC", "FDR", "FC_FDR"))
ggplot(data=df, aes(x=log2FoldChange, y=-log10(adjusted.pvalue_pseudo), colour=Significance)) +
geom_point(alpha=0.4, size=1.75) + theme_bw() +
ggtitle(paste("Feature = ", text, sep = "")) +
xlab("log2 fold change") + ylab("-log10 p-value") +
theme(axis.text.x = element_text(colour="black", size=12, angle=0, hjust=0.5, vjust=0.5, face="bold"),
axis.text.y = element_text(colour="black", size=12, angle=0, hjust=0.5, vjust=0.5, face="bold"),
axis.title.x = element_text(colour="black", size=12, angle=0, hjust=.5, vjust=0.5, face="bold"),
axis.title.y = element_text(colour="black", size=12, angle=90, hjust=.5, vjust=.5, face="bold"),
legend.text=element_text(size=10, face="plain"), legend.title = element_blank(),
plot.title = element_text(hjust = 0.5), legend.position="bottom") +
scale_color_manual(values=c(NS="grey30", FC="forestgreen", FDR="royalblue", FC_FDR="red2"),
labels=c(NS="NS", FC=paste("Log2FC>|", 1, "|", sep=""), FDR=paste("FDR Q<", 0.05, sep=""),
FC_FDR=paste("FDR Q<", 0.05, " & Log2FC>|", 1, "|", sep=""))) +
#Tidy the text labels for a subset of genes
geom_text(data=subset(df, adjusted.pvalue_pseudo<0.05 & abs(log2FoldChange)>=1),
aes(label=rownames(subset(df, adjusted.pvalue_pseudo<0.05 & abs(log2FoldChange)>= 1))),
size=2.25,
#segment.color="black", #This and the next parameter spread out the labels and join them to their points by a line
#segment.size=0.01,
check_overlap=TRUE,
vjust=1.0)+
geom_vline(xintercept=c(-1,1), linetype="dotted") + geom_hline(yintercept=-log10(0.05), linetype="dotted") +
ggsave(filename = paste(prefix, "/", "Feature_", text, "_VolcanoPlot.jpg", sep=""), dpi = 1200, height = 12, width = 12, units = 'cm')
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.