Nothing
R/1_plots.R
In tinyarray: Expression Data Analysis and Visualization
Defines functions
draw_boxplot
draw_venn
draw_volcano
draw_heatmap
draw_pca
Documented in
draw_boxplot
draw_heatmap
draw_pca
draw_venn
draw_volcano
##' draw PCA plots
##'
##' do PCA analysis and warning a PCA plot
##'
##' @param exp A numeric matrix
##' @param group_list A character or factor vector with one label per sample;
##' character input will be converted to factor and unused levels will be
##' dropped.
##' @param color color vector
##' @param addEllipses logical,add ellipses or not
##' @param style plot style,"default","ggplot2"and "3D"
##' @param color.label color legend label
##' @param title plot title
##' @param ... other parameters from \code{fviz_pca_ind}
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 scale_color_manual
##' @importFrom ggplot2 scale_fill_manual
##' @importFrom ggplot2 theme_classic
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 labs
##' @importFrom ggplot2 aes
##' @importFrom ggplot2 stat_ellipse
##' @return a PCA plot according to \code{exp} and grouped by \code{group_list}.
##' @author Xiaojie Sun
##' @export
##' @examplesIf requireNamespace("FactoMineR", quietly = TRUE) && requireNamespace("factoextra", quietly = TRUE)
##' draw_pca(t(iris[,1:4]),iris$Species)
##' #change color
##' draw_pca(t(iris[,1:4]),iris$Species,color = c("#E78AC3", "#A6D854", "#FFD92F"))
##'
##' @examplesIf requireNamespace("FactoMineR", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
##' draw_pca(t(iris[,1:4]),iris$Species,style = "ggplot2")
##'
##' @examplesIf requireNamespace("FactoMineR", quietly = TRUE) && requireNamespace("scatterplot3d", quietly = TRUE)
##' draw_pca(t(iris[,1:4]),iris$Species,style = "3D")
##' @seealso
##' \code{\link{draw_heatmap}};\code{\link{draw_volcano}};\code{\link{draw_venn}}
draw_pca <- function(exp,group_list,
color = c("#2874C5","#f87669","#e6b707","#868686","#92C5DE","#F4A582","#66C2A5","#FC8D62","#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#B3B3B3"),
addEllipses = TRUE,
style = "default",
color.label = "Group",
title = "",
...){
if(!requireNamespace("FactoMineR",quietly = TRUE)) {
stop("Package \"FactoMineR\" needed for this function to work.
Please install it by install.packages('FactoMineR')",call. = FALSE)
}
p1 <- all(apply(exp,2,is.numeric))
if(!p1) stop("exp must be a numeric matrix")
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
dat <- as.data.frame(t(exp))
dat.pca <- FactoMineR::PCA(dat, graph = FALSE)
col = color[1:length(levels(group_list))]
if(style == "default"){
if(!requireNamespace("factoextra",quietly = TRUE)) {
stop("Package \"factoextra\" needed for this function to work.
Please install it by install.packages('factoextra')",call. = FALSE)
}
factoextra::fviz_pca_ind(dat.pca,
geom.ind = "point",
col.ind = group_list,
addEllipses = addEllipses,
palette = col,
legend.title = "Groups",
title = title,
...)
}else if(style == "ggplot2"){
pdat = data.frame(dat.pca[["ind"]][["coord"]],
Group = group_list)
p = ggplot(pdat,aes(Dim.1,Dim.2))+
geom_point(aes(Dim.1,Dim.2,fill = Group),
shape = 21,color = "black")+
scale_color_manual(values = color[1:nlevels(group_list)])+
scale_fill_manual(values = color[1:nlevels(group_list)])+
theme_classic()+
theme(legend.position = "top")+
labs(color = color.label,fill = color.label,title = title)
if(addEllipses) p = p +
stat_ellipse(aes(color = Group,fill = Group),
geom = "polygon",
alpha = 0.3,
linetype = 2)
return(p)
}else if(style == "3D"){
if(!requireNamespace("scatterplot3d",quietly = TRUE)) {
stop("Package \"scatterplot3d\" needed for this function to work.
Please install it by install.packages('scatterplot3d')",call. = FALSE)
}
colors = color[as.numeric(group_list)]
pdat = data.frame(dat.pca[["ind"]][["coord"]],
Group = group_list)
scatterplot3d::scatterplot3d(pdat[,1:3],
color = "black",
pch = 21,
bg = colors,
main = title)
graphics::legend("bottom",col = "black",
legend = levels(group_list),
pt.bg = color[1:nlevels(group_list)], pch = 21,
inset = -0.2, xpd = TRUE,
horiz = TRUE)
}
}
if(getRversion() >= "2.15.1") utils::globalVariables(c(".","Dim.1","Dim.2","Group","pdat"))
##' draw a heatmap plot
##'
##' draw a heatmap plot for an expression matrix grouped by \code{group_list};
##' if \code{scale = TRUE}, the matrix is scaled before plotting.
##'
##' @inheritParams draw_pca
##' @param n A numeric matrix
##' @param scale_before deprecated parameter
##' @param n_cutoff numeric cutoff used after scaling; values outside
##' \code{[-n_cutoff, n_cutoff]} are truncated before plotting.
##' @param annotation_legend logical,show annotation legend or not
##' @param color color for heatmap
##' @param color_an color for column annotation
##' @param legend logical,show legend or not
##' @param show_rownames logical,show rownames or not
##' @param scale logical,scale the matrix or not
##' @param main the title of the plot
##' @param split_column split column by group_list
##' @param show_column_title show column title or not
##' @param ... other parameters from pheatmap
##' @return a heatmap plot according to \code{n} and grouped by
##' \code{group_list}.
##' @author Xiaojie Sun
##' @importFrom pheatmap pheatmap
##' @export
##' @examplesIf requireNamespace("ggplotify", quietly = TRUE) && requireNamespace("pheatmap", quietly = TRUE)
##' #example data
##' exp = matrix(abs(rnorm(60,sd = 16)),nrow = 10)
##' exp[,4:6] <- exp[,4:6]+20
##' colnames(exp) <- paste0("sample",1:6)
##' rownames(exp) <- paste0("gene",1:10)
##' exp[1:4,1:4]
##' group_list = factor(rep(c("A","B"),each = 3))
##' draw_heatmap(exp,group_list)
##' #use iris
##' n = t(iris[,1:4]);colnames(n) = 1:150
##' group_list = iris$Species
##' draw_heatmap(n,group_list)
##' draw_heatmap(n,group_list,color = colorRampPalette(c("green","black","red"))(100),
##' color_an = c("red","blue","pink") )
##' @seealso
##' \code{\link{draw_pca}};\code{\link{draw_volcano}};\code{\link{draw_venn}}
draw_heatmap <- function(n,
group_list,
scale_before = FALSE,
n_cutoff = 3,
legend = FALSE,
show_rownames = FALSE,
annotation_legend=FALSE,
split_column = FALSE,
show_column_title = FALSE,
color = grDevices::colorRampPalette(c("#2fa1dd", "white", "#f87669"))(100),
color_an = c("#2fa1dd","#f87669","#e6b707","#868686","#92C5DE","#F4A582","#66C2A5","#FC8D62","#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#B3B3B3"),
scale = TRUE,
main = NA,...){
dots <- list(...)
cluster_cols <- TRUE
if ("cluster_cols" %in% names(dots)) {
cluster_cols <- isTRUE(dots$cluster_cols)
dots$cluster_cols <- NULL
}
if(!requireNamespace("ggplotify",quietly = TRUE)) {
stop("Package 'ggplotify' needed for this function to work.
Please install it by install.packages('ggplotify')",call. = FALSE)
}
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
if (!cluster_cols) {
reordered <- .tinyarray_order_columns_by_group(n, group_list)
n <- reordered$exp
group_list <- reordered$group_list
}
n = as.data.frame(n)
if(scale_before) {
message("scale_before parameter is deprecated")
scale_before = FALSE
}
p1 <- all(apply(n,2,is.numeric))
if(!p1) stop("non-numeric matrix detected")
heatmap_cluster <- .tinyarray_safe_heatmap_clusters(
nrow_value = nrow(n),
ncol_value = ncol(n),
cluster_rows = TRUE,
cluster_cols = cluster_cols
)
if(split_column){
if(!requireNamespace("circlize",quietly = TRUE)) {
stop("Package 'circlize' needed for this function to work.
Please install it by install.packages('circlize')",call. = FALSE)
}
if(!requireNamespace("ComplexHeatmap",quietly = TRUE)) {
stop("Package 'ComplexHeatmap' needed for this function to work.
Please install it by BiocManager::install('ComplexHeatmap')",call. = FALSE)
}
col_fun = circlize::colorRamp2(c(-n_cutoff, 0, n_cutoff), c(color[1], "white", color[100]))
if(scale){
n = t(scale(t(n)))
n[n > n_cutoff] = n_cutoff
n[n < - n_cutoff] = -n_cutoff
}
m0 = ComplexHeatmap::Heatmap(n,
column_split = group_list,
cluster_rows = heatmap_cluster$cluster_rows,
cluster_columns = heatmap_cluster$cluster_cols)
column_title = suppressWarnings(names(ComplexHeatmap::column_order(m0)))
top_annotation = ComplexHeatmap::HeatmapAnnotation(
cluster = ComplexHeatmap::anno_block(gp = grid::gpar(fill = (color_an[1:nlevels(group_list)])[match(column_title,levels(group_list))]),
labels = column_title,
labels_gp = grid::gpar(col = "white", fontsize = 12)))
if(show_column_title){
p = ComplexHeatmap::Heatmap(n,name = " ",
col = col_fun,
top_annotation = top_annotation,
column_split = group_list,
cluster_rows = heatmap_cluster$cluster_rows,
cluster_columns = heatmap_cluster$cluster_cols,
show_heatmap_legend = legend ,
border = sum(dim(n))< 150,
rect_gp = grid::gpar(col = "grey", lwd = 1),
show_column_names = FALSE,
show_row_names = show_rownames)
}else{
p = ComplexHeatmap::Heatmap(n,name = " ",
col = col_fun,
top_annotation = top_annotation,
column_split = group_list,
cluster_rows = heatmap_cluster$cluster_rows,
cluster_columns = heatmap_cluster$cluster_cols,
show_heatmap_legend = legend ,
border = sum(dim(n))< 150,
rect_gp = grid::gpar(col = "grey", lwd = 1),
show_column_names = FALSE,
show_row_names = show_rownames ,
column_title = NULL)
}
}else{
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n)
col = color_an[1:length(levels(group_list))]
ann_colors = list(group = col)
if(is.null(names(ann_colors$group)))names(ann_colors$group)=levels(group_list)
if(scale) {
scale_row = "row"
breaks = seq(-n_cutoff,n_cutoff,length.out = length(color))
} else {
scale_row = "none"
breaks = NA
}
if(!scale_before){
heatmap_breaks <- if (scale_row == "row") {
breaks
} else {
.tinyarray_safe_heatmap_breaks(n, length(color))
}
p = do.call(
pheatmap,
c(
list(
n,
show_colnames = FALSE,
show_rownames = show_rownames,
scale = scale_row,
color = color,
annotation_col = annotation_col,
annotation_colors = ann_colors,
breaks = heatmap_breaks,
legend = legend,
silent = TRUE,
annotation_legend = annotation_legend,
main = main,
cluster_rows = heatmap_cluster$cluster_rows,
cluster_cols = heatmap_cluster$cluster_cols
),
dots
)
)
p = ggplotify::as.ggplot(p)
}
}
return(p)
}
##' draw a volcano plot
##'
##' warning a volcano plot for Differential analysis result in data.frame format.
##'
##' @param deg a data.frame created by Differential analysis
##' @param pvalue_cutoff Cutoff value of pvalue,0.05 by default.
##' @param logFC_cutoff Cutoff value of logFC,1 by default.
##' @param pkg a integer ,means which Differential analysis packages you used,we support three packages by now, 1,2,3,4 respectively means "DESeq2","edgeR","limma(voom)","limma"; pkg = 4 uses `log2FoldChange` as the default x-axis label.
##' @param adjust a logical value; if TRUE, use \code{padj} and plot
##' \code{-log10(padj)}; otherwise use \code{P.value} and plot
##' \code{-log10(P.value)}.
##' @param change_col column name for the up/down/status label. Default uses
##' `change` when it exists; otherwise the function falls back to the
##' existing automatic classification logic.
##' @param symmetry a logical value ,would you like to get your plot symmetrical
##' @param color color vector
##' @param lab label for x axis in volcano plot; when `NA`, pkg = 4 uses
##' `log2FoldChange`, and the other package codes keep their original labels
##' @param xlab.package whether to use the package name as the x axis name
##' @return a volcano plot according to \code{logFC} and either
##' \code{-log10(P.value)} or \code{-log10(padj)}, depending on \code{adjust}.
##' @author Xiaojie Sun
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 aes
##' @importFrom ggplot2 geom_point
##' @importFrom ggplot2 scale_color_manual
##' @importFrom ggplot2 geom_vline
##' @importFrom ggplot2 geom_hline
##' @importFrom ggplot2 theme_bw
##' @importFrom ggplot2 labs
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 element_text
##' @importFrom ggplot2 scale_x_continuous
##' @importFrom ggplot2 unit
##' @importFrom ggplot2 theme
##' @export
##' @examplesIf requireNamespace("ggplot2", quietly = TRUE)
##' head(deseq_data)
##' draw_volcano(deseq_data)
##' draw_volcano(deseq_data,pvalue_cutoff = 0.01,logFC_cutoff = 2)
##' draw_volcano(deseq_data,color = c("darkgreen", "darkgrey", "#B2182B"))
##' @seealso
##' \code{\link{draw_heatmap}};\code{\link{draw_pca}};\code{\link{draw_venn}}
draw_volcano <- function(deg,lab=NA,xlab.package = TRUE,
pvalue_cutoff = 0.05,logFC_cutoff= 1,
pkg = 1,adjust = TRUE,change_col = "change",
symmetry = FALSE,
color = c("#2874C5", "grey","#f87669")){
if(!is.data.frame(deg)) stop("deg must be a data.frame created by Differential analysis")
if(pvalue_cutoff>0.1)warning("Your pvalue_cutoff seems too large")
if(pvalue_cutoff>=1)stop("pvalue_cutoff will never larger than 1")
pv_col = .tinyarray_pvalue_column(adjust)
if(!adjust){
dat = switch(EXPR = pkg,
v1 = deg[,c(2,5)],
v2 = deg[,c(1,4)],
v3 = deg[,c("logFC", "P.Value")],
v4 = deg[,c("logFC", "P.Value")])
}else{
dat = switch(EXPR = pkg,
v1 = deg[,c(2,6)],
v2 = deg[,c(1,5)],
v3 = deg[,c("logFC", "adj.P.Val")],
v4 = deg[,c("logFC", "adj.P.Val")])
}
colnames(dat)[1:2]=c("logFC",pv_col)
#logFC_cutoff <- with(dat,mean(abs(logFC)) + 2*sd(abs(logFC)) )
if (!is.null(change_col) && length(change_col) == 1L && !is.na(change_col) && nzchar(change_col)) {
if (change_col %in% colnames(deg)) {
dat$change = toupper(as.character(deg[[change_col]]))
} else if (identical(change_col, "change")) {
dat$change = with(dat, ifelse(logFC > logFC_cutoff & dat[[pv_col]] < pvalue_cutoff, "UP",
ifelse(logFC < -logFC_cutoff & dat[[pv_col]] < pvalue_cutoff, "DOWN", "NOT")))
} else {
stop("change_col `", change_col, "` was not found in deg", call. = FALSE)
}
} else if ("change" %in% colnames(deg)) {
dat$change = toupper(as.character(deg$change))
} else {
dat$change = with(dat,ifelse(logFC>logFC_cutoff & dat[[pv_col]]<pvalue_cutoff , "UP",
ifelse(logFC< -logFC_cutoff & dat[[pv_col]]<pvalue_cutoff , "DOWN","NOT")))
}
if(length(lab) == 1L && is.na(lab)) {
if (as.integer(pkg) == 4L) {
lab = "log2FoldChange"
} else if (xlab.package) {
lab = c("DESeq2","edgeR","limma(voom)","limma")[pkg]
} else {
lab = "logFC"
}
}
dat$neglogp <- -log10(pmax(dat[[pv_col]], .Machine$double.xmin))
ylab <- if (adjust) {
"-log10 adjusted P value"
} else {
"-log10 P value"
}
this_tile <- paste0(nrow(dat[dat$change =='DOWN',]),
' Downstream, ',
nrow(dat[dat$change =='UP',]),
' Upstream'
)
p <- ggplot(data = dat,
aes(x = logFC,
y = neglogp)) +
geom_point(alpha=0.4, size=1.75,
aes(color=change)) +
scale_color_manual(values=color)+
geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),linetype=4,col="black",linewidth=0.8) +
geom_hline(yintercept = -log10(pvalue_cutoff),linetype=4,col="black",linewidth=0.8) +
theme_bw()+
labs(title = this_tile , x = lab, y = ylab)+
theme(plot.title = element_text(hjust = 0.5))
ce = ceiling(max(abs(dat$logFC)))
if(symmetry)p <- p+scale_x_continuous(limits = c(-ce,ce),expand = c(0,0))
return(p)
}
utils::globalVariables(c("logFC","P.value","padj","change","neglogp"))
##' draw a venn plot
##'
##' warning a venn plot for deg result created by three packages
##'
##' @inheritParams VennDiagram::venn.diagram
##' @inheritParams VennDiagram::draw.single.venn
##' @param x a list for plot
##' @param color color vector
##' @param reverse logical,reflect the three-set Venn diagram along its central
##' vertical axis of symmetry. Use in combination with rotation
##' to generate all possible set orders
##' @param ... other parameters from venn.diagram
##' @return a venn plot according to \code{x}, \code{y}, and \code{z} using the
##' list names as labels.
##' @author Xiaojie Sun
##' @export
##' @examplesIf requireNamespace("VennDiagram",quietly = TRUE) && requireNamespace("ggplotify",quietly = TRUE) && requireNamespace("cowplot",quietly = TRUE)
##' x = list(Deseq2=sample(1:100,30),edgeR = sample(1:100,30),limma = sample(1:100,30))
##' draw_venn(x,"test")
##' draw_venn(x,"test",color = c("darkgreen", "darkblue", "#B2182B"))
##' @seealso
##' \code{\link{draw_pca}};\code{\link{draw_volcano}};\code{\link{draw_heatmap}}
draw_venn <- function(x,main,
color = c("#2874C5","#f87669","#e6b707","#868686","#66C2A5","#FC8D62","#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#B3B3B3"),
imagetype ="png",
filename=NULL,
lwd=1,
lty=1,
col=color[1:length(x)],
fill=color[1:length(x)],
cat.col=color[1:length(x)],
cat.cex = 1,
cat.dist = -0.15,
rotation.degree = 0,
main.cex = 1,
cex=1,
alpha = 0.1,
reverse=TRUE,
...){
if(as.numeric(grDevices::dev.cur())!=1) grDevices::graphics.off()
if(!requireNamespace("VennDiagram",quietly = TRUE)) {
stop("Package \"VennDiagram\" needed for this function to work.
Please install it by install.packages('VennDiagram')",call. = FALSE)
}
if(!requireNamespace("ggplotify",quietly = TRUE)) {
stop("Package \"ggplotify\" needed for this function to work.
Please install it by install.packages('ggplotify')",call. = FALSE)
}
if(!requireNamespace("cowplot",quietly = TRUE)) {
stop("Package \"cowplot\" needed for this function to work.
Please install it by install.packages('cowplot')",call. = FALSE)
}
if(!is.list(x)) stop("x must be a list")
if(length(x)>7) stop("why do you give me so many elements to compare, I reject!")
temp_device <- NULL
if (is.null(filename)) {
temp_device <- tempfile(fileext = ".png")
grDevices::png(temp_device)
on.exit({
if (grDevices::dev.cur() != 1) {
grDevices::dev.off()
}
if (!is.null(temp_device) && file.exists(temp_device)) {
unlink(temp_device)
}
}, add = TRUE)
}
p = VennDiagram::venn.diagram(x = x,
imagetype =imagetype,
filename=filename,
lwd=lwd,
lty=lty,
col=col,
fill=fill,
cat.col=fill,
cat.cex = cat.cex,
cat.dist = cat.dist,
rotation.degree = rotation.degree,
main = main,
main.cex = main.cex,
cex=cex,
alpha = alpha,
reverse=reverse,
...)
p = ggplotify::as.ggplot(cowplot::as_grob(p))
file.remove(dir(pattern = ("^VennDiagram.*log$")))
return(p)
}
##' draw boxplot for expression
##'
##' draw boxplot for expression
##'
##' @inheritParams draw_pca
##' @param method one of kruskal.test,aov,t.test and wilcox.test
##' @param width width of boxplot and error bar
##' @param sort whether the boxplot will be sorted
##' @param drop whether to discard insignificant values
##' @param pvalue_cutoff if drop = TRUE,genes with p-values below the threshold will be drawn
##' @param xlab title of the x axis
##' @param ylab title of the y axis
##' @param grouplab title of group legend
##' @param p.label whether to show the exact p value in the plot; if
##' `FALSE`, show significance stars
##' @param add_error_bar whether to add error bar
##' @param ... other parameters from stat_compare_means
##' @return a boxplot according to \code{exp} and grouped by \code{group}.
##' @author Xiaojie Sun
##' @importFrom tibble rownames_to_column
##' @importFrom ggplot2 stat_boxplot
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 geom_boxplot
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 aes
##' @importFrom ggplot2 labs
##' @importFrom ggplot2 scale_fill_manual
##' @importFrom ggplot2 element_text
##' @importFrom ggplot2 after_stat
##' @export
##' @examplesIf requireNamespace("tidyr", quietly = TRUE) && requireNamespace("ggpubr", quietly = TRUE)
##' draw_boxplot(t(iris[,1:4]),iris$Species)
##' exp <- matrix(rnorm(60),nrow = 10)
##' colnames(exp) <- paste0("sample",1:6)
##' rownames(exp) <- paste0("gene",1:10)
##' exp[,4:6] = exp[,4:6] +10
##' exp[1:4,1:4]
##' group_list <- factor(rep(c("A","B"),each = 3))
##' draw_boxplot(exp,group_list)
##' draw_boxplot(exp,group_list,color = c("grey","red"))
##' @seealso
##' \code{\link{draw_heatmap}};\code{\link{draw_volcano}};\code{\link{draw_venn}}
draw_boxplot = function(exp,group_list,
method = "kruskal.test",
sort = TRUE,
drop = FALSE,
width = 0.5,
pvalue_cutoff = 0.05,
xlab = "Gene",
ylab = "Expression",
grouplab = "Group",
p.label = FALSE,
add_error_bar = FALSE,
color = c("#2874C5","#f87669","#e6b707","#868686","#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854",
"#FFD92F", "#E5C494", "#B3B3B3"),
...){
if(!requireNamespace("ggpubr",quietly = TRUE)) {
stop("Package \"ggpubr\" needed for this function to work.
Please install it by install.packages('ggpubr')",call. = FALSE)
}
if(!requireNamespace("tidyr",quietly = TRUE)) {
stop("Package \"tidyr\" needed for this function to work.
Please install it by install.packages('tidyr')",call. = FALSE)
}
p0 <- all(apply(exp,2,is.numeric)) & (!is.null(rownames(exp)))
if(!p0) stop("exp must be a numeric matrix with rownames")
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
p1 = method %in% c("kruskal.test","aov","t.test","wilcox.test")
if(!p1) stop("method should be one of kruskal.test,aov,t.test and wilcox.test")
if(nlevels(group_list)>2 & method %in% c("t.test","wilcox.test")) stop("group_list parameter must have exactly 2 group")
if(method=="kruskal.test"){
x = apply(exp, 1, function(x){
stats::kruskal.test(x~group_list)$p.value
})
}else if(method=="aov"){
x = apply(exp, 1, function(x){
summary(stats::aov(x~group_list))[[1]]$`Pr(>F)`[1]
})
}else if(method=="t.test"){
x = apply(exp, 1, function(x){
stats::t.test(x~group_list)$p.value
})
}else if(method=="wilcox.test"){
x = apply(exp, 1, function(x){
stats::wilcox.test(x~group_list)$p.value
})
}
if(sum(x< pvalue_cutoff)==0) {
message("No rows below the threshold,plot all rows")
drop = FALSE
}
if(drop){
x = x[ x < pvalue_cutoff]
exp = exp[rownames(exp) %in% names(x),]
}
if(length(x)>40) message("it seems to be too many rows")
dat = rownames_to_column(as.data.frame(t(exp)),var = "sample")
dat$group = group_list
dat = tidyr::gather(dat,
rows,exp,-sample,-group)
if(sort){
dat$rows = factor(dat$rows,
levels = names(sort(x)),
ordered = TRUE)
}
col = color[1:length(levels(group_list))]
p = ggplot(dat,aes(rows,exp,fill = group))+
geom_boxplot( width = width)+
theme_bw()+
theme(legend.position = "top")+
labs(fill = grouplab,
x = xlab,
y = ylab)+
scale_fill_manual(values = col)
if(add_error_bar) p = p + stat_boxplot(geom ='errorbar', width = width)
if(!p.label){
p = p + ggpubr::stat_compare_means(aes(group = group,
label = ggplot2::after_stat(p.signif)),
method = method,...)
}else{
p = p + ggpubr::stat_compare_means(aes(group = group,
label = ggplot2::after_stat(p.format)),
method = method,...)
}
if(length(x)>10) p = p +
theme(axis.text.x = element_text(vjust = 1,hjust = 1,angle = 80))
return(p)
}
utils::globalVariables(c(".","rows","group","p.signif","p.format"))
##' ggheat
##'
##' draw heatmap plot with annotation by ggplot2
##'
##' @param dat expression matrix for plot
##' @param group group for expression colnames
##' @param cluster logical,cluster in both rows and column or not, default F,now replaced by cluster_rows and cluster_cols.
##' @param cluster_rows logical, if rows (on the plot) should be clustered, default F
##' @param cluster_cols logical, if column (on the plot) should be clustered, default F
##' @param show_rownames logical or NULL. When NULL, the function shows x-axis
##' labels for up to 100 rows, rotates them 45 degrees for 9 to 100 rows, and
##' hides them above 100 rows.
##' @param show_colnames logical,show colnames in plot or not, default T
##' @param rotate_rownames logical or NULL. When NULL, the function uses the
##' automatic rotation rule described in \code{show_rownames}.
##' @param groupname name of group legend
##' @param expname name of exp legend
##' @param fill_mid use median value as geom_tile fill midpoint
##' @param color color for heatmap
##' @param legend_color color for legend
##' @return a ggplot object
##' @author Xiaojie Sun
##' @importFrom pheatmap pheatmap
##' @importFrom tibble rownames_to_column
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 geom_tile
##' @importFrom ggplot2 scale_fill_manual
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 aes
##' @importFrom ggplot2 element_blank
##' @importFrom ggplot2 scale_x_continuous
##' @importFrom ggplot2 element_rect
##' @importFrom ggplot2 scale_fill_gradient2
##' @importFrom ggplot2 scale_x_discrete
##' @importFrom ggplot2 labs
##' @importFrom patchwork wrap_plots
##' @export
##' @examplesIf requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("pheatmap", quietly = TRUE) && requireNamespace("tidyr", quietly = TRUE) && requireNamespace("patchwork", quietly = TRUE)
##' exp_dat = matrix(sample(100:1000,40),ncol = 4)
##' exp_dat[1:(nrow(exp_dat)/2),] = exp_dat[1:(nrow(exp_dat)/2),]-1000
##' rownames(exp_dat) = paste0("sample",1:nrow(exp_dat))
##' colnames(exp_dat) = paste0("gene",1:ncol(exp_dat))
##' group = rep(c("A","B"),each = nrow(exp_dat)/2)
##' group = factor(group,levels = c("A","B"))
##' ggheat(exp_dat,group)
##' ggheat(exp_dat,group,cluster_rows = TRUE)
##' ggheat(exp_dat,group,cluster_rows = TRUE,show_rownames = FALSE,
##' show_colnames = FALSE,groupname = "risk",expname = "expression")
ggheat = function(dat,group,cluster = FALSE,
color = c("#2874C5", "white","#f87669"),
legend_color = c("#2874C5","#f87669","#e6b707","#868686","#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F",
"#E5C494", "#B3B3B3"),
show_rownames = NULL,show_colnames = TRUE,
cluster_rows = FALSE,cluster_cols = FALSE,
groupname = "group",expname = "exp",
fill_mid = TRUE,
rotate_rownames = NULL){
group <- .tinyarray_prepare_group_list_safe(group, arg = "group")
if (is.null(group)) return(NULL)
dat = data.frame(dat)
ph = pheatmap::pheatmap(
t(dat),
silent = TRUE,
breaks = .tinyarray_safe_heatmap_breaks(t(dat), 100L),
color = grDevices::colorRampPalette(c("#2fa1dd", "white", "#f87669"))(100),
cluster_rows = ncol(dat) > 1L,
cluster_cols = nrow(dat) > 1L
)
if(cluster){
cluster_rows = TRUE
cluster_cols = TRUE
}
if(cluster_rows && !is.null(ph$tree_row)){
dat = dat[,ph$tree_row$order]
}
if(cluster_cols && !is.null(ph$tree_col)){
dat = dat[ph$tree_col$order,]
group = group[ph$tree_col$order]
}
dat$group = group
dat2 = tidyr::gather(rownames_to_column(dat,
"samples"),
gene,exp,-group,-samples)
dat2$samples = factor(dat2$samples,levels = rownames(dat))
dat2$gene = factor(dat2$gene,levels = rev(colnames(dat)))
if(!cluster) dat2 = arrange(dat2,group)
ng = length(unique(group))
col = legend_color[1:ng]
names(col) = levels(group)
p1 = ggplot(dat, aes(x = 1:nrow(dat), y = 1)) +
geom_tile(aes(fill = group))+
scale_fill_manual(values = col)+
theme(panel.grid = element_blank(),
panel.background = element_blank(),
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank()) +
scale_x_continuous(expand = c(0,0))+
labs(fill = groupname)
midpoint = ifelse(fill_mid,stats::median(dat2$exp),0)
p2 = ggplot(dat2, aes(samples, gene, fill = exp)) +
geom_tile()+
theme(panel.grid = element_blank(),
panel.background = element_rect(fill = NA),
legend.background = element_rect(fill = NA),
plot.background = element_rect(fill = NA),
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank()
) +
scale_fill_gradient2(low = color[1],
mid = color[2],
high = color[3],
midpoint = midpoint)+
scale_x_discrete(expand = c(0,0))+
labs(fill = expname)
x_n <- nrow(dat)
if (is.null(show_rownames)) {
show_rownames <- x_n <= 100
}
if (is.null(rotate_rownames)) {
rotate_rownames <- x_n > 8 && x_n <= 100
}
if(!show_rownames) {
p2 = p2 + theme(axis.text.x = element_blank())
} else if (rotate_rownames) {
p2 = p2 + theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))
} else {
p2 = p2 + theme(axis.text.x = element_text(vjust = 0.5, hjust = 0.5, angle = 0))
}
if(!show_colnames) p2 = p2 + theme(axis.text.y = element_blank())
p = wrap_plots(p1,p2,nrow = 2,heights = c(1,11))
return(p)
}
utils::globalVariables(c("gene","samples"))
##' draw_tsne
##'
##' draw tsne plot with annotation by ggplot2
##'
##' @inheritParams draw_pca
##' @param perplexity numeric; perplexity parameter for Rtsne
##' @param color.label color legend label
##' @return a ggplot object
##' @author Xiaojie Sun
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 geom_point
##' @importFrom ggplot2 stat_ellipse
##' @importFrom ggplot2 theme_classic
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 aes
##' @export
##' @examplesIf requireNamespace("Rtsne", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
##' exp <- matrix(rnorm(10000),nrow = 50)
##' colnames(exp) <- paste0("sample",1:200)
##' rownames(exp) <- paste0("gene",1:50)
##' exp[1:4,1:4]
##' exp[,1:100] = exp[,1:100]+10
##' group_list <- factor(rep(c("A","B"),each = 100))
##' draw_tsne(exp,group_list)
draw_tsne = function(exp,group_list,perplexity=30,
color = c("#2874C5","#f87669","#e6b707","#868686","#92C5DE", "#F4A582", "#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3",
"#A6D854", "#FFD92F", "#E5C494", "#B3B3B3"),
color.label = "group",
addEllipses = TRUE){
if(!requireNamespace("Rtsne",quietly = TRUE))
stop("Package \"Rtsne\" needed for this function to work.
Please install it by install.packages('Rtsne')",call. = FALSE)
p1 <- all(apply(exp,2,is.numeric))
if(!p1) stop("exp must be a numeric matrix")
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
tsne_out = Rtsne::Rtsne(t(exp),perplexity = perplexity)
pdat = data.frame(tsne_out$Y,group_list)
colnames(pdat) = c("Y1","Y2","group")
p = ggplot(pdat,aes(Y1,Y2))+
geom_point(aes(Y1,Y2,fill = group),
shape = 21,color = "black")+
scale_color_manual(values = color[1:nlevels(group_list)])+
scale_fill_manual(values = color[1:nlevels(group_list)])+
theme_classic()+
theme(legend.position = "top")+
labs(color = color.label,fill = color.label)
if(addEllipses) p = p +
stat_ellipse(aes(color = group,fill = group),
geom = "polygon",
alpha = 0.3,
linetype = 2)
return(p)
}
utils::globalVariables(c("Y1","Y2","group","patient","ri","time","event"))
##' draw_KM
##'
##' draw KM-plot with two or more group
##'
##' @inheritParams draw_pca
##' @param meta survival data with time and event column
##' @param time_col colname of time
##' @param event_col colname of event
##' @param legend.title legend title
##' @param legend.labs character vector specifying legend labels
##' @param ... other parameters from ggsurvplot
##' @return a KM-plot
##' @author Xiaojie Sun
##' @importFrom survival survfit
##' @importFrom survival Surv
##' @importFrom survminer ggsurvplot
##' @export
##' @examplesIf requireNamespace("survival", quietly = TRUE) && requireNamespace("survminer", quietly = TRUE)
##' x = survival::lung
##'draw_KM(meta = x,
##' group_list = x$sex,event_col = "status")
draw_KM = function(meta,
group_list,
time_col = "time",event_col = "event",
legend.title = "Group",
legend.labs = NULL,
color = c("#2874C5","#f87669","#e6b707","#868686","#92C5DE", "#F4A582", "#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3",
"#A6D854", "#FFD92F", "#E5C494", "#B3B3B3"),
...){
p1 <- (time_col %in% colnames(meta)) & (event_col %in% colnames(meta))
if(!p1){
stop("meta data must involved time and event columns")
}else{
z1 = which(colnames(meta)==time_col)
z2 = which(colnames(meta)==event_col)
colnames(meta)[z1] = "time"
colnames(meta)[z2] = "event"
}
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
if(is.null(legend.labs)) legend.labs = levels(group_list)
meta$Group = group_list
sfit <- survival::survfit(Surv(time, event) ~ Group,
data = meta)
p = survminer::ggsurvplot(sfit,pval = T,data = meta,
palette = color[1:nlevels(group_list)],
legend.title = legend.title,
legend.labs = legend.labs,
...)
return(p$plot)
}
##' corscatterplot
##'
##' draw cor scatter plot with density plot by ggplot2
##' @param dat plot data
##' @param x x
##' @param y y
##' @param color_cor color for cor reg.line
##' @param fill_cor fill for cor reg.line
##' @param fill_x fill for top density plot
##' @param fill_y fill for right density plot
##' @param type whether to use a density plot or a histogram plot for the side panel.
##' @param ... other parameters for ggscatter
##' @return a ggplot object
##' @author Xiaojie Sun
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 geom_density
##' @importFrom ggplot2 theme_void
##' @importFrom ggplot2 theme_bw
##' @importFrom ggplot2 geom_histogram
##' @importFrom patchwork plot_layout
##' @importFrom patchwork plot_spacer
##' @export
##' @examplesIf requireNamespace("ggpubr", quietly = TRUE) && requireNamespace("patchwork", quietly = TRUE)
##' corscatterplot(iris,"Sepal.Length","Sepal.Width")
corscatterplot = function(dat,x,y,color_cor = "blue",fill_cor = "lightgray",
fill_x = "#ff820e",fill_y = "#0000fe",
type = "density",...){
if(!requireNamespace("ggpubr",quietly = TRUE))
stop("Package \"ggpubr\" needed for this function to work.
Please install it by install.packages('ggpubr')",call. = FALSE)
p1 <- ggpubr::ggscatter( dat, x = x, y = y,
add = "reg.line", conf.int = TRUE,
add.params = list(color = color_cor, fill = fill_cor),...)+
ggpubr::stat_cor()+
theme_bw()
if(type == "density"){
p2 <- ggplot(dat, do.call(ggplot2::aes, list(x = as.name(x)))) +
geom_density(fill = fill_x) +
theme_void()
p3 <- ggplot(dat, do.call(ggplot2::aes, list(x = as.name(y)))) +
geom_density(fill = fill_y) +
coord_flip() +
theme_void()
}else if( type == "histogram"){
p2 <- ggplot(dat, do.call(ggplot2::aes, list(x = as.name(x)))) +
geom_histogram(fill = fill_x,width = 0.1,color = "black") +
theme_void()
p3 <- ggplot(dat, do.call(ggplot2::aes, list(x = as.name(y)))) +
geom_histogram(fill = fill_y,width = 0.1,color = "black") +
coord_flip() +
theme_void()
}else{
stop("The `type` parameter must be either 'density' or 'histogram'.")
}
p2 + p1 + p3 + plot_spacer() + plot_layout(design = c("AAAAD
BBBBC
BBBBC
BBBBC
BBBBC"))
}
##' corheatmap
##'
##' draw cor heatmap
##' @param exp A numeric matrix
##' @param x genes or cells from exp
##' @param y genes or cells from exp
##' @param color color for heatmap
##' @return a ggplot object
##' @author Xiaojie Sun
##' @importFrom dplyr case_when
##' @importFrom pheatmap pheatmap
##' @importFrom Hmisc rcorr
##' @export
##' @examplesIf requireNamespace("Hmisc", quietly = TRUE) && requireNamespace("pheatmap", quietly = TRUE)
##' x = rownames(exprSet_hub1)[1:3]
##' y = rownames(exprSet_hub1)[4:7]
##' corheatmap(exprSet_hub1,x,y)
corheatmap = function(exp,x,y,color = c("#2fa1dd", "white", "#f87669")){
x = unique(x)
y = unique(y)
y = setdiff(y,x)
if (!length(x) || !length(y)) {
warning("corheatmap requires at least one gene in each set after overlap removal; skipping plot.", call. = FALSE)
return(NULL)
}
da = t(exp[c(x,y),])
rc = rcorr(da)
m = rc$r[seq_len(length(x)), (ncol(da)-length(y)+1):ncol(da), drop = FALSE]
p = rc$P[seq_len(length(x)), (ncol(da)-length(y)+1):ncol(da), drop = FALSE]
tmp = matrix(ifelse(p<=0.0001,"****",
ifelse(p<=0.001,"***",
ifelse(p<=0.01,"**",
ifelse(p<=0.05,"*","")))),
nrow = nrow(p),
ncol = ncol(p))
pheatmap(t(m),
display_numbers = t(tmp),
angle_col = ifelse(length(x)>10,45,0),
color = grDevices::colorRampPalette(color)(100),
breaks = .tinyarray_safe_heatmap_breaks(t(m), 100L),
border_color = "white",
width = 7,
height = 9.1,
treeheight_col = 0,
treeheight_row = 0,
cluster_rows = ncol(m) > 1L,
cluster_cols = nrow(m) > 1L)
}
Try the tinyarray package in your browser
Any scripts or data that you put into this service are public.
tinyarray documentation built on Aug. 2, 2026, 9:07 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 draw_boxplot draw_venn draw_volcano draw_heatmap draw_pca
Documented in draw_boxplot draw_heatmap draw_pca draw_venn draw_volcano
##' draw PCA plots
##'
##' do PCA analysis and warning a PCA plot
##'
##' @param exp A numeric matrix
##' @param group_list A character or factor vector with one label per sample;
##' character input will be converted to factor and unused levels will be
##' dropped.
##' @param color color vector
##' @param addEllipses logical,add ellipses or not
##' @param style plot style,"default","ggplot2"and "3D"
##' @param color.label color legend label
##' @param title plot title
##' @param ... other parameters from \code{fviz_pca_ind}
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 scale_color_manual
##' @importFrom ggplot2 scale_fill_manual
##' @importFrom ggplot2 theme_classic
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 labs
##' @importFrom ggplot2 aes
##' @importFrom ggplot2 stat_ellipse
##' @return a PCA plot according to \code{exp} and grouped by \code{group_list}.
##' @author Xiaojie Sun
##' @export
##' @examplesIf requireNamespace("FactoMineR", quietly = TRUE) && requireNamespace("factoextra", quietly = TRUE)
##' draw_pca(t(iris[,1:4]),iris$Species)
##' #change color
##' draw_pca(t(iris[,1:4]),iris$Species,color = c("#E78AC3", "#A6D854", "#FFD92F"))
##'
##' @examplesIf requireNamespace("FactoMineR", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
##' draw_pca(t(iris[,1:4]),iris$Species,style = "ggplot2")
##'
##' @examplesIf requireNamespace("FactoMineR", quietly = TRUE) && requireNamespace("scatterplot3d", quietly = TRUE)
##' draw_pca(t(iris[,1:4]),iris$Species,style = "3D")
##' @seealso
##' \code{\link{draw_heatmap}};\code{\link{draw_volcano}};\code{\link{draw_venn}}
draw_pca <- function(exp,group_list,
color = c("#2874C5","#f87669","#e6b707","#868686","#92C5DE","#F4A582","#66C2A5","#FC8D62","#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#B3B3B3"),
addEllipses = TRUE,
style = "default",
color.label = "Group",
title = "",
...){
if(!requireNamespace("FactoMineR",quietly = TRUE)) {
stop("Package \"FactoMineR\" needed for this function to work.
Please install it by install.packages('FactoMineR')",call. = FALSE)
}
p1 <- all(apply(exp,2,is.numeric))
if(!p1) stop("exp must be a numeric matrix")
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
dat <- as.data.frame(t(exp))
dat.pca <- FactoMineR::PCA(dat, graph = FALSE)
col = color[1:length(levels(group_list))]
if(style == "default"){
if(!requireNamespace("factoextra",quietly = TRUE)) {
stop("Package \"factoextra\" needed for this function to work.
Please install it by install.packages('factoextra')",call. = FALSE)
}
factoextra::fviz_pca_ind(dat.pca,
geom.ind = "point",
col.ind = group_list,
addEllipses = addEllipses,
palette = col,
legend.title = "Groups",
title = title,
...)
}else if(style == "ggplot2"){
pdat = data.frame(dat.pca[["ind"]][["coord"]],
Group = group_list)
p = ggplot(pdat,aes(Dim.1,Dim.2))+
geom_point(aes(Dim.1,Dim.2,fill = Group),
shape = 21,color = "black")+
scale_color_manual(values = color[1:nlevels(group_list)])+
scale_fill_manual(values = color[1:nlevels(group_list)])+
theme_classic()+
theme(legend.position = "top")+
labs(color = color.label,fill = color.label,title = title)
if(addEllipses) p = p +
stat_ellipse(aes(color = Group,fill = Group),
geom = "polygon",
alpha = 0.3,
linetype = 2)
return(p)
}else if(style == "3D"){
if(!requireNamespace("scatterplot3d",quietly = TRUE)) {
stop("Package \"scatterplot3d\" needed for this function to work.
Please install it by install.packages('scatterplot3d')",call. = FALSE)
}
colors = color[as.numeric(group_list)]
pdat = data.frame(dat.pca[["ind"]][["coord"]],
Group = group_list)
scatterplot3d::scatterplot3d(pdat[,1:3],
color = "black",
pch = 21,
bg = colors,
main = title)
graphics::legend("bottom",col = "black",
legend = levels(group_list),
pt.bg = color[1:nlevels(group_list)], pch = 21,
inset = -0.2, xpd = TRUE,
horiz = TRUE)
}
}
if(getRversion() >= "2.15.1") utils::globalVariables(c(".","Dim.1","Dim.2","Group","pdat"))
##' draw a heatmap plot
##'
##' draw a heatmap plot for an expression matrix grouped by \code{group_list};
##' if \code{scale = TRUE}, the matrix is scaled before plotting.
##'
##' @inheritParams draw_pca
##' @param n A numeric matrix
##' @param scale_before deprecated parameter
##' @param n_cutoff numeric cutoff used after scaling; values outside
##' \code{[-n_cutoff, n_cutoff]} are truncated before plotting.
##' @param annotation_legend logical,show annotation legend or not
##' @param color color for heatmap
##' @param color_an color for column annotation
##' @param legend logical,show legend or not
##' @param show_rownames logical,show rownames or not
##' @param scale logical,scale the matrix or not
##' @param main the title of the plot
##' @param split_column split column by group_list
##' @param show_column_title show column title or not
##' @param ... other parameters from pheatmap
##' @return a heatmap plot according to \code{n} and grouped by
##' \code{group_list}.
##' @author Xiaojie Sun
##' @importFrom pheatmap pheatmap
##' @export
##' @examplesIf requireNamespace("ggplotify", quietly = TRUE) && requireNamespace("pheatmap", quietly = TRUE)
##' #example data
##' exp = matrix(abs(rnorm(60,sd = 16)),nrow = 10)
##' exp[,4:6] <- exp[,4:6]+20
##' colnames(exp) <- paste0("sample",1:6)
##' rownames(exp) <- paste0("gene",1:10)
##' exp[1:4,1:4]
##' group_list = factor(rep(c("A","B"),each = 3))
##' draw_heatmap(exp,group_list)
##' #use iris
##' n = t(iris[,1:4]);colnames(n) = 1:150
##' group_list = iris$Species
##' draw_heatmap(n,group_list)
##' draw_heatmap(n,group_list,color = colorRampPalette(c("green","black","red"))(100),
##' color_an = c("red","blue","pink") )
##' @seealso
##' \code{\link{draw_pca}};\code{\link{draw_volcano}};\code{\link{draw_venn}}
draw_heatmap <- function(n,
group_list,
scale_before = FALSE,
n_cutoff = 3,
legend = FALSE,
show_rownames = FALSE,
annotation_legend=FALSE,
split_column = FALSE,
show_column_title = FALSE,
color = grDevices::colorRampPalette(c("#2fa1dd", "white", "#f87669"))(100),
color_an = c("#2fa1dd","#f87669","#e6b707","#868686","#92C5DE","#F4A582","#66C2A5","#FC8D62","#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#B3B3B3"),
scale = TRUE,
main = NA,...){
dots <- list(...)
cluster_cols <- TRUE
if ("cluster_cols" %in% names(dots)) {
cluster_cols <- isTRUE(dots$cluster_cols)
dots$cluster_cols <- NULL
}
if(!requireNamespace("ggplotify",quietly = TRUE)) {
stop("Package 'ggplotify' needed for this function to work.
Please install it by install.packages('ggplotify')",call. = FALSE)
}
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
if (!cluster_cols) {
reordered <- .tinyarray_order_columns_by_group(n, group_list)
n <- reordered$exp
group_list <- reordered$group_list
}
n = as.data.frame(n)
if(scale_before) {
message("scale_before parameter is deprecated")
scale_before = FALSE
}
p1 <- all(apply(n,2,is.numeric))
if(!p1) stop("non-numeric matrix detected")
heatmap_cluster <- .tinyarray_safe_heatmap_clusters(
nrow_value = nrow(n),
ncol_value = ncol(n),
cluster_rows = TRUE,
cluster_cols = cluster_cols
)
if(split_column){
if(!requireNamespace("circlize",quietly = TRUE)) {
stop("Package 'circlize' needed for this function to work.
Please install it by install.packages('circlize')",call. = FALSE)
}
if(!requireNamespace("ComplexHeatmap",quietly = TRUE)) {
stop("Package 'ComplexHeatmap' needed for this function to work.
Please install it by BiocManager::install('ComplexHeatmap')",call. = FALSE)
}
col_fun = circlize::colorRamp2(c(-n_cutoff, 0, n_cutoff), c(color[1], "white", color[100]))
if(scale){
n = t(scale(t(n)))
n[n > n_cutoff] = n_cutoff
n[n < - n_cutoff] = -n_cutoff
}
m0 = ComplexHeatmap::Heatmap(n,
column_split = group_list,
cluster_rows = heatmap_cluster$cluster_rows,
cluster_columns = heatmap_cluster$cluster_cols)
column_title = suppressWarnings(names(ComplexHeatmap::column_order(m0)))
top_annotation = ComplexHeatmap::HeatmapAnnotation(
cluster = ComplexHeatmap::anno_block(gp = grid::gpar(fill = (color_an[1:nlevels(group_list)])[match(column_title,levels(group_list))]),
labels = column_title,
labels_gp = grid::gpar(col = "white", fontsize = 12)))
if(show_column_title){
p = ComplexHeatmap::Heatmap(n,name = " ",
col = col_fun,
top_annotation = top_annotation,
column_split = group_list,
cluster_rows = heatmap_cluster$cluster_rows,
cluster_columns = heatmap_cluster$cluster_cols,
show_heatmap_legend = legend ,
border = sum(dim(n))< 150,
rect_gp = grid::gpar(col = "grey", lwd = 1),
show_column_names = FALSE,
show_row_names = show_rownames)
}else{
p = ComplexHeatmap::Heatmap(n,name = " ",
col = col_fun,
top_annotation = top_annotation,
column_split = group_list,
cluster_rows = heatmap_cluster$cluster_rows,
cluster_columns = heatmap_cluster$cluster_cols,
show_heatmap_legend = legend ,
border = sum(dim(n))< 150,
rect_gp = grid::gpar(col = "grey", lwd = 1),
show_column_names = FALSE,
show_row_names = show_rownames ,
column_title = NULL)
}
}else{
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n)
col = color_an[1:length(levels(group_list))]
ann_colors = list(group = col)
if(is.null(names(ann_colors$group)))names(ann_colors$group)=levels(group_list)
if(scale) {
scale_row = "row"
breaks = seq(-n_cutoff,n_cutoff,length.out = length(color))
} else {
scale_row = "none"
breaks = NA
}
if(!scale_before){
heatmap_breaks <- if (scale_row == "row") {
breaks
} else {
.tinyarray_safe_heatmap_breaks(n, length(color))
}
p = do.call(
pheatmap,
c(
list(
n,
show_colnames = FALSE,
show_rownames = show_rownames,
scale = scale_row,
color = color,
annotation_col = annotation_col,
annotation_colors = ann_colors,
breaks = heatmap_breaks,
legend = legend,
silent = TRUE,
annotation_legend = annotation_legend,
main = main,
cluster_rows = heatmap_cluster$cluster_rows,
cluster_cols = heatmap_cluster$cluster_cols
),
dots
)
)
p = ggplotify::as.ggplot(p)
}
}
return(p)
}
##' draw a volcano plot
##'
##' warning a volcano plot for Differential analysis result in data.frame format.
##'
##' @param deg a data.frame created by Differential analysis
##' @param pvalue_cutoff Cutoff value of pvalue,0.05 by default.
##' @param logFC_cutoff Cutoff value of logFC,1 by default.
##' @param pkg a integer ,means which Differential analysis packages you used,we support three packages by now, 1,2,3,4 respectively means "DESeq2","edgeR","limma(voom)","limma"; pkg = 4 uses `log2FoldChange` as the default x-axis label.
##' @param adjust a logical value; if TRUE, use \code{padj} and plot
##' \code{-log10(padj)}; otherwise use \code{P.value} and plot
##' \code{-log10(P.value)}.
##' @param change_col column name for the up/down/status label. Default uses
##' `change` when it exists; otherwise the function falls back to the
##' existing automatic classification logic.
##' @param symmetry a logical value ,would you like to get your plot symmetrical
##' @param color color vector
##' @param lab label for x axis in volcano plot; when `NA`, pkg = 4 uses
##' `log2FoldChange`, and the other package codes keep their original labels
##' @param xlab.package whether to use the package name as the x axis name
##' @return a volcano plot according to \code{logFC} and either
##' \code{-log10(P.value)} or \code{-log10(padj)}, depending on \code{adjust}.
##' @author Xiaojie Sun
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 aes
##' @importFrom ggplot2 geom_point
##' @importFrom ggplot2 scale_color_manual
##' @importFrom ggplot2 geom_vline
##' @importFrom ggplot2 geom_hline
##' @importFrom ggplot2 theme_bw
##' @importFrom ggplot2 labs
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 element_text
##' @importFrom ggplot2 scale_x_continuous
##' @importFrom ggplot2 unit
##' @importFrom ggplot2 theme
##' @export
##' @examplesIf requireNamespace("ggplot2", quietly = TRUE)
##' head(deseq_data)
##' draw_volcano(deseq_data)
##' draw_volcano(deseq_data,pvalue_cutoff = 0.01,logFC_cutoff = 2)
##' draw_volcano(deseq_data,color = c("darkgreen", "darkgrey", "#B2182B"))
##' @seealso
##' \code{\link{draw_heatmap}};\code{\link{draw_pca}};\code{\link{draw_venn}}
draw_volcano <- function(deg,lab=NA,xlab.package = TRUE,
pvalue_cutoff = 0.05,logFC_cutoff= 1,
pkg = 1,adjust = TRUE,change_col = "change",
symmetry = FALSE,
color = c("#2874C5", "grey","#f87669")){
if(!is.data.frame(deg)) stop("deg must be a data.frame created by Differential analysis")
if(pvalue_cutoff>0.1)warning("Your pvalue_cutoff seems too large")
if(pvalue_cutoff>=1)stop("pvalue_cutoff will never larger than 1")
pv_col = .tinyarray_pvalue_column(adjust)
if(!adjust){
dat = switch(EXPR = pkg,
v1 = deg[,c(2,5)],
v2 = deg[,c(1,4)],
v3 = deg[,c("logFC", "P.Value")],
v4 = deg[,c("logFC", "P.Value")])
}else{
dat = switch(EXPR = pkg,
v1 = deg[,c(2,6)],
v2 = deg[,c(1,5)],
v3 = deg[,c("logFC", "adj.P.Val")],
v4 = deg[,c("logFC", "adj.P.Val")])
}
colnames(dat)[1:2]=c("logFC",pv_col)
#logFC_cutoff <- with(dat,mean(abs(logFC)) + 2*sd(abs(logFC)) )
if (!is.null(change_col) && length(change_col) == 1L && !is.na(change_col) && nzchar(change_col)) {
if (change_col %in% colnames(deg)) {
dat$change = toupper(as.character(deg[[change_col]]))
} else if (identical(change_col, "change")) {
dat$change = with(dat, ifelse(logFC > logFC_cutoff & dat[[pv_col]] < pvalue_cutoff, "UP",
ifelse(logFC < -logFC_cutoff & dat[[pv_col]] < pvalue_cutoff, "DOWN", "NOT")))
} else {
stop("change_col `", change_col, "` was not found in deg", call. = FALSE)
}
} else if ("change" %in% colnames(deg)) {
dat$change = toupper(as.character(deg$change))
} else {
dat$change = with(dat,ifelse(logFC>logFC_cutoff & dat[[pv_col]]<pvalue_cutoff , "UP",
ifelse(logFC< -logFC_cutoff & dat[[pv_col]]<pvalue_cutoff , "DOWN","NOT")))
}
if(length(lab) == 1L && is.na(lab)) {
if (as.integer(pkg) == 4L) {
lab = "log2FoldChange"
} else if (xlab.package) {
lab = c("DESeq2","edgeR","limma(voom)","limma")[pkg]
} else {
lab = "logFC"
}
}
dat$neglogp <- -log10(pmax(dat[[pv_col]], .Machine$double.xmin))
ylab <- if (adjust) {
"-log10 adjusted P value"
} else {
"-log10 P value"
}
this_tile <- paste0(nrow(dat[dat$change =='DOWN',]),
' Downstream, ',
nrow(dat[dat$change =='UP',]),
' Upstream'
)
p <- ggplot(data = dat,
aes(x = logFC,
y = neglogp)) +
geom_point(alpha=0.4, size=1.75,
aes(color=change)) +
scale_color_manual(values=color)+
geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),linetype=4,col="black",linewidth=0.8) +
geom_hline(yintercept = -log10(pvalue_cutoff),linetype=4,col="black",linewidth=0.8) +
theme_bw()+
labs(title = this_tile , x = lab, y = ylab)+
theme(plot.title = element_text(hjust = 0.5))
ce = ceiling(max(abs(dat$logFC)))
if(symmetry)p <- p+scale_x_continuous(limits = c(-ce,ce),expand = c(0,0))
return(p)
}
utils::globalVariables(c("logFC","P.value","padj","change","neglogp"))
##' draw a venn plot
##'
##' warning a venn plot for deg result created by three packages
##'
##' @inheritParams VennDiagram::venn.diagram
##' @inheritParams VennDiagram::draw.single.venn
##' @param x a list for plot
##' @param color color vector
##' @param reverse logical,reflect the three-set Venn diagram along its central
##' vertical axis of symmetry. Use in combination with rotation
##' to generate all possible set orders
##' @param ... other parameters from venn.diagram
##' @return a venn plot according to \code{x}, \code{y}, and \code{z} using the
##' list names as labels.
##' @author Xiaojie Sun
##' @export
##' @examplesIf requireNamespace("VennDiagram",quietly = TRUE) && requireNamespace("ggplotify",quietly = TRUE) && requireNamespace("cowplot",quietly = TRUE)
##' x = list(Deseq2=sample(1:100,30),edgeR = sample(1:100,30),limma = sample(1:100,30))
##' draw_venn(x,"test")
##' draw_venn(x,"test",color = c("darkgreen", "darkblue", "#B2182B"))
##' @seealso
##' \code{\link{draw_pca}};\code{\link{draw_volcano}};\code{\link{draw_heatmap}}
draw_venn <- function(x,main,
color = c("#2874C5","#f87669","#e6b707","#868686","#66C2A5","#FC8D62","#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#B3B3B3"),
imagetype ="png",
filename=NULL,
lwd=1,
lty=1,
col=color[1:length(x)],
fill=color[1:length(x)],
cat.col=color[1:length(x)],
cat.cex = 1,
cat.dist = -0.15,
rotation.degree = 0,
main.cex = 1,
cex=1,
alpha = 0.1,
reverse=TRUE,
...){
if(as.numeric(grDevices::dev.cur())!=1) grDevices::graphics.off()
if(!requireNamespace("VennDiagram",quietly = TRUE)) {
stop("Package \"VennDiagram\" needed for this function to work.
Please install it by install.packages('VennDiagram')",call. = FALSE)
}
if(!requireNamespace("ggplotify",quietly = TRUE)) {
stop("Package \"ggplotify\" needed for this function to work.
Please install it by install.packages('ggplotify')",call. = FALSE)
}
if(!requireNamespace("cowplot",quietly = TRUE)) {
stop("Package \"cowplot\" needed for this function to work.
Please install it by install.packages('cowplot')",call. = FALSE)
}
if(!is.list(x)) stop("x must be a list")
if(length(x)>7) stop("why do you give me so many elements to compare, I reject!")
temp_device <- NULL
if (is.null(filename)) {
temp_device <- tempfile(fileext = ".png")
grDevices::png(temp_device)
on.exit({
if (grDevices::dev.cur() != 1) {
grDevices::dev.off()
}
if (!is.null(temp_device) && file.exists(temp_device)) {
unlink(temp_device)
}
}, add = TRUE)
}
p = VennDiagram::venn.diagram(x = x,
imagetype =imagetype,
filename=filename,
lwd=lwd,
lty=lty,
col=col,
fill=fill,
cat.col=fill,
cat.cex = cat.cex,
cat.dist = cat.dist,
rotation.degree = rotation.degree,
main = main,
main.cex = main.cex,
cex=cex,
alpha = alpha,
reverse=reverse,
...)
p = ggplotify::as.ggplot(cowplot::as_grob(p))
file.remove(dir(pattern = ("^VennDiagram.*log$")))
return(p)
}
##' draw boxplot for expression
##'
##' draw boxplot for expression
##'
##' @inheritParams draw_pca
##' @param method one of kruskal.test,aov,t.test and wilcox.test
##' @param width width of boxplot and error bar
##' @param sort whether the boxplot will be sorted
##' @param drop whether to discard insignificant values
##' @param pvalue_cutoff if drop = TRUE,genes with p-values below the threshold will be drawn
##' @param xlab title of the x axis
##' @param ylab title of the y axis
##' @param grouplab title of group legend
##' @param p.label whether to show the exact p value in the plot; if
##' `FALSE`, show significance stars
##' @param add_error_bar whether to add error bar
##' @param ... other parameters from stat_compare_means
##' @return a boxplot according to \code{exp} and grouped by \code{group}.
##' @author Xiaojie Sun
##' @importFrom tibble rownames_to_column
##' @importFrom ggplot2 stat_boxplot
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 geom_boxplot
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 aes
##' @importFrom ggplot2 labs
##' @importFrom ggplot2 scale_fill_manual
##' @importFrom ggplot2 element_text
##' @importFrom ggplot2 after_stat
##' @export
##' @examplesIf requireNamespace("tidyr", quietly = TRUE) && requireNamespace("ggpubr", quietly = TRUE)
##' draw_boxplot(t(iris[,1:4]),iris$Species)
##' exp <- matrix(rnorm(60),nrow = 10)
##' colnames(exp) <- paste0("sample",1:6)
##' rownames(exp) <- paste0("gene",1:10)
##' exp[,4:6] = exp[,4:6] +10
##' exp[1:4,1:4]
##' group_list <- factor(rep(c("A","B"),each = 3))
##' draw_boxplot(exp,group_list)
##' draw_boxplot(exp,group_list,color = c("grey","red"))
##' @seealso
##' \code{\link{draw_heatmap}};\code{\link{draw_volcano}};\code{\link{draw_venn}}
draw_boxplot = function(exp,group_list,
method = "kruskal.test",
sort = TRUE,
drop = FALSE,
width = 0.5,
pvalue_cutoff = 0.05,
xlab = "Gene",
ylab = "Expression",
grouplab = "Group",
p.label = FALSE,
add_error_bar = FALSE,
color = c("#2874C5","#f87669","#e6b707","#868686","#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854",
"#FFD92F", "#E5C494", "#B3B3B3"),
...){
if(!requireNamespace("ggpubr",quietly = TRUE)) {
stop("Package \"ggpubr\" needed for this function to work.
Please install it by install.packages('ggpubr')",call. = FALSE)
}
if(!requireNamespace("tidyr",quietly = TRUE)) {
stop("Package \"tidyr\" needed for this function to work.
Please install it by install.packages('tidyr')",call. = FALSE)
}
p0 <- all(apply(exp,2,is.numeric)) & (!is.null(rownames(exp)))
if(!p0) stop("exp must be a numeric matrix with rownames")
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
p1 = method %in% c("kruskal.test","aov","t.test","wilcox.test")
if(!p1) stop("method should be one of kruskal.test,aov,t.test and wilcox.test")
if(nlevels(group_list)>2 & method %in% c("t.test","wilcox.test")) stop("group_list parameter must have exactly 2 group")
if(method=="kruskal.test"){
x = apply(exp, 1, function(x){
stats::kruskal.test(x~group_list)$p.value
})
}else if(method=="aov"){
x = apply(exp, 1, function(x){
summary(stats::aov(x~group_list))[[1]]$`Pr(>F)`[1]
})
}else if(method=="t.test"){
x = apply(exp, 1, function(x){
stats::t.test(x~group_list)$p.value
})
}else if(method=="wilcox.test"){
x = apply(exp, 1, function(x){
stats::wilcox.test(x~group_list)$p.value
})
}
if(sum(x< pvalue_cutoff)==0) {
message("No rows below the threshold,plot all rows")
drop = FALSE
}
if(drop){
x = x[ x < pvalue_cutoff]
exp = exp[rownames(exp) %in% names(x),]
}
if(length(x)>40) message("it seems to be too many rows")
dat = rownames_to_column(as.data.frame(t(exp)),var = "sample")
dat$group = group_list
dat = tidyr::gather(dat,
rows,exp,-sample,-group)
if(sort){
dat$rows = factor(dat$rows,
levels = names(sort(x)),
ordered = TRUE)
}
col = color[1:length(levels(group_list))]
p = ggplot(dat,aes(rows,exp,fill = group))+
geom_boxplot( width = width)+
theme_bw()+
theme(legend.position = "top")+
labs(fill = grouplab,
x = xlab,
y = ylab)+
scale_fill_manual(values = col)
if(add_error_bar) p = p + stat_boxplot(geom ='errorbar', width = width)
if(!p.label){
p = p + ggpubr::stat_compare_means(aes(group = group,
label = ggplot2::after_stat(p.signif)),
method = method,...)
}else{
p = p + ggpubr::stat_compare_means(aes(group = group,
label = ggplot2::after_stat(p.format)),
method = method,...)
}
if(length(x)>10) p = p +
theme(axis.text.x = element_text(vjust = 1,hjust = 1,angle = 80))
return(p)
}
utils::globalVariables(c(".","rows","group","p.signif","p.format"))
##' ggheat
##'
##' draw heatmap plot with annotation by ggplot2
##'
##' @param dat expression matrix for plot
##' @param group group for expression colnames
##' @param cluster logical,cluster in both rows and column or not, default F,now replaced by cluster_rows and cluster_cols.
##' @param cluster_rows logical, if rows (on the plot) should be clustered, default F
##' @param cluster_cols logical, if column (on the plot) should be clustered, default F
##' @param show_rownames logical or NULL. When NULL, the function shows x-axis
##' labels for up to 100 rows, rotates them 45 degrees for 9 to 100 rows, and
##' hides them above 100 rows.
##' @param show_colnames logical,show colnames in plot or not, default T
##' @param rotate_rownames logical or NULL. When NULL, the function uses the
##' automatic rotation rule described in \code{show_rownames}.
##' @param groupname name of group legend
##' @param expname name of exp legend
##' @param fill_mid use median value as geom_tile fill midpoint
##' @param color color for heatmap
##' @param legend_color color for legend
##' @return a ggplot object
##' @author Xiaojie Sun
##' @importFrom pheatmap pheatmap
##' @importFrom tibble rownames_to_column
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 geom_tile
##' @importFrom ggplot2 scale_fill_manual
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 aes
##' @importFrom ggplot2 element_blank
##' @importFrom ggplot2 scale_x_continuous
##' @importFrom ggplot2 element_rect
##' @importFrom ggplot2 scale_fill_gradient2
##' @importFrom ggplot2 scale_x_discrete
##' @importFrom ggplot2 labs
##' @importFrom patchwork wrap_plots
##' @export
##' @examplesIf requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("pheatmap", quietly = TRUE) && requireNamespace("tidyr", quietly = TRUE) && requireNamespace("patchwork", quietly = TRUE)
##' exp_dat = matrix(sample(100:1000,40),ncol = 4)
##' exp_dat[1:(nrow(exp_dat)/2),] = exp_dat[1:(nrow(exp_dat)/2),]-1000
##' rownames(exp_dat) = paste0("sample",1:nrow(exp_dat))
##' colnames(exp_dat) = paste0("gene",1:ncol(exp_dat))
##' group = rep(c("A","B"),each = nrow(exp_dat)/2)
##' group = factor(group,levels = c("A","B"))
##' ggheat(exp_dat,group)
##' ggheat(exp_dat,group,cluster_rows = TRUE)
##' ggheat(exp_dat,group,cluster_rows = TRUE,show_rownames = FALSE,
##' show_colnames = FALSE,groupname = "risk",expname = "expression")
ggheat = function(dat,group,cluster = FALSE,
color = c("#2874C5", "white","#f87669"),
legend_color = c("#2874C5","#f87669","#e6b707","#868686","#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F",
"#E5C494", "#B3B3B3"),
show_rownames = NULL,show_colnames = TRUE,
cluster_rows = FALSE,cluster_cols = FALSE,
groupname = "group",expname = "exp",
fill_mid = TRUE,
rotate_rownames = NULL){
group <- .tinyarray_prepare_group_list_safe(group, arg = "group")
if (is.null(group)) return(NULL)
dat = data.frame(dat)
ph = pheatmap::pheatmap(
t(dat),
silent = TRUE,
breaks = .tinyarray_safe_heatmap_breaks(t(dat), 100L),
color = grDevices::colorRampPalette(c("#2fa1dd", "white", "#f87669"))(100),
cluster_rows = ncol(dat) > 1L,
cluster_cols = nrow(dat) > 1L
)
if(cluster){
cluster_rows = TRUE
cluster_cols = TRUE
}
if(cluster_rows && !is.null(ph$tree_row)){
dat = dat[,ph$tree_row$order]
}
if(cluster_cols && !is.null(ph$tree_col)){
dat = dat[ph$tree_col$order,]
group = group[ph$tree_col$order]
}
dat$group = group
dat2 = tidyr::gather(rownames_to_column(dat,
"samples"),
gene,exp,-group,-samples)
dat2$samples = factor(dat2$samples,levels = rownames(dat))
dat2$gene = factor(dat2$gene,levels = rev(colnames(dat)))
if(!cluster) dat2 = arrange(dat2,group)
ng = length(unique(group))
col = legend_color[1:ng]
names(col) = levels(group)
p1 = ggplot(dat, aes(x = 1:nrow(dat), y = 1)) +
geom_tile(aes(fill = group))+
scale_fill_manual(values = col)+
theme(panel.grid = element_blank(),
panel.background = element_blank(),
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank()) +
scale_x_continuous(expand = c(0,0))+
labs(fill = groupname)
midpoint = ifelse(fill_mid,stats::median(dat2$exp),0)
p2 = ggplot(dat2, aes(samples, gene, fill = exp)) +
geom_tile()+
theme(panel.grid = element_blank(),
panel.background = element_rect(fill = NA),
legend.background = element_rect(fill = NA),
plot.background = element_rect(fill = NA),
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank()
) +
scale_fill_gradient2(low = color[1],
mid = color[2],
high = color[3],
midpoint = midpoint)+
scale_x_discrete(expand = c(0,0))+
labs(fill = expname)
x_n <- nrow(dat)
if (is.null(show_rownames)) {
show_rownames <- x_n <= 100
}
if (is.null(rotate_rownames)) {
rotate_rownames <- x_n > 8 && x_n <= 100
}
if(!show_rownames) {
p2 = p2 + theme(axis.text.x = element_blank())
} else if (rotate_rownames) {
p2 = p2 + theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))
} else {
p2 = p2 + theme(axis.text.x = element_text(vjust = 0.5, hjust = 0.5, angle = 0))
}
if(!show_colnames) p2 = p2 + theme(axis.text.y = element_blank())
p = wrap_plots(p1,p2,nrow = 2,heights = c(1,11))
return(p)
}
utils::globalVariables(c("gene","samples"))
##' draw_tsne
##'
##' draw tsne plot with annotation by ggplot2
##'
##' @inheritParams draw_pca
##' @param perplexity numeric; perplexity parameter for Rtsne
##' @param color.label color legend label
##' @return a ggplot object
##' @author Xiaojie Sun
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 geom_point
##' @importFrom ggplot2 stat_ellipse
##' @importFrom ggplot2 theme_classic
##' @importFrom ggplot2 theme
##' @importFrom ggplot2 aes
##' @export
##' @examplesIf requireNamespace("Rtsne", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)
##' exp <- matrix(rnorm(10000),nrow = 50)
##' colnames(exp) <- paste0("sample",1:200)
##' rownames(exp) <- paste0("gene",1:50)
##' exp[1:4,1:4]
##' exp[,1:100] = exp[,1:100]+10
##' group_list <- factor(rep(c("A","B"),each = 100))
##' draw_tsne(exp,group_list)
draw_tsne = function(exp,group_list,perplexity=30,
color = c("#2874C5","#f87669","#e6b707","#868686","#92C5DE", "#F4A582", "#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3",
"#A6D854", "#FFD92F", "#E5C494", "#B3B3B3"),
color.label = "group",
addEllipses = TRUE){
if(!requireNamespace("Rtsne",quietly = TRUE))
stop("Package \"Rtsne\" needed for this function to work.
Please install it by install.packages('Rtsne')",call. = FALSE)
p1 <- all(apply(exp,2,is.numeric))
if(!p1) stop("exp must be a numeric matrix")
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
tsne_out = Rtsne::Rtsne(t(exp),perplexity = perplexity)
pdat = data.frame(tsne_out$Y,group_list)
colnames(pdat) = c("Y1","Y2","group")
p = ggplot(pdat,aes(Y1,Y2))+
geom_point(aes(Y1,Y2,fill = group),
shape = 21,color = "black")+
scale_color_manual(values = color[1:nlevels(group_list)])+
scale_fill_manual(values = color[1:nlevels(group_list)])+
theme_classic()+
theme(legend.position = "top")+
labs(color = color.label,fill = color.label)
if(addEllipses) p = p +
stat_ellipse(aes(color = group,fill = group),
geom = "polygon",
alpha = 0.3,
linetype = 2)
return(p)
}
utils::globalVariables(c("Y1","Y2","group","patient","ri","time","event"))
##' draw_KM
##'
##' draw KM-plot with two or more group
##'
##' @inheritParams draw_pca
##' @param meta survival data with time and event column
##' @param time_col colname of time
##' @param event_col colname of event
##' @param legend.title legend title
##' @param legend.labs character vector specifying legend labels
##' @param ... other parameters from ggsurvplot
##' @return a KM-plot
##' @author Xiaojie Sun
##' @importFrom survival survfit
##' @importFrom survival Surv
##' @importFrom survminer ggsurvplot
##' @export
##' @examplesIf requireNamespace("survival", quietly = TRUE) && requireNamespace("survminer", quietly = TRUE)
##' x = survival::lung
##'draw_KM(meta = x,
##' group_list = x$sex,event_col = "status")
draw_KM = function(meta,
group_list,
time_col = "time",event_col = "event",
legend.title = "Group",
legend.labs = NULL,
color = c("#2874C5","#f87669","#e6b707","#868686","#92C5DE", "#F4A582", "#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3",
"#A6D854", "#FFD92F", "#E5C494", "#B3B3B3"),
...){
p1 <- (time_col %in% colnames(meta)) & (event_col %in% colnames(meta))
if(!p1){
stop("meta data must involved time and event columns")
}else{
z1 = which(colnames(meta)==time_col)
z2 = which(colnames(meta)==event_col)
colnames(meta)[z1] = "time"
colnames(meta)[z2] = "event"
}
group_list <- .tinyarray_prepare_group_list_safe(group_list)
if (is.null(group_list)) return(NULL)
if(is.null(legend.labs)) legend.labs = levels(group_list)
meta$Group = group_list
sfit <- survival::survfit(Surv(time, event) ~ Group,
data = meta)
p = survminer::ggsurvplot(sfit,pval = T,data = meta,
palette = color[1:nlevels(group_list)],
legend.title = legend.title,
legend.labs = legend.labs,
...)
return(p$plot)
}
##' corscatterplot
##'
##' draw cor scatter plot with density plot by ggplot2
##' @param dat plot data
##' @param x x
##' @param y y
##' @param color_cor color for cor reg.line
##' @param fill_cor fill for cor reg.line
##' @param fill_x fill for top density plot
##' @param fill_y fill for right density plot
##' @param type whether to use a density plot or a histogram plot for the side panel.
##' @param ... other parameters for ggscatter
##' @return a ggplot object
##' @author Xiaojie Sun
##' @importFrom ggplot2 ggplot
##' @importFrom ggplot2 geom_density
##' @importFrom ggplot2 theme_void
##' @importFrom ggplot2 theme_bw
##' @importFrom ggplot2 geom_histogram
##' @importFrom patchwork plot_layout
##' @importFrom patchwork plot_spacer
##' @export
##' @examplesIf requireNamespace("ggpubr", quietly = TRUE) && requireNamespace("patchwork", quietly = TRUE)
##' corscatterplot(iris,"Sepal.Length","Sepal.Width")
corscatterplot = function(dat,x,y,color_cor = "blue",fill_cor = "lightgray",
fill_x = "#ff820e",fill_y = "#0000fe",
type = "density",...){
if(!requireNamespace("ggpubr",quietly = TRUE))
stop("Package \"ggpubr\" needed for this function to work.
Please install it by install.packages('ggpubr')",call. = FALSE)
p1 <- ggpubr::ggscatter( dat, x = x, y = y,
add = "reg.line", conf.int = TRUE,
add.params = list(color = color_cor, fill = fill_cor),...)+
ggpubr::stat_cor()+
theme_bw()
if(type == "density"){
p2 <- ggplot(dat, do.call(ggplot2::aes, list(x = as.name(x)))) +
geom_density(fill = fill_x) +
theme_void()
p3 <- ggplot(dat, do.call(ggplot2::aes, list(x = as.name(y)))) +
geom_density(fill = fill_y) +
coord_flip() +
theme_void()
}else if( type == "histogram"){
p2 <- ggplot(dat, do.call(ggplot2::aes, list(x = as.name(x)))) +
geom_histogram(fill = fill_x,width = 0.1,color = "black") +
theme_void()
p3 <- ggplot(dat, do.call(ggplot2::aes, list(x = as.name(y)))) +
geom_histogram(fill = fill_y,width = 0.1,color = "black") +
coord_flip() +
theme_void()
}else{
stop("The `type` parameter must be either 'density' or 'histogram'.")
}
p2 + p1 + p3 + plot_spacer() + plot_layout(design = c("AAAAD
BBBBC
BBBBC
BBBBC
BBBBC"))
}
##' corheatmap
##'
##' draw cor heatmap
##' @param exp A numeric matrix
##' @param x genes or cells from exp
##' @param y genes or cells from exp
##' @param color color for heatmap
##' @return a ggplot object
##' @author Xiaojie Sun
##' @importFrom dplyr case_when
##' @importFrom pheatmap pheatmap
##' @importFrom Hmisc rcorr
##' @export
##' @examplesIf requireNamespace("Hmisc", quietly = TRUE) && requireNamespace("pheatmap", quietly = TRUE)
##' x = rownames(exprSet_hub1)[1:3]
##' y = rownames(exprSet_hub1)[4:7]
##' corheatmap(exprSet_hub1,x,y)
corheatmap = function(exp,x,y,color = c("#2fa1dd", "white", "#f87669")){
x = unique(x)
y = unique(y)
y = setdiff(y,x)
if (!length(x) || !length(y)) {
warning("corheatmap requires at least one gene in each set after overlap removal; skipping plot.", call. = FALSE)
return(NULL)
}
da = t(exp[c(x,y),])
rc = rcorr(da)
m = rc$r[seq_len(length(x)), (ncol(da)-length(y)+1):ncol(da), drop = FALSE]
p = rc$P[seq_len(length(x)), (ncol(da)-length(y)+1):ncol(da), drop = FALSE]
tmp = matrix(ifelse(p<=0.0001,"****",
ifelse(p<=0.001,"***",
ifelse(p<=0.01,"**",
ifelse(p<=0.05,"*","")))),
nrow = nrow(p),
ncol = ncol(p))
pheatmap(t(m),
display_numbers = t(tmp),
angle_col = ifelse(length(x)>10,45,0),
color = grDevices::colorRampPalette(color)(100),
breaks = .tinyarray_safe_heatmap_breaks(t(m), 100L),
border_color = "white",
width = 7,
height = 9.1,
treeheight_col = 0,
treeheight_row = 0,
cluster_rows = ncol(m) > 1L,
cluster_cols = nrow(m) > 1L)
}
Try the tinyarray package in your browser
Any scripts or data that you put into this service are public.
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.