R/enrichment_analysis.R
In dputhier/dbfmcl: The scigenex package (Single-Cell Informative GENe Explorer)
#################################################################
## Running gprofiler2 to perform enrichment analysis
#################################################################
#' @title
#' Enrichment analysis using GO database of gene clusters from a ClusterSet object
##' @description
#' This function performs enrichment analysis of all gene clusters from a ClusterSet object
#' and store the results in the cluster_annotations slot.
#' @param object an object of class \code{ClusterSet}.
#' @param species a character string indicating the species name,
#' as a concatenation of the first letter of the name (uppercase) and the family name (lowercase),
#' e.g human -> Hsapiens, mouse -> Mmusculus
#' @param ontology a character string indicating the ontology to use for the enrichment analysis. One of "BP", "MF", and "CC" ontology, or "ALL".
#'
#' @return A \code{ClusterSet} object
#' @method enrich_go ClusterSet
#' @export enrich_go
#' @keywords internal
setGeneric("enrich_go",
function(object,
species="Hsapiens",
ontology="all") {
standardGeneric("enrich_go")
})
#' @title
#' Enrichment analysis using GO database of gene clusters from a ClusterSet object
##' @description
#' This function performs enrichment analysis of all gene clusters from a ClusterSet object
#' and store the results in the cluster_annotations slot.
#' @param object an object of class \code{ClusterSet}.
#' @param species a character string indicating the species name,
#' as a concatenation of the first letter of the name (uppercase) and the family name (lowercase),
#' e.g human -> Hsapiens, mouse -> Mmusculus
#' @param ontology a character string indicating the ontology to use for the enrichment analysis. One of "BP", "MF", and "CC" ontology, or "ALL".
#'
#' @return A \code{ClusterSet} object
#' @method enrich_go ClusterSet
#' @export enrich_go
#'
#' @examples
#' load_example_dataset("7871581/files/pbmc3k_medium_clusters")
#' load_example_dataset("8028126/files/pbmc3k_medium_clusters_enr")
#' # pbmc3k_medium_clusters_enr <- enrich_go(pbmc3k_medium_clusters)
#' pbmc3k_medium_clusters_enr
#' is(pbmc3k_medium_clusters_enr@gene_cluster_annotations)
setMethod("enrich_go",
signature(object = "ClusterSet"),
function(object,
species="Hsapiens",
ontology="all") {
## Check format object arg
check_format_cluster_set(object)
if (!species %in% c("Hsapiens", "Mmusculus")){
print_msg("This species name is not supported.",
msg_type = "STOP")}
if (species == "Hsapiens") {
org_db <- org.Hs.eg.db::org.Hs.eg.db
go_name <- "org.Hs.eg.db"
print_msg(msg_type = "INFO",
msg = "Species used : Homo sapiens")
}
if (species == "Mmusculus") {
org_db <- org.Mm.eg.db::org.Mm.eg.db
go_name <- "org.Mm.eg.db"
print_msg(msg_type = "INFO",
msg = "Species used : Mus musculus")
}
for(cluster in unique(names(object@gene_clusters))){
print_msg(paste0("Enrichment analysis for cluster ", cluster),
msg_type = "INFO")
cluster_name = paste0("Cluster_", cluster)
query = object@gene_clusters[[cluster]]
# Remove "MT-" or "mt-" not recognized in AnnotationDbi::select()
query = gsub("MT-", "", query)
query = gsub("mt-", "", query)
# Convert gene id to EntrezId format
suppressMessages(query_entrezid <- AnnotationDbi::select(org_db,
keys = query,
columns = c("ENTREZID", "SYMBOL"),
keytype = "SYMBOL"))
print_msg(msg_type = "DEBUG",
msg = paste0("Cluster ",
cluster,
"\n",
"\t|--Number of gene names converted to EntrezId : ",
length(which(!is.na(query_entrezid[, "ENTREZID"]))),
"\n",
"\t|--Number of gene names not converted to EntrezId : ",
length(which(is.na(query_entrezid[, "ENTREZID"]))),
"\n\t|-- Calling enrichGO."))
# Enrichment analysis using GO database
enrich_go_res <- clusterProfiler::enrichGO(query_entrezid[,"ENTREZID"],
OrgDb = go_name,
ont = ontology,
readable = TRUE)
# Store results in the ClusterSet object
object@gene_cluster_annotations[[cluster]] <- enrich_go_res
}
return(object)
}
)
############################################################################
## Visualization of enrichment analyses results from a ClusterSet object
############################################################################
#' @title
#' Visualization of enrichment analyses results from a ClusterSet object
##' @description
#' Display enrichment analyses results for a ClusterSet object.
#' @param object A \code{ClusterSet} object.
#' @param clusters A vector of cluster id to plot.
#' @param type A vector providing the type of plot.
#' @param nb_terms An integer indicating the number of terms in the plot.
#' @param terms_size An interger indicating the wrap length of terms
#' @param font_size The font size for the x axis.
#'
#' @return A \code{ClusterSet} object
#' @method viz_enrich ClusterSet
#' @export viz_enrich
#' @keywords internal
setGeneric("viz_enrich",
function(object,
clusters = "all",
type = "dotplot",
nb_terms = 20,
terms_size = 50,
font_size=4) {
standardGeneric("viz_enrich")
})
#' @title
#' Visualization of enrichment analyses results from a ClusterSet object
##' @description
#' Display enrichment analyses results for a ClusterSet object.
#' @param object A \code{ClusterSet} object.
#' @param clusters A vector of cluster id to plot.
#' @param type A vector providing the type of plot.
#' @param nb_terms An integer indicating the number of terms in the plot.
#' @param terms_size An interger indicating the wrap length of terms
#' @param font_size The font size for the x axis.
#'
#' @return A \code{ClusterSet} object
#' @method viz_enrich ClusterSet
#' @export viz_enrich
#'
#' @examples
#' load_example_dataset("8028226/files/pbmc3k_medium_clusters_enr_sub")
#' viz_enrich(pbmc3k_medium_clusters_enr_sub, cluster = 1, nb_terms = 5)
setMethod("viz_enrich",
signature(object = "ClusterSet"),
function(object,
clusters = "all",
type = c("dotplot", "barplot"),
nb_terms = 20,
terms_size = 50,
font_size=4) {
## Check format object arg
check_format_cluster_set(object)
if (length(clusters) == 1){
if (clusters == "all"){
clusters <- object@gene_clusters_metadata$cluster_id
}
}
list_plot <- list()
list_dotplot <- list()
list_barplot <- list()
for (cur_cluster in clusters) {
# Check if the current cluster id provided exists
if (!(cur_cluster %in% object@gene_clusters_metadata$cluster_id)) {
stop(paste0("Cluster ", cur_cluster, " doesn't exist."))
}
# Check if there is a result provided by enrich_go function for the current cluster
if(nrow(object@gene_cluster_annotations[[cur_cluster]]@result) == 0){
print_msg(msg_type = "WARNING",
msg = paste0("No functional enrichment analysis results for cluster ", cur_cluster, ".")) #Continue through the next cluster without plotting
} else {
# Print an informative message to announce plot of the results of the current cluster
print_msg(msg_type = "INFO",
msg = paste0("Plot enrichment analysis results for cluster ", cur_cluster))
# Create a ggplot - dotplot
if ("dotplot" %in% type){
if(object@gene_cluster_annotations[[cur_cluster]]@ontology == "GOALL"){
dot_plot <- enrichplot::dotplot(object@gene_cluster_annotations[[cur_cluster]],
split="ONTOLOGY",
showCategory=nb_terms,
label_format = terms_size)
dot_plot <- dot_plot + facet_grid(ONTOLOGY~., scales="free")
dot_plot <- dot_plot + ggtitle(paste0("Gene cluster ", cur_cluster))
dot_plot <- dot_plot + theme(axis.text = element_text(size=font_size))
} else {
dot_plot <- enrichplot::dotplot(object@gene_cluster_annotations[[cur_cluster]],
showCategory=nb_terms,
label_format = terms_size)
dot_plot <- dot_plot + ggtitle(paste0("Gene cluster ", cur_cluster))
dot_plot <- dot_plot + theme(axis.text = element_text(size=font_size))
}
list_dotplot <- append(list_dotplot, list(dot_plot))
}
# Create a ggplot - barplot
if ("barplot" %in% type){
if(object@gene_cluster_annotations[[cur_cluster]]@ontology == "GOALL"){
bar_plot <- graphics::barplot(object@gene_cluster_annotations[[cur_cluster]],
split="ONTOLOGY",
showCategory=nb_terms,
label_format = terms_size)
bar_plot <- bar_plot + facet_grid(ONTOLOGY~., scales="free")
bar_plot <- bar_plot + ggtitle(paste0("Gene cluster ", cur_cluster))
} else {
bar_plot <- graphics::barplot(object@gene_cluster_annotations[[cur_cluster]],
showCategory=nb_terms,
label_format = terms_size)
bar_plot <- bar_plot + ggtitle(paste0("Gene cluster ", cur_cluster))
}
list_barplot <- append(list_barplot, list(bar_plot))
}
list_plot <- append(list_barplot, list_dotplot)
}
}
return(invisible(lapply(list_plot, print)))
}
)
############################################################################
## plot_clust_enrichments()
############################################################################
#' @title Plot Cluster Enrichments
#'
#' @title Display cluster enrichments for gene GO terms.
#'
#' @param object An object of class "ClusterSet".
#' @param stat_shown A character string specifying the statistic to be shown. Must be one of "qvalue", "p.adjust", or "pvalue". Defaults to "qvalue".
#' @param nb_go_term An integer specifying the number of top GO terms to select per cluster. Defaults to 6. Note that some terms may be shared between clusters.
#' @param gradient_palette A vector of colors specifying the gradient palette for the color scale. Defaults to colors_for_gradient("Magma").
#' @param label_fun A function to customize the labels of the GO terms. Defaults to NULL.
#' @param term_order A vector specifying the desired order of the GO terms. Defaults to NULL.
#' @param return_mat Logical. Do not return a ggplot object but the prepared matrix.
#' @param floor A value for flooring the statistic.
#'
#' @return A ggplot object displaying the cluster enrichments for the GO terms.
#'
#' @seealso
#' \code{\link{enrich_go}}
#' \code{\link{colors_for_gradient}}
#'
#' @importFrom ggplot2 ggplot
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 geom_point
#' @importFrom ggplot2 scale_color_gradientn
#' @importFrom ggplot2 theme
#' @importFrom ggplot2 element_blank
#' @importFrom ggplot2 element_rect
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 xlab
#' @importFrom ggplot2 ylab
#' @importFrom ggplot2 guide_colourbar
#' @importFrom ggplot2 facet_grid
#' @importFrom reshape2 melt
#' @importFrom grid unit
#' @export plot_clust_enrichments
#' @examples
#' load_example_dataset("8028226/files/pbmc3k_medium_clusters_enr_sub")
#' plot_clust_enrichments(pbmc3k_medium_clusters_enr_sub)
#' @keywords internal
setGeneric("plot_clust_enrichments",
function(object,
stat_shown=c("qvalue", "p.adjust",
"pvalue"),
nb_go_term=3,
gradient_palette=colors_for_gradient("Ju1"),
label_fun=NULL,
term_order=NULL,
return_mat=FALSE,
floor=NULL) {
standardGeneric("plot_clust_enrichments")
})
#' @title Plot Cluster Enrichments
#'
#' @title Display cluster enrichments for gene GO terms.
#'
#' @param object An object of class "ClusterSet".
#' @param stat_shown A character string specifying the statistic to be shown. Must be one of "qvalue", "p.adjust", or "pvalue". Defaults to "qvalue".
#' @param nb_go_term An integer specifying the number of top GO terms to select per cluster. Defaults to 6. Note that some terms may be shared between clusters.
#' @param gradient_palette A vector of colors specifying the gradient palette for the color scale. Defaults to colors_for_gradient("Magma").
#' @param label_fun A function to customize the labels of the GO terms. Defaults to NULL.
#' @param term_order A vector specifying the desired order of the GO terms. Defaults to NULL.
#' @param return_mat Logical. Do not return a ggplot object but the prepared matrix.
#' @param floor A value for flooring the statistic.
#'
#' @return A ggplot object displaying the cluster enrichments for the GO terms.
#'
#' @seealso
#' \code{\link{enrich_go}}
#' \code{\link{colors_for_gradient}}
#'
#' @importFrom ggplot2 ggplot
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 geom_point
#' @importFrom ggplot2 scale_color_gradientn
#' @importFrom ggplot2 theme
#' @importFrom ggplot2 element_blank
#' @importFrom ggplot2 element_rect
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 xlab
#' @importFrom ggplot2 ylab
#' @importFrom ggplot2 guide_colourbar
#' @importFrom ggplot2 facet_grid
#' @importFrom reshape2 melt
#' @importFrom grid unit
#' @export plot_clust_enrichments
#' @examples
#' load_example_dataset("8028226/files/pbmc3k_medium_clusters_enr_sub")
#' plot_clust_enrichments(pbmc3k_medium_clusters_enr_sub)
setMethod("plot_clust_enrichments",
signature(object = "ClusterSet"),
function(object,
stat_shown=c("qvalue", "p.adjust",
"pvalue"),
nb_go_term=6,
gradient_palette=colors_for_gradient("Magma"),
label_fun=NULL,
term_order=NULL,
return_mat=FALSE,
floor=NULL) {
stat_shown <- match.arg(stat_shown)
print_msg(paste0("Using ", stat_shown, " as ordering statistic"), msg_type = "DEBUG")
print_msg("Checking object format", msg_type = "DEBUG")
check_format_cluster_set(object)
if(is.null(object@gene_cluster_annotations) | length(object@gene_cluster_annotations)==0){
print_msg("Please run enrich_go() first.", msg_type = "STOP")
}
print_msg("Extracting and formating data", msg_type = "DEBUG")
print_msg(paste0("Argument nb_go_term = ", nb_go_term), msg_type = "DEBUG")
result_go <- lapply(object@gene_cluster_annotations, "slot", "result")
order_mat <- function(x, column){x[order(x[,column]), ]}
result_go <- lapply(result_go, order_mat, stat_shown)
add_rank <- function(x, nb_go_term){x$rank <- 1:nrow(x);
x$rank[(nb_go_term + 1):nrow(x)] <- NA; return(x)}
result_go <- lapply(result_go, add_rank, nb_go_term)
for(i in 1:length(result_go)){
result_go[[i]]$cluster <- names(result_go)[i]
}
print_msg("Merging results", msg_type = "DEBUG")
m <- do.call('rbind', result_go)
go_retained <- unique(m$Description[!is.na(m$rank)])
m <- m[m$Description %in% go_retained, ]
if(!is.null(term_order)){
term_order_selected <- intersect(term_order, m$Description)
m <- m[m$Description %in% term_order_selected, ]
}
m$stat <- -log10(m[, stat_shown])
m <- m[order(m$cluster, m$rank), ]
m$is_top <- m$stat
m$is_top[is.na(m$rank)] <- "No"
m$is_top[!is.na(m$rank)] <- "Yes"
m$is_top <- as.factor(m$is_top)
if(return_mat){
return(m)
}
if(!is.null(label_fun)){
m$Description <- sapply(m$Description, label_fun)
}
if(!is.null(term_order)){
desc_levels <- sapply(term_order_selected, label_fun)
}else{
desc_levels <- rev(unique(m$Description))
}
m$Description <- factor(m$Description, levels=desc_levels, ordered = TRUE)
m$cluster <- factor(m$cluster,
levels=as.character(sort(as.numeric(unique(m$cluster)))),
ordered=TRUE)
cluster <- go_term <- stat <- size <- rank <- Counts <- NULL
print_msg("Melting...", msg_type = "DEBUG")
m_melt <- reshape2::melt(m, id.vars = c("Description", "Count", "stat", "cluster", "rank", "is_top"))
colnames(m_melt) <- c("go_term", "Counts", "stat", "cluster", "rank", "is_top", "ID", "value")
m_gg <- m_melt[order(m_melt[,"cluster"]), ]
print_msg("Running ggplot...", msg_type = "DEBUG")
if(!is.null(floor))
m_gg$stat[m_gg$stat > floor] <- floor
ggplot(data=m_gg, mapping=aes(x=cluster, y = go_term, color=stat, size=Counts)) +
ggplot2::theme_bw() +
geom_point(shape=16) +
ggplot2::scale_color_gradientn(colours = gradient_palette,
guide = guide_colourbar(title=paste0("-log10(", stat_shown, ")"), barwidth = 0.75, barheight = 5)) +
ggplot2::theme(
axis.ticks.y = ggplot2::element_blank(),
axis.ticks.x = ggplot2::element_blank(),
panel.spacing = grid::unit(0.0, "lines"),
panel.border = ggplot2::element_blank(),
strip.background.x = ggplot2::element_rect(fill = "#444444", colour = "white"),
strip.text.x = ggplot2::element_text(colour = "white", angle = 0),
axis.text.x = ggplot2::element_blank()
) +
ggplot2::xlab('Clusters') +
ggplot2::ylab('Go terms') +
ggplot2::facet_grid(~cluster, scale="free_x")
})
############################################################################
## plot_markers_to_clusters
############################################################################
#' @title Map cell markers to clusters (hypergeometric/jaccard)
#'
#' @description
#' Use jaccard index and hypergeometric to map a list of cell markers to gene clusters.
#'
#' @param object A ClusterSet object.
#' @param markers A list of cell markers.
#' @param gradient_palette A vector of colors for the gradient palette.
#' @param background The background for hypergeometric test. Typically the list of genes annotated in the considered gene ontology (e.g. BP).
#' @param cell_type_order The order of cell types. Can be "None" or "hclust". Not supported at the moment.
#' @param only_specific Logical value indicating whether to show only specific markers.
#' @param jaccard_scale The scale factor for Jaccard coefficient.
#'
#' @return A plot showing the relationship between markers and clusters.
#'
#' @export
#' @import ggplot2
#' @importFrom reshape2 melt
#' @importFrom grid unit
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 scale_color_gradientn
#' @importFrom ggplot2 element_blank
#' @importFrom ggplot2 element_line
#' @importFrom ggplot2 element_rect
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 xlab
#' @importFrom ggplot2 ylab
#' @importFrom ggplot2 facet_grid
#' @importFrom ggplot2 guide_colourbar
#' @importFrom ggplot2 geom_point
#'
#' @examples
#' load_example_dataset("8028126/files/pbmc3k_medium_clusters_enr")
#' library(clustermole)
#' m <- clustermole::clustermole_markers(species = "hs")
#' markers <- split(m$gene[m$organ == "Blood" & m$species == "Human"],
#' m$celltype[m$organ == "Blood" & m$species == "Human"])
#' plot_markers_to_clusters(pbmc3k_medium_clusters_enr, markers=markers)
#' @keywords internal
setGeneric("plot_markers_to_clusters",
function(object,
markers=NULL,
gradient_palette=colors_for_gradient("Magma"),
background=NULL,
cell_type_order=c("None", "hclust"),
only_specific=FALSE,
jaccard_scale=10) {
standardGeneric("plot_markers_to_clusters")
})
#' @title Map cell markers to clusters (hypergeometric/jaccard)
#'
#' @description
#' Use jaccard index and hypergeometric to map a list of cell markers to gene clusters.
#'
#' @param object A ClusterSet object.
#' @param markers A list of cell markers.
#' @param gradient_palette A vector of colors for the gradient palette.
#' @param background The background for hypergeometric test. Typically the list of genes annotated in the considered gene ontology (e.g. BP).
#' @param cell_type_order The order of cell types. Can be "None" or "hclust". Not supported at the moment.
#' @param only_specific Logical value indicating whether to show only specific markers.
#' @param jaccard_scale The scale factor for Jaccard coefficient.
#'
#' @return A plot showing the relationship between markers and clusters.
#'
#' @export
#' @import ggplot2
#' @importFrom reshape2 melt
#' @importFrom grid unit
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 scale_color_gradientn
#' @importFrom ggplot2 element_blank
#' @importFrom ggplot2 element_line
#' @importFrom ggplot2 element_rect
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 xlab
#' @importFrom ggplot2 ylab
#' @importFrom ggplot2 facet_grid
#' @importFrom ggplot2 guide_colourbar
#' @importFrom ggplot2 geom_point
#' @importFrom stats as.formula
#'
#' @examples
#' load_example_dataset("8028126/files/pbmc3k_medium_clusters_enr")
#' library(clustermole)
#' m <- clustermole_markers(species = "hs")
#' markers <- split(m$gene[m$organ == "Blood" & m$species == "Human"],
#' m$celltype[m$organ == "Blood" & m$species == "Human"])
#' plot_markers_to_clusters(pbmc3k_medium_clusters_enr, markers=markers)
setMethod("plot_markers_to_clusters",
signature(object = "ClusterSet"),
function(object,
markers=NULL,
gradient_palette=c("white", colors_for_gradient("De2")),
background=NULL,
cell_type_order=c("None", "hclust"),
only_specific=FALSE,
jaccard_scale=10) {
cell_type_order <- match.arg(cell_type_order)
if(is.null(background)){
print_msg("WARNING: Running without a background but you should probably provide one.", msg_type = "INFO")
}
print_msg("Checking object format", msg_type = "DEBUG")
check_format_cluster_set(object)
if(length(object@gene_clusters) == 0){
print_msg("No cluster found in the object.", msg_type = "STOP")
}
if(!is.list(markers))
print_msg("The markers argument should be a list.", msg_type = "STOP")
if(is.null(markers) | length(markers)==0){
print_msg("No marker found.", msg_type = "STOP")
}
if(only_specific){
for(i in 1:length(markers)){
cell_marker_curr <- markers[[i]]
all_other_marker <- unique(unlist(cell_marker_curr[(1:length(markers))[-i]]))
cell_marker_curr <- cell_marker_curr[!cell_marker_curr %in% all_other_marker]
markers[[i]] <- cell_marker_curr
}
nb_cell_markers <- lapply(markers, length)
markers <- markers[nb_cell_markers > 0]
}
comparison_jaccard <- compare_genesets(object@gene_clusters,
markers,
background,
stat = "jaccard")
comparison_jaccard_m <- reshape2::melt(comparison_jaccard)
comparison_hyper <- compare_genesets(object@gene_clusters,
markers,
background,
stat = "hypergeom")
comparison_hyper_m <- reshape2::melt(comparison_hyper)
comparisons_m <- comparison_jaccard_m
colnames(comparisons_m) <- c("cluster", "marker", "jaccard")
comparisons_m$pvalue <- comparison_hyper_m$value
comparisons_m$pvalue[comparisons_m$pvalue < 1e-320] <- 1e-320
comparisons_m$pvalue <- -log10(comparisons_m$pvalue)
comparisons_m$cluster <- factor(comparisons_m$cluster)
comparisons_m$marker <- factor(comparisons_m$marker)
m_gg <- comparisons_m[order(comparisons_m[,"cluster"],
comparisons_m[,"marker"],
comparisons_m[,"pvalue"]), ]
cluster <- marker <- pvalue <- jaccard <- NULL
print_msg("Running ggplot...", msg_type = "DEBUG")
ggplot(data=m_gg, mapping=aes(x=cluster, y = marker, color=pvalue, size=jaccard*jaccard_scale)) +
ggplot2::theme_bw() +
geom_point(shape=15) +
ggplot2::scale_color_gradientn(colours = gradient_palette,
guide = guide_colourbar(title=paste0("-log10(pvalue)"), barwidth = 0.75, barheight = 5)) +
ggplot2::theme(
axis.ticks.x = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_line(color="gray"),
panel.spacing = grid::unit(0.0, "lines"),
panel.border = element_rect(size=0.1, color="gray"),
panel.grid = element_blank(),
strip.background.x = ggplot2::element_rect(fill = "#444444", colour = "white"),
strip.text.x = ggplot2::element_text(colour = "white", angle = 0),
axis.text.x = ggplot2::element_blank()
) +
ggplot2::xlab('Clusters') +
ggplot2::ylab('Markers') +
ggplot2::facet_grid(~cluster, scale="free_x")
})
dputhier/dbfmcl documentation built on Dec. 24, 2024, 4:35 p.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.
#################################################################
## Running gprofiler2 to perform enrichment analysis
#################################################################
#' @title
#' Enrichment analysis using GO database of gene clusters from a ClusterSet object
##' @description
#' This function performs enrichment analysis of all gene clusters from a ClusterSet object
#' and store the results in the cluster_annotations slot.
#' @param object an object of class \code{ClusterSet}.
#' @param species a character string indicating the species name,
#' as a concatenation of the first letter of the name (uppercase) and the family name (lowercase),
#' e.g human -> Hsapiens, mouse -> Mmusculus
#' @param ontology a character string indicating the ontology to use for the enrichment analysis. One of "BP", "MF", and "CC" ontology, or "ALL".
#'
#' @return A \code{ClusterSet} object
#' @method enrich_go ClusterSet
#' @export enrich_go
#' @keywords internal
setGeneric("enrich_go",
function(object,
species="Hsapiens",
ontology="all") {
standardGeneric("enrich_go")
})
#' @title
#' Enrichment analysis using GO database of gene clusters from a ClusterSet object
##' @description
#' This function performs enrichment analysis of all gene clusters from a ClusterSet object
#' and store the results in the cluster_annotations slot.
#' @param object an object of class \code{ClusterSet}.
#' @param species a character string indicating the species name,
#' as a concatenation of the first letter of the name (uppercase) and the family name (lowercase),
#' e.g human -> Hsapiens, mouse -> Mmusculus
#' @param ontology a character string indicating the ontology to use for the enrichment analysis. One of "BP", "MF", and "CC" ontology, or "ALL".
#'
#' @return A \code{ClusterSet} object
#' @method enrich_go ClusterSet
#' @export enrich_go
#'
#' @examples
#' load_example_dataset("7871581/files/pbmc3k_medium_clusters")
#' load_example_dataset("8028126/files/pbmc3k_medium_clusters_enr")
#' # pbmc3k_medium_clusters_enr <- enrich_go(pbmc3k_medium_clusters)
#' pbmc3k_medium_clusters_enr
#' is(pbmc3k_medium_clusters_enr@gene_cluster_annotations)
setMethod("enrich_go",
signature(object = "ClusterSet"),
function(object,
species="Hsapiens",
ontology="all") {
## Check format object arg
check_format_cluster_set(object)
if (!species %in% c("Hsapiens", "Mmusculus")){
print_msg("This species name is not supported.",
msg_type = "STOP")}
if (species == "Hsapiens") {
org_db <- org.Hs.eg.db::org.Hs.eg.db
go_name <- "org.Hs.eg.db"
print_msg(msg_type = "INFO",
msg = "Species used : Homo sapiens")
}
if (species == "Mmusculus") {
org_db <- org.Mm.eg.db::org.Mm.eg.db
go_name <- "org.Mm.eg.db"
print_msg(msg_type = "INFO",
msg = "Species used : Mus musculus")
}
for(cluster in unique(names(object@gene_clusters))){
print_msg(paste0("Enrichment analysis for cluster ", cluster),
msg_type = "INFO")
cluster_name = paste0("Cluster_", cluster)
query = object@gene_clusters[[cluster]]
# Remove "MT-" or "mt-" not recognized in AnnotationDbi::select()
query = gsub("MT-", "", query)
query = gsub("mt-", "", query)
# Convert gene id to EntrezId format
suppressMessages(query_entrezid <- AnnotationDbi::select(org_db,
keys = query,
columns = c("ENTREZID", "SYMBOL"),
keytype = "SYMBOL"))
print_msg(msg_type = "DEBUG",
msg = paste0("Cluster ",
cluster,
"\n",
"\t|--Number of gene names converted to EntrezId : ",
length(which(!is.na(query_entrezid[, "ENTREZID"]))),
"\n",
"\t|--Number of gene names not converted to EntrezId : ",
length(which(is.na(query_entrezid[, "ENTREZID"]))),
"\n\t|-- Calling enrichGO."))
# Enrichment analysis using GO database
enrich_go_res <- clusterProfiler::enrichGO(query_entrezid[,"ENTREZID"],
OrgDb = go_name,
ont = ontology,
readable = TRUE)
# Store results in the ClusterSet object
object@gene_cluster_annotations[[cluster]] <- enrich_go_res
}
return(object)
}
)
############################################################################
## Visualization of enrichment analyses results from a ClusterSet object
############################################################################
#' @title
#' Visualization of enrichment analyses results from a ClusterSet object
##' @description
#' Display enrichment analyses results for a ClusterSet object.
#' @param object A \code{ClusterSet} object.
#' @param clusters A vector of cluster id to plot.
#' @param type A vector providing the type of plot.
#' @param nb_terms An integer indicating the number of terms in the plot.
#' @param terms_size An interger indicating the wrap length of terms
#' @param font_size The font size for the x axis.
#'
#' @return A \code{ClusterSet} object
#' @method viz_enrich ClusterSet
#' @export viz_enrich
#' @keywords internal
setGeneric("viz_enrich",
function(object,
clusters = "all",
type = "dotplot",
nb_terms = 20,
terms_size = 50,
font_size=4) {
standardGeneric("viz_enrich")
})
#' @title
#' Visualization of enrichment analyses results from a ClusterSet object
##' @description
#' Display enrichment analyses results for a ClusterSet object.
#' @param object A \code{ClusterSet} object.
#' @param clusters A vector of cluster id to plot.
#' @param type A vector providing the type of plot.
#' @param nb_terms An integer indicating the number of terms in the plot.
#' @param terms_size An interger indicating the wrap length of terms
#' @param font_size The font size for the x axis.
#'
#' @return A \code{ClusterSet} object
#' @method viz_enrich ClusterSet
#' @export viz_enrich
#'
#' @examples
#' load_example_dataset("8028226/files/pbmc3k_medium_clusters_enr_sub")
#' viz_enrich(pbmc3k_medium_clusters_enr_sub, cluster = 1, nb_terms = 5)
setMethod("viz_enrich",
signature(object = "ClusterSet"),
function(object,
clusters = "all",
type = c("dotplot", "barplot"),
nb_terms = 20,
terms_size = 50,
font_size=4) {
## Check format object arg
check_format_cluster_set(object)
if (length(clusters) == 1){
if (clusters == "all"){
clusters <- object@gene_clusters_metadata$cluster_id
}
}
list_plot <- list()
list_dotplot <- list()
list_barplot <- list()
for (cur_cluster in clusters) {
# Check if the current cluster id provided exists
if (!(cur_cluster %in% object@gene_clusters_metadata$cluster_id)) {
stop(paste0("Cluster ", cur_cluster, " doesn't exist."))
}
# Check if there is a result provided by enrich_go function for the current cluster
if(nrow(object@gene_cluster_annotations[[cur_cluster]]@result) == 0){
print_msg(msg_type = "WARNING",
msg = paste0("No functional enrichment analysis results for cluster ", cur_cluster, ".")) #Continue through the next cluster without plotting
} else {
# Print an informative message to announce plot of the results of the current cluster
print_msg(msg_type = "INFO",
msg = paste0("Plot enrichment analysis results for cluster ", cur_cluster))
# Create a ggplot - dotplot
if ("dotplot" %in% type){
if(object@gene_cluster_annotations[[cur_cluster]]@ontology == "GOALL"){
dot_plot <- enrichplot::dotplot(object@gene_cluster_annotations[[cur_cluster]],
split="ONTOLOGY",
showCategory=nb_terms,
label_format = terms_size)
dot_plot <- dot_plot + facet_grid(ONTOLOGY~., scales="free")
dot_plot <- dot_plot + ggtitle(paste0("Gene cluster ", cur_cluster))
dot_plot <- dot_plot + theme(axis.text = element_text(size=font_size))
} else {
dot_plot <- enrichplot::dotplot(object@gene_cluster_annotations[[cur_cluster]],
showCategory=nb_terms,
label_format = terms_size)
dot_plot <- dot_plot + ggtitle(paste0("Gene cluster ", cur_cluster))
dot_plot <- dot_plot + theme(axis.text = element_text(size=font_size))
}
list_dotplot <- append(list_dotplot, list(dot_plot))
}
# Create a ggplot - barplot
if ("barplot" %in% type){
if(object@gene_cluster_annotations[[cur_cluster]]@ontology == "GOALL"){
bar_plot <- graphics::barplot(object@gene_cluster_annotations[[cur_cluster]],
split="ONTOLOGY",
showCategory=nb_terms,
label_format = terms_size)
bar_plot <- bar_plot + facet_grid(ONTOLOGY~., scales="free")
bar_plot <- bar_plot + ggtitle(paste0("Gene cluster ", cur_cluster))
} else {
bar_plot <- graphics::barplot(object@gene_cluster_annotations[[cur_cluster]],
showCategory=nb_terms,
label_format = terms_size)
bar_plot <- bar_plot + ggtitle(paste0("Gene cluster ", cur_cluster))
}
list_barplot <- append(list_barplot, list(bar_plot))
}
list_plot <- append(list_barplot, list_dotplot)
}
}
return(invisible(lapply(list_plot, print)))
}
)
############################################################################
## plot_clust_enrichments()
############################################################################
#' @title Plot Cluster Enrichments
#'
#' @title Display cluster enrichments for gene GO terms.
#'
#' @param object An object of class "ClusterSet".
#' @param stat_shown A character string specifying the statistic to be shown. Must be one of "qvalue", "p.adjust", or "pvalue". Defaults to "qvalue".
#' @param nb_go_term An integer specifying the number of top GO terms to select per cluster. Defaults to 6. Note that some terms may be shared between clusters.
#' @param gradient_palette A vector of colors specifying the gradient palette for the color scale. Defaults to colors_for_gradient("Magma").
#' @param label_fun A function to customize the labels of the GO terms. Defaults to NULL.
#' @param term_order A vector specifying the desired order of the GO terms. Defaults to NULL.
#' @param return_mat Logical. Do not return a ggplot object but the prepared matrix.
#' @param floor A value for flooring the statistic.
#'
#' @return A ggplot object displaying the cluster enrichments for the GO terms.
#'
#' @seealso
#' \code{\link{enrich_go}}
#' \code{\link{colors_for_gradient}}
#'
#' @importFrom ggplot2 ggplot
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 geom_point
#' @importFrom ggplot2 scale_color_gradientn
#' @importFrom ggplot2 theme
#' @importFrom ggplot2 element_blank
#' @importFrom ggplot2 element_rect
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 xlab
#' @importFrom ggplot2 ylab
#' @importFrom ggplot2 guide_colourbar
#' @importFrom ggplot2 facet_grid
#' @importFrom reshape2 melt
#' @importFrom grid unit
#' @export plot_clust_enrichments
#' @examples
#' load_example_dataset("8028226/files/pbmc3k_medium_clusters_enr_sub")
#' plot_clust_enrichments(pbmc3k_medium_clusters_enr_sub)
#' @keywords internal
setGeneric("plot_clust_enrichments",
function(object,
stat_shown=c("qvalue", "p.adjust",
"pvalue"),
nb_go_term=3,
gradient_palette=colors_for_gradient("Ju1"),
label_fun=NULL,
term_order=NULL,
return_mat=FALSE,
floor=NULL) {
standardGeneric("plot_clust_enrichments")
})
#' @title Plot Cluster Enrichments
#'
#' @title Display cluster enrichments for gene GO terms.
#'
#' @param object An object of class "ClusterSet".
#' @param stat_shown A character string specifying the statistic to be shown. Must be one of "qvalue", "p.adjust", or "pvalue". Defaults to "qvalue".
#' @param nb_go_term An integer specifying the number of top GO terms to select per cluster. Defaults to 6. Note that some terms may be shared between clusters.
#' @param gradient_palette A vector of colors specifying the gradient palette for the color scale. Defaults to colors_for_gradient("Magma").
#' @param label_fun A function to customize the labels of the GO terms. Defaults to NULL.
#' @param term_order A vector specifying the desired order of the GO terms. Defaults to NULL.
#' @param return_mat Logical. Do not return a ggplot object but the prepared matrix.
#' @param floor A value for flooring the statistic.
#'
#' @return A ggplot object displaying the cluster enrichments for the GO terms.
#'
#' @seealso
#' \code{\link{enrich_go}}
#' \code{\link{colors_for_gradient}}
#'
#' @importFrom ggplot2 ggplot
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 geom_point
#' @importFrom ggplot2 scale_color_gradientn
#' @importFrom ggplot2 theme
#' @importFrom ggplot2 element_blank
#' @importFrom ggplot2 element_rect
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 xlab
#' @importFrom ggplot2 ylab
#' @importFrom ggplot2 guide_colourbar
#' @importFrom ggplot2 facet_grid
#' @importFrom reshape2 melt
#' @importFrom grid unit
#' @export plot_clust_enrichments
#' @examples
#' load_example_dataset("8028226/files/pbmc3k_medium_clusters_enr_sub")
#' plot_clust_enrichments(pbmc3k_medium_clusters_enr_sub)
setMethod("plot_clust_enrichments",
signature(object = "ClusterSet"),
function(object,
stat_shown=c("qvalue", "p.adjust",
"pvalue"),
nb_go_term=6,
gradient_palette=colors_for_gradient("Magma"),
label_fun=NULL,
term_order=NULL,
return_mat=FALSE,
floor=NULL) {
stat_shown <- match.arg(stat_shown)
print_msg(paste0("Using ", stat_shown, " as ordering statistic"), msg_type = "DEBUG")
print_msg("Checking object format", msg_type = "DEBUG")
check_format_cluster_set(object)
if(is.null(object@gene_cluster_annotations) | length(object@gene_cluster_annotations)==0){
print_msg("Please run enrich_go() first.", msg_type = "STOP")
}
print_msg("Extracting and formating data", msg_type = "DEBUG")
print_msg(paste0("Argument nb_go_term = ", nb_go_term), msg_type = "DEBUG")
result_go <- lapply(object@gene_cluster_annotations, "slot", "result")
order_mat <- function(x, column){x[order(x[,column]), ]}
result_go <- lapply(result_go, order_mat, stat_shown)
add_rank <- function(x, nb_go_term){x$rank <- 1:nrow(x);
x$rank[(nb_go_term + 1):nrow(x)] <- NA; return(x)}
result_go <- lapply(result_go, add_rank, nb_go_term)
for(i in 1:length(result_go)){
result_go[[i]]$cluster <- names(result_go)[i]
}
print_msg("Merging results", msg_type = "DEBUG")
m <- do.call('rbind', result_go)
go_retained <- unique(m$Description[!is.na(m$rank)])
m <- m[m$Description %in% go_retained, ]
if(!is.null(term_order)){
term_order_selected <- intersect(term_order, m$Description)
m <- m[m$Description %in% term_order_selected, ]
}
m$stat <- -log10(m[, stat_shown])
m <- m[order(m$cluster, m$rank), ]
m$is_top <- m$stat
m$is_top[is.na(m$rank)] <- "No"
m$is_top[!is.na(m$rank)] <- "Yes"
m$is_top <- as.factor(m$is_top)
if(return_mat){
return(m)
}
if(!is.null(label_fun)){
m$Description <- sapply(m$Description, label_fun)
}
if(!is.null(term_order)){
desc_levels <- sapply(term_order_selected, label_fun)
}else{
desc_levels <- rev(unique(m$Description))
}
m$Description <- factor(m$Description, levels=desc_levels, ordered = TRUE)
m$cluster <- factor(m$cluster,
levels=as.character(sort(as.numeric(unique(m$cluster)))),
ordered=TRUE)
cluster <- go_term <- stat <- size <- rank <- Counts <- NULL
print_msg("Melting...", msg_type = "DEBUG")
m_melt <- reshape2::melt(m, id.vars = c("Description", "Count", "stat", "cluster", "rank", "is_top"))
colnames(m_melt) <- c("go_term", "Counts", "stat", "cluster", "rank", "is_top", "ID", "value")
m_gg <- m_melt[order(m_melt[,"cluster"]), ]
print_msg("Running ggplot...", msg_type = "DEBUG")
if(!is.null(floor))
m_gg$stat[m_gg$stat > floor] <- floor
ggplot(data=m_gg, mapping=aes(x=cluster, y = go_term, color=stat, size=Counts)) +
ggplot2::theme_bw() +
geom_point(shape=16) +
ggplot2::scale_color_gradientn(colours = gradient_palette,
guide = guide_colourbar(title=paste0("-log10(", stat_shown, ")"), barwidth = 0.75, barheight = 5)) +
ggplot2::theme(
axis.ticks.y = ggplot2::element_blank(),
axis.ticks.x = ggplot2::element_blank(),
panel.spacing = grid::unit(0.0, "lines"),
panel.border = ggplot2::element_blank(),
strip.background.x = ggplot2::element_rect(fill = "#444444", colour = "white"),
strip.text.x = ggplot2::element_text(colour = "white", angle = 0),
axis.text.x = ggplot2::element_blank()
) +
ggplot2::xlab('Clusters') +
ggplot2::ylab('Go terms') +
ggplot2::facet_grid(~cluster, scale="free_x")
})
############################################################################
## plot_markers_to_clusters
############################################################################
#' @title Map cell markers to clusters (hypergeometric/jaccard)
#'
#' @description
#' Use jaccard index and hypergeometric to map a list of cell markers to gene clusters.
#'
#' @param object A ClusterSet object.
#' @param markers A list of cell markers.
#' @param gradient_palette A vector of colors for the gradient palette.
#' @param background The background for hypergeometric test. Typically the list of genes annotated in the considered gene ontology (e.g. BP).
#' @param cell_type_order The order of cell types. Can be "None" or "hclust". Not supported at the moment.
#' @param only_specific Logical value indicating whether to show only specific markers.
#' @param jaccard_scale The scale factor for Jaccard coefficient.
#'
#' @return A plot showing the relationship between markers and clusters.
#'
#' @export
#' @import ggplot2
#' @importFrom reshape2 melt
#' @importFrom grid unit
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 scale_color_gradientn
#' @importFrom ggplot2 element_blank
#' @importFrom ggplot2 element_line
#' @importFrom ggplot2 element_rect
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 xlab
#' @importFrom ggplot2 ylab
#' @importFrom ggplot2 facet_grid
#' @importFrom ggplot2 guide_colourbar
#' @importFrom ggplot2 geom_point
#'
#' @examples
#' load_example_dataset("8028126/files/pbmc3k_medium_clusters_enr")
#' library(clustermole)
#' m <- clustermole::clustermole_markers(species = "hs")
#' markers <- split(m$gene[m$organ == "Blood" & m$species == "Human"],
#' m$celltype[m$organ == "Blood" & m$species == "Human"])
#' plot_markers_to_clusters(pbmc3k_medium_clusters_enr, markers=markers)
#' @keywords internal
setGeneric("plot_markers_to_clusters",
function(object,
markers=NULL,
gradient_palette=colors_for_gradient("Magma"),
background=NULL,
cell_type_order=c("None", "hclust"),
only_specific=FALSE,
jaccard_scale=10) {
standardGeneric("plot_markers_to_clusters")
})
#' @title Map cell markers to clusters (hypergeometric/jaccard)
#'
#' @description
#' Use jaccard index and hypergeometric to map a list of cell markers to gene clusters.
#'
#' @param object A ClusterSet object.
#' @param markers A list of cell markers.
#' @param gradient_palette A vector of colors for the gradient palette.
#' @param background The background for hypergeometric test. Typically the list of genes annotated in the considered gene ontology (e.g. BP).
#' @param cell_type_order The order of cell types. Can be "None" or "hclust". Not supported at the moment.
#' @param only_specific Logical value indicating whether to show only specific markers.
#' @param jaccard_scale The scale factor for Jaccard coefficient.
#'
#' @return A plot showing the relationship between markers and clusters.
#'
#' @export
#' @import ggplot2
#' @importFrom reshape2 melt
#' @importFrom grid unit
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 scale_color_gradientn
#' @importFrom ggplot2 element_blank
#' @importFrom ggplot2 element_line
#' @importFrom ggplot2 element_rect
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 xlab
#' @importFrom ggplot2 ylab
#' @importFrom ggplot2 facet_grid
#' @importFrom ggplot2 guide_colourbar
#' @importFrom ggplot2 geom_point
#' @importFrom stats as.formula
#'
#' @examples
#' load_example_dataset("8028126/files/pbmc3k_medium_clusters_enr")
#' library(clustermole)
#' m <- clustermole_markers(species = "hs")
#' markers <- split(m$gene[m$organ == "Blood" & m$species == "Human"],
#' m$celltype[m$organ == "Blood" & m$species == "Human"])
#' plot_markers_to_clusters(pbmc3k_medium_clusters_enr, markers=markers)
setMethod("plot_markers_to_clusters",
signature(object = "ClusterSet"),
function(object,
markers=NULL,
gradient_palette=c("white", colors_for_gradient("De2")),
background=NULL,
cell_type_order=c("None", "hclust"),
only_specific=FALSE,
jaccard_scale=10) {
cell_type_order <- match.arg(cell_type_order)
if(is.null(background)){
print_msg("WARNING: Running without a background but you should probably provide one.", msg_type = "INFO")
}
print_msg("Checking object format", msg_type = "DEBUG")
check_format_cluster_set(object)
if(length(object@gene_clusters) == 0){
print_msg("No cluster found in the object.", msg_type = "STOP")
}
if(!is.list(markers))
print_msg("The markers argument should be a list.", msg_type = "STOP")
if(is.null(markers) | length(markers)==0){
print_msg("No marker found.", msg_type = "STOP")
}
if(only_specific){
for(i in 1:length(markers)){
cell_marker_curr <- markers[[i]]
all_other_marker <- unique(unlist(cell_marker_curr[(1:length(markers))[-i]]))
cell_marker_curr <- cell_marker_curr[!cell_marker_curr %in% all_other_marker]
markers[[i]] <- cell_marker_curr
}
nb_cell_markers <- lapply(markers, length)
markers <- markers[nb_cell_markers > 0]
}
comparison_jaccard <- compare_genesets(object@gene_clusters,
markers,
background,
stat = "jaccard")
comparison_jaccard_m <- reshape2::melt(comparison_jaccard)
comparison_hyper <- compare_genesets(object@gene_clusters,
markers,
background,
stat = "hypergeom")
comparison_hyper_m <- reshape2::melt(comparison_hyper)
comparisons_m <- comparison_jaccard_m
colnames(comparisons_m) <- c("cluster", "marker", "jaccard")
comparisons_m$pvalue <- comparison_hyper_m$value
comparisons_m$pvalue[comparisons_m$pvalue < 1e-320] <- 1e-320
comparisons_m$pvalue <- -log10(comparisons_m$pvalue)
comparisons_m$cluster <- factor(comparisons_m$cluster)
comparisons_m$marker <- factor(comparisons_m$marker)
m_gg <- comparisons_m[order(comparisons_m[,"cluster"],
comparisons_m[,"marker"],
comparisons_m[,"pvalue"]), ]
cluster <- marker <- pvalue <- jaccard <- NULL
print_msg("Running ggplot...", msg_type = "DEBUG")
ggplot(data=m_gg, mapping=aes(x=cluster, y = marker, color=pvalue, size=jaccard*jaccard_scale)) +
ggplot2::theme_bw() +
geom_point(shape=15) +
ggplot2::scale_color_gradientn(colours = gradient_palette,
guide = guide_colourbar(title=paste0("-log10(pvalue)"), barwidth = 0.75, barheight = 5)) +
ggplot2::theme(
axis.ticks.x = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_line(color="gray"),
panel.spacing = grid::unit(0.0, "lines"),
panel.border = element_rect(size=0.1, color="gray"),
panel.grid = element_blank(),
strip.background.x = ggplot2::element_rect(fill = "#444444", colour = "white"),
strip.text.x = ggplot2::element_text(colour = "white", angle = 0),
axis.text.x = ggplot2::element_blank()
) +
ggplot2::xlab('Clusters') +
ggplot2::ylab('Markers') +
ggplot2::facet_grid(~cluster, scale="free_x")
})
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.