##################################################################################
#' Generate profile heatmap
#'
#' This function generate a profile heatmap using EnrichedHeatmap function. The genes
#' are clustered into different groups based on the profile using k-means clustering
#' algorithm. Optionally, the prior clustering information can be provided through
#' argument geneGroups
#'
#' @param profileMat profile matrix of class \code{normalizedMatrix}
#' @param signalName name of the signal track
#' @param profileColor Color function for profile heatmap generated by \code{colorRamp2}
#' @param columnTitle title to be placed on the column. Default: signalName
#' @param createClusters Logical: whether or not to do kmeans clustering. If FALSE,
#' clusterStorePath must be provided which has the cluster data. Default: FALSE
#' @param numClust number of clusters to generate using k-means clustering. Default: 7
#' @param geneGroups Either a dataframe or a file path with cluster information or NULL.
#' Two columns must be present in this dataframe: \code{geneId, cluster}. In case it is NULL, heatmap
#' rows are not split into different groups.
#' @param showAnnotation Logical: to show or hide top annotation. Default: TRUE
#' @param clusterColor Color assignment object for cluster colors. If nothing provided, colors
#' are generated internally. Default: NULL
#' @param ylimFraction A numeric vector with one or two elements for deciding the \code{ylimit} of top
#' profile heatmap annotation. If the value is single number, it has to be floating point number
#' to extract the quantile and use limit \code{[0, quantile(x)]}. If the value is numeric vector of
#' length 2, the two numbers are used as lower and upper limits for ylimit of top annotation.
#' Default: NULL.
#' @param posLineGpar \code{pos_line_gp} parameter of \code{EnrichedHeatmap}. Default:
#' \code{gpar(lty = 2, alpha = 0.4, lwd = 0.5)}
#' @param rasterPar rasterization parameters as list() with two elements: use, qual. Default:
#' \code{list(use = TRUE, qual = 5)}. \code{rasterPar$use} is used for \code{use_raster} argument and
#' \code{rasterPar$qual} is used for \code{raster_quality} argument of \code{Heatmap} function.
#' @param ... Other arguments for EnrichedHeatmap function
#'
#' @return A list object with following elements:
#' \itemize{
#' \item heatmap: Enriched heatmap
#' \item rowGroupHt: row group heatmap
#' \item legend: Annotation legend for clusters
#' \item cluster: dataframe with cluster information
#' \item clusterColor: cluster color assignment information
#' \item profileColor: Heatmap profile color function generated by colorRamp2
#' \item ylim: ylimit used for the top annotation Y axis limit
#' \item kmeans: if kmeans clustering was performed, a kmeans clustering result object
#' }
#'
#' @export
#'
#' @examples NA
#'
profile_heatmap <- function(profileMat,
signalName,
profileColor,
columnTitle = NULL,
createClusters = FALSE,
numClust = 7,
geneGroups = NULL,
showAnnotation = TRUE,
clusterColor = NULL,
ylimFraction = NULL,
posLineGpar = gpar(lty = 2, alpha = 0.4, lwd = 0.5),
rasterPar = list(use = TRUE, qual = 5),
...){
cat("Generating profile heatmap for sample", signalName, "\n")
if ( missing(columnTitle) || is.null(columnTitle) ) {
columnTitle <- signalName
}
if (! is.null(geneGroups) ){
if( (!is.data.frame(geneGroups)) && isFALSE(createClusters) ){
if(! file.exists(geneGroups) ){
stop("geneGroups argument has to be either NULL or a dataframe or a file which has cluster information")
}
}
}
## hirerchical clustering
# dend = hclust(dist(profileMat), method = "ward.D")
# clusterCut = cutree(dend, numClust)
# clusterData = data.frame(geneId = names(clusterCut), cluster = clusterCut, stringsAsFactors = F, row.names = names(clusterCut))
clusterData <- NULL
lgd <- NULL
clusterHt <- NULL
if( isTRUE(createClusters) ){
## do k-means clustering or read existing clustering information
profileKm <- profile_matrix_kmeans(mat = profileMat,
km = numClust,
clustFile = geneGroups,
name = signalName)
clusterData <- profileKm$geneClusters
} else if( ! is.null(geneGroups) ){
if(is.data.frame(geneGroups)){
## read the clustering information from dataframe
# cat("Using the provided clusterData...\n")
clusterData <- as.data.frame(geneGroups)
rownames(clusterData) <- NULL
} else if(file.exists(geneGroups)) {
## read clustering information from cluster file
# cat("Reading k-means clustering information sample", signalName, "\nNumber of clusters:", numClust, "\n")
clusterData <- suppressMessages(readr::read_tsv(file = geneGroups))
}
} else if( is.null(geneGroups) ){
cat("Heatmap wont be split as geneGroups information is NULL\n")
}
if( !is.null(clusterData) ){
clusterData <- dplyr::select(clusterData, geneId, cluster) %>%
tibble::column_to_rownames(var = "geneId") %>%
as.data.frame()
## make sure that the rows of clusterdata and profile matrix are in same order
## if the order is different, the splitting will be changed
if(! all(rownames(profileMat) == rownames(clusterData) ) ){
stop("Row order for the profile matrix and clustering dataframe does not match")
}
## set colors for clusters. IMP to sort as they will be arranged on heatmap in order
if(is.null(clusterColor)){
clusterNames <- levels(clusterData$cluster)
# cat("Colors for clusters not provided. Generating new colors for clusters: ", clusterNames, "\n")
# clusterColor <- structure(rainbow(length(clusterNames)), names = clusterNames)
clusterColor <- structure(
chipmineEnv$customColors[1:length(clusterNames)],
names = clusterNames)
## if there are too many clusters, use rainbow colors
if(length(clusterNames) > length(chipmineEnv$customColors)){
clusterColor <- structure(rainbow(length(clusterNames)), names = clusterNames)
}
} else{
## modify the clusterColor: remove the cluster_x which is not present in the clusterData
clusterNames <- levels(clusterData$cluster)
clusterColor <- clusterColor[clusterNames]
}
# print(clusterColor)
## annotation legend
lgd <- Legend(at = names(clusterColor), title = "Clusters",
type = "lines", legend_gp = gpar(col = clusterColor))
## use heatmap instead of rowAnnotation for cluster
clusterHt <- ComplexHeatmap::Heatmap(
matrix = as.matrix(clusterData["cluster"]),
name = "clusterAn",
col = clusterColor,
width = unit(0.3, "cm"),
show_row_names = FALSE,
show_heatmap_legend = FALSE
)
} else if( is.null(clusterData)){
## if there is no splitting of the rows, use single color
clusterColor <- "red"
}
## ylimit for the annotation
ylimit <- NULL
if(!is.null(ylimFraction) && length(ylimFraction) == 1){
ylimit <- c(0, quantile(profileMat, ylimFraction, names = FALSE))
} else if(length(ylimFraction) == 2){
ylimit <- ylimFraction
}
## set the top annotation
topAnn <- NULL
if(isTRUE(showAnnotation)){
topAnn <- HeatmapAnnotation(
lines = anno_enriched(gp = gpar(col = clusterColor),
show_error = F,
ylim = ylimit
),
height = unit(2, "cm")
)
}
## axis_names:
axis_names <- NULL
if(isTRUE(attributes(profileMat)$target_is_single_point)){
## when region length == 1
axis_names[1] <- attributes(profileMat)$extend[1] * -1
axis_names[2] <- attributes(profileMat)$target_name
axis_names[3] <- attributes(profileMat)$extend[2]
} else{
## when region length > 1
axis_names[1] <- attributes(profileMat)$extend[1] * -1
axis_names[2] <- "start"
axis_names[3] <- "end"
axis_names[4] <- attributes(profileMat)$extend[2]
}
## generate profile heatmap
ht <- EnrichedHeatmap(
mat = profileMat,
col = profileColor,
name = signalName,
axis_name_rot = 90,
axis_name = axis_names,
column_title = columnTitle,
split = clusterData,
top_annotation = topAnn,
row_title_rot = 0,
pos_line_gp = posLineGpar,
use_raster = rasterPar$use,
raster_quality = rasterPar$qual,
...
)
## reconvert the factor to character
clusterData$cluster <- as.character(clusterData$cluster)
clusterData$geneId <- rownames(clusterData)
returnList <- list(
"heatmap" = ht,
"rowGroupHt" = clusterHt,
"legend" = lgd,
"cluster" = clusterData,
"clusterColor" = clusterColor,
"profileColor" = profileColor,
"ylim" = ylimit
)
if(isTRUE(createClusters)){
returnList[["kmeans"]] <- profileKm$km
}
##returns: profile heatmap, annotation legend, clusterData
return(returnList)
}
##################################################################################
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.